## Solving Partial Differential Equations

In a partial differential equation (PDE), the function being solved for depends on several variables, and the differential equation can include partial derivatives taken with respect to each of the variables. Partial differential equations are useful for modelling waves, heat flow, fluid dispersion, and other phenomena with spatial behavior that changes over time.

### What Types of PDEs Can You Solve with MATLAB?

The MATLAB® PDE solver `pdepe` solves initial-boundary value problems for systems of PDEs in one spatial variable x and time t. You can think of these as ODEs of one variable that also change with respect to time.

`pdepe` uses an informal classification for the 1-D equations it solves:

• Equations with a time derivative are parabolic. An example is the heat equation $\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial {x}^{2}}$.

• Equations without a time derivative are elliptic. An example is the Laplace equation $\frac{{\partial }^{2}u}{\partial {x}^{2}}=0$.

`pdepe` requires at least one parabolic equation in the system. In other words, at least one equation in the system must include a time derivative.

`pdepe` also solves certain 2-D and 3-D problems that reduce to 1-D problems due to angular symmetry (see the argument description for the symmetry constant `m` for more information).

Partial Differential Equation Toolbox™ extends this functionality to generalized problems in 2-D and 3-D with Dirichlet and Neumann boundary conditions.

### Solving 1-D PDEs

A 1-D PDE includes a function u(x,t) that depends on time t and one spatial variable x. The MATLAB PDE solver `pdepe` solves systems of 1-D parabolic and elliptic PDEs of the form

`$c\left(x,t,u,\frac{\partial u}{\partial x}\right)\frac{\partial u}{\partial t}={x}^{-m}\frac{\partial }{\partial x}\left({x}^{m}f\left(x,t,u,\frac{\partial u}{\partial x}\right)\right)+s\left(x,t,u,\frac{\partial u}{\partial x}\right).$`

The equation has the properties:

• The PDEs hold for t0ttf and axb.

• The spatial interval [a, b] must be finite.

• `m` can be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry, respectively. If m > 0, then a ≥ 0 must also hold.

• The coefficient $f\left(x,t,u,\frac{\partial u}{\partial x}\right)$ is a flux term and $s\left(x,t,u,\frac{\partial u}{\partial x}\right)$ is a source term.

• The flux term must depend on the partial derivative u/∂x.

The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix $c\left(x,t,u,\frac{\partial u}{\partial x}\right)$. The diagonal elements of this matrix are either zero or positive. An element that is zero corresponds to an elliptic equation, and any other element corresponds to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if they are mesh points (points where the solution is evaluated). Discontinuities in c and s due to material interfaces are permitted provided that a mesh point is placed at each interface.

#### Solution Process

To solve PDEs with `pdepe`, you must define the equation coefficients for c, f, and s, the initial conditions, the behavior of the solution at the boundaries, and a mesh of points to evaluate the solution on. The function call `sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)` uses this information to calculate a solution on the specified mesh:

Together, the `xmesh` and `tspan` vectors form a 2-D grid that `pdepe` evaluates the solution on.

#### Equations

You must express the PDEs in the standard form expected by `pdepe`. Written in this form, you can read off the values of the coefficients c, f, and s.

In MATLAB you can code the equations with a function of the form

```function [c,f,s] = pdefun(x,t,u,dudx) c = 1; f = dudx; s = 0; end```
In this case `pdefun` defines the equation $\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial {x}^{2}}$. If there are multiple equations, then c, f, and s are vectors with each element corresponding to one equation.

#### Initial Conditions

At the initial time t = t0, for all x, the solution components satisfy initial conditions of the form

`$u\left(x,{t}_{0}\right)={u}_{0}\left(x\right).$`

In MATLAB you can code the initial conditions with a function of the form

```function u0 = icfun(x) u0 = 1; end```
In this case `u0 = 1` defines an initial condition of u0(x,t0) = 1. If there are multiple equations, then `u0` is a vector with each element defining the initial condition of one equation.

#### Boundary Conditions

At the boundary x = a or x = b, for all t, the solution components satisfy boundary conditions of the form

`$p\left(x,t,u\right)+q\left(x,t\right)f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0.$`

q(x,t) is a diagonal matrix with elements that are either zero or never zero. Note that the boundary conditions are expressed in terms of the flux f, rather than the partial derivative of u with respect to x. Also, of the two coefficients p(x,t,u) and q(x,t), only p can depend on u.

In MATLAB you can code the boundary conditions with a function of the form

```function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t) pL = uL; qL = 0; pR = uR - 1; qR = 0; end ```
`pL` and `qL` are the coefficients for the left boundary, while `pR` and `qR` are the coefficients for the right boundary. In this case `bcfun` defines the boundary conditions

`$\begin{array}{l}{u}_{L}\left({x}_{L},t\right)=0\\ {u}_{R}\left({x}_{R},t\right)=1\end{array}$`

If there are multiple equations, then the outputs `pL`, `qL`, `pR`, and `qR` are vectors with each element defining the boundary condition of one equation.

#### Integration Options

The default integration properties in the MATLAB PDE solver are selected to handle common problems. In some cases, you can improve solver performance by overriding these default values. To do this, use `odeset` to create an `options` structure. Then, pass the structure to `pdepe` as the last input argument:

```sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) ```

Of the options for the underlying ODE solver `ode15s`, only those shown in the following table are available for `pdepe`.

Category

Option Name

Error control

Step-size

Event logging

#### Evaluating the Solution

After you solve an equation with `pdepe`, MATLAB returns the solution as a 3-D array `sol`, where `sol(i,j,k)` contains the `k`th component of the solution evaluated at `t(i)` and `x(j)`. In general, you can extract the `k`th solution component with the command `u = sol(:,:,k)`.

