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)