The elliptic solver allows other types of equations to be more
easily implemented. In this section, we show how the parabolic equation
can be reduced to solving elliptic equations. This is done using the
function `parabolic`

.

Consider the equation

$$d\frac{\partial u}{\partial t}-\nabla \cdot \left(c\nabla u\right)+au=f\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{in}\Omega ,$$

with the initial condition

*u*(**x**,0)
= *u*_{0}(**x**)
for **x**∊Ω

where **x** represents a 2-D or
3-D point and there are boundary conditions of the same kind as for
the elliptic equation on ∂Ω.

The heat equation reads

$$\rho C\frac{\partial u}{\partial t}-\nabla \text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(k\nabla u\right)+h\left(u-{u}_{\infty}\right)=f$$

in the presence of distributed heat loss to the surroundings. *ρ* is
the density, *C* is the thermal capacity, *k* is
the thermal conductivity, *h* is the film coefficient, *u*_{∞} is
the ambient temperature, and *f* is the heat source.

For time-independent coefficients, the steady-state solution of the equation is the solution to the standard elliptic equation

–∇ · (*c*∇*u*)
+ *au* = *f*.

Assuming a mesh on Ω and *t* ≥
0, expand the solution to the PDE (as a function of **x**)
in the Finite Element Method basis:

$$u(x,t)={\displaystyle \sum _{i}{U}_{i}(t){\varphi}_{i}(x)}.$$

Plugging the expansion into the PDE, multiplying with a test
function *ϕ _{j}*, integrating
over Ω, and applying Green's formula and the boundary conditions
yield

$$\begin{array}{l}{\displaystyle \sum _{i}{\displaystyle \underset{\Omega}{\int}d{\varphi}_{j}{\varphi}_{i}\text{\hspace{0.17em}}\frac{d{U}_{i}\left(t\right)}{dt}}}\text{\hspace{0.17em}}dx+{\displaystyle \sum _{i}\left({\displaystyle \underset{\Omega}{\int}\left(\nabla {\varphi}_{j}\cdot \left(c\nabla {\varphi}_{i}\right)+a{\varphi}_{j}{\varphi}_{i}\right)\text{\hspace{0.17em}}dx+{\displaystyle \underset{\partial \Omega}{\int}q{\varphi}_{j}{\varphi}_{i}\text{\hspace{0.17em}}ds}}\right){U}_{i}(t)}\\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}={\displaystyle \underset{\Omega}{\int}f{\varphi}_{j}\text{\hspace{0.17em}}dx}+{\displaystyle \underset{\partial \Omega}{\int}g{\varphi}_{j}\text{\hspace{0.17em}}ds}\text{\hspace{1em}}\forall j.\end{array}$$

In matrix notation, we have to solve the *linear, large* and *sparse* ODE
system

$$M\frac{dU}{dt}+KU=F.$$

This method is traditionally called *method of lines* semidiscretization.

Solving the ODE with the initial value

*U _{i}*(0) =

yields the solution to the PDE at each node **x*** _{i}* and
time

–∇ · (*c*∇*u*)
+ *au* = *f* in Ω

with the original boundary conditions, while *M* is
just the mass matrix of the problem

–∇ · (0∇*u*)
+ *du* = 0 in Ω.

When the Dirichlet conditions are time dependent, *F * contains
contributions from time derivatives of *h* and **r**. These derivatives are evaluated by finite
differences of the user-specified data.

The ODE system is ill conditioned. Explicit time integrators
are forced by stability requirements to very short time steps while
implicit solvers can be expensive since they solve an elliptic problem
at every time step. The numerical integration of the ODE system is
performed by the MATLAB^{®} ODE Suite functions, which are efficient
for this class of problems. The time step is controlled to satisfy
a tolerance on the error, and factorizations of coefficient matrices
are performed only when necessary. When coefficients are time dependent,
the necessity of reevaluating and refactorizing the matrices each
time step may still make the solution time consuming, although `parabolic`

reevaluates
only that which varies with time. In certain cases a time-dependent
Dirichlet matrix **h**(*t*)
may cause the error control to fail, even if the problem is mathematically
sound and the solution *u*(*t*)
is smooth. This can happen because the ODE integrator looks only at
the reduced solution *v* with *u* = *Bv* + *ud*.
As **h** changes, the pivoting scheme
employed for numerical stability may change the elimination order
from one step to the next. This means that *B*, *v*,
and *ud* all change discontinuously, although *u* itself
does not.

This example shows how to solve a parabolic equation and to set an initial condition as a variable.

At the MATLAB command prompt, type

`pdetool`

.Draw a rectangle in the PDE app axes.

From the

**Draw**menu, select**Export Geometry Description, Set Formula, Labels**.In the Export dialog box, enter

`gd sf ns`

. Click**OK**.The exported variables are available in the MATLAB workspace.

From the

**Boundary**menu, select**Boundary Mode**.From the

**Boundary**menu, select**Specify Boundary Conditions**.Set the Neumann and Dirichlet boundary conditions. If these conditions are not the same for all the stages, set the conditions accordingly.

From the

**Boundary**menu, select**Export Decomposed Geometry, Boundary Cond's**.In the Export dialog box, enter

`g b`

. Click**OK**.From the

**PDE**menu, select**PDE Mode**.From the

**PDE**menu, select**PDE Specification**.Set the partial differential equation (PDE) coefficients, which are the same for any value of time.

From the

**PDE**menu, select**Export PDE Coefficients.**In the Export dialog box, enter

`c a f d`

. Click**OK**.From the

**Mesh**menu, select**Mesh Mode**.From the

**Mesh**menu, select**Parameters**.Verify the initial mesh, jiggle mesh, and refine mesh values. The mesh is fixed for all stages.

From the

**Mesh**menu, select**Export Mesh**.In the Export dialog box, enter

`p e t`

. Click**OK**.Save the workspace variables into a MAT-file by typing

`save data.mat`

at the MATLAB command prompt.Save the following code as a file:

clear all; close all; load data %For the first stage you need to specify an %initial condition, U0. U0 = 0; %U0 expands to the correct size automatically. %Divide the time range into 4 stages. time = {0:.01:1, 1:.05:3, 3:.1:5, 5:.5:20}; for i = 1:4 U1 = parabolic(U0,time{i},b,p,e,t,c,a,f,d); for j = 1:size(U1,2) H =pdeplot(p,e,t,'xydata',U1(:,j),'zdata',... U1(:,j),'mesh','off'); set(gca,'ZLim',[-80 0]); drawnow end %Reset the initial condition at all points. U0 = U1(:,1); end

This file uses the variables you defined in the MATLAB workspace to solve a parabolic equation in stages. Within this file, you set the initial condition as a variable.

Was this topic helpful?