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

조회 수: 3(최근 30일)
MINATI 29 Apr 2019
편집: MINATI 30 Apr 2019
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
Jan 30 Apr 2019

댓글을 달려면 로그인하십시오.

### 채택된 답변

David Wilson 30 Apr 2019
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')
Resultant plot is as above.
##### 댓글 수: 1표시숨기기 없음
MINATI 30 Apr 2019
Many many thanks David
It worked
Where to put the loop for Gamma G=[0 0.1 0.2 0.3] which will vary the fig

댓글을 달려면 로그인하십시오.

### Community Treasure Hunt

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

Start Hunting!