Fitting an ODE set with time dependent parameters to data

I have this set of ODE that I would like to fit to measured data of two time-dependent parameters Pa and Pe:

Both equations decrease exponentially as a function of time. To solve these equations I understood this needs to be programmed differently by integrating the exponential terms on forehand, while to fitting ODE to data is pretty straight forward.

Is it possible to do both at the same time? I have been breaking my head over this for a few days and I hope anyone here could help me out.

Answer 1

[Star Strider]

You could certainly estimate it by integrating it with ode45, however it has an analytic solution (kindly provided by the Symbolic Math Toolbox):

syms Pe(t) Pa(t) Pa0 Pe0 
a = sym('a', [1 6]);

DEq1 = diff(Pe) == a(1)*exp(a(2)*t)-a(3) * (Pe-Pa+1)
DEq2 = diff(Pa) == a(4)*exp(a(5)*t)-a(6) * (Pe-Pa+1)

Eqs = dsolve(DEq1,DEq2, Pa(0)==Pa0, Pe(0)==Pe0)
Pa = simplify(Eqs.Pa, 'Steps',500)
Pe = simplify(Eqs.Pe, 'Steps',500)
Eqs = matlabFunction([Pa;Pe], 'Vars',{[a Pa0, Pe0],t})

that after a bit of editing resolves to:

Eqs = @(in1,t)[-exp(-t.*(in1(:,3)-in1(:,6))).*((in1(:,6).*exp(t.*(in1(:,3)-in1(:,6)))-(in1(:,1).*in1(:,6).*exp(t.*(in1(:,2)+in1(:,3)-in1(:,6))))./(in1(:,2)+in1(:,3)-in1(:,6))+(in1(:,4).*in1(:,6).*exp(t.*(in1(:,3)+in1(:,5)-in1(:,6))))./(in1(:,3)+in1(:,5)-in1(:,6)))./(in1(:,3)-in1(:,6))-(in1(:,6).*(-in1(:,1).*in1(:,3)+in1(:,2).*in1(:,3)-in1(:,1).*in1(:,5)+in1(:,2).*in1(:,4)+in1(:,1).*in1(:,6)+in1(:,2).*in1(:,5)+in1(:,3).*in1(:,4)-in1(:,2).*in1(:,6)+in1(:,3).*in1(:,5)-in1(:,3).*in1(:,6).*2.0-in1(:,4).*in1(:,6)-in1(:,5).*in1(:,6)-in1(:,7).*in1(:,3).^2-in1(:,7).*in1(:,6).^2+in1(:,8).*in1(:,3).^2+in1(:,8).*in1(:,6).^2+in1(:,3).^2+in1(:,6).^2-in1(:,7).*in1(:,2).*in1(:,3)-in1(:,7).*in1(:,2).*in1(:,5)+in1(:,7).*in1(:,2).*in1(:,6)-in1(:,7).*in1(:,3).*in1(:,5)+in1(:,7).*in1(:,3).*in1(:,6).*2.0+in1(:,7).*in1(:,5).*in1(:,6)+in1(:,8).*in1(:,2).*in1(:,3)+in1(:,8).*in1(:,2).*in1(:,5)-in1(:,8).*in1(:,2).*in1(:,6)+in1(:,8).*in1(:,3).*in1(:,5)-in1(:,8).*in1(:,3).*in1(:,6).*2.0-in1(:,8).*in1(:,5).*in1(:,6)))./((in1(:,3)-in1(:,6)).*(in1(:,2)+in1(:,3)-in1(:,6)).*(in1(:,3)+in1(:,5)-in1(:,6))))-(in1(:,1).*in1(:,5).*in1(:,6).*exp(in1(:,2).*t)-in1(:,2).*in1(:,3).*in1(:,4).*exp(in1(:,5).*t))./(in1(:,2).*in1(:,5).*(in1(:,3)-in1(:,6)))-(in1(:,2).*in1(:,3).*in1(:,4)-in1(:,1).*in1(:,5).*in1(:,6)-in1(:,7).*in1(:,2).*in1(:,3).*in1(:,5)+in1(:,8).*in1(:,2).*in1(:,5).*in1(:,6))./(in1(:,2).*in1(:,5).*(in1(:,3)-in1(:,6)));
               -(in1(:,1).*in1(:,5).*in1(:,6).*exp(in1(:,2).*t)-in1(:,2).*in1(:,3).*in1(:,4).*exp(in1(:,5).*t))./(in1(:,2).*in1(:,5).*(in1(:,3)-in1(:,6)))-(in1(:,3).*exp(-t.*(in1(:,3)-in1(:,6))).*((in1(:,6).*exp(t.*(in1(:,3)-in1(:,6)))-(in1(:,1).*in1(:,6).*exp(t.*(in1(:,2)+in1(:,3)-in1(:,6))))./(in1(:,2)+in1(:,3)-in1(:,6))+(in1(:,4).*in1(:,6).*exp(t.*(in1(:,3)+in1(:,5)-in1(:,6))))./(in1(:,3)+in1(:,5)-in1(:,6)))./(in1(:,3)-in1(:,6))-(in1(:,6).*(-in1(:,1).*in1(:,3)+in1(:,2).*in1(:,3)-in1(:,1).*in1(:,5)+in1(:,2).*in1(:,4)+in1(:,1).*in1(:,6)+in1(:,2).*in1(:,5)+in1(:,3).*in1(:,4)-in1(:,2).*in1(:,6)+in1(:,3).*in1(:,5)-in1(:,3).*in1(:,6).*2.0-in1(:,4).*in1(:,6)-in1(:,5).*in1(:,6)-in1(:,7).*in1(:,3).^2-in1(:,7).*in1(:,6).^2+in1(:,8).*in1(:,3).^2+in1(:,8).*in1(:,6).^2+in1(:,3).^2+in1(:,6).^2-in1(:,7).*in1(:,2).*in1(:,3)-in1(:,7).*in1(:,2).*in1(:,5)+in1(:,7).*in1(:,2).*in1(:,6)-in1(:,7).*in1(:,3).*in1(:,5)+in1(:,7).*in1(:,3).*in1(:,6).*2.0+in1(:,7).*in1(:,5).*in1(:,6)+in1(:,8).*in1(:,2).*in1(:,3)+in1(:,8).*in1(:,2).*in1(:,5)-in1(:,8).*in1(:,2).*in1(:,6)+in1(:,8).*in1(:,3).*in1(:,5)-in1(:,8).*in1(:,3).*in1(:,6).*2.0-in1(:,8).*in1(:,5).*in1(:,6)))./((in1(:,3)-in1(:,6)).*(in1(:,2)+in1(:,3)-in1(:,6)).*(in1(:,3)+in1(:,5)-in1(:,6)))))./in1(:,6)-(in1(:,2).*in1(:,3).*in1(:,4)-in1(:,1).*in1(:,5).*in1(:,6)-in1(:,7).*in1(:,2).*in1(:,3).*in1(:,5)+in1(:,8).*in1(:,2).*in1(:,5).*in1(:,6))./(in1(:,2).*in1(:,5).*(in1(:,3)-in1(:,6)))];

