[MATLAB数学相关] matlab 新方法数值求解非光滑ode时遇到的问题

[wertshw]

大家好！我在用一种新的数值求解方法求解非光滑ode，这种方法叫QSS(量化状态系统)，原先的程序可以求解，可是为了以后方便，我把这个程序矢量化实现，需要像在用ODE45时一样定义包含ode的自定义函数，因为是非光滑ode，需要用到if else来判断条件，可是矢量化后的程序中，if else 语句似乎没有激活，结果中变量x(1)和x(2)应该在19.5和20.5之间来回波动，可是在这里却是在很小的范围内波动，这会不会是因为if else 语句定义在了自定义函数中了，求求各位老师能够帮我看一下问题出在了哪里，我将不胜感激！
第一个程序是原先没有矢量化的程序，复制、粘贴、运行可以得到正确的结果；
第二个程序是矢量化的程序，得到的结果可以和第一个程序对比一下，就会发现它的结果波动很小。
第一个程序：

clc
clear
tic
%为各变量赋初值
delta_q=1e-4;
q1=21;q2=25;q3=17;
x1=q1;x2=q2;x3=q3;
t=0;delta_t=0;
A=zeros(105609,6);
n=0;
C1=550;C2=600;Cw=580;R1=1.8;R2=2.2;R1w=2;R2w=2;Ta=32;p1=6;p2=6;m1=1;m2=1;T_ref=20;pi=3.1415926;
%开始while循环
while (t<3600*12)
    if(x1>T_ref+0.5)
        m1=1;
    elseif(x1<T_ref-0.5)
        m1=0;
    end
    if(x2>T_ref+0.5)
        m2=1;
    elseif(x2<T_ref-0.5)
        m2=0;
    end
     Ta=10*sin(pi*3.4722e-05*t-pi/2)+20;
     TTa=pi*3.4722e-05*10*cos(pi*3.4722e-05*t-pi/2);
    TTTa=-(pi*3.4722e-05)^2*10*sin(pi*3.4722e-05*t-pi/2);
    Dx1=((Ta-q1)/R1+(q3-q1)/R1w-m1*p1)/C1;  
    Dx2=((Ta-q2)/R2+(q3-q2)/R2w-m2*p2)/C2; 
    Dx3=((q1-q3)/R1w+(q2-q3)/R2w)/Cw;
    DDx1=((TTa-Dx1)/R1+(Dx3-Dx1)/R1w)/C1;
    DDx2=((TTa-Dx2)/R2+(Dx3-Dx2)/R2w)/C2;
    DDx3=((Dx1-Dx3)/R1w+(Dx2-Dx3)/R2w);
   %求Δt1和Δt2的值
       delta_t1=sqrt(2*delta_q/abs(DDx1));
       delta_t2=sqrt(2*delta_q/abs(DDx2));
       delta_t3=sqrt(2*delta_q/abs(DDx3));
       
     delta_tmin=min([delta_t1 delta_t2 delta_t3]);


   %比较Δt1和Δt2的大小，进而确定q1和q2谁先跃迁
   if (delta_t1==delta_tmin)
       delta_t=delta_t1;
       t=t+delta_t;
       x1=x1+Dx1*delta_t+0.5*DDx1*delta_t^2;
       x2=x2+Dx2*delta_t+0.5*DDx2*delta_t^2;
       x3=x3+Dx3*delta_t+0.5*DDx3*delta_t^2;
       q1=x1;
       q2=q2+Dx2*delta_t;
       q3=q3+Dx3*delta_t;
   elseif (delta_t2==delta_tmin)
       delta_t=delta_t2;
       t=t+delta_t;
       x1=x1+Dx1*delta_t+0.5*DDx1*delta_t^2;
       x2=x2+Dx2*delta_t+0.5*DDx2*delta_t^2;
       x3=x3+Dx3*delta_t+0.5*DDx3*delta_t^2;
       q2=x2;
       q1=q1+Dx1*delta_t;
       q3=q3+Dx3*delta_t;
   elseif (delta_t3==delta_tmin)
       delta_t=delta_t3;
       t=t+delta_t;
       x1=x1+Dx1*delta_t+0.5*DDx1*delta_t^2;
       x2=x2+Dx2*delta_t+0.5*DDx2*delta_t^2;
       x3=x3+Dx3*delta_t+0.5*DDx3*delta_t^2;
       q3=x3;
       q1=q1+Dx1*delta_t;
       q2=q2+Dx2*delta_t;
   end
   n=n+1;
   A(n,1)=t;
   A(n,2)=x1;
   A(n,3)=x2;
   A(n,4)=x3;
   A(n,5)=m1;
   A(n,6)=m2;
