Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

To resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

Numerical integration with nonlinear least squares curve fitting?

Asked by TP on 15 Jan 2013

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, [], []);

0 Comments

TP

Products

No products are associated with this question.

2 Answers

Answer by Shashank Prasanna on 15 Jan 2013

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

4 Comments

TP on 16 Jan 2013

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 on 16 Jan 2013

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 on 17 Jan 2013

Thank you!

Shashank Prasanna
Answer by Alan Weiss on 17 Jan 2013

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

2 Comments

TP on 17 Jan 2013

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 on 21 Jan 2013

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

Alan Weiss

Contact us