0

I have a Sturm-Liouville operator:

$$\mathcal{Ly(x)}=-\frac{\mathrm{d}}{\mathrm{d}x}\left(p(x)\frac{\mathrm{d}}{\mathrm{d}x}y(x)\right)-q(x)y(x)$$

but I do not yet know the parameter functions $p(x)$ and $q(x)$, nor the arbitrary function $y(x)$. I define the operator in Mathematica:

sl[x_]:=-D[p[x]*D[y[x],x],x]-q[x]*y[x]

What if now I want to specify what $p$, $q$ and $y$ should be? Let's say $p(x)=1$, $q(x)=1$ and $y(x)=\sin(x)$. I do:

res[x_]:=sl[x]/.{p[x]->1,q[x]->0}/.{y[x]->Sin[x]}

To check if this is correct and I can plot res, I do:

Plot[res[x],{x,0,Pi}]

I am getting an empty plot:

enter image description here

but I expected a $-\frac{\mathrm{d}^2}{\mathrm{d}^2x}\left(\sin(x)\right)=\sin(x)$ plot.

What am I doing wrong?

Seemingly reasonable another failed attempts:

1)

sl[p_, q_, y_,x_]:=-D[p[x]*D[y[x],x],x]-q[x]*y[x]
res[x_]:=sl[p[x], q[x], y[x], x]/.{p[x]->1,q[x]->0}/.{y[x]->Sin[x]}

Based on the Currying part of this answer:

sl[x_][p_,q_]:=-D[p[x]*D[y[x],x],x]-q[x]*y[x]
p[x_]:=1
q[x_]:=0
res[x_]:=sl[x][p, q]

Based on comments:

sl[x_]:=((-D[p[x]*D[y[x],x],x]-q[x]*y[x])/.{p[x]->1,q[x]->0})
res[x_]:=sl/.{y[x]->Sin[x]}

Specifying y as well in the first line:

sl[x_]:=((-D[p[x]*D[y[x],x],x]-q[x]*y[x])/.{p[x]->1,q[x]->0,y[x]->Sin[x]})

when plotting it (Plot[sl[x],{x,0,Pi}]), I still get errors and no plot.

zabop
  • 137
  • 6
  • What do you get if you just execute res[x] for an unspecified x? Does the result help you explain the behavior? – Marius Ladegård Meyer Jan 17 '22 at 19:54
  • Thanks for the suggestion! I get just -p'[x]y'[x]-y''[x], which suggests I think that I couldn't pass the value for p properly for some reason, otherwise it would not be shown as p. – zabop Jan 17 '22 at 19:59
  • I tried writing res:= ... according to this answer. Instead of VariationalD[(a[x] - s[x])^2, s[x], x], I have sl[x]. Maybe the problem is that the p and q dependence is not explicitly statet in my res definition, trying to fix that... – zabop Jan 17 '22 at 20:05
  • Yeah, the problem is that in your res definition sl[x] is evaluated before the substitution with /. happens, both when x is a number and when it is a symbol. In my opinion it is sl that should take p and q as input functions, perhaps as pure Functions? – Marius Ladegård Meyer Jan 17 '22 at 20:24
  • Thanks you! I added my latest attempts to the question (number 3 is based on these comments), I haven't found a solution yet but keep working on it... – zabop Jan 17 '22 at 20:31
  • attempt 4 also added. – zabop Jan 17 '22 at 20:34

2 Answers2

1
Clear["Global`*"]

sl[p_, q_, y_, xv_, var_ : x] := 
  (-D[p*D[y, var], var] - q*y) /. var -> xv

p[x_] = 1;
q[x] = 0;
y[x] = Sin[x];
res[xv_] := sl[p[x], q[x], y[x], xv]

Plot[Evaluate@res[x], {x, 0, Pi}]

enter image description here

Another example,

p[x_] = x^2;
q[x] = Sqrt[x];
y[x] = x*Sin[x];
res[xv_] := sl[p[x], q[x], y[x], xv]

Note the equivalent use of the option Evaluated rather than Evaluate

Plot[res[x], {x, 0, Pi}, Evaluated -> True]

enter image description here

Bob Hanlon
  • 157,611
  • 7
  • 77
  • 198
  • Thank you, I indeed get the plots. Here: sl[p_, q_, y_, xv_, var_ : x] := (-D[p*D[y, var], var] - q*y) /. var -> xv what's happening is that we define the right hand side of the formula using var, and then we set var to xv by /. var -> xv, is that correct? Probably var_ : x is also a key element on the left hand side. If you could add some explanation to this or direct me to the right direction researching this, it'd be much appreciated. – zabop Jan 17 '22 at 21:21
  • var_ : x is used to indicate that var is an optional argument with a default value of x. This facilitates handling functions of another variable such as z. The var -> xv insures that the derivatives are taken prior to a numeric value being assigned to the variable. – Bob Hanlon Jan 17 '22 at 21:56
1

Straight replacing the way you have done does not honor derivatives. You can do it this way:

sl[x_] = -D[p[x]*D[y[x], x], x] - q[x]*y[x]

res[x_] = sl[x] /. {p -> (1 &), q[x] -> 0} /. {y -> (Sin[#] &)}
(* Sin[x] *)

And you will get the plot you expect.

If you don't want to use pure functions you can do it this way:

sl[x] /. {p -> Function[{x}, 1], q -> Function[{x}, 0], 
  y -> Function[{x}, Sin[x]]}
(* Sin[x] *)

It is not necessary to make q a function in this case because sl contains no derivatives of q. But maybe you don't know that up front, so it's OK to do so. In any case if you want Mathematica to replace the function as well as its derivatives, you need to tell Mathematica it's dealing with a function. Of course you can always specify p[x], p'[x], etc. separately.

Bill Watts
  • 8,217
  • 1
  • 11
  • 28