Selecting last value of tspan ode45

Question

gorilla3 2017-12-6

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/371408-selecting-last-value-of-tspan-ode45

评论： Jan 2017-12-14

I'm trying to extract results only for the last value of the tspan in ode45. However, I end up with 2 results for tfinal... I don't understand why this is happening.

Here's the code:

PS. The ode45 is ran in a for loop because I want to solve the set of diff eqs for 130 different values of ABP (which is a parameter in the eqs).

clear all 
clc
%%=========== ODE45 ==================   
ABP= linspace(40,170,131);
     for i=1:1:length(ABP)  %change in ABP at every loop 
         sol = ode45(@first,[0 100], [0 0 0 0 0.205 97.6 12 49.67],[],ABP(i)); %(function, timespan, initial condition for xq,xc,xm1,xm,Ca,P1,V_sa,P2)
     end
%%========= FUNCTION ==================
function dvdt = first(t,v,ABP)
xq= v(1); %variables
xc= v(2);
xm= v(3);
xm1= v(4);
Ca=v(5);
P1= v(6);
V_sa= v(7);
P2 = v(8);
   % -----Constants -----
  R_la= 0.4;
  R_sa_b= 5.03;
  R_sv= 1.32;
  R_lv= 0.56;
  P_v= 6;
  V_la=1;
  V_sa_b= 12;
  P_ic= 10;
  k_ven= 0.186;
  P_v1= -2.25;
  V_vn= 28; 
  tau_q= 20; 
  Pa_co2_b= 40; 
  tau_co2= 40; 
  tau1= 2; 
  tau2= 4;
  tau_g= 1; 
  C_a_p=2.87; 
  C_a_n= 0.164;
  g=0.2; 
  E=0.4; 
  K= 0.15; 
  V0= 0.02; 
  q_b=12.5;
  G_q=3; 
  Pa_co2=40; 
  Ca_b= 0.205;
  eps=1;
  u=0.5; 
  Pa_b=100;
  ka=3.68;
