odearguments error for fitting data with system of ODEs

Question

Rosie 2020-9-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/598723-odearguments-error-for-fitting-data-with-system-of-odes

评论： Star Strider 2020-9-30

采纳的回答： Star Strider

fittingdata_t_f.txt

Hi all,

I have some experimental data (attached) that I need to fit with the following model.

$dP/dt = k_nm(t)^{n_c} + k_mM(t) - k_aP(t)^2$

I want to fit the data with M(t), but I need to solve the whole system (obviously). Additionally, there are four parameters that I need to extract from the fit once I am done (

). I've tried to follow a number of previous posts to write the code:

https://www.mathworks.com/matlabcentral/answers/515421-fitting-coupled-ode-s-to-multiple-sets-of-experimental-data

https://www.mathworks.com/matlabcentral/answers/101883-how-do-i-estimate-or-optimize-the-parameters-of-my-ode-system-in-matlab-8-1-r2013a

https://www.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

This is what I have so far, but it's giving me the error below. I plan on plotting the curve resulting from the fit values y (I really only want y(2), or M(t)) and reporting the parameter p values, and I'm not even sure I'll be able to extract all of those results from this script.

Thanks for any help you can give! Let me know if you would like any more information to help figure this out.

p = lsqcurvefit(@MSolve,pguess,t,f);
% t = time
% f = raw data
    function y = MSolve(p,t)
        cond1 = 0;
        cond2 = 0;
        conds = [cond1, cond2];
       
        tspan = [t(1) t(end)];
       
        % y(1) = P(t), y(2) = M(t)
        [T,y] = ode45(@MSolveFit,tspan,conds);
       
        function dydt = MSolveFit(t,y)
            nc = 2; % # monomers for primary nucleation
            m0 = 10; % initial concentration of tau
            mt = m0 - y(2);
           
            % parameters :
            % p(1) = kn
            % p(2) = km
            % p(3) = ka
            % p(4) = kp
           
            % dydt(1) = dPdt(t)
            % dydt(2) = dMdt(T)
            dydt = [0 0];
            dydt(1) = p(1).*mt.^nc + p(2).*y(2) - p(3).*y(1).^2;
            dydt(2) = p(4).*mt.*y(1);          
        end
    end
    

Here's the error:

Error using odearguments (line 93)
MSOLVE/MSOLVEFIT must return a column vector.
Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in Yaofit_chi2>MSolve (line 638)
        [T,Y] = ode45(@MSolveFit,tspan,conds);
Error in lsqcurvefit (line 213)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in Yaofit_chi2 (line 625)
p = lsqcurvefit(@MSolve,pguess,t,f);

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

请先登录，再进行评论。

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

Answer 1

Star Strider 2020-9-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/598723-odearguments-error-for-fitting-data-with-system-of-odes#answer_499297

The ‘dydt’ definition should be:

dydt = [0; 0];

that would define a column vector. The differential equation runs correctrly with that change.

The problem is that you are fitting one column off data (‘f’), however your differential equation returns two columns. This throws:

Error using lsqcurvefit (line 269)
Function value and YDATA sizes are not equal.

That can be solved by using ‘t’ instead of ‘tspan’:

[T,y] = ode45(@MSolveFit,t,conds);

This works:

D = load('fittingdata_t_f.txt');
t = D(:,1);
f = D(:,2);
pguess = rand(4,1)*1E+6;
p = lsqcurvefit(@MSolve,pguess,t,f);
    function y = MSolve(p,t)
        cond1 = 0;
        cond2 = 0;
        conds = [cond1, cond2];
       
        tspan = [t(1) t(end)];
       
        % y(1) = P(t), y(2) = M(t)
        [T,y] = ode45(@MSolveFit,t,conds);
                
        y = y(:,2);
       
        function dydt = MSolveFit(t,y)
            nc = 2; % # monomers for primary nucleation
            m0 = 10; % initial concentration of tau
            mt = m0 - y(2);
           
            % parameters :
            % p(1) = kn
            % p(2) = km
            % p(3) = ka
            % p(4) = kp
           
            % dydt(1) = dPdt(t)
            % dydt(2) = dMdt(T)
            dydt = [0; 0];
            dydt(1) = p(1).*mt.^nc + p(2).*y(2) - p(3).*y(1).^2;
            dydt(2) = p(4).*mt.*y(1);          
        end
    end
figure
plot(t, f, 'p')
hold on
plot(t, MSolve(p,t), '-r')
hold off
grid