Pleas help me to run this simple code

T K
2026 2 17
1 个回答
更新时间：2026 2 17
15 次查看（30 天）
Ran in:

0 个投票

Ran in:
proj()
Error using bvp4c (line 196)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in solution>proj (line 17)
sol= bvp4c(@projfun,@projbc,solinit,options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
function sol= proj
myLegend1 = {};
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.2 0.4 0.6];
for i =1:numel(rr)
a= rr(i);
s=1;h=1;
y0 = [1,0,0,1,1,0,0];
options =bvpset('stats','on','RelTol',1e-4);
m = linspace(0,5);
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))];
end
figure(1)
legend(myLegend1)
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-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+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)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv);
dy(5)=dt;
dy(6)=ddt;
dy(7)=(1/(a*(a5*s^2*(x+2)+a5*h^2*(x+2)-a8*(x+2)^3)*t+a5*s^2*(x+2)^2+a5*h^2*(x+2)-a8*(x+2)^3))*(-(2*a7*(x+2)-2*a*a11*(x+2))*((1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+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)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv))+(s^2+h^2+2*a5*t-a6*(x+2)^2)*ddt-(a*(s^2+h^2+4*a5-a6*(x+2)^2))*t*ddt-(a*(s^2+h^2+a5+a5))*dt*dt-a*(x+2)*(a5*s^2+a5*h^2+a10*(x+2)^2)*dt*ddt-(2*a7*(x+2)-2*a*a11*(x+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)))+2*a*a7*(x+2)*t*((1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+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)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv))+(2*a*a7*(x+2))*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))));
end
function res = projbc(ya,yb)
res = [ya(1);
ya(3);
ya(5)-1;
ya(7);
yb(1);
yb(3);
yb(5);
];
end
end

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

回答（1 个）

Walter Roberson
Walter Roberson 2026-2-17
移动：Walter Roberson 2026-2-17
Ran in:

0 个投票

