数学模型期末作业.docx

重庆大学数学模型

数学模型

课程名称 用打靶法求周期解
学 院 物理学院 年级 2013

学 生 姓 名 和 学 号
l 郑勇 20137022
l 赵茂雄 20136970
l 席明 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
flag=1;
end
end
if 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;
end
end
xvalue=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=0
x01 =

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=0

x01 =

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.0000
3.7843；3.5686；3.3528；3.1370；2.9211；2.7050；2.4889
2.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.7512
9.6267；9.4925；9.3488；9.1958；9.0336；8.8624；8.6823
8.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；0.6073；0.8891；1.1704；1.4512；1.7313
2.0106；2.2890；2.5664；2.8427；3.1177；3.3914；3.6637
3.9345；4.2037；4.4711；4.7368；5.0007；5.2626；5.5226
5.7805；6.0362；6.2898；6.5412；6.7903；7.0370；7.2814
7.5234；7.7629；8.0000
数据分析：周期t=0.01*239=2.39
五，Floquet理论判段稳定性：

由题意可知特征值的模<=1所以周期稳定！
六，误差分析与总结:
1. 下达三个命令是为了选择合适的或比较明显的周期；
2. 由于知识的有限我们采用的是读点计数法来计算的周期t；
3. 对于打靶法的原理及分析我们有借鉴，但更多的是我们的理解；
4. 主要存在误差，读点计数和点的选取；
5. 小组成员查阅了大量的资料及文献；
七，文章检索：
1，基于MATLAB的机械振动系统响应求解，柴保明，王远东，华龙，河北工程大学；

2，求解非线性动力系统周期解的改进打靶法，夏志鹏，郑铁生，“力学与实践”；

3，关于非线性周期系的Floquet理论，史金麟，“数学学报”；

4，Floquet理论在研究周期解中的应用，梁家荣，张棣，“广西师院学报（自然科学版）”；

5，求解非线性动力系统周期解推广的打靶法，李德信，徐健学，“应用力学学报”；附录文献查找：

19