Note: This page has been translated by MathWorks. Please click here

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

**THIS PAGE DESCRIBES
THE LEGACY WORKFLOW.** New features might not be
compatible with the legacy workflow. For the corresponding step in
the recommended workflow, see the recommended examples on the PDE Coefficients page.

This example shows how to write PDE coefficients in character form and in functional form for 2-D 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

The PDE is parabolic,

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:

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

The boundary conditions on the outer boundary (segments 1 through 4) are Dirichlet, with the value , where is time. Suppose the circular boundary (segments 5 through 8) has a generalized Neumann condition, with and .

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 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.

The initial condition is at $ t = 0$.

u0 = 0;

Create the mesh.

generateMesh(model,'GeometricOrder','linear');

Set the time steps for the parabolic solver to 50 steps from time 0 to time 1.

tlist = linspace(0,1,50);

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);

126 successful steps 9 failed attempts 272 function evaluations 1 partial derivatives 35 LU decompositions 271 solutions of linear systems

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

Equivalently, you can write a function for the coefficient `f`

in the syntax described in 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);

126 successful steps 9 failed attempts 272 function evaluations 1 partial derivatives 35 LU decompositions 271 solutions of linear systems

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);

126 successful steps 9 failed attempts 272 function evaluations 1 partial derivatives 35 LU decompositions 271 solutions of linear systems

Was this topic helpful?