This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

Note: This page has been translated by MathWorks. Click here to see
To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

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(x,y)=52+20x.

You can implement this simple interpolation in an anonymous function.

myufunction = @(location,state)52 + 20*location.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 u1(x)=52+20x+10sin(πx3).

The second component satisfies u2(x)=52-20x-10sin(πx3).

Write a function file myufun.m that incorporates these equations in the syntax described in Nonconstant Boundary Conditions.

function bcMatrix = myufun(location,state)
bcMatrix = [52 + 20*location.x + 10*sin(pi*(location.x.^3));
    52 - 20*location.x - 10*sin(pi*(location.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)