The time mesh you specify is used purely for output purposes, and does not affect the internal time steps taken by the solver. However, the spatial mesh you specify can affect the quality and speed of the solution. After solving an equation, you can use `pdeval` to evaluate the solution structure returned by `pdepe` with a different spatial mesh.

### Example: The Heat Equation

An example of a parabolic PDE is the heat equation in one dimension:

`$\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{{\partial }^{2}\mathit{u}}{\partial {\mathit{x}}^{2}}.$`

This equation describes the dissipation of heat for $0\le \mathit{x}\le \mathit{L}$ and $\mathit{t}\ge 0$. The goal is to solve for the temperature $\mathit{u}\left(\mathit{x},\mathit{t}\right)$. The temperature is initially a nonzero constant, so the initial condition is

`$\mathit{u}\left(\mathit{x},0\right)={\mathit{T}}_{0}.$`

Also, the temperature is zero at the left boundary, and nonzero at the right boundary, so the boundary conditions are

`$\mathit{u}\left(0,\mathit{t}\right)=0,$`

`$\mathit{u}\left(\mathit{L},\mathit{t}\right)=1.$`

To solve this equation in MATLAB®, you need to code the equation, initial conditions, and boundary conditions, then select a suitable solution mesh before calling the solver `pdepe`. You either can include the required functions as local functions at the end of a file (as in this example), or save them as separate, named files in a directory on the MATLAB path.

Code Equation

Before you can code the equation, you need to make sure that it is in the form that the `pdepe` solver expects:

`$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)\frac{\partial \mathit{u}}{\partial \mathit{t}}={\mathit{x}}^{-\mathit{m}}\frac{\partial }{\partial \mathit{x}}\left({\mathit{x}}^{\mathit{m}}\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)\right)+\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right).$`

In this form, the heat equation is

`$1\cdot \frac{\partial \mathit{u}}{\partial \mathit{t}}={\mathit{x}}^{0}\frac{\partial }{\partial \mathit{x}}\left({\mathit{x}}^{0}\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)+0.$`

So the values of the coefficients are as follows:

• $\mathit{m}=0$

• $\mathit{c}=1$

• $\mathit{f}=\frac{\partial \mathit{u}}{\partial \mathit{x}}$

• $\mathit{s}=0$

The value of $\mathit{m}$ is passed as an argument to `pdepe`, while the other coefficients are encoded in a function for the equation, which is

```function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end ```

(Note: All functions are included as local functions at the end of the example.)

Code Initial Condition

The initial condition function for the heat equation assigns a constant value for ${\mathit{u}}_{0}$. This function must accept an input for $\mathit{x}$, even if it is unused.

```function u0 = heatic(x) u0 = 0.5; end ```

Code Boundary Conditions

The standard form for the boundary conditions expected by the `pdepe` solver is

`$\mathit{p}\left(\mathit{x},\mathit{t},\mathit{u}\right)+\mathit{q}\left(\mathit{x},\mathit{t}\right)\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$`

Written in this form, the boundary conditions for this problem are

`$\mathit{u}\left(0,\mathit{t}\right)+\left(0\cdot \mathit{f}\right)=0,$`

`$\left(\mathit{u}\left(\mathit{L},\mathit{t}\right)-1\right)+\left(0\cdot \mathit{f}\right)=0.$`

So the values for $\mathit{p}$ and $\mathit{q}$ are

• ${\mathit{p}}_{\mathit{L}}={\mathit{u}}_{\mathit{L}},{\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\text{\hspace{0.17em}}\mathit{q}}_{\mathit{L}}=0.$

• ${\mathit{p}}_{\mathit{R}}={\mathit{u}}_{\mathit{R}}-1,{\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathit{q}}_{\mathit{R}}=0$.

The corresponding function is then

```function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end ```

Select Solution Mesh

Use a spatial mesh of 20 points and a time mesh of 30 points. Since the solution rapidly reaches a steady state, the time points near $\mathit{t}=0$ are more closely spaced together to capture this behavior in the output.

```L = 1; x = linspace(0,L,20); t = [linspace(0,0.05,20), linspace(0.5,5,10)];```

Solve Equation

Finally, solve the equation using the symmetry $\mathit{m}$, the PDE equation, the initial condition, the boundary conditions, and the meshes for $\mathit{x}$ and $\mathit{t}$.

```m = 0; sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);```

Plot Solution

Use `imagesc` to visualize the solution matrix.

```colormap hot imagesc(x,t,sol) colorbar xlabel('Distance x','interpreter','latex') ylabel('Time t','interpreter','latex') title('Heat Equation for \$0 \le x \le 1\$ and \$0 \le t \le 5\$','interpreter','latex')``` Local Functions

```function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end function u0 = heatic(x) u0 = 0.5; end function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end```

### PDE Examples and Files

Several available example files serve as excellent starting points for most common 1-D PDE problems. To explore and run examples, use the Differential Equations Examples app. To run this app, type

`odeexamples`

To open an individual file for editing, type

`edit exampleFileName.m`

To run an example, type

`exampleFileName`

This table contains a list of the available PDE example files.

Example File

Description

`pdex1`

Simple PDE that illustrates the formulation, computation, and plotting of the solution.

`pdex2`

Problem that involves discontinuities.

`pdex3`

Problem that requires computing values of the partial derivative.

`pdex4`

System of two PDEs whose solution has boundary layers at both ends of the interval and changes rapidly for small t.

`pdex5`

System of PDEs with step functions as initial conditions.

 Skeel, R. D. and M. Berzins, "A Method for the Spatial Discretization of Parabolic Equations in One Space Variable," SIAM Journal on Scientific and Statistical Computing, Vol. 11, 1990, pp. 1–32.