Time Dependent Boundary Conditions in PDE Toolbox

9 views (last 30 days)
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);
  1 Comment
Michael
Michael on 7 Aug 2013
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!

Sign in to comment.

Accepted Answer

Bill Greene
Bill Greene on 6 Aug 2013
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
  3 Comments
Michael
Michael on 7 Aug 2013
Edited: Michael on 7 Aug 2013
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
SUBHAM KASHYAP on 9 Jul 2023
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.

Sign in to comment.

More Answers (1)

Bill Greene
Bill Greene on 3 Aug 2013
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
  1 Comment
Michael
Michael on 7 Aug 2013
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.

Sign in to comment.

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!