pdepe
Solve 1-D parabolic and elliptic PDEs
Syntax
Description
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 sol
= pdepe(m
,pdefun
,icfun
,bcfun
,xmesh
,tspan
)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
.
[
also finds where functions of (t,u(x,t)), called event functions, are zero. In the output, sol
,tsol
,sole
,te
,ie
] = pdepe(m
,pdefun
,icfun
,bcfun
,xmesh
,tspan
,options
)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
The equation is defined for at times . The initial condition is defined in terms of the bessel function and its first zero as
Since this problem is in cylindrical coordinates (m = 1
), pdepe
automatically enforces the symmetry condition at . The right boundary condition is
The initial and boundary conditions are selected to be consistent with the analytic solution to the problem,
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
Written in this form, the PDE becomes
With the equation in the proper form you can read off the relevant terms:
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 tox
andt
.dudx
is the partial spatial derivative .The outputs
c
,f
, ands
correspond to coefficients in the standard PDE equation form expected bypdepe
. These coefficients are coded in terms of the input variablesx
,t
,u
, anddudx
.
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 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
Since this problem is in cylindrical coordinates (m = 1
), pdepe
automatically enforces the symmetry condition at , so you do not need to specify a left boundary condition.
For problems posed on the interval , the boundary conditions apply for all and either or . The standard form for the boundary conditions expected by the solver is
Written in this form, the boundary conditions for the partial derivatives of need to be expressed in terms of the flux . So the right boundary condition for this problem is
The boundary function should have the function signature [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
:
The inputs
xl
andul
correspond to and for the left boundary.The inputs
xr
andur
correspond to and for the right boundary.t
is the independent time variable.The outputs
pl
andql
correspond to and for the left boundary ( for this problem).The outputs
pr
andqr
correspond to and for the right boundary ( 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 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 k
th component of the solution 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 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. 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 () 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
The equation is defined for and . The initial condition is
The boundary conditions are
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
The PDE equation is already in this form:
So you can read off the relevant terms:
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 tox
andt
.dudx
is the partial spatial derivative .The outputs
c
,f
, ands
correspond to coefficients in the standard PDE equation form expected bypdepe
. These coefficients are coded in terms of the input variablesx
,t
,u
, anddudx
.
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 is
function u0 = oscic(x) u0 = 1; end
Code Boundary Conditions
Now, write a function that evaluates the boundary conditions
For problems posed on the interval , the boundary conditions apply for all and either or . The standard form for the boundary conditions expected by the solver is
Written in this form, the boundary conditions for this problem become
The boundary function should have the function signature [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)
:
The inputs
xl
andul
correspond to and for the left boundary.The inputs
xr
andur
correspond to and for the right boundary.t
is the independent time variable.The outputs
pl
andql
correspond to and for the left boundary ( for this problem).The outputs
pr
andqr
correspond to and for the right boundary ( 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 topdepe
.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 solutionumesh
. Whenvalue
equals 0, an event occurs.isterminal
specifies whether an event stops the integration. Ifisterminal
is 0, events are logged but the integration does not stop. Ifisterminal
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 at which you want pdepe
to evaluate the solution. For this problem, use a fine mesh of 50 points in the intervals and . 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 bypdepe
.tsol
is a vector of times before a terminal event.tsol
is equal tot
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. Sincevalues = umesh
in the event function,ie
gives the indices ofumesh
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 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 |
---|---|---|
| 1-D Cartesian coordinates with no symmetry | |
| 2-D Cylindrical coordinates with azimuthal angular symmetry | , where (azimuthal symmetry). |
| 3-D Spherical coordinates with azimuthal and zenith angular symmetries | , where (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
The terms in the equation are:
is a flux term.
is a source term.
The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix . 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 t0 ≤ t ≤
tf and a ≤ x ≤
b. The values tspan(1)
and
tspan(end)
correspond to
t0 and
tf, 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 can be rewritten in the form expected by the solver as
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 = t0 and all x, the solution components satisfy initial conditions of the form
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 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
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 boundaryxl = a
, andur
is the approximate solution at the right boundaryxr = b
.pl
andql
are column vectors corresponding to p and q evaluated atxl
. Similarly,pr
andqr
correspond toxr
.The number of elements in the outputs
pl
,ql
,pr
, andqr
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 on the interval . In the form expected by the solver, the boundary conditions are
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
— Option structure
structure array
Option structure, specified as a structure array. Use the odeset
function to create or modify the option structure.
pdepe
supports these options:
In most cases, default values for these options provide satisfactory solutions.
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.
ui = ui
= sol
(:
,:
,i
)
is an approximation to the i
th component of the solution vector
u. The element
ui
(j
,k
) =
sol
(j
,k
,i
)
approximates ui 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
∂ui/∂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 componenti
of the solution at timetspan(j)
and mesh pointsxmesh
, thenpdeval
evaluates the approximation and its partial derivative ∂ui/∂x at the array of pointsxout
and returns them inuout
andduoutdx
:[uout,duoutdx] = pdeval(m,xmesh,uji,xout)
. Thepdeval
function evaluates the partial derivative ∂ui/∂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
Open Example
You have a modified version of this example. Do you want to open this example with your edits?
MATLAB Command
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)