80 mathtut2Open.

nb

Chapter 6: Heat flow and diffusion

■ Section 6.1: Fourier series methods for the heat equation

Example 6.2: An inhomogeneous example

Consider the following IBVP:

∂u ∂2 u
-A = 10-7 , 0 < x < 100, t > 0,
∂t ∂x 2

u (x, 0) = 0, 0 < x < 100,

u (0, t) = 0, t > 0,
u (100, t) = 0, t > 0.

The constant A has value 0.208 cm2 s.

In[544]:= ClearAll[A]
A = 0.208
Out[545]= 0.208

The solution can be written as

∞ nπx
u (x, t) = an (t) sin ( ),
n=1 100

where the coefficient an (t) satisfies the IVP

dan An2 π2
+ an = cn , t > 0,
dt 1002
an (0) = 0.

The values c1 , c2 , c3 , ... are the Fourier sine coefficients of the constant function 10-7 :

In[546]:= $Assumptions = Element[{m, n}, Integers]

ClearAll[c, n, x]
2 / 100 Integrate10^(-7) Sinn Pi x / 100, {x, 0, 100}
Out[546]= m ∈ Integers
nπ 2
Sin 2

Out[548]=
2 500 000 n π
mathtut2Open.nb 81

In[549]:= c[n_] = %
nπ 2
Sin 2

Out[549]=
2 500 000 n π
t 2 π2 (t-s)1002
We compute an (t) by the formula an (t) = ∫0 e-An cn ⅆs.

In[550]:= ClearAll[a, n, x, t]
IntegrateExp-A n^2 Pi^2 (t - s) / 100^2 c[n], {s, 0, t}
2 2 nπ 2
ⅇ-0.000205288 n t -0.000620222 + 0.000620222 ⅇ0.000205288 n t  Sin 2

Out[551]=
n3

In[552]:= a[t_, n_] = %

2 2 nπ 2
ⅇ-0.000205288 n t -0.000620222 + 0.000620222 ⅇ0.000205288 n t  Sin 2

Out[552]=
n3

Now we define the Fourier series solution:

In[553]:= ClearAll[u, x, t, M]
u[x_, t_, M_] := Suma[t, n] Sinn Pi x / 100, {n, 1, M}

We can easily look at some "snapshots" of the solution. For example, we will show the concentration distribution after 10
minutes (600 seconds). Some trial and error may be necessary to determine how many terms in the Fourier series are
required for a qualitatively correct plot. (As discussed in the text, this number decreases as t increases, due to the smooth-
ing of the solution.)

In[555]:= Plot[u[x, 600, 10], {x, 0, 100}, PlotRange → All]

0.00006

0.00005

0.00004

Out[555]= 0.00003

0.00002

0.00001

20 40 60 80 100

The wiggles at the top of the arch suggest that we did not use enough terms in the Fourier series, so we try again with more
terms. (Of course, perhaps those wiggles are really part of the solution. In that case, they will persist when we draw the
graph with more terms.)
82 mathtut2Open.nb

In[556]:= Plot[u[x, 600, 20], {x, 0, 100}, PlotRange → All]

0.00006

0.00005

0.00004

Out[556]= 0.00003

0.00002

0.00001

20 40 60 80 100

In[557]:= Plot[u[x, 600, 40], {x, 0, 100}, PlotRange → All]

0.00006

0.00005

0.00004

Out[557]= 0.00003

0.00002

0.00001

20 40 60 80 100
mathtut2Open.nb 83

In[558]:= Plot[{u[x, 600, 10], u[x, 600, 20], u[x, 600, 40]}, {x, 0, 100}, PlotRange → All]

0.00006

0.00005

0.00004

Out[558]= 0.00003

0.00002

0.00001

20 40 60 80 100

The above graphs suggest that 20 terms is enough for a qualitatively correct graph at t=600.

■ Section 6.4: Finite element methods for the heat equation

