Related to bvp4c matlab code
显示 更早的评论
Dear experts,
Anyone help me,
how to plot isotherm lines, contour plot, surface plot using matlab bvp4c.
suggests any sample codes or any ideas.
Thank you.
this is our code
function hybrid(solver)
if nargin < 1
solver = 'bvp4c';
end
bvpsolver = fcnchk(solver);
xlow=0;
g1=-1;
g2=1;
g3=0.5;
g4=0.5;
global St rad m
St=0.2;
infinity=01;
maxinfinity=15;
solinit=bvpinit(linspace(0,infinity,1001),[St g1 1 g2 0 1 g3 1 g4]);
sol = bvpsolver(@ode,@bc,solinit);
x = sol.x;
y = sol.y;
fprintf('\n');
fprintf('At eta= %g \t %10.8f \t %10.8f \t %10.8f %10.8f.\n',xlow,y(2,1),y(4,1),y(7,1),y(9,1))
for Bnew = infinity:01:maxinfinity
solinit = bvpxtend(sol,Bnew);
sol = bvpsolver(@ode,@bc,solinit);
x = sol.x;
y = sol.y;
fprintf('At eta= %g \t %10.8f \t %10.8f \t %10.8f %10.8f.\n',Bnew,y(2,1),y(4,1),y(7,1),y(9,1))
end
fprintf('%10.8f \t %10.8f \t %10.8f \t %10.8f.\n',y(2,1),y(4,1),y(7,1),y(9,1))
plot(x,y(6,:),'-','Linewidth',1.5)
% hold on
%Cf=rad*(1+rad)^m*((0.32906655)^2+(-1.12865255)^2)^(1/2)
%Nu=-rad*(1+rad)^m*y(7);
%Sh=-rad*(1+rad)^m*y(9)
%Sn=-rad*(1+rad)^m*y(11)
function dydx=ode(x,y)
% rhos1 = 2810; Cps1 = 960; ks1 = 173; Pi1 = 0.1; ss1 = 26.77*(10^6); %AA7075
% rhos2 = 2720; Cps2 = 893.2; ks2 = 222; Pi2 = 0.3; ss2 = 34.83*(10^6); %AA7072
% rhof = 792; Cpf = 2545; kf = 0.2035; sf = 0.5*(10^(-6)); %(Methonal)
Np=2;
alpha=0.15;
eps=0.45;
Re=1;
Pr=1.5;
rad=0.2;
m=1.1;
n=0.1;
Sc=01;
Beta=1;
B=1;
Ha=0.2;
Ec=0.1;
Del=0.01;
h=0.5;
Rp=1.5;
% knf = kf*(ks1+(Np-1)*kf-(Np-1)*Pi1*(kf-ks1))/(ks1+(Np-1)*kf+Pi1*(kf-ks1));
% K =(ks2+(Np-1)*knf+Pi2*(knf-ks2)) /(ks2+(Np-1)*knf-(Np-1)*Pi2*(knf-ks2));%khnf/kf
%
% snf = sf*(ss1+2*sf-2*Pi1*(sf-ss1))/(ss1+2*sf+Pi1*(sf-ss1));
% S = (ss2+2*snf-2*Pi2*(snf-ss2))/(ss2+2*snf+Pi2*(snf-ss2));%shnf/sf
%
% Rh=(1-Pi2)*((1-Pi1)+Pi1*((rhos1*Cps1)/(rhof*Cpf)))+Pi2*((rhos2*Cps2)/(rhof*Cpf));
% Mu=1/((1-Pi1)^(2.5)*(1-Pi2)^(2.5));
% D=1.8;
% RaggRp=3.34;
% Pmax=0.605;
% Pint=(RaggRp)^(D-3);
% Pagg=Pi/Pint;
% rhoagg=(1-Pint)*rhof+Pint*rhos;
% Mu=(1/(Pagg/Pmax))^(-2.5*Pmax);
% Rh=(1-Pagg)+((Pagg*rhoagg)/rhof);
% rhoCpagg=(1-Pint)*(rhof*Cpf)+Pint*rhos*Cps;
% kagg=(kf/4)*((3*Pint-1)*(ks/kf)+(3*(1-Pint)-1)+(((3*Pint-1)*(ks/kf)+(3*(1-Pint)-1))^2+(8*ks/kf))^0.5);
% RhCp=(1-Pagg)+Pagg*(rhoCpagg/(rhof*Cpf));
% K=(kagg+2*kf+2*Pagg*(kagg-kf))/(kagg+2*kf-Pagg*(kagg-kf));
% S=1+(((3*(ss/sf)-1)*Pagg)/((ss/sf)+2)*((ss/sf)-1)*Pagg);
rhos = 2810; Cps = 960; ks = 173; Pi = 0.1; ss = 26.77*(10^6); %AA7075
rhof = 792; Cpf = 2545; kf = 0.2035; sf = 0.5*(10^(-6)); %(Methonal)
Mu=1+2.5*Del+4.5*(1/((h/Rp)*(2+(h/Rp))*(1+(h/Rp))^2));
Rh=(1-Del)+Del*(rhos/rhof);
RhCp=(1-Del)+Del*(rhos*Cps/rhof*Cpf);
K=(ks+2*kf-2*Del*(kf-ks))/(ks+2*kf+Del*(kf-ks));
S=(ss+2*sf-2*Del*(sf-ss))/(ss+2*sf+Del*(sf-ss));
dydx=[y(2)
Rh*(y(1)^2-y(3)^2+m*eps*(x+alpha)*y(1)*y(2)+y(5)*y(2))-(B/Re^2)+S*Ha*y(1)/Mu*Re^((1-n)/(1+n))*(1+rad)^(2*m)*(1+(1/Beta))
y(4)
Rh*(2*y(1)*y(3)+y(5)*y(4)+m*eps*(x+alpha)*y(1)*y(4))+S*Ha*y(3)/Mu*Re^((1-n)/(1+n))*(1+rad)^(2*m)*(1+(1/Beta))
-2*y(1)-m*eps*(x+alpha)*y(2)
y(7)
(y(1)*y(7)*m*eps*(x+alpha)*Pr+Pr*y(5)*y(7)+Pr*(S*RhCp)*Ha*Ec*(y(1)^2+y(3)^2))/(K*RhCp)*Re^((1-n)/(1+n))*(1+rad)^(2*m)
y(9)
(Sc*y(1)*y(9)*m*eps*(x+alpha)+Sc*y(5)*y(9))/(Mu*Re^((1-n)/(1+n))*(1+rad)^(2*m))
];
end
function res=bc(y0,yinf)
res=[y0(1)-St
y0(3)-1
y0(5)
y0(6)-1
y0(8)-1
yinf(1)
yinf(3)
yinf(6)
yinf(8)
];
end
end
回答（1 个）
Syed Sohaib Zafar
2023-10-5
0 个投票
I think You need these kind of graphs.
To plot such graphs, use the command "contourf".
Find more help from the link below.
https://www.mathworks.com/help/matlab/ref/contourf.html
类别
在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
2023-10-5
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
