# Matlab code for system of differential equations (chemical kinetics) fitting to data

35 views (last 30 days)
Brandon Breitling on 29 May 2013
Commented: Renee on 11 Nov 2014
Chemical kinetics fitting as I understand the mathematical process is to: 1)Setup the system of differential equations for the connected processes 2)Use variation of parameters to solve for the coefficients and the differential equations so that everything is in terms of rate constants. 3)Use Runga-Kutta or some other similar method to fit these equations to experimental data to get the rate constants.
Can anyone layout these steps in Matlab code-wise? I think I can even get everything solved for in terms of rate constants. I really just need to know how do I fit a function with parameters to data with code.
If there is an easy way to do this give just data and the system of differentiation equations(unsolved) that would be even better.

Matt Kindig on 29 May 2013
Edited: Matt Kindig on 29 May 2013
To be clear, you want to fit experimental data to a system of differential equations with some number of unknown parameters, correct? If so, your workflow will basically be performing a numerical optimization, with the ODE solution nested inside the optimization, with the following general step. Note that you don't need to have explicitly solved the ODEs to run this approach (and in fact, it is not necessary to even have an analytical solution to them).
1. Set up the system of differential equations in the dx = f(t,x) format required by Matlab's ODE solvers. (I assume since you mention Runge-kutta you have ODEs). See the Documentation for 'ode45' (for example) to see how to do this:
doc ode45
Let's call that function 'myodefcn' and have it take three inputs: your unknown parameters p, the time vector t and state vector x. This function could have the form:
function dx= myodefcn(p, t, x)
%your equations are here, in the form dx= f(t,x)
dx = p(1)*x-2*p(2); %for example
2. Create an objective function that solves the ODE equations and outputs a state variable that can be compared to your data. Calculate the difference between the model data and your experimental data, and sum the squared error (to perform a least-squares fit). Suppose you are modeling a 2-parameter first-order system, and you have experimental data time td and state xd. Such a function might look like:
function SSE = myobjective(td,xd, p, x0)
f = @(t,x) myodefcn(p,t,x); %function to call ode45
[ts,xs] = ode45(f, td, x0); %xs is the state variable predicted by the model
err = xd - xs;
SSE = sum(err.^2); %sum squared-error.
with inputs td (experimental time), xd (experimental data), p (model parameters), and x0 (initial conditions).
3. Use one of Matlab's optimization functions to solve for the parameters that minimize the SSE. If you have the Optimization Toolbox, you have lots of flexibility here, to implement constraints, bounds, etc. on the parameters. If you have basic Matlab, you can still use the fminsearch function:
psolve = fminsearch(@(p) myobjective(td,xd,p,x0), p0);
where p0 is your initial guesses for the unknown parameters p.
This is the general approach. Unfortunately, without more details, I cannot give anything more specific. Let me know if you have any questions.
Renee on 11 Nov 2014
Hi, I dont know if this is ok (to ask a question in the answers section), but I have a very similar problem to the original poster and I am trying to implement your answer and I have a problem. I am confused about passing tdata to my ode solving function. This is the code I have:
c0=[.005;100;.005;100;.1;2.6;.1;2.6;.5;1.16;.2;2.2e10;.05;.01;];
x0 =[20000;75;0;0;0;0;0;0;0;0;.05;0;.05;];
psolve = fminsearch(@(c) myoptim(tdata,xdata,c,x0), c0);
sse function:
function SSE = myoptim(tdata,xdata,c,x0)
f=@(t,x)fitodes_newsplit(c,x);
[ts,xs]=ode45(f,tdata,x0);
err=xd-xs;
SSE=sum(err.^2);
my ode function:
function y = fitodes_newsplit(c,x)
% Define ODE system with parameter c
f = @(t,x,c) [-c(1).*x(1).*x(3)+c(2).*x(4)-c(3).*x(1).*x(5)+c(4).*x(6);
-c(5).*x(2).*x(3)+c(6).*x(5)-c(7).*x(2).*x(4)+c(8).*x(6);
x(12)-c(1).*x(1).*x(3)+c(2).*x(4)-c(5).*x(2).*x(3)+c(6).*x(5)+c(9).*x(9).*x(8);
c(1).*x(1).*x(3)-c(2).*x(4)-c(7).*x(2).*x(4)+c(8).*x(6);
c(5).*x(2).*x(3)-c(6).*x(5)-c(3).*x(1).*x(5)+c(4).*x(6);
c(7).*x(2).*x(4)-c(8).*x(6)+c(3).*x(1).*x(5)-c(4).*x(6)-c(10).*x(6);
c(9).*x(9).*x(8);
c(10).*x(6).*(1-c(11))-c(9).*x(9).*x(8);
c(10).*x(6).*c(11)-c(9).*x(9).*x(8);
c(10).*x(6).*(1-c(11)).*c(12)-x(10);
-x(13).*x(11).*c(13)+x(12).*c(14);
x(13).*x(11).*c(13)-x(12).*c(14);
-x(13).*x(11).*c(13)+x(12).*c(14);];
%init conditions
init =[20000;75;0;0;0;0;0;0;0;0;.05;0;.05;];
%solve
[~,z] = ode23s(@(t,y) f(t,y,c),x,init);
% Extract desired solution component
y = z(:,10);
Here's my errors:
Error using fitodes_newsplit Too many input arguments.
Error in myoptim>@(t,x)fitodes_newsplit(c,t,x) (line 7) f=@(t,x)fitodes_newsplit(c,t,x);
Error in odearguments (line 87) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 113) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in myoptim (line 8) [ts,xs]=ode45(f,tdata,x0);
Error in @(c)myoptim(tdata,xdata,c,x0)
Error in fminsearch (line 190) fv(:,1) = funfcn(x,varargin{:});
Error in optimtest1 (line 5) psolve = fminsearch(@(c) myoptim(tdata,xdata,c,x0), c0);