error in runge kutta -Unable to perform assignment because the left and right sides have a different number of elements.

10 Ansichten (letzte 30 Tage)
Ayush Ranjan
Ayush Ranjan am 12 Dez. 2024 um 10:59
Bearbeitet: Torsten am 12 Dez. 2024 um 12:16
Ran in:
h=10;
numberofsteps=500/h;
k1_1=zeros(numberofsteps,1);
k2_1=zeros(numberofsteps,1);
k3_1=zeros(numberofsteps,1);
k4_1=zeros(numberofsteps,1);
k1_2=zeros(numberofsteps,1);
k2_2=zeros(numberofsteps,1);
k3_2=zeros(numberofsteps,1);
k4_2=zeros(numberofsteps,1);
k1_3=zeros(numberofsteps,1);
k2_3=zeros(numberofsteps,1);
k3_3=zeros(numberofsteps,1);
k4_3=zeros(numberofsteps,1);
k1_4=zeros(numberofsteps,1);
k2_4=zeros(numberofsteps,1);
k3_4=zeros(numberofsteps,1);
k4_4=zeros(numberofsteps,1);
k1_5=zeros(numberofsteps,1);
k2_5=zeros(numberofsteps,1);
k3_5=zeros(numberofsteps,1);
k4_5=zeros(numberofsteps,1);
k1_6=zeros(numberofsteps,1);
k2_6=zeros(numberofsteps,1);
k3_6=zeros(numberofsteps,1);
k4_6=zeros(numberofsteps,1);
k1_7=zeros(numberofsteps,1);
k2_7=zeros(numberofsteps,1);
k3_7=zeros(numberofsteps,1);
k4_7=zeros(numberofsteps,1);
k1_8=zeros(numberofsteps,1);
k2_8=zeros(numberofsteps,1);
k3_8=zeros(numberofsteps,1);
k4_8=zeros(numberofsteps,1);
k1_9=zeros(numberofsteps,1);
k2_9=zeros(numberofsteps,1);
k3_9=zeros(numberofsteps,1);
k4_9=zeros(numberofsteps,1);
k1_10=zeros(numberofsteps,1);
k2_10=zeros(numberofsteps,1);
k3_10=zeros(numberofsteps,1);
k4_10=zeros(numberofsteps,1);
k1_11=zeros(numberofsteps,1);
k2_11=zeros(numberofsteps,1);
k3_11=zeros(numberofsteps,1);
k4_11=zeros(numberofsteps,1);
k1_12=zeros(numberofsteps,1);
k2_12=zeros(numberofsteps,1);
k3_12=zeros(numberofsteps,1);
k4_12=zeros(numberofsteps,1);
k1_13=zeros(numberofsteps,1);
k2_13=zeros(numberofsteps,1);
k3_13=zeros(numberofsteps,1);
k4_13=zeros(numberofsteps,1);
k1_14=zeros(numberofsteps,1);
k2_14=zeros(numberofsteps,1);
k3_14=zeros(numberofsteps,1);
k4_14=zeros(numberofsteps,1);
k1_15=zeros(numberofsteps,1);
k2_15=zeros(numberofsteps,1);
k3_15=zeros(numberofsteps,1);
k4_15=zeros(numberofsteps,1);
k1_16=zeros(numberofsteps,1);
k2_16=zeros(numberofsteps,1);
k3_16=zeros(numberofsteps,1);
k4_16=zeros(numberofsteps,1);
k1_17=zeros(numberofsteps,1);
k2_17=zeros(numberofsteps,1);
k3_17=zeros(numberofsteps,1);
k4_17=zeros(numberofsteps,1);
k1_18=zeros(numberofsteps,1);
k2_18=zeros(numberofsteps,1);
k3_18=zeros(numberofsteps,1);
k4_18=zeros(numberofsteps,1);
k1_19=zeros(numberofsteps,1);
k2_19=zeros(numberofsteps,1);
k3_19=zeros(numberofsteps,1);
k4_19=zeros(numberofsteps,1);
k1_20=zeros(numberofsteps,1);
k2_20=zeros(numberofsteps,1);
k3_20=zeros(numberofsteps,1);
k4_20=zeros(numberofsteps,1);
k1_21=zeros(numberofsteps,1);
k2_21=zeros(numberofsteps,1);
k3_21=zeros(numberofsteps,1);
k4_21=zeros(numberofsteps,1);
k1_22=zeros(numberofsteps,1);
k2_22=zeros(numberofsteps,1);
k3_22=zeros(numberofsteps,1);
k4_22=zeros(numberofsteps,1);
% pathogen load parameters start
kpg=0.008;
kd=10^-4;
um=1.1*10^2;
kpm=2*10^-2;
sm=0.48*2*5;
a=2*10^-3;
kcNA=1000;
b=1.1*10^-3;
kcMA=100;
% pathogen load parameters start
% resting neutrophils and active neutrophils parameters start
knr=2.04;
Ns=2.5*10^6;
kNRd=1.224;
kNAd=0.364;
kNP=0.047;
xNP=100.23;
hNP=2;
kNTNF=10.92;
xNTNF=400;
kN6=10.92;
xN6=393;
kN8=10.64;
xN8=554;
xN10=48;
% resting neutrophils and active neutrophils parameters end
% resting monocytes and active monocytes parameters start
kMR=1;
Minfinty=3*10^4;
kMRd=2/3;
kMP=0.052;
xMP=4.93;
kMTNF=14;
xMTNF=723;
kM6=15.7;
xM6=783.5;
hMP=1;
kMAd=0.52;
xM10=81;
% resting monocytes and active monocytes parameters end
% TNF parameters start
kTNFN=377;
kTNFM=106;
hTNFNM=1.12;
xTNFN=100*6*10^3;
xTNFM=30*3*10^2;
xGTNF=(60/2);
xLTNF=(95*10^-2*2);
xTNF6=250;
xTNF10=235;
kTNF=0.12;
wTNF=4.07;
% TNF parameters end
% IL6 parameters start
k6N=333;
k6M=189;
h6NM=3.2;
x6N=10^5;
x6M=10^3;
k6TNF=53.89;
x6TNF=350;
x66=546;
x610=40;
xL6=(115*10^-2*2);
k6=0.32;
wIL6=3.33;
k6NO=3;
x6N0=0.08;
xG6=(55/2);
% IL6 parameters end
% IL8 parameters start
k8N=17.3*10^2;
x8N=2*10^5;
h8N=1;
k8TNF=17.46;
x8TNF=300;
x810=30;
k8=0.7;
wIL8=25.75;
% IL8 parameters end
% IL10 parameters start
k10N=0.012;
k10M=0.022;
x10N=10^5;
x10M=10^3;
k10TNF=3.8*10^3;
k106=4.3*10^3;
x10TNF=200;
x106=400;
k10=0.293;
wIL10=1.45;
% IL10 parameters end
% Damage parameters start
f=0.07;
T=2510;
kA=2.93;
kdp=15.51;
xdp=142.4;
ud=0.33;
% Damage parameters end
% Nitric Oxide parameters start
nNOTNF=195;
kNONM=1.2*10^-6*1.5;
kNOD=6.3*10^-5*2.5;
nNO10=150;
uno=0.025;
% Nitric Oxide parameters end
% Temperature parameters start
kTTNF=1.2;
kT6=1.6;
kT10=0.4;
kT=0.7;
nTTNF=130;
nT6=140;
nT10=40;
hT=1;
Tb=36.2;
% Temperature parameters end
% Heart rate parameters start
kH=0.46;
HM=190;
BPb=118.8;
Hb=60;
hHP=4;
nHPhBP=14.6;
nHPlBP=22.11;
nHT=0.93;
tauH=0.79;
% Heart rate parameters end
% endotoxin parameters start
kpe=112.3;
npe=100.3;
kde=0.2;
hET=2;
% endotoxin parameters end
% pain perception parameters start
kPTe=2*10^-3;
kPT=0.5;
PTb=781;
% pain perception parameters end
% glucose and glycogen parameters start
sg=109.8;
sIN=0.077;
kLG=38.23;
xING=10.5;
xTNFG=107;
xIL6G=133;
kglucgly=1.74;
kglyglu=1.827;
% glucose and glycogen parameters end
% Lactate parameters start
kGL=0.0046;
%original xNOL start
xNOL=1.5;
%original xNOL end
xETG=100;
uL=0.0913;
Lact0=1.25;
% Lactate parameters end
% Insulin parameters start
alpha=50;
Go=105;
xINTNF=407;
xIN6=333;
kIN=0.0417;
kING=4.9;
wINSULIN=13.5;
% Insulin parameters end
% Respiratory rate parameters start
RRb=20;
kRRlact=100;
nRRlact=1;
nTempRR=0.24;
nPaO2RR=70;
tauRR=0.1;
% Respiratory rate parameters end
% PH parameters start
kPHlact=0.01;
nPaO2Ph=20;
kpH=0.5;
PHbaseline=7.4;
Pao2b=82.5;
% PH parameters end
% PaO2 parameters start
kinfPao2=0.41;
nTempPao2=0.4*1.2;
nPHPao2=7.4;
ko2Pmetab=0.03;
% PaO2 parameters end
P(1)=6*10^4;
NR(1)=10^6;
NA(1)=0;
MR(1)=10^4;
MA(1)=0;
TNF(1)=4.07;
IL6(1)=3.33;
IL8(1)=25.75;
IL10(1)=1.45;
D(1)=0;
NO(1)=0;
Temp(1)=36.2;
HR(1)=60;
ET(1)=0;
PT(1)=781;
Glu(1)=105;
Gly(1)=100;
Lact(1)=1.25;
IN(1)=13.5;
RR(1)=20;
PH(1)=7.4;
Pao2(1)=82.5;
t(1)=0;
tspan=0:h:500;
MAPdata=[115 120 117 112 112.5 113 113 113 113 113 113 113 113];
tsimulation=0:4:48;
if (tspan<=max(tsimulation))
MAP=interp1(tsimulation,MAPdata,tspan,"linear");
else
MAP=interp1(tsimulation,MAPdata,tspan,"linear","extrap");
end
dPdt=@(t,P,NA,MA) kpg*P-kd*P-((sm*P)/(um+kpm*P))-a*HU2(P,kcNA,1)*NA-b*HU2(P,kcMA,1)*MA;
dNRdt=@(t,P,NR,TNF,IL6,IL8,IL10) knr*NR*(1-(NR/Ns))-kNP*HU2(P,xNP,hNP)*NR*(1+kNTNF*HU2(TNF,xNTNF,2))*(1+kN6*HU2(IL6,xN6,2))*(1+kN8*HU2(IL8,xN8,3))*HYX(IL10,xN10,2)-kNRd*NR;
dNAdt=@(t,P,NR,NA,TNF,IL6,IL8,IL10) kNP*HU2(P,xNP,hNP)*NR*(1+kNTNF*HU2(TNF,xNTNF,2))*(1+kN6*HU2(IL6,xN6,2))*(1+kN8*HU2(IL8,xN8,3))*HYX(IL10,xN10,2)-kNAd*NA;
dMRdt=@(t,P,MR,TNF,IL6,IL10) kMR*(1-(MR/Minfinty))*MR-kMP*HU2(P,xMP,hMP)*MR*(1+kMTNF*HU2(TNF,xMTNF,2))*(1+kM6*HU2(IL6,xM6,2))*HYX(IL10,xM10,2)-kMRd*MR;
dMAdt=@(t,P,MR,MA,TNF,IL6,IL10) kMP*HU2(P,xMP,hMP)*MR*(1+kMTNF*HU2(TNF,xMTNF,2))*(1+kM6*HU2(IL6,xM6,2))*HYX(IL10,xM10,2)-kMAd*MA;
dTNFdt=@(t,NA,MA,TNF,IL6,IL10,Glu,Lact) (kTNFN*HU2(NA,xTNFN,hTNFNM)+kTNFM*HU2(MA,xTNFM,hTNFNM))*HU2(Glu,xGTNF,2)*HYX(Lact,xLTNF,1)*HYX(IL6,xTNF6,2)*HYX(IL10,xTNF10,2)-kTNF*(TNF-wTNF);
dIL6dt=@(t,NA,MA,TNF,IL6,IL10,Glu,Lact,NO) (k6N*HU2(NA,x6N,h6NM)+k6M*HU2(MA,x6M,h6NM))*(k6TNF*HU2(TNF,x6TNF,2)+k6NO*HU2(NO,x6N0,2))*HYX(Lact,xL6,1)*HU2(Glu,xG6,2)*HYX(IL6,x66,1)*HYX(IL10,x610,2)-k6*(IL6-wIL6);
dIL8dt=@(t,MA,TNF,IL8,IL10) (k8N*HU2(NA,x8N,h8N)+k8TNF*HU2(TNF-wTNF,x8TNF,3))*HYX(IL10,x810,2)-k8*(IL8-wIL8);
dIL10dt=@(t,NA,MA,TNF,IL6,IL10) (k10N*HU2(NA,x10N,4)+k10M*HU2(MA,x10M,4))*(k10TNF*HU2(TNF,x10TNF,2)+k106*HU2(IL6,x106,4))-k10*(IL10-wIL10);
z=f*(NA+MA)-T;
if (z>0)
u=z;
else
u=0;
end
dDdt=@(t,NA,MA,P,D,IL10) (u/(1+kA*IL10))+kdp*hV(P,xdp)-ud*D;
dNOdt=@(t,NA,MA,TNF,D,NO,IL10) (kNONM*(NA+MA)+kNOD*D)*HU2(TNF-wTNF,nNOTNF,2)*HYX(IL10-wIL10,nNO10,1)-uno*NO;
dTempdt=@(t,Temp,TNF,IL6,IL10) kTTNF*HU2((TNF-wTNF),nTTNF,hT)+kT6*HU2((IL6-wIL6),nT6,hT)-kT10*(1-HYX((IL10-wIL10),nT10,hT))-kT*(Temp-Tb);
Nstep=size(tspan,2);
for j=1:Nstep
if(MAP(j)<=100)
fBP(j)=HU2(BPb-MAP(j),nHPlBP,hHP);
else
fBP(j)=HYX(BPb-MAP(j),nHPhBP,hHP);
end
dHRdt=@(t,HR,Temp) (Hb-HR+(kH*(HM-Hb)*HU2((Temp-Tb),nHT,2)*fBP(j)))/tauH;
end
dETdt=@(t,P,ET) kpe*HU2(P,npe,hET)-kde*ET;
dPTdt=@(t,ET,PT) -kPTe*ET*PT-kPT*(PT-PTb);
dGludt=@(t,Lact,Glu,Gly,TNF,IL6,IN) (sg+kLG*(Lact-Lact0)+(kglyglu*Gly*(1+HU2(TNF-wTNF,xTNFG,2))*(1+HU2(IL6-wIL6,xIL6G,2))))*HYX(IN-wINSULIN,xING,1)-kglucgly*Glu-sIN*Glu*IN;
dGlydt=@(t,Glu,Gly,TNF,IL6,IN) kglucgly*Glu-kglyglu*Gly*(1+HU2(TNF-wTNF,xTNFG,2))*(1+HU2(IL6-wIL6,xIL6G,2))*HYX(IN-wINSULIN,xING,1);
dLactdt=@(t,Glu,NO,Lact,IN,ET) kGL*Glu*(1-HYX(NO,xNOL,2))*(1+HU2(ET,xETG,1))-uL*(Lact-Lact0);
dINdt=@(t,Glu,TNF,IL6,IN) kING*(1+HU2(Glu-Go,alpha,2))*(1-HYX(TNF-wTNF,xINTNF,1))*(1-HYX(IL6-wIL6,xIN6,1))-kIN*(IN-wINSULIN);
dRRdt=@(t,Lact,Temp,RR,Pao2) (RRb-RR+(kRRlact*(1+HU2(Lact,nRRlact,1))*HU2(Temp-Tb,nTempRR,2)*HYX(Pao2,nPaO2RR,1)))/tauRR;
dPHdt=@(t,Lact,PH,Pao2) -kPHlact*PH*Lact*(1-HYX(Pao2b-Pao2,nPHPao2,1))+kpH*(PHbaseline-PH);
dPao2dt=@(t,Lact,Temp,Pao2) -kinfPao2*HU2(Temp-Tb,nTempPao2,2)*HYX(Pao2,nPaO2Ph,1)-ko2Pmetab*(Pao2-Pao2b);
for i=1:numberofsteps
k1_1(i)=dPdt(t(i),P(i),NA(i),MA(i));
k2_1(i)=dPdt(t(i)+.5*h,P(i)+0.5*k1_1(i)*h,NA(i)+0.5*h*k1_1(i),MA(i)+0.5*h*k1_1(i));
k3_1(i)=dPdt(t(i)+.5*h,P(i)+0.5*k2_1(i)*h,NA(i)+0.5*h*k2_1(i),MA(i)+0.5*h*k2_1(i));
k4_1(i)=dPdt(t(i)+h,P(i)+k3_1(i)*h,NA(i)+h*k3_1(i),MA(i)+h*k3_1(i));
P(i+1)=P(i)+(k1_1(i)+2*k2_1(i)+2*k3_1(i)+k4_1(i))*h/6;
if(P(i+1)<1)
P(i+1)=0;
else
P(i+1)=P(i+1);
end
k1_2(i)=dNRdt(t(i),P(i),NR(i),TNF(i),IL6(i),IL8(i),IL10(i));
k2_2(i)=dNRdt(t(i)+.5*h,P(i)+0.5*k1_2(i)*h,NR(i)+0.5*h*k1_2(i),TNF(i)+0.5*h*k1_2(i),IL6(i)+0.5*h*k1_2(i),IL8(i)+0.5*h*k1_2(i),IL10(i)+0.5*h*k1_2(i));
k3_2(i)=dNRdt(t(i)+.5*h,P(i)+0.5*k2_2(i)*h,NR(i)+0.5*h*k2_2(i),TNF(i)+0.5*h*k2_2(i),IL6(i)+0.5*h*k2_2(i),IL8(i)+0.5*h*k2_2(i),IL10(i)+0.5*h*k2_2(i));
k4_2(i)=dNRdt(t(i)+h,P(i)+k3_2(i)*h,NR(i)+h*k3_2(i),TNF(i)+h*k3_2(i),IL6(i)+h*k3_2(i),IL8(i)+h*k3_2(i),IL10(i)+h*k3_2(i));
NR(i+1)=NR(i)+(k1_2(i)+2*k2_2(i)+2*k3_2(i)+k4_2(i))*h/6;
if(NR(i+1)>NR(1))
NR(i+1)=NR;
else
NR(i+1)=NR(i+1);
end
k1_3(i)=dNAdt(t(i),P(i),NR(i),NA(i),TNF(i),IL6(i),IL8(i),IL10(i));
k2_3(i)=dNAdt(t(i)+.5*h,P(i)+0.5*k1_3(i)*h,NR(i)+0.5*h*k1_3(i),NA(i)+0.5*h*k1_3(i),TNF(i)+0.5*h*k1_3(i),IL6(i)+0.5*h*k1_3(i),IL8(i)+0.5*h*k1_3(i),IL10(i)+0.5*h*k1_3(i));
k3_3(i)=dNAdt(t(i)+.5*h,P(i)+0.5*k2_2(i)*h,NR(i)+0.5*h*k2_2(i),NA(i)+0.5*h*k2_2(i),TNF(i)+0.5*h*k2_2(i),IL6(i)+0.5*h*k2_2(i),IL8(i)+0.5*h*k2_2(i),IL10(i)+0.5*h*k2_2(i));
k4_3(i)=dNAdt(t(i)+h,P(i)+k3_3(i)*h,NR(i)+h*k3_3(i),NA(i)+h*k3_3(i),TNF(i)+h*k3_3(i),IL6(i)+h*k3_3(i),IL8(i)+h*k3_3(i),IL10(i)+h*k3_3(i));
NA(i+1)=NA(i)+(k1_3(i)+2*k2_3(i)+2*k3_3(i)+k4_3(i))*h/6;
if(NA(i+1)<1)
NA(i+1)=0;
else
NA(i+1)=NA(i+1);
end
k1_4(i)=dMRdt(t(i),P(i),MR(i),TNF(i),IL6(i),IL10(i));
k2_4(i)=dMRdt(t(i)+.5*h,P(i)+0.5*k1_4(i)*h,MR(i)+0.5*h*k1_4(i),TNF(i)+0.5*h*k1_4(i),IL6(i)+0.5*h*k1_4(i),IL10(i)+0.5*h*k1_4(i));
k3_4(i)=dMRdt(t(i)+.5*h,P(i)+0.5*k2_4(i)*h,MR(i)+0.5*h*k2_4(i),TNF(i)+0.5*h*k2_4(i),IL6(i)+0.5*h*k2_4(i),IL10(i)+0.5*h*k2_4(i));
k4_4(i)=dMRdt(t(i)+h,P(i)+k3_4(i)*h,NR(i)+h*k3_4(i),TNF(i)+h*k3_4(i),IL6(i)+h*k3_4(i),IL10(i)+h*k3_4(i));
MR(i+1)=MR(i)+(k1_4(i)+2*k2_4(i)+2*k3_4(i)+k4_4(i))*h/6;
if(MR(i+1)>MR(1))
MR(i+1)=MR;
else
MR(i+1)=MR(i+1);
end
k1_5(i)=dMAdt(t(i),P(i),MR(i),MA(i),TNF(i),IL6(i),IL10(i));
k2_5(i)=dMAdt(t(i)+.5*h,P(i)+0.5*k1_5(i)*h,MR(i)+0.5*h*k1_5(i),MA(i)+0.5*h*k1_5(i),TNF(i)+0.5*h*k1_5(i),IL6(i)+0.5*h*k1_5(i),IL10(i)+0.5*h*k1_5(i));
k3_5(i)=dMAdt(t(i)+.5*h,P(i)+0.5*k2_5(i)*h,MR(i)+0.5*h*k2_5(i),MA(i)+0.5*h*k2_5(i),TNF(i)+0.5*h*k2_5(i),IL6(i)+0.5*h*k2_5(i),IL10(i)+0.5*h*k2_5(i));
k4_5(i)=dMAdt(t(i)+h,P(i)+k3_5(i)*h,MR(i)+h*k3_5(i),MA(i)+h*k3_5(i),TNF(i)+h*k3_5(i),IL6(i)+h*k3_5(i),IL10(i)+h*k3_5(i));
MA(i+1)=MA(i)+(k1_5(i)+2*k2_5(i)+2*k3_5(i)+k4_5(i))*h/6;
if(MA(i+1)<1)
MA(i+1)=0;
else
MA(i+1)=NA(i+1);
end
k1_6(i)=dTNFdt(t(i),NA(i),MA(i),TNF(i),IL6(i),IL10(i),Glu(i),Lact(i));
k2_6(i)=dTNFdt(t(i)+.5*h,NA(i)+0.5*k1_6(i)*h,MA(i)+0.5*k1_6(i)*h,TNF(i)+0.5*k1_6(i)*h,IL6(i)+0.5*k1_6(i)*h,IL10(i)+0.5*k1_6(i)*h,Glu(i)+0.5*k1_6(i)*h,Lact(i)+0.5*k1_6(i)*h);
k3_6(i)=dTNFdt(t(i)+.5*h,NA(i)+0.5*k2_6(i)*h,MA(i)+0.5*k2_6(i)*h,TNF(i)+0.5*k2_6(i)*h,IL6(i)+0.5*k2_6(i)*h,IL10(i)+0.5*k2_6(i)*h,Glu(i)+0.5*k2_6(i)*h,Lact(i)+0.5*k2_6(i)*h);
k4_6(i)=dTNFdt(t(i)+h,NA(i)+k3_6(i)*h,MA(i)+k3_6(i)*h,TNF(i)+k3_6(i)*h,IL6(i)+k3_6(i)*h,IL10(i)+0.5*k3_6(i)*h,Glu(i)+0.5*k3_6(i)*h,Lact(i)+0.5*k3_6(i)*h);
TNF(i+1)=TNF(i)+(k1_6(i)+2*k2_6(i)+2*k3_6(i)+k4_6(i))*h/6;
k1_7(i)=dIL6dt(t(i),NA(i),MA(i),TNF(i),IL6(i),IL10(i),Glu(i),Lact(i),NO(i));
k2_7(i)=dIL6dt(t(i)+.5*h,NA(i)+0.5*k1_7(i)*h,MA(i)+0.5*k1_7(i)*h,TNF(i)+0.5*k1_7(i)*h,IL6(i)+0.5*k1_7(i)*h,IL10(i)+0.5*k1_7(i)*h,Glu(i)+0.5*k1_7(i)*h,Lact(i)+0.5*k1_7(i)*h,NO(i)+0.5*k1_7(i)*h);
k3_7(i)=dIL6dt(t(i)+.5*h,NA(i)+0.5*k2_7(i)*h,MA(i)+0.5*k2_7(i)*h,TNF(i)+0.5*k2_7(i)*h,IL6(i)+0.5*k2_7(i)*h,IL10(i)+0.5*k2_7(i)*h,Glu(i)+0.5*k2_7(i)*h,Lact(i)+0.5*k2_7(i)*h,NO(i)+0.5*k2_7(i)*h);
k4_7(i)=dIL6dt(t(i)+h,NA(i)+k3_7(i)*h,MA(i)+k3_7(i)*h,TNF(i)+k3_7(i)*h,IL6(i)+k3_7(i)*h,IL10(i)+k3_7(i)*h,Glu(i)+k3_7(i)*h,Lact(i)+k3_7(i)*h,NO(i)+k3_7(i)*h);
IL6(i+1)=IL6(i)+(k1_7(i)+2*k2_7(i)+2*k3_7(i)+k4_7(i))*h/6;
k1_8(i)=dIL8dt(t(i),NA(i),TNF(i),IL10(i),IL8(i));
k2_8(i)=dIL8dt(t(i)+.5*h,NA(i)+0.5*k1_8(i)*h,TNF(i)+0.5*k1_8(i)*h,IL8(i)+0.5*k1_8(i)*h,IL10(i)+0.5*k1_8(i)*h);
k3_8(i)=dIL8dt(t(i)+.5*h,NA(i)+0.5*k2_8(i)*h,TNF(i)+0.5*k2_8(i)*h,IL8(i)+0.5*k2_8(i)*h,IL10(i)+0.5*k2_8(i)*h);
k4_8(i)=dIL8dt(t(i)+h,NA(i)+k3_8(i)*h,TNF(i)+k3_8(i)*h,IL8(i)+k3_8(i)*h,IL10(i)+k3_8(i)*h);
IL8(i+1)=IL8(i)+(k1_8(i)+2*k2_8(i)+2*k3_8(i)+k4_8(i))*h/6;
k1_9(i)=dIL10dt(t(i),NA(i),MA(i),TNF(i),IL6(i),IL10(i));
k2_9(i)=dIL10dt(t(i)+.5*h,NA(i)+0.5*k1_9(i)*h,MA(i)+0.5*k1_9(i)*h,TNF(i)+0.5*k1_9(i)*h,IL6(i)+0.5*k1_9(i)*h,IL10(i)+0.5*k1_9(i)*h);
k3_9(i)=dIL10dt(t(i)+.5*h,NA(i)+0.5*k2_9(i)*h,MA(i)+0.5*k2_9(i)*h,TNF(i)+0.5*k2_9(i)*h,IL6(i)+0.5*k2_9(i)*h,IL10(i)+0.5*k2_9(i)*h);
k4_9(i)=dIL10dt(t(i)+h,NA(i)+k3_9(i)*h,MA(i)+k3_9(i)*h,TNF(i)+k3_9(i)*h,IL6(i)+k3_9(i)*h,IL10(i)+k3_9(i)*h);
IL10(i+1)=IL10(i)+(k1_9(i)+2*k2_9(i)+2*k3_9(i)+k4_9(i))*h/6;
k1_10(i)=dDdt(t(i),NA(i),MA(i),P(i),D(i),IL10(i));
k2_10(i)=dDdt(t(i)+.5*h,NA(i)+0.5*k1_10(i)*h,MA(i)+0.5*k1_10(i)*h,P(i)+0.5*k1_10(i)*h,D(i)+0.5*k1_9(i)*h,IL10(i)+0.5*k1_10(i)*h);
k3_10(i)=dDdt(t(i)+.5*h,NA(i)+0.5*k2_10(i)*h,MA(i)+0.5*k2_10(i)*h,P(i)+0.5*k2_10(i)*h,D(i)+0.5*k2_10(i)*h,IL10(i)+0.5*k2_10(i)*h);
k4_10(i)=dDdt(t(i)+h,NA(i)+k3_10(i)*h,MA(i)+k3_10(i)*h,P(i)+k3_10(i)*h,D(i)+k3_10(i)*h,IL10(i)+k3_10(i)*h);
D(i+1)=D(i)+(k1_10(i)+2*k2_10(i)+2*k3_10(i)+k4_10(i))*h/6;
k1_11(i)=dNOdt(t(i),NA(i),MA(i),TNF(i),D(i),IL10(i),NO(i));
k2_11(i)=dNOdt(t(i)+.5*h,NA(i)+0.5*k1_11(i)*h,MA(i)+0.5*k1_11(i)*h,TNF(i)+0.5*k1_11(i)*h,D(i)+0.5*k1_1(i)*h,IL10(i)+0.5*k1_10(i)*h,NO(i)+0.5*k1_11(i)*h);
k3_11(i)=dNOdt(t(i)+.5*h,NA(i)+0.5*k2_11(i)*h,MA(i)+0.5*k2_11(i)*h,TNF(i)+0.5*k2_11(i)*h,D(i)+0.5*k2_11(i)*h,IL10(i)+0.5*k2_11(i)*h,NO(i)+0.5*k2_11(i)*h);
k4_11(i)=dNOdt(t(i)+h,NA(i)+k3_11(i)*h,MA(i)+k3_11(i)*h,TNF(i)+k3_11(i)*h,D(i)+k3_11(i)*h,IL10(i)+k3_11(i)*h,NO(i)+k3_11(i)*h);
NO(i+1)=NO(i)+(k1_11(i)+2*k2_11(i)+2*k3_11(i)+k4_11(i))*h/6;
k1_12(i)=dTempdt(t(i),TNF(i),IL6(i),IL10(i),Temp(i));
k2_12(i)=dTempdt(t(i)+.5*h,TNF(i)+0.5*k1_12(i)*h,IL6(i)+0.5*k1_12(i)*h,IL10(i)+0.5*k1_12(i)*h,Temp+0.5*k1_12(i)*h);
k3_12(i)=dTempdt(t(i)+.5*h,TNF(i)+0.5*k2_12(i)*h,IL6(i)+0.5*k2_12(i)*h,IL10(i)+0.5*k2_12(i)*h,Temp(i)+0.5*k2_12(i)*h);
k4_12(i)=dTempdt(t(i)+h,TNF(i)+k3_12(i)*h,IL6(i)+k3_12(i)*h,IL10(i)+k3_12(i)*h,Temp(i)+k3_12(i)*h);
Temp(i+1)=Temp(i)+(k1_12(i)+2*k2_12(i)+2*k3_12(i)+k4_12(i))*h/6;
k1_13(i)=dHRdt(t(i),HR(i),Temp(i));
k2_13(i)=dHRdt(t(i)+.5*h,HR(i)+0.5*k1_13(i)*h,Temp+0.5*k1_13(i)*h);
k3_13(i)=dHRdt(t(i)+.5*h,HR(i)+0.5*k2_13(i)*h,Temp+0.5*k2_13(i)*h);
k4_13(i)=dHRdt(t(i)+h,HR(i)+k3_13(i)*h,Temp+k3_13(i)*h);
HR(i+1)=HR(i)+(k1_13(i)+2*k2_13(i)+2*k3_13(i)+k4_13(i))*h/6;
k1_14(i)=dETdt(t(i),P(i),ET(i));
k2_14(i)=dETdt(t(i)+.5*h,P(i)+0.5*k1_14(i)*h,ET+0.5*k1_14(i)*h);
k3_14(i)=dETdt(t(i)+.5*h,P(i)+0.5*k2_14(i)*h,ET+0.5*k2_14(i)*h);
k4_14(i)=dETdt(t(i)+h,P(i)+k3_14(i)*h,ET+k3_14(i)*h);
ET(i+1)=ET(i)+(k1_14(i)+2*k2_14(i)+2*k3_14(i)+k4_14(i))*h/6;
k1_15(i)=dPTdt(t(i),PT(i),ET(i));
k2_15(i)=dPTdt(t(i)+.5*h,PT(i)+0.5*k1_15(i)*h,ET+0.5*k1_15(i)*h);
k3_15(i)=dPTdt(t(i)+.5*h,PT(i)+0.5*k2_15(i)*h,ET+0.5*k2_15(i)*h);
k4_15(i)=dPTdt(t(i)+h,PT(i)+k3_15(i)*h,ET+k3_15(i)*h);
PT(i+1)=PT(i)+(k1_15(i)+2*k2_15(i)+2*k3_15(i)+k4_15(i))*h/6;
k1_16(i)=dGludt(t(i),Glu(i),Gly(i),TNF(i),Lact(i),IL6(i),IN(i));
k2_16(i)=dGludt(t(i)+.5*h,Glu(i)+0.5*k1_16(i)*h,Gly(i)+0.5*k1_16(i)*h,TNF(i)+0.5*k1_16(i)*h,Lact(i)+0.5*k1_16(i)*h,IL6(i)+0.5*k1_16(i)*h,IN(i)+0.5*k1_16(i)*h);
k3_16(i)=dGludt(t(i)+.5*h,Glu(i)+0.5*k2_16(i)*h,Gly(i)+0.5*k2_16(i)*h,TNF(i)+0.5*k2_16(i)*h,Lact(i)+0.5*k2_16(i)*h,IL6(i)+0.5*k2_16(i)*h,IN(i)+0.5*k2_16(i)*h);
k4_16(i)=dGludt(t(i)+h,Glu(i)+k3_16(i)*h,Gly(i)+k3_16(i)*h,TNF(i)+k3_16(i)*h,Lact(i)+k3_16(i)*h,IL6(i)+k3_16(i)*h,IN(i)+k3_16(i)*h);
Glu(i+1)=Glu(i)+(k1_16(i)+2*k2_16(i)+2*k3_16(i)+k4_16(i))*h/6;
k1_17(i)=dGlydt(t(i),Glu(i),Gly(i),TNF(i),IL6(i),IN(i));
k2_17(i)=dGlydt(t(i)+.5*h,Glu(i)+0.5*k1_17(i)*h,Gly(i)+0.5*k1_17(i)*h,TNF(i)+0.5*k1_17(i)*h,IL6(i)+0.5*k1_17(i)*h,IN(i)+0.5*k1_17(i)*h);
k3_17(i)=dGlydt(t(i)+.5*h,Glu(i)+0.5*k2_17(i)*h,Gly(i)+0.5*k2_17(i)*h,TNF(i)+0.5*k2_17(i)*h,IL6(i)+0.5*k2_17(i)*h,IN(i)+0.5*k2_17(i)*h);
k4_17(i)=dGlydt(t(i)+h,Glu(i)+k3_17(i)*h,Gly(i)+k3_17(i)*h,TNF(i)+k3_17(i)*h,IL6(i)+k3_17(i)*h,IN(i)+k3_17(i)*h);
Gly(i+1)=Gly(i)+(k1_17(i)+2*k2_17(i)+2*k3_17(i)+k4_17(i))*h/6;
k1_18(i)=dLactdt(t(i),Glu(i),NO(i),Lact(i),ET(i));
k2_18(i)=dLactdt(t(i)+.5*h,Glu(i)+0.5*k1_18(i)*h,NO(i)+0.5*k1_18(i)*h,Lact(i)+0.5*k1_18(i)*h,ET(i)+0.5*k1_18(i)*h);
k3_18(i)=dLactdt(t(i)+.5*h,Glu(i)+0.5*k2_18(i)*h,NO(i)+0.5*k2_18(i)*h,Lact(i)+0.5*k2_18(i)*h,ET(i)+0.5*k2_18(i)*h);
k4_18(i)=dLactdt(t(i)+h,Glu(i)+k3_18(i)*h,NO(i)+k3_18(i)*h,Lact(i)+k3_18(i)*h,ET(i)+k3_18(i)*h);
Lact(i+1)=Lact(i)+(k1_18(i)+2*k2_18(i)+2*k3_18(i)+k4_18(i))*h/6;
k1_19(i)=dINdt(t(i),Glu(i),TNF(i),IL6(i),IN(i));
k2_19(i)=dINdt(t(i)+.5*h,Glu(i)+0.5*k1_19(i)*h,TNF(i)+0.5*k1_19(i)*h,IL6(i)+0.5*k1_19(i)*h,IN(i)+0.5*k1_19(i)*h);
k3_19(i)=dINdt(t(i)+.5*h,Glu(i)+0.5*k2_19(i)*h,TNF(i)+0.5*k2_19(i)*h,IL6(i)+0.5*k2_19(i)*h,IN(i)+0.5*k2_19(i)*h);
k4_19(i)=dINdt(t(i)+h,Glu(i)+k3_19(i)*h,TNF(i)+k3_19(i)*h,IL6(i)+k3_19(i)*h,IN(i)+k3_19(i)*h);
IN(i+1)=IN(i)+(k1_19(i)+2*k2_19(i)+2*k3_19(i)+k4_19(i))*h/6;
k1_20(i)=dRRdt(t(i),Lact(i),Temp(i),RR(i),Pao2(i));
k2_20(i)=dRRdt(t(i)+.5*h,Lact(i)+0.5*k1_20(i)*h,Temp(i)+0.5*k1_20(i)*h,RR(i)+0.5*k1_20(i)*h,Pao2(i)+0.5*k1_20(i)*h);
k3_20(i)=dRRdt(t(i)+.5*h,Lact(i)+0.5*k2_20(i)*h,Temp(i)+0.5*k2_20(i)*h,RR(i)+0.5*k2_20(i)*h,Pao2(i)+0.5*k2_20(i)*h);
k4_20(i)=dRRdt(t(i)+h,Lact(i)+k3_20(i)*h,Temp(i)+k3_20(i)*h,RR(i)+k3_20(i)*h,Pao2(i)+k3_20(i)*h);
RR(i+1)=RR(i)+(k1_20(i)+2*k2_20(i)+2*k3_20(i)+k4_20(i))*h/6;
k1_21(i)=dPHdt(t(i),Lact(i),PH(i),Pao2(i));
k2_21(i)=dPHdt(t(i)+.5*h,Lact(i)+0.5*k1_21(i)*h,PH(i)+0.5*k1_21(i)*h,Pao2(i)+0.5*k1_21(i)*h);
k3_21(i)=dPHdt(t(i)+.5*h,Lact(i)+0.5*k2_21(i)*h,PH(i)+0.5*k2_21(i)*h,Pao2(i)+0.5*k2_21(i)*h);
k4_21(i)=dPHdt(t(i)+h,Lact(i)+k3_21(i)*h,PH(i)+k3_21(i)*h,Pao2(i)+k3_21(i)*h);
PH(i+1)=PH(i)+(k1_21(i)+2*k2_21(i)+2*k3_21(i)+k4_21(i))*h/6;
k1_22(i)=dPao2dt(t(i),Lact(i),Temp(i),Pao2(i));
k2_22(i)=dPao2dt(t(i)+.5*h,Lact(i)+0.5*k1_22(i)*h,Temp(i)+0.5*k1_22(i)*h,Pao2(i)+0.5*k1_22(i)*h);
k3_22(i)=dPao2dt(t(i)+.5*h,Lact(i)+0.5*k2_22(i)*h,Temp(i)+0.5*k2_22(i)*h,Pao2(i)+0.5*k2_22(i)*h);
k4_21(i)=dPao2dt(t(i)+h,Lact(i)+k3_22(i)*h,Temp(i)+k3_22(i)*h,Pao2(i)+k3_22(i)*h);
Pao2(i+1)=Pao2(i)+(k1_22(i)+2*k2_22(i)+2*k3_22(i)+k4_22(i))*h/6;
t(i+1)=t(i)+h;
Y=[P;NR;NA;MR;MA;TNF;IL6;IL8;I10;D;NO;Temp;HR;ET;PT;Glu;Gly;Lact;IN;RR;PH;Pao2];
k=size(Y,2);
end
Unable to perform assignment because the left and right sides have a different number of elements.
for nm=1:k
figure(nm)
plot(t,Y(:,nm))
end
function qw=HU2(X,n,h)
qw=(X.^h)/(X.^h+n.^h);
end
function qz=HYX(X,n,h)
qz=(n.^h)./(X.^h+n.^h);
end
function h=hV(Z,x)
h=(Z.^6)/(Z.^6+x.^6);
end

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (1)

