# Documentation

## Partial Differential Equations

### Function Summary

#### PDE Solver

This is the MATLAB® PDE solver.

PDE Initial-Boundary Value Problem Solver

Description

Solve initial-boundary value problems for systems of parabolic and elliptic PDEs in one space variable and time.

#### PDE Helper Function

PDE Helper Function

Description

Evaluate the numerical solution of a PDE using the output of `pdepe`.

### Initial Value Problems

`pdepe` solves systems of parabolic and elliptic PDEs in one spatial variable x and time t, 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).$ (10-5)

The PDEs hold for t0ttf and axb. The 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.

In Equation 10-5, f(x,t,u,∂u/∂x) is a flux term and s(x,t,u,∂u/∂x) is a source term. The flux term must depend on ∂u/∂x. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix c(x,t,u,∂u/∂x). The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise 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. Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface.

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).$ (10-6)

At the boundary x = a or x = b, for all t the solution components satisfy a boundary condition 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.$ (10-7)

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

### PDE Solver

#### The PDE Solver

The MATLAB PDE solver, `pdepe`, solves initial-boundary value problems for systems of parabolic and elliptic PDEs in the one space variable x and time t. There must be at least one parabolic equation in the system.

The `pdepe` solver converts the PDEs to ODEs using a second-order accurate spatial discretization based on a fixed set of user-specified nodes. The discretization method is described in [9]. The time integration is done with `ode15s`. The `pdepe` solver exploits the capabilities of `ode15s` for solving the differential-algebraic equations that arise when Equation 10-5 contains elliptic equations, and for handling Jacobians with a specified sparsity pattern. `ode15s` changes both the time step and the formula dynamically.

After discretization, elliptic equations give rise to algebraic equations. If the elements of the initial conditions vector that correspond to elliptic equations are not "consistent" with the discretization, `pdepe` tries to adjust them before beginning the time integration. For this reason, the solution returned for the initial time may have a discretization error comparable to that at any other time. If the mesh is sufficiently fine, `pdepe` can find consistent initial conditions close to the given ones. If `pdepe` displays a message that it has difficulty finding consistent initial conditions, try refining the mesh. No adjustment is necessary for elements of the initial conditions vector that correspond to parabolic equations.

#### PDE Solver Syntax

The basic syntax of the solver is:

```sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) ```
 Note:   Correspondences given are to terms used in Initial Value Problems.

The input arguments are

 `m` Specifies the symmetry of the problem. `m` can be `0` = slab, `1` = cylindrical, or `2` = spherical. It corresponds to m in Equation 10-5. `pdefun` Function that defines the components of the PDE. It computes the terms c, f, and s in Equation 10-5, and has the form```[c,f,s] = pdefun(x,t,u,dudx) ```where `x` and `t` are scalars, and `u` and `dudx` are vectors that approximate the solution u and its partial derivative with respect to x. `c`, `f`, and `s` are column vectors. `c` stores the diagonal elements of the matrix c. `icfun` Function that evaluates the initial conditions. It has the form```u = icfun(x) ```When called with an argument `x`, `icfun` evaluates and returns the initial values of the solution components at `x` in the column vector `u`. `bcfun` Function that evaluates the terms p and q of the boundary conditions. It has the form ```[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t) ```where `ul` is the approximate solution at the left boundary `xl = a` and `ur` is the approximate solution at the right boundary `xr = b`. `pl` and `ql` are column vectors corresponding to p and the diagonal of q evaluated at `xl`. Similarly, `pr` and `qr` correspond to `xr`. When m > 0 and a = 0, boundedness of the solution near x = 0 requires that the f vanish at a = 0. `pdepe` imposes this boundary condition automatically and it ignores values returned in `pl` and `ql`. `xmesh` Vector [`x0`, `x1`, ..., `xn`] specifying the points at which a numerical solution is requested for every value in `tspan`. `x0` and `xn` correspond to a and b, respectively. Second-order approximation to the solution is made on the mesh specified in `xmesh`. Generally, it is best to use closely spaced mesh points where the solution changes rapidly. `pdepe` does not select the mesh in x automatically. You must provide an appropriate fixed mesh in `xmesh`. The cost depends strongly on the length of `xmesh`. When m > 0, it is not necessary to use a fine mesh near x = 0 to account for the coordinate singularity.The elements of `xmesh` must satisfy `x0` < `x1` < ... < `xn`. The length of `xmesh` must be ≥ `3`. `tspan` Vector [`t0`, `t1`, ..., `tf`] specifying the points at which a solution is requested for every value in `xmesh`. `t0` and `tf` correspond to t0 and tf , respectively. `pdepe` performs the time integration with an ODE solver that selects both the time step and formula dynamically. The solutions at the points specified in `tspan` are obtained using the natural continuous extension of the integration formulas. The elements of `tspan` merely specify where you want answers and the cost depends weakly on the length of `tspan`.The elements of `tspan` must satisfy `t0` < `t1` < ... < `tf`. The length of `tspan` must be ≥ `3`.

