数值分析课程设计.doc

数值分析课程设计报 告 专业： 信息与计算科学 学号： 091101215 学生姓名： 彭旭芳 指导教师： 缪红益 一、 问题提出：用不同的方法于下面的方程，挑出你满意的算法。1） y=(y+xy)/x,y(1)=1,计算x=2.2） y=0.5(x+4y)-x,y(2)=-1,计算到x=4.方法：1. 欧拉方法 2. 三阶龙格库塔公式二、理论基础下面主要讨论一阶常微分方程的初值问题，其一般形式为： （1）我们知道，只要函数f (x,y)适当光滑譬如关于y 满足利普希茨(Lipschitz)条件f (x,y) f (x,) L y理论上就可以保证初值问题(1)的解y=y(x)存在并且唯一。所谓数值解法， 就是求问题( 1 ) 在某些离散点a=xx x =b的近似值y,y,yy的方法。y,y,yy就称为问题(1)的数值解。h = xx 成为到x到x的步长，我们为了方便取为常量h 。欧拉方法将微分方程离散化，用向前差商代替微分y(x)，代入(1)中的微分方程，可得=f(x,y(x) (n=1,2,3)化简可得y (x)= y(x)+ f(x,y(x)h (n=1,2,3)如果用y近似y(x)代入上式便可得到y(x)的近似值y，计算式为：y = y + f(x, y)h (2)这样问题(1)的近似解可通过求解下面的差分初值问题： (n=1,2,3) (3) 得到，按(3)式由初值y可逐次求出y,yy。Eule 方法就是用差分方程初值问题(3)的解来近似微分方程初值问题(1)的解。即由公式 (3)算出y(x) 的近似值。这组公式求问题(1)的数值解就是著名的欧拉(Euler)公式。改进的欧拉方法用数值积分方法离散化问题(1)，两端积分可得y -y =f(x,y(x)dx (n=1,2,3)对右端积分使用梯形公式可得,f(x,y(x)dx f(x, y(x)+f(x,y(x)再用y, y代替y(x)，y(x)则得计算公式y =y+f(x, y)+f(x,y) (7)很明显可以注意到(7)式为隐式形式。改进的欧拉方法是先用欧拉公式求y(x)的一个近似值，称为预测值，然后用梯形公式进行矫正并求得近似值y。即 (8)改进的欧拉方法的误差估计由(7)式可知，y(x)- y(x)- y(x)+ y (x)=hy(x)+y(x)- y (x)+ y(x)+ hy”(x)+=-+ O(h)所以改进的欧拉方法是二阶的，精度较欧拉方法要高，实用性更加广泛。龙格-库塔方法的基本原理是：用方程中函数在前一节点上取值的线性组合构造一个表示的近似值公式，从而避免了求时用的高阶导数。该方法有多种推导方法，这里用数值积分法推导。为此，先将微分方程略加变形得出：，对该式两边在相邻两节点和之间求积分，移项得出： (3-15)采用不同的方法计算式(3-15)中定积分，便可得出数值积分的不同近似结果。如果用矩形求积公式计算，代入式(3-15)就可得出，这个结果和由Taylor公式得出的Euler公式结果完全一样。三阶龙格-库塔公式若用抛物线(Simpson)求积公式计算式(3-15)中的积分，则有：式中，而和都是未知的，可以用欧拉格式估算和。类似二阶龙格-库塔公式的推导，并令：，把它们代入式(3-15)则得出三阶龙格-库塔公式，它的局部截断误差为。 （3-17）二、 实验内容用改进的欧拉方法解y=(y+xy)/x,y(1)=1,计算到x=2.结果如下x = Columns 1 through 3 1.000000000000000 1.100000000000000 1.200000000000000 Columns 4 through 6 1.300000000000000 1.400000000000000 1.500000000000000 Columns 7 through 9 1.600000000000001 1.700000000000001 1.800000000000001 Columns 10 through 11 1.900000000000001 2.000000000000001y = Columns 1 through 3 1.000000000000000 1.214049586776860 1.463023084659781 Columns 4 through 6 1.753893874434661 2.095715940460330 2.500440717359425 Columns 7 through 9 2.984154699578947 3.569011699297371 4.286357858166035 Columns 10 through 11 5.182000314097829 6.325542138921840用改进的欧拉方法解y=0.5(x+4y)-x,y(2)=-1,计算到x=4结果如下x = Columns 1 through 3 2.000000000000000 2.200000000000000 2.400000000000000 Columns 4 through 6 2.600000000000001 2.800000000000001 3.000000000000001 Columns 7 through 9 3.200000000000001 3.400000000000001 3.600000000000001 Columns 10 through 11 3.800000000000002 4.000000000000002y = Column 1-1.000000000000000 Column 2-1.420000000000000 + 0.043588989435407i Column 3-1.864815806729992 + 0.154578681815478i Column 4-2.335513215091207 + 0.302190870011942i Column 5-2.832810379104140 + 0.480424194155551i Column 6-3.357061940013564 + 0.686119295671758i Column 7-3.908505582487285 + 0.917276040308307i Column 8-4.487324281537418 + 1.172501943336878i Column 9-5.093669390660537 + 1.450768159759122i Column 10-5.727671332019924 + 1.751282584504018i Column 11-6.389445354304673 + 2.073416594846941i用三阶龙格-库塔公式解y=(y+xy)/x,y(1)=1,计算到x=2.结果如下x = y= 1.000000000000000 1.000000000000 1.100000000000000 1. 215807374731411 1.200000000000000 1. 467336379130580 1.300000000000000 1. 761926153294213 1.400000000000000 2. 109135358348288 1.500000000000000 2. 521678923847024 1.600000000000000 3. 016829085059781 1.700000000000000 3. 618638982351298 1.800000000000000 4. 361611090780754 1.900000000000000 5. 297035294574241 2.000000000000000 6. 505701359717193用三阶龙格-库塔公式解y=0.5(x+4y)-x,y(2)=-1,计算到x=4结果如下x = y = 2.0000 -1.0000 2.2000 -1.4135 + 0.0599i 2.4000 -1.8539 + 0.1721i 2.6000 -2.3213 + 0.3200i 2.8000 -3.3377 + 0.7035i 3.0000 -3.8870 + 0.9343i 3.2000 -3.8870 + 0.9343i 3.4000 -4.4639 + 1.1889i 3.6000 -5.0685 + 1.4666i 3.8000 -5.7008 + 1.7665i 4.0000 -6.3611 + 2.0879i 三、 结果分析解第一个方程根据实验数据画出如下表格，取四位有效数字x1.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.0000Y欧拉1.00001.21401.46301.75392.09572.50042.98423.56904.28645.18206.3255Y龙格1.00001.21581.46731.76192.10912.52173.01683.61864. 36165.29706.5057解第二个方程根据实验数据画出如下表格，其中复数只取实部，取五位有效数字X2.00002.20002.40002.60002.80003.00003.20003.40003.60003.80004.0000Y欧拉-1.0000-1.4200-1.8648-2.3355-2.8328-3.3571-3.9085-4.4873-5.0936-5.7276-6.3894Y龙格-1.0000-1.4135-1.8539 -2.3213-2.8354-3.3377-3.8870-4.4639-5.0685-5.7008-6.3611通过对欧拉法，改进的欧拉法和龙格-库塔方法的研究，我们发现龙格-库塔方法的精度远高于欧拉法和改进的欧拉法。改进的欧拉法比欧拉法精度高的原因在于，它在确定平均斜率是，多取了一个斜率的值。如果我们在【x,x】上多取几个斜率值，然后对它们作线性组合得到平均斜率，则有可能得到精度更高的计算方法。而这就是龙格-库塔方法的基本思想。龙格-库塔方法是欧拉方法思想的提高，属于精度较高的单步法。我们从表中发现，它们的值还是比较接近的。五、参考文献1数值分析，李庆扬，王能超，易大义，清华大学出版社（第四版），2001。2数值方法，关治，陆金甫，清华大学出版社，2006。3数值分析与实验学习指导，蔡大用，清华大学出版社，2001。4数值分析与实验，薛毅，2005，北京工业大学出版社。5 数值计算原理M，王高熊，周之铭，朱思铭等。清华大学出版社，20006 数值分析M李庆扬，王能超，易大义清华大学出版社，施普林格出版社，20017 数值计算原理M李庆扬，关治，白峰杉清华大学出版社，2000附录改进的欧拉方法代码如下：function x,y=eulerpro(fun,x0,xfinal,y0,n);if nargin5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i);y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;endfunction x,y=eulerpro(fun,x0,xfinal,y0,n);if nargin5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i);y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end三阶龙格-库塔代码如下：function y = DELGKT3_kuta(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2); K3 = Funval(f,varvec,x(i-1)+h y(i-1)-h*K