Please help me to run the surface figure from the attached code

6 vues (au cours des 30 derniers jours)
Tarek
Tarek le 23 Juil 2025
Réponse apportée : Star Strider le 23 Juil 2025
%Array indices must be positive integers or logical values.
%Error in Tarek7_2025 (line 66)
% N=(kt/kf)*(1+Rd)*d*R(1/(muf/rhof)^2)*(sol.y(7,:).^2)+((mut/kf)*(1/Ti)*(1/at^2))*R(d*(muf/rhof)*(sol.y(2,:).^2)+e(sol.y(4,:).^2)+R*d*p*(sol.y(1,:).^2));
%code
function sol= proj
clc;clf;clear;
%Relation of base fluid
rhof=997.1*10^-3;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;
alfaf=kf/(rhof*cpf);
bef=21*10^-5;
ky=muf/rhof;
disp('ky');disp((muf/rhof));
%sigf=0.05*10^-8;
%Ag
ph1=0.01;
rho1=10500*10^-3;
cp1=235*10^4;
k1=429*10^5;be1=21*10^-5;
%sig1=0.74*10^-2;
%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=400*10^5;
%sig2=5.96*10^-1;
be2=1.67*10^-5;
%Alumina
ph3=0.01;
rho3=3970*10^-3;
cp3=765*10^4;
k3=40*10^5;
be3=0.85*10^-5;
%sig3=3.5*10^-1;
%Relation of ternary hyprid
kn=kf*((k3+2*kf-2*ph3*(kf-k3))/(k3+2*kf+ph3*(kf-k3)));
kh=kn*((k2+2*kn-2*ph2*(kn-k2))/(k2+2*kn+ph2*(kn-k2)));
kt=kh*((k1+2*kh-2*ph1*(kh-k1))/(k1+2*kh+ph1*(kh-k1)));
mut= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5);
rhot=(1-ph1)*((1-ph2)*((1-ph3)+ph3*(rho3/rhof))+ph2*(rho2/rhof))+ph1*(rho1/rhof);
%vt=rhot*cpt
vt =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*cp3)/(rhof*cpf)))+ph2*((rho2*cp2)/(rhof*cpf)))+ph1*((rho1*cp1)/(rhof*cpf));
%disp('vt');disp(vt);
%vb=rho*betb
vb =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*be3)/(rhof*bef)))+ph2*((rho2*be2)/(rhof*bef)))+ph1*((rho1*be1)/(rhof*bef));
%disp('vb');disp(vb);disp(ky);
myLegend1 = {};myLegend2 = {};
rr = [1 2 3 4]
numR = numel(rr);
m = linspace(0,1);
a=0.09;b=0.001;p=.001/((1-0.01)*(mut/muf)*(rhof/rhot));
Ec=10;Gr=0.55;at=0.09;gamma=pi/4;
Rd=1;
prf=6.9;d=0.0001;e=0.001;Rd=0.8;
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1];options =bvpset('stats','on','RelTol',1e-5);
Z = zeros(numR, length(m));
for i = 1:numR
R= rr(i);
if i==1
solinit = bvpinit(m, y0);
else
guess = @(x)interp1(sol.x,(sol.y).',x);
solinit = bvpinit(sol.x,guess);
end
sol = bvp4c(@projfun,@projbc, solinit, options);
N=(kt/kf)*(1+Rd)*d*R(1/(muf/rhof)^2)*(sol.y(7,:).^2)+((mut/kf)*(1/Ti)*(1/at^2))*R(d*(muf/rhof)*(sol.y(2,:).^2)+e(sol.y(4,:).^2)+R*d*p*(sol.y(1,:).^2));
Z(i,:) = interp1(sol.x,N,m);
end
[X, Y] = meshgrid(m, rr);
shading interp;
surf(X, Y, Z);
xlabel('eta');
ylabel('Gr');
zlabel('NG');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;
function dy= projfun(~,y)
dy= zeros(8,1);
% alignComments
E = y(1);
dE = y(2);
F = y(3);
dF= y(4);
w = y(5);
dw=y(6);
t = y(7);
dt = y(8);
dy(1) = dE;
dy(2) = (((rhot/mut)*(a*(muf/rhof)^0.5*(E*F+E^2)+a*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+Gr*a*(muf/rhof)*sin(gamma)*(vb/(rhof*bef))*t)));
dy(3) = dF;
dy(4) = (((rhot/mut)*(b*(muf/rhof)^0.5*(F^2+F*E)+(muf/rhof)*b^0.5*a^(1.5)*dF)));
dy(5) =-(a*F+b*E);
dy(6) = (((rhot/mut)*((muf/rhof)^0.5*w*dw+Gr*b*(muf/rhof)*cos(gamma)*(vb/(rhof*bef))*t)));
dy(7) = dt;
dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*((vt/(rhof*cpf))*(muf/rhof)^0.5*w*dt-(mut/muf)*Ec*1*dw^2) ;
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5);
ya(6);
ya(7)+1-(1/0.9)*ya(8);
yb(1)-0.01;
yb(3);
yb(7);
% yb(7);
];
end