The output argument `sol` is a three-dimensional array, such that

• `sol`(`:`,`:`,`k`) approximates component `k` of the solution u.

• `sol`(`i`,`:`,`k`) approximates component `k` of the solution at time `tspan`(`i`) and mesh points `xmesh(:)`.

• `sol`(`i,j,k`) approximates component `k` of the solution at time `tspan`(`i`) and the mesh point `xmesh(j)`.

#### PDE Solver Options

For more advanced applications, you can also specify as input arguments solver options and additional parameters that are passed to the PDE functions.

 `options` Structure of optional parameters that change the default integration properties. This is the seventh input argument.```sol = pdepe(m,pdefun,icfun,bcfun,... xmesh,tspan,options) ```See Integrator Options for more information.

### Integrator 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 defaults. You do this by supplying `pdepe` with one or more property values in an options structure.

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

Use `odeset` to create the options structure. Only those options of the underlying ODE solver shown in the following table are available for `pdepe`. The defaults obtained by leaving off the input argument `options` are generally satisfactory. Integrator Options tells you how to create the structure and describes the properties.

PDE Properties

Properties Category

Property Name

Error control

`RelTol`, `AbsTol`, `NormControl`

Step-size

`InitialStep`, `MaxStep`

### Examples

#### Single PDE

Solving the Equation.  This example illustrates the straightforward formulation, solution, and plotting of the solution of a single PDE

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

This equation holds on an interval 0 ≤ x ≤ 1 for times t ≥ 0. At t = 0, the solution satisfies the initial condition

$u\left(x,0\right)=\mathrm{sin}\pi x.$

At x = 0 and x = 1, the solution satisfies the boundary conditions

$\begin{array}{l}u\left(0,t\right)=0\\ \pi {e}^{-t}+\frac{\partial u}{\partial x}\left(1,t\right)=0.\end{array}$

 Note:   The file, `pdex1.m``pdex1.m`, contains the complete code for this example. It contains all the required functions coded as local functions. To see the code in an editor, type `edit pdex1` at the command line. To run it, type `pdex1` at the command line. See PDE Solver Syntax for more information.
1. Rewrite the PDE. Write the PDE in 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).$

This is the form shown in Equation 10-5 and expected by `pdepe`. See Initial Value Problems for more information. For this example, the resulting equation is

${\pi }^{2}\frac{\partial u}{\partial t}={x}^{0}\frac{\partial }{\partial x}\left({x}^{0}\frac{\partial u}{\partial x}\right)+0$

with parameter m = 0 and the terms

$\begin{array}{l}c\left(x,t,u,\frac{\partial u}{\partial x}\right)={\pi }^{2}\\ f\left(x,t,u,\frac{\partial u}{\partial x}\right)=\frac{\partial u}{\partial x}\\ s\left(x,t,u,\frac{\partial u}{\partial x}\right)=0.\end{array}$

2. Code the PDE. Once you rewrite the PDE in the form shown above (Equation 10-5) and identify the terms, you can code the PDE in a function that `pdepe` can use. The function must be of the form

```[c,f,s] = pdefun(x,t,u,dudx) ```

where `c`, `f`, and `s` correspond to the c, f, and s terms. The code below computes `c`, `f`, and `s` for the example problem.

