數(shù)學(xué)模型期末作業(yè)

1、數(shù)學(xué)模型課程名稱 用打靶法求周期解 學(xué) 院 物理學(xué)院 年級 2013 學(xué) 生 姓 名 和 學(xué) 號l 鄭勇 20137022 l 趙茂雄 20136970l 席明 20136952 開 課 時 間 2014 至 2015 學(xué)年第 1 學(xué)期一，問題引入： 二，原理分析：1，打靶法原理：打靶法的實質(zhì)是把微分方程的邊值問題，轉(zhuǎn)化成初值問題求解；求非線性動力系統(tǒng)的周期解就是，求如下形式的邊值問題：T為響應(yīng)的最小周期.。在某一給定的初始值表示初值迭代的次數(shù))下具體計算過程如下:首先根據(jù)方法可得由此可求得若（為允許誤差)則為所求非線性系統(tǒng)的周期解.若不然令,使用如下Newotn迭代關(guān)系求得下次計算需要的初值

2、其中，微分算子在此關(guān)鍵要求出的值。令且方程兩邊對求導(dǎo)可得矩陣的初值問題其中，再次使用方法即可求得的值；代入式(4)即可求得下次計算需要的初值。在對式(5)使用Runge一Kutta方法求步的值時但在第1次的Runge一Kutta方法中,我們并沒有求的式（7）需要的的值，其中的計算通常采用線性插值方法求得，即，或者重新縮短步長計算獲得，但是在計算精度上與Runge一Kutta方法不匹配,且會增加計算量,必然在以后的計算中產(chǎn)生誤差.為此我們只需要在第2次使用Runge-Kutta方法時直接對（2）式兩邊求導(dǎo)即可得第2次使用Runge-Kutta方法所需的值，具體過程如下：式（2）兩邊對求導(dǎo)可得下面

3、求，在式（3）中令則有顯然由于都是第1次使用Runge-Kutta方法求得的，而所以第步的值可表示為顯然使用該種方法計算Runge-Kutta所需要的參數(shù)時候不需要的值，避免了誤差的產(chǎn)生。2，F(xiàn)loquet理論：若系統(tǒng)(2)的特征根的模小于l,則系統(tǒng)(1)存在唯一的w一周期解。三，理論推理：的周期解的公式推導(dǎo)如下：令則現(xiàn)使用打靶法求解如下：令T為響應(yīng)的最小周期。根據(jù)runge-kutta方法推導(dǎo)如下Rnge-kutta的一般形式為：確定了階數(shù)之后，再通過Taylor展開，比較兩邊系數(shù)的方法，確定各待定系數(shù):.將展開式代入，得到局部截斷誤差：要使得方法是二階的，則局部截斷誤差應(yīng)該為三階小量，即：

4、。三個方程，四個未知數(shù)，所以其解不唯一，可令a=1, .則積分公式為：。類似前面的推導(dǎo)，可以導(dǎo)出各種四階的runge-kutta公式，它們的局部截斷誤差滿足。下面列出最常見的一個：從而根據(jù)次類似的推理，并將上式的迭代的次數(shù)n換成i,Y變成x.從而有下述結(jié)論由此可求得若（為允許誤差)則為所求非線性系統(tǒng)的周期解.T為最小響應(yīng)周期。Ww重點分析和結(jié)論：對于常微分方程經(jīng)過我們的推導(dǎo)和研究，我們發(fā)現(xiàn)等式x+(4+6cost)x-3x2+x+2x3=0.01cost;可改寫為x+(4+6cost)x-0.01cost=3x2-x-2x3=A,A為常數(shù)。x+(4+6cost)x-0.01cost=A有周期解

