94数值分析实验报告

石家庄经济学院2015/2016学年第1学期数值分析实验报告班级学号姓名指导教师１（一）实验题目：二分法求非线性方程根一、程序功能编程实现二分法求解非线性方程组的根。二、实验算例选择计算：01)(3xxxf的根。三、算法步骤1：事前执行次数估计步骤2：计算a,b的中点,判断是否跳出步骤3:a=m或b=m四、重要标识符说明a,b表示区间[a,b]f表示原函数五、程序运行实例六、源程序f=inline('x^3-x-1');a=1.0;b=1.5;c=(b-a)/2;２k=0;err=1e-5;y0=f(a);whilecepsa0=a;b0=b;x1=(a0+b0)/2;iff(x1)==0x1elseiff(x1)0a1=a0;b1=x1;elsea1=x1;b1=b0;endc=(b1-a1)/2;k=k+1;a=a1;b=b1;endx=(a+b)/2（一）实验题目：牛顿法求非线性方程根一、程序功能编程实现牛顿法求解非线性方程组的根。二、实验算例选择计算：01)(3xxxf的根。三、算法步骤1：求出迭代公式，原函数Function1，导数Function2。步骤2:迭代。步骤3:判断X(k+1)-X(k)与e的大小，决定是否终止。四、重要标识符说明k代表迭代次数，x代表迭代初值五、程序运行实例３六、源程序clearallError=1e-5;x=1;fork=1:10xk=x;x=x-(x^3-x-1)/(3*x^2-1);if(abs(xk-x)=Error)break;endend七、个人实验总结初步了解了matlab的用法。分别用二分法和牛顿法解非线性方程可清晰看出牛顿法求解要比二分简单的多且效率高。（二）实验题目：高斯列主元素消去法解方程组一、程序功能编程实现高斯列主元素消去法解方程组。４二、实验算例选择x1+2x2+x3=2-2x1-2x2-x3=-32x1-3x2-2x3=-1三、算法步骤1：选主元及从增广矩阵第一列中选择绝对值最大的数所在的行与第一行互换，并将新矩阵第二第三行的第一列数化0；步骤2：从新矩阵的第2,3行第二列中选择绝对值最大的数换到第二行，第三行二列的数化0；步骤3：回代，分别求出x3，x2，x1四、重要标识符说明A,b代表方程的系数矩阵nb的长度abs求绝对值A(k:n,k)取矩阵A中第k列n行符合条件的数Disp输出五、程序运行实例六、源程序A=[121;-2-2-1;2-3-2];b=[2-3-1];５disp('原方程为AX=b：')Abdisp('------------------------')n=length(b);fork=1:n-1[main,index]=max(abs(A(k:n,k)));index=index+k-1;ifabs(main)epsdisp('列元素太小！！');break;elseifindexktemp=A(k,:);A(k,:)=A(index,:);A(index,:)=temp;endfori=k+1:nm(i,k)=A(i,k)/A(k,k);A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);b(i)=b(i)-m(i,k)*b(k);endenddisp('消元后所得到的三角阵是')Ab(n)=b(n)/A(n,n);fori=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)))/A(i,i);endclearx;x=b;disp('AX=b的解x是')x６(二)实验题目：道立特三角分解法解方程组一、程序功能编程实现三角分解法解方程组。二、实验算例选择x1+2x2+x3=2-2x1-2x2-x3=-32x1-3x2-2x3=-1三、算法步骤1：计算L下三角矩阵；步骤2：计算U；步骤3：Ly=b,Ux=y求出x,y。四、重要标识符说明RA代表A的秩RB代表B的秩五、程序运行实例六、源程序A=[121;-2-2-1;2-3-2];b=[2;-3;-1];B=[Ab];n=length(b);RA=rank(A);RB=rank(B);７zhica=RB-RA;ifzhica0disp('请注意：因为RA~=RB，所以此方程组无解.')return;endifRA==RBifRA==ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.')X=zeros(n,1);forp=1:n-1fork=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);enddisp('AX=b的解x是')Xelsedisp('请注意：因为RA=RBn，所以此方程组有无穷多解.')endend七、个人实验总结高斯列主元消元法就是在高斯消元法的基础上选取每列最大值最为对角线上的元素，以减小误差；而三角分解法需求出上三角和下三角然后解出X,Y并不用选主元。但前者数值较为稳定。８（三）实验题目：雅可比解方程组一、程序功能编程实现雅可比法解方程组二、实验算例选择A=[-4,1,1,1;1,-4,1,1,1;1,1,-4,1;1,1,1,-4]B=[1;1;1;1];三、算法步骤1：求出矩阵的范数，判断迭代公式是否收敛。步骤2:代入雅克比的矩阵迭代公式x0迭代，得到解四、重要标识符说明A为方程组的系数矩阵b为方程组的常数向量X0为迭代初值Eps为误差限T为最大迭代次数，超过则认为无解。五、程序运行实例六、源程序functionjacobiclc;clear;A=[4-10-100;-14-10-10;0-1400-1;-1004-10;0-10-14-1;00-10-14];９B=[0;5;0;6;-2;6];Err_user=0.01;N=50;[m,n]=size(A);X=zeros(n,1);k=1;whilek=NXk=X;fori=1:nforj=1:nifi~=jAX(j)=A(i,j)*Xk(j);endendSum_AX=sum(AX);AX=0;X(i)=(B(i)-Sum_AX)/A(i,i);endE=max(abs(Xk-X));ifEErr_userbreak;endk=k+1;enddisp(X)disp(k)（三）实验题目：高斯—赛德尔解方程组一、程序功能编程实现雅可比法解方程组二、实验算例选择A=[-4,1,1,1;1,-4,1,1,1;1,1,-4,1;1,1,1,-4]B=[1;1;1;1]三、算法步骤1：求出矩阵的范数，判断迭代公式是否收敛。１０步骤2:代入高斯-赛德尔的矩阵迭代公式x0迭代，得到解四、重要标识符说明A为方程组的系数矩阵b为方程组的常数向量X0为迭代初值Eps为误差限T为最大迭代次数五、程序运行实例六、源程序functionx=GS(A,b,x0,eps,t)ifnargin==0;A=[-4111;1-411;11-41;111-4];b=[1111]';x0=[0000]';eps=1e-6;m=200;elseifnargin==3;eps=1e-6;m=200;elseifnargin3error('输入有误');return;elseifnargin==5m=t;endD=diag(diag(A));１１L=-tril(A,-1);U=-triu(A,1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;whilenorm(x-x0)=epsx0=x;x=B*x0+f;n=n+1;if(n=m)disp('迭代次数过多，可能不收敛');return;end;end七、个人实验总结高斯赛德尔是雅克比的一种变形，用新数据代替旧数据，加快了迭代速度。（四）实验题目：拉格朗日插值多项式一、程序功能用拉格朗日插值多项式的方法求xi的近似解二、实验算例选择xk0.50.70.91.11.31.51.71.9Sinxk0.47940.64420.78330.89120.96360.99750.99170.9463求函数x在0.60.81.0处的近似解１２三、算法步骤1：给定xi的值，计算dxnum=xi-X,即xi-x(1),xi-x(2),...dxden=x(i)-x(1:n)即xi-xk...步骤2：将步一代入计算公式Pn(x)=sinX0*L0(x)+...+sinXn*Ln(x)步骤3：计算结果四、重要标识符说明X代表函数中的自变量Y代表在X作用下的因变量Xi是插值点Prod是连乘函数五、程序运行实例六、源程序X=[0.50.70.91.11.31.51.71.9];Y=[0.47940.64420.78330.89120.96360.99750.99170.9463];xi=0.8;n=length(X);L=zeros(size(Y));fori=1:ndisp(i);dxnum=xi-X;dxden=X(i)-X(1:n);fork=1:nifi==kdisp(k)dxnum(i)=1;dxden(i)=1;１３lnum=prod(dxnum(1:n));lden=prod(dxden(1:n));L(i)=lnum/ldenendendendyi=sum(Y.*L)（四）实验题目：牛顿差值多项式一、程序功能用牛顿插值多项式的方法求xi的近似解二、实验算例选择xk0.50.70.91.11.31.51.71.9Sinxk0.47940.64420.78330.89120.96360.99750.99170.9463求函数x在0.60.81.0处的近似解三、算法步骤1：计算dxnum=xi-X,即xi-x(1),xi-x(2),...步骤2：构造差商表即D（i,j）步骤3：构造一个七次牛顿插值多项式，计算出结果四、重要标识符说明Xi是函数中的自变量Y代表在X作用下的因变量X是插值点Prod连乘函数Formatlong定义数据格式，使结果更精确五、程序运行实例１４六、源程序functionnewIntclc;clear;formatlong;Xi=[0.50.70.91.11.31.51.71.9];Y=[0.47940.64420.78330.89120.96360.99750.99170.9463];x=0.8;dx=x-Xi;n=length(Xi);yi(1)=Y(1);D(:,1)=Y(:);forj=2:nfori=j:nD(i,j)=(D(i,j-1)-D(i-1,j-1))/(Xi(i)-Xi(i-j+1));endyi(j)=D(j,j)*prod(dx(1:j-1));endPn=sum(yi);disp(Pn);七、个人实验总结拉格朗日插值多项式形式对称,计算较方便,但一旦有更新数据就必须重新计算。而牛顿插值就算增加节点也不需重新计算。１５（五）实验题目：曲线拟合的最小二乘法一、程序功能编程实现用最小二乘法做曲线拟合二、实验算例选择y=sinx三、算法步骤1：计算一组x,y的数据;步骤2：求出法方程步骤3：用递推公式拟合曲线方程四、重要标识符说明x是插入点横坐标y为所求函数M是拟合次数五、程序运行实例六、源程序pi=3.1415926;x=(pi/180)*linspace(0,2*180,2*18);y=sin(x);M=4;[m,n]=size(x);w=ones(1,n);１６c=zeros(1,n);a=zeros(1,n);b=zeros(1,n);p=zeros(M,n);p0=1;c(1)=sum(w.*y)/sum(w);a(1)=sum(w.*x)/sum(w);p(1,:)=p0*(x-a(1));c(2)=sum((w.*y).*p(1,:))/sum(w.*p(1