# How to solve Sturm-Liouville problems in SageMath?

I've conducted DuckDuckGo searches for "Sturm-Liouville sagemath" and received no usable results, and as a result I am here to ask, how can one solve Sturm-Liouville problems with SageMath? I know how to do so numerically with Julia, GNU Octave and Scilab, using matrices, but is there any way to solve at least some of them analytically in SageMath (e.g. $\frac{d^2 y}{dx^2}-xy = \lambda y$ has solutions $y_n=C\cdot \textit{Ai}(x+\lambda_n)$, where $\lambda_n$ are the zeros of the Airy $\textit{Ai}(x)$ function, if $y(0)=y(\infty)=0$)? Failing this, is there a way to numerically solve them in SageMath (without resorting to its interfaces with the languages I just mentioned)?

If it helps, here's the code I use in GNU Octave to solve the aforementioned SLE:

# Standard values are for transformation x=(y-L)/(y+L), y in 0, inf
# *lin are for x=2/L*y-1 (linear transformation)
clear all
format long e;
N=5000;
perc=0.266
Nfrag=round(perc*N);
n=0:N;
# 75 at N=1k; 162 for N=3k; 225 for N=5k;
L=226;
t=pi+n'*pi/N;
x=cos(t);
y=L*(cot(t/2)).^2;
ysub=y(2:N);
ylintrans = L/2*(x+1);
# Chebyshev polys of the first kind
T         = cos(t*n);
# Now for arrays that do not include the endpoints
xsub      = x(2:N);
Tsub      = T(2:N,:);
Usub      = diag(1./sqrt(1-xsub.^2))*sin(acos(xsub)*n);
dTsub     = Usub*diag(n);
dT        = [-(n.^2).*(-1).^(n); dTsub ; (n).^2];
dTL       = diag(2*L./((y+L).^2))*dT;
dTlin     = 2/L*dT;
D1        = dTL/T;
D1lin     = dTlin/T;
D2        = D1^2;
D2lin     = D1lin^2;
H         = D2-diag(y);
Hlin      = D2lin-diag(ylintrans);
H         = H(2:N,2:N);
Hlin      = Hlin(2:N,2:N);
[Y, Lam]  = eig(H);
[Ylin, Lamlin] = eig(Hlin);
Lam       = diag(Lam);
Lamlin    = diag(Lamlin);
[Lam, IX] = sort(Lam, 'ascend');
[Lamlin, IXlin] = sort(Lamlin, 'ascend');
Y         = Y(:,IX);
Ylin      = Ylin(:,IXlin);
aizerochk = abs(airy(0,Lam));
aizerochklin = abs(airy(0,Lamlin));
L
errrms    = sqrt(aizerochk(1:Nfrag)'*aizerochk(1:Nfrag)/(Nfrag+1))
errrmslin    = sqrt(aizerochklin(1:Nfrag)'*aizerochklin(1:Nfrag)/(Nfrag+1))
figure(1)
plot(ysub(1:3*N/5), Y(1:3*N/5,1));
figure(2)
plot(ysub(1:3*N/5), Y(1:3*N/5,2));
figure(3)
plot(ysub(1:3*N/5), Y(1:3*N/5,Nfrag));


. It uses a Chebyshev spectral method on the extrema grid, and then analyses the results.

edit retag close merge delete

I think that anything which can be computed in Julia, Octave or Scilab can also be computed in SageMath. So I would suggest to post the code you use to solve the problem in one of these languages, e.g. in Octave. In this way, you could get help to translate it to SageMath.

( 2019-07-24 04:41:53 -0600 )edit

Thanks, shall do.

( 2019-07-24 05:01:37 -0600 )edit

From airy_ai?:

   The Airy Ai function operatorname{Ai}(x) is (along with
operatorname{Bi}(x)) one of the two linearly independent standard
solutions to the Airy differential equation f''(x) - x f(x) = 0. It
is defined by the initial conditions:

operatorname{Ai}(0)=frac{1}{2^{2/3}
Gamma(frac{2}{3})},  operatorname{Ai}'(0)=-frac{
1}{2^{1/3}Gamma(frac{1}{3})}.


But Sage doesn't use it when solving the ODE...

( 2019-07-24 08:30:22 -0600 )edit

Sort by » oldest newest most voted

Sage can give us a symbolic solution in terms of Bessel functions:

sage: y=function("y")
sage: var("lambda_")
lambda_
sage: E=diff(y(x),x,2)-x*y(x)==lambda_*y(x)
sage: S=desolve(E,y(x),ivar=x, contrib_ode=True);S
[y(x) == 3/2*(2/3)^(2/3)*(2/3*I*(lambda_ + x)^(3/2))^(1/3)*_K2*sqrt(lambda_ + x)*bessel_I(1/3, 2/3*(lambda_ + x)^(3/2))/((lambda_ + x)^(3/2))^(1/3) + _K1*sqrt(lambda_ + x)*bessel_Y(1/3, 2/3*I*(lambda_ + x)^(3/2))]
sage: var("_K1, _K2")


wich is better seen/understood via \LaTeX:

$$y\left(x\right) = \frac{3 \left(\frac{2}{3}\right)^{\frac{2}{3}} \left(\frac{2}{3} i {\left(\lambda + x\right)}^{\frac{3}{2}}\right)^{\frac{1}{3}} K_{2} \sqrt{\lambda + x} I_{\frac{1}{3}}(\frac{2}{3} {\left(\lambda + x\right)}^{\frac{3}{2}})}{2 {\left({\left(\lambda + x\right)}^{\frac{3}{2}}\right)}^{\frac{1}{3}}} + K_{1} \sqrt{\lambda + x} Y_{\frac{1}{3}}(\frac{2}{3} i {\left(\lambda + x\right)}^{\frac{3}{2}})$$

I am not qualified to check this solution.

more

Incorrect, because it does not take the boundary conditions into account, nor does it take into account that lambda is an eigenvalue, not just any old parameter. The solution only satisfies the boundary conditions when lambda is a zero of the Airy Ai(x) function.

( 2019-07-24 05:32:11 -0600 )edit

Where are stated those boundary conditions ? They may be accounted for by using the ics argument of desolve (thus defining $K_1$ and _$K_2$).

The proposed solution is a solution of the differential equation, not necessarily the solution of the problem.

( 2019-07-24 08:28:09 -0600 )edit

If you look at my question the conditions are: $y(0)=y(\infty)=0$ (granted, I added that like 30 mins before you published this answer, so maybe you missed it). I didn't really expect desolve to work on this, otherwise my DuckDuckGo search probably would have detected a mention of its use for this purpose (SLEs are a different kettle of fish to usual ODEs in terms of how they are solved, so I would have expected a special mention of them, if desolve could solve them). Nonetheless, I gave it a go, here's the code I used and I placed the error in a comment: https://gist.github.com/fusion809/b3e....

( 2019-07-24 08:48:24 -0600 )edit