How to draw Fig. 1 from the attached pdf with this code

3 次查看（过去 30 天）

MINATI 2019-4-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code

编辑： MINATI 2019-4-30

724547.pdf

function main
Pr=1; G=0.1;
% phi=input('phi=');   %%0,.05, .1, .15, .2
phi=0.0;
rhof=997.1;Cpf=4179;kf=0.613;           %for    WATER
rhos=6320;Cps=531.8;ks=76.5;               %for    CuO
a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks));  %%%%Knf
xa=0;xb=6;
solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
sol=bvp4c(@ode,@bc,solinit);
xint=linspace(xa,xb,101);
sxint=deval(sol,xint);
figure(1)
plot(xint,(1-phi)^-2.5*sxint(3,:),'-','Linewidth',1.5);   %for f''(0)/(1-phi)^2.5  vs  phi
xlabel('\eta'); 
ylabel('f''(0)/(1-phi)^2.5');
hold on
function res=bc(ya,yb)
res=[ya(1); ya(2)-1-G*ya(3); ya(4)-1; yb(2); yb(4)];  
end
 function dydx=ode(x,y)    
dydx=[y(2);  y(3);  a1*(y(2)^2-y(3)*y(1));  y(5);  -A*Pr*a2*y(1)*y(5)];              
end
end

http://doi.org/10.1155/2013/724547

[EDITED, Jan, Attachment added].

11 个评论
显示 9更早的评论隐藏 9更早的评论

David Wilson 2019-4-30

If I understand correctly, the BVP you are trying to solve has BCs at infinity. You have chosen 6 (& the paper uses 8), so you might like to validate that approximation.

What exactly is the problem ?

My solution is for the Pr=6.2, \gamma=0.1. This seems to follow your Fig 1.

Jan 2019-4-30

@MINATI: See https://en.wikipedia.org/wiki/Stiff_equation

请先登录，再进行评论。

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

采纳的回答

David Wilson 2019-4-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/459075-how-to-draw-fig-1-from-the-attached-pdf-with-this-code#answer_372949

I didn't bother draw the other 3 lines, but you just ned to make the necessary changes to gamma for that.

If you run something like what you had originally, you only want the fist point of f''().

Pr=6.2; G=0.1;
% phi=input('phi=');   %%0,.05, .1, .15, .2
phi=0.0;
rhof=997.1;Cpf=4179;kf=0.613;           %for    WATER
rhos=6320;Cps=531.8;ks=76.5;               %for    CuO
a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks));  %%%%Knf
BCres= @(ya,yb) ... 
       [ya(1); ya(2)-1-G*ya(3); ya(4)-1; yb(2); yb(4)];  
fODE = @(x,y) ... 
       [y(2);  y(3);  a1*(y(2)^2-y(3)*y(1));  y(5);  -A*Pr*a2*y(1)*y(5)];  
xa=0;xb=8;
solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
sol=bvp4c(fODE,BCres,solinit);
xint=linspace(xa,xb,101);
sxint=deval(sol,xint);
figure(1)
plot(xint,(1-phi)^-2.5*sxint(3,:),'-','Linewidth',1.5);   %for f''(0)/(1-phi)^2.5  vs  phi
xlabel('\eta'); 
ylabel('f''(0)/(1-phi)^2.5');

Now you have to re-run the above, but change phi over the range given in the Fig.

xa=0;xb=8;
phiv = [0:0.04:0.2]'; 
p = []; % collect points here 
for i=1:length(phiv)
    phi = phiv(i); 
    a1=((1-phi)^2.5)*(1-phi+phi*(rhos/rhof));
    a2=(1-phi+phi*((rhos*Cps)/(rhof*Cpf)));
    A=(ks+2*kf+phi*(kf-ks))/(ks+2*kf-2*phi*(kf-ks));  %%%%Knf
    
    fODE = @(x,y) ... 
       [y(2);  y(3);  a1*(y(2)^2-y(3)*y(1));  y(5);  -A*Pr*a2*y(1)*y(5)];  
    solinit=bvpinit(linspace(xa,xb,101),[0 1 0 1 0]);
    sol=bvp4c(fODE,BCres,solinit);
    
    p(i,1) = (1-phi)^-2.5*sxint(3,1)
end 
plot(phiv, p,'o-')
xlabel('\phi'); ylabel('f''''(0) & stuff')