# Documentation

### This is machine translation

Translated by
Mouse over text to see original. Click the button below to return to the English verison of the page.

## Parabolic Equations

### Reducing Parabolic Equations to Elliptic Equations

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

Partial Differential Equation Toolbox™ solves equations of the form

`$m\frac{{\partial }^{2}u}{\partial {t}^{2}}+d\frac{\partial u}{\partial t}-\nabla ·\left(c\nabla u\right)+au=f$`

When the m coefficient is 0, but d is not, the documentation refers to the equation as parabolic, whether or not it is mathematically in parabolic form.

A parabolic problem is to solve the equation

with the initial condition

u(x,0) = u0(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}}·\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

–∇ · (cu) + 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\left(x,t\right)=\sum _{i}{U}_{i}\left(t\right){\varphi }_{i}\left(x\right)$`

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}\sum _{i}\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+\sum _{i}\left(\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+\underset{\partial \Omega }{\int }q{\varphi }_{j}{\varphi }_{i}\text{\hspace{0.17em}}ds\right){U}_{i}\left(t\right)\\ \text{ }\text{ }\text{ }\text{ }=\underset{\Omega }{\int }f{\varphi }_{j}\text{\hspace{0.17em}}dx+\underset{\partial \Omega }{\int }g{\varphi }_{j}\text{\hspace{0.17em}}ds\text{ }\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

Ui(0) = u0(xi)

yields the solution to the PDE at each node xi and time t. Note that K and F are the stiffness matrix and the right-hand side of the elliptic problem

–∇ · (cu) + 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.

Was this topic helpful?

Get trial now