Selecting last value of tspan ode45
10 views (last 30 days)
Show older comments
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 Comments
Jan
on 6 Dec 2017
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."
Accepted Answer
Torsten
on 6 Dec 2017
Edited: Torsten
on 6 Dec 2017
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.
More Answers (1)
Jan
on 6 Dec 2017
Edited: Jan
on 6 Dec 2017
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 Comments
Jan
on 14 Dec 2017
@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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!