# Initial conditions for PDEPE

3 views (last 30 days)
Matterazzi on 31 Mar 2023
Commented: Matterazzi on 4 Apr 2023
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

Torsten on 1 Apr 2023
Edited: Torsten on 1 Apr 2023
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

Bill Greene on 31 Mar 2023
Edited: Bill Greene on 31 Mar 2023
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
.
Matterazzi on 4 Apr 2023
I am slapping myself now for being so blind. Thank you very much @Torsten.

### Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

R2020a

### Community Treasure Hunt

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

Start Hunting!