2022年数值分析实验报告 .pdf

《2022年数值分析实验报告 .pdf》由会员分享，可在线阅读，更多相关《2022年数值分析实验报告 .pdf（23页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数值分析实验报告姓名：张献鹏学号：173511038专业：冶金工程班级：重冶二班精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 23 页数值分析实验报告目录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小组元的误差

2、. 84二分法求方程的根 . 94.1问题背景 . 94.2数学模型 . 94.3计算方法 . 94.4二分法的收敛性 . 105雅可比迭代求解方程组. 115.1问题背景 . 115.2数学模型 . 115.2.1理论基础 . 115.2.2实例 . 11精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 23 页数值分析实验报告5.3计算方法 . 125.4收敛性分析 . 136Romberg求积法 . 136.1问题背景 . 136.2数学模型： . 146.2.1理论基础 . 146.2.2实例 . 146.3计算方法 . 146.

3、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数学模型 . 198.4误差分析 . 20精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 23 页数值分析实验报告1 1拉格朗日插值1.1问题背景对于函数211)(xxf,55x求拉格朗日插值。10n，把)(xf和插值多项式的曲线

4、画在同一张图上进行比较，观察数值积分中的Lagrange插值。1.2数学模型取等距差值节点 ?=-5+10? /n，? =0，1，.，n，构造 n 次 lagrange插值多项式：?= 11 + ?2?=0?+1(?)(? - ?2)?+1(?)当 n=10 时，十次插值多项式L10(x)以及函数 f(x)的图像可以由 Matlab 画出。1.3计算方法f.m：function f= f( x )f=1./(1+x.2);end Lagrange.mfunction y=Lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:m z=x(i); s

5、=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); endend s=p*y0(k)+s; end y(i)=s;End 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 23 页数值分析实验报告2 拉格朗日插值的曲线：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) hold onplot(x0,y1,r) 运行这个文件可以得到如下

6、拉格朗日图形曲线：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问题背景用复化 Simpson公式计算定积分 ?2= ?21?的近似值，要求误差限 =1/2 10-7，利用其余项对算法做出步长的事前估

7、计；并将计算结果与精确值进行比较。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 23 页数值分析实验报告3 2.2数学模型将积分区间 a,b分为 n 等分，h=(b-a)/n，xk=a+ k h，k=0，1，n。在每个子区间 xk，xk+1上用 Simpson公式可得： ? (x)dx = ?(?)?6?-1?=0? ?(?) + 4? (?+12) + ?(?+1)?-1?=0其中 xk+1/2=xk+1/2h。Sn? =?6 ? (?) + 4?(?+12) + ? ( ?+1) =?-1?=0?6? (? ) + 4 ? (?+

8、12) + 2 ? (?) + ?(?)?-1?=1?-1?=0此式即为复化 Simpson公式。设 f(x) C4a,b，由 Simpson公式的误差有?= ? - ?= -190(?2)5?(4)(?)?-1?=0，?,?+1。则复化 Simpson公式的余项为：?= -?-?2880?4?(4)(?)，? ，?复化 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：clc clear syms

9、 x% 定义自变量 xf=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 （）函数求的是表达式的最小值，且要求精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 23 页数值分析实验报告4 表达式带引号，故取负号，一边求最大值e=5*10(-8) % 精度要求值a=1 % 积分下限

10、b=2 % 积分上限x1=fminbnd(f3,1,2) % 求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的 x值for n=2:1000000 % 求等分数 n Rn=-(b-a)/180*(b-a)/(2*n)4*f2(x1) % 计算余项if abs(Rn)0(i=0，1，. ，n)则求积公式是稳定的。3矩阵的 LU 分解3.1问题背景矩阵的 LU 分解主要用来求解线性方程组或者计算行列式。在使用初等行变换法求解线性方程组的过程中，系数矩阵的变化情况如下：A=12-1310-1-1-2经过 E12(-3)、E13(1)、E23(1/5)可得到 12-10-5300-12/5。

11、由上可知： 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)精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 23 页数值分析实验报告6 A = 12-1310-1-1-2=100310-1-1/51 12-10-5300-12/5=LU如果方阵 A 可以分解成单位下三角矩阵L 与上三角矩阵 U 的乘积，则式 A=LU

