数值分析实验报告

数值分析实验报告学院：姓名：实验一一、实验名称：线性方程组的直接解法二、实验目的：1、理解Gauss列主元素法和LU分解法的基本原理2、掌握Gauss列主元素法和LU分解法3、编写MATLAB程序实现Gauss列主元素法和LU分解法三、实验步骤及内容：使用Gauss列主元素法和LU分解法求解下列方程组：0.61.31.150.0020.036.05.40.296.05.696.00.5321321321xxxxxxxxx1．基本计算公式：（一）Gauss列主元素法：主要包括消元和回代步骤：[1]消元过程：主要是把原方程组化为三角方阵的方程组也即kkikikaal/kkjikkijkijalaa1kkikkikblbb1i与Gauss消元法的区别在于在每次消元前寻找该列绝对值最大的元素作为主元，然后进行消元。[2]回代过程：按变量的逆序逐步回代可得到方程组的解：有nnnnnnabx/kkklklkklkkkaxabx/)(n1）（,1,22,-n1,-n根据以上公式可以求出方程的解。（二）LU分解法：其基本思想是：对系数矩阵A进行分解，ULA,分别得到L为下三角矩阵，U为上三角矩阵，最后通过两个三角方程组yUxL和by，最后得出方程组的解。2、基本程序设计：列主元素法消元法的MATLAB基本程序如下所示：function[x,det,flag]=Gauss(A,b)%求线性方程组的列主元Gauss消元法，其中%A为方程组的系数矩阵；%b为方程组的右端项；%x为方程组的解；%det为系数矩阵A的行列式的值；%flag为指标向量，flag='failure'表示计算失败，flag=‘OK’表示计算成功。A=[5.00.966.5;2.04.50.36;0.501.13.1];b=[0.960.0206.0];[n,m]=size(A);nb=length(b)%当方程组行与列的维数不相等时，停止计算，并输出出错信息ifn~=merror('TherowandcolumnsofmatrixAmustbeepual!');return;end%当方程组与右端项的维数不匹配时，停止计算，并输出错误信息ifn~=nberror('ThecolumnsofAmustbeequalthelengthofb!')return;end%开始计算，先赋初值flag='OK';det=1;x=zeros(n,1);fork=1:n-1%选主元max1=0;fori=k:nifabs(A(i,k))max1max1=abs(A(i,k));r=i;endendifmax11e-10flag='failure';return;end%交换两行ifrkforj=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;end%消元过程fori=k+1:nm=A(i,k)/A(k,k);forj=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);%回代过程ifabs(A(n,n))1e-10flag='failure';return;endfork=n:-1:1forj=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);endx(k)=b(k)/A(k,k);endLU分解法求解线性方程组程序代码：functiony=isexist(A,b)[m,n]=size(A);[m1,n1]=size(b);ifm~=m1error('您的输入有误!!!');return;endr=rank(A);s=rank([A,b]);ifr==s&r==ny=1;elseifr==s&rny=Inf;elsey=0;endfunctionx=LUSeparation(A,b)flag=isexist(A,b);ifflag==0disp('该线性方程组无解!!!');x==[];return;elser=rank(A);[m,n]=size(A);[L,U,p]=lu(A);b=p*b;y(1)=b(1);ifm1fori=2:my(i)=b(i)-L(i,1:i-1)*y(1:i-1)';endendy=y';x0(r)=y(r)/U(r,r);ifr1fori=r-1:-1:1x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i);endendx0=x0';ifflag==1x=x0;return;elseformatrat;Z=null(A,'r');[m2,n2]=size(Z);x0(r+1:n)=0;fori=1:n2;t=sym(char([10748+i]));k(i)=t;endx=x0;fori=1:n2x=x+k(i)*Z(:,i);endendend四、实验结果：由列主元素消元法解的方程组的解为：在命令窗口中输入：LU分解法运行结果：五、总结：列主元素法相比高斯消元法而言，其避免了小主元的出现的问题，抑制了舍入误差的增长，使得计算结果相比而言误差较小。但与全主元素法相比，它的精度稍低于全主元素法，但是计算简单，工作量大为减少，所以列主元素法是求解小型的线性方程组的非常好的方法之一。三角分解法作为求解方程组直接方法之一，由于系数举证分解时与右端项无关，因而计算系数矩阵相同而右端项不同的方程组时，用三角求解法更为简单。实验二一、实验名称：线性方程组的求解的间接方法二、实验目的：1、理解Jacobi迭代法和Gauss-Seideld迭代法的基本原理2、掌握Jacobi迭代法和Gauss-Seideld迭代法3、编写MATLAB程序实现Jacobi迭代法和Gauss-Seideld迭代法的求解三、实验步骤及内容：试用（1）Jacobi迭代法和（2）Gauss-Seideld迭代法，求解下列的线性方程组。4258321072210321321321xxxxxxxxx选取初始向量为T)0,0,0,0,0(0x。1、基本原理：（1）Jacobi迭代法计算公式为iinijjkjijiikiabixaax111若用系数矩阵A来表示则有：A-1D-IBb1Dg其中),,,,(D11332211aaaadiag（2）Gauss-Seideld迭代法矩阵表示法为：bxk1-1-1kL-DUL-Dx）（）（2、Jacobi迭代法求解线性方程组的MATLAB程序如下所示：clearfprintf('Jacobi')x1_(1)=0;x2_(1)=0;x3_(1)=0;fori=1:9x1_(i+1)=7.2+0.1*x2_(i)+0.2*x3_(i);x2_(i+1)=8.3+0.1*x1_(i)+0.2*x3_(i);x3_(i+1)=8.4+0.2*x1_(i)+0.2*x2_(i);endx=[x1_',x2_',x3_']3、Gauss-Seideld迭代法求解线性方程组的MATLAB程序代码如下：clearfprintf('gauss-seidel迭代法')x1_(1)=0;x2_(1)=0;x3_(1)=0;fori=1:9x1_(i+1)=7.2+0.1*x2_(i)+0.2*x3_(i);x2_(i+1)=8.3+0.1*x1_(i+1)+0.2*x3_(i);x3_(i+1)=8.4+0.2*x1_(i+1)+0.2*x2_(i+1);endx=[x1_',x2_',x3_']四、实验结果：Jacobi迭代法运行结果：Gauss-Seideld迭代法运行结果：五、总结：由上面的计算明显可以看出，用两种方法解的最后的结果几乎是一致的，在误差容许的范围之内。但是用Gauss-Seide迭代法比Jacobi迭代法的效果好，说明了在此条件下Gauss-Seide迭代法比Jacobi迭代法的收敛速度要快。这也是在有些场合下应用较为广泛的原因。实验三一、实验名称：插值多项式的收敛性二、实验目的：1．理解插值的基本原理；2．掌握多项式插值的概念、存在唯一性；3．编写MATLAB程序实现Lagrange插值和Newton插值，验证Runge现象、分析插值多项式的收敛性。三、实验内容及步骤：1．已知数据如下：ix0.20.40.60.81.0()ifx0.980.920.810.640.38（1）用MATLAB语言编写按Langrage插值法和Newton插值法计算插值的程序，对以上数据进行插值；（2）利用MATLAB在第一个图中画出离散数据及插值函数曲线。2．给定函数f(x)=11+25x2,x∈[−1,1]，利用上题编好的Langrage插值程序（或Newton插值程序），分别取3个，5个、9个、11个等距节点作多项式插值，分别画出插值函数及原函数()fx的图形，以验证Runge现象、分析插值多项式的收敛性。2、Lagrange插值法（1）Lagrange插值法的基本思想：步骤1：构造01,,,nxxxL处的插值基函数01(),(),,()nlxlxlxL，其中，插值节点ix处的插值基函数()ilx为011011()()()()()()()()()iiniiiiiiinxxxxxxxxlxxxxxxxxx-+-+----=----LLLL；步骤2：以iy作为()ilx的系数，使得()iiylx通过插值点(,)iixy；步骤3：把所有的()iiylx线性叠加，得到通过所有插值点(,),0,1,,iixyin=L的插值函数0()()nniiiLxylx==å。（2）、Lagrange插值伪代码：给定n个插值点0011(,),(,),,(,)nnxyxyxyL的情况下，求插值函数()nLx在点t处的函数值。/*输入参数*y=(y0,y1,…,yn);被插函数f(x)在插值节点处的函数值*t求插值函数Ln(x)在t处的函数值*返回值插值函数Ln(x)在t处的函数值*/procedureLagrangeresult0;fori=1tonli(t)1;forj=1tonifi≠jli(t)li(t)*(t-xi)/(xi-xj);endifendforresultresult+yi*li(t);endforreturnresult;endprocedure（3）Lagrange插值子程序lagr1:functiony=lagr1(x0,y0,x)%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;End3、Newton插值算法公式:001001011()()[,]()[,,,]()()()nnnNxfxfxxxxfxxxxxxxxxLLL余项式()()()nnRxfxNx)(],,,,[010ininxxxxxxf0(1)()()(1)!niifnxxn其中(,)ab与x有关.（1）Newton插值伪代码：/*输入参数*x=(x0,x1….,xn),插值节点*y=(y0,y1,…,yn);被插函数f(x)在插值节点处的函数值*t求插值函数Pn(x)在t处的函数值*返回值插值函数Pn(x)在t处的函数值*/procedureNewtonforj=0tond1jyj;endforforj=1tonfori=jtondij(di,j-1-di-1,j-1)/(xi-xi-j+1);endforendforresultd11;temp1;fori=1ton