Now we show how to use the backward Euler method with the finite element method to (approximately) solve the heat
equation. Since the backward Euler method is implicit, it is necessary to solve an equation at each step. This makes it
difficult to write a general-purpose program implementing backward Euler, and we do not attempt to do so. Instead, the
command beuler (defined below) applies the algorithm to the system of ODEs
dα
M + Kα = f (t), t > 0,
dt
α (0) = α0,

which is the result of applying the finite element method to the heat equation.

In[559]:= ClearAll[M, K, f, a0, n, dt, beuler]

beuler[M_, K_, f_, a0_, n_, dt_] :=
Modulei, L, U = Tablei * dt, Null, i, 0, n;
U[[1, 2]] = a0;
L = M + dt * K;
DoUi + 1, 2 = LinearSolveL, M.Ui, 2 + dt * fUi + 1, 1, i, 1, n;
U

To solve a specific problem, we have to compute the mass matrix M, the stiffness matrix K, the load vector f, and the initial
data α0 . The techniques should by now be familiar.

Example 6.8

An 100 cm iron bar, with ρ=7.88 gcm3 , c=0.437 J/(g K), and κ=0.836 W/(cm K), is chilled to an initial temperature of 0
degrees, and then heated internally with both ends maintained at 0 degrees. The heat source is described by the following
function:
84 mathtut2Open.nb

In[561]:= ClearAll[F, x, t]
F[x_, t_] = 10^(-8) t x (100 - x)^2
t (100 - x)2 x
Out[562]=
100 000 000

The temperature distribution is the solution of the IBVP

∂u ∂2 u
ρc -κ = F (x, t), 0 < x < 100, t > 0,
∂t ∂x2
u (x, 0) = 0, 0 < x < 100,
u (0, t) = 0, t > 0,
u (100, t) = 0, t > 0.

We begin by defining the mesh and the constants:

In[563]:= ClearAll[n, h, k, p, c]
n = 100
h = 100. / n
k = 0.836
p = 7.88
c = 0.437
Out[564]= 100
Out[565]= 1.
Out[566]= 0.836
Out[567]= 7.88
Out[568]= 0.437

Next we define the stiffness matrix (for this constant coefficient problem, there is no need to perform any integrations---we
already know the entries in the stiffness matrix):

In[569]:= ClearAllK, i, j
K = Table0.0, i, 1, n - 1, j, 1, n - 1;
DoKi, i = 2 k / h, i, 1, n - 1
DoKi, i + 1 = -k / h;
Ki + 1, i = Ki, i + 1, i, 1, n - 2

Similarly, we know the entries of the mass matrix:

In[573]:= ClearAllM, i, j
M = Table0.0, i, 1, n - 1, j, 1, n - 1;
DoMi, i = 2 h p c / 3, i, 1, n - 1
DoMi, i + 1 = h p c / 6;
Mi + 1, i = Mi, i + 1, i, 1, n - 2

The load vector f = f (t) (which is a vector-valued function) can be computed with the Table command:
mathtut2Open.nb 85

In[577]:= ClearAllphi1, phi2, x, i

phi1x_, i_ = x - i - 1 * h  h
phi2x_, i_ = -x - i + 1 * h  h
Out[578]= 1. (-1. (-1 + i) + x)
Out[579]= 1. (1. (1 + i) - x)

In[580]:= ClearAll[f, t]
f[t_] = TableIntegrateF[x, t] × phi1x, i, x, i * h - h, i * h +
IntegrateF[x, t] × phi2x, i, x, i * h, i * h + h, i, 1, n - 1;

Finally, before invoking beuler, we need the initial vector:

In[582]:= ClearAll[a0]
a0 = Table0.0, i, 1, n - 1;

Now we choose the time step and the number of time steps, and invoke the backward Euler method:

In[584]:= ClearAll[dt, steps, U]

dt = 2.
steps = 180. / dt
U = beuler[M, K, f, a0, steps, dt];
Out[585]= 2.
Out[586]= 90.

