IF condition in the function containing the PDE - pdepe solver

[Patrice Perreault]

Hi,

I'm using pdepe to find optimal parameters for kinetic constants while performing curve fitting with lsqcurvefit. Everything works perfectly (exact same [c,f,s] as below but without if conditions) as long as I don't use "if" conditions in the function containing the pde ([c,f,s]).

The problem is the following :
I need to solve the mole balance (convection + diffusion terms), and following the height, some reactive rate terms (empty reactor tube, followed by a small reaction zone, then some more empty reactor tube).

I though of using a
if z >= under the reactive zone
if z <= reactive zone
c = [1;1;1;1];
f = [D_z*dCdz(1); D_z*dCdz(2); 0; 0];
s = [-(C(1)*dUdz + U_z_z*dCdz(1)) - k1*C(1)*Theta(C(3),C(4))*(1-eps)/eps;
-(C(2)*dUdz + U_z_z*dCdz(2)) + k2*C(3)*C(4)^2*(1-eps)/eps; ...
(k1*C(1)*Theta(C(3),C(4)) - k2*C(3)*C(4)^2)*R*T/P; ...
-k2*C(3)*C(4)^2*R*T/P];
else % Non reactive zone
c = [1;1;1;1];
f = [D_z*dCdz(1); D_z*dCdz(2); 0; 0];
s = [-(C(1)*dUdz + U_z_z*dCdz(1));
-(C(2)*dUdz + U_z_z*dCdz(2)); ...
0; ...
0];
end
c = [1;1;1;1];
f = [D_z*dCdz(1); D_z*dCdz(2); 0; 0];
s = [-(C(1)*dUdz + U_z_z*dCdz(1));
-(C(2)*dUdz + U_z_z*dCdz(2)); ...
0; ...
];
end

The only other z-dependant term is the velocity U_z_z = f(z).

Any suggestions?

Complete code (the one that works) :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [c,f,s] = pdefun_lsq(z,t,C,dCdz,k)
% c(x,t,u,du/dx)*du/dt = x^(-m)*d/dx[x^m*f(x,t,u,du/dx)] + s(x,t,u,du/dx)
% C(1) = y_CO
% C(2) = y_CO2
% C(3) = Theta_CO*
% C(4) = Theta_O*



L_reac = 0.30; % m
V_gaz = 40; % Nml/min
h_lit = 0.01; % m
D_reac = 7E-3; % m
P = 101325; % Pa
T = 813; % K
R = 8.314; % Pa*m^3/mol/K
eps = 0.45; %

D_z = 8.7E-5; % m^2/s

U_z_0 = V_gaz*1E-6/60*4/(pi*D_reac^2);


U_z_z = U_z_0*(-2E-05*(z-.28).^4 + 0.0014*(z-.28).^3 - 0.0402*(z-.28).^2 ...
+ 0.5635*(z-.28) + 1.3);


dUdz = U_z_0*(4*-2E-05*(z-.28).^3 + 3*0.0014*(z-.28).^2 ...
- 2*0.0402*(z-.28) + 0.5635);

]

k1 = k(1);
k2 = k(2);

c = [1; 1; 1; 1];
f = [D_z*dCdz(1); D_z*dCdz(2); 0; 0];
s = [-(C(1)*dUdz + U_z_z*dCdz(1)) - k1*C(1)*Theta(C(3),C(4))*(1-eps)/eps; ...
-(C(2)*dUdz + U_z_z*dCdz(2)) + k2*C(3)*C(4)^2*(1-eps)/eps; ...
(k1*C(1)*Theta(C(3),C(4)) - k2*C(3)*C(4)^2)*R*T/P; ...
-k2*C(3)*C(4)^2*R*T/P];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Error message :
??? Error using ==> feval
Output argument "c" (and maybe others) not assigned during call to
"C:\Users\Patrice\Documents\MATLAB\Modélisation_CO\PDE - ajustement de
courbe\Ajustement-CO-CO2-archive\3
paramètres\pdefun_lsq_ifz.m>pdefun_lsq_ifz".

Error in ==> pdepe at 259
c = feval(pde,xi(1),t(1),U,Ux,varargin{:});

Error in ==> fun_pdepe_lsq at 18
sol = pdepe(m,@pdefun_lsq_ifz,@icfun_lsq,@bcfun_lsq,z,t,[],k);

Error in ==> lsqncommon at 152
fuser = funfcn{3}(xargin{:},varargin{:});

Error in ==> lsqcurvefit at 186
[x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...

Error in ==> PDE_solve_lsq at 55
[K,resnorm,residual,exitflag] = lsqcurvefit(@fun_pdepe_lsq,k,t,ydata,
...

Caused by:
Failure in initial user-supplied objective function evaluation.
LSQCURVEFIT cannot continue.

Thank you
Patrice

[Torsten]

I suspect that your if-statement does not cover the complete z-range
for which pdepe requires the values for c, f and s.

Best wishes
Torsten.

[Bruno Luong]

"Patrice Perreault" <patrice....@polymtl.ca> wrote in message <jl33h6$lhf$1...@newscl01ah.mathworks.com>...

>
> I though of using a
> if z >= under the reactive zone

[ snip ...]
> end
>

What happens if your function is called with z < under the reactive zone ?

Your code is also (Matlab) syntactically incorrect, which make hard to follow.

Bruno