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!
MATLAB: Curve fit toolbox help system ODE and excel
curve fittingexcelodesystemtoolbox
Related Solutions
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);figureplot(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];%
figureplot(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
figureplot(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.
‘What am I doing wrong?’
First, not posting the actual problem, that being:
Warning: Failure at t=7.149052e+00. Unable to meet integration tolerances withoutreducing the step size below the smallest value allowed (1.421085e-14) at time t. so the size discrepancy is due to ode45 not being able to estimate the values for all the time elements.
Second, setting ‘c0’ to have only 3 elements when 5 are required. Fixing two of the parameters within the objective function is not appropriate. A better approach is to constrain them using the parameter bounds vectors.
Note that ‘c0’ needs to have realistic values. All the elements of ‘Sdata’ are between 1 and 17, so the last 3 elements of ‘c0’ in the slightly revised code below needs to reflect that.
I changed the code slightly to reflect this, and also to estimate the initial conditions as part of the other parameters. I leave you to experiment with it:
Time = [0 1 2 3 4 5 6 7 8 9];Sdata = [3 17 9 4 8 2 1 9 2 4];% Sdata = Sdata';
figureplot(Time, Sdata)gridc0 = rand(8,1)*10;[c,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Fit,c0,Time,Sdata);function S = Fit(c, t, ~)y0 = c(6:end); [T,Sv] = ode45(@DifEq, t, y0); function dS = DifEq(t, y) ydot = zeros(3,1); ydot(1) = (c(1)*y(1).^2) + (c(2)*(y(2))) - (c(3)*y(3)); ydot(2) = (c(3)*y(3)) + (c(4)*y(1)); ydot(3) = (c(1)*y(1).^2) - (c(5)*y(2)) + (c(3)*y(3)); dS = ydot; endS = Sv(:,3);endThe problem appears to be that the differential equations may not correctly describe the system, or that the initial parameter estimates are inappropirate. Also, the system appears to be stiff, so ode15s would be an appropriate solver.
Estimating the parameters using the Global Optimization Toolbox ga (genetic algorithm) function would likely provide the best results, since lsqcurvefit requires that the initial parameter estimates are reasonably close to the actual optimal values.
Best Answer