数值分析实验报告43278

1、数值分析实验报告姓名：张献鹏学号：173511038专业：冶金工程班级：重冶二班数值分析实验报告目录1拉格朗日插值 .11.1问题背景 .11.2数学模型 .11.3计算方法 .11.4数值分析 .22复化辛普森求积公式 .22.1问题背景 .22.2数学模型 .32.3计算方法 .32.4数值分析 .53矩阵的 LU 分解 .53.1问题背景 .53.2数学模型 .63.2.1理论基础 .63.2.2实例 .63.3计算方法 .63.4小组元的误差 .84二分法求方程的根 .94.1问题背景 .94.2数学模型 .94.3计算方法 .94.4二分法的收敛性 .105雅可比迭代求解方程组 .1

2、15.1问题背景 .115.2数学模型 .115.2.1理论基础 .115.2.2实例 .11数值分析实验报告5.3计算方法 .125.4收敛性分析 .136Romberg 求积法 .136.1问题背景 .136.2数学模型： .146.2.1理论基础 .146.2.2实例 .146.3计算方法 .146.4误差分析 .157幂法.167.1问题背景 .167.2数学模型 .167.2.1理论基础 .167.2.2实例 .177.3计算方法 .177.4误差分析 .188改进欧拉法 .188.1问题背景 .188.2数学模型 .188.2.1理论基础 .188.2.2实例 .188.3数学模型

3、 .198.4误差分析 .20数值分析实验报告1 拉格朗日插值1.1问题背景1f ( x), 5 x 5求拉格朗日插值。 n 10 ，把 f (x) 和插值多项式的对于函数1 x2曲线画在同一张图上进行比较，观察数值积分中的Lagrange 插值。1.2数学模型取等距差值节点 ?=-5+10?/n，?=0，1， .，n，构造 n 次 lagrange插值多项式：?1?(?)?= ?+122?1 + ? (?- ?)? (?)?=0?+1当 n=10 时，十次插值多项式L10(x)以及函数 f(x)的图像可以由 Matlab 画出。1.3计算方法f.m ：functionf= f( x )f=1

4、./(1+x.2);endLagrange.mfunctiony=Lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s;End1数值分析实验报告拉格朗日插值的曲线：x=-5:1:5;y=1./(1+x.2);x0=-5:0.001:5;y0=Lagrange(x,y,x0);y1=1./(1+x0.2);plot(x0,y0,'b')holdon

5、plot(x0,y1,'r')运行这个文件可以得到如下拉格朗日图形曲线：1.4数值分析L10(x)的断误差 R10(x)= f(x)-L10(x)在区间 -5，5的两端非常大。例如，L10(4.8)=1.80438，而 f(4.8)=0.04160。这种现象称之为龙格现象。不管n 取多大， Runge 现象依然存在。因此，对函数作插值多项式时，必须小心处理，不能认为差值节点取得越多，差值余项就越小。此外，当节点增多时，舍入误差的影响不能低估。为了克服高次插值的不足，应采用分段低次插值。2 复化辛普森求积公式2.1问题背景22?-7的近似值，要求误差限× ，用复化 Si

6、mpson 公式计算定积分 ?= =1/2?101利用其余项对算法做出步长的事前估计；并将计算结果与精确值进行比较。2数值分析实验报告2.2数学模型将积分区间 a,b 分为 n 等分，h=(b-a)/n，xk=a+ k h，k=0，1， n。在每个子区间 xk， xk+1 上用 Simpson公式可得：?-1?-1()dx =( )+ 4?(?1) + ?(?+1)?x?(?)? ?=06?=0?+2其中 xk+1/2=xk+1/2h。 ?-1()()?( )?-1?-1()S?+ 4?(?1 =?(?1) + 2+ ?(?)? =? ?) + ? + 4?n6?+2?+16?+2?=0?=1

7、?=0此式即为复化 Simpson 公式。设 f(x)C4a,b，由 Simpson 公式的误差有?= ?- ? =?-11?5()(?) ，? ?,?+1。-()?4?=0902则复化 Simpson 公式的余项为：?-?4( 4)?= - 2880 ? ? (?) ，? ?，?复化 Simpson 公式四阶收敛。2.3计算方法程序 1(求 f(x) 的 n 阶导数 )：syms xf=x*exp(x)%定义函数 f （ x）n=input(' 输入所求导数阶数:' )f2=diff(f,x,n)%求 f(x)的 n阶导数程序 2：clcclearsyms x%定义自变量 x

