Ode system solution with unknown constant
显示 更早的评论
I'm trying to obtain and plot a function that adjusts to my data. I have a table with time, position and velocity, and an ode function which I have rewritten as a system:
ydot(1) = y(2)
ydot(2) = -a*y(2) - b*y(1) + c*cos(wt)
a,b,c and w are unknown parameters but I can estimate their aproximate values.
How can I solve the problem?
采纳的回答
Star Strider
2020-12-7
You have already seen and posted a Comment to Parameter Estimation for a System of Differential Equations and that (or similar threads) are what I would direct you to.
For your application, the ‘kinetics’ objective function in that solution would be:
function Y = kinetics(theta,t)
a = theta(1);
b = theta(2);
c = theta(3);
w = theta(4);
y0 = theta(5:6);
[t,yv] = ode45(@Difeq,t,y0);
function ydot = Difeq(t,y)
ydot = zeros(2,1);
ydot(1) = y(2);
ydot(2) = -a*y(2) - b*y(1) + c*cos(w*t);
end
Y = yv;
end
Note that here the code will estimate the intial conditions as well, so ‘theta0’ must be a 6-element vector.
I obviously cannot test this, so I am posting it as UNTESTED CODE. It should work, although it may require a bit of editing.
6 个评论
ana ramirez de las heras
2020-12-7
Star Strider
2020-12-7
My pleasure!
If I had your data, I could test it.
ana ramirez de las heras
2020-12-8
Star Strider
2020-12-9
My pleasure!
I could not get a good fit with lsqcurvefit, so I went with ‘Plan B’ and wrote ga (genetic algorithm) code for it. I added a phase term to the cos call (entirely appropriate, and does not assume that the cos term always begins at 0) and while the fit improved slightly, it did so only for position (matching the frequency and phase, however not the amplitude), not velocity, that remains difficult (essentially impossible) to fit.
Using this ‘solver’ code with the added phase term:
%Calcular los parametros que se ajustan a los datos
function Y=solver(theta,t)
a = theta(1);
b = theta(2);
c = theta(3);
w = theta(4);
p = theta(5); % Phase
y0 = theta(6:7);
[t,yv] = ode45(@Difeq,t,y0);
function ydot = Difeq(t,y)
ydot = zeros(2,1);
ydot(1) = y(2);
ydot(2) = -a*y(2) - b*y(1) + c*cos(w*t + p);
end
Y = yv;
end
likely the best parameter set ga could estimate were:
Rate Constants:
Theta(1) = 4.22996
Theta(2) = 6.33071
Theta(3) = 19.42851
Theta(4) = 2.01714
Theta(5) = 5.17042
Theta(6) = 0.69420
Theta(7) = 2.37232
with a fitness value of 63.9.
If you have the Global Optimization Toolbox, I will post my ga code for your problem.
Star Strider
2020-12-9
My pleasure!
I discovered that the columns of ‘y’ are incorrect, and should be reversed, the correct orientation being:
y=[v p];
since the derivatrive is the first column of ‘Y’ and its integral is the second column. (I was too tired yesterday to detect that error in your code.)
Make that change, and the fit is excellent:
Rate Constants:
Theta(1) = 0.04702
Theta(2) = 9.99643
Theta(3) = 20.76383
Theta(4) = 1.99939
Theta(5) = 3.17972
Theta(6) = 0.94858
Theta(7) = 1.39551
The fitness value (norm of the residuals) for this run was typical at 13.6113.
I plotted the two variables separately, since it is difficult to see them when all are plotted together.
Star Strider
2020-12-9
As always, my pleasure!
I typically let the optimisation function calculate everything, in order not to constrain it beyond any obvious limits (for example constraining ‘w’ and ‘p’ to 
here).
here). 更多回答（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 
