Chapter 3
Second-Order Partial Differential Equations

3.1 Equations with Constant Coefficients

The operations carried out in Equations (3.24)-(3.29), which involve a straightforward but somewhat lengthy sequence of algebraic steps, can be performed easily in Mathematica.  We first define the function u as a function of xi[x,y] and eta[x,y]:

f=u[xi[x,y],eta[x,y]];

We now carry out the operations required to obtain (3.29).  The original differential equation is transformed into the new variables, with the derivatives carried out using the chain rule automatically.

Eqn=Expand[a D[f,{x,2}] + b D[f,x,y] +c D[f,{y,2}]];

The coefficients A, B and C of the second-order derivatives in the transformed variables are then obtained using the command Coefficient.  For example,

A=Coefficient[Eqn,
     Derivative[2,0][u][xi[x,y],eta[x,y]]]

Similarly, for obtaining B and C we use the commands

B=Coefficient[Eqn,
     Derivative[1,1][u][xi[x,y],eta[x,y]]]

Unprotect[C]

C=Coefficient[Eqn,
     Derivative[0,2][u][xi[x,y],eta[x,y]]]

Protect[C]

Finally, the quantity B^2-4AC is constructed, and the result in (3.28) is obtained:

Factor[B^2-4 A C]

3.2 Boundary Conditions of Second-Order Equations

The appearance of the condition Delta≠0 in the evaluation of the higher-order terms in (3.32) can be demonstrated explicitly using Mathematica.  The kth-order terms in (3.32) involve the k+1 distinct partial derivatives of u with respect to x and y, since the assumed analyticity of u implies that the order of differentiation is
unimportant.  The function u and the second-order terms in the differential
equation are first defined:

x=x0[s];
y=y0[s];
f=u[x,y];

Eqn=a[s] D[f,{x,2}] + b[s] D[f,x,y] + c[s] D[f,{y,2}];

The k+1 equations for the kth-order derivatives of u are obtained from differentiating the derivatives of order k-1 of f with respect to s, which yields k equations, and differentiating (3.4) k-2 times.  The resulting matrix of coefficients, M[k] is constructed as a Table:

M[k_]:=Table[If[i<k+1,
     Coefficient[D[f,s,{x,k-i},{y,i-1}],
     Derivative[k+1-j,j-1][u][x,y]],
     Coefficient[D[Eqn,{x,k-2}],
     Derivative[k+1-j,j-1][u][x,y]]],
     {i,1,k+1},{j,1,k+1}]

The argument k has been defined to correspond to the order of the partial derivatives that are required.  The coefficients contained in M can be displayed as a matrix using MatrixForm and the determinant of M can be evaluated using Det.  Examples are shown for the second and fourth-order terms.  In both cases, the determinant is proportional to (3.37).  This can be seen from the structure of M as a function of k and the evaluation of the determinant using minors.   In this way, the condition for obtaining solutions of kth-order derivatives in the expansion (3.32) can be extended to all orders.

MatrixForm[M[2]]

Simplify[Det[M[2]]]

MatrixForm[M[4]]

Simplify[Det[M[4]]]

3.4 Hyperbolic Equations

Example 3.1

The influence of the two terms on the form of the solution in (3.64) can be seen by
considering two cases:

            f[x]=Exp[-x^2],     g[x]=0                                                        (3.68a)
               
            f[x]=0,      g[x]=Exp[-x^2]                                                       (3.68b)
               
The profiles of these two solutions with v=1 that have been generated with Plot are shown in Figure 3.2.

For the initial conditions in (3.68a), the initial form is seen to break up into two travelling waves moving in opposite directions with unit speed, i.e. along the characteristics.  The solution for the initial conditions (3.68b) evolves by the value at each point rising to the value 0.886 with a 'wavefront' that also moves with unit speed.

A direct correspondence between the time development of these solutions and the
characteristics can be more easily seen by using ContourPlot, as shown in
Figure 3.3.  For both initial conditions in (3.68), the contours of the solutions
eventually coincide with the characteristics.

Figures 3.2 and 3.3

The left panel of Figure 3.2 is obtained by first constructing the solution using the d'Alembert form and then using Plot to display the solutions for t=0, t=1, t=2, and t=3.

Clear[x]

u[x_,t_]:=(1/2)(Exp[-(x + t)^2]+Exp[-(x - t)^2])

Do[Plot[u[x,k],{x,-5,5},AxesLabel->{x,u},
     PlotRange->{0,1}],{k,0,3}]

The top panel of Figure 3.3 is obtained with ContourPlot:

ContourPlot[u[x,t],{x,-5,5},{t,0,6},PlotRange->All,
     Contours->20,ContourShading->False,
     PlotPoints->60]

The right panel of Figure 3.2 is generated in a similar fashion to the left panel:

