Using the Equations We Derived in General

It is better to work the series solution out by hand and determine the underlying function using Mathematica whenever possible. However, if one wanted to create truncated solutions automatically using Mathematica and the equations we worked out in class, one could. Here is one way to do it.

Differential Equations

Mathematica's DSolve

Here is a differential equation, defined with coefficients left unspecified:

L[y_,P_,Q_,R_]:=P*D[y,{x,2}]+Q*D[y,{x,1}]+R*y

We can use this to construct any differential equation we want:

L[y[x],1,0,-1]==0

DSolve[L[y[x],1,0,-1]==0,y[x],x]

DSolve[L[y[x],(1-x),x,x]==0,y[x],x]

Notice the special functions which arise. This is very common with variable coefficient linear differential equations.

Example 1: Module to Find Solutions About Ordinary Point

Here is a module that does the series solution about an ordinary point, by using the recursion relation we worked out in general. Notice we need to use Mathematica to determine the coefficients of the Taylor series expansions for p(x) and q(x).

SSOP[P_,Q_,R_,mx_,x0_,x_]:=Module[{a,p,q,k,solution,y},

p[k_]:=p[k]=(1/(k!))D[(Q/P),{x,k}]/.x->x0;
q[k_]:=q[k]=(1/(k!))D[(R/P),{x,k}]/.x->x0;

a[0]=c1;
a[1]=c2;
a[n_]:=a[n]=-(1/((n)(n-1)))((∑_(k=0)^(n-2))(((n-2)-k+1)p[k]a[(n-2)-k+1]+q[k] a[(n-2)-k]));
solution=(∑_(n=0)^mx)a[n](((x-x0))^n);
Return[Collect[solution,{c1,c2}]]]

We can use the module to solve differential equations:

SSOP[1,0,1,12,0,x]

Compare this with Mathematica's solution:

Clear[y]
DSolve[L[y[x],1,0,1]==0,y[x],x]
Series[Cos[x],{x,0,10}]
Series[Sin[x],{x,0,10}]

We can solve a wide variety of problems with this module, as long as we are expanding about an ordinary point, and check that they are really solutions by substituting back into the differential equation. Here we do two solutions, keeping more terms in the second:

sol1[x_]=SSOP[1,x,(x^2)-1,10,2,x]
Simplify[L[sol1[x],1,x,(x^2)-1]]

sol2[x_]=SSOP[1,x,(x^2)-1,12,2,x]
Simplify[L[sol2[x],1,x,(x^2)-1]]

Plot[{sol1[x]/.c1->1/.c2->2,sol2[x]/.c1->1/.c2->2},{x,0,4},PlotRange->{{0,5},{-25,5}},PlotStyle->{Thick},ImageSize->{400,Automatic}]

Plot[{Abs[((sol1[x]-sol2[x])/(sol1[x]))]/.c1->1/.c2->2},{x,0,4},PlotRange->{{0,5},{0,1}},PlotStyle->{Thick},ImageSize->{400,Automatic}]

From the above sketches, we could be confident of our solution truncated at (((x-2))^10) for 0.5 < x < 3.0.

It would be difficult, however, to pull the pattern out of our series and write it in closed summation form, but we could certainly try! If you are using a truncated series, you definitely need to consider where your solution is valid, since their might be significant truncation error as well as error from being outside the interval of convergence.

Example 2: Module to Find First Solution About Regular Singular Point

Here is a module that does the series solution about a regular singular point, by using the recursion relation we worked out in general and reduction of order (since reduction of order will always find your second solution). Notice we use Mathematica to determine the coefficients in the Taylor series expansion of (x-(Subscript[x, 0]))p(x) and (((x-(Subscript[x, 0])))^2)q( x).

Clear[y,sol]

SSRSP[P_,Q_,R_,mx_,x0_,x_]:=Module[{a,p,q,r1,solution1,solution2,k,n},

p[k_]:=p[k]=(1/(k!))D[(x-x0)*Q/P,{x,k}]/.x->x0;
q[k_]:=q[k]=(1/(k!))D[(((x-x0))^2)*R/P,{x,k}]/.x->x0;

r1=Simplify[(1/2) (1-p[0]+(Sqrt[1-2 p[0]+((p[0])^2)-4 q[0]]))];

a[0]=1;
a[n_]:=a[n]=Simplify[-(1/((n+r1)(n+r1-1)+p[0](n+r1)+q[0]))(∑_(k=1)^n)((n-k+r1)p[k]+q[k])a[n-k]];
solution1=(∑_(n=0)^mx)a[n](((x-x0))^(n+r1));

solution2=Collect[(∫Normal[Series[((E^(-∫(Q/P)\[DifferentialD]x))/(((solution1))^2)),{x,x0,mx}]]\[DifferentialD]x) *solution1,x];

Return[Collect[c1*solution1+c2*solution2,{c1,c2}]]
]

