数值分析实验报告(共23页).docx

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

1、精选优质文档-倾情为你奉上数值分析实验报告姓名：张献鹏学号：专业：冶金工程班级：重冶二班专心-专注-专业目录1拉格朗日插值1.1问题背景对于函数,求拉格朗日插值。，把和插值多项式的曲线画在同一张图上进行比较，观察数值积分中的Lagrange插值。1.2数学模型取等距差值节点xk=-5+10k/n，k=0，1，.，n，构造n次lagrange插值多项式：Ln=i=0n11+xi2wn+1(x)(x-xi2)wn+1(xi)当n=10时，十次插值多项式L10(x)以及函数f(x)的图像可以由Matlab画出。1.3计算方法f.m：function f= f( x )f=1./(1+x.2);end

2、Lagrange.mfunction y=Lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:m z=x(i); s=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); end end s=p*y0(k)+s; end y(i)=s;End拉格朗日插值的曲线：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 on plot(x0,y

3、1,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问题背景用复化Simpson公式计算定积分e2=12xexdx的近似值，要求误差限=1/210-7，利用其余项对算法做

4、出步长的事前估计；并将计算结果与精确值进行比较。2.2数学模型将积分区间a,b分为n等分，h=(b-a)/n，xk=a+ k h，k=0，1，n。在每个子区间xk，xk+1上用Simpson公式可得：abfxdx=k=0n-1f(x)dxh6k=0n-1fxk+4fxk+12+f(xk+1)其中xk+1/2=xk+1/2h。Snf=h6k=0n-1fxk+4fxk+12+fxk+1=h6fa+4k=0n-1fxk+12+2k=1n-1fxk+f(b)此式即为复化Simpson公式。设f(x)C4a,b，由Simpson公式的误差有Rsn=I-Sn=k=0n-1-190h25f4(k)，kxk,

5、xk+1。则复化Simpson公式的余项为：Rsn=-b-a2880h4f4(k)，ka，b复化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 %定义自变量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=-

6、(4*exp(x) + x*exp(x) %因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值e=5*10(-8) %精度要求值 a=1 %积分下限b=2 %积分上限x1=fminbnd(f3,1,2) %求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值for n=2: %求等分数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分解主要用来求解线性方程组或者计算行列式。在使用初等行变换法求解线性方程组的过

7、程中，系数矩阵的变化情况如下：A=12-1310-1-1-2经过E12(-3)、E13(1)、E23(1/5)可得到12-10-5300-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)A=12-1310-1-1-2=-1-1/5112-10-5300-12/5=LU如果方阵A可以分解成单位下三角矩阵L与上三角矩阵U的乘积，则式A=LU称为A的LU分解或三角分解。3.2数学模型3.2.1

8、理论基础矩阵的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=a110，2=a11a12a21a22，3=a11a1nan1ann则存在惟一的主对角线上元素全为1的下三角阵L与惟一的上三角阵U，使得A=LU。3.

9、2.2实例将矩阵进行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;endfor 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)

10、=A(i,j)-s; end for 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、ULU输入一个方阵，输出结果如下：3.4小组元的误差例如线性方程组Ax = b，其中：A=0111，b=10，可得理论解x=-11。如果系数矩阵被扰动成=10-20111，可手算知：L=1010-201，U=10-20101-10-20。若上述过程在计算机中进行，由浮点数运算规则可知，两数相加时，大数吃掉小数，则计算机中产生的矩阵为：A=，L=L，U=10-2010-10-20这时会发

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

12、将ak+1,jk(k)，an.jk(k)化为零。继续这个过程，直至将矩阵A化为行阶梯形。4二分法求方程的根4.1问题背景在科学研究与工程计算中，常遇到方程（组）求根问题。若干个世纪以来，工程师和数学家花了大量时用于探索求解方程（组），研究各种各样的方程求解方法。对于方程f(x)=0，当f(x)为线性函数时，称f(x)=0为线性方程；当f(x)为非线性函数时，称式f(x)=0为非线性方程。对于线性方程（组）的求解，理论与数值求法的成果丰富；对于非线性方程的求解，由于f(x)的多样性，尚无一般的解析解法。当f(x)为非线性函数时，若f(x)=0无解析解，但如果对任意的精度要求，设计迭代方程，数值计

13、算出方程的近似解，则可以认为求根的计算问题已经解决，至少能够满足实际要求。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;endplot(y,.-);grid on在命令窗口下执行：实验结果如下：可以得到迭代区间中点数列分布及图像，数值如下：ans =0.5，0.75，0.625，0.6875，0.65625，0.，0.。0.，

14、0.，0.，0.，0.，0.，0.，0.，0.，0.依据题目要求的精度，则需做十七次，由实验数据知x=0.即为所求的根。4.4二分法的收敛性二分法算法简单，容易理解，且总是收敛的；但是二分法收敛速度太慢，浪费时间，并且不能求复根跟偶数重根。5雅可比迭代求解方程组5.1问题背景在自然科学和工程技术中有很多问题的解决常常归结为解线性方程组，例如电学中的网络问题，化学中的配平方程式模型问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，非线性方程组求解问题，用差分法或者有限元法解常微分方程、偏微分方程边值问题等都导致求解线性方程组。在实践中，通常采用数值方法来讨论线性方程