u[x_,t_]=(1/2)Integrate[Exp[-s^2],{s,x-t,x+t}]

Do[Plot[u[x,k],{x,-5,5},AxesLabel->{x,u},
     PlotRange->{0,1}],{k,0,3}]

ContourPlot[u[x,t],{x,-5,5},{t,0,6},PlotRange->All,
     Contours->20,ContourShading->False,
     PlotPoints->60]

3.5 Parabolic Equations

Example 3.2

Mathematica can be used to investigate solutions of the diffusion equation for a wide range of initial conditions by using the symbolic capabilities to carry out the differentiations and the summation in the expansion in (3.79).  This is not the most efficient way of solving the diffusion equation, but it does provide an instance where the power series method provides an explicit exact solution of this equation by using only elementary operations.

Consider the initial-value problem with f[x]=Exp[-x^2].  The solution (3.79) is obtained by first defining the initial function f, performing the required derivatives of f and assembling the terms in (3.79) as a partial sum using Sum:

f[x_]:=Exp[-x^2]

a[x_,k_]:=D[f[s],{s,2k}]/.s->x

u[x_,y_,n_]:=Sum[(1/k!)a[x,k]y^k,{k,0,n}]

The profiles of this solution for y=0, 0.05, 0.10 and 0.15 obtained by retaining the required number of terms to achieve an accuracy of 1% are shown in Figure 3.4.  With increasing time, the solution is seen to decrease in value near the origin and the profile broadens.  This is a characteristic feature of the diffusion equation and is related to the behavior of the solution in (3.81) as a function of the quantity k (Problem 19).

Figure 3.4

This example illustrates just how many terms are required to obtain good convergence in this solution even for short times, i.e. t<1.  In fact, the diagram in Figure 3.4 was not actually generated with the power series method.  To have done so would have required too many terms in the expansion in (3.79) to obtain convergence over the time span shown.  Instead, the solution was obtained from the Green's function of the diffusion equation (Chapter 7).  For t=0, 0.05, and 0.10, and 0.15, there is no problem in obtaining the solution to 1% accuracy, but the diagram is not as revealing as that in Figure 3.4.

The plots below show just how slow the convergence of the power series solution is.  The solutions below have been constructed by retaining the first 5, the first 10, the first 20, and the first 50 terms in the expansion in (3.79). For each series below, the solution is shown for t=0.05, 0.1, 0.15, 0.20 and 0.25.  Even after the inclusion of 50 terms in the series, the solution clearly has not converged at t=0.25.

u5[x_,y_]=u[x,y,5];

Do[Plot[u5[x,j 0.05],{x,-3,3},PlotRange->{0,1},
     AxesLabel->{x,u}],{j,1,5}]

u10[x_,y_]=u[x,y,10];

Do[Plot[u10[x,j 0.05],{x,-3,3},PlotRange->{0,1},
     AxesLabel->{x,u}],{j,1,5}]

u20[x_,y_]=u[x,y,20];

Do[Plot[u20[x,j 0.05],{x,-3,3},PlotRange->{0,1},
     AxesLabel->{x,u}],{j,1,5}]

u50[x_,y_]=u[x,y,50];

Do[Plot[u50[x,j 0.05],{x,-3,3},PlotRange->{0,1},
     AxesLabel->{x,u}],{j,1,5}]

3.6 Elliptic Equations

Example 3.5/Figure 3.5

The top panel in Figure 3.5 is obtained by applying ContourPlot to the d'Alembertain form of the solution given in (3.88) with the f and g given in (3.89).  Thus,

f[x_]:=1/(1+x^2)

u[x_,y_]=Together[(1/2)(f[x+I y]+f[x-I y])]

ContourPlot[u[x,y],{x,-2,2},{y,-2,2},PlotPoints->80,
     Contours->20,ContourShading->False]

The top panel is generated with Plot3D:

Plot3D[u[x,y],{x,-2,2},{y,-2,2},PlotRange->{-2,2},
     Shading->False,PlotPoints->60,ClipFill->None,
     BoxRatios->{1,1,1}]

Example 3.6/Figure 3.6

The diagram produced in the book is not correct, but has been corrected in the second printing. The correct panels are shown below.  The solutions in (3.98) and (3.99) are combined into a single expression by using an If statement.  ContourPlot and Plot3D are then used to obtain the two panels of the correct figure.

r[x_,y_]:=Sqrt[x^2+y^2];

u[x_,y_]:=If[r[x,y]<1,x,x/(r[x,y])^(2)]

ContourPlot[u[x,y],{x,-2,2},{y,-2,2},PlotPoints->80,
     Contours->40,ContourShading->False]

Plot3D[u[x,y],{x,-2,2},{y,-2,2},PlotRange->{-1,1},
     Shading->False,PlotPoints->60,
     ViewPoint->{-1.3,2.4,1.5},BoxRatios->{1,1,1}]