Pleas help me to run this simple code

T K
2026 2월 17
1 답변
업데이트 시간: 2026 2월 17
조회 수: 13 (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에서 Programming에 대해 자세히 알아보기

태그

질문:

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