I am getting an error of 'Subscript indices must either be real positive integers or logicals' while solving a runge kutta problem

3 ビュー (過去 30 日間)
clf;clc;
%parameters
s=0.01;
t_final=10;
global Re Pr R phi b eps Mn Bi alp rhof cpf kf rhos cps ks
Re= 2.5;
Pr= 6.8;
R= 1;
phi= 0.1;
b=1.0;
eps=0.1;
Mn= 1.0;
Bi= 1.0;
alp= 1;
rhof= 997.1;
cpf= 4179.0;
kf= 0.613;
rhos= 8933.0;
cps= 385.0;
ks= 400;
phi1= (1-phi)^2.5*((1-phi)+phi*(rhos/rhof));
phi2= (1-phi)+phi((rhos*cps)/(rhof*cpf));
k = (ks+2*kf-2*phi*(kf-ks))/(ks+2*kf+phi*(kf-ks));
%initial conditions
t(1)=0;
h(1)=1;
%function handle
f=@(t,h) -(1+b)*h+eps*((1+b^2)*phi1*Re+(1+b)*((1-phi)^2.5)*Mn^2)*h^3/3+...
2*eps*alp*(1-phi)^2.5*(k/(k+Bi*h))*h^2-(2/3)*eps^2*(1-phi)^2.5*alp*Pr*(Re*phi2/k+R)*...
((2*k*Bi^2*h^3*h_t)/(k+Bi*h)^3+(2*(1+b)*Bi*k*h^3)/(k+Bi*h)^3+((1+b)*(3*k+Bi*h)*k*h^2)/(k+Bi*h)^2)*h^2/2
% RK4 loop
for i=1:ceil(t_final/s)
t(i+1)=t(i)+s;
k1= f(t(i) , h(i));
k2= f(t(i)+0.5*s , h(i)+0.5*k1*s);
k3= f(t(i)+0.5*s , h(i)+0.5*k2*s);
k4= f(t(i)+s , h(i)+k3*s);
h(i+1)= h(i)+s/6*(k1+2*k2+2*k3+k4);
end

採用された回答

Birdman
Birdman 2018 年 3 月 20 日
Change this line
phi2= (1-phi)+phi((rhos*cps)/(rhof*cpf));
to
phi2= (1-phi)+phi*((rhos*cps)/(rhof*cpf));
  2 件のコメント
naygarp
naygarp 2018 年 3 月 20 日
Thanks, how silly of me.
naygarp
naygarp 2018 年 3 月 20 日
How would I add another for loop if I want to vary the value of 'phi' in the above code and obtain 3 graphs for each value of 'phi'. It would be a nested for loop. say like
qq=[0 0.05 0.1];
for phi=qq

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeNumerical Integration and Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by