end
tt=0:3600*12;
yy=10*sin(pi*3.4722e-05*tt-pi/2)+20;
figure
  plot(A(:,1),A(:,2),'--','LineWidth',1.25);hold on
  plot(A(:,1),A(:,3),'-.','LineWidth',1.25);
  plot(A(:,1),A(:,4),'-','LineWidth', 1.25);
  plot(tt,yy,'-k');
  xticks(0:3600:13*3600); %这样x轴会每隔10显示一个刻度
  xticklabels({'6','7','8','9','10','11','12','13','14','15','16','17','18'});%为说明效果，省略了部分内容，写代码时应该与xticks对应写够10+1个label
  xlabel('时刻（h）','FontSize',12,'FontWeight','bold','Color','k');
  ylabel('摄氏度（℃）','FontSize',12,'FontWeight','bold','Color','k');
  grid on
   legend('room1','room2','wall','室外温度')
% figure
%    plot(A(:,1),A(:,5));hold on
%    plot(A(:,1),A(:,6));
toc

第二个程序：

% QSS1矢量化
clc
clear
tic
%为各变量赋初值
delta_q=[1; 1; 1]*1e-4;
q=[21; 25; 17];
x=q;
t=0;
A=zeros(107666,length(q)+1);
delta_x=zeros(length(q),1);
n=0;
%开始while循环
while (t<3600*12)
    
    ff=Dx(t,q);
    aa=ff(1:length(q));
    aaa=ff(length(q)+1:end);
   %求Δt1和Δt2的值
    delta_t=sqrt(2*delta_q./abs(aaa));
   %比较Δt1和Δt2的大小，进而确定q1和q2谁先跃迁
    cc=find(delta_t==min(delta_t));
   %存在多个最小值相等，取第一个最小值
    if length(cc)>1
        cc=cc(1);
    end
    dd=find(delta_t~=min(delta_t));
    t=t+delta_t(cc);
    x=x+aa.*delta_t(cc)+0.5*aaa.*delta_t(cc)^2;
    q(cc)=x(cc);
    q(dd)=q(dd)+aa(dd).*delta_t(cc);
   %赋值 循环
    n=n+1;
    A(n,1)=t;
    A(n,2)=x(1);
    A(n,3)=x(2);
    A(n,4)=x(3);
end
toc

% % 绘图
hold all;
LineWidth=1.25;
tt=0:3600*12;
yy=10*sin(pi*3.4722e-05*tt-pi/2)+20;
plot(A(:,1),A(:,2),'--','LineWidth',LineWidth);
plot(A(:,1),A(:,3),'-.','LineWidth',LineWidth);
plot(A(:,1),A(:,4),'-','LineWidth', LineWidth);
plot(tt,yy,'-k');
xticks(0:3600:13*3600); %这样x轴会每隔10显示一个刻度
xticklabels({'6','7','8','9','10','11','12','13','14','15','16','17','18'});%为说明效果，省略了部分内容，写代码时应该与xticks对应写够10+1个label
  xlabel('时刻（h）','FontSize',12,'FontWeight','bold','Color','k');
  ylabel('摄氏度（℃）','FontSize',12,'FontWeight','bold','Color','k');
grid on; grid minor;
% title({'QSS量子=',string(delta_q(1)) });
% legend
% legend('QSS求解');
%比较误差

function f=Dx(t,x)
    C1=550;C2=600;Cw=580;R1=1.8;R2=2.2;R1w=2;R2w=2;
    p1=6;p2=6;m1=1;m2=1;T_ref=20;pi=3.1415926;
    if(x(1)>T_ref+0.5)
        m1=1;
    elseif(x(1)<T_ref-0.5)
        m1=0;
    end
    if(x(2)>T_ref+0.5)
        m2=1;
    elseif(x(2)<T_ref-0.5)
        m2=0;
    end
     Ta=10*sin(pi*3.4722e-05*t-pi/2)+20;
     TTa=pi*3.4722e-05*10*cos(pi*3.4722e-05*t-pi/2);
    f=zeros(6,1);
    f(1)=((Ta-x(1))/R1+(x(3)-x(1))/R1w-m1*p1)/C1;
    f(2)=((Ta-x(2))/R2+(x(3)-x(2))/R2w-m2*p2)/C2; 
    f(3)=((x(1)-x(3))/R1w+(x(2)-x(3))/R2w)/Cw;
    f(4)=((TTa-f(1))/R1+(f(3)-f(1))/R1w)/C1;
    f(5)=((TTa-f(2))/R2+(f(3)-f(2))/R2w)/C2;
    f(6)=((f(1)-f(3))/R1w+(f(2)-f(3))/R2w);
end