8、f=inline('x*exp(x)', 'x' )%定义函数 f(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=inline('(4*exp(x) + x*exp(x)', 'x' )%定义 f(x)的四阶导数，输入程序1里求出的f2 即可f3= '-(4*exp(x) + x*exp(x)'%因 fminbnd （）函数求的是表达式的最小值，且要求3数值分析实验报告表达式带引号，故取负号，一边求最大值e=5*10(-8)%精度要求值a=1%积分下限b=2%积分上限x1=fminbnd(f3,1,2

9、)%求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的 x值for n=2:1000000%求等分数 nRn=-(b-a)/180*(b-a)/(2*n)4*f2(x1)%计算余项ifabs(Rn)<e%用余项进行判断break% 符合要求时结束endendh=(b-a)/n%求 hSn1=0Sn2=0for k=0:n-1%求两组连加和xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b)%因 Sn2多加了 k=0时的值，故减去 f(a)z=exp(2)R=Sn-z%

10、求已知值与计算值的差fprintf(' 用 Simpson公式计算的结果Sn=' )disp(Sn)fprintf(' 等分数 n=' )disp(n)fprintf(' 已知值与计算值的误差R=' )disp(R)运行结果为：4数值分析实验报告2.4数值分析误差分析：在上述计算中，若采用复化梯形公式，则可以知，其与精确值的误差为-8-2.8300 ×e，等分数为 n=7019；而采用复化 Simpson 公式知与精确值的误差为2.7284 ×e8，等分数为 n=24。故与复化梯形公式相比，复化Simpson 公式误差相对较小

11、。?( )?复化 Simpson 公式的余项是：收敛性分析：若 lim = ( )? ?=0?，?= -?-? 4 ( 4)? ? (?) ，? ?，?2880?可以看出误差是 h4 阶，实际上若 f(x) C(a,b)，lim?，因此复化 Simpson? = ?(?)? ?公式是收敛的。稳定性分析：由于求积公式中Ai >0(i=0，1， .，n)则求积公式是稳定的。3 矩阵的 LU 分解3.1问题背景矩阵的 LU 分解主要用来求解线性方程组或者计算行列式。在使用初等行变换法求12-1解线性方程组的过程中，系数矩阵的变化情况如下：A= 310 经过 E12(-3)、-1-1-212-1

12、E13(1)、 E23(1/5)可得到 0-53 。00-12/5由上可知： E23(1/5)E13(1)E12(-3)A=U其中 U 就是上面矩阵 A经过行变换后的上三角矩阵， Eij 表示将 i行元素与 j 行元素互换的初等矩阵； Eij (k)表示将 i 行元素的 k 倍加到 j 行上。因此： A=E12(3)E13(-1)E23(-1/5)5数值分析实验报告12-110012-1A=310= 310 0-53 =LU-1-1-2-1-1/5100-12/5如果方阵 A 可以分解成单位下三角矩阵L 与上三角矩阵U 的乘积，则式 A=LU 称为 A 的 LU 分解或三角分解。3.2数学模型

13、理论基础矩阵的 LU 分解在求解线性方程组时将十分简便。 如对线性方程组 Ax=b，设 A=LU 是其 LU 分解。我们先求解方程组 Ly=b。由于 L 是下三角矩阵， 则解向量 y 可以通过依次求出其分量 y1，y2， ， yn 而求出，再求解方程组 Ux=y。解向量 x 可以通过该方程组依次求出分量 xn，xn-1， ， x2，x1 而快速得出。于是由两个方程组 Ux=y，Ly=b 的求解而给出 LUx=L y=b=A x 的解。若矩阵 A 非奇异，则 A 能分解为 LU 的充分必要条件是A 的顺序主子行列式不为0。?= ? 0，?= = 11?1?1112 ， ，? 1112?21?22