```function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = 0; ```
3. Code the initial conditions function. You must code the initial conditions in a function of the form

```u = icfun(x) ```

The code below represents the initial conditions in the function `pdex1ic`.

```function u0 = pdex1ic(x) u0 = sin(pi*x); ```
4. Code the boundary conditions function. You must also code the boundary conditions in a function of the form

```[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t) ```

The boundary conditions, written in the same form as Equation 10-7, are

and

The code below evaluates the components p(x,t,u) and q(x,t) of the boundary conditions in the function `pdex1bc`.

```function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = pi * exp(-t); qr = 1; ```

In the function `pdex1bc`, `pl` and `ql` correspond to the left boundary conditions (x = 0), and `pr` and `qr` correspond to the right boundary condition (x = 1).

5. Select mesh points for the solution. Before you use the MATLAB PDE solver, you need to specify the mesh points (t,x) at which you want `pdepe` to evaluate the solution. Specify the points as vectors `t` and `x`.

The vectors `t` and `x` play different roles in the solver (see PDE Solver). In particular, the cost and the accuracy of the solution depend strongly on the length of the vector `x`. However, the computation is much less sensitive to the values in the vector `t`.

This example requests the solution on the mesh produced by 20 equally spaced points from the spatial interval [0,1] and five values of `t` from the time interval [0,2].

```x = linspace(0,1,20); t = linspace(0,2,5); ```
6. Apply the PDE solver. The example calls `pdepe` with `m = 0`, the functions `pdex1pde`, `pdex1ic`, and `pdex1bc`, and the mesh defined by `x` and `t` at which `pdepe` is to evaluate the solution. The `pdepe` function returns the numerical solution in a three-dimensional array `sol`, where `sol(i,j,k)` approximates the `k`th component of the solution, uk, evaluated at `t(i)` and `x(j)`.

```m = 0; sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); ```

This example uses `@` to pass `pdex1pde`, `pdex1ic`, and `pdex1bc `as function handles to `pdepe`.

 Note:   See the `function_handle` (@), `func2str`, and `str2func` reference pages, and the `@` section of MATLAB Programming Fundamentals for information about function handles.
7. View the results. Complete the example by displaying the results:

1. Extract and display the first solution component. In this example, the solution u has only one component, but for illustrative purposes, the example "extracts" it from the three-dimensional array. The surface plot shows the behavior of the solution.

``` u = sol(:,:,1); surf(x,t,u) title('Numerical solution computed with 20 mesh points') xlabel('Distance x') ylabel('Time t') ```

2. Display a solution profile at tf , the final value of t. In this example, tf  = t = 2.

``` figure plot(x,u(end,:)) title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)') ```

Evaluating the Solution.  After obtaining and plotting the solution above, you might be interested in a solution profile for a particular value of `t`, or the time changes of the solution at a particular point `x`. The `k`th column `u(:,k)` (of the solution extracted in step 7) contains the time history of the solution at `x(k)`. The `j`th row `u(j,:)` contains the solution profile at `t(j)`.

Using the vectors `x` and `u(j,:)`, and the helper function `pdeval`, you can evaluate the solution `u` and its derivative ∂u/∂x at any set of points `xout`

```[uout,DuoutDx] = pdeval(m,x,u(j,:),xout) ```

The example `pdex3``pdex3` uses `pdeval` to evaluate the derivative of the solution at `xout = 0`. See `pdeval` for details.

#### System of PDEs

This example illustrates the solution of a system of partial differential equations. The problem is taken from electrodynamics. It has boundary layers at both ends of the interval, and the solution changes rapidly for small t.

The PDEs are

$\begin{array}{l}\frac{\partial {u}_{1}}{\partial t}=0.024\frac{{\partial }^{2}{u}_{1}}{\partial {x}^{2}}-F\left({u}_{1}-{u}_{2}\right)\\ \frac{\partial {u}_{2}}{\partial t}=0.170\frac{{\partial }^{2}{u}_{1}}{\partial {x}^{2}}+F\left({u}_{1}-{u}_{2}\right)\end{array}$

where F(y) = exp(5.73y) – exp(–11.46y). The equations hold on an interval 0 ≤ x ≤ 1 for times t ≥ 0.

