pdepe

Solve 1-D parabolic and elliptic PDEs

Syntax

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

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

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

Description

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) solves a system of parabolic and elliptic PDEs with one spatial variable x and time t. At least one equation must be parabolic. The scalar m represents the symmetry of the problem (slab, cylindrical, or spherical). The equations being solved are coded in pdefun, the initial value is coded in icfun, and the boundary conditions are coded in bcfun. The ordinary differential equations (ODEs) resulting from discretization in space are integrated to obtain approximate solutions at the times specified in tspan. The pdepe function returns values of the solution on a mesh provided in xmesh.

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) also uses the integration settings defined by options, which is created using the odeset function. For example, use the AbsTol and RelTol options to specify absolute and relative error tolerances.

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) also finds where functions of (t,u(x,t)), called event functions, are zero. In the output, te is the time of the event, sole is the solution at the time of the event, and ie is the index of the triggered event. tsol is a column vector of times specified in tspan, prior to the first terminal event.

For each event function, specify whether the integration is to terminate at a zero and whether the direction of the zero-crossing matters. Do this by setting the Events option of odeset to a function, such as @myEventFcn, and creating a corresponding function: [value,isterminal,direction] = myEventFcn(m,t,xmesh,umesh). The xmesh input contains the spatial mesh and umesh is the solution at the mesh points.

Examples

Solve Heat Equation in Cylindrical Coordinates

Solve the heat equation in cylindrical coordinates using pdepe, and plot the solution.

In cylindrical coordinates with angular symmetry the heat equation is

$$\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{1}{\mathit{x}}\frac{\partial }{\partial \mathit{x}}\left(\mathit{x}\frac{\partial \mathit{u}}{\partial \mathit{x}}\right).$$

The equation is defined for $0\le \mathit{x}\le 1$ at times $\mathit{t}\ge 0$. The initial condition is defined in terms of the bessel function ${\mathit{J}}_0\left(\mathit{x}\right)$ and its first zero $\mathit{n}=2.404825557695773$ as

$$\mathit{u}\left(\mathit{x},0\right)={\mathit{J}}_0\left(\mathrm{nx}\right).$$

Since this problem is in cylindrical coordinates (m = 1), pdepe automatically enforces the symmetry condition at $\mathit{x}=0$. The right boundary condition is

$$\mathit{u}\left(1,\mathit{t}\right)={\mathit{J}}_0\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^2\mathit{t}}.$$

The initial and boundary conditions are selected to be consistent with the analytic solution to the problem,

$$\mathit{u}\left(\mathit{x},\mathit{t}\right)={\mathit{J}}_0\left(\mathrm{nx}\right){\mathit{e}}^{-{\mathit{n}}^2\mathit{t}}.$$

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 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

$$\frac{\partial \mathit{u}}{\partial \mathit{t}}={\mathit{x}}^{-1}\frac{\partial }{\partial \mathit{x}}\left({\mathit{x}}^1\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)$$

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

• $\mathit{m}=1$

• $\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=1$

• $\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] = heatcyl(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] = heatcyl(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 Conditions

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 = heatic(x).

The corresponding function for $\mathit{u}\left(\mathit{x},0\right)={\mathit{J}}_0\left(\mathrm{nx}\right)$ is

function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end

Code Boundary Conditions

Now, write a function that evaluates the boundary condition

$\mathit{u}\left(1,\mathit{t}\right)={\mathit{J}}_0\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^2\mathit{t}}.$

Since this problem is in cylindrical coordinates (m = 1), pdepe automatically enforces the symmetry condition at $\mathit{x}=0$, so you do not need to specify a left boundary condition.

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

Written in this form, the boundary conditions for the partial derivatives of $\mathit{u}$ need to be expressed in terms of the flux $\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)$. So the right boundary condition for this problem is

$$\mathit{u}\left(1,\mathit{t}\right)-{\mathit{J}}_0\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^2\mathit{t}}+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$$

