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

    %% Solve PDE with Coefficients in Functional Form
%%
% This example shows how to write PDE coefficients in character form and in
% functional form for 2-D geometry.
%% Geometry
% The geometry is a rectangle with a circular hole. Create a PDE model
% container, and incorporate the geometry into the container.
model = createpde(1);

% Rectangle is code 3, 4 sides,
% followed by x-coordinates and then y-coordinates
R1 = [3,4,-1,1,1,-1,-.4,-.4,.4,.4]';
% Circle is code 1, center (.5,0), radius .2
C1 = [1,.5,0,.2]';
% Pad C1 with zeros to enable concatenation with R1
C1 = [C1;zeros(length(R1)-length(C1),1)];
geom = [R1,C1];
% Names for the two geometric objects
ns = (char('R1','C1'))';
% Set formula
sf = 'R1 - C1';
% Create geometry
gd = decsg(geom,sf,ns);

% Include the geometry in the model
geometryFromEdges(model,gd);
% View geometry
pdegplot(model,'EdgeLabels','on')
xlim([-1.1 1.1])
axis equal
%% PDE Coefficients
% The PDE is parabolic,
%
% $$d \frac{\partial u}{\partial t} - \nabla \cdot \left( c \nabla u \right) + au= f$$
%
% with the following coefficients:
%
% * d = 5
% * a = 0
% * f is a linear ramp up to 10, holds at 10, then ramps back down to 0:
%
% $$
% u = 10 \cdot
% \left\{
%  \begin{array}{ll}
% 10t & 0 \leq t \leq 0.1 \\ 1 & 0.1 \leq t \leq 0.9 \\ 10-10t & 0.9 \leq t \leq 1
% \end{array}  \right.
% $$
%
% * $c = 1 + x^2 + y^2$
%
% Write a function for the |f| coefficient.
%  function f = framp(t)
%  
%  if t <= 0.1
%      f = 10*t;
%  elseif t <= 0.9
%      f = 1;
%  else
%      f = 10-10*t;
%  end
%  f = 10*f;
%  end
%% Boundary Conditions
% The boundary conditions on the outer boundary (segments 1 through 4) are
% Dirichlet, with the value $u \left( x,y \right) = t \left( x - y
% \right)$, where $t$ is time. Suppose the circular boundary (segments 5
% through 8) has a generalized Neumann condition, with $q = 1$ and $g = x^2
% + y^2$.
myufun = @(region,state)state.time*(region.x - region.y);
mygfun = @(region,state)(region.x.^2 + region.y.^2);
applyBoundaryCondition(model,'edge',1:4,'u',myufun,'Vectorized','on');
applyBoundaryCondition(model,'edge',5:8,'q',1,'g',mygfun,'Vectorized','on');
%%
% The boundary conditions are the same as in
% <http://www.mathworks.com/help/pde/ug/boundary-conditions-by-writing-functions.html?searchHighlight=Boundary%20Conditions%20for%20Scalar%20PDE#btj13_7
% Boundary Conditions for Scalar PDE>. That description uses the older
% function form for specifying boundary conditions, which is no longer
% recommended. This description uses the recommended object form.
%% Initial Conditions
% The initial condition is $u \left( x,y \right) = 0$ at $ t = 0$.
u0 = 0;
%% Mesh
% Create the mesh.
generateMesh(model);
%% Parabolic Solution Times
% Set the time steps for the parabolic solver to 50 steps from time 0 to
% time 1.
tlist = linspace(0,1,50);
%% Solution
% Solve the parabolic PDE.
d = 5;
a = 0;
f = 'framp(t)';
c = '1 + x.^2 + y.^2';
u = parabolic(u0,tlist,model,c,a,f,d);
%%
% View an animation of the solution.
for tt = 1:size(u,2) % number of steps
    pdeplot(model,'XYData',u(:,tt),'ZData',u(:,tt),'ColorMap','jet')
    axis([-1 1 -1/2 1/2 -1.5 1.5 -1.5 1.5]) % use fixed axis
    title(['Step ' num2str(tt)])
    view(-45,22)
    drawnow
    pause(.1)
end
%% Alternative Coefficient Syntax
% Equivalently, you can write a function for the coefficient |f| in the
% syntax described
% in <http://www.mathworks.com/help/pde/ug/scalar-coefficients-in-function-form.html?searchHighlight=Specify%202-D%20Scalar%20Coefficients%20in%20Function%20Form#responsive_offcanvas
% Specify 2-D Scalar Coefficients in Function Form>.
%
%  function f = framp2(p,t,u,time)
%  
%  if time <= 0.1
%      f = 10*time;
%  elseif time <= 0.9
%      f = 1;
%  else
%      f = 10 - 10*time;
%  end
%  f = 10*f;
%  end
%%
% Call this function by setting
f = @framp2;
u = parabolic(u0,tlist,model,c,a,f,d);
%%
% You can also write a function for the coefficient |c|, though it is more
% complicated than the character formulation.
%  function c = cfunc(p,t,u,time)
%  
%  % Triangle point indices
%  it1 = t(1,:);
%  it2 = t(2,:);
%  it3 = t(3,:);
%  
%  % Find centroids of triangles
%  xpts = (p(1,it1) + p(1,it2) + p(1,it3))/3;
%  ypts = (p(2,it1) + p(2,it2) + p(2,it3))/3;
%  
%  c = 1 + xpts.^2 + ypts.^2;
%  end
%%
% Call this function by setting
c = @cfunc;
u = parabolic(u0,tlist,model,c,a,f,d);