Selecting last value of tspan ode45

조회 수: 9(최근 30일)
gorilla3
gorilla3 2017년 12월 6일
댓글: 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
gorilla3
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.

댓글을 달려면 로그인하십시오.

채택된 답변

Torsten
Torsten 2017년 12월 6일
편집: 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
gorilla3
gorilla3 2017년 12월 14일
Wonderful, thanks a lot Torsten!

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Jan
Jan 2017년 12월 6일
편집: Jan 2017년 12월 6일
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
Jan
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.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by