The boundary function should have the function signature [pl,ql,pr,qr] = heatbc(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}}_{\mathit{L}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and ${\mathit{q}}_{\mathit{L}}\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}}_{\mathit{R}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and ${\mathit{q}}_{\mathit{R}}\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] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
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 25 equally spaced points in the spatial interval [0,1] for both x and t.

x = linspace(0,1,25);
t = linspace(0,1,25);

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 = 1;
sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

pdepe returns the solution in a 3-D array sol, where sol(i,j,k) approximates the kth 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 25-by-25 matrix, but in general you can extract the kth 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. Since the problem is posed with cylindrical coordinates on a disc, the x values show the temperature on the disc some distance from the center, and the t values show how the temperature at a particular location changes over time.

surf(x,t,u)
xlabel("x")
ylabel("t")
zlabel("u(x,t)")
view([150 25])

Plot the temperature change at the center ($\mathit{x}=0$) of the disc.

plot(t,sol(:,1))
xlabel("Time")
ylabel("Temperature u(0,t)")
title("Temperature change at center of disc")

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] = heatcyl(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
%----------------------------------------------
function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end
%----------------------------------------------
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
end
%----------------------------------------------

Solve Oscillatory PDE with Event Logging

Solve a partial differential equation and use an event function to log zero-crossings in the oscillatory solution.

Consider the equation

$$\frac{1}{\mathit{x}}\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{\partial }{\partial \mathit{x}}\left(\frac{1}{\mathit{t}}\mathit{u}\right).$$

The equation is defined for $0\le \mathit{x}\le 1$ and $\mathit{t}\ge \mathrm{0}.1$. The initial condition is

$$\mathit{u}\left(\mathit{x},0.1\right)=1.$$

The boundary conditions are

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

$$\mathit{u}\left(1,\mathit{t}\right)=\cos \left(\pi \mathit{t}\right).$$

Additionally, the zero-crossings of the solution are of interest.

To solve this equation in MATLAB, you need to code the equation, initial conditions, boundary conditions, and event function, 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).$$

The PDE equation is already in this form:

$\frac{1}{\mathit{x}}\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{\partial }{\partial \mathit{x}}\left(\frac{1}{\mathit{t}}\mathit{u}\right).$

So 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)=\frac{1}{\mathit{x}}$$

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

• $$\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] = oscpde(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] = oscpde(x,t,u,dudx)
c = 1/x;
f = u/t;
s = 0;
end

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

Code Initial Conditions

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 = oscic(x).

The corresponding function for $\mathit{u}\left(\mathit{x},0.1\right)=1$ is

function u0 = oscic(x)
u0 = 1;
end

Code Boundary Conditions

Now, write a function that evaluates the boundary conditions

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

$$\mathit{u}\left(1,\mathit{t}\right)=\cos \left(\pi \mathit{t}\right).$$

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

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

$$\mathit{u}\left(0,\mathit{t}\right)-1+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0,$$

$$\mathit{u}\left(1,\mathit{t}\right)-\cos \left(\pi \mathit{t}\right)+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$$

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

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

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

• t is the independent time variable.

• The outputs pl and ql correspond to ${\mathit{p}}_{\mathit{L}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and ${\mathit{q}}_{\mathit{L}}\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}}_{\mathit{R}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and ${\mathit{q}}_{\mathit{R}}\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] = oscbc(xl,ul,xr,ur,t)
pl = ul - 1;
ql = 0;
pr = ur - cos(pi*t);
qr = 0;
end

Code Event Function

Use an event function to log the zero-crossings of the solution in the integration. The event function has a function signature of [value,isterminal,direction] = pdevents(m,t,xmesh,umesh):

• m is the coordinate symmetry specified as the first input to pdepe.

• t is the current time (a scalar).

• xmesh is the spatial mesh.

• umesh contains the solution at the mesh points.

• value is an equation of interest, usually expressed in terms of the solution umesh. When value equals 0, an event occurs.

