Documentation Center

  • Trial Software
  • Product Updates

Scalar PDE Functional Form and Calling Syntax

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

The geometry is a rectangle with a circular hole.

 Code for generating the figure

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:

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

The boundary conditions are the same as in Boundary Conditions for Scalar PDE.

 Boundary conditions

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

After running the code for creating the geometry, 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);

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.

b = @pdebound;
d = 5;
a = 0;
f = 'framp(t)';
c = '1+x.^2+y.^2';
u = parabolic(0,tlist,b,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

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(0,tlist,b,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(0,tlist,b,p,e,t,c,a,f,d);

Related Examples

More About

Was this topic helpful?