Parameter Estimation for Differential Equation according to data measured

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

Answer 1

[Star Strider]

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.

[Netanel shachak]

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]

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.)