Numerical integration with nonlinear least squares curve fitting?

I have a system which can be modelled by coupled differential equations:

1. dS/dt = kt * (A - S) - ka * S * (Rmax - R)+ kd * R;
2. dR/dt = ka * S * (Rmax - R) - kd * R;

Therefore, variables are: S, R and unknown parameters are: kt, ka, kd and Rmax.

I need to find the unknown parameters and simulate a curve that will fit my data. I have measured R as a function of time experimentally, which is the only data I have. I know the value of A. How do I carry out numerical integration and use nonlinear least squares curve fitting on my data?

Here is something I tried, but the calculation goes on for hours until I have to abort it manually.

1st m-file:

function S = NumInt(U, t)

%       dS/dt = kt * (A - S)- ka * S * (Rmax - R)+ kd * R;
%       dR/dt = (ka * S * (Rmax - R) - kd * R;
%       Variables:  y(1) = S,  y(2) = R;
%       Parameters: kt = U(1),  ka = U(2),  kd = U(3),  Rmax = U(4)

y0 = rand(2,1);

A= 50;

   [T,Y] = ode45(@Int, t, y0);
   function dy = Int(t, y)
   dy = zeros(2,1);
   dy(1) = U(1) * (A - y(1)) - U(2) * y(1) * (U(4)- y(2) + U(3) * y(2));
   dy(2) = U(2) * y(1) * (U(4) - y(2)) - U(3) * y(2);

   end

end

2nd m-file:

U0 = rand(4,1) * 10;
[U] = lsqcurvefit(@NumInt, U0, t, y, [], []);

Answer 1

[Shashank Prasanna]

Since you are trying to fit a curve to an ODE, this page in the documentation has good pointers on the best solver to choose and how to go about it:

http://www.mathworks.com/help/optim/ug/optimizing-a-simulation-or-ordinary-differential-equation.html

[Shashank Prasanna]

Since the optimization just minimizes the objective function, make sure that you write your objective function such that it finds the sum of squared errors of the difference between the data and the fitted model for each iteration

norm((ydata-y)^2)

[TP]

Thanks for directing me to the ode page. I am now using ode15s which seems to have solved the calculation time issue.

However, now I am getting the following error:

Error using lsqcurvefit (line 247)
Function value and YDATA sizes are incommensurate.

Error in testcall (line 4)
[U] = lsqcurvefit('NumInt',U0,t,y,[],[]);

My y data is 283 x 1, t is 283 x 1.

I don't know how to solve this. Can you please help me?

[Shashank Prasanna]

I haven't seen you entire code, but this link should help you:

http://www.mathworks.com/support/solutions/en/data/1-19B8E/index.html

[TP]

Thank you!

Answer 2

[Alan Weiss]

Your function NumInt returns a solution structure S. Your lsqcurvefit call takes (xdata,ydata) arguments as t and y. Did you do some internal conversions to ensure that the correct data is passed between your functions?

Also, I hope you realize that the previous comment about using norm((ydata-y)^2) is incorrect. The note in the lsqcurvefit function reference page explains this very important point.

Alan Weiss

MATLAB mathematical toolbox documentation

[TP]

Thanks Alan!

(a) There is indeed an issue with the correct data being passed between the functions. Could you please redirect me/show me how to do it correctly?

(b) How can I get the output [T,Y] as variables?

Here is an example dataset and the entire code:

Data

t = [0:0.5:58]'
y =

    2.7000
    0.7000
         0
    0.1000
    0.4000
    1.0000
    1.4000
   15.8000
   45.0000
   52.4000
   56.1000
   59.0000
   61.1000
   63.0000
   64.3000
   65.6000
   67.0000
   68.4000
   69.8000
   70.8000
   71.9000
   72.9000
   73.7000
   74.8000
   75.5000
   76.2000
   77.0000
   77.9000
   78.5000
   79.4000
   79.9000
   80.6000
   81.0000
   81.7000
   82.1000
   82.6000
   83.0000
   83.6000
   84.1000
   83.9000
   84.5000
   85.1000
   85.6000
   86.0000
   86.1000
   86.5000
   86.9000
   87.0000
   87.6000
   87.9000
   88.1000
   88.4000
   88.4000
   88.7000
   88.9000
   89.3000
   89.3000
   89.5000
   89.5000
   89.8000
   89.9000
   89.9000
   90.4000
   90.2000
   90.2000
   90.4000
   90.6000
   90.4000
   90.6000
   90.7000
   90.7000
   90.9000
   90.9000
   90.9000
   91.1000
   91.0000
   90.9000
   90.9000
   90.5000
   90.7000
   90.7000
   90.9000
   90.8000
   91.1000
   90.9000
   91.0000
   90.8000
   91.0000
   91.0000
   90.9000
   90.8000
   91.0000
   90.8000
   90.8000
   91.0000
   90.9000
   91.3000
   91.2000
   91.0000
   90.9000
   91.0000
   91.0000
   91.1000
   90.7000
   90.6000
   90.6000
   90.7000
   90.9000
   90.8000
   90.5000
   90.9000
   90.7000
   90.6000
   90.9000
   90.6000
   90.8000
   90.6000

(1) NumInt.m

function S = NumInt(U,t);
% To solve the differential equations:
%
%       dS/dt = kt * (A - S)- ka * S * (Rmax - R)+ kd * R;
%       dR/dt = (ka * S * (Rmax - R) - kd * R;
%
%       Variables:  dydt(1) = S,  dydt(2) = R;
%       Parameters: kt = U(1),  ka = U(2),  kd = U(3),  Rmax = U(4);%
%       Known parameter: A = injected concentration%
%       Collected data: change is response R (y) vs time t (t)

 A = 50;
 xvalues = [1:1:117];

    [T,Y] = ode15s(@Int,xvalues,[0 0]);

   figure(1); plot(T,Y(:,1),'r'); %this is S
   figure(2); plot(T,Y(:,2),'b'); %this is R

   function dydt = Int(~,y);
      dydt=zeros(2,1);
      dydt(1) = U(1) * (A-y(1)) - U(2) * y(1) * (U(4)-y(2)+U(3) * y(2));
      dydt(2) = U(2) * y(1) * (U(4)-y(2))-U(3) * y(2);
   end

 end

(2) myfit.m

U0 = [0.1 1000 0.3 200];
options = optimset('MaxFunEvals',1e5,'MaxIter',1e5);

    [U,Resnorm,Rsd] = lsqcurvefit(@NumInt,U0,t,y,0,[],options);

I am quite new to matlab and any help will be much appreciated.

[Alan Weiss]

The function reference page gives the syntax you want:

 [T,Y] = ode15s(@Int,xvalues,[0 0]);

Instead of having NumInt return S, it should return [T,Y]:

 function [T,Y] = NumInt(U,t);

This outputs the times T and values Y of the calculated solution points of the ODE.

Use these values, which I suppose become the xdata and ydata for lsqcurvefit, and you should be able to continue.

You might also want to set DiffMinChange to a higher value than default, perhaps 1e-4, as suggested here.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation