why do i get a wrong answers when using the if else sentences in the defined-function?

1 view (last 30 days)

Commented: 汉武 沈 on 2 Feb 2021
hi all! I use a new method to solve ode numerically, called QSS, it doesn't matter if you can't understand the method. My problem is that when I vectorized my codes which means I have to defined the ode in the user defined-function, bacause the ode is nonsmooth, I have to use the 'if else' sentence in the defined-function, however the results is not right compared to the original codes(the original codes puts the if else sentence in the loop), can anyone tell me how could that happen? Why dosen't the if else sentence that put in the function work? I really need someone to help me fix the problem. Thank you all:)
Below is my codes.First one is the original codes, second one is the vectorized codes.
Besides, both codes work well when solving the smooth odes.
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','ambient temperature')
% figure
% plot(A(:,1),A(:,5));hold on
% plot(A(:,1),A(:,6));
toc
second one
% 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

Thank you for your comment.I am still confused.The first above codes is the original one, which compare for equilty works fine, the codes is okay. The second one failed, it doesn't compare for equilty. My problem in the second codes is that I put the 'if else' in the defiined-function which will be passed to the main function just like we did when we use the ode suits, however it seemed that the 'if else' sentence didn't work, the variables x(1)、x(2) should fluctuate between 20.5 and 19.5, if you run the first codes, you will get the right answers.Is there something wrong when 'if else' sentences be put in the defined-function? As you can see the first codes also catains the exactly same 'if else' sentences, but it works well.