Initial conditions for PDEPE
조회 수: 3 (최근 30일)
이전 댓글 표시
Hi community,
I am trying to solve 9 PDEs using pdepe routine for which the initial condition for all these equations is zero, Ci(t=0) = 0. My solution mesh is as:
% Solution mesh
x = linspace(0,1,50); % space
t = linspace(0,100,100); % time
I attempted to code the initial conditions as:
% icfun_sr() defines the initial conditions
function u0 = icfun_sr(x)
u0 = zeros(9,1);
end
which I can then pass to solve the equations as:
% Solve equations
m=0;
sol = pdepe(m, @pdefun_sr, @icfun_sr, @bcfun_sr, x, t, [], A, B, D);
The error thrown as
Error using sr_code>icfun_sr
Too many input arguments.
Error in pdepe (line 229)
temp = feval(ic,xmesh(1),varargin{:});
Error in sr_code (line 32)
sol = pdepe(m, @pdefun_sr, @icfun_sr, @bcfun_sr, x, t, [], A, B, D);
Full code below:
close all
clear all
clc
global A B D Ncomp
% Constants
Rtime = 2;
u = 0.7;
L = 1;
D_ax = 1.5e-3;
Bo = u*L/D_ax;
density = 2000;
poros = 0.45;
Ncomp = 9;
% mesh
x = linspace(0,1,50); % space
t = linspace(0,100,100); % time
% Other parameters;
A = -1/Rtime;
B = 1/Rtime*Bo;
D = density/poros;
% Solve equation
m=0;
pdef = @(x, t, u, dudx) pdefun_sr(x, t, u, dudx);
sol = pdepe(m, pdef, @icfun_sr, @bcfun_sr, x, t);
% Extract solutions
for i = 1:Ncomp
u(i) = sol(:,:,i);
end
waterfall(x,t,u(1)), map=[0 0 0]; colormap(map), view(15,60)
function [c,f,s] = pdefun_sr(x, t, u, dudx)
global A B D Ncomp
c = ones(Ncomp,1); % diagonal terms
k = ones(12,1);
u(10) = 6; % or u(10) = 6 - exp(-10*t);
f = [0; % diffusion term
0;
0;
0;
A*u(5) + B*dudx(5);
0;
A*u(7) + B*dudx(7);
A*u(8) + B*dudx(8);
A*u(9) + B*dudx(9)
];
% source term
s = [k(7).*u(3) - k(1).*u(1);
k(9).*u(6) - k(2).*u(2);
k(12).*u(10) - (k(3) + k(5) + k(8) + k(7)).*u(3); % fails at u(10)
k(10).*u(6) + k(11).*u(5) - k(6).*u(4) - k(4).*u(4);
D*( k(6).*u(4) - k(11).*u(5) );
k(8).*u(3) - (k(9) + k(10)).*u(6);
D*( k(1).*u(1) );
D*( k(2).*u(2) + k(3).*u(3) );
D*( k(5).*u(3) + k(4).*u(4) )
];
end
% icfun_sr() defines the initial conditions
function u0 = icfun_sr(x)
global Ncomp
u0 = zeros(Ncomp,1);
end
function [pl,ql,pr,qr] = bcfun_sr(xl,ul,xr,ur,t,A,B,D)
% bcfun_sr() defines the initial conditions
global Ncomp
pl = zeros(Ncomp,1).*ul;
ql = zeros(Ncomp,1);
pr = zeros(Ncomp,1);
qr = ones(Ncomp,1); % or ones(dudx)*Ncomp times
end
댓글 수: 0
채택된 답변
Torsten
2023년 4월 1일
편집: Torsten
2023년 4월 1일
Maybe you could add a little diffusivity to the components with spurious oscillations (those with f=0). If this is not possible: I gave you two alternatives.
global A B D Ncomp
% Constants
Rtime = 2;
u = 0.7;
L = 1;
D_ax = 1.5e-3;
Bo = u*L/D_ax;
density = 2000;
poros = 0.45;
Ncomp = 9;
% mesh
x = linspace(0,1,500); % space
t = linspace(0,100,100); % time
% Other parameters;
A = -1/Rtime;
B = 1/Rtime*Bo;
D = density/poros;
% Solve equation
m=0;
sol = pdepe(m, @pdefun_sr, @icfun_sr, @bcfun_sr, x, t);
u = sol(40,:,6);
plot(x,u)
function [c,f,s] = pdefun_sr(x, t, u, dudx)
global A B D Ncomp
c = ones(Ncomp,1); % diagonal terms
k = ones(12,1);
u10 = 6; % or u10 = 6 - exp(-10*t);
f = [0; % diffusion term
0;
0;
0;
B*dudx(5);
0;
B*dudx(7);
B*dudx(8);
B*dudx(9)
];
% source term
s = [k(7).*u(3) - k(1).*u(1);
k(9).*u(6) - k(2).*u(2);
k(12).*u10 - (k(3) + k(5) + k(8) + k(7)).*u(3); % fails at u(10)
k(10).*u(6) + k(11).*u(5) - k(6).*u(4) - k(4).*u(4);
A*dudx(5) + D*( k(6).*u(4) - k(11).*u(5) );
k(8).*u(3) - (k(9) + k(10)).*u(6);
A*dudx(7) + D*( k(1).*u(1) );
A*dudx(8) + D*( k(2).*u(2) + k(3).*u(3) );
A*dudx(9) + D*( k(5).*u(3) + k(4).*u(4) )
];
end
% icfun_sr() defines the initial conditions
function u0 = icfun_sr(x)
global Ncomp
u0 = zeros(Ncomp,1);
end
function [pl,ql,pr,qr] = bcfun_sr(xl,ul,xr,ur,t,A,B,D)
% bcfun_sr() defines the initial conditions
global Ncomp
pl = ul;
ql = zeros(Ncomp,1);
pr = zeros(Ncomp,1);
qr = ones(Ncomp,1); % or ones(dudx)*Ncomp times
end
댓글 수: 0
추가 답변 (1개)
Bill Greene
2023년 3월 31일
편집: Bill Greene
2023년 3월 31일
The problem is that you are using an old, undocumented form of the of the pdepe function
to pass additional parameters to the functions called by pdepe. If you use this form, you must
include the extra parameters in ALL of the called functions, e.g.
function u0 = icfun_sr(x,A,B,C)
The reason that this old form of pdepe is undocumented is that anonymous functions provide a
superior way of passing extra parameters to called functions. In your example, you could define
the anonymous function
pdef = @(x, t, u, dudx) pdefun_sr(x, t, u, dudx, A, B, D);
and call pdepe
sol = pdepe(m, pdef, @icfun_sr, @bcfun_sr, x, t);
But even after changing your code to either of the two solutions, you will still get errors due
to the variable k being undefined in function pdefun_sr
.
댓글 수: 16
Torsten
2023년 4월 4일
Why is k equal to initial guess k0?
Because you overwrite the parameters supplied by lsqcurvefit in this line:
k = ones(12,1)*1E-2;
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!