ODEs with Structs in codes

Hello all, I'm new to Matlab and trying to estimate two parameters p : p(1) and p(2) in ode45. I've followed the code by Star Strider's 'Igor_Moura'' function and created my own lsqcurvefit function; however, it got some issues in my case.

Is it a problem that i put the "struts" as "para" in 'odefun' function and in the 'XP' function as input?

Why not put the parameter 'p' also as an output in 'odefun' function?

%Data
p0=[0,0]
lb=[0,0]
ub=[10,200]

%Parameters Estimation 
[p]=lsqcurvefit(@XP,p0,xdata,Ydata,lb,ub)

Here's the 'heatpara' function. I'm not sure if the struts are ok to be put as an input in 'odefun' function.

function Yfit=XP(p,xlength)
para.n = 4;                                    
para.Y_0 = linspace(562,569,para.n);          %%%inlet temperature versus r_i radius
Y0 = zeros(input.n,1); %x0 is an matrix for 10*1
Y0(1:para.n) = para.Y_0;  
[x, Y] = ode45(@(x,Y) odefcn(x,Y,para), xlength, Y0, para);
    function ddx= odefun(x,T,para)
        
        ddx(input.n) = A.*(-(p(2)).*(Y(para.n)-para.Y2)-B.*para.p(1);
            
        ddx=ddx';
    end
Yfit=Y;
end

Any advice or help would be greatly appreciated!

Answer 1

[Torsten]

What is the dimension of Tdata ? And what are your Tdata ?

%%%input as an array of structs for fixed value variables
input.v_z = 0.93;                   % m/s
input.n = 4;                      % number of control volume
%input.T_0r = linspace(562,569,input.n);         % K
input.T_k = 593;                   % K  (T_wall)
input.d_t=0.014;                   % m    %%%reactor radius
input.rho = 0.9;                   % kg/m^3
input.c_P = 1300;                  % J/kg/K
input.rho_cat = 2200;              % kg/m^3
input.w_A_0r = 0.03;
input.k_01 = 5.05.*10.^2;          % mol/kg/s
input.E_A1 = 63.1.*1000;           % J/mol
input.n_1 = 0.85;                  % -
input.k_02 = 2.2.*10.^4;           % mol/kg/s
input.E_A2 = 85.*1000;             % J/mol
input.n_2 = 1.2;                   % -
input.R = 8.314;                   % J/mol/K
input.delta_r_h1 = 3125.*1000;     % J/mol
input.delta_r_h2 = 3410.*1000;     % J/mol
input.epsilon_cat = 0.443;         % -
z_length=linspace(0,0.84,13);       % like tspan

% Discretization of radius
input.r_Ri = linspace(0,input.d_t,input.n+1)';     % Position of Edge point (face)
input.delta_r = diff(input.r_Ri)';                 % Radius of Middle Points
input.delta_r = input.delta_r(1)';
input.r_i = input.delta_r./2+input.r_Ri(1:end-1);
%Data
%parameter to be estimated with the right value
% p(1) as lambda_r = 0.65;             % W/m/K
% p(2) as k_W = 160;                   % W/m^2/K
p0=[0,0];
lb=[0,0];
ub=[10,200];
zdata = z_length;
Tdata= [289.6  290.2 296.5	320
294.6	295.4	298.1	320
298.2   298.8	302.1   320
301.1	302.3	307.4 	320
302.0	303.2	308.3   320
304.9	304.6	308.6	320
306.9	306.7   309.6	320
307.7	308.4	311.0	320
308.8	309.2	314.5   320
310.0	311.7	315.1	320
312.2	310.6	315.6	320
312.2	311.2   315.8	320
308.8 	306.7	312.1 	320]+273.15;
input.T_0r = Tdata(1,:);
%Para Estimation
[p, Rsd]=lsqcurvefit(@(p,z)heatpara(p,z,input),p0,zdata,Tdata,lb,ub)

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

p = 1×2

1.4551 200.0000

Rsd = 1.8409e+03

function Tfit=heatpara(p,z_length,input)
T0 = input.T_0r;  %%%%%Inlet Temperature
[z, Tfit] = ode15s(@odefun,z_length, T0);
  function dTdz= odefun(z,T)
  C = p(1)./(input.rho.*input.c_P.*input.v_z);
  D = 1./(input.rho.*input.c_P.*input.v_z);
  E = (1-input.epsilon_cat).*input.rho_cat;
  dTdz(1) = C.*((input.r_i(1)+input.delta_r./2).*(T(2)-T(1))./(input.r_i(1).*input.delta_r.^2))...
    +D.*((-input.delta_r_h1).*E.*input.k_01.*exp(-input.E_A1./(input.R.*T(1))).*(input.w_A_0r).^input.n_1...
    +(-input.delta_r_h2).*E.*input.k_02.*exp(-input.E_A2./(input.R.* T(1))).*(input.w_A_0r).^input.n_2);
  dTdz(2:input.n-1) = C.*((input.r_i(2:input.n-1)+input.delta_r./2).*(T(3:input.n)-T(2:input.n-1))./(input.r_i(2:input.n-1).*input.delta_r.^2)...
    -(input.r_i(2:input.n-1)-input.delta_r./2).*(T(2:input.n-1)-T(1:input.n-2))./(input.r_i(2:input.n-1).*input.delta_r.^2))...
    +D.*((-input.delta_r_h1).*E.*input.k_01.*exp(-input.E_A1./(input.R.*T(2:input.n-1))).*((input.w_A_0r)).^input.n_1...
    +(-input.delta_r_h2).*E.*input.k_02.*exp(-input.E_A2./(input.R.*T(2:input.n-1))).*((input.w_A_0r)).^input.n_2);
  dTdz(input.n) = D.*(-(p(2)).*(T(input.n)-input.T_k).*(input.r_i(input.n)+input.delta_r./2)./(input.r_i(input.n).*input.delta_r)...
    -p(1).*(input.r_i(input.n)-input.delta_r./2).*(T(input.n)-T(input.n-1))./(input.r_i(input.n).*input.delta_r.^2))...
    +D.*((-input.delta_r_h1).*E.*input.k_01.*exp(-input.E_A1./(input.R.*T(input.n)))*(input.w_A_0r).^input.n_1...
    +(-input.delta_r_h2).*E.*input.k_02.*exp(-input.E_A2./(input.R.*T(input.n))).*(input.w_A_0r).^input.n_2);
  dTdz=dTdz.';
  end
end

[Torsten]

Code adapted.

And don't use ode45 for your system ; you need a stiff integrator (ode15s).

[uki71319]

Dear @Torsten, thank you very very much, the code finally works!!

But i still have one questions is: Why is it better to use ode15s intead of ode45 for this equation?

My original PDE equation is dTdz=C*(1/r)(d(r*dTdr)dr)+D*source term

Thanks in advance!!!

[Torsten]

Q2: Why is it better to use ode15s intead of ode45 for this equation?

My original PDE equation is dTdz=C*(1/r)(d(r*dTdr)dr)+D*source term

Thought the original PDE solutuon for simulating the temperature profile is using ode45, so i continued using ode45 solver for estimating parameter before.

I already gave the answer. Discretized reaction-diffusion equations like yours are highly stiff in general. That's why you need ODE15S as a stiff solver. It should solve your ODEs much faster than ODE45 does.