Please help me to run this simple code

Question

T K on 25 Dec 2025 at 14:02

0
Link

Direct link to this question

https://uk.mathworks.com/matlabcentral/answers/2181971-please-help-me-to-run-this-simple-code

Edited: Torsten on 26 Dec 2025 at 1:16

Open in MATLAB Online

% Error

Array indices must be positive integers or logical values.

Error in Untitled2 (line 45)

kh=kf*((k1+2*kf-2*ph1(kf-k1))/(k1+2*kf+ph1*(kf-k1)));

% code
proj()
Array indices must be positive integers or logical values.

Error in solution>proj (line 47)
kh=kf*((k1+2*kf-2*ph1(kf-k1))/(k1+2*kf+ph1*(kf-k1)));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
function sol = proj
    clc; clf; clear;
    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;
sigf=0.05*10^-8;
ky=muf/rhof;
%Titanium oxide
ph4=0.01;
rho4=4250*10^-3;
cp4=711*10^4;
k4=8.953*10^5;
sig4=2.6*10*10^-1;
%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
kh=kf*((k1+2*kf-2*ph1(kf-k1))/(k1+2*kf+ph1*(kf-k1)));
kh=kn*((k2+2*kn-2*ph2(kn-k2))/(k2+2*kn+ph2*(kn-k2)));
kt=kh*((k3+2*kh-2*ph3(kh-k3))/(k3+2*kh+ph3*(kh-k3)));
kT=kt*((k4+2*kt-2*ph4*(kt-k4))/(k4+2*kt+ph4*(kt-k4)));
muT= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5*(1-ph4)^2.5);
rhoT=rhof*((1-ph4)*((1-ph3)*((1-ph2)*((1-ph1)+ph1*((ph1*rho1*(1/rhof))))+(ph2*rho2*(1/rhof))))+(ph3*rho3*(1/rhof))+ph4*rho4*(1/rhof));
sign = sigf*(1+(3*(sigf-1)*ph1)/((sigf+2)-(sigf-1)*ph1));
sigh = sign*(1+(3*(sign-1)*ph2)/((sign+2)-(sign-1)*ph2));
sigt = sigh*(1+(3*(sigh-1)*ph3)/((sigh+2)-(sigh-1)*ph3));
sigT = sigt*(1+(3*(sigt-1)*ph4)/((sigt+2)-(sigt-1)*ph4));
%vt=rhot*cpt
VT=(rhof*cpf)*((1-ph4)*((1-ph3)*((1-ph2)*((1-ph1)+ph1*((ph1*cp1*(1/rhof*cpf))))+(ph2*cp2*(1/rhof*cpf))))+(ph3*cp3*(1/rhof*cpf))+ph4*cp4*(1/rhof*cpf));
%disp('vt');disp(vt);
%vb=rho*betb
vT =(1-ph1)*((1-ph2)*((1-ph3)+ph3*((rho3*be3)/(rhof*bef)))+ph2*((rho2*be2)/(rhof*bef)))+(1-ph1)*ph1*((rho1*be1)/(rhof*bef));  
 myLegend1 = {};myLegend2 = {};
rr = [0.1 0.3 0.5 0.7]
for i =1:numel(rr)
    Rd= rr(i)
M=0.5;
    R=1;Pr=6.9;
    m = linspace(0,1);