deltaCa_p=2.87;
deltaCa_n=0.164;
% =========== System of diff eq =========================================
dxq= (-xq+G_q*( ( ( (P1-P2)/(R_sv+0.5*(R_sa_b*(V_sa_b/V_sa)^2))) -q_b) /q_b) )/tau_q;
dxc=(-xc +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
dxm=(eps*u-tau2*xm1-xm)/tau1^2;
dxm1=xm1;
sum_xqxcxm=xm+xc-xq; %sum of feedback mechanisms
if t==100 % Last value of tspan!
    disp(sum_xqxcxm) %should get 14 displayed values (because i=14) but I get 28!
end
   if (ABP==40)||(ABP==41) ||(ABP==42) ||(ABP==43) ||(ABP==44) ||(ABP==45) ||(ABP==46) ||(ABP==47)||(ABP==48)||(ABP==49)||(ABP==50)||(ABP==51)||(ABP==52)
      delta=deltaCa_p; %because sum_xqxcxm >0 for ABP=40
elseif (ABP==53)||(ABP==54) %add other values later
     delta=deltaCa_n; %sum<0
 end
dCa=0.5*delta*(1-tanh(2*sum_xqxcxm/delta)^2);
dP1= 1/Ca * ((ABP-P1)/(R_la+0.5*(R_sa_b*(V_sa_b/V_sa)^2)) - (P1-P2)/(R_sv+0.5*(R_sa_b*(V_sa_b/V_sa)^2)) -dCa*(P1-P_ic));
dV_sa= dCa*(P1-P_ic) + Ca*dP1;
dP2=1/(1/(k_ven*(P2-P_ic-P_v1))) *((P1-P2)/(R_sv+0.5*(R_sa_b*(V_sa_b/V_sa)^2)) -(P2-P_v)/R_lv);
dvdt=[dxq;dxc;dxm;dxm1;dCa;dP1;dV_sa;dP2];
%%==== print to file====
Rsa= R_sa_b*(V_sa_b/V_sa)^2;
CBF= (P1-P2)/(Rsa*0.5 + R_sv);
CBF_ch= (CBF-q_b)/q_b *100;
if t==100 %save only final solution
    fileID=fopen('results.txt','a');
     fprintf(fileID,' %-5s %-2s %-5s %-5s %-5s %-5s %-5s %-5s %-5s %-5s\n','xq','xc','xm','xm1','Ca','P1','V_sa','P2', 'CBF','CBFchange'); %how to write this only as a header and not repeat every time?
fprintf(fileID,' %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f\n\n',xq,xc,xm,xm1,Ca,P1,V_sa,P2,CBF,CBF_ch)
fclose(fileID);
end
end

2 个评论
显示无隐藏无

Jan 2017-12-6

Note: The brute clearing of everything "clear all" is rarely useful, but wastes a lot of time with reloading all functions from the disk.

Where do you "end up with 2 results" wand what do you expect instead? Please explain what "this" is in: "I don't understand why this is happening."

gorilla3 2017-12-6

if you look in the code, you will see that I put a comment next to disp(sum_xqxcxm). When I display it only 1 value should be displayed, instead I get 2 values.

请先登录，再进行评论。

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

Answer 1

Torsten 2017-12-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/371408-selecting-last-value-of-tspan-ode45#answer_294967

编辑：Torsten 2017-12-6

Function "first" is called several times for the same t-values.

Better call the function "first" once more after ode45 has finished (i.e. after sol = ode45(...)). Or use OutputFcn.

Best wishes

Torsten.

8 个评论
显示 6更早的评论隐藏 6更早的评论

gorilla3 2017-12-14

在 MATLAB Online 中打开

I have tried using your method but I keep getting errors:

ABP= linspace(40,170,131);
t=[0 100];
       for i=1:1:length(ABP)  %change in ABP at every loop        
           [T,Y] = ode45(@first,t, [0 0 0 0 0.205 62.516 12 13.201],[],ABP(i)); %(function, timespan, initial condition for xq,xc,xm1,xm,Ca,P1,V_sa,P2)
           t=T(end);
             y=Y(end,:);
             dy=first2(t,y);
       end
  %%========= FUNCTION ==================
  function dvdt = first(t,v,ABP)
xq= v(1); %variables
xc= v(2);
xm= v(3);
xm1= v(4);
Ca=v(5);
P1= v(6);
Vsa= v(7);
P2 = v(8);
     % -----Constants -----
    Rla= 0.4;
    Rsab= 5.03;
    Rsv= 1.32;
    Rlv= 0.56;
    Pv= 6;
    Vla=1;
    Vsab= 12;
    Pic= 10;
    kven= 0.186;
    Pv1= -2.25;
    Vvn= 28; 
    tauq= 20; 
    Pcb= 40; 
    tauc= 40; 
    tau1= 2; 
    tau2= 4;
    taug= 1; 
    g=0.2; 
    E=0.4; 
    K= 0.15; 
    V0= 0.02; 
  %   qb=12.5;
  qb=12.5;
    Gq=3; 
    Pc=40; 
    Ca_b= 0.205;
    eps=1;
    u=0.5; 
    Pab=100;
    ka=3.68;
deltap=2.87;
deltan=0.164;
% =========== System of diff eq =========================================
Rsa= Rsab/((Vsa/Vsab)^2);
q=(P1-P2)/(Rsv+0.5*Rsa);
Cv=1/(kven*(P2-Pic-Pv1));
fn=0.3+3*tanh(Pc/Pcb -1.1);
dxq= 1/tauq *(-xq+Gq*((q-qb)/qb));
dxc=1/tauc*(-xc +fn);
dxm=1/tau2*(-tau1^2*xm1-xm);
dxm1=xm1;
x=xm+xc-xq; %sum of feedback mechanisms
dx=dxm+dxc-dxq;
if x>=0
    delta=deltap;
else
    delta=deltan;
end
dCa=dx* 1/(cosh(2*x/delta))^2;
dP1=1/Ca *((ABP-P1)/(Rla+Rsa*0.5) -q - dCa*(P1-Pic));
dVsa= dCa*(P1-Pic)+Ca*dP1;
dP2=1/Cv *(q-(P2-Pv)/Rlv);
dvdt=[dxq;dxc;dxm;dxm1;dCa;dP1;dVsa;dP2];
end
%%========= FUNCTION 2 ==================
function dy = first2(t,v,ABP)
xq= v(1); %variables
xc= v(2);
xm= v(3);
xm1= v(4);
Ca=v(5);
P1= v(6);
Vsa= v(7);
P2 = v(8);
     % -----Constants -----
    Rla= 0.4;
    Rsab= 5.03;
    Rsv= 1.32;
    Rlv= 0.56;
    Pv= 6;
    Vla=1;
    Vsab= 12;
    Pic= 10;
    kven= 0.186;
    Pv1= -2.25;
    Vvn= 28; 
    tauq= 20; 
    Pcb= 40; 
    tauc= 40; 
    tau1= 2; 
    tau2= 4;
    taug= 1; 
    g=0.2; 
    E=0.4; 
    K= 0.15; 
    V0= 0.02; 
  %   qb=12.5;
  qb=12.5;
    Gq=3; 
    Pc=40; 
    Ca_b= 0.205;
    eps=1;
    u=0.5; 
    Pab=100;
    ka=3.68;
deltap=2.87;
deltan=0.164;
% =========== System of diff eq =========================================
Rsa= Rsab/((Vsa/Vsab)^2);
q=(P1-P2)/(Rsv+0.5*Rsa);
Cv=1/(kven*(P2-Pic-Pv1));
fn=0.3+3*tanh(Pc/Pcb -1.1);
dxq= 1/tauq *(-xq+Gq*((q-qb)/qb));
dxc=1/tauc*(-xc +fn);
dxm=1/tau2*(-tau1^2*xm1-xm);
dxm1=xm1;
x=xm+xc-xq; %sum of feedback mechanisms
dx=dxm+dxc-dxq;
if x>=0
    delta=deltap;
else
    delta=deltan;
end
dCa=dx* 1/(cosh(2*x/delta))^2;
dP1=1/Ca *((ABP-P1)/(Rla+Rsa*0.5) -q - dCa*(P1-Pic));
dVsa= dCa*(P1-Pic)+Ca*dP1;
dP2=1/Cv *(q-(P2-Pv)/Rlv);
dy=[dxq;dxc;dxm;dxm1;dCa;dP1;dVsa;dP2];
qch= (q-qb)/qb *100;
figure(1)
  if t==100
      Vv= (1/kven)*log(P2-Pic-Pv1)+Vvn;
       if t==100 %save only final solution
    fileID=fopen('venousV.txt','a');
fprintf(fileID,' %4.3f\n',Vv);
fclose(fileID);
fileID=fopen('Vsa.txt','a');
fprintf(fileID,' %4.3f\n',Vsa);
fclose(fileID);
       end
  end

I'm making the plots only in first2 now, as suggested.

The error is:
Not enough input arguments.
Error in loop1412>first2 (line 158)
dP1=1/Ca *((ABP-P1)/(Rla+Rsa*0.5) -q - dCa*(P1-Pic));
Error in loop1412 (line 12)
           dy=first2(t,y);

Torsten 2017-12-14

编辑：Torsten 2017-12-14

在 MATLAB Online 中打开

dy = first2(t,y,ABP(i));

And open your output files in the main program.

And remove the if-statement concerning t=100 in first2 ; they are no longer needed.

Best wishes

Torsten.

gorilla3 2017-12-14

Wonderful, thanks a lot Torsten!

请先登录，再进行评论。

Answer 2

Jan 2017-12-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/371408-selecting-last-value-of-tspan-ode45#answer_294968

编辑：Jan 2017-12-6

在 MATLAB Online 中打开

Do not use the parameter as 5th input, because this is deprecated for over 17 years now. See http://www.mathworks.com/matlabcentral/answers/1971.

A bold guess: You overwrite sol in each iteration. Here is one way to collect them instead:

solC = cell(1, length(ABP));
yIni = [0 0 0 0 0.205 97.6 12 49.67];
for k = 1:length(ABP)
   sol{k} = ode45(@(t, y) first(t, y, ABP(k)), [0 100], yIni);
end

Or maybe you want this:

yR = zeros(length(ABP), 8);
for k = 1:length(ABP)
   [T, Y]   = ode45(@(t, y) first(t, y, ABP(k)), [0 100], yIni);
   yR(k, :) = Y(end, :);
end

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

gorilla3 2017-12-13

Here's a better explanation of the issue:

https://uk.mathworks.com/matlabcentral/answers/372804-for-loop-ode-45-changing-parameter-values

Jan 2017-12-14

@gorilla3: As you have found out, "sol" is a typo. Call it "solC" instead. I think, you should be able to fix this by your own, and if the result is provided in a variable called "sol", you should be able to find it there also.

请先登录，再进行评论。

Selecting last value of tspan ode45

2 个评论
显示无隐藏无

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（1 个）

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

另请参阅

类别

标签

Community Treasure Hunt

Selecting last value of tspan ode45

2 个评论 显示 无隐藏 无

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（1 个）

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

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

8 个评论
显示 6更早的评论隐藏 6更早的评论

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