Finally, we set up a table of data and display the result at time t=10. Note that U[[steps+1,1]] is the final time (in this case,
t = 180) and U[[steps+1,2]] is the vector of nodal values of the piecewise linear function approximating u(x, t) for this
value of t.
86 mathtut2Open.nb

In[588]:= ClearAllX, i, u, dat

X = Tablei * h, i, 0, n;
u = Flatten[{0, U[[steps + 1, 2]], 0}];
dat = Transpose[{X, u}];
ListPlot[dat]

4
Out[592]=
3

20 40 60 80 100

Chapter 8: First-order PDEs and the Method of Characteristics

■ Section 8.1: The simplest PDE and the method of characteristics

When solving PDEs in two variables, it is sometimes desirable to graph the solution as a function of two variables (that is,
as a surface), rather than plotting snapshots of the solution. This is particularly appropriate when neither variable is time.

The Plot3D function plots a function of two variables:

In[593]:= ClearAll[u, x, y]
u[x_, y_] = 1 / (1 + (x + y / 2)^2)
1
Out[594]=
y 2
1 + x + 2 
mathtut2Open.nb 87

In[595]:= Plot3D[u[x, y], {x, -5, 5}, {y, 0, 10}]

Out[595]=

If you would like to see the surface from a different angle, you can click on the figure and rotate it by moving the mouse
(i.e. put the pointer on the figure, hold down the left mouse button, and move the mouse).

■ Section 8.2: First-order quasi-linear PDEs

The purpose of the method of characteristics is to reduce a PDE to a family of ODEs. Mathematica has a function called
DSolve that will solve many ordinary differential equations symbolically. We will illustrate the use of DSolve in the
context of Example 8.6 from the text.

The method of characteristics reduces the PDE in Example 8.6 to the IVP
dv 1
= v2 , v (0) = .
dt 1 + s2

DSolve will solve this problem as follows:

In[596]:= ClearAll[v, t, s]
DSolve[{v'[t] ⩵ v[t]^2, v[0] == 1 / (1 + s^2)}, v[t], t]
1
Out[597]= v[t] → 
1 + s2 - t

If no initial condition is given, DSolve will return the general solution of the ODE:
88 mathtut2Open.nb

In[598]:= DSolve[v'[t] ⩵ v[t]^2, v[t], t]

1
Out[598]= v[t] → 
-t - C[1]

DSolve will solve a system of ODEs (when possible). Here is the system of ODEs from Example 8.7:
dx
= v, x (0) = s,
dt
dy
= y, y (0) = 1,
dt
dv
= x, v (0) = 2 s.
dt

The solution is given as follows:

In[599]:= ClearAll[x, y, v, t, s, sol]
sol = DSolve[{x'[t] == v[t], y'[t] == y[t], v'[t] == x[t],
x[0] == s, y[0] == 1, v[0] == 2 s}, {x[t], y[t], v[t]}, t]
1 1
Out[600]= v[t] → ⅇ-t 1 + 3 ⅇ2 t  s, x[t] → ⅇ-t -1 + 3 ⅇ2 t  s, y[t] → ⅇt 
2 2

The solutions can be defined for further use as follows:

In[601]:= x[t_] = (x[t] /. sol[[1]])

y[t_] = (y[t] /. sol[[1]])
v[t_] = (v[t] /. sol[[1]])
1
Out[601]= ⅇ-t -1 + 3 ⅇ2 t  s
2
Out[602]= ⅇt
1
Out[603]= ⅇ-t 1 + 3 ⅇ2 t  s
2

Chapter 11: Problems in multiple spatial dimensions

■ Section 11.2: Fourier series on a rectangular domain

Fourier series calculations on a rectangular domain proceed in almost the same fashion as in one-dimensional problems.
The key difference is that we must compute double integrals and double sums in place of single integrals and single sums.
Fortunately, Mathematica makes this easy. We do not need to learn any new commands, since a double integral over a
rectangle can be computed as an iterated integral.

As an example, we compute the Fourier double sine series of the following function f on the unit square: