# Selecting last value of tspan ode45

조회 수: 25(최근 30일)
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표시숨기기 이전 댓글 수: 1
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 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표시숨기기 이전 댓글 수: 7
gorilla3 2017년 12월 14일
Wonderful, thanks a lot Torsten!

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

### 추가 답변(1개)

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표시숨기기 이전 댓글 수: 2
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.

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

### 범주

Find more on Physics in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by