Ask Your Question
1

How to solve Sturm-Liouville problems in SageMath?

asked 2019-07-24 00:17:12 -0600

Fusion809 gravatar image

updated 2019-07-24 08:04:24 -0600

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))
ain=quad("airynorm",0,100);
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 flag offensive close merge delete

Comments

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.

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

Thanks, shall do.

Fusion809 gravatar imageFusion809 ( 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...

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2019-07-24 08:30:22 -0600 )edit

1 answer

Sort by ยป oldest newest most voted
0

answered 2019-07-24 05:21:54 -0600

Emmanuel Charpentier gravatar image

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.

edit flag offensive delete link more

Comments

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.

Fusion809 gravatar imageFusion809 ( 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.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 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....

Fusion809 gravatar imageFusion809 ( 2019-07-24 08:48:24 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2019-07-24 00:17:12 -0600

Seen: 75 times

Last updated: Jul 24