I am unable to obtain multiple graphs by varying a parameter in a for loop while solving a Runge-Kutta based problem

Question

naygarp 2018-3-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/389612-i-am-unable-to-obtain-multiple-graphs-by-varying-a-parameter-in-a-for-loop-while-solving-a-runge-kut

编辑： Torsten 2018-3-23

While solving a Runge-Kutta problem, I wanted to obtain 3 graphs varying the values of 'phi' in a 'for' loop, but I am unable to obtain it. Please, can someone correct it?

 clear;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 h_t
Re= 2.5;
Pr= 6.8;
R= 1;
pp=[0 0.05 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;
h_t=2;
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;
 h0(1)=1;
 %function handle
 f=@(t,h0) -(1+b)*h0+eps*((1+b^2)*phi1*Re+(1+b)*((1-phi)^2.5)*Mn^2)*h0^3/3+...
 2*eps*alp*(1-phi)^2.5*(k/(k+Bi*h0))*h0^2-(2/3)*eps^2*(1-phi)^2.5*alp*Pr*(Re*phi2/k+R)*...
     ((2*k*Bi^2*h0^3*h_t)/(k+Bi*h0)^3+(2*(1+b)*Bi*k*h0^3)/(k+Bi*h0)^3+((1+b)*(3*k+Bi*h0)*k*h0^2)/(k+Bi*h0)^2)*h0^2/2;
 % RK4 loop
 for i=1:numel(pp);
     phi=pp(i);
  for i=1:ceil(t_final/s);
     t(i+1)=t(i)+s;
     k1= f(t(i)       ,  h0(i));
     k2= f(t(i)+0.5*s ,  h0(i)+0.5*k1*s);
     k3= f(t(i)+0.5*s ,  h0(i)+0.5*k2*s);
     k4= f(t(i)+s     ,    h0(i)+k3*s);
     h0(i+1)= h0(i)+s/6*(k1+2*k2+2*k3+k4);
 fprintf('\n%16.3f%16.3f',t,h0);
 plot(t,h0);hold on
  end
  end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2018-3-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/389612-i-am-unable-to-obtain-multiple-graphs-by-varying-a-parameter-in-a-for-loop-while-solving-a-runge-kut#answer_311042

在 MATLAB Online 中打开

   clear;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 h_t
  Re= 2.5;
  Pr= 6.8;
  R= 1;
  pp=[0 0.05 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;
  h_t=2;
  for j=1:numel(pp)
    phi=pp(j);
    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(j,1)=1;
    %function handle
    f=@(t,y) -(1+b)*y+eps*((1+b^2)*phi1*Re+(1+b)*((1-phi)^2.5)*Mn^2)*y^3/3+...
             2*eps*alp*(1-phi)^2.5*(k/(k+Bi*y))*y^2-(2/3)*eps^2*(1-phi)^2.5*alp*Pr*(Re*phi2/k+R)*...
             ((2*k*Bi^2*y^3*h_t)/(k+Bi*y)^3+(2*(1+b)*Bi*k*y^3)/(k+Bi*y)^3+((1+b)*(3*k+Bi*y)*k*y^2)/(k+Bi*y)^2)*y^2/2;
   % RK4 loop
    for i=1:ceil(t_final/s);
       t(i+1)=t(i)+s;
       k1= f(t(i)       ,  h(j,i));
       k2= f(t(i)+0.5*s ,  h(j,i)+0.5*k1*s);
       k3= f(t(i)+0.5*s ,  h(j,i)+0.5*k2*s);
       k4= f(t(i)+s     ,  h(j,i)+k3*s);
       h(j,i+1)= h(j,i)+s/6*(k1+2*k2+2*k3+k4);
    end
  end
  plot(t,h(1,:),t,h(2,:),t,h(3,:))

And you should preallocate t and h before entering the first for-loop.

Best wishes

Torsten.

5 个评论
显示 3更早的评论隐藏 3更早的评论

naygarp 2018-3-22

在 MATLAB Online 中打开

will I have to add the function handle 'f=@(t,h)' after in the 2nd for loop

    for i=1:ceil(t_final/s);
       t(i+1)=t(i)+s;
       k1= f(t(i)       ,  h(j,i));
       k2= f(t(i)+0.5*s ,  h(j,i)+0.5*k1*s);
       k3= f(t(i)+0.5*s ,  h(j,i)+0.5*k2*s);
       k4= f(t(i)+s     ,  h(j,i)+k3*s);
       h(j,i+1)= h(j,i)+s/6*(k1+2*k2+2*k3+k4);
       h_t=-(1-b)*h(j,i+1);
       f=@(t,y) -(1+b)*y+eps*((1+b^2)*phi1*Re+(1+b)*((1-phi)^2.5)*Mn^2)*y^3/3+...
             2*eps*alp*(1-phi)^2.5*(k/(k+Bi*y))*y^2-(2/3)*eps^2*(1-phi)^2.5*alp*Pr*(Re*phi2/k+R)*...
             ((2*k*Bi^2*y^3*h_t)/(k+Bi*y)^3+(2*(1+b)*Bi*k*y^3)/(k+Bi*y)^3+((1+b)*(3*k+Bi*y)*k*y^2)/(k+Bi*y)^2)*y^2/2;
    end