Please help me to run this simple code

43 次查看（过去 30 天）
T K
T K 2026-2-13，18:18
评论： dpb 2026-2-13，22:28
Ran in:
proj()
rr = 1×3
0 0.3000 0.6000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Array indices must be positive integers or logical values.

Error in solution>proj/projfun (line 48)
dy(2)=(1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a(s^2+a1*h^2)*t-a2*s*h*a(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in bvparguments (line 96)
testODE = ode(x1,y1,odeExtras{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in bvp4c (line 119)
bvparguments(solver_name,ode,bc,solinit,options,varargin);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in solution>proj (line 19)
sol= bvp4c(@projfun,@projbc,solinit,options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
function sol= proj
clc;clf;clear;
myLegend1 = {};
myLegend2 = {};
k0=386; ce=3.831*10^2; mu=38.6*10^9;alfat=1.78*10^-5; rho=89.54*10^2; lamda=77.6*10^9;taw=0.5;Tnot=2.93*10^2;
c0=sqrt((lamda+2*mu)/(rho)); Betanot=(3*lamda+2*mu)*alfat; a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);
rr = [0 0.3 0.6 ]
for i =1:numel(rr)
a= rr(i);
s=0.5;h=0.1;
y0 = [1,0,0,0,1,0,0];
options =bvpset('stats','on','RelTol',1e-4);
m = linspace(0,10);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
figure(1)
plot(sol.x,(sol.y(5,:)))
grid on,hold on
myLegend1{i}=['alfa= ',num2str(rr(i))];
figure(2)
plot(sol.x,(sol.y(6,:)))
title('Temperature')
grid on,hold on
myLegend2{i}=['alfa= ',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(7,1);
v = y(1);
dv = y(2);
u=y(3);
du=y(4);
t=y(5);
dt=y(6);
ddt=y(7);
dy(1) = dv;
dy(2)=(1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a(s^2+a1*h^2)*t-a2*s*h*a(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h));
dy(3)=du;
dy(4) =(1/(a*(s^2+a1*h^2-(2*x+2+1).^2))-(s^2+a1*h^2-(2*x+2+1).^2))*((a2*s*h-a*a2*s*h*t)*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a(s^2+a1*h^2)*t-a2*s*h*a(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+(((a*2*(x+1)*t)-2*(x+1)-(a*s^2+a*a1*h^2)*t)*du)+(2*a*a3*s*t-a3*s)*dt-(a*a4*s*h+a*a1*s*h)*dt*dv);
dy(5)=dt;
dy(6)=ddt;
dy(7)=(1/(a*(a5*s^2*2*(x+1+1)+a5*h^2*2*(x+1+1))*t-(a5*s^2*2*(x+1+1)+a5*h^2*2*(x+1)-a8*(2*(x+1+1))^3)))*((s^2+h^2+2*a5*h*2-a6*(2*(x+1)).^2-3*a8*2*(x+1+1)*2*(x+1))*ddt);
end
end
function res= projbc(ya,yb)
res= [ya(1);
ya(3);
ya(5)-1;
ya(7);
yb(1);
yb(3);
yb(5);
];
end
  3 个评论
显示 1更早的评论
dpb
dpb about 3 hours 前
编辑：dpb about 3 hours 前
@T K
dy(2)=(1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a(s^2+a1*h^2)*t-a2*s*h*a(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h));
Format your code by wrapping long lines and appropriate spacing, etc., so it's possible to read...there's almost no chance of being able to parse those lines as written without missing something...making shorter, temporary variables wouldn't necessarily be a bad idea, either as one way to reduce the apparent complexity.
Stephen23
Stephen23 about 3 hours 前
And get rid of cargo-cult code like this: clc;clf;clear;
What do you think clearing an empty function workspace will do ?

请先登录，再进行评论。

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

采纳的回答

Stephen23
Stephen23 about 3 hours 前
Ran in:
You're missing the multiplication operator * in several places. Here's the corrected code with the missing * operators added:
proj()
The solution was obtained on a mesh of 100 points. The maximum residual is 4.453e-06. There were 1303 calls to the ODE function. There were 25 calls to the BC function. The solution was obtained on a mesh of 100 points. The maximum residual is 3.046e-06. There were 1502 calls to the ODE function. There were 26 calls to the BC function. The solution was obtained on a mesh of 100 points. The maximum residual is 5.754e-06. There were 1701 calls to the ODE function. There were 27 calls to the BC function.
ans = struct with fields:
solver: 'bvp4c' x: [0 0.1010 0.2020 0.3030 0.4040 0.5051 0.6061 0.7071 0.8081 0.9091 1.0101 1.1111 1.2121 1.3131 1.4141 1.5152 1.6162 1.7172 1.8182 1.9192 2.0202 2.1212 … ] (1×100 double) y: [7×100 double] yp: [7×100 double] stats: [1×1 struct]
function sol= proj
myLegend1 = {};
myLegend2 = {};
k0=386; ce=3.831*10^2; mu=38.6*10^9;alfat=1.78*10^-5; rho=89.54*10^2; lamda=77.6*10^9;taw=0.5;Tnot=2.93*10^2;
c0=sqrt((lamda+2*mu)/(rho)); Betanot=(3*lamda+2*mu)*alfat; a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);
rr = [0 0.3 0.6 ];
for i =1:numel(rr)
a= rr(i);
s=0.5;h=0.1;
y0 = [1,0,0,0,1,0,0];
options =bvpset('stats','on','RelTol',1e-4);
m = linspace(0,10);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
figure(1)
plot(sol.x,(sol.y(5,:)))
grid on,hold on
myLegend1{i}=['alfa= ',num2str(rr(i))];
figure(2)
plot(sol.x,(sol.y(6,:)))
title('Temperature')
grid on,hold on
myLegend2{i}=['alfa= ',num2str(rr(i))];
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
hold on
function dy = projfun(x,y)
dy = zeros(7,1);
v = y(1);
dv = y(2);
u = y(3);
du = y(4);
t = y(5);
dt = y(6);
ddt = y(7);
dy(1) = dv;
dy(2)=(1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h));
dy(3)=du;
dy(4) =(1/(a*(s^2+a1*h^2-(2*x+2+1).^2)-(s^2+a1*h^2-(2*x+2+1).^2)))*((a2*s*h-a*a2*s*h*t)*((1/((((s^2+a1*h^2-(2*x+2+1).^2)*(1-a*t))*((h^2+a1*s^2-(2*(x+1+1)).^2)*(1-a*t)))-((a2*s*h)*(1-a*t))^2))*(((a2*s*h*a*2*(x+1)+2*(x+1)*a*a2*s*h)*t-a2*s*h*a*2*(x+1)*t^2-a2*s*h*2*(x+1)+(a*a2*s*h*a*(s^2+a1*h^2)*t-a2*s*h*a*(s^2+a1*h^2))*dt)*du+((s^2+a1*h^2-2*(x+1))*a*(h^2+a1*s^2)-(a*2*(x+1)*a*(a4*s*h+a1*s*h))+((a*a2*s*h)*a*(a4*s*h+a1*s*h)-(a*(s^2+a1*h^2-2*(x+1)))*(a*(h^2+a1*s^2))))*dt*dv+(a2*s*h*2*a*a3*s+a3*s*a*a2*s*h+2*a*a3*h*(s^2+a1*h^2-4*(x+1+1)^2)+a*(s^2+a1*h^2-(2*(x+1+1))^2)*(a3*h))*t*dt-(a2*s*h*a3*s)*dt-((s^2+a1*h^2-(2*(x+1))^2)*a3*h)*dt-dt*t^2*(a*a2*s*h*2*a*a3*s+a*(s^2+a1*h^2-(2*(x+1))^2)*2*a*a3*h)))+(((a*2*(x+1)*t)-2*(x+1)-(a*s^2+a*a1*h^2)*t)*du)+(2*a*a3*s*t-a3*s)*dt-(a*a4*s*h+a*a1*s*h)*dt*dv);
dy(5)=dt;
dy(6)=ddt;
dy(7)=(1/(a*(a5*s^2*2*(x+1+1)+a5*h^2*2*(x+1+1))*t-(a5*s^2*2*(x+1+1)+a5*h^2*2*(x+1)-a8*(2*(x+1+1))^3)))*((s^2+h^2+2*a5*h*2-a6*(2*(x+1)).^2-3*a8*2*(x+1+1)*2*(x+1))*ddt);
end
function res = projbc(ya,yb)
res = [ya(1);
ya(3);
ya(5)-1;
ya(7);
yb(1);
yb(3);
yb(5);
];
end
end
  2 个评论
dpb
dpb about 1 hour 前
@Stephen23 has far more patience (and better eyesight, undoubtedly) than I... <vbg>

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by