y0 = [0,1,0,1,0,0,0,0];options =bvpset('stats','on','RelTol',1e-5);
   solinit = bvpinit(m,y0);
      sol= bvp4c(@projfun,@projbc,solinit,options);
          figure(1)
    plot(sol.x,(sol.y(1,:)))
    grid on,hold on
    myLegend1{i}=['n= ',num2str(rr(i))];
    
    figure(2)
 plot(sol.x,(sol.y(7,:)))
    grid on,hold on
    myLegend2{i}=['n= ',num2str(rr(i))];
    
    i=i+1;
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)                                                                                                                                                                                                                                                                               
    hold on
    function dy = projfun(x, y)
      dy= zeros(8,1);
        % alignComments
        f = y(1);
        df = y(2);
         g = y(3);
        dg= y(4);
        h= y(5);
        dh = y(6);
         t = y(7);
       dt=y(8);
               dy(1) = df; 
        dy(2) =(1/(1+x^2))*(3*x*df+R*((muT/muf))*(f^2-f*df+h*df-g^2+(sigT/sigf)*M*f));
         dy(3) = dg;
     dy(4) = (1/(1+x^2))*(3*x*dg+R*((muT/muf))*(2*f*g-f*dg+h*dg+(sigT/sigf)*M*g));
         dy(5) =dh ;
      dy(6) =(1/(1+x^2))*(2*x*dh-h+R*((muT/muf))*(f*dh-h*dh-f*h));
         dy(7) =dt;
       dy(8)=(1/(1+x^2+(4/3)*Rd*(1/(kT/kf))))*(((x-4)*dt)-4*t-Pr*R*((vT/(rhof*cpf))*(kf/kT))*(x*f*dt-2*f*t-h*dt));
                     
    end
end
function res= projbc(ya,yb)
res= [ya(1)-0.1;
  ya(3)-1;
  ya(5)-0.1;
    ya(7)-1;
    
    yb(1)-0.1;
        yb(3);
        yb(5);
     yb(7);
       ];
       end

6 Comments
Show 4 older commentsHide 4 older comments

T K on 25 Dec 2025 at 19:30

Edited: Torsten on 26 Dec 2025 at 1:11

Open in MATLAB Online

Iam added multiplication sign but the code not run.Pleae help me

% Error using bvp4c (line 251)

Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in proj (line 73)

sol= bvp4c(@projfun,@projbc,solinit,options);

% code

proj()

rr = 1×4

0.1000 0.3000 0.5000 0.7000

Rd = 0.1000

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.240883, while requested accuracy is 1e-05.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 9.13308e-94.

The solution was obtained on a mesh of 2674 points. The maximum residual is 2.409e-01. There were 107407 calls to the ODE function. There were 354 calls to the BC function.

Rd = 0.3000

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.240883, while requested accuracy is 1e-05.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 1.43969e-92.

The solution was obtained on a mesh of 2674 points. The maximum residual is 2.409e-01. There were 107374 calls to the ODE function. There were 354 calls to the BC function.

Rd = 0.5000

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.240883, while requested accuracy is 1e-05.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 1.26505e-91.

The solution was obtained on a mesh of 2674 points. The maximum residual is 2.409e-01. There were 107369 calls to the ODE function. There were 354 calls to the BC function.

Rd = 0.7000

Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.240883, while requested accuracy is 1e-05.

Warning: Results may be inaccurate because of an ill-conditioned Jacobian whose reciprocal condition number is 9.33478e-91.

The solution was obtained on a mesh of 2674 points. The maximum residual is 2.409e-01. There were 107331 calls to the ODE function. There were 354 calls to the BC function.

ans = struct with fields:

solver: 'bvp4c' x: [0 0.0011 0.0022 0.0034 0.0045 0.0056 0.0067 0.0079 0.0090 0.0101 0.0112 0.0123 0.0135 0.0146 0.0157 0.0168 0.0180 0.0191 0.0202 0.0213 0.0224 0.0236 … ] (1×892 double) y: [8×892 double] yp: [8×892 double] stats: [1×1 struct]

function sol = proj

clc; clf; clear;

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;

sigf=0.05*10^-8;

ky=muf/rhof;

%Titanium oxide

ph4=0.01;

rho4=4250*10^-3;

cp4=711*10^4;

k4=8.953*10^5;

sig4=2.6*10*10^-1;

%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*((k1+2*kf-2*ph1*(kf-k1))/(k1+2*kf+ph1*(kf-k1)));

kh=kn*((k2+2*kn-2*ph2*(kn-k2))/(k2+2*kn+ph2*(kn-k2)));

kt=kh*((k3+2*kh-2*ph3*(kh-k3))/(k3+2*kh+ph3*(kh-k3)));

Km=kt*((k4+2*kt-2*ph4*(kt-k4))/(k4+2*kt+ph4*(kt-k4)));

