Visualizing PDE Solutions

Using separation of variables, we obtain solutions of PDEs in the form of summations that may not convey much intuitive meaning. This notebook contains some sample code that may be helpful for both obtaining solutions and presenting them in understandable ways.

Example: Vibrating String (1D Wave Equation)

In class, we derived the governing partial differential equation (Subscript[u, tt])=(c^2)(Subscript[u, xx])with boundary conditions for fixed ends, u(0,t)=0 and u(L,t)=0, and specified initial displacement, u(x,0)=f(x). (For simplicity, let's assume that the string has unit length and is displaced and released from rest so that (Subscript[u, t])(x,0)=0.)

We note that if we look for solutions of the form u(x,t)=X(x)T(t), then the PDE separates into two ODEs, X''(x)-k X(x)=0 and (Overscript[T, ¨])(t)-k (c^2) T(t)=0, and the homogeneous boundary conditions separate into X(0)=0 and X(1)=0. Let's solve the ODEs, starting with the X-equation (which has the homogeneous conditions). If we want to call the functions X and Y, then we are about the violate the rule of thumb of using names that start with lowercase letters, so let's check first to make sure that other definitions do not exist for  X and Y (and L):

?X

?Y

?L

No other information was found, so we can proceed. Start by finding solutions of the X-equation:

DSolve[{X''[x]-k*X[x]==0, X[0]==0, X[1]==0},
X,x]

If we plug in the ODE with both BCs, we get the result that X is a function that returns 0, i.e. a trivial solution. This is one of those FREQUENT cases when some human intelligence is essential in guiding the system toward the desired solution. Recall that to get non-trivial solutions, we needed to specify that the separation constant was negative; so let's rename k -> -(p^2)and let's leave off the right end conditionfor the moment.

DSolve[
{X''[x]+p^2*X[x]==0, X[0]==0},
X,x]

Now we see that X(x)~ Sin[p x] and to get non-trivial solutions we need to find values of p satisfying Sin[p ] = 0:

Solve[ Sin[p] == 0, p]

Now take the suggestion and use Reduce:

Reduce[ Sin[p]==0, p]

which is a complicated way of telling us that p=n π, where n is an integer. We can now enumerate a set of spatial solutions (eigenfunctions or mode shapes):

(Subscript[X, n_]) [x_] := Sin[ p[n]*x]

(Subscript[X, n])[x]

?p

p[n_] := n*Pi

Now solve the T-equation. (Recall that we are assuming zero initial velocity.) Let's also chose (Subscript[T, n])[0]==1and include the arbitrary constant in the summation later:

tSolution = DSolve[{(Subscript[T, n])''[t]+p[n]^2*c^2*(Subscript[T, n])[t]==0,(Subscript[T, n])[0]==1, (Subscript[T, n])'[0]==0},
(Subscript[T, n]), t]

(Subscript[T, n])[t] /.  tSolution

(Subscript[u, n_])[x_,t_]:= (Subscript[X, n])[x]*(Subscript[T, n])[t] /. Flatten[tSolution]

(Subscript[u, n])[x,t]

u[x_, t_, m_] := Sum[(Subscript[b, n])* (Subscript[u, n])[x,t], {n,1,m}]

Test the definition of u:

u[x,t,2]

Now we need to plug in for the b's the Fourier coefficients for the odd periodic extension of the initial displacement (Piecewise):

f[x_] := (Piecewise[{{x, 0≤x≤0.1}, {0.2-x, 0.1<x≤0.2}, {0, 0.2<x<0.8}, {x-0.8, 0.8≤x≤0.9}, {1-x, 0.9<x≤1.0}}])

Plot[f[x], {x,0,1}, PlotStyle -> AbsoluteThickness[3]]

coefficientRule[i_] := (Subscript[b, i]) -> 2*Integrate[f[x]*Sin[i*Pi*x],{x,0,1}]

coefficientRule[1]

coefficientRuleList = Table[ coefficientRule[i], {i,1,40}]

u[x,t,40] /. coefficientRuleList

We need to specify a numerical value for c in order to be able to plot. Let's let c=1.

soln[x_, t_] :=0.012259351126204434` Cos[c π t] Sin[π x]+0.` Cos[2 c π t] Sin[2 π x]+0.030035087543061837` Cos[3 c π t] Sin[3 π x]+0.` Cos[4 c π t] Sin[4 π x]+0.03242277876554809` Cos[5 c π t] Sin[5 π x]+0.` Cos[6 c π t] Sin[6 π x]+0.02124924823309842` Cos[7 c π t] Sin[7 π x]+0.` Cos[8 c π t] Sin[8 π x]+0.006033334951166327` Cos[9 c π t] Sin[9 π x]+0.` Cos[10 c π t] Sin[10 π x]-0.004038844058218781` Cos[11 c π t] Sin[11 π x]+0.` Cos[12 c π t] Sin[12 π x]-0.006161024635632086` Cos[13 c π t] Sin[13 π x]+0.` Cos[14 c π t] Sin[14 π x]-0.0036025309739497876` Cos[15 c π t] Sin[15 π x]+0.` Cos[16 c π t] Sin[16 π x]-0.0009353487470157665` Cos[17 c π t] Sin[17 π x]+0.` Cos[18 c π t] Sin[18 π x]-0.00003395942140222835` Cos[19 c π t] Sin[19 π x]+0.` Cos[20 c π t] Sin[20 π x]+0.000027798982145588287` Cos[21 c π t] Sin[21 π x]+0.` Cos[22 c π t] Sin[22 π x]+0.0005109939279537932` Cos[23 c π t] Sin[23 π x]+0.` Cos[24 c π t] Sin[24 π x]+0.0012969111506219236` Cos[25 c π t] Sin[25 π x]+0.` Cos[26 c π t] Sin[26 π x]+0.0014282759443372051` Cos[27 c π t] Sin[27 π x]+0.` Cos[28 c π t] Sin[28 π x]+0.00058109409161055` Cos[29 c π t] Sin[29 π x]+0.` Cos[30 c π t] Sin[30 π x]-0.0005085329147184938` Cos[31 c π t] Sin[31 π x]+0.` Cos[32 c π t] Sin[32 π x]-0.0009561186073662283` Cos[33 c π t] Sin[33 π x]+0.` Cos[34 c π t] Sin[34 π x]-0.0006616893625622059` Cos[35 c π t] Sin[35 π x]+0.` Cos[36 c π t] Sin[36 π x]-0.0001974549217586242` Cos[37 c π t] Sin[37 π x]+0.` Cos[38 c π t] Sin[38 π x]-8.060059912034473`*^-6 Cos[39 c π t] Sin[39 π x]+0.` Cos[40 c π t] Sin[40 π x]/. c-> 1

Before trying to plot, make sure your function evaluates to something numerical:

soln[0.3, 2.4]

Visualizing the solution

Plot

Plot[soln[x,0], {x,-2,2}, PlotRange -> {-0.125,0.125}, PlotStyle -> AbsoluteThickness[3]]

Options[Plot]

Plot3D

Plot3D[soln[x,t], {x,0,1}, {t,0,2},AxesLabel->{X,T,U},
SphericalRegion -> True]

Plot3D[soln[x,t], {x,0,1}, {t,0,2}, PlotPoints -> 55, AxesLabel -> {"X", "T", "u[x,t]"},ViewPoint->{0.102, 0.000, 3.441}]

ContourPlot

ContourPlot[soln[x,t], {x,0,1}, {t,0,2}]

DensityPlot

DensityPlot[soln[x,t], {x,0,1}, {t,0,2}]

DensityPlot[soln[x,t], {x,0,1}, {t,0,2}, Mesh -> False, PlotPoints -> 50]

Sequence of images: Animation and GraphicsArray

Options[Plot]

Table[Plot[soln[x,t], {x,0,1},PlotStyle -> AbsoluteThickness[3], PlotRange -> {-0.125, 0.125}, Ticks-> False], {t,0,2.0, 0.05}]

wavePlots = %;

Length[wavePlots]

Partition[wavePlots, 6];

Show[GraphicsGrid[%]]

ListAnimate[wavePlots]