MATLAB Answers

odearguments error for fitting data with system of ODEs

17 views (last 30 days)
Rosie
Rosie on 23 Sep 2020 at 19:43
Commented: Star Strider on 30 Sep 2020 at 15:33
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
dM/dt = k_pm(t)P(t)
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 (k_n, k_m, k_a, and k_p). I've tried to follow a number of previous posts to write the code:
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 Comments

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 23 Sep 2020 at 20:15
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
However I cannot get a good fit to your data, probably because I have no idea what you are doing or what the initial parammeter estimates should be.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by