An ode system solution as nlinfit function

5 次查看（过去 30 天）
Andrea Ristori
Andrea Ristori 2020-7-22
评论： Star Strider 2020-7-22
Hi everyone,
I am trying to use one of the solutions of a differential equations system to fit a set of data using nlinfit and nlparci since I found an example online. In particular, I think I am doing the iteration: solve the system -> obtain the solution -> use this solution to try to fit the data -> get better parameters from nlinfit and nlparci -> use this parameter to solve the system
and so on.
I would like to know a couple of thing:
  1. Does nlinfit solve the system every iteration, in order to obtain the function, or not?
  2. If not, how can I use the optimized parameters obtained from nlparci to solve again the system and to go on with the iteration?
Thank you for any help.
Here is my code:
Gamma_0 = 1.2;
gamma_0 = 0.3;
gamma_c_0 = 5.1;
c_0 = 6.5;
A_0 = 1;
t0_0 = -0.6;
y0_0 = 0;
parguess = [Gamma_0, gamma_0, gamma_c_0, c_0, A_0, t0_0, y0_0];
options = statset();
options.MaxIter = 1000;
[pars, resid, J] = nlinfit(tdata,ydata,@myfunc,parguess,options);
alpha = 0.05;
pars_co = nlparci(pars, resid, 'jacobian', J, 'alpha', alpha);
Gamma = pars(1);
gamma = pars(2);
gamma_c = pars(3);
c = pars(4);
A = pars(5);
t0 = pars(6);
y0 = pars(7);
parst = transpose(pars);
for i=1:length(pars)
errors(1,i) = max(abs(parst(i)-pars_co(i,1)),abs(pars(i)-pars_co(i,2)));
errors(2,i) = abs(errors(1,i)/parst(i))*100;
end
tfit = transpose(linspace(tdata(1),tdata(end),1000));
fit = myfunc(pars, tfit);
figure();
plot(tfit,fit,'b',tdata,ydata,'r');
legend('fit','data','FontSize',12);
figure();
semilogy(tfit,fit,'b',tdata,ydata,'r');
legend('fit','data','FontSize',12);
%%%%
function uscita = myfunc(pars, tdata)
Gamma = pars(1);
gamma = pars(2);
gamma_c = pars(3);
c = pars(4);
A = pars(5);
t0 = pars(6);
y0 = pars(7);
k = 760;
function dydt = ode(t,y)
dydt = zeros(8,1);
dydt(1)=-y(1)*gamma_c*2*(exp(-gamma_c*(t+t0)))+y(4)*gamma+y(8)*k; %p00
dydt(2)=y(1)*gamma_c*(exp(-gamma_c*(t+t0)))-y(2)*c*gamma_c*(exp(-gamma_c*(t+t0)))+y(5)*gamma; %p10
dydt(3)=y(1)*gamma_c*(exp(-gamma_c*(t+t0)))-y(3)*c*gamma_c*(exp(-gamma_c*(t+t0)))+y(6)*gamma; %p01
dydt(4)=-y(4)*(gamma+Gamma+gamma_c+2*(exp(-gamma_c*(t+t0))))+(y(2)+y(3))*c*gamma_c*(exp(-gamma_c*(t+t0)))+y(7)*gamma+y(8)*Gamma; %p11
dydt(5)=y(4)*gamma_c*(exp(-gamma_c*(t+t0)))-y(5)*(gamma+c*gamma_c*(exp(-gamma_c*(t+t0)))); %p21
dydt(6)=y(4)*gamma_c*(exp(-gamma_c*(t+t0)))-y(6)*(gamma+c*gamma_c*(exp(-gamma_c*(t+t0)))); %p12
dydt(7)=(y(5)+y(6))*c*gamma_c*(exp(-gamma_c*(t+t0)))-y(7)*gamma; %p22
dtdy(8)=y(4)*Gamma-y(8)*(k+Gamma); %pcm
end
y00 = 1;
y10 = 0;
y01 = 0;
y11 = 0;
y21 = 0;
y12 = 0;
y22 = 0;
ycm = 0;
[time,y] = ode45(@ode,tdata,[y00,y10,y01,y11,y21,y12,y22,ycm]);
uscita = A*y(:,4)+y0;
end

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

回答（1 个）

Star Strider
Star Strider 2020-7-22
‘Does nlinfit solve the system every iteration, in order to obtain the function, or not?’
It solves it at every iteration. It updates the ‘pars’ parameters and continues iterating until it converges on an acceptable parameter set. The nlparci results will be the parameter confidence intervals for the converged parameter set.
  4 个评论
显示 2更早的评论
Andrea Ristori
Andrea Ristori 2020-7-22
With "So my code is ok and it should work, is it?" I meant that the process is correct and it does the iteration I need to do.
Regarding y0, do you mean that is not re-defined in the sytem of equations? If it is so, yes I wanted only to intoduce a constant value to the fitting function (and an amplitute, A), to have an additional parameter to better fit the data.
Thank you again.
Star Strider
Star Strider 2020-7-22
As I wrote previously, I have no idea.
I have no idea what you are doing. If using ‘y0’ there is correct, then use it.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Computations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by