清空筛选条件
清空筛选条件

Parameter Estimation for Differential Equation according to data measured

7 次查看（过去 30 天）
Netanel shachak
Netanel shachak 2018-10-11
评论： Star Strider 2018-10-12
Thanks in advance to all supporters, this is the first time a ask a question - that's because there are answers to most of the problems i had.
After reviewing the proposed solution :
  • i tried to apply "Star Strider" ‘Igor_Moura’ function to another ode problem i have but get Error "using lsqcurvefit"
  • i did a modification to your code to solve the eq:
  • A,B,C are parameters i look for ([A B C] inside the code are theta(1:3). [y,y'] are c(1:2), dcdt(2) is [y''] and x is t)
function C=kinetics(theta,t)
c0=[293;0]; %initial condition
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(2,1);
dcdt(1)=c(2);
dcdt(2)=theta(3)-theta(2).*c(1)-theta(1).*t.*c(2);
dC=dcdt;
end
C=Cv;
end
global t c
load('t.mat'); t=t*1000;
load('c.mat');
theta0=[1;1;1]; % initial guess
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
  • unfortunately, I encounter a problem of dimensions :
Error using lsqcurvefit (line 248)
Function value and YDATA sizes are not equal.
  • at line 199 :
  • initVals.F get the size :
  • so at line 248 :
  • i got the err because YDATA length is 18896X1 so their sizes are not equal.
any idea?? thanks Netanel

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2018-10-11
编辑：Star Strider 2018-10-11
I will now delete your Comment from the other post.
I found the problem, and your code now works. Your ‘kinetics’ function should be:
function C=kinetics(theta,t)
c0=[0;293]; %initial condition
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(2,1);
dcdt(1)=c(2);
dcdt(2)=theta(3)-theta(2).*c(1)-theta(1).*t.*c(2);
dC=dcdt;
end
C=Cv(:,2);
end
Your parameters are very difficult to estimate. After failing to get a good fit with lsqcurvefit, I decided to let ga (Genetic Algorithm) estimate them. That has now been going for 11 hours, and does not appear close to converging.
EDIT — (11 Oct 2018 at 20:16)
The best parameters I can estimate (unconstrained genetic algorithm) are:
Rate Constants:
Theta(1) = 98063.79853
Theta(2) = 99458.01337
Theta(3) = 97298.39633
They do not come close to approximating your data, and have a residual norm of 1543.5.
I suspect that your model is not appropriate to your data.
  2 个评论
Netanel shachak
Netanel shachak 2018-10-12
Thanks Professor!
Strange though.
i got different theta (same code above):
Rate Constants:
Theta(1) = 15.53345
Theta(2) = -15.55330
Theta(3) = 198.99352
it looks not bad:
but i do not understand basic things:
as i know - ode45 returns:
[T,Cv]=ode45(@DifEq,t,c0);
while The first column of Cv corresponds to c(1) means y(x) from the original equation , and the second column to c(2) means y'(x). if so - why you correct the code :
C=Cv(:,2);
that's means the function "kinetics" returns the first derivative of y?
the same confusion with the initial condition: you correct to :
c0=[0;293]; %initial condition
while i thought that c0(1) is y(0) and c0(2) is y'(0). my initial condition (as you can see at the plot) is nearly y(0)=293;
i expect to : c0=[293;0]
last thing: the homogeneous equetion is:
i insert a step input f(t) to my system above (real experiment-measuring temperature using Hot-Wire technic, that's why the initial y(0) is 293 Kelvin (room temperature) and i insert my prob to environment with 355 Kelvin:
f(t)=355; means y(x--->inf)=355 you can see that the result at the plot agree with that BUT the coefficient theta(3) which is corresponds to my f(t) is fount to be
Theta(3) = 198.99352
any idea?
Thank you very much for your willingness to help
Netanel
Star Strider
Star Strider 2018-10-12
As always, my pleasure!
‘... that's means the function "kinetics" returns the first derivative of y?’
No. It is returning the integrated solution for ‘y’. The first derivative would be ‘Cv(:,1)’.
‘... while i thought that c0(1) is y(0) and c0(2) is y'(0).’
You simply have them reversed.
‘... any idea?’
Unfortunately, no. (If I think of something, I will post back here.)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Plot Settings 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by