www.gusucode.com > pde 案例源码 matlab代码程序 > pde/pdedemo1.m

    %% Poisson's Equation on a Unit Disk
% This example shows how to numerically solve a Poisson's equation using
% the |solvepde| function in Partial Differential Equation Toolbox(TM).
%
% The particular PDE is
%
% $$ -\Delta u = 1$$
%
% on the unit disk with zero-Dirichlet boundary conditions. The exact
% solution is 
%
% $$ u(x,y) = \frac{1-x^2-y^2}{4}. $$
%
% For most partial differential equations, the exact solution is not known.
% In this example, however, we can use the known, exact solution to show
% how the error decreases as the mesh is refined.


% Copyright 1994-2015 The MathWorks, Inc.

%% Problem Definition
% The following variables will define our problem:
%%
% * |g|: A specification function that is used by |geometryFromEdges|. For
% more information, please see the documentation pages for |circleg| and
% the documentation section
% <matlab:helpview(fullfile(docroot,'toolbox','pde','helptargets.map'),'pde_geometry_fcn');
% Create Geometry Using a Geometry Function>.
% * |c|, |a|, |f|: The coefficients and inhomogeneous term.
g = @circleg;
c = 1;
a = 0;
f = 1;

%% PDE Coefficients and Boundary Conditions
% Plot the geometry and display the edge labels for use in the boundary
% condition definition.
figure 
pdegplot(g,'EdgeLabels','on'); 
axis equal
% Create a PDE Model with a single dependent variable
numberOfPDE = 1;
pdem = createpde(numberOfPDE);
% Create a geometry entity
geometryFromEdges(pdem,g);
% Specify PDE Coefficients
specifyCoefficients(pdem,'m',0,'d',0,'c',c,'a',a,'f',f);
% The solution is zero at all four outer edges of the circle
applyBoundaryCondition(pdem,'dirichlet','Edge',(1:4),'u',0);

%% Generate Initial Mesh
% The function |generateMesh| takes a geometry specification function and
% returns a discretization of that domain. The |'Hmax'| option lets the
% user specify the maximum edge length. In this case, because the domain is
% a unit disk, a maximum edge length of one creates a coarse discretization.
hmax = 1;
generateMesh(pdem,'Hmax',hmax);
figure
pdemesh(pdem); 
axis equal

%% Refinement
% We repeatedly refine the mesh until the infinity-norm of the error vector
% is less than $10^{-3}$. 
% 
% For this domain, each refinement halves the |'Hmax'| option.
er = Inf;
while er > 0.001
    hmax =hmax/2;
    generateMesh(pdem,'Hmax',hmax);
    result = solvepde(pdem);
    msh = pdem.Mesh;
    u = result.NodalSolution;
    exact = (1-msh.Nodes(1,:).^2-msh.Nodes(2,:).^2)'/4;
    er = norm(u-exact,'inf');
    fprintf('Error: %e. Number of nodes: %d\n',er,size(msh.Nodes,2));
end


%% Plot Final Mesh
figure
pdemesh(pdem); 
axis equal

%% Plot FEM Solution
figure
pdeplot(pdem,'XYData',u,'ZData',u,'ColorBar','off')

displayEndOfDemoMessage(mfilename)