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

1 次查看（过去 30 天）
汉武 沈
汉武 沈 2021-2-2
评论： 汉武 沈 2021-2-2
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
  2 个评论
Walter Roberson
Walter Roberson 2021-2-2
http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F
Do not compare calculated numbers for equality.
Also, even at the best of times, all three of those if/elseif tests are going to fail. I suspect that you are looking to test whether the values are in particular ranges, which requires a range test rather than an equality test.
if value <= first_boundary
do first range stuff
elseif value <= second_boundary
do second boundary stuff
elseif value <= third_boundary
do third boundary stuff
else
do stuff when you are further than any bound
end
汉武 沈
汉武 沈 2021-2-2
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by