ODE and Data fitting
显示 更早的评论
I have 3 differential equations:
The initial condition for all the equations at timepoint=0 equal to 0. After solving this equaitons, I want to a data fitting on the equation below in order to find the nest fitted A, B, and C.
Can someone help me on this, because I face several errors when I tried to do it.
Maybe the best way is to find the u, v, and w functions and then place them in equations x and y; but how can I do it? is there any way that ode45 give the formula of u,v, and w?
采纳的回答
Star Strider
2020-1-8
0 个投票
All I can do is to point you in the direction of successful efforts to do what you want.
One such is: Parameter Estimation for a System of Differential Equations and another is of course Monod kinetics and curve fitting.
We can’t help without seeing your code.
5 个评论
Ehsan Homaee
2020-1-9
Star Strider
2020-1-9
编辑：Star Strider
2020-1-9
First, do not use global variables! They create code that is almost impossible to debug.
Second, you apparently have three differential equations (I am going on your original Question, since I cannot decipher your code) and you are fitting one variable. I have no idea which one that is.
My ‘Monod kinetics’ code does something similar, so you can use that as a prototype.
Star Strider
2020-1-9
The updated post clarified everything. This runs and gives a reasonably decent fit, however it may be necessary to use the optimoptions function to provide an options structure to allow lsqcurvefit a higher function evaluation limit. I assume the parameters are allowed to become negative (not always permitted), so I did not constrain them, since you did not in your code.
Try this:
function Crodata1_Fit
function y = Objfcn(b,t)
x0 = b(end-2:end);
[tv,xv] = ode45(@ODEs, t, x0);
function dY = ODEs(T,Y) % Y(1) = v; Y(2) = u, Y(3) = w
dY(1,:) = Y(1).*(-1.0./2.0e+1)+1.0./(Y(1).^2.*2.5e+1+Y(2).^2.*2.5e+1+2.5e+1);
dY(2,:) = (Y(3).^2.*(3.0./5.0e+1))./(Y(3).^2+1.0)+(Y(2).^2.*(2.0./2.5e+1))./(Y(1).^2+Y(2).^2+1.0)+Y(2)./2.5e+1;
dY(3,:) = Y(3).*(-3.0./2.5e+1)+1.0./(Y(1).^2.*2.5e+1+Y(2).^2.*2.5e+1+2.5e+1);
end
y = 1./(b(1) + b(2).*xv(:,2).^2 + b(3).*xv(:,1).^2);
end
tdata=[0 0.5 1 2 3 5 7.5 10 20 30 60 120 180].';
Crodata1=[0.037373477
2.938766958
2.428716866
4.673151154
4.039997314
7.286068277
5.229206437
6.722726692
1.809539197
1.975990902
1.366188561
1.103164147
0.29209409];
B0 = [rand(3,1); zeros(3,1)];
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(B)
fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
end
tv = linspace(min(tdata), max(tdata));
Yfit = Objfcn(B, tv);
figure(1)
plot(tdata, Crodata1, 'p')
hold on
hlp = plot(tv, Yfit);
hold off
grid
xlabel('Time')
ylabel('Unknown Unit')
end
I run all this code inside of a specific function I use for such purposes, so to make life easier, I named this one ‘Crodata1_Fit’. Save it as a separate function and run it by calling:
Crodata1_Fit
from a script or your Command Window.
Make appropriate changes to get the result you want.
One run produced this plot:
and these parameters:
B(1) = 0.14574
B(2) = 0.00926
B(3) = 2.80074
B(4) = -0.32271
B(5) = -0.02885
B(6) = -1.57333
Other optimisers (specifically the ga funciton) will likely produce even better results.
NOTE —
I used this code to derive the ‘ODEs’ equations:
syms u(t) v(t) w(t) t T Y
DEs = [diff(u)==0.08*u^2/(1+u^2+v^2)+0.06*w^2/(1+w^2)+0.04*u; diff(v)==0.04*1/(1+u^2+v^2)-0.05*v; diff(w)==0.04*1/(1+u^2+v^2)-0.12*w];
[VF,Sbs] = odeToVectorField(DEs)
ODEs = matlabFunction(VF, 'Vars',{T,Y})
It is easier (and more accurate) to let the Symbolic Math Toolbox do the ‘heavy lifting’.
Star Strider
2020-1-13
Ehsan Homaee’s Answer moved here —
Thanks for your response. Why the output has 6 diffeent B? Since the objective function only has 3 parameters,what is the other 3 parameters?
Furthermore, if I want to restrict the parameters to only possitive values, how can I do it?
Star Strider
2020-1-13
As always, my pleasure!
The first 3 ‘B’ parameters correspond to ‘A’, ‘B’, and ‘C’. The last 3 are the initial conditions of the system of differential equations, since I always let the optimization function estimate them as well.
To restrict the first 3 to positive values, the lsqcurvefit call changes to:
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1, [zeros(1,3), -Inf(1,3)]);
I would not constrain the initial conditions. However if you want to constrain them as well to be positive:
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1,zeros(1,6));
I suspect the ga function would come up with the best parameter set, since it searches the entire parameter space and is not usually affected by local minima.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
