数值分析实验报告

实验一、误差分析一、实验目的1．通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；2．通过上机计算，了解误差、绝对误差、误差界、相对误差界的有关概念；3．通过上机计算，了解舍入误差所引起的数值不稳定性。二．实验原理误差问题是数值分析的基础，又是数值分析中一个困难的课题。在实际计算中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。同时，由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法的好坏会影响到数值结果的精度。三．实验内容对20,,2,1,0n，计算定积分105dxxxynn.算法1：利用递推公式151nnyny,20,,2,1n,取100182322.05ln6ln51dxxy.算法2：利用递推公式nnyny515111,,19,20n.注意到1010202010201051515611261dxxdxxxdxx,取008730.0)12611051(20120y.:四.实验程序及运行结果程序一：t=log(6)-log(5);n=1;y(1)=t;fork=2:1:20y(k)=1/k-5*y(k-1);n=n+1;endyy=0.0884y=0.0581y=0.0431y=0.0346y=0.0271y=0.0313y=-0.0134y=0.1920y=-0.8487y=4.3436y=-21.6268y=108.2176y=-541.0110y=2.7051e+003y=-1.3526e+004y=6.7628e+004y=-3.3814e+005y=1.6907e+006y=-8.4535e+006y=4.2267e+007程序2：y=zeros(20,1);n=1;y1=(1/105+1/126)/2;y(20)=y1;fork=20:-1:2y(k-1)=1/(5*k)-(1/5)*y(k);n=n+1;end运行结果：y=0.08840.05800.04310.03430.02850.02120.01880.01690.01540.01410.01300.01200.01120.01050.00990.00930.0089实验二、插值法一、实验目的1、理解插值的基本概念，掌握各种插值方法，包括拉格朗日插值和牛顿插值等，注意其不同特点；2、通过实验进一步理解并掌握各种插值的基本算法。二．实验原理插值法是函数逼近的一种重要方法，它是数值积分、微分方程数值解等数值计算的基础与工具，其中多项式插值是最常用和最基本的方法。拉格朗日插值多项式的优点是表达式简单明确，形式对称，便于记忆，它的缺点是如果想要增加插值节点，公式必须整个改变，这就增加了计算工作量。而牛顿插值多项式对此做了改进，当增加一个节点时只需在原牛顿插值多项式基础上增加一项，此时原有的项无需改变，从而达到节省计算次数、节约存储单元、应用较少节点达到应有精度的目的。三．实验内容1、已知函数表ix0.561600.562800.564010.56521iy0.827410.826590.825770.82495用二次拉格朗日插值多项式求5635.0x时的函数近似值。2、已知函数表ix0.40.550.650.80.9iy0.410750.578150.696750.888111.02652用牛顿插值多项式求)596.0(3N和)895.0(4N。四．实验程序及运行结果1.程序如下：function[y,R]=lagranzi(X,Y,x,M)n=length(X);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;q1=1.0;c1=1.0;forj=1:nifj~=kp=p*(z-X(j))/(X(k)-X(j));endq1=abs(q1*(z-X(j)));c1=c1*j;ends=p*Y(k)+s;endy(i)=s;endR=M*q1/c1;在MATLAB工作窗口输入程序：x=0.5635;M=1;X=[0.56280,0.56401];Y=[0.82659,0.82577];[y,R]=lagranzi(X,Y,x,M)实验结果如下：y=0.8261R=1.7850e-0072.程序如下：function[y,R]=newcz(X,Y,x,M)n=length(X);m=length(x);fort=1:mz=x(t);A=zeros(n,n);A(:,1)=Y';s=0.0;p=1.0;q1=1.0;c1=1.0;forj=2:nfori=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));endq1=abs(q1*(z-X(j-1)));c1=c1*j;endC=A(n,n);q1=abs(q1*(z-X(n)));fork=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endy(k)=polyval(C,z);endR=M*q1/c1;在MATLAB工作窗口输入程序：三次牛顿差值程序如下：X=[0.4,0.55,0.6];Y=[0.41075,0.57815,0.69675];[y,R]=newcz(X,Y,x,M)实验结果如下：y=0.6861R=6.0107e-006四次牛顿差值程序如下：x=0.895;X=[0.55,0.65,0.8,0.9];Y=[0.57815,0.69675,0.88811,1.02652];实验结果如下：y=1.0194R=1.6729e-006实验三、解线性方程组的直接法一、实验目的1．了解求线性方程组的直接法的有关理论和方法；2．会编制列主元消去法、LU分解法的程序；二．实验原理解线性方程组的直接法是指经过有限步运算后能求得方程组精确解的方法。但由于实际计算中舍入误差是客观存在的，因而使用这类方法也只能得到近似解。目前较实用的直接法是古老的高斯消去法的变形，即主元素消去法及矩阵的三角分解法。引进选主元的技巧是为了控制计算过程中舍入误差的增长，减少舍入误差的影响。一般说来，列主元消去法及列主元三角分解法是数值稳定的算法，它具有精确度较高、计算量不大和算法组织容易等优点，是目前计算机上解中、小型稠密矩阵方程组可靠而有效的常用方法。三．实验内容1.用列主元高斯消去法求解方程组035.3643.5072.1835.2137.2623.4712.3347.1183.1555.3304.2101.0321321321xxxxxxxxx.2．用矩阵直接三角分解法求解方程组bAx，其中381265973274581221A，4911427b.四．实验程序及运行结果1.主程序function[RA,RB,n,X]=liezhu(A,b)B=[Ab];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;ifzhica0,disp('请注意：因为RA~=RB，所以此方程组无解.')returnendifRA==RBifRA==ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.')X=zeros(n,1);C=zeros(1,n+1);forp=1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;fork=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);forq=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RBn，所以此方程组有无穷多解.')endend在MATLAB工作窗口输入程序：A=[0.1012.3043.555;-1.3473.7124.623;-2.8351.0725.643];b=[1.183;2.137;3.035];[RA,RB,n,x]=liezhu(A,b)运行结果RA=3RB=3n=3x=-0.39820.01380.33512.主程序functionX=LUjfcz(A,b)[n,n]=size(A);X=zeros(n,1);Y=zeros(n,1);C=zeros(1,n);r=1:n;forp=1:n-1[max1,j]=max(abs(A(p:n,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;g=r(p);r(p)=r(j+p-1);r(j+p-1)=g;ifA(p,p)==0disp('A是奇异阵,方程组无唯一解');break;endfork=p+1:nH=A(k,p)/A(p,p);A(k,p)=H;A(k,p+1:n)=A(k,p+1:n)-H*A(p,p+1:n);endendY(1)=b(r(1));fork=2:nY(k)=b(r(k))-A(k,1:k-1)*Y(1:k-1);endX(n)=Y(n)/A(n,n);fori=n-1:-1:1X(i)=(Y(i)-A(i,i+1:n)*X(i+1:n))/A(i,i);end在MATLAB工作窗口输入程序：A=[12-128;547-2;-3795;6-12-83];b=[27;4;11;49];X=LUjfcz(A,b)运行结果X=3.0000-2.00001.00005.0000实验四、解线性方程组的迭代法一、实验目的1、熟悉迭代法的有关理论和方法；2、会编制雅可比迭代法、高斯-塞德尔迭代法的程序；3、注意所用方法的收敛性及其收敛速度问题。二．实验原理解线性方程组的迭代法是用某种极限过程去逐步逼近线性方程组精确解的方法，即是从一个初始向量出发，按照一定的迭代格式产生一个向量序列，使其收敛到方程组的解。迭代法的优点是所需计算机存储单元少，程序设计简单，原始系数矩阵在计算过程中始终不变等。但迭代法存在收敛性及收敛速度问题。迭代法是解大型稀疏矩阵方程组的重要方法。三.实验内容1.用改进欧拉法解2)0(50'2yxxyy.2．用龙格-库塔方法求解2)0(50'2yxxyy四.实验程序及运行结果1.主程序functionX=jacdd(A,b,X0,P,wucha,max1)[nm]=size(A);forj=1:ma(j)=sum(abs(A(:,j)))-2*(abs(A(j,j)));endfork=1:max1kforj=1:mX(j)=(b(j)-A(j,[1:j-1,j+1:m])*X0([1:j-1,j+1:m]))/A(j,j);endX,djwcX=norm(X'-X0,P);xdwcX=djwcX/(norm(X',P)+eps);X0=X';X1=A\b;if(djwcXwucha)&(xdwcXwucha)disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：')returnendendif(djwcXwucha)&(xdwcXwucha)disp('请注意：雅可比迭代次数已经超过最大迭代次数max1')enda,X=X;jX=X1',在MATLAB工作窗口输入程序：A=[12-2;111;221];b=[7;2;5];X=jacdd(A,b,X0,P,wucha,max1)运行结果X=12-12.主程序functionX=gsdddy(A,b,X0,P,wucha,max1)D=diag(diag(A));U