where ‘in1’ are the parameters ‘a(1)’ to ‘a(6)’ in order, and the intital conditions ‘Pa(0)’ and ‘Pe(0)’ (also in that order) as a row vector, so the initial parameter estimates must be a row vector or the function will throw an error. (This is a pecularity of the functions created by the matlabFunction function.) The data you are fitting must be in row vectors as well, since they will be fitting [Pa; Pe] as row vectors, with ‘Pa’ the first row and ‘Pe’ the second row. I would use lsqcurvefit with these.

To clarify, unless I am missing something, these are time dependent equations, and do not appear to have time-dependent parameters, since the parameters do not appear to vary with time.

[Jos Rozema]

Thank you very much! After some trying I managed to make it work. I should have thought of finding an analytical solution first...

I will be trying some other ODEs on the data as well. I don't suppose I should expect them all to be analytical as well. In your opinion, what would be the best approach to handle those cases?

[Star Strider]

As always, my pleasure!

I would then use the numeric ODE solvers, and one of the parameter estimation functions.

A recent example using the ga function to estimate the parameters is Parameter Estimation for a System of Differential Equations and a similar older one using lsqcurvefit is Parameter Estimation for a System of Differential Equations (same title, similar approach except for the optimisation function).

You will likely need to adapt these to your data and differential equations. I will provide what help I can.

[Jos Rozema]

Thank you for the suggestions. That will be a useful starting point. :)

[Star Strider]

As always, my pleasure!