lsqnonlin fitting parameters with multiple equations (ODE) and dependency between parameters
16 views (last 30 days)
Show older comments
For an assignment, I'm trying to fit a set of four ODEs with experimental data. I think I need to use lsqnonlin, but I'm completely lost on how to implement this properly with ODEs and on top of that, the parameters are dependent on each other. I've worked with lsqnonlin before, but this is outside my abilities. The goal is as follows:
I have 4 conditions:
- The concentration of O2- (red line) should be 0 (or close to 0 like 0.001) after 1 sec
- the ratio CO/C2H4 should be 0.62 after 1 sec
- the ratio H2/C2H4 should be 0.64 after 1 sec
- the conversion of CH4 (CH4 (at t=0) - CH4 (at t=1))/CH4 (at t=0)) should be 0.405
as you can see, these are all boundary conditions at 1 sec time. I have 5 parameters that should be fitted, which are constants (k1 to k4) and the initial condition for CH4. Here is the code:
% Parameters to fit
k1=1;
k2=1;
k3=1;
k4=1;
CH4IN=0.076;
% Time interval and initial conditions
tmin = 0;
tmax = 100;
t_interval = [tmin,tmax];
init_cond = [CH4IN,0.076,0,0]'; %initial amount of respectively CH4, O2-, C2H4, H2
% Solution
[t,y] = ode45(@(t,Y) SOEfcn(t,Y,k1,k2,k3,k4) , t_interval , init_cond);
% Calculate and plot fluxes
% Ethane, CO en H2O cant be calculated directly, thus this must be
% Calculated in the following manner:
rEthane=k2.*y(:,3).*y(:,4).^2.*y(:,2);
Ethane = cumtrapz(rEthane);
rCO = 2*k3.*y(:,3).*y(:,2).^4;
CO = cumtrapz(rCO);
rH2O = k1.*y(:,1).^2.*y(:,2)+k2.*y(:,3).*y(:,4).^2.*y(:,2)+2.*k3.*y(:,3).*y(:,2).^4+k4.*y(:,4).*y(:,2);
H2O = cumtrapz(rH2O);
plot(t,y(:,1),'b',t,y(:,2),'r',t,y(:,3),'g',t,y(:,4),'y',t,Ethane,'v',t,CO,'p',t,H2O,'x');
legend('CH4','O2-','C2H4','H2','C2H6','CO','H2O')
xlabel('time (s)')
ylabel('concentration')
with SOEfcn (in ode45):
function dYdt = SOEfcn(t,Y,k,k2,k3,k4)
%Y(1) is CH4, Y(2) is oxide (O2-), Y(3) is ethylene (C2H4), Y(4) is H2
dYdt = [-2*k(1)*Y(1)^2*Y(2); %methane rate
-k(1)*Y(1)^2*Y(2)-k2*Y(3)*Y(4)^2*Y(2)-4*k3*Y(3)*Y(2)^4-k4*Y(4)*Y(2); %oxide rate
k(1)*Y(1)^2*Y(2)-k2*Y(3)*Y(4)^2*Y(2)-k3*Y(3)*Y(2)^4; %ethylene rate
%k2*Y(3)*Y(4)^2*Y(2); %ethane rate
%2*k3*Y(3)*Y(2)^4; %CO rate
k(1)*Y(1)^2*Y(2)-2*k2*Y(3)*Y(4)^2*Y(2)-k4*Y(4)*Y(2);]; %H2 rate
%k1*Y(1)^2*Y(2)+k2*Y(3)*Y(4)^2*Y(2)+2*k3*Y(3)*Y(2)^4+k4*Y(4)*Y(2);];%H2O rate
end
Just to be clear, the ethane, CO and H2 rate couldn't be solved within the ode45, so I had to calculate it differently.
Does anyone know how to fit the parameters to the given boundary conditions (I think with lsqnonlin, but if you know another method that's also fine)?
4 Comments
Torsten
on 25 May 2021
Take a look at the examples provided under
de.mathworks.com/help/optim/ug/fit-differential-equation-ode.html?lang=en),
Accepted Answer
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!