【求助】遇到一个十分奇怪的毛病关于 elseif

[Jozgoo]

为了实现在t 在 [400,414] [900,910] 区间
常数Q 为0
t在其余情况 Q都为5.3*10^(-6)

Q设初始值也设为 5.3*10^(-6);
t=0:0.1:2000

写了下面一段程序

if((t >= 400) && (t <= 414))
Q = 0;
elseif ((t >=900) && (t <= 910))
Q = 0;
else
Q = 5.3*10^(-6);
end;

结果中间elseif语句完全没有效果， 在t= [900,910] 区间中 Q仍然=5.3*10^(-6)
而 在t = [400, 414]区间中 Q是=0的

但是我改到 t= [700, 710] 区间中, Q在这时候却可以变成 0了
这时候 t 在 两个区间[400, 414] [700, 710] 中 Q值都变成了0

if ((t >= 400) && (t <= 414))
Q = 0;
elseif ((t >=700) && (t <= 710))
Q = 0;
else
Q = 5.3*10^(-6);
end;

我还试过下面这些 始终都是同样问题，900-910 就不行，700-710就可以

if((t >= 400) & (t <= 414))
Q = 0;
elseif ((t >=900) & (t <= 910))
Q = 0;
else
Q = 5.3*10^(-6);
end;

if (t >= 400) & (t <= 414)
Q = 0;
elseif (t >=900) & (t <= 910)
Q = 0;
else
Q = 5.3*10^(-6);
end;

if t >= 400 & t <= 414
Q = 0;
elseif t >=900 & t <= 910
Q = 0;
else
Q = 5.3*10^(-6);
end;

[Jozgoo]

为什么 完全相同的程序，900-910 elseif 语句就没作用 换成700-710 却可以了？

这是什么毛病啊，

[fanxing39]

请把调用的命令都发上来,看下. 你现在的这些程序,我运行了下,没有问题.
我调用的命令是:
t=905;
if((t >= 400) && (t <= 414))
Q = 0;
elseif ((t >=900) && (t <= 910))
Q = 0;
else
Q = 5.3*10^(-6);
end

程序输出中 Q=0

[Jozgoo]

程序一共有两个.m 文件，一个是主程序scriptEx5.m 另外一个是调用的 EquationEx5.m
EquationEx5.m 是需要求解的函数方程式部分， if, elseif 语句也在这部分程序里

下面是程序原件，可以直接分别拷贝为scriptEx5.m 和EquationEx5.m 运行
运行主程序scriptEx5, 会有图显示出来，在Q=0的时候，曲线会有突起，表示温度有升高
400-414 部分始终正常，900-910 就没变化， 但是将EquationEx5.m 中900-910 改成 700-710
就可以了，这时候可以观察到400-414 和700-710 两个突起

手动将常数Q 改为0 可以在任何位置观察到温度升高的现象， 但用elseif 语句却只有700-710
部分比较正常

下面是主程序 ScriptEx5.m 中的代码

clear all, close all, clc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Localization of unstable steady state using Q = 0 m^3/s (2e)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Model parameters
Q = 5.3*10^(-6);
V = 95*10^(-6);
CA_0 = 800;
CB_0 = 1200;

MyA = -1;
MyB = -2;
CA_s = 0;
CB_s = 0;
T_s = 300;
UA = 0;
Rho = 1000;
Cp = 4186.3;
T_a = 293.0;

DHR = -553.0*10^3;
A = 17.692;
B = 8500;
AlphaA = 1;
AlphaB = 1;
T_0 = 275.0;
CC_0 = 0;
CD_0 = 0;
CE_0 = 0;
MyC = 0.5;
MyD = 0.5;
MyE = 2;
CC_s = 0;
CD_s = 0;
CE_s = 0;

Parameters = [Q V CA_0 CB_0 MyA MyB CA_s CB_s T_s UA Rho Cp T_a DHR A B AlphaA AlphaB T_0 CC_0 CD_0 CE_0 MyC MyD MyE CC_s CD_s CE_s];


% Initial values of state variables
state0 = [CA_s CB_s T_s CC_s CD_s CE_s];


% Setting up simulation time in seconds
tstart1 = 0;
tend1 = 1500;

t = [tstart1: 0.1: tend1];

% Simulation
[time1,State1] = ode15s(@EquationEx5,t,state0,[],Parameters);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure (1)
plot(time1,State1(:,3),'-k')
legend('Q = 5.3*10^-^6 or 0 m^3/s ')
title('T(t) in a CSTR for finding the unstable steady state','FontSize',14)
xlabel('Time [s]','FontSize',11)
ylabel('T [K]','FontSize',11)
% axis([0 100 0 1400])

[Jozgoo]

主程序ScriptEx5.m 调用的EquationEx5.m (ode15s函数中调用) 如下

function dstate = EquationEx5(t,state,pars)

CA = state(1); %
CB = state(2); %
T = state(3); %
CC = state(4); %
CD = state(5); %
CE = state(6); %

%%%%%%%%%%%%%%%%%%
Q = pars (1);
V = pars(2);
CA_0 = pars(3);
CB_0 = pars(4);

MyA = pars(5);
MyB = pars(6);
CA_s = pars(7);
CB_s = pars(8);
T_s = pars(9);
UA = pars(10);
Rho = pars(11);
Cp = pars(12);
T_a = pars(13);

