资源描述：
《常微分方程初值问题数值解法的比较.doc》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。
1、常微分方程初值问题数值解法的比较数值计算实践—课程设计报告课题名称常微分方程初值问题数值解法的比较完成时间2013-1-17姓名班级学号成绩一.实验目的及内容1实验目的：（1）了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究；（2）熟练掌握欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析；（3）比较欧拉法、改进欧拉法及龙格-库塔法，能够选择合适的方法进行问题的研究计算；2实验内容：求微分方程（欧拉法求解）求微分方程（改进欧拉法求解）求微分方程（龙格-库塔求解）根据实验的结果进行分析，了解一般方法的的优缺点，稳定性，收敛性以及截断误差的分析，针
2、对相应问题拿出有效方法得出最优的结果。二．相关背景知识介绍以及初值问题稳定性的研究：在科学与工程问题中，常微分方程表述物理量的变化规律，应用非常广泛，比如，天体运动的轨迹，机器人控制，化学反应过程的描述和控制以及电路瞬态过程分析等。这些问题中要求解随时间变化的物理量，即位置函数表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。考虑一阶常微分方程的初值问题，如果存在实数使得则称f关于y满足利普希茨条件，L称为利普希茨常数。15对于常微分方程初值问题，考虑初值的扰动是问题的解发生偏差的情形。若时的偏差被控制在有界范围内，则称该初值问题是稳定的，否
3、则该初值问题不稳定的。特别地，若时的偏差收敛于零，则称该初值问题是渐进稳定的。对于初值问题稳定性的研究，易知其准确解为，假定初值经过扰动后变为，对于扰动后的解为因此带来的扰动误差为，因此考虑时的值，它取决于。易知，若，则原问题是稳定的；若，则原问题是不稳定的；若，则原问题是渐进稳定的。实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面的讨论数值解法时这常常是默认条件。1．欧拉法：依据：积分曲线上一点的切线斜率等于函数值。方法：推进法，初始点出发，依照方向场在改点的方向推进到向前欧拉法的得到:(1)将在处泰勒展开取h的线性部分，得（2）将初值问题中得导数
4、用向前差商来代替有，因此（3）将两边同时对x的区间上积分对右端用左矩形公式得，此方法称向前欧拉法，也叫显示欧拉法。（4）对右端用右矩形公式得，也叫隐式欧拉法。15误差分析：1.称为计算时的局部截断误差；2.如果数值方法的局部截断误差为，那么称这种数值方法的阶数是p，其实p为非负整数。通常情况下，步长h越小，p越高，则局部截断误差越小，计算精初泰勒展开有则有可见欧拉方法是一阶方法，精度不是很高。1．改进欧拉方法：梯形公式：对右端用梯形公式得+显然梯形公式是隐式公式。改进欧拉公式：先用欧拉公式求的一个初步的近似值，成为预测值，预测值的精度可能达不到要求，在用梯形
5、公式将他校正一次，记为，这个结果成为校正值。预测：校正：误差分析：记为改进欧拉公式在处的截断误差，记因此有，表示在出的局部截断误差。由此得，梯形公式的局部截断误差为，因此改进欧拉的截断误差为，可见改进欧拉的方法是二阶方法，改进欧拉方法优于欧拉方法。2．龙格—库塔法:根据拉格朗日微分中值定理，，记得到，这样，只要给出一种计算的算法，就可以得到相应的计算公式。欧拉公式可以写为15改进欧拉公式可以写成因此推出一般的推出广式称为p阶龙格-库塔方法，简称p阶R-K方法。因为这里的均为常数。因为给定的系数不唯一，因此这里的常数有无穷多个解，下面是特殊情况下和一般情况下得
6、结果（一阶龙格-库塔）当r=1时，这就是欧拉法。（二阶龙格-库塔）当r=2时，，这就是改进欧拉法。三阶和四阶龙格-库塔也只是在一般情况下得结果。三阶四阶其局部截断误差是：15三．程序代码欧拉法代码如下：（1）向前欧拉法：functiony=Euler1(fun,x0,y0,xN,N)%fun为一阶微分方程，x0，y0为初始条件，xN为取值范围的一个端点，N为区间个数x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;%h为区间步长for(n=1:N)x(n+1)=x(n)+h;y(n+1)=y(
7、n)+h*feval(fun,x(n),y(n));%根据向前欧拉公式计算y值endT=[x',y']（2）向后欧拉法：functiony=Euler2(fun,x0,y0,xN,N)%fun为一阶微分方程，x0，y0为初始条件，xN为取值范围的一个端点，N为区间个数x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;%h为区间步长for(n=1:N)x(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n));for(k=1:3)z1=y(n)+h*feval(fu
8、n,x(n+1),z0);if(abs(z1-z0)
