数学模型期末作业.docx

重庆大学数学模型数学模型课程名称 用打靶法求周期解 学 院 物理学院 年级 2013 学 生 姓 名 和 学 号l 郑勇 20137022 l 赵茂雄 20136970l 席明 20136952 开 课 时 间 2014 至 2015 学年第 1 学期一，问题引入： 二，原理分析：1，打靶法原理：打靶法的实质是把微分方程的边值问题，转化成初值问题求解；求非线性动力系统的周期解就是，求如下形式的边值问题：T为响应的最小周期.。在某一给定的初始值表示初值迭代的次数)下具体计算过程如下:首先根据方法可得由此可求得若（为允许误差)则为所求非线性系统的周期解.若不然令,使用如下Newotn迭代关系求得下次计算需要的初值其中，微分算子在此关键要求出的值。令且方程两边对求导可得矩阵的初值问题其中，再次使用方法即可求得的值；代入式(4)即可求得下次计算需要的初值。在对式(5)使用Runge一Kutta方法求步的值时但在第1次的Runge一Kutta方法中,我们并没有求的式（7）需要的的值，其中的计算通常采用线性插值方法求得，即，或者重新缩短步长计算获得，但是在计算精度上与Runge一Kutta方法不匹配,且会增加计算量,必然在以后的计算中产生误差.为此我们只需要在第2次使用Runge-Kutta方法时直接对（2）式两边求导即可得第2次使用Runge-Kutta方法所需的值，具体过程如下：式（2）两边对求导可得下面求，在式（3）中令则有显然由于都是第1次使用Runge-Kutta方法求得的，而所以第步的值可表示为显然使用该种方法计算Runge-Kutta所需要的参数时候不需要的值，避免了误差的产生。2，Floquet理论：若系统(2)的特征根的模小于l,则系统(1)存在唯一的w一周期解。三，理论推理：的周期解的公式推导如下：令则现使用打靶法求解如下：令T为响应的最小周期。根据runge-kutta方法推导如下Rnge-kutta的一般形式为：确定了阶数之后，再通过Taylor展开，比较两边系数的方法，确定各待定系数:.将展开式代入，得到局部截断误差：要使得方法是二阶的，则局部截断误差应该为三阶小量，即：。三个方程，四个未知数，所以其解不唯一，可令a=1, .则积分公式为：。类似前面的推导，可以导出各种四阶的runge-kutta公式，它们的局部截断误差满足。下面列出最常见的一个：从而根据次类似的推理，并将上式的迭代的次数n换成i,Y变成x.从而有下述结论由此可求得若（为允许误差)则为所求非线性系统的周期解.T为最小响应周期。Ww重点分析和结论：对于常微分方程经过我们的推导和研究，我们发现等式x+(4+6cost)x-3x2+x+2x3=0.01cost;可改写为x+(4+6cost)x-0.01cost=3x2-x-2x3=A,A为常数。x+(4+6cost)x-0.01cost=A有周期解，3x2-x-2x3=A为常数解，故方程中3x2+x+2x3部分对整个方程的结的周期无影响。故求解x+(4+6cost)x-3x2+x+2x3=0.01cost的周期解时，我们可通过求解方程x+(4+6cost)x-0.01cost=0的周期解来实现。四，程序实现：在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); if 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 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: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); end 下达命令一： 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); x01=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=0.01; dbf(f,x01,x0u,x01,x0u,h,1e-6); y=ans(:,2); x=ans(:,1); plot(x,y,-r)得图：同时可由命令三得数据：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；-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.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.9485；-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；-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.7004-2.3888 -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.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；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.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.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；-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.8031；-0.5212；-0.2391 0.0431；0.3253