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

6 Ansichten (letzte 30 Tage)
Tarek
Tarek am 23 Jul. 2025
Beantwortet: Star Strider am 23 Jul. 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

Melden Sie sich an, um diese Frage zu beantworten.

Akzeptierte Antwort

Star Strider
Star Strider am 23 Jul. 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
.

Weitere Antworten (0)

Kategorien

Mehr zu Eigenvalue Problems finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by