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

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.

Answer 1

[Matt Kindig]

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.

[Brandon Breitling]

So here would be my functions:

function dydt = Bigode(t,Y,k) dydt =zeros(4,1); dydt(3) = k(1)*Y(1) + k(2)*Y(2)-k(3)*Y(3); dydt(1) = -k(1)*Y(1); dydt(2) = -k(2)*Y(2);

function SSE = myobjective(td,Yd, k, Y0) f = @(t,Y) Bigode(k,t,Y); %function to call ode45 [ts,Ys] = ode45(f, td, Y0); %xs is the state variable predicted by the model err = Yd - Ys; SSE = sum(err.^2); %sum squared-error.

And my time t and concentration Y variables with initial guesses for the parameters k0

t=[0 15 30 45 60 75 90 105 120 ]

Y=[100.00 54.46 50.66 39.86 33.33 29.22 25.40 24.45 22.08 ]

k0=[0.05 0.05 0.05 ]

I'm getting a bunch of errors with this, can you help me? I'm pretty new to matlab.

[Matt Kindig]

Can you edit your comment above to format your code properly first? Highlight the code portions, and click the {} Code button above to do the formatting. Also, include the code that actually calls your optimization (e.g., your fminsearch, fmincon, or other call).

Second, can you post your errors?

Third, have you tried running each function separately? For example, first call Bigode() with some sample values of t, Y, and k-- does this throw errors? Second, call myobjective() with your data set--does it return a number? This kind of debugging is part of the fun of Matlab :)

[Brandon Breitling]

So here are the two .m file functions I made:

function dydt = Bigode(t,Y,k) dydt =zeros(4,1); dydt(3) = k(1)*Y(1) + k(2)*Y(2)-k(3)*Y(3); dydt(1) = -k(1)*Y(1); dydt(2) = -k(2)*Y(2);

function SSE = myobjective(td,Yd, k, Y0) f = @(t,Y) Bigode(k,t,Y); %function to call ode45 [ts,Ys] = ode45(f, td, Y0); %xs is the state variable predicted by the model err = Yd - Ys; SSE = sum(err.^2); %sum squared-error.

And my time t and concentration Y variables with initial guesses for the parameters k0

t=[0 15 30 45 60 75 90 105 120 ]

Y=[63.5271   36.4729         0
 54.9296    0.0000   17.4348
 47.4956    0.0000   17.2703
 41.0677    0.0000    6.6331
 35.5098    0.0000    3.0565
 30.7040    0.0000    1.4077
 26.5487    0.0000    2.3618
 22.9557    0.0000    0.8042
 19.8489         0    2.4638 ]

k0=[0.05 0.05 0.05 ]

And the code I ran to solve the system of differential equations:

psolve = fminsearch(@(p) myobjective(td,Yd,k,Y0), k0);

Here are my errors:

Undefined function or variable 'td'.

Error in @(k)myobjective(td,Yd,k,Y0)

Error in fminsearch (line 191) fv(:,1) = funfcn(x,varargin{:});

Running Bigode with the variables t Y and k0 gives:

ans =

    -3.1764
   -2.7465
    3.5481
         0

Running myobjective with

myobjective(t,Y, k0, Y)

gives:

Attempted to access Y(2); index out of bounds because numel(Y)=1.

Error in Bigode (line 3) dydt(3) = k(1)*Y(1) + k(2)*Y(2)-k(3)*Y(3);

Error in myobjective>@(t,Y)Bigode(k,t,Y) (line 2) f = @(t,Y) Bigode(k,t,Y); %function to call ode45

Error in odearguments (line 88) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 114) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in myobjective (line 3) [ts,Ys] = ode45(f, td, Y0); %xs is the state variable predicted by the model

[Matt Kindig]

First, I would rename your experimental variables from 't' and 'Y' to something else, as you are already using those variables in your ODE solver. Your call for myobjective() calls them 'td' and 'Yd'; use those names instead.

Also, where is Y0 defined?

P.S.: Your code still isn't formatted the way that's preferred on this forum. Please insert line breaks between each statement, such that the code is formatted the same way as in Matlab itself. This is to make it easier for us to copy-paste your code into Matlab, as well as identify syntax errors, comments, etc.

[Matt Kindig]

Another issue is that your anonymous function f in myobjective.m has the arguments switched from what Bigode.m expects. This is the correct line:

f = @(t,Y) Bigode(t,Y,k);  %note that this is the same order as Bigode.m

[Renee]

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:

xdata = xlsread('50nm split.xlsx','B2:B1201');
tdata = xlsread('50nm split.xlsx','A2:A1201');
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);