## Numerical integration with nonlinear least squares curve fitting?

### TP (view profile)

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

## Products

No products are associated with this question.

### Shashank Prasanna (view profile)

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

TP

### TP (view profile)

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.

Shashank Prasanna

### Shashank Prasanna (view profile)

on 16 Jan 2013

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

TP

on 17 Jan 2013

Thank you!

### Alan Weiss (view profile)

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

TP

### TP (view profile)

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

### Alan Weiss (view profile)

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

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