# How to define specify parameter functions of a Sturm-Liouville operator & plot its effect on specific functions?

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:

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]


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.

• What do you get if you just execute res[x] for an unspecified x? Does the result help you explain the behavior? Jan 17, 2022 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. Jan 17, 2022 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... Jan 17, 2022 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? Jan 17, 2022 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... Jan 17, 2022 at 20:31

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.

• Thank you. I indeed get the plot. Would you be able to add some explanation or direct me in the right direction to understand why using p -> (1 &) works, but p[x]->1 does not? Jan 17, 2022 at 21:16
• Because p[x]->1 will properly replace p[x] in sl[x], but it will not replace p'[x]. Jan 17, 2022 at 21:48
• You don't have that problem with q[x] because there are no derivatives of q[x] in sl[x]. Jan 17, 2022 at 21:50
• Thank you, that clears it up. Jan 17, 2022 at 21:54
• You are welcome. Jan 17, 2022 at 22:27
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}]


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]


• 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. Jan 17, 2022 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. Jan 17, 2022 at 21:56