Automating computations for PDEs

with Mathematica

We’ve seen lots of cool, computer generated illustrations of solutions of PDEs this semester. With the test just around the corner, I’d like to show you a tool that will ease your algebraic workload. When it comes to computer algebra, there’s definitely nothing more powerful than Wolfram Research’s Mathematica.

Assumptions

The following cell will set the assumption that m and n are positive integers, affecting the results of algebraic integration. For example, cos(n π) will be simplified to (((-1))^n).

$Assumptions={m∈Integers, m>0,n∈Integers, n>0};

Computing a full Fourier series

Let’s try to find the full Fourier series of f(x) = (e^x) over the interval [-L,L], where L=1. I guess it should look like

((Subscript[a, 0])/2)+(∑_(n=1)^∞)((Subscript[a, n])cos(n π x/L)+(Subscript[b, n])sin(n π x/L)),

where

(Subscript[b, n]) = (2/L)(∫_0^L)f(x)sin(n π x/L)dx   and   (Subscript[b, n]) = (2/L)(∫_0^L)f(x)sin(n π x/L)dx.

Here’s how to automate this process with Mathematica:

L = 1;
f[x_]= Exp[x];
a[0] = Integrate[f[x],{x,-L,L}]/L
a[n_] = Integrate[f[x]*Cos[n*Pi*x/L],{x,-L,L}]/L
b[n_] = Integrate[f[x]*Sin[n*Pi*x/L],{x,-L,L}]/L

Let’s take a look:

F[x_] = a[0]/2 + Sum[a[n]*Cos[n*Pi*x/L]+b[n]*Sin[n*Pi*x/L],{n,1,12}];
Plot[{f[x],F[x]},{x,-L,L}]

That looks pretty good!

Computing a Fourier sine series

Now let’s try to find the Fourier sine series of f(x) = 1-x+(x^2)-(x^3) over the interval [0,L], where L=2. I guess it should look like

(∑_(n=1)^∞)(Subscript[b, n])sin(n π x/L),

where

(Subscript[b, n]) = (2/L)(∫_0^L)f(x)sin(n π x/L)dx.

Here’s how to automate this with Mathematica:

L = 2;
f[x_]= 1-x+x^2-x^3;
b[n_] = 2Integrate[f[x]*Sin[n*Pi*x/L],{x,0,L}]/L

Let’s take a look:

F[x_] = Sum[b[n]*Sin[n*Pi*x/L],{n,1,42}];
Plot[{f[x],F[x]},{x,0,L}]

That looks good, too!

Application

Of course, you should understand that this series arises in the solution to the heat problem

(Subscript[u, t])=(Subscript[u, x x]), u(0,t)=u(2,t)=0, u(x,0)=1-x+(x^2)-(x^3).

In fact the solution should look like

u(x,t) = (∑_(n=1)^∞)(Subscript[b, n])(e^(-(n^2)(π^2)t/(L^2)))sin(n π x/L).

Let’s take a look at that too:

Clear[u]

u[x_,t_] = Sum[
b[n]Exp[-n^2*Pi^2*t/L^2]*Sin[n*Pi*x/L],
{n,1,42}];
u[x_,0.0] = f[x];
Manipulate[Quiet@Plot[u[x,t],{x,0,L},
PlotRange -> {-6,2}],{t,0.0,2.0}]

Computing a Fourier cosine series

Next, let’s try to find the Fourier cosine series of

f(x) = (Piecewise[{{x-1, if x<3}, {x-3, if x>2}}])

over the interval [0,L], where L=4. I guess it should look like

((Subscript[a, 0])/2)+(∑_(n=1)^∞)(Subscript[a, n])cos(n π x/L),

where

(Subscript[a, n]) = (2/L)(∫_0^L)f(x)cos(n π x/L)dx.

Here’s how to automate this with Mathematica:

L = 4;
f[x_]= If[x<2,x-1,x-3];
b[0]  = 2Integrate[f[x],{x,0,L}]/L
b[n_] = 2Integrate[f[x]*Cos[n*Pi*x/L],{x,0,L}]/L //Simplify

Let’s take a look:

F[x_] = Sum[b[0]/2+b[n]*Cos[n*Pi*x/L],{n,1,42}];
Plot[{f[x],F[x]},{x,0,L}]

That looks good, too!

Application

Of course, you should understand that this series arises in the solution to the heat problem

