1面向方程的数值积分方法仿真
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验１面向方程的数值积分方法仿真
1.实验目的
运用CSS01.C仿真程序解题,培养阅读及修改仿真程序的能力,学习并了解仿真程序的结构及特点.通过实验,加深理解4阶龙格-库塔法的原理及其稳定域.
2.实验设备：装有BC语言的PC机一台
3.实验内容
修改CSS01.C仿真程序,对如下系统进行仿真.
(1)线性定常系统
= + u
y=
= , u=1(t)
a)、对龙格库塔法进行分析：它是一种数值积分法，也就是微分方程初值问题数值计算法，是对初值微分方程的离散化求解。对于数值积分法我们常用的是欧拉法以及二阶和四阶龙格库塔法其原理分别如下：
对于形如 = 的微分方程
用欧拉法仿真其迭代公式为 ，其中h为仿真步长（下同）；
二阶龙格库塔法其迭代公式为yk+1= yk+h/2*(k1+k2)，其中 ，
Y[1]=k[1][1]+k[1][2]；
= ——g[1]，（计算k1），
g[1]*h/2——k[1][2]，（把k1*h/2存储到k[1][2]）
第二次计算：
Y[1]=k[1][1]+k[1][3]；
= ——g[1]，（计算k2），
g[1]*h/2——k[1][3]，（把k2*h/2存储到k[1][3]）
下面对上述非线性微分方程组做稳定性分析.
首先求平衡点,令 得(0,0)和 .若只考虑 的稳定性,则可作变量代换，使原点为平稳点,令 代入(1)式得 设 , ,去掉二次项使之线性化,得 即 (2)
可以求出上式系数行列式的特征根为 .故线性方程组的通解为 (3)
由上式可以看出:
(a) 的线性化解,既不增长也不衰减,而是连续振动.这意味着平稳是亚稳定的,这是一种广义稳定( 意义下);在平衡点领域显示出稳定的有界循环.
所以可得如程序中的计算公式：
Y[1]= k[1][1]+(2*k[1][2]+4* k[1][3]+2* k[1][4]+h*g[1]）；
可见该仿真程序为四阶龙格库塔法的数值积分仿真程序。
以上为只是一个微分方程的情况，若是微分方程组其过程也大概如此。只是将数组y、g和k的其它部分也相应赋值或做运算即可。其实是把g[1],y[1]分别用g[i],y[i]替代；k[1][j]用k[i][j]替代即可（j=1、2、3、4）。
从对非线性微分方程组的直接讨论也可知道轨线是一族以平衡点为中心越来越扩展的封闭曲线.封闭曲线对应着(1)式的,周期解,记周期为T, 分别为在一个周期内x(t)和y(t)的平均值,则 其中y(t)为周期函数,y(T)=y(0).
对上式应用劳斯判据，有
令劳斯表中 行首项为0得 根据 行的系数，得如下辅助方程： 代入 并令 解出交点坐标 。
所以当超过这两点系统进入不稳定状态。
（2）实验内容----非线性系统（弱肉强食模型）：
(1)ቤተ መጻሕፍቲ ባይዱ
其中: ， ， ， ， ， 。
a)、理论分析：设兔子数量为 ,老虎的数量为 ,将它们置于一个特定的环境里, 将以自然增长率 增长,即 .但是老虎以兔子为食,致使兔子的增长率降低,设降低程度与老虎数量 成正比,于是相对增长率为 .常数 反映了老虎捕食兔子的能力.如果没有兔子,老虎就无法生存,设老虎的自然死亡率为 ,则 .兔子为老虎提供了食物,致使老虎死亡率降低,或者说为老虎提供了增长的条件.设增长率与兔子的数量 成正比,于是老虎的相对增长率 .常数 反映了兔子对老虎的供养能力.所以形成上面的模型.
通过作以上的修改，等程序运行后分别输入T1=6.6,T2=0.01 N1=3 N2=0
N3=661J8=1Input initial values of state variables=0回车0回车0回车
实验结果：由仿真结果可以看出：该模型是一个系统阶跃响应实例。当t趋于无穷时y(t)的理论值趋于1若取t为k*0.01(k=1,2,3,4,5,6…)计算y(t)的值并与仿真结果中的y[1]进行比较，发现其值大致相等，所以该仿真是正确的，合理的。(程序见附录xianxin)
四阶龙格库塔法其迭代公式为： ，
其中 ，k2=f(tk+h/2,yk+k1*h/2)，k3=f(tk+h/2,yk+k2*h/2)，k4=f(tk+h,yk+k3*h)；
本程序中大致分以下四次计算：
第一次计算：
g[1]*h/2——k[1][2]，y[1]——k[1][1]，把y[1]的值赋给k[1][1]
(b)令 ,则振动周期T=
(c)兔子和老虎的总数的振动相位差为 ,见图
兔子和老虎数量变化图示
(4)从解中消去时间t,得
这条轨线在总体相空间的图形如图
轨线在总体相空间的图形（横虚线表示 ）
(5)可定出在相空间运动的方向,将 代入线性化方程可得
因此轨线在总体相空间的图形四个区域中有I: II: III: ;IV: .
第三次计算：
= ——g[1]，（计算k3），
Y[1]=k[1][1]+k[1][4]；
k[1][4]=g[1]*h；（把k3*h存储到k[1][4]）
第四次计算：
= ——g[1]，（计算k4），
此时把k1,k2,k3,k4的值代入下式：
即为 ，
由于h*k1=2*k[1][2]，2*h*k2=4*k[1][3]，2*h*k3=2*k[1][4]，h*k4=h*g[1]
实验小结：再用根轨迹法分析：
1）、确定实轴上的根轨迹，实轴上[0，—200.0165]区域必为根轨迹
2）、确定渐进线（根轨迹） 故有三条轨迹渐进线 ， ， （k分别为0、1、２）
３）、求分离点：由 有 既转化为
（舍去）设分离点处开环增益为 把 代入 中得 当 时临界阻尼 当 时为欠阻尼
4）、确定根轨迹与虚轴交点：本题闭环特征方程式为：
b)、计算 理论表达式，首先引入状态变量:
则有
则 利用分式部分法，取拉氏反变换有
c)、程序实现：在CSS01.C仿真程序基础上,增加 u； 将 再修改为 在 中的将100改为800，再在 中将101改为801最后在输入程序块 程序改为；
在输出转换 程序块中程序改为；
以下确定仿真步长：由上面方程可得系统传递函数为 其开环惯量为 ， 理论上仿真步长h<2 (a为惯量极点)，则有h<2 =0.01，可取h为0.01。由于程序中存储的仿真点数最多为800个，所以可取总时间为8。
1.实验目的
运用CSS01.C仿真程序解题,培养阅读及修改仿真程序的能力,学习并了解仿真程序的结构及特点.通过实验,加深理解4阶龙格-库塔法的原理及其稳定域.
2.实验设备：装有BC语言的PC机一台
3.实验内容
修改CSS01.C仿真程序,对如下系统进行仿真.
(1)线性定常系统
= + u
y=
= , u=1(t)
a)、对龙格库塔法进行分析：它是一种数值积分法，也就是微分方程初值问题数值计算法，是对初值微分方程的离散化求解。对于数值积分法我们常用的是欧拉法以及二阶和四阶龙格库塔法其原理分别如下：
对于形如 = 的微分方程
用欧拉法仿真其迭代公式为 ，其中h为仿真步长（下同）；
二阶龙格库塔法其迭代公式为yk+1= yk+h/2*(k1+k2)，其中 ，
Y[1]=k[1][1]+k[1][2]；
= ——g[1]，（计算k1），
g[1]*h/2——k[1][2]，（把k1*h/2存储到k[1][2]）
第二次计算：
Y[1]=k[1][1]+k[1][3]；
= ——g[1]，（计算k2），
g[1]*h/2——k[1][3]，（把k2*h/2存储到k[1][3]）
下面对上述非线性微分方程组做稳定性分析.
首先求平衡点,令 得(0,0)和 .若只考虑 的稳定性,则可作变量代换，使原点为平稳点,令 代入(1)式得 设 , ,去掉二次项使之线性化,得 即 (2)
可以求出上式系数行列式的特征根为 .故线性方程组的通解为 (3)
由上式可以看出:
(a) 的线性化解,既不增长也不衰减,而是连续振动.这意味着平稳是亚稳定的,这是一种广义稳定( 意义下);在平衡点领域显示出稳定的有界循环.
所以可得如程序中的计算公式：
Y[1]= k[1][1]+(2*k[1][2]+4* k[1][3]+2* k[1][4]+h*g[1]）；
可见该仿真程序为四阶龙格库塔法的数值积分仿真程序。
以上为只是一个微分方程的情况，若是微分方程组其过程也大概如此。只是将数组y、g和k的其它部分也相应赋值或做运算即可。其实是把g[1],y[1]分别用g[i],y[i]替代；k[1][j]用k[i][j]替代即可（j=1、2、3、4）。
从对非线性微分方程组的直接讨论也可知道轨线是一族以平衡点为中心越来越扩展的封闭曲线.封闭曲线对应着(1)式的,周期解,记周期为T, 分别为在一个周期内x(t)和y(t)的平均值,则 其中y(t)为周期函数,y(T)=y(0).
对上式应用劳斯判据，有
令劳斯表中 行首项为0得 根据 行的系数，得如下辅助方程： 代入 并令 解出交点坐标 。
所以当超过这两点系统进入不稳定状态。
（2）实验内容----非线性系统（弱肉强食模型）：
(1)ቤተ መጻሕፍቲ ባይዱ
其中: ， ， ， ， ， 。
a)、理论分析：设兔子数量为 ,老虎的数量为 ,将它们置于一个特定的环境里, 将以自然增长率 增长,即 .但是老虎以兔子为食,致使兔子的增长率降低,设降低程度与老虎数量 成正比,于是相对增长率为 .常数 反映了老虎捕食兔子的能力.如果没有兔子,老虎就无法生存,设老虎的自然死亡率为 ,则 .兔子为老虎提供了食物,致使老虎死亡率降低,或者说为老虎提供了增长的条件.设增长率与兔子的数量 成正比,于是老虎的相对增长率 .常数 反映了兔子对老虎的供养能力.所以形成上面的模型.
通过作以上的修改，等程序运行后分别输入T1=6.6,T2=0.01 N1=3 N2=0
N3=661J8=1Input initial values of state variables=0回车0回车0回车
实验结果：由仿真结果可以看出：该模型是一个系统阶跃响应实例。当t趋于无穷时y(t)的理论值趋于1若取t为k*0.01(k=1,2,3,4,5,6…)计算y(t)的值并与仿真结果中的y[1]进行比较，发现其值大致相等，所以该仿真是正确的，合理的。(程序见附录xianxin)
四阶龙格库塔法其迭代公式为： ，
其中 ，k2=f(tk+h/2,yk+k1*h/2)，k3=f(tk+h/2,yk+k2*h/2)，k4=f(tk+h,yk+k3*h)；
本程序中大致分以下四次计算：
第一次计算：
g[1]*h/2——k[1][2]，y[1]——k[1][1]，把y[1]的值赋给k[1][1]
(b)令 ,则振动周期T=
(c)兔子和老虎的总数的振动相位差为 ,见图
兔子和老虎数量变化图示
(4)从解中消去时间t,得
这条轨线在总体相空间的图形如图
轨线在总体相空间的图形（横虚线表示 ）
(5)可定出在相空间运动的方向,将 代入线性化方程可得
因此轨线在总体相空间的图形四个区域中有I: II: III: ;IV: .
第三次计算：
= ——g[1]，（计算k3），
Y[1]=k[1][1]+k[1][4]；
k[1][4]=g[1]*h；（把k3*h存储到k[1][4]）
第四次计算：
= ——g[1]，（计算k4），
此时把k1,k2,k3,k4的值代入下式：
即为 ，
由于h*k1=2*k[1][2]，2*h*k2=4*k[1][3]，2*h*k3=2*k[1][4]，h*k4=h*g[1]
实验小结：再用根轨迹法分析：
1）、确定实轴上的根轨迹，实轴上[0，—200.0165]区域必为根轨迹
2）、确定渐进线（根轨迹） 故有三条轨迹渐进线 ， ， （k分别为0、1、２）
３）、求分离点：由 有 既转化为
（舍去）设分离点处开环增益为 把 代入 中得 当 时临界阻尼 当 时为欠阻尼
4）、确定根轨迹与虚轴交点：本题闭环特征方程式为：
b)、计算 理论表达式，首先引入状态变量:
则有
则 利用分式部分法，取拉氏反变换有
c)、程序实现：在CSS01.C仿真程序基础上,增加 u； 将 再修改为 在 中的将100改为800，再在 中将101改为801最后在输入程序块 程序改为；
在输出转换 程序块中程序改为；
以下确定仿真步长：由上面方程可得系统传递函数为 其开环惯量为 ， 理论上仿真步长h<2 (a为惯量极点)，则有h<2 =0.01，可取h为0.01。由于程序中存储的仿真点数最多为800个，所以可取总时间为8。