12、 称为 A 的 LU 分解或三角分解。3.2数学模型3.2.1理论基础矩阵的 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=Ly=b=Ax 的解。若矩阵 A 非奇异，则 A 能分解为 LU 的充分必要条件是A 的顺序主子行列式不为0。?1= ?110，?2= ?11?12?2

13、1?22，?3= ?11?1?1?则存在惟一的主对角线上元素全为1的下三角阵 L 与惟一的上三角阵 U， 使得 A=LU 。3.2.2实例将矩阵 1020302045803080171进行 LU 分解。3.3计算方法程序：clear allclc A=input( 请输入一个方阵 );% 输入一个 n阶方阵n,n=size(A); L=zeros(n,n); U=zeros(n,n); for i=1:n %将L的主对角线元素赋值1 L(i,i)=1; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 23 页数值分析实验报告7 endf

14、or j=1:n %求矩阵 U的第一行元素 U(1,j)=A(1,j); endfor k=2:n %求矩阵 L的第一列元素 L(k,1)=A(k,1)/U(1,1); endfor i=2:n %求L、U矩阵元素for j=i:n s=0; for t=1:i-1 s=s+L(i,t)*U(t,j); end U(i,j)=A(i,j)-s; endfor k=i+1:n r=0; for t=1:i-1 r=r+L(k,t)*U(t,i); end L(k,i)=(A(k,i)-r)/U(i,i); endend%输出矩阵 L、UL U输入一个方阵，输出结果如下：精选学习资料 - - -

15、- - - - - - 名师归纳总结 - - - - - - -第 10 页，共 23 页数值分析实验报告8 3.4小组元的误差例如线性方程组 Ax = b，其中： A = 0111，b=10，可得理论解 x=-11。如 果 系 数 矩 阵 被 扰 动 成 ? = 10-20111 ， 可 手 算 知 ： ?= 1010-201 ， ?=10-20101 - 10-20。若上述过程在计算机中进行， 由浮点数运算规则可知， 两数相加时，大数吃掉小数，则计算机中产生的矩阵为：?= ?，?= ?，?= 10-2010-10-20这时会发现 ?= 10-20110，且据?x=b 解出的理论解 x=01

16、，明显不再等于前面的理论解。这说明 LU 分解是稳定的，但是将LU 分解用到解线性方程组上是不稳定的。究其原因，是因为 ?中的第一个主元 10-20太小，导致第二个主元中的1 与值 10-20相差悬殊，出现大数吃小数。为了避免上述危害， 引入一种选主元手段， 即在消去的过程中，通过适当的选主元，避免放大数据误差。 常用的选主元技术就是列选主元法(除此之外还有全选主元法、 对角选主元法和随机选主元法等)：对 mn 阶矩阵 A，在确定第 k 个主元 ?(?)(?k)时，先从该列自主元位置 (k,jk)至列尾的所有元素中选择绝对值最大的元素，与?(?)交换，然后将 ?+1,?(?)， ?.?(?)化

17、为零。继续这个过程，直至将矩阵A 化为行阶梯形。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 23 页数值分析实验报告9 4二分法求方程的根4.1问题背景在科学研究与工程计算中，常遇到方程（组）求根问题。若干个世纪以来，工程师和数学家花了大量时用于探索求解方程（组） ， 研究各种各样的方程求解方法。 对于方程f(x)=0，当 f(x)为线性函数时， 称 f(x)=0 为线性方程； 当 f(x)为非线性函数时， 称式 f(x)=0为非线性方程。对于线性方程（组）的求解，理论与数值求法的成果丰富；对于非线性方程的求解，由于 f(x)的多

18、样性，尚无一般的解析解法。 当 f(x)为非线性函数时， 若 f(x)=0无解析解，但如果对任意的精度要求，设计迭代方程，数值计算出方程的近似解，则可以认为求根的计算问题已经解决，至少能够满足实际要求。4.2数学模型使用二分法求方程x3+x-1=0 在0,1内的近似根（误差 er xk=(a+b)/2; fx=subs(ff,x,xk); fa=subs(ff,x,a); k=k+1;if fx=0 y(k)=xk;break;elseif fa*fx0 b=xk;else a=xk;end y(k)=xk;end精选学习资料 - - - - - - - - - 名师归纳总结 - - - -

