Please help me to run this simple code

Question

T K 2025-12-25

0
链接

此问题的直接链接

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

编辑： Torsten 2025-12-26

在 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 个评论
显示 4更早的评论隐藏 4更早的评论

T K 2025-12-25

编辑：Torsten 2025-12-26

在 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 2025-12-26

编辑：Torsten 2025-12-26

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.

请先登录，再进行评论。

请先登录，再回答此问题。

Please help me to run this simple code

6 个评论
显示 4更早的评论隐藏 4更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Please help me to run this simple code

6 个评论 显示 4更早的评论隐藏 4更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论