• isterminal specifies whether an event stops the integration. If isterminal is 0, events are logged but the integration does not stop. If isterminal is 1, then the integration stops when an event occurs.

• direction specifies the direction of zero crossing. If 1, only zero-crossings with a positive slope trigger an event. If -1, the zero-crossings must have a negative slope. If 0, then any zero-crossing triggers an event.

At each time in the integration, the solver calls the event function to check for zero crossings. To log all zero crossings, value should look for changes in sign in the solution vector umesh. Specify isterminal and direction as vectors of zeros with the same size, since the events in this example are not terminal and the zero-crossings can occur with any slope.

The event function for this problem is

function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
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. For this problem, use a fine mesh of 50 points in the intervals $0\le \mathit{x}\le 1$ and $0.1\le \mathit{t}\le 1$. A fine mesh gives good resolution of the oscillatory solution.

x = linspace(0,1,50);
t = linspace(0.1,pi,50);

Solve Equation

Finally, solve the equation using the symmetry m, the PDE equation, the initial conditions, the boundary conditions, the event function, and the meshes for x and t. Use odeset to create an options structure that references the events function, and pass in the structure as the last input argument to pdepe. Specify five output arguments to return information from both the event function and the solver:

• sol is the solution computed by pdepe.

• tsol is a vector of times before a terminal event. tsol is equal to t when no events are terminal.

• sole is the solution at the time of each event.

• te is the time of each event.

• ie is the index of each event. Since values = umesh in the event function, ie gives the indices of umesh that triggered an event at each time step.

m = 0;
options = odeset(Events=@pdevents);
[sol,tsol,sole,te,ie] = pdepe(m,@oscpde,@oscic,@oscbc,x,t,options);

Extract the solution as a matrix u.

u = sol(:,:,1);

Plot Solution

Create a surface plot of the solution and view the plot from above.

surf(x,t,u)
view(2)

Plot the points where events occurred, with a surface $\mathit{f}\left(\mathit{x},\mathit{t}\right)=0$ for reference. The output index vector ie is useful to pick out the event locations. The expression x(ie)' gives the x-values where events occurred, and the expression sole(x==x(ie)') gives the corresponding solution values.

view([39 30])
xlabel("x")
ylabel("t")
zlabel("u(x,t)")
hold on
plot3(x(ie)',te,sole(x==x(ie)'),'r*')
surf(x,t,zeros(size(u)),EdgeColor="flat")
hold off

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] = oscpde(x,t,u,dudx)
c = 1/x;
f = u/t;
s = 0;
end
%----------------------------------------------
function u0 = oscic(x)
u0 = 1;
end
%----------------------------------------------
function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)
pl = ul - 1;
ql = 0;
pr = ur - cos(pi*t);
qr = 0;
end
%----------------------------------------------
function [value, isterminal, direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end
%----------------------------------------------

Input Arguments

m — Symmetry constant

0 | 1 | 2

Symmetry constant, specified as one of the values in this table. m specifies the type of problem in the equations specified by pdefun. Once you rewrite the equations in the form expected by pdepe, you can read the value of m from the equation.

Value	Description	Laplacian of flux
0	1-D Cartesian coordinates with no symmetry	$$\Delta f=\frac{{\partial }^{2}f}{\partial x^2}$$
1	2-D Cylindrical coordinates with azimuthal angular symmetry	$$\Delta f=\frac{1}{\rho }\frac{\partial }{\partial \rho }\left(\rho \frac{\partial f}{\partial \rho }\right)+\frac{1}{{\rho }^2}\frac{{\partial }^{2}f}{\partial {\phi }^2}$$, where $$\frac{\partial f}{\partial \phi }=0$$ (azimuthal symmetry).
2	3-D Spherical coordinates with azimuthal and zenith angular symmetries	$$\Delta f=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial f}{\partial r}\right)+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial f}{\partial \theta }\right)+\frac{1}{r^2{\sin }^2\theta }\frac{{\partial }^{2}f}{\partial {\phi }^2}$$, where $$\frac{\partial f}{\partial \theta }=\frac{\partial f}{\partial \phi }=0$$ (azimuthal and zenith angular symmetries).

