ode45 for loop stuck at 1st loop

2 次查看（过去 30 天）
gorilla3
gorilla3 2017-11-29
评论： Torsten 2017-11-29
I'm trying to solve an ode45 for 130 different values of ABP (a parameter in the equations), hence I decided to make a for loop around the solver but the solver is stuck in the first value of ABP (ABP(i)=40, i=1) Could you please tell me how to make it shift to i=2, hence ABP=41 ? Here's the full code so you can see what's happening by running it.
clear all
clc
%%=========== ODE45 ==================
ABP= linspace(40,170,131) %arterial pressure
for i=1:1:length(ABP) %change in ABP at every time step
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,xm,xm1,Ca,P1,V_sa,P2 ; options ; call new value of ABP)
end
%%========= FUNCTION ==================
function dvdt = first(t,v,ABP)
xq= v(1);
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/V_sa_b)^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
%----- AMPLITUDE OF COMPLIANCE -----
if (ABP==40)
delta=deltaCa_n; %because sum_xqxcxm <0 for ABP=40
elseif (ABP==41)
delta=deltaCa_n;
end
%----- ARTERIAL COMPLIANCE -----
if ABP==100 %baseline condition
Ca=Ca_b
elseif (ABP<100)
Ca= Ca_b+0.5*delta*tanh(2*sum_xqxcxm/delta)
elseif(ABP>100)
Ca= Ca_b+0.5*delta*tanh(2*sum_xqxcxm/delta)
end
dCa=0.5*delta*(1- ((cosh(4*(xm+xc-xq)/delta))-1)/(cosh(4*(xm+xc-xq)) +1));
dP1= 1/Ca * ((ABP-P1)/(R_la+0.5*(R_sa_b/(V_sa/V_sa_b)^2)) - (P1-P2)/(R_sv+0.5*(R_sa_b/(V_sa/V_sa_b)^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/V_sa_b)^2)) -(P2-P_v)/R_lv);
dvdt=[dxq;dxc;dxm;dxm1;dCa;dP1;dV_sa;dP2]
%Calculate CBF
Rsa= R_sa_b*(V_sa_b/V_sa)^2;
CBF= (P1 -P2)/(Rsa*0.5 + R_sv);
figure(1)
xlabel('ABP')
ylabel('CBF')
title('Cerebral blood flow dependence on arterial blood pressure')
plot(ABP,CBF)
hold on
end
  6 个评论
显示 4更早的评论
gorilla3
gorilla3 2017-11-29
Thanks but that still doesn't solve the issue of the for loop not moving forward. I get the same output as in the original code I posted above.
Torsten
Torsten 2017-11-29
Before you use the loop, you should be almost sure that the solver manages to integrate your system for the parameters given.
So first test for single values of ABP what the problem is. A good starting point would thus be abp = 40.
Best wishes
Torsten.

请先登录，再进行评论。

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

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by