Note: This page has been translated by MathWorks. Please click here

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

Solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D

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

| A parameter corresponding to the symmetry of the problem. |

| A handle to a function that defines the components of the PDE. |

| A handle to a function that defines the initial conditions. |

| A handle to a function that defines the boundary conditions. |

| A vector [ |

| A vector [ |

| Some options of the underlying ODE solver are available
in |

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

solves
initial-boundary value problems for systems of parabolic and elliptic
PDEs in the one space variable * x* and time

`pdefun`

, `icfun`

,
and `bcfun`

are function handles. See Create Function Handle for
more information. The ordinary differential equations (ODEs) resulting
from discretization in space are integrated to obtain approximate
solutions at times specified in `tspan`

. The `pdepe`

function
returns values of the solution on a mesh provided in `xmesh`

.Parameterizing Functions explains how to provide additional
parameters to the functions `pdefun`

, `icfun`

,
or `bcfun`

, if necessary.

`pdepe`

solves 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)$$ | (1-3) |

The PDEs hold for *t*_{0} ≤ * t* ≤

In Equation 1-3, * f* (

For * t* =

$$u(x,{t}_{0})={u}_{0}(x)$$ | (1-4) |

For all * t* and either

$$p(x,t,u)+q(x,t)f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0$$ | (1-5) |

Elements of * q* are either identically zero
or never zero. Note that the boundary conditions are expressed in
terms of the flux

In the call `sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)`

:

`m`

corresponds to.*m*`xmesh(1)`

and`xmesh(end)`

correspond toand*a*.*b*`tspan(1)`

and`tspan(end)`

correspond to*t*_{0}and.*t*_{f}`pdefun`

computes the terms,*c*, and*f*(Equation 1-3). It has the form*s*[c,f,s] = pdefun(x,t,u,dudx)

The input arguments are scalars

`x`

and`t`

and vectors`u`

and`dudx`

that approximate the solutionand its partial derivative with respect to*u*, respectively.*x*`c`

,`f`

, and`s`

are column vectors.`c`

stores the diagonal elements of the matrix(Equation 1-3).*c*`icfun`

evaluates the initial conditions. It has the formu = 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`

evaluates the termsand*p*of the boundary conditions (Equation 1-5). It has the form*q*[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

`ul`

is the approximate solution at the left boundary`xl`

=and*a*`ur`

is the approximate solution at the right boundary`xr`

=.*b*`pl`

and`ql`

are column vectors corresponding toand*p*evaluated at*q*`xl`

, similarly`pr`

and`qr`

correspond to`xr`

. When> 0 and*m*= 0, boundedness of the solution near*a*= 0 requires that the flux*x*vanish at*f*= 0.*a*`pdepe`

imposes this boundary condition automatically and it ignores values returned in`pl`

and`ql`

.

`pdepe`

returns the solution as a multidimensional
array `sol`

. * u_{i}* =

`ui`

= `sol`

(`:`

,`:`

,`i`

)
is an approximation to the `i`

th component of the
solution vector `ui`

(`j`

,`k`

)
= `sol`

(`j`

,`k`

,`i`

)
approximates `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}*/∂

`xmesh`

. See `pdeval`

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

solves
as above with default integration parameters replaced by values in `options`

,
an argument created with the `odeset`

function.
Only some of the options of the underlying ODE solver are available
in `pdepe`

: `RelTol`

, `AbsTol`

, `NormControl`

, `InitialStep`

,
and `MaxStep`

. The defaults obtained by leaving off
the input argument `options`

will generally be satisfactory.
See `odeset`

for details.

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

with
the `'Events'`

property in `options`

set
to a function handle `Events`

, solves as above while
also finding where event functions `g(t,u(x,t))`

are
zero. For each function you specify whether the integration is to
terminate at a zero and whether the direction of the zero crossing
matters. Three column vectors are returned by `events`

:
`[value,isterminal,direction] = events(m,t,xmesh,umesh)`

. `xmesh`

contains
the spatial mesh and `umesh`

is the solution at the
mesh points. Use `pdeval`

to evaluate the solution
between mesh points. For the I-th event function, `value(i)`

is
the value of the function, `ISTERMINAL(I) = 1`

if
the integration is to terminate at a zero of this event function and
0 otherwise. `direction(i) = 0`

if all zeros are
to be computed (the default), +1 if only zeros where the event function
is increasing, and -1 if only zeros where the event function is
decreasing. Output `tsol`

is a column vector of
times specified in `tspan`

, prior to first terminal
event. `SOL(j,:,:)`

is the solution at `T(j)`

. `TE`

is
a vector of times at which events occur. `SOLE(j,:,:)`

is
the solution at `TE(j)`

and indices in vector `IE`

specify
which event occurred.

If `UI = SOL(j,:,i)`

approximates component
i of the solution at time `TSPAN(j)`

and mesh points `XMESH`

, `pdeval`

evaluates the approximation and
its partial derivative ∂* u_{i}*/∂

`XOUT`

and returns them in `UOUT`

and `DUOUTDX: [UOUT,DUOUTDX] = PDEVAL(M,XMESH,UI,XOUT)`

is
evaluated here rather than the flux. The flux is continuous, but
at a material interface the partial derivative may have a jump.x |

**Example 1.** This example illustrates
the straightforward formulation, computation, and plotting of the
solution of a single PDE.

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

This equation holds on an interval 0 ≤ * x* ≤
1 for times

The PDE satisfies the initial condition

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

and boundary conditions

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

It is convenient to use local functions to place all the functions
required by `pdepe`

in a single function.

function pdex11 m = 0; x = linspace(0,1,20); t = linspace(0,2,5); sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % Extract the first solution component as u. u = sol(:,:,1); % A surface plot is often a good way to study a solution. surf(x,t,u) title('Numerical solution computed with 20 mesh points.') xlabel('Distance x') ylabel('Time t') % A solution profile can also be illuminating. figure plot(x,u(end,:)) title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)') % -------------------------------------------------------------- function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = 0; % -------------------------------------------------------------- function u0 = pdex1ic(x) u0 = sin(pi*x); % -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = pi * exp(-t); qr = 1;

