Demo #6
Linear ODE's with Variable Coefficients

In this Demo, we will have a look at Mathematica's ability to solve homogeneous, linear ODE's with non-constant coefficients. In the most general case, a second order, linear ODE with variable coefficients, i.e.

y''[x] + p[x] y'[x] + q[x] y[x] == 0

has two linearly independent solutions. However, these solutions are often not easily expressed in terms of simple functions with which you are familiar. These solutions can always be developed in power series, but the easiest way to learn about a given ODE is to simply try to "look it up" using Mathematica's DSolve function!

To start, let's feed DSolve with a well-studied equation, namely Airy's differential equation:

airy = y''[x] - x y[x] == 0

DSolve[ airy , y[x],x]

So, Mathematica is educated about Airy's equation!  It says that the general solution is a linear combination of two functions, AiryAi[x] and AiryBi[x], that it seems to know about. (Note that these functions are usually denoted Ai[x] and Bi[x]).

Let's now plot the two Airy functions so that we understand their behavior:

Plot[ AiryAi[x], {x,-10,3}]

So, the Ai Airy function decays rapidly at positive x, but oscillates and decays more slowly at negative x.

Plot[ AiryBi[x],{x, -10,2}]

This other solution of the Airy equation evidently blows up for positive x, but oscillates out of
phase with Ai[x] for negative x.

To verify that the two functions are indeed linearly independent (LI), let's form the Wronskian:

w[x_]:=  AiryAi[x] AiryBiPrime[x] - AiryBi[x] AiryAiPrime[x]

We now plot the Wronskian to verify that it is nonvanishing and confirm linear independence:

Plot[ w[x],{x,-10,2},PlotRange -> {0,1}]

w[0.0]

Indeed, it is a positive constant, approximately 0.31831, apparently over the whole real x axis! LI is thus demonstrated!

Let's now see what Mathematica knows about a second very important linear, homogeneous ODE known as Bessel's equation:

beseq = x^2 y''[x] + x y'[x] + (x^2 - n^2 ) y[x] == 0

Here "n" is just a constant parameter in the equation.

DSolve[ beseq ,y[x],x]

Indeed, Mathematica knows about Bessel equation! DSolve tells us that the General solution is a linear combination
of two "Bessel" functions with + or - values of the parameter n, BesselJ[n,x]. Note that Mathematica writes Sqrt[x^2] for x!

In the special case of n an integer, the two solutions are not linearly independent, so we have to consider it separately. For example, n=0:

\(\(\(
\)\(beseq2\  = \ beseq\  /. \ n -> 0\)\)\)

\(\(\(
\)\(DSolve[\ beseq2\ , y[x], x]\)\)\)

Thus, we see two different linearly independent functions for n=0: BesselJ[0,x] and BesselY[0,x].More generally,
BesselJ[n,x] and BesselY[n,x] are linearly independent solutions of Bessel's equation regardless of whether n is
an integer or not. Thus, it is safer to use BesselY[n,x] rather than BesselJ[-n,x] as the second linearly independent
solution of Bessel's equation.

Let's now plot these up for the special case of n=0 to get a first sense for how Bessel functions behave:

Plot[BesselJ[0,x] ,{x,0,10}]

The "J" function thus oscillates and decays.

Plot[BesselY[0,x] ,{x,0,10}]

So, the "Y" function diverges to negative values at x=0 and oscillates and decays for large x.

Clearly, we could continue this exercise of running through all the "classical" ODE's that Mathematica knows about. Instead, however, let's now throw DSolve a non-classical one:

tough = y''[x] + Tanh[x] y'[x] + x y[x] == 0

DSolve[ tough ,y[x],x]

Thus, when Mathematica is asked to analytically solve a "tough" ODE, i.e. one whose solutions can evidently
not be expressed in terms of simple functions, it simply gives up. Besides approximate analytical methods, or power
series methods, the only choice in such a case is to seek a numerical solution (once initial conditions or boundary
conditions have been specified and any inhomogeneous term restored).