ode45在规定区间内提前停止运算！！！

飞翼

function roll_numerical_prediction

y0=[0.1;0.1;0.1;0.1;0];

V0=0.3;

k1=1;k2=1;k3=0.6;

u1=0.09;u2=0.2;u3=0.09;

P1=0.3;P2=0.5;P3=0.5;P4=1;

a1=0.15;a2=0.09;a3=0.07;

options=odeset;options.RelTol=1e-5;options.AbsTol=1e-6;

w=1.2;

T=2*pi/w;

h=T/100;

t0=0;tf=h;

for j=1:2072

   j

    t0=t0+(j-1)*h;

    tf=tf+(j-1)*h;

    if (y0(4)-V0)<0||((-k3*y0(1)-2*u3*y0(2)-k2*y0(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*y0(4)-3*V0*a3*P4*y0(4)^2-a3*P4*y0(4)^3+P3-a2*V0*P4+a3*P4*V0^3+a1*P4)<0&&(-k3*y0(1)-2*u3*y0(2)-k2*y0(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*y0(4)-3*V0*a3*P4*y0(4)^2-a3*P4*y0(4)^3+P3-a2*V0*P4+a3*P4*V0^3-a1*P4)<0&&(y0(4)-V0)==0)

      

        [t,dx]=ode45(@sb1,[t0:h/2:tf],y0,options);

       T(j)=t(end);

       X1(j)=dx(end,1);

       X2(j)=dx(end,2);

       X3(j)=dx(end,3);

       X4(j)=dx(end,4);

       X5(j)=dx(end,5);

       y0=[X1(j);X2(j);X3(j);X4(j);X5(j)];

           if (X4(j)-V0)==0

              y0=[X1(j);X2(j);X3(j);X4(j);X5(j)];

           end

           if (X4(j)-V0)>0&&(X4(j-1)-V0)<0

              a=(V0-X4(j-1))/(X4(j)-V0);

              X1(j)=(X1(j-1)+a*X1(j))/(1+a);

              X2(j)=(X2(j-1)+a*X2(j))/(1+a);

              X3(j)=(X3(j-1)+a*X3(j))/(1+a);

              X4(j)=(X4(j-1)+a*X4(j))/(1+a);

              y0=[X1(j);X2(j);X3(j);X4(j);X5(j)];

           end

           continue;

    end

    if (y0(4)-V0)>0||((-k3*y0(1)-2*u3*y0(2)-k2*y0(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*y0(4)-3*V0*a3*P4*y0(4)^2-a3*P4*y0(4)^3+P3-a2*V0*P4+a3*P4*V0^3+a1*P4)>0&&(-k3*y0(1)-2*u3*y0(2)-k2*y0(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*y0(4)-3*V0*a3*P4*y0(4)^2-a3*P4*y0(4)^3+P3-a2*V0*P4+a3*P4*V0^3-a1*P4)>0&&(y0(4)-V0)==0)

       [t,dx]=ode45(@sb2,[t0:h/2:tf],y0,options);

       T(j)=t(end);

       X1(j)=dx(end,1);

       X2(j)=dx(end,2);

       X3(j)=dx(end,3);

       X4(j)=dx(end,4);

       X5(j)=dx(end,5);

       y0=[X1(j);X2(j);X3(j);X4(j);X5(j)];

       if (X4(j)-V0)==0

              y0=[X1(j);X2(j);X3(j);X4(j);X5(j)];

           end

       if (X4(j)-V0)<0&&(X4(j-1)-V0)>0

          a=(V0-X4(j-1))/(X4(j)-V0);

          X1(j)=(X1(j-1)+a*X1(j))/(1+a);

          X2(j)=(X2(j-1)+a*X2(j))/(1+a);

          X3(j)=(X3(j-1)+a*X3(j))/(1+a);

          X4(j)=(X4(j-1)+a*X4(j))/(1+a);

          y0=[X1(j);X2(j);X3(j);X4(j);X5(j)];

       end

      continue;

    end

    if (y0(4)-V0)==0&&((-k3*y0(1)-2*u3*y0(2)-k2*y0(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*y0(4)-3*V0*a3*P4*y0(4)^2-a3*P4*y0(4)^3+P3-a2*V0*P4+a3*P4*V0^3+a1*P4)>=0&&(-k3*y0(1)-2*u3*y0(2)-k2*y0(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*y0(4)-3*V0*a3*P4*y0(4)^2-a3*P4*y0(4)^3+P3-a2*V0*P4+a3*P4*V0^3-a1*P4)<=0)

       [t,dx]=ode45(@sb,[t0:h/2:tf],y0,options);

       T(j)=t(end);

       X1(j)=dx(end,1);

       X2(j)=dx(end,2);

       X3(j)=dx(end,3);

       X4(j)=dx(end,4);

       X5(j)=dx(end,5);

       y0=[X1(j);X2(j);X3(j);X4(j);X5(j)];

       continue;

    end

end

figure(2)

plot(X3(1:end),X4(1:end),'k-')

X1(end-5:end)

X2(end-5:end)

X3(end-5:end)

X4(end-5:end)

X5(end-5:end)

function dx=sb1(t,x)

k1=1;k2=1;k3=0.6;

u1=0.09;u2=0.2;u3=0.09;

P1=0.3;P2=0.5;P3=0.5;P4=1;

a1=0.15;a2=0.09;a3=0.07;

V0=0.3;w=1.2;

dx=[x(2);-(k1+k3)*x(1)-2*(u1+u3)*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*x(4)-3*V0*a3*P4*x(4)^2-a3*P4*x(4)^3+P2*cos(x(5))+P1+P3-a2*V0*P4+a3*P4*V0^3+a1*P4;x(4);-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*x(4)-3*V0*a3*P4*x(4)^2-a3*P4*x(4)^3+P3-a2*V0*P4+a3*P4*V0^3+a1*P4;w];

function dx=sb2(t,x)

k1=1;k2=1;k3=0.6;

u1=0.09;u2=0.2;u3=0.09;

P1=0.3;P2=0.5;P3=0.5;P4=1;

a1=0.15;a2=0.09;a3=0.07;

V0=0.3;w=1.2;

dx=[x(2);-(k1+k3)*x(1)-2*(u1+u3)*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*x(4)-3*V0*a3*P4*x(4)^2-a3*P4*x(4)^3+P2*cos(x(5))+P1+P3-a2*V0*P4+a3*P4*V0^3-a1*P4;x(4);-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*x(4)-3*V0*a3*P4*x(4)^2-a3*P4*x(4)^3+P3-a2*V0*P4+a3*P4*V0^3-a1*P4;w];

function dx=sb(t,x)

k1=1;k2=1;k3=0.6;

u1=0.09;u2=0.2;u3=0.09;

P1=0.3;P2=0.5;P3=0.5;P4=1;

a1=0.15;a2=0.09;a3=0.07;

V0=0.3;w=1.2;

dx=[1/2*((-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P3-a2*V0*P4+a3*P4*V0^3+a1*P4)*x(2)-(-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P3-a2*V0*P4+a3*P4*V0^3-a1*P4)*x(2))/a1/P4;

    1/2*((-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P3-a2*V0*P4+a3*P4*V0^3+a1*P4)*((-k1-k3)*x(1)-(2*u1+2*u3)*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P2*cos(x(5))+P1+P3-a2*V0*P4+a3*P4*V0^3-a1*P4)-(-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P3-a2*V0*P4+a3*P4*V0^3-a1*P4)*((-k1-k3)*x(1)-(2*u1+2*u3)*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P2*cos(x(5))+P1+P3-a2*V0*P4+a3*P4*V0^3+a1*P4))/a1/P4;

    1/2*((-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P3-a2*V0*P4+a3*P4*V0^3+a1*P4)*V0-(-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P3-a2*V0*P4+a3*P4*V0^3-a1*P4)*V0)/a1/P4;

    0;

    1/2*((-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P3-a2*V0*P4+a3*P4*V0^3+a1*P4)*w-(-k3*x(1)-2*u3*x(2)-k2*x(3)-(2*u2+3*V0^2*a3*P4-a2*P4)*V0-3*V0*a3*P4*V0^2-a3*P4*V0^3+P3-a2*V0*P4+a3*P4*V0^3-a1*P4)*w)/a1/P4];

for语句运行一次代表在步长h上积分一次，j取100为运行一个周期，每次积分用的方程用if 语句进行选择，共三个方程。奇怪的是j的取值从2069开始往后，X1到X4的值将一直保持不变，似乎方程不再参数运算，X1X2X3X4最后6个取值为：

X1=-0.0454   -0.0547   -0.0555   -0.0555   -0.0555  -0.0555

X2=-0.1923   -0.1632   -0.1601   -0.1601   -0.1601  -0.1601

X3=0.1781    0.1932    0.1949    0.1949    0.1949  0.1949

X4=0.2770    0.2979     0.3         0.3          0.3        0.3

plot图上也是画到标注为红色的点后就停止了，希望各位帮帮忙

   

meiyongyuandeze

我看红色数据之后的数据变化很小，你们可以再循环中设置一个收敛残差，残差可以使用没两步解之间的差。

飞翼

不是数据变化小，而是根本就不变了，因此在figure图上的效果就是画到某个点之后就停止了，因为后面的点都一样，是同一个点

meiyongyuandeze

那还是可以用判定语句强制跳出计算的！

飞翼

谢谢！嗯，基于这种情况确实可以跳出运算，但是我的目的是想让它继续运算下去，而不是跳出。
实际情况应该是状态变量的值能够持续更新，一直到规定的时间结束，这里取j=1:2072是我为了看出从哪里开始停止运算而取的，其实我最初取的是j=1:5000，我就是不知道为啥j取从2069往后，状态变量的值保持不变

meiyongyuandeze

那这个应该不是你编程的问题吧，可能是不是算法的的原因！

飞翼

不知道呢，不过还是谢谢了！