%置初值
%置仿真的时间范围
%求微分方程
%绘制第一条状态响应曲线（黑色、实线） %保持在同一个图形窗口绘图
%绘制第二条状态响应曲线（红色、点划线）
%标注图例 %标识X轴的名称 %标识Y轴的名称 %加网格线
%图形标题 子程序：
function xdot=Appl_1_2_func(t, x) %函数头 xdot1(1)=- x(1) *(x(1)*x(1)+x(2)*x(2))+x(2); xdot1(2)=-x(1)-x(2) (x(1)*x(1)+x(2)*x(2)); xdot=xdot1'; 主程序：
clear x0=[10,10];%置初值
tspan=[0,120]; %置仿真的时间范围
[t,x]=ode45('Appl_1_2_func', tspan, x0); %求微分方程
plot(t, x(:,1), 'k', 'LineWidth',2); %绘制第一条状态响应曲线（黑色、实线）
hold %保持在同一个图形窗口绘图
plot(t, x(:,2), 'r-.', 'LineWidth',2); %绘制第二条状态响应曲线（红色、点划线） legend('x(1)', 'x(2)'); %标注图例 xlabel('t'); %标识X轴的名称
ylabel('状态向量'); %标识Y轴的名称 grid; %加网格线
title('系统的状态响应'); %图形标题 3、（本题15分）某地区某病菌传染的系统动力学模型为
?x ?x1(0)?6201??0.001x1x2?
? ??x2(0)?10?x2?0.001x1x2?0.072x2
?
?x?x3(0)?70 3?0.072x2??
式中， x1表示可能受到传染的人数，x2 表示已经被传染的病人数，x3表示已治愈的人数。试用ode45编程，对其进行仿真研究，并在同一个图形窗口用不同的线型（实线、虚线和冒号线）绘制状态响应曲线，仿真时间范围：0~30天。 图形要求：
图形标题：病菌传染模型的状态响应；X轴的名称：“time(天) t0=0, tf=30”， Y轴的名称：“x(人):x1(0)=620,x2(0)=10;x3(0)=70”。
子程序是描述微分方程组的M函数，函数名为“fun2_4”， 函数的输入变量分别为： t, x，输出列向量为：xdot。 主程序：
% 置状态变量初值 % 置仿真时间区间
% 调用ode45求仿真解
% 用不同的线型绘制仿真结果曲线 % 对横轴进行标注 % 对纵轴进行标注
%标注图例，其中第一条第、第二条和第三条分别为x1、x2和x3 %加网格线
%图形标题 子程序：
%函数头
% 第一个微分方程 % 第二个微分方程 % 第三个微分方程 xdot=xdot1'; 主程序： clear
x0=[620,10,70]; % 置状态变量初值 tspan=[0,30]; % 置仿真时间区间 [t,x]=ode45('fun2_4',tspan,x0); % 调用ode45求仿真解
plot(t,x(:,1), t,x(:,2),‘g--’,t,x(:,3),‘r:’); % 用不同的线型绘制仿真结果曲线
xlabel('time(天) t0=0, tf=30'); % 对横轴进行标注 ylabel('x(人):x1(0)=620,x2(0)=10;x3(0)=70'); % 对纵轴进行标注
legend('x1','x2','x3'); %标注图例，其中第一条第、第二条和第三条分别为x1、x2和x3
grid; %加网格线
title('病菌传染模型的状态响应'); %图形标题 子程序：
function xdot=fun2_4(t,x) %函数头
xdot1(1)=-0.001*x(1)*x(2); % 第一个微分方程 xdot1(2)=0.001*x(1)*x(2)-0.072*x(2); % 第二个微分方程 xdot1(3)=0.072*x(2); % 第三个微分方程 xdot=xdot1';
计算题
1、用二阶龙格—库塔法求解方程对数值稳定性的影响。
1y'??y,??0，分析对计算步长h有何限制，说明h
?h?y?y?(k1?k2)k?k?12?1? 解： ?k1??yk
??11?k??(y?ykh)k?2???h2yk?1?yk(1??2)?2? 稳定系统最终渐进收敛。 得到
h 系统稳定则
h21??2?1?2?h 计算得0?h?2?。
h的选取不能超出上述范围，否则系统不稳定。
2、（本题15分）已知y?y?t,y(0)?1，取计算步距h=0.1，试分别用欧拉法、四阶龙格—库塔法求t=h时的y值，并将求得的y值与精确解y(t)?2e?1?t比较，并说明造成差异的原因。 解：（1） 欧拉法：
yn?1?yn?(yn?tn)h
y1?1?(1?0)?0.1?1.1 （5分） （2） 四阶龙格—库塔法： yn?1?yn?t.h(k1?2k2?2k3?k4) 6?k1?yn?tn??k2?yn?hk1?tn?h?22 ?
hh?k3?yn?k2?tn??22?k?y?hk?t?hn3n?4 k1=1，k2=1.1，k3=1.105，k4=1.2105
y1?1.1103 （5分）
y(0.1)=1.1103 （2分）
计算结果产生差异是由于两种方法的精度不一样，RK4方法精度更高。（3分）
3、（本题10分）设Ty(t)?y(t)?k，试分别用欧拉法、二阶龙格—库塔法求y(t)的差分方程，如果步长h大于2T将会产生什么结果？试说明其原因。
欧拉法：
.hkhym?1?(1?)ym? （4分）
TTRK2法：
ym?1hh2khkh2?(1??)ym?? （4分）
T2T2T2T2显然，当h?2T时，数值解将发散。系统的特征值???出稳定性范围。（2分）
1，若h?2T，则?h?2，超T4、（本题15分）已知y??y?t,y(0)?1，，取计算步距h=0.1，试分别用欧拉法、四阶龙格—库塔法求t=h时的y值，并说明造成差异的原因。
解：被求函数y的导函数y?f(t,y)??y?t,y(0)?1，以下分别用两种方法求解
（1） 欧拉法 由欧拉法的递推公式
得：
.2.2 （5分）
（2） 四阶龙格—库塔法 RK4的递推公式为：
其中
由已知条件，
，
，由
递推出
时的值
（5分）
（3）计算结果产生差异是由于两种方法的精度不一样，RK4方法精度更高。（5分）
5、（本题15分）已知微分方程及其初值：
取计算步距h=0.2，试用四阶龙格—库塔法计算y(0.4)的近似值，至少保留四位小数。 解： 此处f (t，y)＝8－3y，四阶龙格――库塔法公式为
其中 ?1＝f (tk，yk)；?2＝f (tk+0.5h，yk+0.5h?1)； ?3＝f (tk+0.5h，yk+0.5h?2)；?4＝f (tk+h，yk+h?3)
其中 ?1＝8－3yk；?2＝5.6－2.1yk； ?3＝6.32－2.37yk；?4＝4.208－1.578yk
＝1.2016＋0.5494yk (k＝0，1，2，…)
当x0＝0，y0＝2，
控制系统数字仿真习题 doc
