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

    %% Solve PDEs with Nonconstant Boundary Conditions
%%
% This example shows how to write functions for a nonconstant boundary
% condition specification.
%% Geometry
% All the specifications use the same geometry, which is a rectangle with a
% circular hole.

% 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
g = decsg(geom,sf,ns);

% Create geometry model
model = createpde;

% Include the geometry in the model and view the geometry
geometryFromEdges(model,g);
pdegplot(model,'EdgeLabels','on')
xlim([-1.1 1.1])
axis equal
%% Scalar Problem
%
% * Edge 3 has Dirichlet conditions with value 32.
% * Edge 1 has Dirichlet conditions with value 72.
% * Edges 2 and 4 have Dirichlet conditions that linearly interpolate between edges 1 and 3.
% * The circular edges (5 through 8) have Neumann conditions with |q = 0|,
% |g = -1|.
applyBoundaryCondition(model,'dirichlet','edge',3,'u',32);
applyBoundaryCondition(model,'dirichlet','edge',1,'u',72);
applyBoundaryCondition(model,'neumann','edge',5:8,'g',-1); % q = 0 by default
%%
% Edges 2 and 4 need functions that perform the linear interpolation. Each
% edge can use the same function that returns the value $u \left( x,y
% \right) = 52 + 20x$.
%
% You can implement this simple interpolation in an anonymous function.
myufunction = @(region,state)52 + 20*region.x;
%%
% Include the function for edges 2 and 4. To help speed the solver, allow a
% vectorized evaluation.
applyBoundaryCondition(model,'dirichlet','edge',[2,4],'u',myufunction,'Vectorized','on');
%%
% Solve an elliptic PDE with these boundary conditions, using the
% parameters |c = 1|, |a = 0|, and | f = 10|. Because the shorter
% rectangular side has length 0.8, to ensure that the mesh is not too
% coarse choose a maximum mesh size |Hmax = 0.1|.
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',10);
generateMesh(model,'Hmax',0.1);
results = solvepde(model);
u = results.NodalSolution;
pdeplot(model,'XYData',u)
%% System of PDEs
% Suppose that the system has |N = 2|.
%
% * Edge 3 has Dirichlet conditions with values |[32,72]|.
% * Edge 1 has Dirichlet conditions with values |[72,32]|.
% * Edges 2 and 4 have Dirichlet conditions that interpolate between the
% conditions on edges 1 and 3, and include a sinusoidal variation.
% * Circular edges (edges 5 through 8) have |q = 0| and |g = -10|.
model = createpde(2);
geometryFromEdges(model,g);

applyBoundaryCondition(model,'dirichlet','edge',3,'u',[32,72]);
applyBoundaryCondition(model,'dirichlet','edge',1,'u',[72,32]);
applyBoundaryCondition(model,'neumann','edge',5:8,'g',[-10,-10]);
%%
% The first component of edges 2 and 4 satisfies the equation $u_1 \left( x
% \right) = 52 + 20x + 10 \sin \left( \pi x^3 \right)$.
%
% The second component satisfies $u_2 \left( x
% \right) = 52 - 20x - 10 \sin \left( \pi x^3 \right)$.
%
% Write a function file |myufun.m| that incorporates these equations in the
% syntax from "Nonconstant Boundary Conditions" in
% <http://www.mathworks.com/help/pde/ug/steps-to-specify-a-boundary-conditions-object.html
% Specify Boundary Conditions>.
%
%  function bcMatrix = myufun(region,state)
%  bcMatrix = [52 + 20*region.x + 10*sin(pi*(region.x.^3));
%      52 - 20*region.x - 10*sin(pi*(region.x.^3))]; % OK to vectorize
%  end
%%
% Include this function in the edge 2 and edge 4 boundary condition.
applyBoundaryCondition(model,'dirichlet','edge',[2,4],'u',@myufun,'Vectorized','on');
%%
% Solve an elliptic PDE with these boundary conditions, with the parameters
% |c = 1|, |a = 0|, and |f = (10,-10)|. Because the shorter rectangular
% side has length 0.8, to ensure that the mesh is not too coarse choose a
% maximum mesh size |Hmax = 0.1|.
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',[10;-10]);
generateMesh(model,'Hmax',0.1);
results = solvepde(model);
u = results.NodalSolution;

subplot(1,2,1)
pdeplot(model,'XYData',u(:,1),'ZData',u(:,1),'ColorBar','off')
view(-9,24)
subplot(1,2,2)
pdeplot(model,'XYData',u(:,2),'ZData',u(:,2),'ColorBar','off')
view(-9,24)