In this example, the PDE, initial condition, and boundary conditions
are coded in local functions `pdex1pde`

, `pdex1ic`

,
and `pdex1bc`

.

The surface plot shows the behavior of the solution.

The following plot shows the solution profile at the final value
of `t`

(i.e., `t = 2`

).

**Example 2.** This example illustrates
the solution of a system of PDEs. The problem has boundary layers
at both ends of the interval. The solution changes rapidly for small * t*.

The PDEs are

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

where * F*(

This equation holds on an interval 0 ≤ * x* ≤
1 for times

The PDE satisfies the initial conditions

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

and boundary conditions

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

In the form expected by `pdepe`

, the equations
are

$$\left[\begin{array}{l}1\\ 1\end{array}\right]\cdot *\frac{\partial}{\partial t}\left[\begin{array}{l}{u}_{1}\\ {u}_{2}\end{array}\right]=\frac{\partial}{\partial x}\left[\begin{array}{l}0.024(\partial {u}_{1}\text{/}\partial x)\\ 0.170(\partial {u}_{2}\text{/}\partial x)\end{array}\right]+\left[\begin{array}{r}-F({u}_{1}-{u}_{2})\\ F({u}_{1}-{u}_{2})\end{array}\right]$$

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}{l}1\\ 0\end{array}\right]\cdot *\left[\begin{array}{l}0.024(\partial {u}_{1}\text{/}\partial x)\\ 0.170(\partial {u}_{2}\text{/}\partial x)\end{array}\right]=\left[\begin{array}{l}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}{l}0\\ 1\end{array}\right]\cdot *\left[\begin{array}{l}0.024(\partial {u}_{1}\text{/}\partial x)\\ 0.170(\partial {u}_{2}\text{/}\partial x)\end{array}\right]=\left[\begin{array}{l}0\\ 0\end{array}\right]$$

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, the example must select the
output times accordingly. There are boundary layers in the solution
at both ends of [0,1], so the example places mesh points near

`0`

and `1`

to
resolve these sharp changes. Often some experimentation is needed
to select a mesh that reveals the behavior of the solution.function pdex4 m = 0; 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]; sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t); 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') % -------------------------------------------------------------- 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]; % -------------------------------------------------------------- function u0 = pdex4ic(x); u0 = [1; 0]; % -------------------------------------------------------------- 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];

In this example, the PDEs, initial conditions, and boundary
conditions are coded in local functions `pdex4pde`

, `pdex4ic`

,
and `pdex4bc`

.

The surface plots show the behavior of the solution components.

The arrays

`xmesh`

and`tspan`

play different roles in`pdepe`

.– The`tspan`

`pdepe`

function performs the time integration with an ODE solver that selects both the time step and formula dynamically. The elements of`tspan`

merely specify where you want answers and the cost depends weakly on the length of`tspan`

.– Second order approximations to the solution are made on the mesh specified in`xmesh`

`xmesh`

. Generally, it is best to use closely spaced mesh points where the solution changes rapidly.`pdepe`

does*not*select the mesh inautomatically. You must provide an appropriate fixed mesh in*x*`xmesh`

. The cost depends strongly on the length of`xmesh`

. When> 0, it is not necessary to use a fine mesh near*m*= 0 to account for the coordinate singularity.*x*The time integration is done with

`ode15s`

.`pdepe`

exploits the capabilities of`ode15s`

for solving the differential-algebraic equations that arise when Equation 1-3 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.

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

Was this topic helpful?