数值分析.doc

数值分析实验报告信息与计算科学02 李志刚 100920351.设,计算,，哪个接近。算法分析：(matlab语言)用Taylor级数计算。计算的误差为1=(其中介于0到-5),计算的误差为2=(其中介于0到5)，由于计算的误差较小，2的上界近似等于1的下界，故2误差较小。也可以由交错级数收敛得慢得出误差较大。程序如下：function lizhigangx=-5;he=1;jc=1;n=24;for k=1:n; jc=jc*k; he=he+1/jc*xk;endP=vpa(he,16)P =0.006737963097703459function lzgx=5;he=1;jc=1;n=24;for k=1:n; jc=jc*k; he=he+1/jc*xk;endP=vpa(he,16);Q=1/PQ = 0.006737947000163260结果如下： 故用计算的结果更精确。2.将不选列主元和列主元的Gauss消去法编成程序，并求解下面的84阶方程组。最后将你的结果与方程精确解比较，并谈谈你对Gauss消去法的看法。算法分析：(Java语言)在不选列主元的Gauss消去法中可能会出现小主元，这就会引起其他元素数量级的剧增和舍入误差的增长，导致计算结果不可靠。而选列主元的Gauss消去法有利于控制误差的增长，使计算结果更可靠。(其精确解为x都为1)不选列主元的Gauss程序如下:public class 第二题Guauss消去法 /构造一个处理增广矩阵的方法 public static double func(doubleA) final int n=A.length;for(int i=0;in-1;i+)for(int ii=i+1;iin;ii+)double rate=Aiii/Aii;for(int jj=i;jjn+1;jj+)Aiijj=Aiijj-rate*Aijj; return A; public static void main(String args) final int n=84;double Ab=new doublenn+1; /增广矩阵double b=new doublen;for(int i=0;in;i+) /初始化矩阵AAbii=6.0;for(int i=1;in;i+)Abii-1=8.0;for(int i=1;in;i+)Abi-1i=1.0;Ab0n=7; /初始化增广矩阵中的向量bAbn-1n=14;for(int i=1;in-1;i+)Abin=15.0;func(Ab);for(int i=0;i=0;i-)double sum=0;for(int j=i+1;jn;j+)sum+=Abij*xj;xi=(bi-sum)/Abii;for(int i=0;in;i+) /显示结果System.out.print(x+(i+1)+=);System.out.print(xi+t);if(i+1)%3=0) System.out.println();结果如下：x1=1.0000000000000000x2=1.0000000000000002x3=0.9999999999999996x4=1.0000000000000009x5=0.9999999999999982x6=1.0000000000000036x7=0.9999999999999929x8=1.0000000000000142x9=0.9999999999999716x10=1.0000000000000568x11=0.9999999999998863x12=1.0000000000002274x13=0.9999999999995453x14=1.0000000000009095x15=0.999999999998181x16=1.000000000003638x17=0.9999999999927245x18=1.000000000014551x19=0.999999999970898x20=1.000000000058204x21=0.9999999998835918x22=1.0000000002328164x23=0.9999999995343671x24=1.0000000009312657x25=0.9999999981374685x26=1.000000003725063x27=0.9999999925498742x28=1.0000000149002517x29=0.9999999701994966x30=1.0000000596010068x31=0.9999998807979864x32=1.0000002384040272x33=0.9999995231919456x34=1.0000009536161087x35=0.9999980927677825x36=1.000003814464435x37=0.99999237107113x38=1.00001525785774x39=0.9999694842845201x40=1.0000610314309597x41=0.9998779371380806x42=1.0002441257238388x43=0.9995117485523225x44=1.0009765028953548x45=0.9980469942092913x46=1.0039060115814138x47=0.9921879768371866x48=1.01562404632557x49=0.9687519073490876x50=1.0624961853009154x51=0.8750076294018072x52=1.2499847411818346x53=0.500030517694535x54=1.9999389643781136x55=-0.999877927824961x56=4.99975585192486x57=-6.999511688949468x58=16.99902331829793x59=-30.99804639819183x60=64.99609184276756x61=-126.9921798710706x62=256.9843444842836x63=-510.96862793713626x64=1024.9370117485487x65=-2046.873046994202x66=4096.742187976823x67=-8190.468751907319x68=16383.875007629336x69=-32764.50003051746x70=65531.00012207008x71=-131055.0004882807x72=262097.00195312407x73=-524127.0078124981x74=1048001.0312499963x75=-2094975.1249999925x76=4185857.499999985x77=-8355328.99999997x78=1.664512899999994E7x79=-3.3028126999999E7x80=6.50077449999997E7x81=-1.2582143899999955E8x82=2.34866688999999E8x83=-4.0262860699998E8x84=5.368381449999981E8(大概从处误差开始扩大,到了后面误差急剧增大。)选列主元的Gauss程序如下:public class 第二题列主元Guauss消去法 /构造一个用列主元Gauss消元法处理增广矩阵A(n*(n+1)的方法public static double func(doubleA) final int n=A.length;for(int i=0;in-1;i+)double tem1=Aii;int tem2=0;for(int p=i+1;ptem1)tem1=Api;tem2=p;if(tem2!=0)for(int q=i;qn+1;q+)double temm=Aiq;Aiq=Atem2q;Atem2q=temm;for(int ii=i+1;iin;ii+)double rate=Aiii/Aii;for(int jj=i;jjn+1;jj+)Aiijj=Aiijj-rate*Aijj; return A; public static void main(String args) final int n=84;double Ab=new doublenn+1; /增广矩阵double b=new doublen;for(int i=0;in;i+) /初始化矩阵AAbii=6.0;for(int i=1;in;i+)Abii-1=8.0;for(int i=1;in;i+)Abi-1i=1.0;Ab0n=7; /初始化增广矩阵中的向量bAbn-1n=14;for(int i=1;in-1;i+)Abin=15.0;func(Ab);for(int i=0;i=0;i-)double sum=0;for(int j=i+1;jn;j+)sum+=Abij*xj;xi=(bi-sum)/Abii;for(int i=0;i1.0e-3) x0=x1; x1=Bj*x0+gj; jacobitimes=jacobitimes+1;endx0jacobitimes %Gauss迭代x1=zeros(4,1); %构造初始向量x0=ones(4,1);gausstimes=0; %Gauss迭代次数while(norm(x1-x0)1.0e-3) x0=x1; x1=Bg*x0+gg; gausstimes=gausstimes+1;endx0gausstimes结果如下：Jacobi迭代： Gauss迭代： 从中可以看出，高斯-赛德尔迭代的收敛速度较快，仅为10次。4观察高次插值多项式的龙格现象。 给定函数。取等距结点，构造牛顿插值多项式和及三次样条插值函数。分别将三种插值多项式与的曲线画在同一个坐标系上进行比较。分析：（本题采用matlab语言）构造牛顿多项式时主要是构造差商矩阵；计算时，由于是等距结点，故简化了计算。程序如下:x5=linspace(-1,1,6); %6个插值结点y5=1./(1+25*x5.2); %6个在插值结点处的函数值F=zeros(6,7); %构造差商矩阵F，5*6阶矩阵F(:,1)=rot90(x5,3); %初始化F第一列F(:,2)=rot90(y5,3); %初始化F第二列 for j=3:7 %完整初始化F for i=(j-1):6 F(i,j)=(F(i,j-1)-F(i-1,j-1)/(F(i,1)-F(i-(j-2),1); endend p=; %用p向量存储F的关键系数for i=1:6 p(1,i)=F(i,i+1);end p=vpa(p,8); syms x;syms y;y=p(1);for i=2:6 j=2; q=1; while(j=i) q=q*(x-F(j-1,1); j=j+1; end y=y+p(i)*q;endsyms N5N5=vpa(simplify(y),3) x0=linspace(-1,1,200);y0=subs(y,x,x0);plot(x0,y0,b) %做N5图像 x=linspace(-1,1,1000);y=1./(1+25*x.2);hold onplot(x,y,r*,markersize,1) %画原函数图像grid on text(0.1,1,1/(1+25x2)图像);text(0.05,0.6,图像); %作N10，其基本算法同上x10=linspace(-1,1,11); %11个插值结点y10=1./(1+25*x10.2); %11个在插值结点处的函数值F=zeros(11,12); %构造差商矩阵F，11*12阶矩阵F(:,1)=rot90(x10,3); %初始化第一列F(:,2)=rot90(y10,3); %初始化第二列 for j=3:12 %完整初始化F for i=(j-1):11 F(i,j)=(F(i,j-1)-F(i-1,j-1)/(F(i,1)-F(i-(j-2),1); endend p=; %用p向量存储F的关键系数for i=1:11 p(1,i)=F(i,i+1);end p=vpa(p,8); syms x;syms y;y=p(1);for i=2:11 j=2; q=1; while(j=i) q=q*(x-F(j-1,1); j=j+1; end y=y+p(i)*q;endsyms n10N10=vpa(simplify(y),3) x0=linspace(-1,1,200);y0=subs(y,x,x0);plot(x0,y0) %作N10图像text(0.6,1.5,N10图像); %以下是做三次样条插值函数sx=linspace(-1,1,11); %等距取11个插值结点sy=1./(1+25*sx.2);h=2/10;u=0.5;v=1-u; A=; %三弯矩方程系数矩阵A=2.*eye(9,9);b=0.5.*ones(8,1);A=A+diag(b,1)+diag(b,-1); d=; %Ax=d,构造difor i=1:9 d(i,1)=6*F(i+2,4);end %求M0和M10syms x;y=1/(1+25*x2);M0=subs(diff(y,x,2),x,-1);M10=subs(diff(y,x,2),x,1); d(1,1)=d(1,1)-u*M0;d(9,1)=d(9,1)-v*M10; M=inv(A)*d;M=M0;M;M10; for i=1:10 t=linspace(sx(i),sx(i+1),100); s=(sx(i+1)-t).3/(6*h)*M(i)+(t-sx(i).3/(6*h)*M(i+1)+(sy(i)-h2/6*M(i). *(sx(i+1)-t)./h+(sy(i+1)-h2/6*M(i+1)*(t-sx(i)./h; syms S; syms x; S=(sx(i+1)-x)3/(6*h)*M(i)+(x-sx(i)3/(6*h)*M(i+1)+(sy(i)-h2/6*M(i). *(sx(i+1)-x)/h+(sy(i+1)-h2/6*M(i+1)*(x-sx(i)/h; S=vpa(expand(S),3) grid