muT= muf/((1-ph1)^2.5*(1-ph2)^2.5*(1-ph3)^2.5*(1-ph4)^2.5);

rhoT=rhof*((1-ph4)*((1-ph3)*((1-ph2)*((1-ph1)+ph1*((ph1*rho1*(1/rhof))))+(ph2*rho2*(1/rhof))))+(ph3*rho3*(1/rhof))+ph4*rho4*(1/rhof));

sign = sigf*(1+(3*(sigf-1)*ph1)/((sigf+2)-(sigf-1)*ph1));

sigh = sign*(1+(3*(sign-1)*ph2)/((sign+2)-(sign-1)*ph2));

sigt = sigh*(1+(3*(sigh-1)*ph3)/((sigh+2)-(sigh-1)*ph3));

sigT = sigt*(1+(3*(sigt-1)*ph4)/((sigt+2)-(sigt-1)*ph4));

%vt=rhot*cpt

vt=(rhof*cpf)*((1-ph4)*((1-ph3)*((1-ph2)*((1-ph1)+ph1*((ph1*cp1*(1/rhof*cpf))))+(ph2*cp2*(1/rhof*cpf))))+(ph3*cp3*(1/rhof*cpf))+ph4*cp4*(1/rhof*cpf));

myLegend1 = {};myLegend2 = {};

rr = [0.1 0.3 0.5 0.7]

for i =1:numel(rr)

Rd= rr(i)

M=0.5;

R=1;Pr=6.9;

m = linspace(0,1);

y0 = [0,1,0,1,0,0,0,0];options =bvpset('stats','on','RelTol',1e-5);

solinit = bvpinit(m,y0);

sol= bvp4c(@projfun,@projbc,solinit,options);

figure(1)

plot(sol.x,(sol.y(1,:)))

grid on,hold on

myLegend1{i}=['n= ',num2str(rr(i))];

figure(2)

plot(sol.x,(sol.y(7,:)))

grid on,hold on

myLegend2{i}=['n= ',num2str(rr(i))];

i=i+1;

end

figure(1)

legend(myLegend1)

hold on

figure(2)

legend(myLegend2)

hold on

function dy = projfun(x, y)

dy= zeros(8,1);

% alignComments

f = y(1);

df = y(2);

g = y(3);

dg= y(4);

h= y(5);

dh = y(6);

t = y(7);

dt=y(8);

dy(1) = df;

dy(2) =(1/(1+x^2))*(3*x*df+R*((muT/muf))*(f^2-f*df+h*df-g^2+(sigT/sigf)*M*f));

dy(3) = dg;

dy(4) = (1/(1+x^2))*(3*x*dg+R*((muT/muf))*(2*f*g-f*dg+h*dg+(sigT/sigf)*M*g));

dy(5) =dh ;

dy(6) =(1/(1+x^2))*(2*x*dh-h+R*((muT/muf))*(f*dh-h*dh-f*h));

dy(7) =dt;

dy(8)=(1/(1+x^2+(4/3)*Rd*(1/(Km/kf))))*(((x-4)*dt)-4*t-Pr*R*((vt/(rhof*cpf))*(kf/Km))*(x*f*dt-2*f*t-h*dt));

end

function res= projbc(ya,yb)

res= [ya(1)-0.1;

ya(3)-1;

ya(5)-0.1;

ya(7)-1;

yb(1)-0.1;

yb(3);

yb(5);

yb(7);

];

end

Torsten on 26 Dec 2025 at 1:15

Edited: Torsten on 26 Dec 2025 at 1:16

The warning message of the recent MATLAB release is different (see above).

MATLAB cannot solve your equations. Maybe there are mistakes in the expressions or in the parameters, maybe the equations are just very difficult to solve - we cannot tell.

Sign in to comment.

Sign in to answer this question.

Please help me to run this simple code

6 Comments
Show 4 older commentsHide 4 older comments

Answers (0)

See Also

Categories

Tags

Community Treasure Hunt

Please help me to run this simple code

6 Comments Show 4 older commentsHide 4 older comments

Answers (0)

See Also

Categories

Tags

Community Treasure Hunt

6 Comments
Show 4 older commentsHide 4 older comments