curve fit toolbox help system ODE and excel?

Hi, I am trying to curve fit using the toolbox to a function which is a system of ODE. my data corresponds to one equation of the function. how can i specify this? Also, the tutorial seems to say you should put your data into the command line, but mine is in excel, is there some easier way? Thanks!

Answer 1

[Matt Tearle]

If I interpret your problem correctly, you have a system of ODEs with a parameter -- y' = f(t,y;c) -- and an Excel spreadsheet that has data for t and y3 (or whatever -- one of the components of y). You're trying to do a curve fit to determine c.

That's fun!

Here's an example that does that with the classic mass-spring system. I'm making some t and y data; you could read this from Excel with xlsread. You need to do this only once at the beginning, then you have the data to fit to.

% ---------------------------------------------------------------
% Make some data -- replace all this with a call to XLSREAD
m = 0.1*rand;
k = 4*rand;
w = 0.5*sqrt(4*k-m*m);
t = linspace(0,2*pi)';
y = exp(m*t/2).*cos(w*t) + 0.02*randn(size(t));
% ---------------------------------------------------------------
% Do curvefitting
cfit = lsqcurvefit(@myode,rand(2,1),t,y)
% See results
plot(t,y,'o',t,myode(cfit,t))

And here's the curve-fitting function:

function y = myode(c,x)
% Define ODE system with parameter c
f = @(t,y,c) [y(2);-c(1)*y(1)+c(2)*y(2)];
% Solve ODE system with given initial condition
[~,z] = ode45(@(t,y) f(t,y,c),x,[1;0]);
% Extract desired solution component
y = z(:,1);

[Renee]

Thanks! So does "c" work like "y", in that c(1) is the first parameter etc? if so, how can I specify specific bounds for c(1) vs c(2) etc?

[Matt Tearle]

Yes, the parameters are specified as elements of a vector. You can pass vectors of bounds to lsqcurvefit. For example, if you have 3 parameters and c1 has to be between 0 and 5, and c3 has to be negative,

lb = [0; -Inf; -Inf];
ub = [5; Inf; 0];
cfit = lsqcurvefit(@myfun,c0,t,y,lb,ub)

[Renee]

Thanks! That solved one of the errors I was getting. Now it says "Dimensions of matrices being concatenated are not consistent." Do you happen to know what its referring to? There's 9 ODEs, and I have one initial condition for each ODE. Beyond that I am lost as to what dimensions its talking about.

    function y = splitodes(c,x)
% Define ODE system with parameter c
f = @(t,y,c) [-c(1) .*y(1) .*y(2) + c(2) .* y(6); - c(1) .* y(1) .*y(2) + c(2) .* y(6);-c(3) .*y(6).*y(3) + c(4) .*y(7) - c(5).*y(3).*y(8) + c(6) .*y(9);-c(7) .*y(4).*y(6) + c(8) .*y(8) - c(9) .*y(4).*y(7) + c(10) .*y(9);c(11) *c(12)*c(13) .* y(9) - y(5);c(1).*y(1).*y(2) - c(2).*y(6) -c(7) .*y(6).*y(4)+c(8).*y(8)-c(3).*y(3).*y(6) + c(4).*y(8) + c(11).*y(9);c(3).*y(6).*y(3)- c(4).*y(7)-c(9).*y(4).*y(7) + c(10).*y(9);-c(5).*y(8).*y(3)+c(6).*y(9) + c(7).*y(4).*y(6)-c(8).*y(8);c(5).*y(8).*y(3)-c(6).*y(9) + c(9).*y(4).*y(7)-c(10).*y(9) -c(11).*y(9)];
%error
options = odeset('RelTol',1e-6);
%init conditions
init = [1;0;0;0;0;0;0;0;0;];
%solve
[~,z] = ode45(@(t,y) f(t,y,c),x,init,options);
% Extract desired solution component
y = z(:,5);

[Matt Tearle]

This is weird, but be careful with your whitespace. MATLAB is pretty forgiving with spaces, BUT...

y(9) - c(11) is the same as y(9)-c(11). However, y(9) -c(11) is taken as two elements: y(9) and -c(11).

So, see how you have c(10).*y(9) -c(11).*y(9) right at the end of your definition of f? MATLAB's interpreting that as two equations concatenated horizontally (which is not OK, because the whole thing is supposed to be one column vector). Similarly, you have c(2).*y(6) -c(7) .*y(6) in the middle.

Clean those up and you should be right.

[Renee]

Thanks! One more question lol, is there a way I can have an equation to constrain sets of parameters to?

[Matt Tearle]

Yes... but it gets hairier. lsqcurvefit accepts only bound constraints. If you want other kinds of constraints, you'll need to go to something like fmincon. In this case, you'll need to write an objective function that is the error in the fit ( lsqcurvefit takes care of that for you).

[Renee]

Well basically my data is a logarithmic shape, and the output of lsqcurvefit (despite varying my initial param values) is consistently linear, and usually intersects the data at the origin and some other place. I dont think its taking nearly long enough to evaluate a system with 15 param and 11 ode, because it usually gives output in a few seconds. I added these options: help?

options = optimoptions(@lsqcurvefit,'TolFun',1.0000e-12,'TolX',1.0000e-12,'ScaleProblem','Jacobian','MaxIter',1500)

[Matt Tearle]

In general, nonlinear fitting problems with a lot of parameters can be very touchy.

Are you able to post any of your actual data? Without being able to reproduce it, it will be hard to figure out the problem. If you can, I'd recommend posting a new question (so you get more fresh eyes on it) with the code and the data included.

[Sean de Wolski]

I'd advise against decreasing tolerances like that. All it does it make the solver take longer. Instead, consider using multiple start points, MultiStart in Global Optim Toolbox does this for you.

doc multistart

[Renee]

thank you i took your advice. and to sean, i also said this in the post i just made but I am pretty new to matlab and dont have that toolbox (yet), I tried to use the curve fitting toolbox and it didnt seem to like my ode so I'm concerned it wont work (I never see examples with systems of ode and like I said I dont know what I'm doing to be able to reason it through).

[Matt Tearle]

If you make a new question, BTW, can you link it here, so I don't miss it? Thanks.

Also, have you tried just evaluating the ODE with some random parameters to see (a) how long it takes to run, and (b) what the results look like? That might help give a sense of what's happening with the curve-fitting part. Maybe there's a problem with how the ODEs are being evaluated; maybe the problem is with the curve-fitting part. Hard to tell without splitting them apart and checking.

[Renee]

hopefully this link will work: new post I'm working on the random numbers thing.

[Renee]

The output with random numbers looks the same.