curve fit with ode bad fit ??

Hi, I'm a grad student with minimal matlab experience trying to complete this curve fit. My data has a logarithmic shape and is (hopefully) described by one ODE out of a system of 11 ODEs with 15 parameters. The only info I know is initial conditions of ODEs and some of the ratios of the parameters (so individual parameter values can vary greatly).

I have gotten my output to look logarithmic by manual manipulation of parameters (fig 1; data = blue, model = green circles), but when I try to get it to curve fit the output is always linear, usually intersecting in 2 places but sometimes not at all (fig 2: data = blue, model = green). (dont even know how its getting 10^22 with the bounds I specified)

I have heard of using global search as a possible solution, but I have the curve fitting toolbox and it doesnt like my system of ODE (probably because I dont know what im doing), so if someone can tell me that the global search toolbox will definitely work with a system of ODE then i'll give it a shot (because its going to take me forever lol).

curve fit code:

y = xlsread('dataforfitting.xlsx','B2:B1201');
t = xlsread('dataforfitting.xlsx','A2:A1201');
%init values for params
c0=[7000000000;.0000001;50000000;.000001;30000;.001;260000;.0001;.5;.05;.0000001;2;.1;450000000000000000000;.5];
%bounds for params
lb = [6000000000;.0000001;40000000;.000001;20000;.0001;200000;.00001;.4;.01;.00000001;1;.01;300000000000000000000;.1];
ub = [8000000000;.000001;500000000;.00001;300000;.01;300000;.001;.6;.1;.000001;5;1;550000000000000000000;1];
%curve fit
options = optimoptions(@lsqcurvefit,'ScaleProblem','Jacobian','MaxIter',20000,'MaxFunEvals',20000);
cfit = lsqcurvefit(@fitodes,c0,t,y,lb,ub,options);
% See results
plot(t,y,'o',t,fitodes(cfit,t))

ode code:

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

[Sean de Wolski]

First, your data has logarithmic behavior so I would advise using ode15s. Perhaps the optimizer is passing values that cause ode45 to fail.

[Renee]

I used ode15s for my manual "fit", but it doesnt seem to make a difference with lsqcurvefit (says its hitting the tolerance levels)

Answer 1

[Matt Tearle]

Ah. This may have been partly my fault for mixing notation (lots of t and x and y all floating around). A good check is to see what happens if you do

plot(t,y,t,fitodes(c0,t))

And you will see that you get the same graph that blows up to 10^21 right away. That's because of the muddled notation in fitodes. When doing a curve fit, the function should be of the form y = f(x,c) where x is the independent variable and c is the vector of parameters. In your case, the independent variable is the times in t.

In fitodes the independent variable is called x. That's fine, names don't matter, but... in the definition of the ODEs, f is defined as a function of t, y, and c. But the actual function definitions use x! (And recall that the x values here are actually the t values from your spreadsheet.) Confused yet? Bottom line: instead of solving an ODE like dy/dt = y, it's solving dy/dt = t.

The "best" solution would be to change all the "x"s in the ODE definition to "y"s. But that would be tedious. The much easier solution (although one that will leave you with some highly confusing code), is to simply change the names of the dummy variables in the definition of the function:

f = @(t,x,c) [c(1).*x(4)-c(2).*x(1)...
        ^

Having made that change, I ran with the c0 values (as above) and got a graph that looks not entirely unlike the data. However, to get it to run in any sensible time, I used ode23s with the default options.

It also barfed a lot of singularity warnings. That's not surprising. I noticed as I was looking at the code that there's some hideous scaling going on. My inner fluid dynamicist was twitching at the sight of parameters that are O(10^-9) along with other parameters that are O(10^20). Is there any way you can nondimensionalize your equations? This isn't a MATLAB issue, BTW. Any numerical solution is going to have a very hard time with equations spanning 30 orders of magnitude!

But, ignoring all that, I managed to get something that is almost working:

data = xlsread('dataforfitting.xlsx','A2:B1201');
t = data(:,1);
y = data(:,2);
figure
plot(t,y)

wstate = warning('query','MATLAB:nearlySingularMatrix');
warning('off','MATLAB:nearlySingularMatrix');

%init values for params
c0=[7000000000;.0000001;50000000;.000001;30000;.001;260000;.0001;.5;.05;.0000001;2;.1;450000000000000000000;.5];
%
figure
plot(t,y,t,fitodes(c0,t))
%
%bounds for params
lb = [6000000000;.0000001;40000000;.000001;20000;.0001;200000;.00001;.4;.01;.00000001;1;.01;300000000000000000000;.1];
ub = [8000000000;.000001;500000000;.00001;300000;.01;300000;.001;.6;.1;.000001;5;1;550000000000000000000;1];
%curve fit
cfit = lsqcurvefit(@fitodes,c0,t,y,lb,ub);%,options);
% See results
figure
plot(t,y,'o',t,fitodes(cfit,t))

warning(wstate.state,'MATLAB:nearlySingularMatrix');

with

function y = fitodes(c,x)
% Define ODE system with parameter c
f = @(t,x,c) [c(1).*x(4)-c(2).*x(1).*x(3)+c(3).*x(6)-c(4).*x(1).*x(5);
    c(5).*x(5)-c(6).*x(2).*x(3)+c(7).*x(6)-c(8).*x(4).*x(2);
    c(9).*c(10).*x(10).^c(12).*x(10).^c(12)-c(9).*c(11).*x(11)+c(5).*x(5)-c(6).*x(2).*x(3)+c(1).*x(4)-c(2).*x(1).*x(3);
    -c(1).*x(4)+c(2).*x(1).*x(3)+c(7).*x(6)-c(8).*x(4);
    c(6).*x(2).*x(3)-c(5).*x(5)+c(3).*x(6)-c(4).*x(1).*x(5);
    -c(7).*x(6)+c(8).*x(4)-c(3).*x(6)+c(4).*x(5).*x(1)-c(13).*x(6);
    c(14)*c(15).*x(6)-x(7);
    c(13).*x(6)*c(15);
    c(13).*x(6)*(1-c(15));
    -c(10).*x(10).^c(12).*x(10).^c(12)+c(11).*x(11);
    c(10).*x(10).^c(12).*x(10).^c(12)-c(11).*x(11);];
%init conditions
init =[20000;75;0;0;0;0;0;0;0;.45;0;];
%solve
[~,z] = ode23s(@(t,y) f(t,y,c),x,init);
% Extract desired solution component
y = z(:,7);

I say "almost working" because it hits the function eval limit and produces something that isn't a very good fit. But it's at least ballpark. So, it's on the right path. You can try tweaking options now, but I really think that the scaling is going to make your life miserable.

[Renee]

I couldn't improve on your solution so I'm going to try optimization now lol. I'm assuming by scaling you mean the large differences between my parameters? I can get them within a few orders of magnitude of each other but I'm always going to have one extremely large one unless I normalize my data somehow, but that wouldnt really work for my end goals (I have to compare data sets with increasing values)

new question