Solve Single PDE

This example shows how to formulate, compute, and plot the solution to a single PDE.

Consider the partial differential equation

`${\pi }^{2}\frac{\partial \text{\hspace{0.17em}}\mathit{u}}{\partial \text{\hspace{0.17em}}\mathit{t}}=\frac{{\partial }^{2}\mathit{u}}{\partial {\mathit{x}}^{2}}.$`

The equation is defined on the interval $0\le \mathit{x}\le 1$ for times $\mathit{t}\ge 0$. At $\mathit{t}=0$, the solution satisfies the initial condition

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

Also, at $\mathit{x}=0$ and $\mathit{x}=1$, the solution satisfies the boundary conditions

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

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

To solve this equation in MATLAB®, you need to code the equation, the initial conditions, and the 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 done here), 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 rewrite it in a form that the `pdepe` solver expects. The standard form that `pdepe` expects is

$\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)$.

Written in this form, the PDE becomes

${\pi }^{2}\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$.

With the equation in the proper form you can read off the relevant terms:

`$\mathit{m}=0$`

`$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)={\pi }^{2}$`

`$\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\frac{\partial \mathit{u}}{\partial \mathit{x}}$`

`$\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0$`

Now you can create a function to code the equation. The function should have the signature `[c,f,s] = pdex1pde(x,t,u,dudx)`:

• `x` is the independent spatial variable.

• `t` is the independent time variable.

• `u` is the dependent variable being differentiated with respect to `x` and `t`.

• `dudx` is the partial spatial derivative $\partial \mathit{u}/\partial \mathit{x}$.

• The outputs `c`, `f`, and `s` correspond to coefficients in the standard PDE equation form expected by `pdepe`. These coefficients are coded in terms of the input variables `x`, `t`, `u`, and `dudx`.

As a result, the equation in this example can be represented by the function:

```function [c,f,s] = pdex1pde(x,t,u,dudx) c = pi^2; f = dudx; s = 0; end ```

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

Code Initial Condition

Next, write a function that returns the initial condition. The initial condition is applied at the first time value `tspan(1)`. The function should have the signature `u0 = pdex1ic(x)`.

The corresponding function is

```function u0 = pdex1ic(x) u0 = sin(pi*x); end ```

Code Boundary Conditions

Now, write a function that evaluates the boundary conditions. For problems posed on the interval $\mathit{a}\le \mathit{x}\le \mathit{b}$, the boundary conditions apply for all $\mathit{t}$ and either $\mathit{x}=\mathit{a}$ or $\mathit{x}=\mathit{b}$. The standard form for the boundary conditions expected by the 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.$`

Rewrite the boundary conditions in this standard form and read off the coefficient values.

For $\mathit{x}=0$, the equation is

`$\mathit{u}\left(0,\mathit{t}\right)=0\text{\hspace{0.17em}}\to \text{\hspace{0.17em}}\mathit{u}+0\cdot \frac{\partial \mathit{u}}{\partial \mathit{x}}=0.$`

The coefficients are:

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

• $\mathit{q}\left(0,\mathit{t}\right)=0$

For $\mathit{x}=1$, the equation is

`$\pi {\mathit{e}}^{-\mathit{t}}+\frac{\partial \mathit{u}}{\partial \mathit{x}}\left(1,\mathit{t}\right)=0\text{\hspace{0.17em}}\to \pi {\mathit{e}}^{-\mathit{t}}+1\cdot \frac{\partial \mathit{u}}{\partial \mathit{x}}\left(1,\mathit{t}\right)=0.$`

The coefficients are:

• $\mathit{p}\left(1,\mathit{t},\mathit{u}\right)=\pi {\mathit{e}}^{-\mathit{t}}$

• $\mathit{q}\left(1,\mathit{t}\right)=1$

Since the boundary condition function is expressed in terms of $\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)$, and this term is already defined in the main PDE function, you do not need to specify this piece of the equation in the boundary condition function. You need only specify the values of $\mathit{p}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and $\mathit{q}\left(\mathit{x},\mathit{t}\right)$ at each boundary.

The boundary function should use the function signature `[pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)`:

• The inputs `xl` and `ul` correspond to $\mathit{u}$ and $\mathit{x}$ for the left boundary.

• The inputs `xr` and `ur` correspond to $\mathit{u}$ and $\mathit{x}$ for the right boundary.

• `t` is the independent time variable.

• The outputs `pl` and `ql` correspond to $\mathit{p}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and $\mathit{q}\left(\mathit{x},\mathit{t}\right)$ for the left boundary ($\mathit{x}=0$ for this problem).

• The outputs `pr` and `qr` correspond to $\mathit{p}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and $\mathit{q}\left(\mathit{x},\mathit{t}\right)$ for the right boundary ($\mathit{x}=1$ for this problem).

The boundary conditions in this example are represented by the function:

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

Select Solution Mesh

Before solving the equation you need to specify the mesh points $\left(\mathit{t},\mathit{x}\right)$ 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. In particular, the cost and 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`.

For this problem, use a mesh with 20 equally spaced points in 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);```

Solve Equation

Finally, solve the equation using the symmetry `m`, the PDE equation, the initial conditions, the boundary conditions, and the meshes for `x` and `t`.

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

`pdepe` returns the solution in a 3-D array `sol`, where `sol(i,j,k)` approximates the `k`th component of the solution ${\mathit{u}}_{\mathit{k}}$ evaluated at `t(i)` and `x(j)`. The size of `sol` is `length(t)`-by-`length(x)`-by-`length(u0)`, since `u0` specifies an initial condition for each solution component. For this problem, `u` has only one component, so `sol` is a 5-by-20 matrix, but in general you can extract the `k`th solution component with the command `u = sol(:,:,k)`.

Extract the first solution component from `sol`.

`u = sol(:,:,1);`

Plot Solution

Create a surface plot of the solution.

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

The initial condition and boundary conditions of this problem were chosen so that there would be an analytical solution, given by

`$\mathit{u}\left(\mathit{x},\mathit{t}\right)={\mathit{e}}^{-\mathit{t}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{sin}\left(\pi \mathit{x}\right).$`

Plot the analytical solution with the same mesh points.

```surf(x,t,exp(-t)'*sin(pi*x)) title('True solution plotted with 20 mesh points') xlabel('Distance x') ylabel('Time t')```

Now, compare the numerical and analytical solutions at ${\mathit{t}}_{\mathit{f}}$, the final value of $\mathit{t}$. In this example ${\mathit{t}}_{\mathit{f}}=2$.

```plot(x,u(end,:),'o',x,exp(-t(end))*sin(pi*x)) title('Solution at t = 2') legend('Numerical, 20 mesh points','Analytical','Location','South') xlabel('Distance x') ylabel('u(x,2)')```

Local Functions

Listed here are the local helper functions that the PDE solver pdepe calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

```function [c,f,s] = pdex1pde(x,t,u,dudx) % Equation to solve c = pi^2; f = dudx; s = 0; end %---------------------------------------------- function u0 = pdex1ic(x) % Initial conditions u0 = sin(pi*x); end %---------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) % Boundary conditions pl = ul; ql = 0; pr = pi * exp(-t); qr = 1; end %----------------------------------------------```