Connectez-vous pour répondre à cette question.

Réponse acceptée

Star Strider
Star Strider le 23 Juil 2025
Ran in:
The error is because you are missing a multiplication operator
dy(2) = (((rhot/mut)*(a*(muf/rhof)^0.5*(E*F+E^2)+a*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+R*a*(muf/rhof)*sin(gamma)*(vb/(rhof*bef))*t)));
HERE ^
and
dy(6) = (((rhot/mut)*((muf/rhof)^0.5*w*dw+R*b*(muf/rhof)*cos(gamma)*(vb/(rhof*bef))*t)));
HERE ^
and
N=(kt/kf)*(1+Rd)*d*R*(1/(muf/rhof)^2)*(sol.y(7,:).^2)+((mut/kf)*(1/Ti)*(1/at^2))*R*(d*(muf/rhof)*(sol.y(2,:).^2)+e*(sol.y(4,:).^2)+R*d*p*(sol.y(1,:).^2));
HERE ^ ^ ^
I inserted them (there are 3 in 'N').
The result otherwise appears to be appropriate.
Running your code --
proj
ky 0.0100
rr = 1×4
1 2 3 4
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
coe 0.8935 The solution was obtained on a mesh of 285 points. The maximum residual is 2.085e-06. There were 15848 calls to the ODE function. There were 185 calls to the BC function. The solution was obtained on a mesh of 308 points. The maximum residual is 1.901e-06. There were 15138 calls to the ODE function. There were 105 calls to the BC function. The solution was obtained on a mesh of 352 points. The maximum residual is 4.438e-06. There were 13545 calls to the ODE function. There were 80 calls to the BC function. The solution was obtained on a mesh of 352 points. The maximum residual is 7.622e-06. There were 9852 calls to the ODE function. There were 54 calls to the BC function.
ans = struct with fields:
solver: 'bvp4c' x: [0 0.0076 0.0152 0.0227 0.0303 0.0370 0.0438 0.0513 0.0589 0.0665 0.0741 0.0816 0.0892 0.0968 0.1044 0.1120 0.1195 0.1271 0.1347 0.1423 0.1498 0.1574 … ] (1×352 double) y: [8×352 double] yp: [8×352 double] stats: [1×1 struct]
%rray indices must be positive integers or logical values.
%Error in Tarek7_2025/projfun (line 92)
% dy(2) = (((rhot/mut)*(a*(muf/rhof)^0.5*(E*F+E^2)+a*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+R*a(muf/rhof)*sin(gamma)*(vb/(rhof*bef))*t)));
%code
function sol= proj
clc;clf;clear;
%Relation of base fluid
rhof=997.1*10^-3;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;
alfaf=kf/(rhof*cpf);
bef=21*10^-5;
ky=muf/rhof;
disp('ky');disp((muf/rhof));
%sigf=0.05*10^-8;
%Ag
ph1=0.01;
rho1=10500*10^-3;
cp1=235*10^4;
k1=429*10^5;be1=21*10^-5;
%sig1=0.74*10^-2;
%copper
ph2=0.01;
rho2=8933*10^-3;
cp2=385*10^4;
k2=400*10^5;
%sig2=5.96*10^-1;
be2=1.67*10^-5;
%Alumina
ph3=0.01;
rho3=3970*10^-3;
cp3=765*10^4;
k3=40*10^5;
be3=0.85*10^-5;
%sig3=3.5*10^-1;
%Relation of ternary hyprid
kn=kf*((k3+2*kf-2*ph3*(kf-k3))/(k3+2*kf+ph3*(kf-k3)));
kh=kn*((k2+2*kn-2*ph2*(kn-k2))/(k2+2*kn+ph2*(kn-k2)));
kt=kh*((k1+2*kh-2*ph1*(kh-k1))/(k1+2*kh+ph1*(kh-k1)));
mut= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5);
rhot=(1-ph1)*((1-ph2)*((1-ph3)+ph3*(rho3/rhof))+ph2*(rho2/rhof))+ph1*(rho1/rhof);
%vt=rhot*cpt
vt =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*cp3)/(rhof*cpf)))+ph2*((rho2*cp2)/(rhof*cpf)))+ph1*((rho1*cp1)/(rhof*cpf));
%disp('vt');disp(vt);
%vb=rho*betb
vb =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*be3)/(rhof*bef)))+ph2*((rho2*be2)/(rhof*bef)))+ph1*((rho1*be1)/(rhof*bef));
%disp('vb');disp(vb);disp(ky);
myLegend1 = {};myLegend2 = {};
rr = [1 2 3 4]
numR = numel(rr);
m = linspace(0,1);
a=0.001;b=0.001;p=.001/((1-0.01)*(mut/muf)*(rhof/rhot));
Ec=10;at=0.01;gamma=pi/4;
prf=6.9;d=0.009;e=0.1;Rd=0.45;
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
disp('coe');disp((mut/muf)*(rhof/rhot));
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1];options =bvpset('stats','on','RelTol',1e-5);
Z = zeros(numR, length(m));
for i = 1:numR
R= rr(i);
if i==1
solinit = bvpinit(m, y0);
else
guess = @(x)interp1(sol.x,(sol.y).',x);
solinit = bvpinit(sol.x,guess);
end
sol = bvp4c(@projfun,@projbc, solinit, options);
N=(kt/kf)*(1+Rd)*d*R*(1/(muf/rhof)^2)*(sol.y(7,:).^2)+((mut/kf)*(1/Ti)*(1/at^2))*R*(d*(muf/rhof)*(sol.y(2,:).^2)+e*(sol.y(4,:).^2)+R*d*p*(sol.y(1,:).^2));
Z(i,:) = interp1(sol.x,N,m);
end
[X, Y] = meshgrid(m, rr);
shading interp;
surf(X, Y, Z);
xlabel('eta');
ylabel('Gr');
zlabel('NG');
title('Variation of velocity with Grashof number,Gr in 3D' );
grid on
shading flat;
colorbar;
function dy= projfun(~,y)
dy= zeros(8,1);
% alignComments
E = y(1);
dE = y(2);
F = y(3);
dF= y(4);
w = y(5);
dw=y(6);
t = y(7);
dt = y(8);
dy(1) = dE;
dy(2) = (((rhot/mut)*(a*(muf/rhof)^0.5*(E*F+E^2)+a*(muf/rhof)*w*dE-(mut/muf)*(rhof/rhot)*p*(1-0.01)*E+R*a*(muf/rhof)*sin(gamma)*(vb/(rhof*bef))*t)));
dy(3) = dF;
dy(4) = (((rhot/mut)*(b*(muf/rhof)^0.5*(F^2+F*E)+(muf/rhof)*b^0.5*a^(1.5)*dF)));
dy(5) =-(a*F+b*E);
dy(6) = (((rhot/mut)*((muf/rhof)^0.5*w*dw+R*b*(muf/rhof)*cos(gamma)*(vb/(rhof*bef))*t)));
dy(7) = dt;
dy(8)=prf*(1/(kt/kf))*(1/(1+((prf*Rd)/((kt/kf)))))*((vt/(rhof*cpf))*(muf/rhof)^0.5*w*dt-(mut/muf)*Ec*1*dw^2) ;
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5);
ya(6);
ya(7)+1-(1/0.9)*ya(8);
yb(1)-0.01;
yb(3);
yb(7);
% yb(7);
];
end
.

Plus de réponses (0)

Catégories

En savoir plus sur Eigenvalue Problems dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by