Please help me to run this simple code

39 次查看（过去 30 天）
T K
T K 2026-1-16，20:27
编辑： Torsten 2026-1-17，21:15
Ran in:
%ErrorError using surf (line 71)
X, Y, Z, and C cannot be complex.
Error in Untitled (line 66)
surf(X, Y, Z);
proj
ky 0.0100 The solution was obtained on a mesh of 754 points. The maximum residual is 2.042e-06. There were 39912 calls to the ODE function. There were 274 calls to the BC function. The solution was obtained on a mesh of 192 points. The maximum residual is 8.109e-06. There were 6379 calls to the ODE function. There were 110 calls to the BC function.
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.184664, while requested accuracy is 1e-05.
The solution was obtained on a mesh of 2674 points. The maximum residual is 1.847e-01. There were 85426 calls to the ODE function. There were 341 calls to the BC function.
Error using surf (line 71)
X, Y, Z, and C cannot be complex.

Error in solution>proj (line 69)
surf(X, Y, Z);
^^^^^^^^^^^^^^
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;
disp('ky');
disp((muf/rhof));
%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);
sigt = sigf + (3 * ((ph1 * sig1 + ph2 * sig2) - sigf * (ph1 + ph2)) /(((ph1 * sig1 + ph2 * sig2) / (sigf * (ph1 + ph2))) + 2 - ((ph1 * sig1 + ph2 * sig2) / sigf) + (ph1 + ph2)));
%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)))+(1-ph1)*ph1*((rho1*be1)/(rhof*bef));
myLegend1 = {};myLegend2 = {};
rr = [0.5 1 1.5];
numn = numel(rr);
m = linspace(0, 7);
R=3;M=2;Qh=0.1;
Rd=1;Pr=6.9;
y0 = [0,1,0,1,0,0,0,0];options =bvpset('stats','on','RelTol',1e-5);
solinit = bvpinit(m,y0);
Z = zeros(numn, length(m));
% Solve the BVP for each Prf
for i = 1:numn
n = rr(i);
solinit = bvpinit(m, y0);
sol = bvp4c(@projfun, @projbc, solinit, options);
Z(i, :) = deval(sol,m,1); % Store the z-axis data
end
[X, Y] = meshgrid(m, rr);
figure;
surf(X, Y, Z);
xlabel('x');
ylabel('Prf');
zlabel('Solution y(6,:)');
title('Surface Plot of Solution');
grid on;
i=i+1;
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
hold on
function dy = projfun(~, 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) =(R^((1-n)/2)*(muf/mut)*(rhot/rhof)*(f^2-g^2+R*h*df+(sigt/sigf)*M*f))^(1/n);
dy(3) = dg;
dy(4) = (R^((1-n)/2)*(muf/mut)*(rhot/rhof)*(f*g+R*dg+(sigt/sigf)*M*g))^(1/n);
dy(5) =dh ;
dy(6) =(R^((-n))*((muf/mut)*(rhot/rhof)*(2*h*f+R*h*dh)-2*h))^(1/n);
dy(7) =dt;
dy(8)=(1/(1+(4/3)*Rd*(1/(kt/kf))))*(1/R)*(-t-Qh*Pr*(kf/kt)*t+Pr*((vt/(rhof*cpf))*(kf/kt)*(f*t+R*h*dt)));
end
end
function res= projbc(ya,yb)
res= [ya(1)-0.3;
ya(3)-1;
ya(5)-0.1;
ya(7)-1;
yb(1)-0.3;
yb(3);
yb(5);
yb(7);
];
end

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

回答（1 个）

Torsten
Torsten 2026-1-16，20:52
编辑：Torsten 2026-1-17，21:15
Your solution for rr(3) becomes complex-valued.
For complex-valued functions, you can plot e.g. their real part, their imaginary part or their absolute value:
surf(X,Y,real(Z))
surf(X,Y,imag(Z))
surf(X,Y,abs(Z))
You try to plot the solution y(1,:). Why do you use a legend that tells you plot y(6,:) ?
Do you really mean
dy(6) =(R^((-n))*((muf/mut)*(rhot/rhof)*(2*h*f+R*h*dh)-2*h))^(1/n);
If I look at dy(2) and dy(4), it could be
dy(6) =(R^((-n)*((muf/mut)*(rhot/rhof)*(2*h*f+R*h*dh)-2*h)))^(1/n);

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by