(Subscript[u, t])=(Subscript[u, x x]), (Subscript[u, x])(0,t)=(Subscript[u, x])(4,t)=0, u(x,0)= (Piecewise[{{x-1, if x<3}, {x-3, if x>2}}])

In fact the solution should look like

u(x,t) =((Subscript[a, 0])/2)+ (∑_(n=1)^∞)(Subscript[a, n])(e^(-(n^2)(π^2)t/(L^2)))cos(n π x/L).

Let’s take a look at that too:

Clear[u];
u[x_,t_] = b[0]/2+Sum[
b[n]Exp[-n^2*Pi^2*t/L^2]*Cos[n*Pi*x/L],
{n,1,42}];
u[x_,0.0] = f[x];
Manipulate[Quiet@Plot[u[x,t],{x,0,L},
PlotRange -> {-1,1}],{t,0.0,3.0}]

Sturm-Liouville problems

Example 1

Consider the Sturm-Liouville problem

-y''=λ y.

with boundary conditions y(0)=y'(1)=0. This is about the simplest modification of our most basic eigenvalue problem and wee worked through a problem like this by hand last time.

The general solution

To find the eigenvalues and eigenfunctions, the first step is to find the general solution.

Clear[y];
y[x_]=DSolveValue[-y''[x]==λ*y[x],y[x],x]

Finding λ

We use the boundary values to find possible values for λ.  Now, we need y'(0)=0, which implies that (Subscript[c, 1])=0.

y[0]==0

Next, we need y'(1)=0, already knowing that (Subscript[c, 1])=0.

y'[1]==0 /. C[1]->0

The only way this happens is if (Sqrt[λ])=(2n-1)π/2 for some n.  Thus our solutions to the SLP have the form

y(x)=(Subscript[c, n])sin((2n-1)π x/2)

Using the determinant

If y(x)=(Subscript[c, 1])cos((Sqrt[λ])x)+(Subscript[c, 2])sin((Sqrt[λ])x), then

y'(x)=-(Subscript[c, 1])(Sqrt[λ])sin((Sqrt[λ])x)+(Subscript[c, 2])(Sqrt[λ])cos((Sqrt[λ])x).

Thus, the boundary conditions y(0)=y'(1)=0 can be written

(Subscript[c, 1])+0=0
-(Sqrt[λ])sin((Sqrt[λ]))(Subscript[c, 1])+(Subscript[c, 2])(Sqrt[λ])cos((Sqrt[λ]))=0

or, as the matrix equation

(({{1, 0}, {-(Sqrt[λ])sin((Sqrt[λ])), (Sqrt[λ])cos((Sqrt[λ]))}}))(({{(Subscript[c, 1])}, {(Subscript[c, 2])}}))=(({{0}, {0}})).

This has a non-trivial solution iff the determinant of the matrix is zero.  The determinant is just (Sqrt[λ])cos((Sqrt[λ])), which puts us where we were before.

Checking orthogonality

The  following computation indicates that these functions are orthogonal for m≠n.

y[n_,x_] = Sin[(2n-1)Pi*x/2];
Integrate[y[m,x]*y[n,x],{x,0,1}]

And the norm-squared of each function is 1/2.

Integrate[((y[n,x])^2),{x,0,1}]

Note that we’ll have to account for the 1/2 in approximations.

A PDE

Let’s put this in the context of a PDE problem.  Suppose we’d like to solve (Subscript[u, t])=(Subscript[u, x x]), subject to u(0,t)=u'(1,t)=0 and u(x,0)=x(2-x).

c[n_] = 2Integrate[x(2-x)y[n,x],{x,0,1}]

Manipulate[
Plot[{x(2-x),Sum[c[n]*Exp[-t((((2n-1)Pi/2))^2)]y[n,x],{n,1,numTerms}]},
{x,0,1},PlotRange -> {0,1.2}],
{t,0,1}, {numTerms,1,10,1}]

Example II

Consider the Sturm-Liouville problem

-y''=λ y.

with boundary conditions y(0)=y(1)+y'(1)=0 - a mild variation of a text book problem.

The general solution

We certainly have the same general solution.

Clear[y];
y[x_]=DSolveValue[-y''[x]==λ*y[x],y[x],x]

Finding λ

We again have y(0)=(Subscript[c, 1])+0=0, which implies (Subscript[c, 1])=0.  
Next, we have

y(1)+y'(1)=(Subscript[c, 1])cos((Sqrt[λ]))+(Subscript[c, 2])sin((Sqrt[λ]))-(Sqrt[λ])sin((Sqrt[λ]))(Subscript[c, 1])+(Subscript[c, 2])(Sqrt[λ])cos((Sqrt[λ]))=0.

