Initial conditions for PDEPE

Question

Matterazzi 2023-3-31

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1938429-initial-conditions-for-pdepe

评论： Matterazzi 2023-4-4

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-4-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1938429-initial-conditions-for-pdepe#answer_1206379

编辑：Torsten 2023-4-1

在 MATLAB Online 中打开

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;

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Bill Greene 2023-3-31

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1938429-initial-conditions-for-pdepe#answer_1205679

编辑：Bill Greene 2023-3-31

在 MATLAB Online 中打开

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 个评论
显示 14更早的评论隐藏 14更早的评论

Matterazzi 2023-4-4

在 MATLAB Online 中打开

sr_data.txt

@Torsten, I am not sure if I should formulate this as a different question but I will ask here since it is a continuation of the problem above. Why is k equal to initial guess k0? I have tried several optimization approaches but nothing seems to solve the problem.

close all

clear

clc

% This problem does not solve because the function outpu does not change at all

% when your optimization parameters change a tiny amount.

% Experimental data

Y=load('sr_data.txt'); % t C16O18O H218O CH3C16O18OH C18O2

tdata = Y(1:end,1);

c9 = Y(1:end,2);

c8 = Y(1:end,3);

c5 = Y(1:end,4);

c7 = Y(1:end,5);

c = [c5,c7,c8,c9]; % 36x4

k0 = [.02; .03; .01; .06; .03; .002; .4; .5; .01; .05; .1; .1];

% Three options to minimize

% lsqcurvefit returns error: Initial point is a local minimum

% lsqnonlin returns error: 'Exiting due to infeasibility: 11 lower bounds exceed the corresponding upper bounds.'

% patternsearch: from the web problem???

options = optimoptions('lsqcurvefit','FiniteDifferenceStepSize',5e-3, 'Algorithm','trust-region-reflective');

lb = zeros(12,1);

%lb = [];

ub = [];

[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,tdata,c,lb,ub,options)

Initial point is a local minimum. Optimization completed because the size of the gradient at the initial point is less than the value of the optimality tolerance.

k = 12×1

0.0200 0.0300 0.0100 0.0600 0.0300 0.0020 0.4000 0.5000 0.0100 0.0500

Rsdnrm = 9.5548e+04

Rsd = 36×4

-0.0024 -0.0030 -0.0046 0 -0.0058 0.0433 3.7644 3.7552 -0.0061 0.2084 7.9363 7.9050 -0.0051 0.4801 11.6359 11.5685 -0.0052 0.8454 14.9631 14.8392 -0.0179 1.2656 17.8793 17.7765 -0.0336 1.7379 20.5297 20.4439 -0.0593 2.2617 22.9251 22.7920 -0.1002 2.8260 25.0457 24.9382 -0.1277 3.4086 26.9769 26.8719

ExFlg = 1

OptmInfo = struct with fields:

firstorderopt: 0 iterations: 0 funcCount: 13 cgiterations: 0 algorithm: 'trust-region-reflective' stepsize: 1 message: 'Initial point is a local minimum.↵↵Optimization completed because the size of the gradient at the initial point ↵is less than the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The final point is the initial point.↵The first-order optimality measure, 0.000000e+00, is less than↵options.OptimalityTolerance = 1.000000e-06.' bestfeasible: [] constrviolation: []

Lmda = struct with fields:

lower: [12×1 double] upper: [12×1 double] eqlin: [] ineqlin: [] eqnonlin: [] ineqnonlin: []

Jmat =

All zero sparse: 144×12

Cfit = kinetics(k, tdata);

% Plotting

plot(tdata, c, 'x')

hold on

plot(tdata, Cfit, '-',"LineWidth",1.2)

function C = kinetics(k,t)

global A B D Ncomp

% Constants

tau = 11.8;

u = .9;

L = .3;

D_ax = 1.5e-3;

Bo = 8.9;%u*L/D_ax;

rho_B = 2000;

eps_B = 0.45;

Ncomp = 9;

% solution mesh

x = linspace(0,1,20); % space

t = linspace(0,100,36); % time

% Other parameters;

A = -1/tau;

B = 1/tau*Bo;

D = rho_B/eps_B;

% Solve equation

m=0;

options=odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7);

sol = pdepe(m, @pdefun_sr, @icfun_sr, @bcfun_sr, x, t,options);

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)*1E-2;

u10 = 6;

f = [B*1e-7*dudx(1); % diffusion term

B*1e-7*dudx(2);

B*1e-7*dudx(3);

B*1e-7*dudx(4);

B*dudx(5);

B*1e-7*dudx(6);

B*dudx(7);

B*dudx(8);

B*dudx(9)

];

s = [k(7).*u(3) - k(1).*u(1); % source term

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

% Initial conditions

function u0 = icfun_sr(x)

u0 = zeros(Ncomp,1);

end

% Boundary conditions

function [pl,ql,pr,qr] = bcfun_sr(xl,ul,xr,ur,t,A,B,D)

pl = ul;

ql = zeros(Ncomp,1);

pr = zeros(Ncomp,1);

qr = ones(Ncomp,1); % or ones(dudx)*Ncomp times

end

% Extract modeled data used for fitting experimental data (at reactor exit)

u5 = sol(:,20,5);

u7 = sol(:,20,7);

u8 = sol(:,20,8);

u9 = sol(:,20,9);

C = [u5 u7 u8 u9]; % very good

end

Torsten 2023-4-4

在 MATLAB Online 中打开

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;

Matterazzi 2023-4-4

I am slapping myself now for being so blind. Thank you very much @Torsten.

请先登录，再进行评论。

Initial conditions for PDEPE

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

16 个评论
显示 14更早的评论隐藏 14更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Initial conditions for PDEPE

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

16 个评论 显示 14更早的评论隐藏 14更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

16 个评论
显示 14更早的评论隐藏 14更早的评论