5、，3x2-x-2x3=A為常數(shù)解，故方程中3x2+x+2x3部分對整個方程的結(jié)的周期無影響。故求解x+(4+6cost)x-3x2+x+2x3=0.01cost的周期解時，我們可通過求解方程x+(4+6cost)x-0.01cost=0的周期解來實現(xiàn)。四，程序?qū)崿F(xiàn)：在MATLAB中編輯公式如下：function ys=dbf(f,a,b,alfa,beta,h,eps)ff=(x,y)y(2),f(y(1),y(2),x); xvalue=a:h:b; n=length(xvalue) ; s0=a-0.01;x0=alfa,s0;flag=0;y0=rk4(ff,a,x0,h,a,b); i

6、f abs(y0(1,n)-beta)<=eps flag=1; y1=y0; else s1=s0+1; x0=alfa,s1; y1=rk4(ff,a,x0,h,a,b); if abs(y1(1,n)-beta)<=eps flag=1; endendif flag=1 while abs(y1(1,n)-beta)>eps s2=s1-(y1(1,n)-beta)*(s1-s0)/(y1(1,n)-y0(1,n); x0=alfa,s2; y2=rk4(ff,a,x0,h,a,b); s0=s1; s1=s2;y0=y1;y1=y2; endendxvalue=a:h

7、:b;yvalue=y1(1,:); ys=xvalue',yvalue' function x=rk4(f,t0,x0,h,a,b) t=a:h:b; m=length(t); t(1)=t0; x(:,1)=x0; for i=1:m-1 L1=f(t(i),x(:,i); L2=f(t(i)+h/2,x(:,i)'+(h/2)*L1); L3=f(t(i)+h/2,x(:,i)'+(h/2)*L2); L4=f(t(i)+h,x(:,i)'+h*L3); x(:,i+1)=x(:,i)'+(h/6)*(L1+2*L2+2*L3+L4); e

8、nd 下達命令一：> f=(x,y,z)(0.01*cos(z)-(4+6*cos(z)*x);>> x01=0x01 = 0>> x0u=2;>> alfa=0;>> beta=2;>> h=0.01;>> dbf(f,x01,x0u,x01,x0u,h,1e-6);>> y=ans(:,2)；>> x=ans(:,1)；>> plot(x,y,'-r')得圖：命令二：f=(x,y,z)(0.01*cos(z)-(4+6*cos(z)*x);>> x0

9、1=0x01 = 0>> x0u=10;>> alfa=0;>> beta=2;>> h=0.01;>> dbf(f,x01,x0u,x01,x0u,h,1e-6);>> y=ans(:,2);>> x=ans(:,1);>> plot(x,y,'-r')得圖：命令三：f=(x,y,z)(0.01*cos(z)-(4+6*cos(z)*x);>> x01=4;>> x0u=8;>> alfa=0;>> beta=2;>> h

10、=0.01;>> dbf(f,x01,x0u,x01,x0u,h,1e-6);>> y=ans(:,2);>> x=ans(:,1);>> plot(x,y,'-r')得圖：同時可由命令三得數(shù)據(jù)：4.00003.7843；3.5686；3.3528；3.1370；2.9211；2.7050；2.48892.2727；2.0564；1.8399；1.6234；1.4068；1.1901；0.9733 0.7564；0.5394；0.3225；0.1054；-0.1116；-0.3286；-0.5456 -0.7625；-0.9793

11、；-1.1960；-1.4126；-1.6290；-1.8451 -2.0610；-2.2766；-2.4918；-2.7066；-2.9210；-3.1349 -3.3482；-3.5610；-3.7730；-3.9844；-4.1949；-4.4046 -4.6133；-4.8211；-5.0278；-5.2334；-5.4377；-5.6407 -5.8424；-6.0426；-6.2413；-6.4383；-6.6336；-6.8271；-7.0186 -7.2082；-7.3957；-7.5810；-7.7639；-7.9445；-8.1226 -8.2981；-8.4708；-8.

12、6408；-8.8078；-8.9718；-9.1326 -9.2902；-9.4444；-9.5952；-9.7423；-9.8857；-10.0253 -10.1610；-10.2926；-10.4201；-10.5433；-10.6621 -10.7763；-10.8860；-10.9909；-11.0910；-11.1862；-11.2763 -11.3612；-11.4408；-11.5150；-11.5838；-11.6470；-11.7045 -11.7562；-11.8020；-11.8419；-11.8757；-11.9033；-11.9247 -11.9398；-11.94

13、85；-11.9508；-11.9465；-11.9356；-11.9180 -11.8937；-11.8626；-11.8247；-11.7799；-11.7282；-11.6695 -11.6039；-11.5312；-11.4515；-11.3648；-11.2710；-11.1701 -11.0622；-10.9472；-10.8251；-10.6961；-10.5600；-10.4170 -10.2670；-10.1102；-9.9465；-9.7761；-9.5990；-9.4152 -9.2248；-9.0280；-8.8248；-8.6153；-8.3996；-8.1778；-

14、7.9501 -7.7166；-7.4774；-7.2326；-6.9824；-6.7270；-6.4665 -6.2010；-5.9308；-5.6561；-5.3769；-5.0935；-4.8062；-4.5150 -4.2202；-3.9221；-3.6208；-3.3166；-3.0097； -2.0753；-1.7600；-1.4433；-1.1253；-0.8064；-0.4869；-0.1669 0.1533；0.4733；0.7929；1.1118；1.4297；1.7464；2.0616 2.3749；2.6862；2.9950；3.3012；3.6045；3.9045；4

15、.2010 4.4937；4.7824；5.0666；5.3463；5.6210；5.8906；6.1547 6.4132；6.6657；6.9120；7.1518；7.3850；7.6112；7.8302 8.0418；8.2458；8.4420；8.6302；8.8101；8.9816；9.1445 9.2986；9.4438；9.5798；9.7066； 9.8240；9.9318；10.0300 10.1184；10.1969；10.2655；10.3240；10.3724；10.4106；10.4386 10.4563；10.4636；10.4607；10.4474；10.4237；

16、10.3898；10.3455 10.2909；10.2262；10.1512；10.0662；9.9711；9.8661；9.75129.6267；9.4925；9.3488；9.1958；9.0336；8.8624；8.68238.4935；8.2963；8.0907；7.8771；7.6556；7.4264；7.1898 6.9461；6.6954；6.4380；6.1742；5.9042；5.6284；5.3470 5.0602；4.7684；4.4719；4.1709；3.8658；3.5569；3.2444 2.9288；2.6102；2.2891；1.9657；1.6404；1.

17、3135；0.9853 0.6562；0.3264；-0.0038；-0.3339；-0.6637；-0.9928 -1.3210；-1.6479；-1.9732；-2.2966；-2.6178；-2.9366 -3.2524；-3.5652；-3.8746；-4.1804；-4.4821；-4.7797 -5.0727；-5.3610；-5.6443；-5.9223；-6.1948；-6.4615 -6.7223；-6.9768；-7.2250；-7.4666；-7.7013；-7.9290 -8.1496；-8.3628；-8.5684；-8.7664；-8.9566；-9.1387 -9

18、.3128；-9.4787；-9.6363；-9.7854；-9.9260；-10.0579 -10.1812；-10.2957；-10.4014；-10.4983；-10.5862；-10.6652 -10.7353；-10.7963；-10.8484；-10.8915；-10.9256；-10.9508 -10.9670；-10.9744；-10.9729；-10.9626；-10.9435；-10.9158 -10.8795；-10.8347；-10.7814；-10.7197；-10.6499；-10.5718 -10.4857；-10.3917；-10.2900；-10.1805；-

19、10.0635；-9.9391 -9.8074；-9.6686；-9.5228；-9.3702；-9.2110；-9.0453 -8.8732；-8.6950；-8.5107；-8.3207；-8.1249；-7.9238 -7.7173；-7.5057；-7.2891；-7.0678；-6.8420；-6.6117 -6.3773；-6.1388；-5.8965；-5.6505；-5.4011；-5.1483 -4.8925；-4.6338；-4.3723；-4.1083；-3.8418；-3.5732 -3.3026；-3.0301；-2.7559；-2.4801；-2.2031；-1.9248 -1.6455；-1.3654；-1.0845；-0.803