Selecting last value of tspan ode45
조회 수: 12 (최근 30일)
이전 댓글 표시
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."
채택된 답변
추가 답변 (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
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!