Torsten
Torsten am 12 Dez. 2024 um 12:05
Bearbeitet: Torsten am 12 Dez. 2024 um 12:16
You always have to use the variables at the ith time step, not the complete array in the assignments.
I didn't check your complete code, but below are some instances (marked in bold).
Further you have to start your loop over the time steps later - I corrected that in your code (see above).
If I were you, I'd use one of the MATLAB adaptive ODE integrators (e.g. ode45, ode15s) instead of a self-written fixed-step classical Runge-Kutta code.
k1_12(i)=dTempdt(t(i),TNF(i),IL6(i),IL10(i),Temp(i));
k2_12(i)=dTempdt(t(i)+.5*h,TNF(i)+0.5*k1_12(i)*h,IL6(i)+0.5*k1_12(i)*h,IL10(i)+0.5*k1_12(i)*h,Temp+0.5*k1_12(i)*h);
k3_12(i)=dTempdt(t(i)+.5*h,TNF(i)+0.5*k2_12(i)*h,IL6(i)+0.5*k2_12(i)*h,IL10(i)+0.5*k2_12(i)*h,Temp(i)+0.5*k2_12(i)*h);
k4_12(i)=dTempdt(t(i)+h,TNF(i)+k3_12(i)*h,IL6(i)+k3_12(i)*h,IL10(i)+k3_12(i)*h,Temp(i)+k3_12(i)*h);
Temp(i+1)=Temp(i)+(k1_12(i)+2*k2_12(i)+2*k3_12(i)+k4_12(i))*h/6;
k1_13(i)=dHRdt(t(i),HR(i),Temp(i));
k2_13(i)=dHRdt(t(i)+.5*h,HR(i)+0.5*k1_13(i)*h,Temp+0.5*k1_13(i)*h);
k3_13(i)=dHRdt(t(i)+.5*h,HR(i)+0.5*k2_13(i)*h,Temp+0.5*k2_13(i)*h);
k4_13(i)=dHRdt(t(i)+h,HR(i)+k3_13(i)*h,Temp+k3_13(i)*h);
HR(i+1)=HR(i)+(k1_13(i)+2*k2_13(i)+2*k3_13(i)+k4_13(i))*h/6;
k1_14(i)=dETdt(t(i),P(i),ET(i));
k2_14(i)=dETdt(t(i)+.5*h,P(i)+0.5*k1_14(i)*h,ET+0.5*k1_14(i)*h);
k3_14(i)=dETdt(t(i)+.5*h,P(i)+0.5*k2_14(i)*h,ET+0.5*k2_14(i)*h);
k4_14(i)=dETdt(t(i)+h,P(i)+k3_14(i)*h,ET+k3_14(i)*h);
ET(i+1)=ET(i)+(k1_14(i)+2*k2_14(i)+2*k3_14(i)+k4_14(i))*h/6;
k1_15(i)=dPTdt(t(i),PT(i),ET(i));
k2_15(i)=dPTdt(t(i)+.5*h,PT(i)+0.5*k1_15(i)*h,ET+0.5*k1_15(i)*h);
k3_15(i)=dPTdt(t(i)+.5*h,PT(i)+0.5*k2_15(i)*h,ET+0.5*k2_15(i)*h);
k4_15(i)=dPTdt(t(i)+h,PT(i)+k3_15(i)*h,ET+k3_15(i)*h);
PT(i+1)=PT(i)+(k1_15(i)+2*k2_15(i)+2*k3_15(i)+k4_15(i))*h/6;

Kategorien

Mehr zu Runge Kutta Methods finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by