19、- - -第 12 页，共 23 页数值分析实验报告10 plot(y,.-);grid on在命令窗口下执行：实验结果如下：可以得到迭代区间中点数列分布及图像，数值如下：ans =0.5，0.75，0.625，0.6875，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二分法

20、的收敛性二分法算法简单， 容易理解，且总是收敛的； 但是二分法收敛速度太慢， 浪费时间，并且不能求复根跟偶数重根。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 23 页数值分析实验报告11 5雅可比迭代求解方程组5.1问题背景在自然科学和工程技术中有很多问题的解决常常归结为解线性方程组，例如电学中的网络问题，化学中的配平方程式模型问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，非线性方程组求解问题，用差分法或者有限元法解常微分方程、偏微分方程边值问题等都导致求解线性方程组。在实践中，通常采用数值方法来

21、讨论线性方程组Ax=b 的解，其中迭代法是一种重要方法。5.2数学模型5.2.1理论基础雅可比迭代法：首先将方程组中的系数矩阵A 分解成三部分，即： A = L+D+U，其中 D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。之后确定迭代格式，X(k+1) = B*X( k) +f， （这里 表示的是上标，括号内数字即迭代次数） ，其中 B 称为迭代矩阵，雅克比迭代法中一般记为 J（k= 0,1，.） ，再选取初始迭代向量X(0)，开始逐次迭代。设 Ax= b，其中 A=D+L+U 为非奇异矩阵，且对角阵D 也非奇异，则当迭代矩阵J的谱半径 (J)1 时，雅克比迭代法收敛。对于 Ax=b，A

22、= 0?21?0?1?2?30 +?11?22?+0?12?13?1?0?22?2?0?3?0=L+D+U因为?0(Jacob假设)所以 D 可逆。故有：(L+D+U) x=bDx=-(L+U) x+bx=-D-1(L+U)+D-1b5.2.2实例用雅可比迭代法解方程组：43035-10-14 ?1?2?3=2430-24精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 23 页数值分析实验报告12 5.3计算方法雅可比迭代法程序：function x,k,index=Jacobi(A,b,ep,it_max)if nargin4 it