When m > 0, pdepe requires the left spatial boundary a to be ≥ 0, and it imposes the left boundary condition automatically to account for the coordinate discontinuity at the origin. In that case, pdepe ignores any specified left boundary conditions in bcfun.

pdefun — PDE function for equations

function handle

PDE function for equations, specified as a function handle that defines the coefficients of the PDE.

pdefun encodes the coefficients for 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 terms in the equation are:

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

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

• 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 all other elements correspond to parabolic equations. 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 those values of x are mesh points.

• Discontinuities in c and s due to material interfaces are permitted provided that a mesh point is placed at each interface.

The PDEs hold for t₀ ≤ t ≤ t_f and a ≤ x ≤ b. The values tspan(1) and tspan(end) correspond to t₀ and t_f, while xmesh(1) and xmesh(end) correspond to a and b. 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 must be ≥ 0.

Coding the Equation

Write a function pdefun that calculates the values of the coefficients c, f, and s. Use the function signature

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

The input arguments to pdefun are scalars x and t and vectors u and dudx. The vector u approximates the solution u, and dudx approximates its partial derivative with respect to x. The output arguments c, f, and s are column vectors with a number of elements equal to the number of equations. c stores the diagonal elements of the matrix c.

Example

The heat equation $$\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}$$ can be rewritten in the form expected by the solver as

$$1·\frac{\partial u}{\partial t}=x^0\frac{\partial }{\partial x}\left(x^0\frac{\partial u}{\partial x}\right).$$

In this form you can read off the value of the coefficients as m = 0, c = 1, f = ∂u/∂x, and s = 0. The function for this equation is:

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

Data Types: function_handle

icfun — Initial condition function

function handle

Initial condition function, specified as a function handle that defines the initial condition.

For t = t₀ and all x, the solution components satisfy initial conditions of the form

$$u\left(x,t_0\right)=u_0\left(x\right).$$

Coding the Initial Condition

Write a function icfun that defines the initial conditions. Use the function signature:

function u0 = icfun(x)

When called with an argument x, the function icfun evaluates and returns the initial values of the solution components at x in the column vector u0. The number of elements in u0 is equal to the number of equations.

Example

The constant initial condition $$u\left(x,0\right)=1/2$$ corresponds to this function:

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

Data Types: function_handle

bcfun — Boundary condition function

function handle

Boundary condition function, specified as a function handle that defines the boundary conditions.

For all t and either x = a or x = b, 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.$$

Elements of q are either zero or never zero. Note that the boundary conditions are expressed in terms of the flux f rather than ∂u/∂x. Also, of the two coefficients, only p can depend on u.

Coding the Boundary Conditions

Write a function bcfun that defines the terms p and q of the boundary conditions. Use the function signature

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

• 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 q evaluated at xl. Similarly, pr and qr correspond to xr.

• The number of elements in the outputs pl, ql, pr, and qr is equal to the number of equations.

When m > 0 and a = 0, boundedness of the solution near x = 0 requires that the flux f vanish at a = 0. pdepe imposes this boundary condition automatically, and it ignores values returned in pl and ql.

Example

Consider the boundary conditions $$u\left(0,t\right)=0,u\left(1,t\right)=1$$ on the interval $$0\le x\le 1$$. In the form expected by the solver, the boundary conditions are

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

In this form it is clear that both q terms are zero. The p terms are nonzero to reflect the values of u. The corresponding function is

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

Data Types: function_handle

xmesh — Spatial mesh

vector

Spatial mesh, specified as a vector [x0 x1 ... xn] specifying the points at which a numerical solution is requested for every value in tspan. Second-order approximations to the solution are made on the mesh specified in xmesh. pdepe does not select the mesh in x automatically. You must provide an appropriate fixed mesh in xmesh.

