Accelerating the pace of engineering and science

# Documentation

## Scalar PDE Functional Form and Calling Syntax

This example shows how to write PDE coefficients in string form and in functional form.

Geometry

The geometry is a rectangle with a circular hole.

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:

$f=10*\left\{\begin{array}{ll}10t\hfill & 0\le t\le 0.1\hfill \\ 1\hfill & 0.1\le t\le 0.9\hfill \\ 10-10t\hfill & 0.9\le t\le 1\hfill \end{array}$

• c = 1 +.x2 + y2

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

Boundary conditions

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.

Suppose the boundary conditions on the outer boundary (segments 1 through 4) are Dirichlet, with the value u(x,y) = t(x – y), where t is time. Suppose the circular boundary (segments 5 through 8) has a generalized Neumann condition, with q = 1 and g = x2 + y2.

```myufun = @(problem,region,state)state.time*(region.x - region.y);
mygfun = @(problem,region,state)(region.x.^2 + region.y.^2);
pg = pdeGeometryFromEdges(gd); % create geometry container
bc1 = pdeBoundaryConditions(pg.Edges(1:4),'u',myufun,'Vectorized','on');
bc2 = pdeBoundaryConditions(pg.Edges(5:8),'q',1,'g',mygfun,'Vectorized','on');
problem = pde(); % scalar problem
problem.BoundaryConditions = [bc1,bc2]; % all boundary conditions```

Initial Conditions

The initial condition is u(x,y) = 0 at t = 0.

`u0 = 0;`

Mesh

Create the mesh, refine it twice, and jiggle it once.

```[p,e,t] = initmesh(gd);
[p,e,t] = refinemesh(gd,p,e,t);
[p,e,t] = refinemesh(gd,p,e,t);
p = jigglemesh(p,e,t);```

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,problem,p,e,t,c,a,f,d);```

View an animation of the solution.

```for tt = 1:size(u,2) % number of steps
pdeplot(p,e,t,'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 Scalar PDE 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;```

Call this function by setting

```f = @framp2;
u = parabolic(u0,tlist,problem,p,e,t,c,a,f,d);```

You can also write a function for the coefficient c, though it is more complicated than the string 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;```

Call this function by setting

```c = @cfunc;
u = parabolic(u0,tlist,problem,p,e,t,c,a,f,d);```