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!
采纳的回答
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)
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
3 个评论
Torsten
2022-10-24
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.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 General PDEs 的更多信息
Translated by 