Here is the solution our module returns:

sol[x_]=SSRSP[x^2,x,(x-1),6,0,x]

Compare with the solution Mathematica finds:

DSolve[L[y[x],x^2,x,(x-1)]==0,y[x],x]

Since we did not obtain the pattern in the recursion relations, we do not have an easy way to verify our truncated series solutions match these functions (they might not, since we have found a different fundamental set of solutions). Let’s substitute our solutions back into the DE to see if they work.

L[(x-((x^2)/3)+((x^3)/24)-((x^4)/360)+((x^5)/8640)-((x^6)/302400)+((x^7)/14515200)),x^2,x,(x-1)]//FullSimplify

Not zero, but it is O((x^8)) which is the order of x we truncated at. So this is what we expect.

The second solution also works, but it is a bit more complicated to see how well it works just by substituting back into the DE due to the logarithm in the solution.

L[(-(1/2)-(1/(2 x))+((x^8)/1254113280)-((10890589 (x^9))/118513704960000)+((1320881 (x^10))/79009136640000)-((14757733 (x^11))/11377315676160000)+((313573 (x^12))/5333116723200000)-((393719 (x^13))/227546313523200000)+((42797 (x^14))/1023958410854400000)+(x^2) ((19/432)-((Log[x])/12))+(x^4) ((163/129600)-((Log[x])/1440))+(x^6) ((37/17418240)-((Log[x])/1209600))+(x^7) (-(1/19906560)+((Log[x])/58060800))+(x^5) (-(401/6220800)+((Log[x])/34560))+(x^3) (-(23/1728)+((Log[x])/96))+x ((29/144)+((Log[x])/4))),x^2,x,(x-1)]//FullSimplify

The above quantity should be close to zero in the interval where our truncated series is a good approximation to the solution to the DE.

Initial Value Problems

Mathematica's DSolve

Here is a differential equation, defined with coefficients left unspecified:

L[y_,P_,Q_,R_]:=P*D[y,{x,2}]+Q*D[y,{x,1}]+R*y

We can use this to construct any differential equation we want:

L[y[x],1,x,2]==0

DSolve[{L[y[x],1,x,2]==0,y[0]==1,y'[0]==2},y[x],x]

Notice the special functions which arise. This is very common with variable coefficient linear differential equations.

Example 1: Module to Find Solutions About Ordinary Point

Here is a module that does the series solution about an ordinary point, by using the recursion relation we worked out in general. Notice we need to use Mathematica to determine the coefficients of the Taylor series expansions for p(x) and q(x).

SSOPIVP[P_,Q_,R_,mx_,x0_,x_,y0_,y0p_]:=Module[{a,p,q,k,solution,y},

p[k_]:=p[k]=(1/(k!))D[(Q/P),{x,k}]/.x->x0;
q[k_]:=q[k]=(1/(k!))D[(R/P),{x,k}]/.x->x0;

a[0]=y0;
a[1]=y0p;
a[n_]:=a[n]=-(1/((n)(n-1)))((∑_(k=0)^(n-2))(((n-2)-k+1)p[k]a[(n-2)-k+1]+q[k] a[(n-2)-k]));
solution=(∑_(n=0)^mx)a[n](((x-x0))^n);
Return[solution]]

We can solve a wide variety of problems with this module, as long as we are expanding about an ordinary point. Here we do two solutions, keeping more terms in the second:

sol1[x_]=SSOPIVP[1,x,2,200,2,x,1,-3];

sol2[x_]=SSOPIVP[1,x,2,220,2,x,1,-3];

Plot[{sol1[x],sol2[x]},{x,-15,15},PlotRange->{{-12,12},{-6,6}},PlotStyle->{Thick},ImageSize->{400,Automatic}]

Plot[{Abs[((sol1[x]-sol2[x])/(sol1[x]))]},{x,-12,12},PlotRange->{{-12,12},{0,1}},PlotStyle->{Thick},ImageSize->{400,Automatic}]

From the above sketches, we could be confident of our solution truncated at (((x-2))^200) for -5 < x < 9.

Example 2: Module to Find First Solution About Regular Singular Point

This would be a pain to set up, so why bother? With Ordinary points, the initial conditions were easy to include. Here, not so much.