The PDE Toolbox is a tool to solve partial differential equations (PDE)
by making it easy to input the 2-D domain, specify the PDE coefficients
and boundary conditions, and numerically solve a finite element
discretization using piecewise linear elements. Problems can be completely
specified and solved within a graphical user interface (GUI) called
* pdetool* or the GUI can be used to specify only some of the data
such as the domain, boundary conditions, and mesh description. These
can then be exported to the main MATLAB workspace for use with user-defined
numerical algorithms.

From the **Options** menu, select *Axes Limits...* to set the
displayed (x,y) area (the default is x=-1.5,1.5 and y=-1,1), select
*Grid* to place a grid on the displayed area, and *Snap* to
automatically snap to grid points (this makes it easier to draw simple
regions). Select *Grid Spacing* if you wish to
change the grid spacing (the default is 0.5 in x and 0.2 in y).

Domains are constructed by the addition and subtraction of primitive domains such as rectangles, polygons, and ellipses. To draw a circle, whose center will be at the initial mouse position, click on the button with the ellipse with "+", then drag the mouse with held right mouse button until the desired radius is obtained. A similar procedure is used for rectangles. Clicking on the rectangle button without the "+" places the left corner of the rectangle at the intial position of the mouse. To draw a polygon, click on the polygon button and then form each side by dragging the mouse with held left button.

To form the desired domain, type in a formula in the Set formula box (above the displayed region) to combine the primitives by addition or subtraction.

To save this domain for later use, select
*Export geometry description* from the **Draw** menu and
click OK in the box. This exports the geometry data, set formulas, and labels
with the names *gd sf ns* (unless you change them). As we shall
see later, this form of the geometry (called the Constructive Solid Geometry
model) must be processed (into the "decomposed geometry") before it can
be used in several MATLAB commands. This may be done after exiting
*pdetool* in the main MATLAB work space by typing
* g = decsg(gd,sf,ns);*. However, it is easier to
wait until the boundary conditions are specified and then export
both the boundary conditions and decomposed geometry simultaneously.

Enter boundary mode by clicking on the "d Omega" button or selecting
*Boundary Mode* from the **Boundary** menu. Select one boundary
segment by clicking on it, or several boundary segments by holding down the
*Shift* key and clicking on the desired segments, or all the boundary
segments by selecting *Select All* from the **Edit** menu. To
input the boundary conditions, for this segment or group of segments, select
*Specify Boundary Conditions...* from the **Boundary** menu.
Dirichlet boundary conditions of the type h*u=r for the segments chosen are
specified by inputting the values of h and r. Note the default is h=1 and
r=0. Neumann boundary conditions are input by first clicking on the Neumann
box and then specifying the coefficients g and q. The general form of the
boundary condition appears at the top of the box (n is the unit outward normal
and c is a coefficient in the PDE).

To save the boundary conditions in a form that can be used in a MATLAB
program, select *Export Decomposed Geometry, Boundary Conds* from
the **Boundary** menu and click OK in the box. This exports the decomposed
geometry data and boundary data with the names *g b*
(unless you change them).

From the **Mesh** menu, select *Parameters* to (optionally) choose
a desired maximum edge size, where `Inf` gives the coarsest possible
mesh. To generate an initial mesh, click on the triangle button
or select *Initialize Mesh* from the **Mesh** menu.

To save the mesh in a form that can be used in a MATLAB
program, select *Export Mesh* from the **Mesh** menu and click
OK in the box. This exports the triangle vertices,
triangle edges, and triangle ordering with the names: *p e t*
(unless you change them). Note that if you wish to compare the computed
and exact solution at the vertices or compute various norms of the error,
this information will be needed.

Click on the subdivided triangle or select *Refine Mesh* from
the **Mesh** menu to get a refined mesh. The default refinement method
is regular. To change to longest edge, select *Parameters* and
then *longest* from the **Mesh** menu. The numerical solution
of the PDE can be done adaptively by selecting *Solve Parameters*
and then *Adaptive Mode* from the **Solve** menu.
Enter the desired maximum number of triangles or maximum number
of refinement steps.

Select *Solve PDE* from the **Solve** menu. To output the
solution to the main MATLAB workspace, select *Export Solution*.
Various types of plots are available. Select *Parameters*
from the **Plot** menu.

If you have solved the boundary value problem using the GUI and
have exported the data described above, you can now process this
data in the main MATLAB workspace by exiting the GUI and returning
to the MATLAB prompt >>. First, you may wish to save the decomposed
geometry *g* in a file in your directory named, for example,
*prob1g.m* so that it may be input to a future problem. This is done
by typing fid = wgeom(g, 'prob1g'); To save the boundary data *b* in a
file in your directory named, for example, *prob1b.m* so that it may be
input to a future problem, type fid = wbound(b, 'prob1b');

In fact, with these two files in your directory, it is fairly easy
to solve a simple boundary value problem without using the GUI.
To do so, first place a mesh on the geometry defined in the file
prob1g.m by typing *[p,e,t] = initmesh('prob1g');* The initmesh
command also allows the specification of certain parameters. For example,
*[p,e,t] = initmesh('prob1g', 'Hmax',1);* specifies that the
maximum edge size is 1. If you replace 1 by inf, you get the coarsest
mesh. This mesh may be viewed by typing *pdemesh(p,e,t)*.
To get a finer mesh, type *[p,e,t] = refinemesh('prob1g',p,e,t);*
The coefficients of the partial differential equation can then
be specified directly. For example, the generic scalar equation has the form
-div(c grad u) + a u = f. Hence to solve the problem with c=1, a=0,
and f = 8-16(x^{2} + y^{2}), type
c= '1'; a= '0'; f = '8 - 16.*(x.^2 + y.^2)';
To obtain the approximate solution of the boundary value problem
with boundary conditions specified in the file prob1b.m, type
*u = assempde('prob1b',p,e,t,c,a,f);* You may then view a plot of the
solution u by typing *pdesurf(p,t,u)*

There are several MATLAB commands which are useful in obtaining information
about the approximate solution and in computing errors (for model problems)
where the exact solution is known. Typing
*[ux,uy] = pdegrad(p,t,u);* produces the gradient of u on each triangle
(the first component is in the vector ux and the second component in uy).
Typing *[area,a1,a2,a3]=pdetrg(p,t);* produces a vector (area) containing
the area of each triangle along with three other geometric quantities.
Note that the L^{2} norm of the gradient of u can then be
computed by typing
*gradnorm = sqrt(sum(area.*(ux.^2 + uy.^2)));*
The values of u can also be obtained at non-vertex points by using the
function *tri2grid*. For example, typing
*v = tri2grid(p,t,u,x,y)* computes the values of u over the grid defined
by the vectors x and y (assuming u is already defined at the vertices).