The elements of xmesh must satisfy x0 < x1 < ... < xn. Generally, it is best to use closely spaced mesh points where the solution changes rapidly. The length of xmesh must be >=3, and the computational cost of the solution depends strongly on the length of xmesh.

For problems with discontinuities, you should place mesh points at the discontinuities so that the problem is smooth within each subinterval. However, when m > 0 , it is not necessary to use a fine mesh near x = 0 to account for the coordinate singularity.

Example: xmesh = linspace(0,5,25) uses a spatial mesh of 25 points between 0 and 5.

Data Types: single | double

tspan — Time span of integration

vector

Time span of integration, specified as a vector [t0 t1 ... tf] specifying the points at which a solution is requested for every value in xmesh. The pdepe function performs the time integration with an ODE solver that selects both the time step and the formula dynamically. The elements of tspan merely specify where you want answers, and the computational cost of the solution therefore depends only weakly on the length of tspan.

The elements of tspan must satisfy t0 < t1 < ... < tf. The length of tspan must be >=3.

Example: tspan = linspace(0,5,5) uses five time points between 0 and 5.

Data Types: single | double

options — Options structure

structure array

Options structure, specified as a structure array. Use the odeset function to create or modify the options structure. You can set any of the options in odeset, but pdepe overrides the following options when solving the PDE:

• Mass

• MStateDependence

• MvPattern

• MassSingular

• Jacobian

• JPattern

• Vectorized

Example: options = odeset(RelTol=1e-5) specifies a relative error tolerance of 1e-5.

Data Types: struct

Output Arguments

sol — Solution array

3-D array

Solution array, returned as a 3-D array.

pdepe returns the solution as a multidimensional array. u_i = ui = sol(:,:,i) is an approximation to the ith component of the solution vector u. The element ui(j,k) = sol(j,k,i) approximates u_i at (t,x) = (tspan(j),xmesh(k)).

ui = sol(j,:,i) approximates component i of the solution at time tspan(j) and mesh points xmesh(:). Use pdeval to compute the approximation and its partial derivative ∂u_i/∂x at points not included in xmesh. See pdeval for details.

tsol — Solution times

column vector

Solution times, returned as a column vector of times specified in tspan that are prior to the first terminal event.

sole — Solution at time of events

array

Solution at time of events, returned as an array. The event times in te correspond to solutions returned in sole, and ie specifies which event occurred.

te — Time of events

column vector

Time of events, returned as a column vector. The event times in te correspond to solutions returned in sole, and ie specifies which event occurred.

ie — Index of triggered event function

column vector

Index of triggered event function, returned as a column vector. The event times in te correspond to solutions returned in sole, and ie specifies which event occurred.

Tips

• If uji = sol(j,:,i) approximates component i of the solution at time tspan(j) and mesh points xmesh, then pdeval evaluates the approximation and its partial derivative ∂u_i/∂x at the array of points xout and returns them in uout and duoutdx: [uout,duoutdx] = pdeval(m,xmesh,uji,xout). The pdeval function evaluates the partial derivative ∂u_i/∂x rather than the flux. The flux is continuous, but at a material interface the partial derivative may have a jump.

Algorithms

The time integration is done with the ode15s solver. pdepe exploits the capabilities of ode15s for solving the differential-algebraic equations that arise when the PDE contains elliptic equations, and for handling Jacobians with a specified sparsity pattern.

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.

References

[1] 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.

Extended Capabilities

Thread-Based Environment
Run code in the background using MATLAB® backgroundPool or accelerate code with Parallel Computing Toolbox™ ThreadPool.

Version History

Introduced before R2006a

R2025a: Specify output function and other odeset options

You can now specify all odeset options, including the output function, when using the pdepe function. However, pdepe overrides some options when defining the PDE to solve. See options for more details.

Previously, pdepe supported only the RelTol, AbsTol, NormControl, InitialStep, MaxStep, and Events options.