15、组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=0a210an1an2an30+a11a22ann+0a

16、12a13a1n0a22a2n0a3n0=L+D+U因为ann0(Jacob假设)所以D可逆。故有：(L+D+U)x=bDx=-(L+U)x+bx=-D-1(L+U)+D-1b5.2.2实例用雅可比迭代法解方程组：43035-10-14x1x2x3=2430-245.3计算方法雅可比迭代法程序：function x,k,index=Jacobi(A,b,ep,it_max)if nargin4 it_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_max

17、for i=1:n if 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); end if 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,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

18、(T(k+1,j+1)-T(k,j);endI=T(k+1,k+1);step=k;输入：可以得到结果：6.4误差分析Romberg求积公式是具有八阶精度的算法，收敛且稳定，比梯形公式，Simpson公式以及Cotes公式收敛的快。龙贝格公式的余项为：Rm,kf=abf(x)dx-Tm,k=-B2m+22m+1m+2k2m!b-a2m+3f2m+2()Romberg积分公式高速有效，易于编程，适合于计算机计算。但他有一个主要的缺点是，每当把积分区间对分后，就要对被积函数f(x)计算它在新分点处的值，而这些函数值的个数是成倍增加的。7幂法7.1问题背景工程及物理中的许多实际问题需要计算矩阵的特征

19、值及特征向量，对于给定的n阶实矩阵A，当n很小时用传统的方法是可以的，但当n稍大时，计算工作量将以惊人的速度增大，并且由于计算存在误差，用方程(I-A)x=0求解十分困难。在实际应用中，有的时候不一定需要求出矩阵的全部特征值和特征向量，例如只需要求出按模最大的或最小的特征值。求这类特征值的方法，通常采用迭代法，其中两种是幂法和反幂法。7.2数学模型7.2.1理论基础设An有n个线性相关的特征向量v1，v2，vn，对应的特征值l1，l2，ln，满足|l1| |l2| |ln| ，因为v1，v2，vn为Cn的一组基，所以任给x(0)0，所以有：若a1 0，则因知，当k充分大时A(k)x(0) l1

20、ka1v1 = cv1属1的特征向量另一方面，记max(x)=xi，其中|xi| = |x|，则当k充分大时，若a1=0，则因舍入误差的影响，会有某次迭代向量在v1方向上的分量不为0，迭代下去可求得1及对应特征向量的近似值。7.2.2实例用幂法计算矩阵A=3-10-绝对值最大的特征值及相应的特征向量，取精度要求为=10-4。7.3计算方法幂法Matlab程序：function 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;m

21、1=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在命令窗口输入和输出如下图所示：7.4误差分析幂法程序可以用来计算矩阵绝对值最大的特征值及相应的特征向量。幂法的缺点是开始的时候并不知道矩阵是否有单一的主特征值，也不知道如何选择v0以保证它关于矩阵特征向量的表达中包含一个与主特征值相关的非零特征向量。采用幂法和反幂法，求矩阵的最大和最小特征值，从原理上看，这两种方法都是迭代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度

22、上。同时，原点位移m的选取也影响收敛的速度，但原点位移m0的适当选取依赖于对矩阵A的大致了解。8改进欧拉法8.1问题背景在工程和科学技术的实际问题中，常需求解微分方程，但常微分方程中往往只有少数较简单和典型的常微分方程（例如线性常系数常微分方程等）可求出其解析解,对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程的求解困难就更不用说了。大多数情况下，常微分方程只能用近似方法求解。这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等;另一类是数值解法，它给出方程在一些离散点上的近似值。8.2数学模型8.2.1理论基础如果在计算中，先用Euler公式求得一个初步的近

23、似值yn+1（称为预测值），再用梯形公式将它校正一次，就称为预测-校正公式，即改进的Euler方法，其迭代格式为预测：校正：8.2.2实例求解初值问题dydx=-y+x+1y|x=0=1，已知精确解为：y=x + e-x，h=0.1。8.3数学模型改进欧拉法程序：clearclc X0=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)+x(i)

24、+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);end y=vpa(y,9)plot(x,y,*)hold onplot(XX,YY)运行结果：改进欧拉法图像：8.4误差分析改进的Euler法与Euler法相比较，改进的Euler法的计算精度更高，相对误差也比较小。因此在求解微分方程的数值解时，改进的Euler法优于Euler法。改进的Euler法步长h越小则计算精度越高，相对误差较小。因此，计算能力允许的范围内，选取步长越小可以得到更加精确的结果。在利用改进的Euler法过程中，前一次迭代的结果会对下一轮求得的数值产生影响。因此，一旦上一轮迭代所得的结果有偏差，下一轮结果的偏差将大于上一轮的偏差。因此会导致伴随迭代次数的增加而产生更大的偏差。