[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])