Time Dependent Boundary Conditions in PDE Toolbox

I'm struggling to get my time dependent boundary condition function working. I went back to the beginning and copied the example straight out of the documentation for "Boundary Conditions for Scalar PDE".

There appears to be two issues with this example. First, the critical step of generating the mesh is noticeably absent. The boundary condition function, as shown, only accounts for the eight edges of the basic geometry. However, once the initmesh function is performed, the edge matrix becomes much larger (all the edges of the individual triangular elements) and the edge boundary labels in the boundary condition function (1-8) are no longer valid.

Second, and probably most important, is that the "time" variable is an empty matrix causing errors in the calculations.

Has anyone had any success with this method?

Thanks - Michael

Here is the script I'm using:

% Time Variant Boundary Conditions
clear
clc

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

% View geometry
figure(1)
pdegplot(gd,'edgeLabels','on')
xlim([-1.1 1.1])
axis equal

[p, e, t] = initmesh(gd,'Hmax',Inf);
% [p, e, t] = refinemesh(gd,p,e,t);
% p = jigglemesh(p,e,t);

figure(2)
pdemesh(p,e,t)
xlim([-1.1 1.1])
axis equal

b = @BC_fun;
u0 = 0;
c = 1;
a = 0;
f = 10;
d = 1;
tlist = 0:.01:1;

u = parabolic(u0,tlist,b,p,e,t,c,a,f,d);

And here is the function file:

function [ qmatrix, gmatrix, hmatrix, rmatrix ] = BC_fun( p, e, u, time )
% PDE Boundary Condition Function

ne = size(e,2); % number of edges
qmatrix = zeros(1,ne);
gmatrix = qmatrix;
hmatrix = zeros(1,2*ne);
rmatrix = hmatrix;

for k = 1:ne
    x1 = p(1,e(1,k)); % x at first point in segment
    x2 = p(1,e(2,k)); % x at second point in segment
    xm = (x1 + x2)/2; % x at segment midpoint
    y1 = p(2,e(1,k)); % y at first point in segment
    y2 = p(2,e(2,k)); % y at second point in segment
    ym = (y1 + y2)/2; % y at segment midpoint
    switch e(5,k)
        case {1,2,3,4} % rectangle boundaries
            hmatrix(k) = 1;
            hmatrix(k+ne) = 1;
            rmatrix(k) = time*(x1 - y1);
            rmatrix(k+ne) = time*(x2 - y2);
        otherwise % same as case {5,6,7,8}, circle boundaries
            qmatrix(k) = 1;
            gmatrix(k) = xm^2 + ym^2;
    end

end

Here is the error returned:

In an assignment  A(I) = B, the number of elements in B and I must be the same.
Error in BC_fun (line 21)
            rmatrix(k) = time*(x1 - y1);
Error in pdeefxpd (line 10)
  [q,g,h,r]=feval(bl,p,e,u,time);
Error in pdeexpd (line 40)
[q,g,h,r]=pdeefxpd(p,e,u,time,bl);
Error in pdeODEInfo/setN (line 178)
        [q,g,h,r]=pdeexpd(self.p,self.e(:,1),self.b);
Error in pdeODEInfo/checkFuncDepen (line 56)
      self=self.setN(u0);
Error in pdeParabolicInfo (line 14)
      obj=obj.checkFuncDepen(u0, tlist);
Error in parabolic (line 41)
  pdeprb=pdeParabolicInfo(u0,tlist,b,p,e,t,c,a,f,d);
Error in TimeVaryBC (line 46)
u = parabolic(u0,tlist,b,p,e,t,c,a,f,d);

[Michael]

I've continued working this problem, and I think I am starting to understand the real issue that I am having. The function PARABOLIC in the Partial Differential Equation Toolbox can't handle discontinuous functions with respect to time. Example:

function [ f ] = f_coeff( p, t, u, time )

f = zeros(1,size(t,2));

for j = 1:size(t,2)
    if t(4,j) == 2
        f(j) = 50*(time >= 5*pi);
    end
end

end

Back to the drawing board!

Answer 1

[Bill Greene]

Michael,

Before commenting on your f-function, let me make sure I understand the situation. This issue is independent of your original post, right? I assume you have resolved the original problem?

As far as your f_coeff function, I don't know of any reason why PDE Toolbox would necessarily have problems with coefficients that are discontinuous with respect to time. What problem are you seeing? One problem I see immediately is that your original geometry (your script above) has only one subdomain. Your f_coeff function is setting f(j) to be non-zero only for a non-existent subdomain 2.

Another issue concerns how PDE Toolbox detects whether a coefficient is time dependent or not. As a test, the code passes time=NaN and expects the function to return a result that contains at least one NaN. This is documented here: http://www.mathworks.com/help/pde/ug/pde-coefficients.html but is easy to overlook.

I suggest adding code like this to the top of your f_coeff function:

if(isnan(time))
f = NaN*ones(1,size(t,2));
end

Bill

[Michael]

I added a subdomain 2 by changing the formula to:

   sf = 'R1+C1';

Now that I've solved the first issue, I'll look at your latest suggestion.

[Michael]

What I found worked was the following:

     if(isnan(time))
          f = NaN*ones(1,size(t,2));
     else
          % my coeffcient conditions
     end

I feel like I'm taking leaps and bounds toward my end goal now. Thanks Bill for all your help.

Michael

[SUBHAM KASHYAP]

Hey,

I want to include seismic motions in the pde calculations. This essentially requires defining time dependent acceleration or displacement boundary conditions at the fixed support. How can I possibly do it using MATLAB pde toolbox.

Answer 2

[Bill Greene]

Hi Michael,

It looks like you may have encountered a bug in the code. A work-around is to replace the line in your script that sets the initial conditions:

u0 = 0;

with these lines:

np = size(p,2);
u0 = zeros(np,1);

that explicitly set the initial conditions at all points in the mesh to zero.

Regards,

Bill

[Michael]

Thanks Bill this did allow the code to run. I wasn't making the connection to the initial condition becasue it would run fine with a scalar initial condition (u0 = 0;) until I included the time variable in the boundary condition function.