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

Question

Tarek on 23 Jul 2025

0
Link

Direct link to this question

https://ch.mathworks.com/matlabcentral/answers/2178664-please-help-me-to-run-the-surface-figure-from-the-attached-code

Answered: Star Strider on 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

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

Sign in to answer this question.

Answer 1

Star Strider on 23 Jul 2025

0
Link

Direct link to this answer

https://ch.mathworks.com/matlabcentral/answers/2178664-please-help-me-to-run-the-surface-figure-from-the-attached-code#answer_1568325

Open in MATLAB Online

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

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

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

.

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

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

0 Comments
Show -2 older commentsHide -2 older comments

Accepted Answer

0 Comments
Show -2 older commentsHide -2 older comments

More Answers (0)

See Also

Categories

Tags

Community Treasure Hunt

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

0 Comments Show -2 older commentsHide -2 older comments

Accepted Answer

0 Comments Show -2 older commentsHide -2 older comments

More Answers (0)

See Also

Categories

Tags

Community Treasure Hunt

0 Comments
Show -2 older commentsHide -2 older comments

0 Comments
Show -2 older commentsHide -2 older comments