ODE and Data fitting

Question

Ehsan Homaee 2020-1-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/499509-ode-and-data-fitting

评论： Star Strider 2020-1-13

采纳的回答： Star Strider

I have 3 differential equations:

The initial condition for all the equations at timepoint=0 equal to 0. After solving this equaitons, I want to a data fitting on the equation below in order to find the nest fitted A, B, and C.

Can someone help me on this, because I face several errors when I tried to do it.

Maybe the best way is to find the u, v, and w functions and then place them in equations x and y; but how can I do it? is there any way that ode45 give the formula of u,v, and w?

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

请先登录，再进行评论。

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

Answer 1

Star Strider 2020-1-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/499509-ode-and-data-fitting#answer_409290

All I can do is to point you in the direction of successful efforts to do what you want.

One such is: Parameter Estimation for a System of Differential Equations and another is of course Monod kinetics and curve fitting.

We can’t help without seeing your code.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Star Strider 2020-1-9

在 MATLAB Online 中打开

The updated post clarified everything. This runs and gives a reasonably decent fit, however it may be necessary to use the optimoptions function to provide an options structure to allow lsqcurvefit a higher function evaluation limit. I assume the parameters are allowed to become negative (not always permitted), so I did not constrain them, since you did not in your code.

Try this:

    function Crodata1_Fit
        function y = Objfcn(b,t)
            x0 = b(end-2:end);
            
            [tv,xv] = ode45(@ODEs, t, x0);
            
            function dY = ODEs(T,Y)                     % Y(1) = v;  Y(2) = u,  Y(3) = w
                dY(1,:) = Y(1).*(-1.0./2.0e+1)+1.0./(Y(1).^2.*2.5e+1+Y(2).^2.*2.5e+1+2.5e+1);
                dY(2,:) = (Y(3).^2.*(3.0./5.0e+1))./(Y(3).^2+1.0)+(Y(2).^2.*(2.0./2.5e+1))./(Y(1).^2+Y(2).^2+1.0)+Y(2)./2.5e+1;
                dY(3,:) = Y(3).*(-3.0./2.5e+1)+1.0./(Y(1).^2.*2.5e+1+Y(2).^2.*2.5e+1+2.5e+1);
            end
            y = 1./(b(1) + b(2).*xv(:,2).^2 + b(3).*xv(:,1).^2);
        end
        
        tdata=[0	0.5	1	2	3	5	7.5	10	20	30	60	120	180].';
        Crodata1=[0.037373477
            2.938766958
            2.428716866
            4.673151154
            4.039997314
            7.286068277
            5.229206437
            6.722726692
            1.809539197
            1.975990902
            1.366188561
            1.103164147
            0.29209409];
        
        B0 = [rand(3,1); zeros(3,1)];
        
        [B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1);
        
        fprintf(1,'\tRate Constants:\n')
        for k1 = 1:length(B)
            fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
        end
        
        tv = linspace(min(tdata), max(tdata));
        Yfit = Objfcn(B, tv);
        
        figure(1)
        plot(tdata, Crodata1, 'p')
        hold on
        hlp = plot(tv, Yfit);
        hold off
        grid
        xlabel('Time')
        ylabel('Unknown Unit')
    end

I run all this code inside of a specific function I use for such purposes, so to make life easier, I named this one ‘Crodata1_Fit’. Save it as a separate function and run it by calling:

Crodata1_Fit

from a script or your Command Window.

Make appropriate changes to get the result you want.

One run produced this plot:

and these parameters:

		B(1) =  0.14574
		B(2) =  0.00926
		B(3) =  2.80074
		B(4) = -0.32271
		B(5) = -0.02885
		B(6) = -1.57333

Other optimisers (specifically the ga funciton) will likely produce even better results.

NOTE —

I used this code to derive the ‘ODEs’ equations:

syms u(t) v(t) w(t) t T Y 
DEs = [diff(u)==0.08*u^2/(1+u^2+v^2)+0.06*w^2/(1+w^2)+0.04*u; diff(v)==0.04*1/(1+u^2+v^2)-0.05*v; diff(w)==0.04*1/(1+u^2+v^2)-0.12*w];
[VF,Sbs] = odeToVectorField(DEs)
ODEs = matlabFunction(VF, 'Vars',{T,Y})

It is easier (and more accurate) to let the Symbolic Math Toolbox do the ‘heavy lifting’.

Star Strider 2020-1-13

Ehsan Homaee’s Answer moved here —

Thanks for your response. Why the output has 6 diffeent B? Since the objective function only has 3 parameters,what is the other 3 parameters?

Furthermore, if I want to restrict the parameters to only possitive values, how can I do it?

Star Strider 2020-1-13

在 MATLAB Online 中打开

As always, my pleasure!

The first 3 ‘B’ parameters correspond to ‘A’, ‘B’, and ‘C’. The last 3 are the initial conditions of the system of differential equations, since I always let the optimization function estimate them as well.

To restrict the first 3 to positive values, the lsqcurvefit call changes to:

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1, [zeros(1,3), -Inf(1,3)]);

I would not constrain the initial conditions. However if you want to constrain them as well to be positive:

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@Objfcn,B0,tdata,Crodata1,zeros(1,6));

I suspect the ga function would come up with the best parameter set, since it searches the entire parameter space and is not usually affected by local minima.

请先登录，再进行评论。

ODE and Data fitting

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

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

ODE and Data fitting

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

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

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

5 个评论
显示 3更早的评论隐藏 3更早的评论