23、_max=100;endif nargin3 ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while k=it_maxfor i=1:nif abs(A(i,i)1e-10 index=0;return;end y(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i)/A(i,i);endif norm(y-x,inf)eps k=k+1; h=h/2; Q=0;for i=1:M x=a+h*(2*i-1); Q=Q+subs(sym(f),findsym(sym(f),x);end T(k+1

24、,1)=T(k,1)/2+h*Q; M=2*M;for j=1:k T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j)/(4j-1);end tol=abs(T(k+1,j+1)-T(k,j);endI=T(k+1,k+1);step=k;输入：可以得到结果：6.4误差分析Romberg求积公式是具有八阶精度的算法，收敛且稳定，比梯形公式，Simpson公精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 18 页，共 23 页数值分析实验报告16 式以及 Cotes公式收敛的快。龙贝格公式的余项为：?,?(? ) = ?(?)?

25、- ?,?= -?2?+22(?+1 )( ?+2?)2?!(? - ? )2?+3?(2?+2 )(?)Romberg积分公式高速有效，易于编程，适合于计算机计算。但他有一个主要的缺点是，每当把积分区间对分后，就要对被积函数?(?)计算它在新分点处的值，而这些函数值的个数是成倍增加的。7幂法7.1问题背景工程及物理中的许多实际问题需要计算矩阵的特征值及特征向量，对于给定的 n 阶实矩阵 A，当 n 很小时用传统的方法是可以的，但当n稍大时，计算工作量将以惊人的速度增大， 并且由于计算存在误差， 用方程 (I-A)x=0求解十分困难。 在实际应用中， 有的时候不一定需要求出矩阵的全部特征值和特

26、征向量，例如只需要求出按模最大的或最小的特征值。求这类特征值的方法，通常采用迭代法，其中两种是幂法和反幂法。7.2数学模型7.2.1理论基础设 An 有 n 个线性相关的特征向量v1，v2，vn，对应的特征值1，2，n，满足|1| |2| |n| ，因为 v1，v2， vn为 Cn 的一组基，所以任给x(0)0，niiivax1)0(，所以有：)()(21111111)0(niiikiknikkiiniikiniiikkvavavavAavaAxA若 a1 0，则因11i知，当 k 充分大时 A(k)x(0)1ka1v1 = cv1属 1的特征向量另一方面，记 max(x)=xi，其中 |xi

27、| = |x|，则当 k 充分大时，111111111111111)0(1)0()max()max()max()max()max()max(vavavavaxAxAkkkkkk若 a1=0，则因舍入误差的影响，会有某次迭代向量在v1 方向上的分量不为0，迭代下去可求得 1及对应特征向量的近似值。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 19 页，共 23 页数值分析实验报告17 7.2.2实例用幂法计算矩阵?= 3-10-132023绝对值最大的特征值及相应的特征向量，取精度要求为=10-4。7.3计算方法幂法 Matlab 程序：funct

28、ion m,u,index=pow(A,ep,it_max)if nargin3 it_max=100;endif nargin2 ep=1e-5;endn=length(A);u=ones(n,1);index=0;k=0;m1=0;while k=it_max v=A*u;,i=max(abs(v); m=v(i);u=v/m;if abs(m-m1)ep index=1;break;end m1=m;k=k+1;end 在命令窗口输入和输出如下图所示：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 20 页，共 23 页数值分析实验报告18

29、7.4误差分析幂法程序可以用来计算矩阵绝对值最大的特征值及相应的特征向量。幂法的缺点是开始的时候并不知道矩阵是否有单一的主特征值，也不知道如何选择 v0以保证它关于矩阵特征向量的表达中包含一个与主特征值相关的非零特征向量。采用幂法和反幂法，求矩阵的最大和最小特征值，从原理上看，这两种方法都是迭代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。同时，原点位移 m 的选取也影响收敛的速度，但原点位移m0的适当选取依赖于对矩阵 A 的大致了解。8改进欧拉法8.1问题背景在工程和科学技术的实际问题中，常需求解微分方程，但常微分方程中往往只有少数较简单和典型的常微分方程 （例如

30、线性常系数常微分方程等）可求出其解析解 ,对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程的求解困难就更不用说了。大多数情况下，常微分方程只能用近似方法求解。 这种近似解法可分为两大类：一类是近似解析法， 如级数解法、 逐次逼近法等 ;另一类是数值解法， 它给出方程在一些离散点上的近似值。8.2数学模型8.2.1理论基础如果在计算中，先用Euler 公式求得一个初步的近似值yn+1（称为预测值），再用梯形公式将它校正一次，就称为预测-校正公式，即改进的Euler 方法，其迭代格式为预测：-1),(nnnnyxhfyy校正：)(),(21,11nnnnnnyxfyxfhyy8.

31、2.2实例求解初值问题 ?= -? + ? + 1?|?=0= 1，已知精确解为：y=x + e-x，h=0.1。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 21 页，共 23 页数值分析实验报告19 8.3数学模型改进欧拉法程序：clearclcX0=0;X1=1;n=10;h=1/n;y(1)=1;x(1)=X0;X=X0:h:X1;Y=X+exp(-X);Y=vpa(Y,9); %精确解XX=X0:0.0001:X1;YY=XX+exp(-XX);for i=1:n x(i+1)=x(i)+h; y(i+1)=y(i)+h*(-y(i)+

32、x(i)+1);endfor i=1:n y(i+1)=y(i)+0.5*h*(-y(i)+x(i)+1)+(-y(i+1)+x(i+1)+1);endy=vpa(y,9)plot(x,y,*)hold onplot(XX,YY)运行结果：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 22 页，共 23 页数值分析实验报告20 改进欧拉法图像：8.4误差分析改进的 Euler 法与 Euler 法相比较，改进的 Euler 法的计算精度更高，相对误差也比较小。因此在求解微分方程的数值解时，改进的Euler 法优于 Euler 法。改进的 Euler 法步长 h 越小则计算精度越高，相对误差较小。因此，计算能力允许的范围内，选取步长越小可以得到更加精确的结果。在利用改进的 Euler 法过程中，前一次迭代的结果会对下一轮求得的数值产生影响。因此，一旦上一轮迭代所得的结果有偏差，下一轮结果的偏差将大于上一轮的偏差。因此会导致伴随迭代次数的增加而产生更大的偏差。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 23 页，共 23 页