Documentation |
This topic describes the solution of some elliptic PDE problems. The last problem, a minimal surface problem, is nonlinear and illustrates the use of the nonlinear solver. The problems are solved using both the PDE app and command-line functions. The topics include:
On this page… |
---|
This example shows how to solve a simple elliptic PDE in the form of Poisson's equation on a unit disk.
The problem formulation is
–ΔU = 1 in Ω, U = 0 on ∂Ω,
where Ω is the unit disk. In this case, the exact solution is
$$U\left(x,y\right)=\frac{1-{x}^{2}-{y}^{2}}{4},$$
so the error of the numeric solution can be evaluated for different meshes.
With the PDE app started, perform the following steps using the generic scalar mode:
Using some of the Option menu features, add a grid and turn on the "snap-to-grid" feature. Draw a circle by clicking the button with the ellipse icon with the + sign, and then click-and-drag from the origin, using the right mouse button, to a point at the circle's perimeter. If the circle that you create is not a perfect unit circle, double-click the circle. This opens a dialog box where you can specify the exact center location and radius of the circle.
Enter the boundary mode by clicking the button with the ∂Ω icon. The boundaries of the decomposed geometry are plotted, and the outer boundaries are assigned a default boundary condition (Dirichlet boundary condition, u = 0 on the boundary). In this case, this is what we want. If the boundary condition is different, double-click the boundary to open a dialog box through which you can enter and display the boundary condition.
To define the partial differential equation, click the PDE button. This opens a dialog box, where you can define the PDE coefficients c, a, and f. In this simple case, they are all constants: c = 1, f = 1, and a = 0.
Click the button or select Initialize Mesh from the Mesh menu. This initializes and displays a triangular mesh.
Click the button or select Refine Mesh from the Mesh menu. This causes a refinement of the initial mesh, and the new mesh is displayed.
To solve the system, just click the = button. The toolbox assembles the PDE problem and solves the linear system. It also provides a plot of the solution. Using the Plot Selection dialog box, you can select different types of solution plots.
To compare the numerical solution to the exact solution, select the user entry in the Property pop-up menu for Color in the Plot Selection dialog box. Then input the MATLAB^{®} expression u-(1-x.^2-y.^2)/4 in the user entry edit field. You obtain a plot of the absolute error in the solution.
You can also compare the numerical solution to the exact solution by entering some simple command-line-oriented commands. It is easy to export the mesh data and the solution to the MATLAB main workspace by using the Export options from the Mesh and Solve menus. To refine the mesh and solve the PDE successively, simply click the refine and = buttons until the desired accuracy is achieved. (Another possibility is to use the adaptive solver.)
This example shows how to solve Poisson's equation using command-line functions. First you must create a function that parameterizes the 2-D geometry--in this case a unit circle. The circleg.m file returns the coordinates of points on the unit circle's boundary. The file conforms to the file format described on the reference page for pdegeom. You can display the file by typing type circleg.
Also, you need a function that describes the boundary condition. This is a Dirichlet boundary condition where u = 0 on the boundary. The circleb1.m file provides the boundary condition. The file conforms to the file format described on the reference page for pdebound. You can display the file by typing type circleb1.
Now you can start working at the command line:
[p,e,t] = initmesh('circleg','Hmax',1); % create mesh error = []; err = 1; while err > 0.001, % run until error <= 0.001 [p,e,t] = refinemesh('circleg',p,e,t); % refine mesh u = assempde('circleb1',p,e,t,1,0,1); % solve equation exact = -(p(1,:).^2+p(2,:).^2-1)/4; err = norm(u-exact',inf); % compare with exact solution error = [error err]; % keep history of err end pdesurf(p,t,u-exact') % plot error
pdedemo1 performs all the previous steps.
This example shows how to solve a simple scattering problem, where you compute the waves reflected from an object illuminated by incident waves. For this problem, assume an infinite horizontal membrane subjected to small vertical displacements U. The membrane is fixed at the object boundary.
We assume that the medium is homogeneous so that the wave speed is constant, c.
Note Do not confuse this c with the parameter c appearing in Partial Differential Equation Toolbox™ functions. |
When the illumination is harmonic in time, we can compute the field by solving a single steady problem. With
U(x,y,t) = u(x,y)e^{–iωt},
the wave equation
$$\frac{{\partial}^{2}U}{\partial {t}^{2}}-{c}^{2}\Delta U=0$$
turns into
–ω^{2}u – c^{2}Δu = 0
or the Helmholtz's equation
–Δu – k^{2}u = 0,
where k, the wave number, is related to the angular frequency ω, the frequency f, and the wavelength λ by
$$k=\frac{\omega}{c}=\frac{2\pi f}{c}=\frac{2\pi}{\lambda}.$$
We have yet to specify the boundary conditions. Let the incident wave be a plane wave traveling in the direction $$\overrightarrow{a}$$ = (cos(a), sin(a)):
$$V(x,y,t)={e}^{i\left(k\overrightarrow{a}\cdot \overrightarrow{x}-\omega t\right)}=v(x,y){e}^{-i\omega t},$$
where
$$v(x,y)={e}^{ik\overrightarrow{a}\cdot \overrightarrow{x}}.$$
u is the sum of v and the reflected wave r,
u = v + r.
The boundary condition for the object's boundary is easy: u = 0, i.e.,
r = –v(x,y)
For acoustic waves, where v is the pressure disturbance, the proper condition would be
$$\frac{\partial u}{\partial n}=0.$$
The reflected wave r travels outward from the object. The condition at the outer computational boundary should be chosen to allow waves to pass without reflection. Such conditions are usually called nonreflecting, and we use the classical Sommerfeld radiation condition. As $$\left|\overrightarrow{x}\right|$$ approaches infinity, r approximately satisfies the one-way wave equation
$$\frac{\partial r}{\partial t}+c\overrightarrow{\xi}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla r=0,$$
which allows waves moving in the positive ξ-direction only (ξ is the radial distance from the object). With the time-harmonic solution, this turns into the generalized Neumann boundary condition
$$\overrightarrow{\xi}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla r=ikr.$$
For simplicity, let us make the outward normal of the computational domain approximate the outward ξ-direction.
You can now use the PDE app to solve this scattering problem. Using the generic scalar mode, start by drawing the 2-D geometry of the problem. Let the illuminated object be a square SQ1 with a side of 0.1 units and center in [0.8 0.5] and rotated 45 degrees, and let the computational domain be a circle C1 with a radius of 0.45 units and the same center location. The Constructive Solid Geometry (CSG) model is then given by C1-SQ1.
For the outer boundary (the circle perimeter), the boundary condition is a generalized Neumann condition with q = –ik. The wave number k = 60, which corresponds to a wavelength of about 0.1 units, so enter -60i as a constant q and 0 as a constant g.
For the square object's boundary, you have a Dirichlet boundary condition:
$$r=-v\left(x,y\right)=-{e}^{ik\overrightarrow{a}\text{\hspace{0.17em}}\cdot \text{\hspace{0.17em}}\overrightarrow{x}}.$$
In this problem, the incident wave is traveling in the –x direction, so the boundary condition is simply
r = –e^{–ikx}.
Enter this boundary condition in the Boundary Condition dialog box as a Dirichlet condition: h = 1, r = -exp(-i*60*x). The real part of this is a sinusoid.
For sufficient accuracy, about 10 finite elements per wavelength are needed. The outer boundary should be located a few object diameters from the object itself. An initial mesh generation and two successive mesh refinements give approximately the desired resolution.
Although originally a wave equation, the transformation into a Helmholtz's equation makes it—in the Partial Differential Equation Toolbox context, but not strictly mathematically—an elliptic equation. The elliptic PDE coefficients for this problem are c = 1, a = -k^{2} = -3600, and f = 0. Open the PDE Specification dialog box and enter these values.
The problem can now be solved, and the solution is complex. For a complex solution, the real part is plotted and a warning message is issued.
The propagation of the reflected waves is computed as
Re(r(x,y)e^{–iωt}),
which is the reflex of
$$\mathrm{Re}\left({e}^{i\left(k\overrightarrow{a}\cdot \overrightarrow{x}-\omega t\right)}\right).$$
To see the whole field, plot
$$\mathrm{Re}\left(\left(r(x,y)+{e}^{ik\overrightarrow{a}\cdot \overrightarrow{x}}\right){e}^{-i\omega t}\right).$$
The reflected waves and the "shadow" behind the object are clearly visible when you plot the reflected wave.
To make an animation of the reflected wave, the solution and the mesh data must first be exported to the main workspace. Then make a script file or type the following commands at the MATLAB prompt:
h = newplot; hf = get(h,'Parent'); set(hf,'Renderer','zbuffer') axis tight, set(gca,'DataAspectRatio',[1 1 1]); axis off M = moviein(10,hf); maxu = max(abs(u)); colormap(cool) for j = 1:10, ur = real(exp(-j*2*pi/10*sqrt(-1))*u)); pdeplot(p,e,t,'xydata',ur,'colorbar','off','mesh','off'); caxis([-maxu maxu]); axis tight, set(gca,'DataAspectRatio',[1 1 1]); axis off M(:,j) = getframe; end movie(hf,M,50);
pdedemo2 contains a full command-line implementation of the scattering problem.
This example shows how to solve a nonlinear problem for this equation:
$$-\nabla \cdot \left(\frac{1}{\sqrt{1+{\left|\nabla u\right|}^{2}}}\nabla u\right)=0$$
where the coefficients c, a, and f do not depend only on x and y, but also on the solution u.
The problem geometry is a unit disk, specified as Ω = {(x, y) | x^{2} + y^{2} ≤ 1}, with u = x^{2} on ∂Ω.
This nonlinear and cannot be solved with the regular elliptic solver. Instead, the nonlinear solver pdenonlin is used.
This example show how to solve this minimal surface problem using both the PDE app and command-line functions.
Make sure that the application mode in the PDE app is set to Generic Scalar. The problem domain is simply a unit circle. Draw it and move to the boundary mode to define the boundary conditions. Use Select All from the Edit menu to select all boundaries. Then double-click a boundary to open the Boundary Condition dialog box. The Dirichlet condition u = x^{2} is entered by typing x.^2 into the r edit box. Next, open the PDE Specification dialog box to define the PDE. This is an elliptic equation with
$$c=\frac{1}{\sqrt{1+{\left|\nabla u\right|}^{2}}},\text{}a=0,\text{and}f=0.$$
The nonlinear c is entered into the c edit box as
1./sqrt(1+ux.^2+uy.^2)
Initialize a mesh and refine it once.
Before solving the PDE, select Parameters from the Solve menu and check the Use nonlinear solver option. Also, set the tolerance parameter to 0.001.
Click the = button to solve the PDE. Use the Plot Selection dialog box to plot the solution in 3-D (check u and continuous selections in the Height column) to visualize the saddle shape of the solution.
This example shows how to solve the minimal surface problem using command-line functions. The files circleg and circleb2 contain the geometry specification and boundary condition functions, respectively.
g = 'circleg'; b = 'circleb2'; c = '1./sqrt(1+ux.^2+uy.^2)'; rtol = 1e-3; [p,e,t] = initmesh(g); [p,e,t] = refinemesh(g,p,e,t); u = pdenonlin(b,p,e,t,c,0,0,'Tol',rtol); pdesurf(p,t,u)
You can also run this example by typing pdedemo3.
This example shows how to perform one-level domain decomposition for complicated geometries, where you can decompose this geometry into the union of more subdomains of simpler structure. Such structures are often introduced by the PDE app.
Assume now that Ω is the disjoint union of some subdomains Ω_{1}, Ω_{2}, . . . , Ωn. Then you could renumber the nodes of a mesh on Ω such that the indices of the nodes of each subdomain are grouped together, while all the indices of nodes common to two or more subdomains come last. Since K has nonzero entries only at the lines and columns that are indices of neighboring nodes, the stiffness matrix is partitioned as follows:
$$K=\left(\begin{array}{ccccc}{K}_{1}& 0& \cdots & 0& {B}_{1}^{T}\\ 0& {K}_{2}& \cdots & 0& {B}_{2}^{T}\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0& 0& \cdots & {K}_{n}& {B}_{n}^{T}\\ {B}_{1}& {B}_{2}& \cdots & {B}_{n}& C\end{array}\right)$$
while the right side is
$$F=\left(\begin{array}{c}{f}_{1}\\ {f}_{2}\\ \vdots \\ {f}_{n}\\ {f}_{c}\end{array}\right)$$
The Partial Differential Equation Toolbox function assempde can assemble the matrices K_{j}, B_{j}, f_{j}, and C separately. You have full control over the storage and further processing of these matrices.
Furthermore, the structure of the linear system
Ku = F
is simplified by decomposing K into the partitioned matrix.
Now consider the geometry of the L-shaped membrane. You can plot the geometry of the membrane by typing
pdegplot('lshapeg')
Notice the borders between the subdomains. There are three subdomains. Thus the matrix formulas with n = 3 can be used. Now generate a mesh for the geometry:
[p,e,t] = initmesh('lshapeg'); [p,e,t] = refinemesh('lshapeg',p,e,t); [p,e,t] = refinemesh('lshapeg',p,e,t);
So for this case, with n = 3, you have
$$\left(\begin{array}{cccc}{K}_{1}& 0& 0& {B}_{1}^{T}\\ 0& {K}_{2}& 0& {B}_{2}^{T}\\ 0& 0& {K}_{3}& {B}_{3}^{T}\\ {B}_{1}& {B}_{2}& {B}_{3}& C\end{array}\right)\left(\begin{array}{c}{u}_{1}\\ {u}_{2}\\ {u}_{3}\\ {u}_{c}\end{array}\right)=\left(\begin{array}{c}{f}_{1}\\ {f}_{2}\\ {f}_{3}\\ {f}_{c}\end{array}\right),$$
and the solution is given by block elimination:
$$\begin{array}{l}(C-{B}_{1}{K}_{1}^{-1}{B}_{1}^{T}-{B}_{2}{K}_{2}^{-1}{B}_{2}^{T}-{B}_{3}{K}_{3}^{-1}{B}_{3}^{T}){u}_{c}={f}_{c}-{B}_{1}{K}_{1}^{-1}{f}_{1}-{B}_{2}{K}_{2}^{-1}{f}_{2}-{B}_{3}{K}_{3}^{-1}{f}_{3}\\ {u}_{1}={K}_{1}^{-1}({f}_{1}-{B}_{1}^{T}{u}_{c})\\ \cdots \end{array}$$
In the following MATLAB solution, a more efficient algorithm using Cholesky factorization is used:
time = []; np = size(p,2); % Find common points c = pdesdp(p,e,t); nc = length(c); C = zeros(nc,nc); FC = zeros(nc,1); [i1,c1] = pdesdp(p,e,t,1);ic1 = pdesubix(c,c1); [K,F] = assempde('lshapeb',p,e,t,1,0,1,time,1); K1 = K(i1,i1);d = symamd(K1);i1 = i1(d); K1 = chol(K1(d,d));B1 = K(c1,i1);a1 = B1/K1; C(ic1,ic1) = C(ic1,ic1)+K(c1,c1)-a1*a1'; f1 = F(i1);e1 = K1'\f1;FC(ic1) = FC(ic1)+F(c1)-a1*e1; [i2,c2] = pdesdp(p,e,t,2);ic2 = pdesubix(c,c2); [K,F] = assempde('lshapeb',p,e,t,1,0,1,time,2); K2 = K(i2,i2);d = symamd(K2);i2 = i2(d); K2 = chol(K2(d,d));B2 = K(c2,i2);a2 = B2/K2; C(ic2,ic2) = C(ic2,ic2)+K(c2,c2)-a2*a2'; f2 = F(i2);e2 = K2'\f2;FC(ic2) = FC(ic2)+F(c2)-a2*e2; [i3,c3] = pdesdp(p,e,t,3);ic3 = pdesubix(c,c3); [K,F] = assempde('lshapeb',p,e,t,1,0,1,time,3); K3 = K(i3,i3);d = symamd(K3);i3 = i3(d); K3 = chol(K3(d,d));B3 = K(c3,i3);a3 = B3/K3; C(ic3,ic3) = C(ic3,ic3)+K(c3,c3)-a3*a3'; f3 = F(i3);e3 = K3'\f3;FC(ic3) = FC(ic3)+F(c3)-a3*e3; % Solve u = zeros(np,1); u(c) = C\ FC; u(i1) = K1\(e1-a1'*u(c1)); u(i2) = K2\(e2-a2'*u(c2)); u(i3) = K3\(e3-a3'*u(c3));
The problem can also be solved by typing
% Compare with solution not using subdomains [K,F] = assempde('lshapeb',p,e,t,1,0,1);u1 = K\F; norm(u-u1,'inf') pdesurf(p,t,u)
You can run this entire example by typing pdedemo4.