DHR = pars(14);
A = pars(15);
B = pars(16);
AlphaA = pars(17);
AlphaB = pars(18);
T_0 = pars(19);
CC_0 = pars(20);
CD_0 = pars(21);
CE_0 = pars(22);
MyC = pars(23);
MyD = pars(24);
MyE = pars(25);
CC_s = pars(26);
CD_s = pars(27);
CE_s = pars(28);

%%%%%%%%%%%%%%%%%%


if((t >= 400) && (t <= 414))
Q = 0;
elseif ((t >= 900) && (t <= 910))
Q = 0;
else
Q = 5.3*10^(-6);
end

% Balances

dCAdt = Q/V*(CA_0 - CA) + MyA*(exp(A - B/T)*CA^AlphaA*CB^AlphaB);

dCBdt = Q/V*(CB_0 - CB) + MyB*(exp(A - B/T)*CA^AlphaA*CB^AlphaB);

dTdt = Q/V*(T_0 - T) + UA/(V*Rho*Cp)*(T_a - T) + (-DHR)/(Rho*Cp)*(exp(A - B/T)*CA^AlphaA*CB^AlphaB);

dCCdt = Q/V*(CC_0 - CC) + MyC*(exp(A - B/T)*CA^AlphaA*CB^AlphaB);

dCDdt = Q/V*(CD_0 - CD) + MyD*(exp(A - B/T)*CA^AlphaA*CB^AlphaB);

dCEdt = Q/V*(CE_0 - CE) + MyE*(exp(A - B/T)*CA^AlphaA*CB^AlphaB);


dstate = [dCAdt dCBdt dTdt dCCdt dCDdt dCEdt]';

[Jozgoo]

发现是什么毛病了，

因为ode15s 这个solver不精确，

根本没有运行 t在 900-910 之间的数据，直接跳过去了

汗啊

ode45 精度会相对高些， 但其实也同样可能出现类似问题，
虽然设了 t = [0: 0.1: 2000]; 但实际上t值跨度可能达到几十甚至上百。

看来matlab 这软件还需要完善啊

[fanxing39]

正在思考中,确实很奇怪哦
这个是刚性微分方程的求解,对于那par 28个参数是个什么东东,程序里好象没有全用
不防把你要求解的微分方程附上来,大家一起看下? 如何

[fanxing39]

哈哈,经过你的提醒,我明白了, 调整下" 相对误差容许上限" 就可以了
原程序中"ode45"命令前面先进行求解参数的修改,使得求解保证足够的精度
插入的命令如下:
options=odeset('RelTol',1e-7)%默认是0.001,现在改为10^(-7)
[time1,State1] = ode45(@EquationEx5,t,state0,options,Parameters);

[Jozgoo]

作者: fanxing39

哈哈,经过你的提醒,我明白了, 调整下" 相对误差容许上限" 就可以了
原程序中"ode45"命令前面先进行求解参数的修改,使得求解保证足够的精度
插入的命令如下:
options=odeset('RelTol',1e-7)%默认是0.0...

hehe 多谢阿 以前没注意这个问题 t=[0:0.1:2000] 本来以为可以设置 取t值间隔0.1
结果根本没效果

看了一下matlab的帮助文件，这个RelTol 还有 AbsTol 指的
似乎都是相对于函数值来说的

在函数值相对平缓的区域， 计算的过程中，自变量取值跨度会相当的大，不知道有什么
办法才可以直接控制自变量的精度

看帮助文件里面写可以直接设置 t = [t1,t2,...,tf]， 在中间直接指定要计算的自变量值
但我试了一下，发现不行 如果我指定 t = [0,401,405,414,901,905,909,1500]
整个函数竟然变成直的了，全乱套了

[fanxing39]

作者: Jozgoo

hehe 多谢阿 以前没注意这个问题 t=[0:0.1:2000] 本来以为可以设置 取t值间隔0.1
结果根本没效果

看了一下matlab的帮助文件，这个RelTol 还有 AbsTol 指的
似乎都是相对于函数值来说的

在函数值相对平缓的区域， 计算的过程中，自变量取值跨度会相...

RelTol 为相对误差上限,默认值为0.001(即0.1%的相对误差),在一些特殊微分方程求解中,为了保证较高的精度,还应当再适当减小该值.
AbsTol 为一个向量,其分量表示每个状态变量容许的绝对误差,其默认值为10^(-6).当然可以自由设置其值,以改变求解精度
MaxStep 为求解方程最大容许的步长
这些都可以由options来更改,方法同上面的那个程序.
微分方程的求解,你上面用的都是数值逼近的方法,每个方法各有优缺点,建议你参阅下相关书籍.
总的来说求解过程要注意三个方面
1.选择适当的步长. 如果步长太大,误差很大,你上面出现的问题主要是由于步长过大导致的.但是如果步长太小,又会产生较大的累积误差.
2. 改进近似算法精度.比较成功的方法有RUNGE-KUTTA法,Adams法等,ode45就是综合了这两种方法.
3.采用变步长的方法. 前面说的"适当"地选择步长,这本身就是个模糊的概念,如果适当地选择步长取决于经验.事实上,很多种方法都容许变步长的求解,如果误差较小时,可自动地增加步长,而误差较大时再自动减小步长,从而精确,有效地求解给出的常微分方程的初值问题.

以上言论摘自 "高等应用数学的MALTAB求解" 一书,作者 薛定宇 陈阳泉
我手上没有这本书的电子版,如果有谁有的话,不妨分享下. 我觉得买一本更好,这本书应该说是matlab方面的经典书籍了.