面向方程的数值积分方法仿真(共12页).doc

《面向方程的数值积分方法仿真(共12页).doc》由会员分享，可在线阅读，更多相关《面向方程的数值积分方法仿真(共12页).doc（12页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上实验面向方程的数值积分方法仿真1. 实验目的 运用CSS01.C仿真程序解题,培养阅读及修改仿真程序的能力,学习并了解仿真程序的结构及特点.通过实验,加深理解4阶龙格-库塔法的原理及其稳定域.2. 实验设备：装有BC语言的PC机一台3. 实验内容修改CSS01.C仿真程序,对如下系统进行仿真.(1) 线性定常系统 =+uy=, u=1(t)a)、对龙格库塔法进行分析：它是一种数值积分法，也就是微分方程初值问题数值计算法，是对初值微分方程的离散化求解。对于数值积分法我们常用的是欧拉法以及二阶和四阶龙格库塔法其原理分别如下：对于形如=的微分方程用欧拉法仿真其迭代公式为 ，

2、其中h为仿真步长（下同）；二阶龙格库塔法其迭代公式为 yk+1= yk+h/2*(k1+k2)，其中， 四阶龙格库塔法其迭代公式为：， 其中，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)；本程序中大致分以下四次计算：第一次计算：g1*h/2k12，y1k11，把y1的值赋给k11Y1= k11+ k12；=g1 ， （计算k1），g1*h/2k12，（把k1*h/2存储到k12）第二次计算：Y1= k11+ k13；=g1 ， （计算k2），g1*h/2k13，（把k2*h/2存储到k13）第三次计算：=

3、g1，（计算k3）， Y1= k11+ k14；k14= g1*h；（把k3*h存储到k14）第四次计算：=g1，（计算k4），此时 把k1,k2,k3,k4 的值代入下式： 即为 ，由于h*k1=2*k12，2*h*k2=4*k13，2*h*k3=2*k14，h*k4=h*g1所以可得如程序中的计算公式：Y1= k11+(2*k12+4* k13+2* k14+h*g1）；可见该仿真程序为四阶龙格库塔法的数值积分仿真程序。以上为只是一个微分方程的情况，若是微分方程组其过程也大概如此。只是将数组y、g和k的其它部分也相应赋值或做运算即可。其实是把g1,y1分别用gi,yi替代；k1j用kij替

4、代即可（j=1、2、3、4）。b)、计算理论表达式，首先引入状态变量: 则有 则利用分式部分法，取拉氏反变换有 c)、程序实现：在CSS01.C仿真程序基础上,增加 u； 将再修改为 在 中的将100改为800，再在中将101改为801 最后在输入程序块程序改为；在输出转换程序块中程序改为； 以下确定仿真步长：由上面方程可得系统传递函数为其开环惯量为，理论上仿真步长h2(a为惯量极点)，则有h2=0.01，可取h为0.01。由于程序中存储的仿真点数最多为800个，所以可取总时间为8。通过作以上的修改，等程序运行后分别输入T1=6.6,T2=0.01 N1=3 N2=0N3=661 J8=1 I

5、nput 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)的值并与仿真结果中的y1进行比较，发现其值大致相等，所以该仿真是正确的，合理的。(程序见附录xianxin)实验小结：再用根轨迹法分析：1）、确定实轴上的根轨迹，实轴上0，200.0165区域必为根轨迹2）、确定渐进线（根轨迹） 故有三条轨迹渐进线 ， ， （ k 分别为0、1、 ）、求分离点：由有既转化为（舍去）设分离点处开环增益为

6、把代入中得当时临界阻尼当时为欠阻尼4）、确定根轨迹与虚轴交点：本题闭环特征方程式为：对上式应用劳斯判据，有令劳斯表中行首项为0 得根据行的系数，得如下辅助方程：代入 并令解出交点坐标。 所以当超过这两点系统进入不稳定状态。（2）实验内容-非线性系统（弱肉强食模型）： (1) 其中: ，。 a)、 理论分析： 设兔子数量为,老虎的数量为,将它们置于一个特定的环境里, 将以自然增长率增长,即.但是老虎以兔子为食,致使兔子的增长率降低,设降低程度与老虎数量成正比,于是相对增长率为.常数反映了老虎捕食兔子的能力.如果没有兔子,老虎就无法生存,设老虎的自然死亡率为,则.兔子为老虎提供了食物,致使老虎死亡

7、率降低,或者说为老虎提供了增长的条件.设增长率与兔子的数量成正比,于是老虎的相对增长率.常数反映了兔子对老虎的供养能力.所以形成上面的模型. 下面对上述非线性微分方程组做稳定性分析.首先求平衡点,令得(0,0)和.若只考虑的稳定性,则可作变量代换，使原点为平稳点,令代入(1)式得设,去掉二次项使之线性化,得即 (2)可以求出上式系数行列式的特征根为.故线性方程组的通解为 (3)由上式可以看出:(a)的线性化解,既不增长也不衰减,而是连续振动.这意味着平稳是亚稳定的,这是一种广义稳定(意义下);在平衡点领域显示出稳定的有界循环.(b) 令,则振动周期T=(c) 兔子和老虎的总数的振动相位差为,见

8、图兔子和老虎数量变化图示 (4) 从解中消去时间t,得这条轨线在总体相空间的图形如图轨线在总体相空间的图形 （横虚线表示）(5) 可定出在相空间运动的方向,将代入线性化方程可得因此轨线在总体相空间的图形四个区域中有I:II:III:;IV: .从对非线性微分方程组的直接讨论也可知道轨线是一族以平衡点为中心越来越扩展的封闭曲线.封闭曲线对应着(1)式的,周期解,记周期为T,分别为在一个周期内x(t)和y(t)的平均值,则 其中y(t)为周期函数,y(T)=y(0).同理,.对于周期性变化的x(t)和y(t),用它们的平均值表示其大小.上面的分析表明,兔子的平均数量取决于老虎方程中的参数s和b,而

9、老虎的平均数量取决于兔子方程中的参数r和a当兔子的自然增长率r下降时，老虎的数量将减少，而当老虎捕食兔子的能力提高时，对兔子没有影响，只使老虎减少另一方面，老虎死亡率s上升将导致兔子增多，而兔子对老虎的供养能力b的提高，会导致兔子减少。b)、程序实现：同样在CSS01.C仿真程序基础上作如下修改，先将中的100改为800，再在中将101改为801最后在系统模型输入程序块中程序改为；在输入转换程序块中程序改为； 因为在这个封闭的生态系统里，循环周期较长，故可取仿真总时间为天，步长为天，即T1=6400,T2=10 N1=2 N2=0N3=640 J8=1 Input initial values

10、 of state variables=12000回车 600回车 0回车作图程序： y 实验结果：仿真结果与理论分析图形基本一样如前面两图所示：兔子和老虎的数量各自在趋于一个动态平衡。即在一定的值上下来回振动，这是一种广义的平衡状态。它们之间数量关系近似为一个椭圆 （程序见附录feixian1）4、实验小结：对数值积分的稳定性分析：一个本来是稳定的系统，利用数值积分法进行仿真的时候常常会得出不稳定的结论。造成这种现象的原因往往是计算的字长选取的太大，当步选取的太大的时候，数值积分法会使得各种误差传递出去，从而引起计算不稳定。而误差的来源一般是以下几个：初始误差（在实际计算中给出的初值通常会不

11、太准确）；由于计算机字长限制造成的舍入误差：在某一步长下产生的截断误差。而这些误差都可能在计算中向下传播。数值解的稳定性，是指在拢动（初始误差、舍入、截断误差等效影响下，其计算过程中的累积误差不会随计算步数的增加而无限增长。不同的数值法对应着不同的差分递推公式。一个数值是否稳定取决于该差分方程的特征根是否满足稳定性要求。下面以一阶显式亚当姆斯方法来讨论数值积分的稳定性分析方法假定系统的微分方程为，（）其递推公式为（）假定（k=0,1,2,.）为它的一个仿真解，另外设为其准确解，则有（）（）（）得整理后为其特征方程为为了使误差序列不随k的增加而增加，必须要求其特征根在单位圆内，即，它所对应的区域就是一阶显式亚当姆斯算法的稳定域。因此，如果系统方程，运用一阶显式亚当姆斯算法时，则只有当才能保证计算是稳定的。专心-专注-专业