matlab 新方法数值求解非光滑ode时遇到的问题
 

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

[/CODE]