Parameter Estimation for Differential Equation according to data measured
6 views (last 30 days)
Show older comments
Netanel shachak
on 11 Oct 2018
Commented: Star Strider
on 12 Oct 2018
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
0 Comments
Accepted Answer
Star Strider
on 11 Oct 2018
Edited: Star Strider
on 11 Oct 2018
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 Comments
Star Strider
on 12 Oct 2018
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.)
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!