Since (Subscript[c, 1])=0, this reduces to

(Subscript[c, 2])sin((Sqrt[λ]))+(Subscript[c, 2])(Sqrt[λ])cos((Sqrt[λ]))=0  or  (Sqrt[λ])=-tan((Sqrt[λ])).

Thus, our solutions have the specific form

y(x)=c sin((Sqrt[(Subscript[λ, n])])x)

where (Subscript[λ, n]) is one of the solutions to (Sqrt[λ])=-tan((Sqrt[λ])).  Note that the Sturm-Liouville theory guarantees the existence of infinitely many such solutions.  This equation is not solvable in closed form.  We can numerically approximate some solutions within a particular range, though.

intersections =λ/. NSolve[{-Tan[(Sqrt[λ])]==(Sqrt[λ]),0<λ<220},λ]

We can even understand the asymptotic behavior.  The following plot of y=(Sqrt[λ]) together with y=-tan((Sqrt[λ])) indicates why (Sqrt[λ]) asymptotically approaches the roots of cos(x).

tanPart=Plot[-Tan[(Sqrt[λ])],{λ,0,220},
Exclusions -> Cos[(Sqrt[λ])]==0, ExclusionsStyle -> Dashed,PlotRange -> {-5,15}];
sqrtPart=Plot[Sqrt[λ],{λ,0,220},
PlotStyle -> Directive[ColorData[1,2]]];
Show[{tanPart,sqrtPart, Graphics[{PointSize[Medium],
Point[{#,(Sqrt[#])}]& /@ intersections}]}]

Checking orthogonality

Again, our solutions have the specific form

y(x)=c sin((Sqrt[(Subscript[λ, n])])x)

where (Subscript[λ, n]) is one of the solutions to (Sqrt[λ])=-tan((Sqrt[λ])); the first few such values are

λns =λ/. NSolve[{-Tan[(Sqrt[λ])]==(Sqrt[λ]),0<λ<1000},λ]

We can run a numerical check for orthogonality:

Clear[y];
y[n_,x_] := Sin[(Sqrt[λns[[n]]])x];
Table[NIntegrate[
y[m,x]*y[n,x],{x,0,1}],
{n,1,Length[λns]},{m,1,Length[λns]}]//Chop//MatrixForm//Quiet

Here’s what these things look like:

Plot[Evaluate[Table[y[n,x]/3,{n,1,5}]],{x,0,1},
AspectRatio->Automatic]

Note that, on the right, we should have y(1)+y'(1)=0.

Approximating a function

Let’s use this to approximate the function f(x)=((x(1-x))^2).  Note that this function satisfies the boundary conditions y(0)=0, y(1)+y'(1)=0, which is kinda fun.

Clear[n];
c[n_] :=c[n]= ((NIntegrate[x(((1-x))^2)y[n,x],{x,0,1}])/(NIntegrate[((y[n,x])^2),{x,0,1}]));

Manipulate[
Plot[{x(((1-x))^2),Sum[c[n]*y[n,x],{n,1,numTerms}]},{x,0,1},
PlotRange -> {0,0.16}],
{numTerms,1,8,1}]

Example III - Airy functions

We now consider the Sturm-Liouville problem

-y''+x y=λ y.

with boundary conditions y(0)=y(1)=0.  This will lead to some odd functions called the Airy functions.

The general solution

To find the eigenvalues and eigenfunctions, the first step is to ask Mathematica to find the general solution for us.

Remove[y];
y[x_]=DSolveValue[-y''[x]+x*y[x]==λ*y[x],y[x],x]

The solution is expressed in terms of strange entities called Airy functions. These were discovered in the late 1830s by George Airy in is study of optics. They have found their way into  many applications since then.

The Airy functions

The Airy functions are most simply  described as the solutions of x''(t)=t x(t):

DSolveValue[x''[t]==t*x[t],x[t],t]

Note that this differential equation looks a lot like the one for the trig functions, but we have the variable coefficient x, rather than a constant. Thus, the solutions should look somewhat trigonometric for x<0 but exponential for x>0.

Plot[{AiryAi[x],AiryBi[x]},{x,-12,3}]

Finding λ

The boundary conditions state that we would like the solution to be zero at both endpoints.  Of course, we're not interested in the trivial solution, so the question is - for which values of λ are there constants, not both zero, so that the boundary conditions are simultaneously satisfied. We can express this in our matrix formulation as follows.

(({{Ai(-λ), Bi(-λ)}, {Ai(1-λ), Bi(1-λ)}}))(({{(Subscript[c, 1])}, {(Subscript[c, 2])}}))=(({{0}, {0}})).

Of course, this system has a non-trivial solution precisely when the matrix

M(λ)=(({{Ai(-λ), Bi(-λ)}, {Ai(1-λ), Bi(1-λ)}}))

is singular.  Again, there is a simple test for singularity, namely the determinant must be zero.  Thus, let's define the matrix and the determinant as functions of λ and plot the value of the determinant in terms of λ.

M[λ_] = (({{AiryAi[-λ], AiryBi[-λ]}, {AiryAi[1-λ], AiryBi[1-λ]}}));
det[λ_] = Det[M[λ]]
Plot[det[λ],{λ,0,400}]

The roots are eigenvalues of the problem.  Here are all the eigenvalues less than 2000.

eigenvalues = λ/.NSolve[det[λ]==0&&0<λ<2000,λ,
WorkingPrecision->30];
N[eigenvalues]

Now, for each of the eigenvalues, there corresponds a choice of (Subscript[c, 1]) and (Subscript[c, 2]) such that (Subscript[c, 1])Ai(-λ)+(Subscript[c, 2])Bi(-λ)=0 and (Subscript[c, 1])Ai(1-λ)+(Subscript[c, 2])Bi(1-λ)=0.  Put another way, the vector ⟨(Subscript[c, 1]),(Subscript[c, 2])⟩ should map to zero under multiplication by the matrix M(λ).  Such a vector is said to lie in the null-space of the matrix. The function (Subscript[y, λ])(x)=(Subscript[c, 1])Ai(x-λ)+(Subscript[c, 2])Bi(x-λ) is then the corresponding eigenfunction.  The NullSpace command is a bit delicate numerically, which is why we requested high precision eigenvalues.  Thus, we now compute the null spaces, toss away the extra precision, define the eigenfunctions, normalize, and visualize.

nullVectors = N[First/@Table[NullSpace[M[λ]],{λ,eigenvalues}]];
eigenvalues=N[eigenvalues];
(Subscript[y, k_])[x_]:=nullVectors[[k]].{AiryAi[x-eigenvalues[[k]]],
AiryBi[x-eigenvalues[[k]]]};
Do[normSquared[k]=NIntegrate[(((Subscript[y, k])[x])^2),{x,0,1}],
{k,1,Length[eigenvalues]}]
Plot[Evaluate@Table[(Subscript[y, k])[x]/Sqrt[normSquared[k]],{k,1,8}],
{x,0,1}]

If this all worked, then the following computation should yield an identity matrix.

Table[NIntegrate[(((Subscript[y, k])[x])/(Sqrt[normSquared[k]]))(((Subscript[y, j])[x])/(Sqrt[normSquared[j]])),
{x,0,1}],{k,1,8},{j,1,8}]//Quiet//Chop//MatrixForm

Approximating a function

OK, let's try to approximate a function.

f[x_] = (1-Abs[1-2x])/2;
Do[(Subscript[c, k])=NIntegrate[((f[x](Subscript[y, k])[x])/(normSquared[k])),{x,0,1}],
{k,1,14}];
m=Length[eigenvalues];
approx[x_] = (∑_(k=1)^m)(Subscript[c, k])(Subscript[y, k])[x];
Plot[{f[x],approx[x]},{x,0,1}]

Necessary numerics

You should understand that there are integrals that can’t be evaluated in closed form.

Integrate[Cos[Cos[x]],x]

Thus, there are Fourier series whose coefficients can’t be evaluated in closed form.

L = Pi;
f[x_]= Cos[Cos[x]];
b[n_] = 2Integrate[f[x]*Sin[n*Pi*x/L],{x,0,L}]/L

Nonetheless these integrals can still be evaluated numerically:

Clear[b];
b[n_]:= b[n]= 2 NIntegrate[f[x]*Sin[n*Pi*x/L], {x, 0, L}]/L;
Table[b[n],{n,1.,10}]//Quiet//Chop

and we can use these decimal numbers to construct an approximate Fourier series:

Clear[F];
F[x_Real] := Sum[b[n]*Sin[n*Pi*x/L],{n,1,42}];
Plot[{f[x],F[x]},{x,0,L}]//Quiet