# Solve 2D Sturm-Liouville equation

I'm trying to solve the 2D SL equation for an analytic solution $F(x,t)$, with the boundary conditions $F(0,t) = F(100,t) = 0$ and $F(x, 0) = \delta(x-1)$. $$\frac{\partial F}{\partial t} = (\frac{x+1}{2})\frac{\partial^2 F}{\partial x^2} + (0.01x+1)\frac{\partial F}{\partial x} + 0.01F.$$

I tried doing

    pde2 = 0.01 F[x, t] + (0.01 x + 1) D[F[x, t], x] + ((x + 1)/2) D[
F[x, t], {x, 2}] - D[F[x, t], t] == 0
soln = DSolve[pde2, F[x, t], {x, t}]


but this does not work, and just echoes the second command.

I tried helping Mathematica by plugging in the ansatz $F(x,t) = y(x)e^{t}$, and we can then simplify it to:

 pde = -0.99 y[x] + (0.01 x + 1) D[y[x], x] + ((x + 1)/2) D[
y[x], {x, 2}] == 0
soln = DSolve[pde, y[x], {x}]


This gives us

{{y[x] ->
E^(-0.02 (1. + 1. x))
C[1] HypergeometricU[100.98, 1.98, 0.02 + 0.02 x] +
E^(-0.02 (1. + 1. x))
C[2] LaguerreL[-100.98, 0.98, 0.02 + 0.02 x]}}


However, once I add initial conditions, such as

soln = DSolve[{pde, y[0] == 0}, y[x], {x}]


the method no longer works, and prints out an empty set for solutions. Any help is appreciated.

• What're you trying to find? The analytic solution, or numeric solution? If the former, do you need the general solution, the particular solution with one initial condition defined in infinite domain, or the particular solution with one initial condition and two explicit boundary condition defined in a bounded domain? – xzczd Feb 2 '18 at 6:51
• Sorry if the problem wasn't clear - I want the analytic solution. I updated the question. – Alex Feb 2 '18 at 7:29

Analytically, using separation of variables, you take F[x, t] = X[x] T[t]. This would yield

D[T[t],t]/T[t] = ((x+1)/2X[x]) X''[x] + ((0.01x + 1)/X[x]) X'[x] + 0.01


So if you take

T[t] = Exp[m t]


then

((x+1)/2X[x]) X''[x] + ((0.01x + 1)/X[x]) X'[x] + 0.01 - m = 0


so the result is (replacing X[x] with y[x] in the above):

pde = (-m + 0.01) y[x] + (0.01 x + 1) D[y[x], x] + ((x + 1)/2) D[
y[x], {x, 2}] == 0
soln = DSolve[{pde, y[0] == 0}, y[x], {x}]

(*{{y[x] -> (0.980199 E^(-0.02 x)
C[1] (HypergeometricU[0.98 + 100. m, 1.98,
0.02 + 0.02 x] LaguerreL[-0.98 - 100. m, 0.98, 0.02] -
1. HypergeometricU[0.98 + 100. m, 1.98,
0.02] LaguerreL[-0.98 - 100. m, 0.98, 0.02 + 0.02 x]))/
LaguerreL[-0.98 - 100. m, 0.98, 0.02]}}*)


note that C[1] is determined by F'[0,0] initial condition. Remember that in x coordinate, you have 2nd order differential equation which requires 2 boundary or initial values to fully determine the answer.

All this is equivalent to taking ansatz:

F[x_,t_] := y[x] Exp[m t]


which is the most general solution you can have for this equation. Note that m must be pure imaginary to have a solution. So F[x,t] = y[x] Exp[I m t] is an oscillatory function in time. That is why you won't get any answer with m = 1 or any other integer or real numbers. However trying 'm = 0. + 1. I' would yield:

(*{{y[x] ->
0.980199 E^(-0.02 x) (C[1] HypergeometricU[0.98 + 100. I, 1.98,
0.02 + 0.02 x] + (4.35948*10^65 - 8.61911*10^65 I) C[
1] LaguerreL[-0.98 - 100. I, 0.98, 0.02 + 0.02 x])}}*)


Assuming

y'[0] = 1


the solution is

pde = (-m + 0.01) y[x] + (0.01 x + 1) D[y[x], x] + ((x + 1)/2) D[y[x],
{x, 2}] == 0
soln = DSolve[{pde, y[0] == 0, y'[0] == 1}, y[x], {x}]

(*{{y[x] -> -((0.5 E^(-0.02 x) (1. HypergeometricU[0.98 + 100. m, 1.98,
0.02 + 0.02 x] LaguerreL[-0.98 - 100. m, 0.98, 0.02] -
1. HypergeometricU[0.98 + 100. m, 1.98,
0.02] LaguerreL[-0.98 - 100. m, 0.98,
0.02 + 0.02 x]))/(-0.01 HypergeometricU[0.98 + 100. m,
1.98, 0.02] LaguerreL[-1.98 - 100. m, 1.98, 0.02] +
0.0098 HypergeometricU[1.98 + 100. m, 2.98,
0.02] LaguerreL[-0.98 - 100. m, 0.98, 0.02] +
1. m HypergeometricU[1.98 + 100. m, 2.98,
0.02] LaguerreL[-0.98 - 100. m, 0.98, 0.02]))}}*)


so

Plot[Re[(y[x] /. soln[[1]]) /. m -> I], {x, 0, 10}, Frame -> True,
FrameLabel -> {"x", "Re[y[x]"}]
Plot[Im[(y[x] /. soln[[1]]) /. m -> I], {x, 0, 10}, Frame -> True,
FrameLabel -> {"x", "Im[y[x]"}]