14、3?1?则存在惟一的主对角线上元素全为1 的下三角阵 L 与惟一的上三角阵 U，使得 A=LU 。实例102030将矩阵 204580 进行 LU 分解。30801713.3计算方法程序：clearallclcA=input(' 请输入一个方阵');% 输入一个 n阶方阵n,n=size(A);L=zeros(n,n);U=zeros(n,n);fori=1:n %将 L的主对角线元素赋值1L(i,i)=1;6数值分析实验报告endforj=1:n%求矩阵 U的第一行元素U(1,j)=A(1,j);endfork=2:n%求矩阵 L的第一列元素L(k,1)=A(k,1)/U(1

15、,1);endfori=2:n%求 L、U矩阵元素forj=i:ns=0;fort=1:i-1s=s+L(i,t)*U(t,j);endU(i,j)=A(i,j)-s;endfork=i+1:nr=0;fort=1:i-1r=r+L(k,t)*U(t,i);endL(k,i)=(A(k,i)-r)/U(i,i);endend%输出矩阵 L、 ULU输入一个方阵，输出结果如下：7数值分析实验报告3.4小组元的误差例如线性方程组 Ax = b ，其中： A = 011-111， b=0 ，可得理论解 x=1 。如果系数矩阵被扰动成?= 10-20110，?=11，可手算知：?= 10-201 10

16、-201 -20 。01- 10若上述过程在计算机中进行， 由浮点数运算规则可知， 两数相加时，大数吃掉小数，则计算机中产生的矩阵为：?， ?10-201= ? = ?，? = ?0-10-2010-20 0，明显不再等这时会发现 ?=1 ，且据 ? ? x=b解出的理论解 x =110于前面的理论解。这说明 LU 分解是稳定的，但是将LU 分解用到解线性方程组上是不稳定的。究其原因，是因为 ?中的第一个主元10-20 太小，导致第二个主元中的1 与值 10-20 相差悬殊，出现大数吃小数。为了避免上述危害，引入一种选主元手段， 即在消去的过程中，通过适当的选主元，避免放大数据误差。 常用的选

17、主元技术就是列选主元法 (除此之外还有全选主元法、 对角选主元法和随机选主元法等 )：(?)对 m× n 阶矩阵 A ，在确定第 k 个主元 ?(? k)时，先从该列自主元位置 (k,jk) 至?列尾的所有元素中选择绝对值最大的元素，(?)(?)(?)化为零。与?交换，然后将 ?， ?+1,?.?继续这个过程，直至将矩阵A 化为行阶梯形。8数值分析实验报告4 二分法求方程的根4.1问题背景在科学研究与工程计算中，常遇到方程（组）求根问题。若干个世纪以来，工程师和数学家花了大量时用于探索求解方程（组），研究各种各样的方程求解方法。对于方程f(x)=0，当 f(x)为线性函数时， 称 f

18、(x)=0 为线性方程； 当 f(x)为非线性函数时， 称式 f(x)=0 为非线性方程。对于线性方程（组）的求解，理论与数值求法的成果丰富；对于非线性方程的求解，由于 f(x)的多样性，尚无一般的解析解法。 当 f(x)为非线性函数时，若 f(x)=0 无解析解，但如果对任意的精度要求，设计迭代方程，数值计算出方程的近似解，则可以认为求根的计算问题已经解决，至少能够满足实际要求。4.2数学模型使用二分法求方程x3+x-1=0 在0,1内的近似根（误差 <10-5）。二分法：二分法是最简单的求根方法，它是利用连续函数的零点定理，将含根区间逐次减半缩小，取区间的中点构造收敛点列 xk 来逼

19、近根 x。4.3计算方法二分法程序代码：functiony=erfen1(m,n,er)symsxxka=m;b=n;k=0;ff=x3+x-1;whileb-a>erxk=(a+b)/2;fx=subs(ff,x,xk);fa=subs(ff,x,a);k=k+1;iffx=0y(k)=xk;break;elseiffa*fx<0b=xk;elsea=xk;endy(k)=xk;end9数值分析实验报告plot(y,'.-');gridon在命令窗口下执行：实验结果如下：可以得到迭代区间中点数列分布及图像，数值如下：ans =0.5，0.75，0.625，0.68

20、75，0.65625， 0.671875， 0.6796875。 0.68359375，0.68164062， 0.68261719，0.68212891， 0.68237305，0.68225098，0.68231201，0.68234253， 0.68232727，0.6823349依据题目要求的精度，则需做十七次，由实验数据知x=0.6823349 即为所求的根。4.4二分法的收敛性二分法算法简单， 容易理解，且总是收敛的；但是二分法收敛速度太慢， 浪费时间，并且不能求复根跟偶数重根。10数值分析实验报告5 雅可比迭代求解方程组5.1问题背景在自然科学和工程技术中有很多问题的解决常常归结

