Optimizing parameters in ODE

Question

Bibek Dhami 2019-9-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/481023-optimizing-parameters-in-ode

评论： Star Strider 2019-9-24

Hi I have a set of experimental data. I want to fit this experimental data to first order differential equation of the form dy/dt = -a*n-b*n^2-c*n^3 to optimize the value of constants a,b and c. Can anyone help in this regards? I am new to matlab as this question might be too simple for others. Thanks in advance.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Star Strider 2019-9-18

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/481023-optimizing-parameters-in-ode#answer_392510

在 MATLAB Online 中打开

This is a simple, separable differential equation that you can likely solve by hand.

Using the Symbolic Math Toolbox:

syms a b c n y(t) y0 
DEqn = diff(y) == -a*n-b*n^2-c*n^3;
Eqn = dsolve(DEqn, y(0)==y0)
fcn = matlabFunction(Eqn, 'Vars',{[a,b,c],t,n,y0})

produces:

Eqn =
y0 - t*(c*n^3 + b*n^2 + a*n)
 
fcn =
  function_handle with value:
    @(in1,t,n,y0) y0-t.*(in1(:,1).*n+in1(:,2).*n.^2+in1(:,3).*n.^3)

or more conveniently:

fcn = @(in1,t,n,y0) y0-t.*(in1(:,1).*n+in1(:,2).*n.^2+in1(:,3).*n.^3);

with ‘in1’ corresponding to [a,b,c] in that order. Supply values for ‘n’ and ‘y0’, then present it to the nonlinear parameter estimation function of your choice as:

objfcn = @(in1,t) fcn(in1,t,n,y0)

Or, since it is ‘linear in the parameters’ you can re-write it as a design matrix and use linear methods such as mldivide,\ to solve it as well.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Star Strider 2019-9-19

在 MATLAB Online 中打开

It would have been helpful to know that yesterday. Your differential equation is nonnlinear, so like most such, there is no symbolic solution to it. You will have to integrate it numerically in order to fit it to your data.

Try this:

    function N = objfcn(b,t)
        n0 = b(4);
        DifEq = @(t,n) -b(1).*n - b(2).*n.^2 - b(3).*n.^3;
        [t,N] = ode45(DifEq,t,n0);
    end
t = (0:9)';                                                     % Create Data
data = randn(size(t));                                          % Create Data
B0 = ones(1,4);
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@objfcn,B0,t,data);
prmnam = {'a','b','c','ic'};
fprintf(1,'\tParameters:\n')
for k = 1:length(B)
    fprintf(1, '\t\t%s\t = %23.15E\n', prmnam{k}, B(k))
end
tv = linspace(min(t), max(t));
Cfit = objfcn(B, tv);
figure(1)
plot(t, data, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Value')

This ran without error when I experimented with it.

Save the ‘objfcn’ function file to a separate file called objfcn.m on your MATLAB user path.

Remember that ‘data’ (the vector you want to fit) must be a column vector. Transpose it to be such if necessary.

You may need to experiment with different values of ‘B0’ (the initial parameter estimates) to get a good fit to your data. Nonnlinear parameter estimation is very sensitive to the initial parameter estimates, so several attempts may be necessary if you do not already have a good idea of what your parameters should be.

Bibek Dhami 2019-9-24

Thank you Star Strider for your code. It helped me a lot though I am stuck in the initial condition of parameters to exactly fit it.

Star Strider 2019-9-24

As always, my pleasure!

My code estimates the initial condition as well, estimating it as ‘b(4)’, in the printed results as ‘ic’. So an initial estimate for it shoulld be ‘B0(4)’.

请先登录，再进行评论。