Preface
This section presents basic material about parabolic equations and their main applications to heat or diffusion equations. Our presentation mostly devoted to examples.
Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VI of the course APMA0340
Introduction to Linear Algebra with Mathematica
Glossary
Parabolic equations
Most partial differential equations are of three basic types: elliptic, hyperbolic, and parabolic. In this section, we discuss the only one type of partial differential equations (PDEs for short)---parabolic equations and its most important applications: heat transfer equations and diffussion equations. Parabolic equations are often called diffusion equations because they describe the diffusion and convection of some substance (such as heat). The dependent variable usually represents the density of the substance. These equations require initial conditions (what the initial concentration of the substance is) as well as boundary conditions (to specify, for instance, whether the substance can cross the boundary or not).
The most general second order linear partial differential equation with constant coefficients
\[
\sum_{i,j=1}^n a_{ij}\, \frac{\partial^2 u}{\partial x_i \,\partial x_j}
+ \sum_{i=1}^n b_i \,\frac{\partial u}{\partial x_i} +cu=d,
\]
may be placed in the form
\[
\lambda_1 \frac{\partial^2 u}{\partial x_1^2} +
\lambda_2 \frac{\partial^2 u}{\partial x_2^2} + \cdots +
\lambda_n \frac{\partial^2 u}{\partial x_n^2} + \cdots = d .
\]
If any of the λi are zero,
then the given equation is parabolic at some point.
Heat Transfer Equations
Heat flows by conduction when different points in a body are at different temperatures. In general, the temperature may vary both in space (x, y, and z) and time (t). The analytic investigation of the heat conduction is basically a study of the space-time variation of the temperature and determination of the equation:
\[
T = T(x,y,z,t) ,
\]
which is the mathematical expression of the temperature field. The above equation represents a set of temperatures at all points of the space at any given time.
The temperature field expressed by the above equation is referred to as transient or
non-steady temperature field. If the temperature in the space does not change with
time, the temperature field is a function of the space coordinates only and it is
referred to as a steady state. Mathematically, it is expressed as
\[
T = T(x,y,z) \qquad \iff \qquad \frac{\partial T}{\partial t} =0 .
\]
The temperature fields represented by these equations are three-dimensional
fields since they are function of three coordinates.
A temperature field, which is a function of two coordinates say x and y, is known as two-dimensional and is described by the equation
\[
T = T(x,y) \qquad \iff \qquad \frac{\partial T}{\partial z} =\frac{\partial T}{\partial t} =0 .
\]
Similarly, a temperature distribution is said to be one-dimensional if and only if
\[
T = T(x) \qquad \mbox{or} \qquad T= T(x,t) .
\]
The thermal conduction means the energy transformation between parts of
continuum by the transfer of kinetic energy between
particles or groups of particles at the atomic level.
It is described by an empirical relationship, known as the Fourier law (1822):
\[
{\bf q} = - k \, \nabla T ,
\]
where q is the thermal flux (the amount of heat conducted across a unit cross-sectional area in unit time, with units W/m2), k is the thermal conductivity with units J/(m*s*K), and
\( \nabla T \) is the temperature gradient.
The above law can be written in another form:
\[
{\bf Q} = - \left( k \cdot A \right) \nabla T ,
\]
where Q is the heat flow rate by conduction, measured in Watts (W); k is the thermal conductivity or heat conduction coefficient, with units W/(m*K); A is the cross sectional area normal to flow, with unites m2, \( \nabla T \) is the temperature gradient, with units K/m. This equation removes the impediment
of having to measure q directly and reduces to measurement of local temperatures (with thermometers or thermocouples).
Heat flux is defined as the amount of heat transferred per unit area per unit time from or to a surface.
A major source of heat loss from a house is through the windows. Calculate the rate of heat flux through a glass window of dimensions 1 m x 0.80 m with width 3 mm, if the temperatures at the inner and outer surfaces are 17℃ and 15℃, respectively. Calculate the heat flux through this window.
Solution: At this point, we know the temperatures at the surfaces of material. These temperatures are given also by conditions inside the house and outside the house. In this case, heat flows by conduction through the glass from the higher inside temperature to the lower outside temperature. We use the Fourier’s law of thermal conduction equation
\[
{\bf q} = - k \, \nabla T ,
\]
where q is the local heat flux density [W/m2]; k is the glass conductivity, which we take as 0.96 W/(m K); and ΔT is approximation to the gradient. The heat flux then be
\[
q \approx 0.96\,\left[\mbox{W}/\mbox{m}/\mbox{K} \right] \times 2 \left[ \mbox{K} \right] /3.0 \times 10^{-2} \left[ \mbox{m}\right] = 640 \left[ \mbox{W}/\mbox{m}^2 \right] .
\]
The total heat loss through this window will be
\[
Q = q\cdot A = 640 \times 1.0 \times 0.8 = 512\,\mbox{Watts}.
\]
The above Fourier equations are sufficient to describe heat conduction under steady-state conditions, that is, where the temperature at each point in the conducting medium is invariant and the flux is constant in time and space. To obtain the second law of heat conduction, the principle of energy conservation in the form of the continuity equation is invoked, which states that, in the absence of any sources or sinks of heat, the time rate of change in heat content of a volume element of the conducting medium must equal the change of flux with distance:
\[
\rho\,c_m \,\frac{\partial T}{\partial t} = -\nabla \cdot {\bf q} ,
\]
where ρ is mass density, cm specific heat capacity per unit mass (called simply ‘specific heat’ and defined as the change in heat content of a unit mass of the body per unit change in temperature). Combining the above equation with the Fourier law, we obtain
\[
\rho\,c_m \,\frac{\partial T}{\partial t} = \nabla \cdot \left( k\,\nabla T \right) .
\]
In one-dimensional form it becomes
\[
\rho\,c_m \,\frac{\partial T}{\partial t} = \frac{\partial}{\partial x} \left( k\, \frac{\partial T}{\partial x} \right) .
\]
Sometimes there is need to account for the possible occurrence of heat sources or sinks in the realm where heat flow takes place. Heat sources include such phenomena as organic matter decomposition, wetting of initially dry soil material, passing electricity, nuclear fission, and condensation of water vapor. Heat sinks are generally associated with evaporation. Combining all these sources and sinks into a single term S, we can rewrite the last equation:
\[
\rho\,c_m \,\frac{\partial T}{\partial t} = \frac{\partial}{\partial x} \left( k\, \frac{\partial T}{\partial x} \right) + S(x,t).
\]
If different points of the body having the same temperature are joined, we obtain an isothermal surface. Intersection of such isothermal surfaces by a plane gives a family of isotherms on this plane as shown below.
The temperature gradient is a vector normal to the isothermal surface. Mathematically, it is the derivative in this direction
\[
\nabla T = {\bf n} \,\frac{\partial T}{\partial {\bf n}} .
\]
where n is the unit vector normal to the isothermal surface and \( \partial T / \partial {\bf n} \) is the temperature
derivative along the normal. Projections of the vector \( \nabla T \) can be made on the
coordinate axes ox, oy, and oz. Accordingly, the components of the normal heat flow vector q (or thermal flux) are
\[
q_x = -k\,\frac{\partial T}{\partial x} , \quad q_y = -k\,\frac{\partial T}{\partial y} , \quad q_z = -k\,\frac{\partial T}{\partial z} .
\]
These equations state that the heat flow rate in, for instance, the x direction is directly
proportional to the temperature gradient and the cross sectional area A
normal to the heat flow. The minus sign indicates that the heat flow is positive in the direction of decreasing temperature.
Example: Consider the initial value problem
\[
\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \qquad (-\infty < x < \infty , \ 0< t ), \quad u(x,0) = f(x) .
\]
Assuming that the unknown function u(x, t) and its derivative approach zero at infinity, we apply the Fourier transformation to the given problem. So we multiply both sides by \( e^{{\bf j} x\xi} \) and integrated the result, we obtain (after integration by parts) the ODE subject to the initial condition:
\[
\frac{{\text d} u^F}{{\text d}t} + \alpha\,\xi^2 u^F = 0 , \qquad u^F (0) = f^F (\xi ) ,
\]
where
\[
u^F = ℱ\left[ u(x,t) \right] = \int_{-\infty}^{\infty} u(x,t)\, e^{{\bf j} x\xi}\, {\text d}x \qquad\mbox{and} \qquad f^F = ℱ\left[ f(x) \right] = \int_{-\infty}^{\infty} f(x)\, {\text d}x
\]
are Fourier transforms of u(x, t) and f(x), respectively. The above initial value problem is easy to solve:
\[
u^F = f^F (\xi )\, e^{-\alpha \xi^2 t} .
\]
Then taking inverse Fourier transform, we obtain the required solution
\[
u(x,t) = ℱ^{-1}\left[ f^F (\xi )\, e^{-\alpha \xi^2 t} \right] = ℱ^{-1}\left[ f^F (\xi )\, G(\xi , t) \right] ,
\]
where G is the Green's function. Its inverse Fourier transform is
\begin{align*}
g(x,t) &= ℱ^{-1}\left[ G(\xi , t) \right] = ℱ^{-1}\left[ e^{-\alpha \xi^2 t} \right]
\\
&= \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-\alpha \xi^2 t} e^{-{\bf j}x\xi} \, {\text d} \xi = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-\alpha \xi^2 t- {\bf j}x\xi} \, {\text d} \xi
\\
&= \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-\alpha t \left( \xi + {\bf j}x/2\alpha t \right)^2 - x^2 /(4\alpha t)} \, {\text d} \xi
\\
&= e^{- x^2 /(4\alpha t)} \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-\alpha t \left( \xi + {\bf j}x/2\alpha t \right)^2} \, {\text d} \xi = e^{- x^2 /(4\alpha t)} \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-\alpha t s} \, {\text d} s
\\
&= \frac{1}{\sqrt{4\alpha \pi t}} \, e^{- x^2 /(4\alpha t)} .
\end{align*}
With the Green's function at hand, we find the solution to the initial value problem
\[
u(x,t) = \frac{1}{\sqrt{4\alpha \pi t}} \int_{-\infty}^{\infty} e^{- \frac{(x-s)^2}{4\alpha \,t}}\,f(s)\, {\text d} s .
\]
Assuming[A > 0, Integrate[Exp[-A*s^2], {s, -Infinity, Infinity}]]
Sqrt[\[Pi]]/Sqrt[A]
■ Example: Consider the initial value problem
\[
\frac{\partial u}{\partial t} = \alpha\, \frac{\partial^2 u}{\partial x^2} \qquad (-\infty < x < \infty , \ 0< t ), \quad u(x,0) = f(x) .
\]
Assuming that the unknown function u(x, t) and its derivative approach zero at infinity, we apply the Laplace transformation to the given problem. So we multiply both sides by \( e^{-\lambda\,t} \) and integrated the result with respect to t, we obtain (after integration by parts) the ODE subject to the radiation condition:
\[
\frac{{\text d}^2 u^L}{{\text d}x^2} = \frac{\lambda}{\alpha}\, u^L - \frac{f(x)}{\alpha} , \qquad \lim_{x\to \pm \infty}\,u^L (x,\lambda ) = 0 ,
\]
where
\[
u^L = {\cal L}\left[ u(x,t) \right] = \int_{0}^{\infty} u(x,t)\, e^{-\lambda\,t}\, {\text d}t
\]
is the Laplace transforms of u(x, t). The above ordinary differential equation is easy to solve using the variation of parameters:
\[
u^L = A(x)\, e^{\omega\,x} + B(x)\, e^{-\omega\,x} , \qquad \omega^2 = \frac{\lambda}{\alpha} ,
\]
where unknown functions A(x) and B(x) are solutions of the following system of equations
\[
\begin{split}
A'(x)\, e^{\omega\,x} + B' (x)\, e^{-\omega\,x} &= 0 ,
\\
\omega\,A'(x)\, e^{\omega\,x} - \omega\,B' (x)\, e^{-\omega\,x} &= -\frac{f(x)}{\alpha} .
\end{split}
\]
We use Mathematica to solve the above system to obtain
Solve[{a*Exp[omega*x] + b*Exp[-omega*x] == 0,
omega*a*Exp[omega*x] - omega* b*Exp[-omega*x] == f[x]/alpha}, {a,
b}]
{{a -> (E^(-omega x) f[x])/(2 alpha omega),
b -> -((E^(omega x) f[x])/(2 alpha omega))}}
\[
\begin{split}
A'(x) &= -\frac{f(x)}{2\alpha\,\omega}\, e^{\omega\,x} \qquad\Longrightarrow \qquad A(x) = \frac{1}{2\,\sqrt{\lambda\alpha}}\, \int_{-x}^{\infty} e^{\omega\,\xi} \,f(\xi )\,{\text d}\xi ,
\\
B'(x) &= \frac{f(x)}{2\alpha\,\omega}\, e^{-\omega\,x} \qquad\Longrightarrow \qquad B(x) = \frac{1}{2\,\sqrt{\lambda\alpha}}\, \int_{-\infty}^x e^{\omega\,\xi} \,f(\xi )\,{\text d}\xi ,
\end{split}
\]
This yields the solution
\[
u^L (x,\lambda ) = \frac{1}{2\,\sqrt{\lambda\alpha}}\, \int_{-\infty}^{\infty} e^{-\omega\left\vert x-\xi \right\vert} \,f(\xi )\,{\text d}\xi .
\]
Then taking inverse Laplace transform, we obtain the required solution
\[
u(x,t) = {\cal L}^{-1} \left[ u^L \right] = \frac{1}{2\,\sqrt{\pi\alpha\,t}}\,\int_{-\infty}^{\infty} e^{-\left\vert x-\xi \right\vert^2 /(4\alpha\,t)} \,f(\xi )\,{\text d}\xi
\]
because the inverse Laplace transform is
\[
{\cal L}^{-1} \left[ \frac{1}{\sqrt{\lambda}}\, e^{- A\,\sqrt{\lambda}} \right] = \frac{1}{\sqrt{\pi\,t}}\, e^{-A^2 /(4t)} .
\]
Assuming[a > 0,
InverseLaplaceTransform[Exp[-a*Sqrt[s]]/Sqrt[s], s, t]]
ConditionalExpression[E^(-(a^2/(4 t)))/(Sqrt[\[Pi]] Sqrt[t]), a > 0]
■ Example: The initial value problem
\begin{align*}
&\mbox{Heat equation:} \qquad &u_t &= \alpha\,u_{xx} \qquad\mbox{for } x > 0 \mbox{ and } t > 0, \\
&\mbox{Initial condition:} \qquad &u(x,0) &= 0 , \qquad x > 0, \\
&\mbox{Regularity condition:} \qquad & \mbox{no condition} &
\end{align*}
The function
\[
u(x,t) = \sum_{n\ge 0} \frac{g^{(n)} (t)}{(2n)!}\, x^{2n} ,
\]
where
\[
g(t) = \begin{cases}
e^{-1/t^2} , & \ t> 0, \\
0, & \ t\le 0,
\end{cases}
\]
is a formal solution of the given IVP. The function g(t) is infinitely differentiable for −∞ < t < t. The series for the sum-function u(x, t) converges (together with all its derivatives) on compact (= bounded and closed) subsets of the xt-plane. Moreover, there exists a constant θ, 0 < θ < 1, such that
\[
\left\vert u(x,t) \right\vert \le \exp \left\{ \frac{x^2}{\theta\,t} - \frac{1}{2 t^2} \right\} , \qquad -\infty < x < \infty , \quad t> 0.
\]
Note that
\[
\frac{x^2}{\theta\,t} - \frac{1}{2 t^2} \le \frac{x^4}{2\theta^2} , \qquad -\infty < x < \infty , \quad t> 0.
\]
■
A hyperbolic model for heat conduction
\[
\frac{1}{C^2} \,\frac{\partial^2 \theta}{\partial t^2} + \frac{1}{\alpha} \,\frac{\partial \theta}{\partial t} = \nabla^2 \theta .
\]
In this equation, C is called the speed of second sound (i.e. the fictitious quantum particles, photons), which was predicted by Tisza (1938) and Landau (1941). The equation is known as the hyperbolic heat conduction (HHC) equation. Mathematically, it is the same as the telegrapher's equation, which is derived from Maxwell’s equations of electrodynamics. The main reason of this model is to overcome instantaneous change in temperature, θ. The above telegrapher's equation is attributed to Cattàneo, Vernotte, and Chester.
Nonlinear Parabolic Equations
Example: Consider nonlinear diffusion equation
\[
\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \left( u^2 \frac{\partial u}{\partial x} \right) \qquad \Longleftrightarrow \qquad
u_t = u^2 u_{xx} + 2y\,u_x ,
\]
subject to the initial condition
\[
u(x,0) = \frac{x+h}{2\sqrt{c}} ,
\]
where c > 0 and h are arbitrary constants.
The exact solution of the given equation is (according to C. Chun et al)
\[
u(x,0) = \frac{x+h}{2\sqrt{c-t}} ,
\]
Assuming that the solution is represented by a convergent power series
\[
u(x,t) = \sum_{i\ge 0} \sum_{j\ge 0} a_{i,j} x^i t^j ,
\]
we substitute it into the diffusion equation. Using formulas
\begin{align*}
u_t &= \sum_{i\ge 0} \sum_{j\ge 0} \left( j+1 \right) a_{i,j+1} x^i t^j ,
\\
u_x &= \sum_{i\ge 0} \sum_{j\ge 0} \left( i+1 \right) a_{i+1,j} x^i t^j ,
\\
u_{xx} &= \sum_{i\ge 0} \sum_{j\ge 0} \left( i+1 \right) \left( i+2 \right)
a_{i+2,j} x^i t^j ,
\\
u^2 &= \sum_{i\ge 0} \sum_{j\ge 0} \left[ \sum_{r=0}^i \sum_{s=0}^j a_{r,s} a_{i-r,j-s} \right] x^i t^j ,
\\
u^2 u_{xx} &= \sum_{i\ge 0} \sum_{j\ge 0} \left[ \sum_{r=0}^i \sum_{s=0}^j
\left( i-p+1 \right) \left( i-p+2 \right) a_{i-p+2,j-q} \left(
\sum_{r=0}^p \sum_{s=0}^q a_{r,s} a_{p-i,q-j} \right) \right] x^i t^j ,
\\
u\,u_x^2 &= \sum_{i\ge 0} \sum_{j\ge 0} \left[ \sum_{r=0}^i \sum_{s=0}^j
a_{i-p.j-q} \left( \sum_{r=0}^p \sum_{S=0}^Q \left( r+1 \right) \left( p-r+1
\right) a_{r+1,s} a_{p-r+1, q-s} \right) \right] x^i t^j ,
\end{align*}
we obtain the recurrence relation:
\begin{align*}
a_{i, j+1} &= \frac{1}{j+1} \left[ \sum_{q=0}^j \sum_{p=0}^i \left( i-p+1 \right) \left( i-p+2 \right) a_{i-p+2, j-q} \left( \sum_{r=0}^i \sum_{s=0}^j a_{r,s} a_{p-r, q-s} \right) + \right.
\\
&\quad \left. 2\sum_{q=0}^j \sum_{p=0}^i a_{i-p, j-q} \left(
\sum_{r=0}^p \sum_{s=0}^q \left( r+1 \right) \left( p-r+1 \right) a_{r+1, s} a_{p-r+1,q-s} \right) \right] ,
\end{align*}
where
\[
a_{i,0} =
\begin{cases}
\frac{h}{2\sqrt{c}} , & \ \mbox{if } i=0, \\
\frac{1}{2\sqrt{c}} , & \ \mbox{if } i=1, \\
0, & \ \mbox{ otherwise}.
\end{cases}
\]
By applying the above recurrence relations for several values of
i and j, the polynomial approximation for
u(x,t) is obtained
\[
u(x,t) = \frac{h}{2\sqrt{c}} + \frac{ht}{4\,c^{3/2}} + \frac{3h\,t^2}{16\,c^{5/2}} + \frac{5h\,t^3}{32\,c^{7/2}} + \frac{x}{2\sqrt{c}} +
\frac{xt}{c^{3/2}} + \frac{3t^2 x}{16\,c^{5/2}} + \frac{5t^3 x}{12\,c^{7/2}} + \cdots .
\]
■ Example: Consider the nonlinear PDE
\[
u_t = u^2 u_{xx} + u^2 + u_x +u , \qquad u(x,0) = \cos x .
\]
The true solution of the above initial value problem is
\( u(x,t) = e^t \cos (x+t) . \)
D[u[x, t], t] - (u[x, t])^2*D[u[x, t], x, x] - (u[x, t])^3 -
1*D[u[x, t], x] - u[x, t]
Suppose that the solution is represented as a power series:
\[
u(x,t) = \sum_{i\ge 0} \sum_{j\ge 0} a_{i,j} x^i t^j .
\]
Then
\begin{align*}
u_t &= \sum_{i\ge 0} \sum_{j\ge 0} \left( j+1 \right) a_{i,j+1} x^i t^j ,
\\
u_x &= \sum_{i\ge 0} \sum_{j\ge 0} \left( i+1 \right) a_{i+1,j} x^i t^j ,
\\
u_{xx} &= \sum_{i\ge 0} \sum_{j\ge 0} \left( i+1 \right) \left( i+2 \right)
a_{i+2,j} x^i t^j ,
\\
u^2 &= \sum_{i\ge 0} \sum_{j\ge 0} \left[ \sum_{r=0}^i \sum_{s=0}^j a_{r,s} a_{i-r, j-s} \right] x^i t^j .
\end{align*}
Therefore,
\[
u^2 u_{xx} = \sum_{i\ge 0} \sum_{j\ge 0} \left[ \sum_{p=0}^i \sum_{q=0}^j \left( i-1 \right) \left( i-p+2 \right) a_{i-p+2,j-q} \left( \sum_{r=0}^p \sum_{s=0}^q a_{r,s} a_{p-r,q-s} \right) \right] x^i t^j ,
\]
and
\[
u^3 = \sum_{i\ge 0} \sum_{j\ge 0} \left[ \sum_{r=0}^i \sum_{s=0}^j \left( \sum_{p=0}^r \sum_{q=0}^s i-1 a_{p,q}a_{r-p, s-q} \right) a_{i-r,j-s} \right] x^i t^j .
\]
Substituting the above series into the given partial equation, the following recurrence equation is obtained
\begin{align*}
a_{i, j+1} &= \frac{1}{j+1} \left[ \sum_{q=0}^j \sum_{p=0}^i \left( i-p+1 \right) \left( i-p+2 \right) a_{i-p+2, j-q} \left( \sum_{r=0}^i \sum_{s=0}^j a_{r,s} a_{p-r, q-s} \right) \left( \sum_{r=0}^p \sum_{s=0}^q a_{r,s} a_{p-r,q-s} \right) \right.
\\
&\quad + \left. \left( \sum_{s=0}^j \sum_{r=0}^i \left( \sum_{p=0}^i \sum_{q=0}^s a_{p, q} a_{r-p, q-s} \right) a_{i-r,j-s} \right) +
\left( i+1 \right) a_{r+1, j} + a_{i,j} \right) \right] .
\end{align*}
The approximate solution is
\[
u(x,t) = 1 - \frac{t^2}{2} -xt + \frac{x t^3}{6} - \frac{x^2}{2} + \frac{x^2 t^2}{4} + \frac{x^3 t}{6} - \frac{x^3 t^3}{36} + \cdots .
\]
■
- Chun, C., Jafari, H., Kim, Y., Numerical method for the wave and the nonlinear diffusion equations with the homotopy perturbation method, Computers & Mathematics with Applications, 2009, Volume 57, Issue 7, pp. 1226--1231. https://doi.org/10.1016/j.camwa.2009.01.013
- Grigoryan, V, Parial Differential Equations, 2010, University of California, Santa Barbara.
- Nuseir, A.S., and Al-Hasoon, A., Power series solutions for nonlinear systems of partial differential equations, Applied Mathematical Sciences, Vol. 6, 2012, no. 104, 5147 - 5159
Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions