1面向方程的数值积分方法仿真

实验１面向方程的数值积分方法仿真1.实验目的运用CSS01.C仿真程序解题,培养阅读及修改仿真程序的能力,学习并了解仿真程序的结构及特点.通过实验,加深理解4阶龙格-库塔法的原理及其稳定域.2.实验设备：装有BC语言的PC机一台3.实验内容修改CSS01.C仿真程序,对如下系统进行仿真.(1)线性定常系统321xxx=00600120000120321xxx+60000uy=001Txxx321)0()0()0(321xxx=000,u=1(t)a)、对龙格库塔法进行分析：它是一种数值积分法，也就是微分方程初值问题数值计算法，是对初值微分方程的离散化求解。对于数值积分法我们常用的是欧拉法以及二阶和四阶龙格库塔法其原理分别如下：对于形如dtdy=),(ytf的微分方程用欧拉法仿真其迭代公式为),(*1kkkkytfhyy，其中h为仿真步长（下同）；二阶龙格库塔法其迭代公式为yk+1=yk+h/2*(k1+k2)，其中),(1kkytfk，);*,(21hkyhtfkkk四阶龙格库塔法其迭代公式为：)432221(*6/1kkkkhyykk，其中),(1kkytfk，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]Y[1]=k[1][1]+k[1][2]；dtdy=),(ytf——g[1]，（计算k1），g[1]*h/2——k[1][2]，（把k1*h/2存储到k[1][2]）第二次计算：Y[1]=k[1][1]+k[1][3]；dtdy=),(ytf——g[1]，（计算k2），g[1]*h/2——k[1][3]，（把k2*h/2存储到k[1][3]）第三次计算：dtdy=),(ytf——g[1]，（计算k3），Y[1]=k[1][1]+k[1][4]；k[1][4]=g[1]*h；（把k3*h存储到k[1][4]）第四次计算：dtdy=),(ytf——g[1]，（计算k4），此时把k1,k2,k3,k4的值代入下式：)432221(*6/1kkkkhyykk即为6/)4*3**22**21*(1khkhkhkhyykk，由于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]所以可得如程序中的计算公式：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）。b)、计算)(ty理论表达式，首先引入状态变量:yx1dtdyxx122223dtydxx则有ydtdyx202222320020xdtdydtydxuydtdydtyddtyd600600400022022336004000220600)(23ssssG)6125.61)(05.01)(005.01(1)15125.0)(9832.19)(0165.200(600600400220600)(23ssssssssss则)6125.61)(05.01)(005.01(1)(sssssY利用分式部分法，取拉氏反变换有sAsAsAsA6125.6105.01005.0143211)(01sssYA70165.200210203.4)005.01()(ssYA49832.1931024.4)05.01()(sssYA62158.6)6125.61)((15125.04sssYAttteeety0165.20079832.19415125.010204.41024.462158.61)(0tc)、程序实现：在CSS01.C仿真程序基础上,增加floatu；int1u将]6][102[w再修改为]6][801[w在100323(}10033({nnwhilennifdo中的将100改为800，再在);101;1(iiifor中将101改为801最后在输入程序块)(modsub程序改为；;*600]1[*600]3[];3[]2[*200]2[];2[]1[*20]1[];1[]22[;]21[uygyygyygyyuy在输出转换)(repsub程序块中程序改为；];21[]1[ya];22[]2[ya];2[]3[ya];3[]4[ya以下确定仿真步长：由上面方程可得系统传递函数为600)20)(200(600)(ssssG其开环惯量为201s，2001s理论上仿真步长h2|1|a(a为惯量极点)，则有h2|2001|=0.01，可取h为0.01。由于程序中存储的仿真点数最多为800个，所以可取总时间为8。通过作以上的修改，等程序运行后分别输入T1=6.6,T2=0.01N1=3N2=0N3=661J8=1Inputinitialvaluesofstatevariables=0回车0回车0回车实验结果：由仿真结果可以看出：该模型是一个系统阶跃响应实例。当t趋于无穷时y(t)的理论值趋于1若取t为k*0.01(k=1,2,3,4,5,6…)计算y(t)的值并与仿真结果中的y[1]进行比较，发现其值大致相等，所以该仿真是正确的，合理的。(程序见附录xianxin)实验小结：再用根轨迹法分析：1）、确定实轴上的根轨迹，实轴上[0，—200.0165]区域必为根轨迹2）、确定渐进线（根轨迹）3mn故有三条轨迹渐进线333.73a60)12(mnka，180，300（k分别为0、1、２）３）、求分离点：由011niipd有0200120101ddd既转化为0400044032dd7373865722.9112ds92928.1363s（舍去）设分离点处开环增益为k把7373865722.9112ds代入0400022023ksss中得13969.19013k当13969.19013k时临界阻尼当13969.19013k时为欠阻尼4）、确定根轨迹与虚轴交点：本题闭环特征方程式为：04000220*3ksss对上式应用劳斯判据，有*0*1*230220400022022040001kskskss令劳斯表中1s行首项为0得880000*k根据2s行的系数，得如下辅助方程：0220*2ks代入880000*k并令jws解出交点坐标2455532.63w。所以当超过这两点系统进入不稳定状态。（2）实验内容----非线性系统（弱肉强食模型）：)()()()()()()()(tytbxtsydttdytytaxtrxdttdx(1)其中:001.0r，6102a，12000)0(x，01.0s，6101b，600)0(y。a)、理论分析：设兔子数量为)(tx,老虎的数量为)(ty,将它们置于一个特定的环境里,)(tx将以自然增长率r增长,即rxx.但是老虎以兔子为食,致使兔子的增长率降低,设降低程度与老虎数量)(ty成正比,于是相对增长率为ayrrx.常数a反映了老虎捕食兔子的能力.如果没有兔子,老虎就无法生存,设老虎的自然死亡率为s,则syy.兔子为老虎提供了食物,致使老虎死亡率降低,或者说为老虎提供了增长的条件.设增长率与兔子的数量)(tx成正比,于是老虎的相对增长率bxsry.常数b反映了兔子对老虎的供养能力.所以形成上面的模型.下面对上述非线性微分方程组做稳定性分析.首先求平衡点,令,0,0yx得(0,0)和),(arbs.若只考虑),(arbs的稳定性,则可作变量代换，使原点为平稳点,令)()(),()(tartytbstx代入(1)式得)()(arbdtdbsadtd设bs,ar,去掉二次项使之线性化,得abrdtdbasdtd即00abrbasdtddtd(2)可以求出上式系数行列式的特征根为rsi.故线性方程组的通解为)sin()()cos()(00rsitrsit(3)由上式可以看出:(a),的线性化解,既不增长也不衰减,而是连续振动.这意味着平稳是亚稳定的,这是一种广义稳定(Kolmogorov意义下);在平衡点领域显示出稳定的有界循环.(b)令2rsT,则振动周期T=rs2(c)兔子和老虎的总数的振动相位差为090,见图兔子和老虎数量变化图示(4)从解中消去时间t,得12020这条轨线在总体相空间的图形如图轨线在总体相空间的图形（横虚线表示r）(5)可定出在相空间运动的方向,将)()(),()(tartytbstx代入线性化方程可得)()(bsxabrdtdyarybasdtdx因此轨线在总体相空间的图形四个区域中有I:;0,0yxII:;0,0yxIII:0,0yx;IV:0,0yx.从对非线性微分方程组的直接讨论也可知道轨线是一族以平衡点为中心越来越扩展的封闭曲线.封闭曲线对应着(1)式的,周期解,记周期为T,yx,分别为在一个周期内x(t)和y(t)的平均值,则bsdtyysbTdttxTxT)(11)(10其中y(t)为周期函数,y(T)=y(0).同理,ary.对于周期性变化的x(t)和y(t),用它们的平均值yx,表示其大小.上面的分析表明,兔子的平均数量取决于老虎方程)(bxsyy中的参数s和b,而老虎的平均数量取决于兔子方程)(ayrxx中的参数r和a．当兔子的自然增长率r下降时，老虎的数量将减少，而当老虎捕食兔子的能力提高时，对兔子没有影响，只使老虎减少．另一方面，老虎死亡率s上升将导致兔子增多，而兔子对老虎的供养能力b的提高，会导致兔子减少。b)、程序实现：同样在CSS01.C仿真程序基础上作如下修改，先将100323(}10033({nnwhilennifdo中的100改为800，再在);101;1(iiifor中将101改为801最后在系统模型输入程序块)(modsub中程序改为；];2[*]1[*000001.0]2[*01.0]2[];2[*]1[*000002.0]1[*001.0]1[yyygyyyg在输入转换)(repsub程序块中程序改为；];1[]1[ya];2[]2[ya];1[]3[ga];2[]4[ga因为在这个封闭的生态系统里，循环周期较长，故可取仿真总时间为６４００天，步长为１０天，即T1=6400,T2=10N1=2N2=0N3=640J8=1Inputinitialvaluesofstatevariables=12000回车600回车0回车作图程序：void()blot{double;1;0ji);\\\\:,mod,&(&bgibcdeggdrive