Ran in:
Your equations are going to +/- infinity in the last component.
Below, max_dy is printed out each time max(abs(dy)) increases
proj()
max dy increased to -1.03787e+09 at component #7 max dy increased to -1.04008e+09 at component #7 max dy increased to -1.0901e+09 at component #7 max dy increased to -1.17508e+09 at component #7 max dy increased to -1.31614e+09 at component #7 max dy increased to -1.56033e+09 at component #7 max dy increased to -2.03538e+09 at component #7 max dy increased to -3.24771e+09 at component #7 max dy increased to -1.14454e+10 at component #7 max dy increased to 1.51966e+10 at component #7 max dy increased to 1.62624e+10 at component #7 max dy increased to -9.60812e+13 at component #7 max dy increased to -2.15787e+14 at component #7 max dy increased to -3.60165e+14 at component #7 max dy increased to -5.29562e+14 at component #7 max dy increased to -7.23281e+14 at component #7 max dy increased to -9.3935e+14 at component #7 max dy increased to -1.17413e+15 at component #7 max dy increased to -1.42209e+15 at component #7 max dy increased to -1.67524e+15 at component #7 max dy increased to -1.92399e+15 at component #7 max dy increased to -2.14796e+15 at component #7 max dy increased to -2.22608e+15 at component #7 max dy increased to -2.48603e+15 at component #7 max dy increased to -2.64374e+15 at component #7 max dy increased to 3.73381e+15 at component #7 max dy increased to 6.47483e+15 at component #7 max dy increased to 9.99645e+15 at component #7 max dy increased to 1.44491e+16 at component #7 max dy increased to 2.00057e+16 at component #7 max dy increased to 2.68613e+16 at component #7 max dy increased to 3.52405e+16 at component #7 max dy increased to 4.53984e+16 at component #7 max dy increased to 5.76213e+16 at component #7 max dy increased to 7.22314e+16 at component #7 max dy increased to 8.95989e+16 at component #7 max dy increased to 1.10133e+17 at component #7 max dy increased to 1.343e+17 at component #7 max dy increased to 1.62614e+17 at component #7 max dy increased to 1.95671e+17 at component #7 max dy increased to 2.34109e+17 at component #7 max dy increased to 2.78674e+17 at component #7 max dy increased to 3.30185e+17 at component #7 max dy increased to 3.89547e+17 at component #7 max dy increased to 4.57799e+17 at component #7 max dy increased to 5.36075e+17 at component #7 max dy increased to 6.25684e+17 at component #7 max dy increased to 7.28032e+17 at component #7 max dy increased to 8.44769e+17 at component #7 max dy increased to 9.77667e+17 at component #7 max dy increased to 1.12873e+18 at component #7 max dy increased to 1.30023e+18 at component #7 max dy increased to 1.49472e+18 at component #7 max dy increased to 1.71502e+18 at component #7 max dy increased to 1.96437e+18 at component #7 max dy increased to 2.24637e+18 at component #7 max dy increased to 2.56499e+18 at component #7 max dy increased to 2.92487e+18 at component #7 max dy increased to 3.3311e+18 at component #7 max dy increased to 3.78963e+18 at component #7 max dy increased to 4.3069e+18 at component #7 max dy increased to 4.89047e+18 at component #7 max dy increased to 5.54868e+18 at component #7 max dy increased to 6.2911e+18 at component #7 max dy increased to 7.12886e+18 at component #7 max dy increased to 8.07401e+18 at component #7 max dy increased to 9.14104e+18 at component #7 max dy increased to 1.0346e+19 at component #7 max dy increased to 1.17075e+19 at component #7 max dy increased to 1.32462e+19 at component #7 max dy increased to 1.49867e+19 at component #7 max dy increased to 1.69571e+19 at component #7 max dy increased to 1.91888e+19 at component #7 max dy increased to 2.17194e+19 at component #7 max dy increased to 2.45911e+19 at component #7 max dy increased to 2.78533e+19 at component #7 max dy increased to 3.15632e+19 at component #7 max dy increased to 3.57857e+19 at component #7 max dy increased to 4.06003e+19 at component #7 max dy increased to 4.60941e+19 at component #7 max dy increased to 5.2373e+19 at component #7 max dy increased to 5.9558e+19 at component #7 max dy increased to 6.77929e+19 at component #7 max dy increased to 7.72436e+19 at component #7 max dy increased to 8.8113e+19 at component #7 max dy increased to 1.00628e+20 at component #7 max dy increased to 1.15065e+20 at component #7 max dy increased to 1.31758e+20 at component #7 max dy increased to 1.51085e+20 at component #7 max dy increased to 1.73516e+20 at component #7 max dy increased to 1.99599e+20 at component #7 max dy increased to 2.30006e+20 at component #7 max dy increased to 2.65529e+20 at component #7 max dy increased to 3.07146e+20 at component #7 max dy increased to 3.56019e+20 at component #7 max dy increased to 4.13585e+20 at component #7 max dy increased to 4.81598e+20 at component #7 max dy increased to 5.62216e+20 at component #7 max dy increased to 6.58104e+20 at component #7 max dy increased to 7.72611e+20 at component #7 max dy increased to 9.09916e+20 at component #7 max dy increased to 1.07528e+21 at component #7 max dy increased to 1.27549e+21 at component #7 max dy increased to 1.51921e+21 at component #7 max dy increased to 4.7741e+22 at component #7 max dy increased to 6.71942e+22 at component #7 max dy increased to 9.10873e+22 at component #7 max dy increased to 1.19593e+23 at component #7 max dy increased to 1.52414e+23 at component #7 max dy increased to 1.88675e+23 at component #7 max dy increased to 2.2652e+23 at component #7 max dy increased to 2.62798e+23 at component #7 max dy increased to 2.91515e+23 at component #7 max dy increased to 2.99859e+23 at component #7 max dy increased to -5.75886e+23 at component #7 max dy increased to -1.14348e+24 at component #7 max dy increased to -1.95779e+24 at component #7 max dy increased to -3.09227e+24 at component #7 max dy increased to -4.63926e+24 at component #7 max dy increased to -6.71261e+24 at component #7 max dy increased to -9.45247e+24 at component #7 max dy increased to -1.30304e+25 at component #7 max dy increased to -1.76537e+25 at component #7 max dy increased to -2.35782e+25 at component #7 max dy increased to -3.10957e+25 at component #7 max dy increased to -4.05812e+25 at component #7 max dy increased to -5.24738e+25 at component #7 max dy increased to -6.72817e+25 at component #7 max dy increased to -8.56113e+25 at component #7 max dy increased to -1.08236e+26 at component #7 max dy increased to -1.35987e+26 at component #7 max dy increased to -1.69926e+26 at component #7 max dy increased to -2.1122e+26 at component #7 max dy increased to -2.61423e+26 at component #7 max dy increased to -3.22037e+26 at component #7 max dy increased to -3.95277e+26 at component #7 max dy increased to -4.83392e+26 at component #7 max dy increased to -5.89055e+26 at component #7 max dy increased to -7.15744e+26 at component #7 max dy increased to -8.67043e+26 at component #7 max dy increased to -1.04794e+27 at component #7 max dy increased to -1.2631e+27 at component #7 max dy increased to -1.51972e+27 at component #7 max dy increased to -1.82441e+27 at component #7 max dy increased to -2.18602e+27 at component #7 max dy increased to -2.61536e+27 at component #7 max dy increased to -3.12451e+27 at component #7 max dy increased to -3.72798e+27 at component #7 max dy increased to -4.44332e+27 at component #7 max dy increased to -5.2913e+27 at component #7 max dy increased to -6.29377e+27 at component #7 max dy increased to -7.48353e+27 at component #7 max dy increased to -8.89177e+27 at component #7 max dy increased to -1.05643e+28 at component #7 max dy increased to -1.2545e+28 at component #7 max dy increased to -1.48981e+28 at component #7 max dy increased to -1.76901e+28 at component #7 max dy increased to -2.10081e+28 at component #7 max dy increased to -2.4962e+28 at component #7 max dy increased to -2.96602e+28 at component #7 max dy increased to -3.52729e+28 at component #7 max dy increased to -4.19687e+28 at component #7 max dy increased to -4.99824e+28 at component #7 max dy increased to -5.95495e+28 at component #7 max dy increased to -7.10367e+28 at component #7 max dy increased to -8.48394e+28 at component #7 max dy increased to -1.01411e+29 at component #7 max dy increased to -1.21408e+29 at component #7 max dy increased to -1.45518e+29 at component #7 max dy increased to -1.74677e+29 at component #7 max dy increased to -2.10016e+29 at component #7 max dy increased to -2.52833e+29 at component #7 max dy increased to -3.05063e+29 at component #7 max dy increased to -3.68578e+29 at component #7 max dy increased to -4.46245e+29 at component #7 max dy increased to -5.4126e+29 at component #7 max dy increased to -6.57881e+29 at component #7 max dy increased to -8.01162e+29 at component #7 max dy increased to -9.78251e+29 at component #7 max dy increased to -1.19656e+30 at component #7 max dy increased to -1.46725e+30 at component #7 max dy increased to -1.80408e+30 at component #7 max dy increased to -2.22272e+30 at component #7 max dy increased to -2.74642e+30 at component #7 max dy increased to -3.40196e+30 at component #7 max dy increased to -4.22666e+30 at component #7 max dy increased to -5.26523e+30 at component #7 max dy increased to -6.58025e+30 at component #7 max dy increased to -8.24736e+30 at component #7 max dy increased to -1.03717e+31 at component #7 max dy increased to -1.3088e+31 at component #7 max dy increased to -1.65747e+31 at component #7 max dy increased to -2.10686e+31 at component #7 max dy increased to -2.68927e+31 at component #7 max dy increased to -3.44745e+31 at component #7 max dy increased to -4.43946e+31 at component #7 max dy increased to -5.74656e+31 at component #7 max dy increased to -7.47928e+31 at component #7 max dy increased to -7.14289e+46 at component #7 max dy increased to -4.43626e+54 at component #7 max dy increased to -5.19348e+54 at component #7 max dy increased to 3.19614e+63 at component #7 max dy increased to 3.56268e+75 at component #7 max dy increased to -1.02917e+79 at component #7 max dy increased to -9.74423e+83 at component #7 max dy increased to -6.39195e+86 at component #7 max dy increased to -1.08617e+96 at component #7 max dy increased to -2.50327e+139 at component #7 max dy increased to -2.5451e+147 at component #7 max dy increased to -1.82551e+148 at component #7 max dy increased to -1.82551e+148 at component #7 max dy increased to -4.68542e+152 at component #7 max dy increased to 6.2098e+163 at component #7 max dy increased to 2.86595e+179 at component #7 max dy increased to 7.25633e+239 at component #7 max dy increased to -Inf at component #2
Error using bvp4c (line 196)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in solution>proj (line 17)
sol= bvp4c(@projfun,@projbc,solinit,options);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
function sol= proj
myLegend1 = {};
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.2 0.4 0.6];
for i =1:numel(rr)
a= rr(i);
s=1;h=1;
y0 = [1,0,0,1,1,0,0];
options =bvpset('stats','on','RelTol',1e-4);
m = linspace(0,5);
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))];
end
figure(1)
legend(myLegend1)
hold on
function dy = projfun(x,y)
global max_dy
if isempty(max_dy); max_dy = 0; end
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-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+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)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv);
dy(5)=dt;
dy(6)=ddt;
dy(7)=(1/(a*(a5*s^2*(x+2)+a5*h^2*(x+2)-a8*(x+2)^3)*t+a5*s^2*(x+2)^2+a5*h^2*(x+2)-a8*(x+2)^3))*(-(2*a7*(x+2)-2*a*a11*(x+2))*((1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+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)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv))+(s^2+h^2+2*a5*t-a6*(x+2)^2)*ddt-(a*(s^2+h^2+4*a5-a6*(x+2)^2))*t*ddt-(a*(s^2+h^2+a5+a5))*dt*dt-a*(x+2)*(a5*s^2+a5*h^2+a10*(x+2)^2)*dt*ddt-(2*a7*(x+2)-2*a*a11*(x+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)))+2*a*a7*(x+2)*t*((1/(a*(s^2+a1*h^2-(x+1+1)^2)*t-(s^2+a1*h^2-(x+1+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)))+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dv))+(2*a*a7*(x+2))*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))));
[~, pos] = max(abs(dy));
T = dy(pos(1));
if abs(T) > abs(max_dy)
max_dy = T;
fprintf('max dy increased to %g at component #%d\n', max_dy, pos(1));
end
end
function res = projbc(ya,yb)
res = [ya(1);
ya(3);
ya(5)-1;
ya(7);
yb(1);
yb(3);
yb(5);
];
end
end

3 个评论
显示 1更早的评论 隐藏 1更早的评论

T K
T K 2026-2-17
编辑：T K 2026-2-17
@ Walter Roberson. Thank you for help me. This error is eliminated by the disappearance of terms for the component that causes infinity in the description of the latter equation?
T K
T K 2026-2-17
What do you mean component #7 ?
Walter Roberson
Walter Roberson 2026-2-17
Component #7 refers to dy(7)
Component #2 refers to dy(2)
"This error is eliminated by the disappearance of terms for the component that causes infinity in the description of the latter equation?"
No, not likely. If you remove those components, then you will be working with a model that does not match the equations.
You need to verify that the implementing equations are correct.
You could consider initially implementing the equations using the Symbolic Toolbox, and then converting the equations to numeric form using odeFunction

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Language Fundamentals 的更多信息

标签

提问:

T K
T K
2026-2-17

评论：

Walter Roberson
Walter Roberson
2026-2-17

Community Treasure Hunt

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

Start Hunting!

Translated by