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