Example From Class (1-x) y''[x]+x y'[x]-y[x]=0

Mathematica Solution: Infinite series (not so good)

Here is the differential equation we solved by hand in class

Clear[y,a,x]
L[y_]:=(1-x)*D[y,{x,2}]+x*D[y,{x,1}]-y
L[y[x]]==0

Let's substitute an infinite series about (Subscript[x, 0])=0  into the differential equation, and see what we can learn:

yassumed[x_]:= Sum[a[i]x^i,{i,0,Infinity}]
yassumed[x]

Substitute into the differential equation:

deq1=L[yassumed[x]]

Collect powers of x

deq2 = Expand[deq1]
deq3 = Collect[deq2,x]

OK, so we are kind of stuck. Getting the recurssion relation is easy to do by hand, but difficult to do in Mathematica (there is a way, but it uses some Mathematica commands that aren't very easy to understand. Since I want us to always understand what Mathematica is doing, we should get the recursion relation by hand (we did it in class, it wasn't difficult, although it may be a bit tedious).

Mathematica Solution: Truncated series (better, but we can't recognize the pattern doing this)

Here is the differential equation we solved by hand in class

Clear[y,a,x]
L[y_]:=(1-x)*D[y,{x,2}]+x*D[y,{x,1}]-y
L[y[x]]==0

Let's substitute a truncated series about (Subscript[x, 0])=0  into the differential equation, and see what we can learn. The series will have n terms.

n=5
yassumed[x_]:= Sum[a[i]x^i,{i,0,n}]
yassumed[x]

Substitute into the differential equation:

deq1=L[yassumed[x]]

Collect powers of x

deq2 = Expand[deq1]
deq3 = Collect[deq2,x]

Strip off the coefficients, to get equations the a[n] must satisfy, then solve the equations recursively. The a[0] and a[1] are arbitrary, which is why we do not include them in the list of variables to solve for. We will need n-1 equations for the n-1 unknown a[2],a[3],...,a[n-1],a[n]. We get these equations from the first n-1 coefficients, which means taking i from zero to n-2.

eqs=Table[Coefficient[deq3,x,i]==0,{i,0,n-2}]

vars=Table[a[i],{i,2,n}]

sol=Solve[eqs,vars]

Let's substitute these back into our solution:

ygeneral[x_]=yassumed[x]/.Flatten[sol]

Now, we can collect the solution in a[0] and a[1]:

ygeneral[x_]=Collect[ygeneral[x],{a[0],a[1]}]

And finally, get a fundamental set of solutions:

y1[x_]=Coefficient[ygeneral[x],a[0]]
y2[x_]=Coefficient[ygeneral[x],a[1]]

Notice that we get the series solution to any order we like (we can simply go back and increase n if we want. The problem with this is we might have difficulty recognizing the pattern from what Mathematica returns.

The solution y2[x] is not an infinite series, therefore it converges for all x and we don't have to worry about it any more.

The solution y1[x] is an infinite series, so we should be worried about how accurate we are with it when we truncate the infinte series.

Let's check our work:

Simplify[L[y1[x]]==0]
Simplify[L[y2[x]]==0]

These should both be identically zero--the first one isn't because we have truncated the series. Keeping more terms should get a solution that is more accurate over the radius of convergence of the series. Since the series is truncated, we don't have access to the general form of a[n], and we can't use the ratio test to determine the radius of convergence. In this case, the best thing to do is to get a numerical estimate of how accurate your solution is over a given domain.

Here is a plot of the first solution, so we can estimate how accurate it is.

Plot[{y1[x],y1[x]-x^5/120},{x,-3,4},PlotStyle->{Thick}, ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

It might be more informative to plot the absolute error:

Plot[{(y1[x]-(y1[x]-x^5/120))/y1[x]},{x,-3,4},PlotStyle->{Thick}, ImageSize->{400,Automatic},AxesLabel->{"x","y"},PlotRange->{All,All}]

From the above we can decide if the solution is accurate enough for what we want to do. If you want more accuracy, choose a larger n right at the beginning (keep more terms in the assumed solution).

Mathematica  to Assist Hand Calculation:  (best, patterns are recognized, closed form of solution is found!)

Here is the differential equation we solved by hand in class

Clear[y,a,x,n,ygeneral,c1,c2]
L[y_]:=(1-x)*D[y,{x,2}]+x*D[y,{x,1}]-y
L[y[x]]==0

Let's use a hybrid process of solution, where we calculate the recursion relation by hand, and use Mathematica to assist us along the way. This truly is the best of both worlds.

The recursion relation we worked out in class was:

(n+2)(n+1)a[n+2]-(n+1)n a[n+1]+(n-1)a[n]==0, n=0,1,2,3,...

Mathematica can easily work with recursion relations, but we need to rewrite it slightly. If we let m=n+2, the recursion relation becomes:

recursion=(n+2)(n+1)a[n+2]-(n+1)n*a[n+1]+(n-1)a[n]==0/.n->m-2

Solve[recursion,a[m]]

Mathematica can now get the terms for you easily (I've simply copy/pasted the above and made some syntax changes):

a[0]=c1;
a[1]=c2;
a[m_]:=a[m]=((3 a[-2+m]-m a[-2+m]+2 a[-1+m]-3 m a[-1+m]+(m^2) a[-1+m])/((-1+m) m))

Table[{m,a[m]},{m,0,6}]//TableForm

The problem now is how do we recognize the pattern? That's the thing we want to be able to do at this stage. Mathematica has multiplied everything out, and we typically don't do that when we are looking for patterns. As we saw in class, the pattern was easy to pick up when this problem was done by hand. Here, since the pattern is basically a factorial, the pattern can be decerned as well from the Mathematica output--for more complicated patterns, searching for the pattern by hand may be the way to go! Mathematica can by used to double check your calculations.

ygeneral[x_]:=  a[1](x)+a[0](1+ Sum[1/i!x^i,{i,2,Infinity}])
ygeneral[x]

Wow! What happened? Mathematica is pretty good at some things, and this is one of them. Let's investigate:

1+Sum[1/i!x^i,{i,2,Infinity}]

This is where Mathematica can really help you--when it comes to turning the infinite summation into a known function! You'll be surprised at some of the bizarre looking things that turn into known functions.

Now, we can collect the solution in a[0] and a[1]:

ygeneral[x]

ygeneral[x_]=Collect[ygeneral[x],{c1,c2}]

y1[x_]=Coefficient[ygeneral[x],c1]
y2[x_]=Coefficient[ygeneral[x],c2]

Check

Simplify[L[y1[x]]==0]
Simplify[L[y2[x]]==0]

We rock. Here is a plot of the two solutions that form our fundamental set of solutions.

Plot[{y1[x],y2[x]},{x,-3,4},PlotStyle->{Thick}, ImageSize->{400,Automatic},AxesLabel->{"x","y"}]

Aside: RSolve to solve Recursion Relations

Mathematica has a neat command to solve recursion relations, called RSolve. It is very similar to DSolve. Let's try it!

Clear[a,m]
RSolve[a[m]== ((3 a[-2+m]-m a[-2+m]+2 a[-1+m]-3 m a[-1+m]+(m^2) a[-1+m])/((-1+m) m)),a[m],m]

Pretty cool. But alas, since we don't know how RSolve works, this would be like using DSolve to solve a differential equation. Wonderful that it's there, we can use it to check our work, but that bastard Barry won't let us use it in this course since the point of the course is to understand how things work!