The solution u satisfies the initial conditions

$\begin{array}{l}{u}_{1}\left(x,0\right)\equiv 1\\ {u}_{2}\left(x,0\right)\equiv 0\end{array}$

and boundary conditions

$\begin{array}{l}\frac{\partial {u}_{1}}{\partial x}\left(0,t\right)\equiv 0\\ {u}_{2}\left(0,t\right)\equiv 0\\ {u}_{1}\left(1,t\right)\equiv 1\\ \frac{\partial {u}_{2}}{\partial x}\left(1,t\right)\equiv 0.\end{array}$

 Note:   The file, `pdex4.m``pdex4.m`, contains the complete code for this example. It contains all the required functions coded as local functions. To see the code in an editor, type `edit pdex4` at the command line. To run it, type `pdex4` at the command line.
1. Rewrite the PDE. In the form expected by `pdepe`, the equations are

The boundary conditions on the partial derivatives of u have to be written in terms of the flux. In the form expected by `pdepe`, the left boundary condition is

$\left[\begin{array}{c}0\\ {u}_{2}\end{array}\right]+\left[\begin{array}{c}1\\ 0\end{array}\right]\cdot \ast \left[\begin{array}{c}0.024\left(\partial {u}_{1}/\partial x\right)\\ 0.170\left(\partial {u}_{2}/\partial x\right)\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$

and the right boundary condition is

$\left[\begin{array}{c}{u}_{1}-1\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ 1\end{array}\right]\cdot \ast \left[\begin{array}{c}0.024\left(\partial {u}_{1}/\partial x\right)\\ 0.170\left(\partial {u}_{2}/\partial x\right)\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$

2. Code the PDE. After you rewrite the PDE in the form shown above, you can code it as a function that `pdepe` can use. The function must be of the form

```[c,f,s] = pdefun(x,t,u,dudx) ```

where `c`, `f`, and `s` correspond to the c, f , and s terms in Equation 10-5.

```function [c,f,s] = pdex4pde(x,t,u,DuDx) c = [1; 1]; f = [0.024; 0.17] .* DuDx; y = u(1) - u(2); F = exp(5.73*y)-exp(-11.47*y); s = [-F; F]; ```
3. Code the initial conditions function. The initial conditions function must be of the form

```u = icfun(x) ```

The code below represents the initial conditions in the function `pdex4ic`.

```function u0 = pdex4ic(x); u0 = [1; 0]; ```
4. Code the boundary conditions function. The boundary conditions functions must be of the form

```[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t) ```

The code below evaluates the components p(x,t,u) and q(x,t) (Equation 10-7) of the boundary conditions in the function `pdex4bc`.

```function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t) pl = [0; ul(2)]; ql = [1; 0]; pr = [ur(1)-1; 0]; qr = [0; 1]; ```
5. Select mesh points for the solution. The solution changes rapidly for small t. The program selects the step size in time to resolve this sharp change, but to see this behavior in the plots, output times must be selected accordingly. There are boundary layers in the solution at both ends of [0,1], so mesh points must be placed there to resolve these sharp changes. Often some experimentation is needed to select the mesh that reveals the behavior of the solution.

```x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2]; ```
6. Apply the PDE solver. The example calls `pdepe` with `m = 0`, the functions `pdex4pde`, `pdex4ic`, and `pdex4bc`, and the mesh defined by x and t at which `pdepe` is to evaluate the solution. The `pdepe` function returns the numerical solution in a three-dimensional array `sol`, where `sol(i,j,k)` approximates the `k`th component of the solution, μk, evaluated at `t(i)` and `x(j)`.

```m = 0; sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t); ```
7. View the results. The surface plots show the behavior of the solution components.

```u1 = sol(:,:,1); u2 = sol(:,:,2); figure surf(x,t,u1) title('u1(x,t)') xlabel('Distance x') ylabel('Time t') ```

```figure surf(x,t,u2) title('u2(x,t)') xlabel('Distance x') ylabel('Time t') ```

The following additional examples are available. Type

`edit examplename`

to view the code and

`examplename`

to run the example.

Example Name

Description

`pdex1`

Simple PDE that illustrates the straightforward 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