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

Learn moreOpportunities for recent engineering grads.

Apply Today
Asked by TP on 15 Jan 2013

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

- dS/dt = kt * (A - S) - ka * S * (Rmax - R)+ kd * R;
- 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, [], []);

*No products are associated with this question.*

Answer by Shashank 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:

Show 1 older comment

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 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

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

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

## 0 Comments