21、为解线性方程组，例如电学中的网络问题，化学中的配平方程式模型问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，非线性方程组求解问题，用差分法或者有限元法解常微分方程、偏微分方程边值问题等都导致求解线性方程组。在实践中，通常采用数值方法来讨论线性方程组Ax=b 的解，其中迭代法是一种重要方法。5.2数学模型理论基础雅可比迭代法：首先将方程组中的系数矩阵A 分解成三部分，即： A = L+D+U ，其中 D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。之后确定迭代格式， X( k+1) = B*X( k) +f，（这里 表示的是上标，括号内数字即迭代次数） ，其中

22、B 称为迭代矩阵，雅克比迭代法中一般记为 J（ k= 0,1，.），再选取初始迭代向量 X(0) ，开始逐次迭代。设 Ax= b，其中 A=D+L+U 为非奇异矩阵，且对角阵 D 也非奇异，则当迭代矩阵 J的谱半径 (J)<1 时，雅克比迭代法收敛。对于 Ax=b ，?0?012131?110?22?2?A = ?210 +?+220?3? ?1 ?2 ?30? ?0 ? =L+D+U因为 ? 0(Jacob假设 )?所以 D 可逆。故有：(L+D+U) x=bDx=-(L+U) x+bx=-D-1(L+U)+D -1b实例用雅可比迭代法解方程组：430?124 35-1 ?302=0-

23、14?-24311数值分析实验报告5.3计算方法雅可比迭代法程序：functionx,k,index=Jacobi(A,b,ep,it_max)ifnargin<4it_max=100;endifnargin<3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;whilek<=it_maxfori=1:nifabs(A(i,i)<1e-10index=0;return;endy(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i)/A(i,i);endifnorm(y-x,inf

24、)<epbreak;endk=k+1;x=y;end在命令窗口输入的命令和结果如下图：12数值分析实验报告5.4收敛性分析由上面运行的结果可知方程组的解为4.2000，2.4000，-5.4000，迭代次数为 39，由index=1 表示计算成功。雅克比迭代法的优点明显，计算公式简单，每迭代一次只需计算一次矩阵和向量的乘法，且计算过程中原始矩阵A 始终不变，比较容易并行计算。 然而这种迭代方式收敛速度较慢，而且占据的存储空间较大，所以工程中一般不直接用雅克比迭代法，而用其改进方法。与逐次超松弛迭代法相比，雅可比迭代法收敛速度相对较慢，逐次超松弛迭代法收敛速度较快。由逐次超松弛迭代法求出的

25、方程组的数值解与该方程组的精确解十分接近。并且离散化后线性方程组的逐次超松弛迭代法的精确性较高，故相对于雅可比迭代法，逐次超松弛迭代法更加广泛地应用于实际，可以用逐次超松弛迭代法求解高阶稀疏线性方程组。6Romberg 求积法6.1问题背景复化求积方法对于提高精度是行之有效的方法， 但复化公式的一个主要缺点在于要事先估计出部长。若步长过大，则精度难于保证；若步长过小，则计算量又不会太大。而用复化公式的截断误差来估计步长，其结果是步长往往过小， 而且f'' (x)和f(4) ( x)在13数值分析实验报告区间 a, b 上的上界 M 的估计是较为困难的。在实际计算中通常采用变步长

26、的方法，即把步长逐次分半（也就是把步长二等分），直到达到某种精度为止，这种方法就是Romberg 积分法的思想。6.2数学模型：6.2.1理论基础记：TN(0)TN ，将区间 N 等分的梯形值。 TN(1)SN ，将区间 N 等分的 Simpson；TN(2)C N ，将区间 N 等分的 Cotes。 TN(3)RN ，将区间 N 等分的 Romberg。由其可构造一个序列 TN(k) ，次序列称为 Romberg序列，并满足如下递推关系：(0)b a(0)1(0)b aNb aT1 f (a)f (b), T2 N2TN2Nf (a (2 k 1),2k 12Nk4k T2(Nk 1)TN(

27、k 1), k1,2,LTN4k1以上递推公式就是Romberg 积分递推公式。实例1.24利用 Romberg 积分法的思想，用龙贝格公式数值积分求函数? ?的积分。-1.26.3计算方法龙贝格公式积分程序：functionI,step=romberg(f,a,b,eps)%romberg.m为用龙贝格求积分%f 为被积函数%a.b为积分区间的上下限%eps 为积分结果精度%I 为积分结果；step为积分的子区间数if(nargin=3)eps=1.0e-4;end ;M=1;tol=10;14数值分析实验报告k=0;T=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(