MATLAB PDE Toolbox Commands

What does the MATLAB PDE Toolbox do?

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.

Starting the PDE Toolbox

To start the PDE Toolbox graphical user interface, first start MATLAB, by typing matlab at the unix prompt and then typing pdetool when the MATLAB prompt ">>" appears.

Specifying the Domain

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.

Specifying the Boundary Conditions

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

Specifying the PDE

From the toolbar or by selecting Application from the Options menu, choose an application type (the default is a generic scalar problem of the form -div(c*grad(u))+ a*u =f). From the PDE menu, select PDE specification and then choose a problem type (the default is Elliptic). Enter the coefficients c and a of the PDE and the right hand side function f. The coefficients can be functions of x, y and for nonlinear problems also of u, ux, uy : Use x,y,u,ux,uy in your expression, noting that MATLAB thinks of these as a vector representing these quantities at a fixed set of points. Hence to express the function (xy) in this notation, use x.*y instead of x*y, since you want elementwise multiplication to produce a vector whose components are the values of (xy) at the fixed set of points.

Specifying the Initial Mesh

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.

Mesh Refinement

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.

Solving the boundary value problem and plotting the solution

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.

Working in the MATLAB workspace

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(x2 + y2), 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)

Post-Processing the Solution

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 L2 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).

Visualization Commands

If geom.m is the name of a geometry M-file, then the PDE geometry described in this file is displayed by typing (at the MATLAB prompt) pdegplot('geom'). To display the triangular mesh given by the mesh data p, e, t, type pdemesh(p,e,t) The command pdesurf(p,t,u) plots a 3-D surface using the node or triangle data of u. More general plots are possible using the pdeplot command.