数值分析实验报告.docx

数值分析实验报告 姓名： XXX 学号： XXXXXXXX 学院： XXXXXXXX老师： XXXX 实验一一、实验内容用雅克比迭代法和高斯塞德尔迭代法求解课本例3.1，设置精度为10-6。二、实验公式 1、雅克比迭代法雅克比迭代法的基本思想：设方程组的系数矩阵的对角线元素，根据方程组推导出一个迭代公式，然后将任意选取的一初始向量代入迭代公式，求出，再以代入同一迭代公式，求出，如此反复进行，得到向量序列。当收敛时，其极限即为原方程组的解。2、高斯塞德尔迭代法：在雅可比（Jacobi）迭代法中，如果当新的分量求出后，马上用它来代替旧的分量，则可能会更快地接近方程组的准确解。基于这种设想构造的迭代公式称为高斯-塞德尔(Gauss-Seidel)迭代法。算法可相应地从雅可比（Jacobi）迭代法改造得到。3、两种迭代法的公式（1）雅可比迭代法：xi(k+1)=1/aii(bi-aijxj(k) （2）高斯-塞德尔迭代法：xi(k+1)=1/aii(bi-aijxj(k+1)-aijxj(k)三、实验程序1、雅克比迭代法的程序format long; %设置为长格式显示,显示15位小数A=8 -3 2;4 11 -1;6 3 12; %输入系数矩阵b=20 33 36; e=0.000001; %控制误差n=max(size(A); %测定维数for m=1:n if A(m,m)=0 对角元为零，不能求解 return endendx=zeros(n,1); %设置初始解k=0; %预设迭代次数为0kend=50; %最大迭代次数为50r=1; %前后项之差的无穷范数，初始值设为1while ke%达到预定精度或迭代超过50次推出计算 x0=x; %几下前次近似解 for m=1:n s=0; for l=1:m-1 s=s+A(m,l)*x0(l); end for l=m+1:n s=s+A(m,l)*x0(l); end x(m)=b(m)/A(m,m)-s/A(m,m); endr=norm(x-x0,inf); %重新计算前后项之差的无穷范数k=k+1;endif kkend 迭代不收敛，失败else 求解成功 x kendans =求解成功x = 2.999999988436388 1.999999749136996 0.999999871142208k =162、高斯塞德尔迭代法的源程序只要在上面的雅克比迭代法的程序中作一处更改便可，即将上面的第21行的程序s=s+A(m,l)*x0(l);改为s=s+A(m,l)*x(l);其结果如下：x = 2.999999983670503 2.000000019116444 1.000000003385637k = 9四、实验分析由实验可见对此题高斯赛德尔迭代法比雅可比迭代法的收敛速度快。因为由高斯赛德尔迭代法和雅克比迭代法的原理可知，雅克比迭代法只是利用了前次迭代的结果，而高斯赛德尔迭代法也利用了当前迭代过程中已经算出来的新的分量，即当新的分量求出后，马上用它来代替旧的分量，则可能会更快地接近方程组的准确解。因此一般情况下高斯赛德尔迭代法比雅可比迭代法的收敛速度快。实验二一、实验内容例5.1 判别方程的实根存在区间，要求区间长度不大于1，然后用二分法求出最小正根的近似值，精度。二、实验原理对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在a,b上连续，f(a)f(b)0，且f(x)在a,b内仅有一个实根x*，取区间中点c，若f(c)e，则c恰为其根，否则根据f(a)f(c) ezplot(x3-3*x+1);axis(-2 2 -10 10);grid由上图可知其根分布区间大概分布在(-2,-1),(0,1)和（1,2）中。编写二分法求解最小正根的近似值程序如下: format long;f=inline(x3-3*x+1);a=0;b=1;Eps=1E-5;for k=1:50A(k)=a;B(k)=b;ya=feval(f,a);yb=feval(f,b);temp=(a+b)/2;X(k)=temp;yt=feval(f,temp);F(k)=yt;if abs(yt)Eps break;endif yt*ya0 a=a;b=temp;elseif yt*yb B=;%清空变量Bx0=pi/6 pi/4 pi/3 pi/2; %初始数据y0=sin(x0);x=pi/18 2*pi/9; %插值点n=max(size(x0);for k=1:2l=1;y=y0(1);B(l,k)=y;s=1; %一次因子的乘积,预设为1dx=y0; %差商for i=1:n-1 dx0=dx; for j=1:n-i dx(j)=(dx0(j+1)-dx0(j)/(x0(i+j)-x0(j); end df=dx(1); s=s*(x(k)-x0(i); %一次因子乘积 y=y+s*df; %计算各次Newton插值的值 l=l+1; B(l,k)=y;%结果保存在变量B中endendB=1:l,B;disp(迭代过程中间结果);disp(k sin(pi/18) sin(pi/4.5);disp(B);迭代过程中间结果k sin(pi/18) sin(pi/4.5) 1.0000 0.5000 0.5000 2.0000 0.2239 0.6381 3.0000 0.1489 0.64344.0000 0.1659 0.6429四、结果分析牛顿插值可以在计算过程中，根据精度要求逐步增加节点，且计算量小。而Newton差商差值多项式有一个显著的优点，即它的每一项都是按x的指数作升幂排列，这样当需要增加节点提高差值多项式次数时，可以充分利用前面已经计算出的结果。如显然,k次牛顿插值公式仅在k-1阶牛顿插值公式后面增加一项(x- x0)L(x- xk-1) f(x0,x1 ,L,xk)，作为对k-1阶牛顿插值公式的补偿或修正，提高了插值精度。即当插值节点增加后,仅在原有计算结果的基础上增加一些项，增加节点前的所有计算结果均起作用，不会造成计算的浪费。克服了拉格朗日插值的弱点。实验四一、实验内容例 7.1 设有如下数据:13456789101054211234利用最小二乘法求该组数据的多项式拟合曲线.二、画图判断曲线的大致形状 x0=1:10;x0(2)= y0=10 5 4 2 1 1 2 3 4 figure,plot(x0,y0,x),title(数据观察图);hold on,plot(x0,y0,or);axis(0 11 0 11);x0 = 1 3 4 5 6 7 8 9 10y0 =10 5 4 2 1 1 2 3 4得到的图形如下：观察到图形大概为一条抛物线，因此它的最高次数为2。构造其正规方程的程序如下:三、 实验程序： format long;x0=1:10;x0(2)=;%初始数据%x为数据点的横坐标,y为数据点的纵坐标x=x0;y0=10 5 4 2 1 1 2 3 4;y=y0;m=2;%最高次数为2n=length(x);b=zeros(1,m+1);f=zeros(n,m+1);%f为正规方程的系数,初始为0for k=1:m+1 f(:,k)=x.(k-1);enda=f*f;b=f*y;%解方程,得到多项式由高到低的系数所构成的向量cc=ab;c=flipud(c);disp(c);c=c;%以下是对多项式系数进行重新计算,得到拟合的效果图x=0:0.1:11;F=c(1)*x.2+c(2)*x+c(3);plot(x0,y0, ok,LineWidth,2),grid on;hold on;plot(x,F,-r,LineWidth,1);title(拟合后的效果图);axis(0 11 0 11); 0.267570664629486 -3.605309396485843 13.459663865546169拟合后的效果如下：实验五一、实验内容例 9.1 设有微分方程: 分别用 Euler法、 Euler预测校正法作数值计算，并比较其计算结果（取h=0.1 ）。二、实验程序 %Euler法F=y-2*x/y;a=0;b=1;h=0.1;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;for i=2:n+1 x=X(i-1); y=Y(i-1); Y(i)=Y(i-1)+eval(F)*h;end%Euler预测-校正公式Y1=zeros(1,n+1);Y1(1)=1;for i=2:n+1 x=X(i-1); y=Y1(i-1); ty=Y1(i-1)+eval(F)*h;Y1(i)=Y1(i-1)+h/2*eval(F);x=X(i); y=ty; Y1(i)=Y1(i)+h/2*eval(F);end%准确解temp=;f=dsolve(Dy=y-2*x/y,y(0)=1,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 Euler法 Euler预测-校正公式 准确值);disp(X,Y,Y1,df);%画图观察效果figure;plot(X,df,k-,X,Y,-r,X,Y1,.-b);grid on;title(Euler法和Euler预测-校正法解常微分方程);legend(准确值,Euler法,Euler预测-校正法); 步长 Euler法 Euler预测-校正公式 准确值 0 1.000000000000000 1.000000000000000 1.000000000000000 0.100000000000000 1.100000000000000 1.095909090909091 1.095445115010332 0.200000000000000 1.191818181818182 1.184096569242997 1.183215956619923 0.300000000000000 1.277437833714722 1.266201360875776 1.264911064067352 0.400000000000000 1.358212599560289 1.343360151483998 1.341640786499874 0.500000000000000 1.435132918657796 1.416401928536909 1.414213562373095 0.600000000000000 1.508966253566332 1.485955602415668 1.483239697419133 0.700000000000000 1.580338237655217 1.552514091326145 1.549193338482967 0.800000000000000 1.649783431047711 1.61647