数值分析计算实习题

《数值分析》计算实习题姓名：学号：班级：第二章1、程序代码Clear;clc;x1=[0.20.40.60.81.0];y1=[0.980.920.810.640.38];n=length(y1);c=y1(:);forj=2:n%求差商fori=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsymsxdfd;df(1)=1;d(1)=y1(1);fori=2:n%求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i-1)*df(i);endP4=vpa(sum(d),5)%P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1,'variational');%调用三次样条函数q=pp.coefs;q1=q(1,:)*[(x-.2)^3;(x-.2)^2;(x-.2);1];q1=vpa(collect(q1),5)q2=q(1,:)*[(x-.4)^3;(x-.4)^2;(x-.4);1];q2=vpa(collect(q2),5)q3=q(1,:)*[(x-.6)^3;(x-.6)^2;(x-.6);1];q3=vpa(collect(q3),5)q4=q(1,:)*[(x-.8)^3;(x-.8)^2;(x-.8);1];q4=vpa(collect(q4),5)%求解并化简多项式2、运行结果P4=0.98*x-0.3*(x-0.2)*(x-0.4)-0.625*(x-0.2)*(x-0.4)*(x-0.6)-0.20833*(x-0.2)*(x-0.4)*(x-0.8)*(x-0.6)+0.784q1=-1.3393*x^3+0.80357*x^2-0.40714*x+1.04q2=-1.3393*x^3+1.6071*x^2-0.88929*x+1.1643q3=-1.3393*x^3+2.4107*x^2-1.6929*x+1.4171q4=-1.3393*x^3+3.2143*x^2-2.8179*x+1.86293、问题结果4次牛顿差值多项式4()Px=0.98*x-0.3*(x-0.2)*(x-0.4)-0.625*(x-0.2)*(x-0.4)*(x-0.6)-0.20833*(x-0.2)*(x-0.4)*(x-0.8)*(x-0.6)+0.784三次样条差值多项式()Qx00.10.20.30.40.50.60.70.80.910.40.50.60.70.80.911.1已知的点牛顿插值多项式三次样条函数323232321.33930.803570.407141.04,[0.2,0.4]1.33931.60710.889291.1643,[0.4,0.6]1.33932.41071.69291.4171,[0.6,0.8]1.33933.21432.81791.8629,[0.8,1.0]xxxxxxxxxxxxxxxx第三章1、程序代码Clear;clc;x=[00.10.20.30.50.81];y=[10.410.50.610.912.022.46];p1=polyfit(x,y,3)%三次多项式拟合p2=polyfit(x,y,4)%四次多项式拟合y1=polyval(p1,x);y2=polyval(p2,x);%多项式求值plot(x,y,'c--',x,y1,'r:',x,y2,'y-.')p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。y3=polyval(p3,x);plot(x,y,'c--',x,y1,'r:',x,y2,'y-.',x,y3,'k--')%画出四种拟合曲线2、运行结果p1=-6.622112.8147-4.65910.9266p2=2.8853-12.334816.2747-5.29870.9427p3=3.1316-1.24000.73563、问题结果三次多项式拟合P1=32-6.622112.81474.65910.9266xxx四次多项式拟合P2=4322.885312.334816.27475.29870.9427xxxx二次多项式拟合P3=23.13161.24000.7356xx第四章1、程序代码1)建立函数文件f.m:functiony=fun(x);y=sqrt(x)*log(x);2)编写程序：a.利用复化梯形公式及复化辛普森公式求解：Clear;clc;h=0.001;%h为步长,可分别令h=1,0.1,0.01,0.001n=1/h;t=0;s1=0;s2=0;fori=1:n-1t=t+f(i*h);endT=h/2*(0+2*t+f(1));T=vpa(T,7)%梯形公式00.10.20.30.40.50.60.70.80.9100.511.522.53数据曲线三次多项式拟合四次多项式拟合二次多项式拟合fori=0:n-1s1=s1+f(h/2+i*h);endfori=1:n-1s2=s2+f(i*h);endS=h/6*(0+4*s1+2*s2+f(1));S=vpa(S,7)%辛普森公式a’复化梯形公式和复化辛普生公式程序代码也可以是：Clear;clc;x=0:0.001:1;%h为步长,可分别令h=1,0.1,0.01,0.001y=sqrt(x).*log(x+eps);T=trapz(x,y);T=vpa(T,7)（只是h=1的运行结果不一样，T=1.110223*10^(-16)，而其余情况结果都相同）Clear;clc;f=inline('sqrt(x).*log(x)',x);S=quadl(f,0,1);S=vpa(S,7)b.利用龙贝格公式求解：Clear;clc;m=14;%m+1即为二分次数h=2;fori=1:mh=h/2;n=1/h;t=0;forj=1:n-1t=t+f(j*h);endT(i)=h/2*(0+2*t+f(1));%梯形公式endfori=1:m-1forj=m:i+1T(j)=4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1);%通过不断的迭代求得T（j)，即T表的对角线上的元素。endendvpa(T(m),7)2、运行结果T=-0.4443875S=-0.4444345ans=-0.44444143、问题结果a.利用复合梯形公式及复合辛普森公式求解：步长h10.10.010.001梯形求积T=0[1.110223*10^(-16)]-0.4170628-0.4431179-0.4443875辛普森求积S=-0.3267527-0.4386308-0.4441945-0.4444345b.利用龙贝格公式求解：通过15次二分，得到结果：-0.4444414第五章1、程序代码（1）LU分解解线性方程组：Clear;clc;A=[10-701-32.099999625-15-12102];b=[8;5.900001;5;1];[m,n]=size(A);L=eye(n);U=zeros(n);flag='ok';fori=1:nU(1,i)=A(1,i);endforr=2:nL(r,1)=A(r,1)/U(1,1);endfori=2:nforj=i:nz=0;forr=1:i-1z=z+L(i,r)*U(r,j);endU(i,j)=A(i,j)-z;endifabs(U(i,i))epsflag='failure'return;endfork=i+1:nm=0;forq=1:i-1m=m+L(k,q)*U(q,i);endL(k,i)=(A(k,i)-m)/U(i,i);endendLUy=L\b;x=U\ydetA=det(L*U)（2）列主元消去法：functionx=gauss(A,b);A=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];[n,n]=size(A);x=zeros(n,1);Aug=[A,b];%增广矩阵fork=1:n-1[piv,r]=max(abs(Aug(k:n,k)));%找列主元所在子矩阵的行rr=r+k-1;%列主元所在大矩阵的行ifrktemp=Aug(k,:);Aug(k,:)=Aug(r,:);Aug(r,:)=temp;endifAug(k,k)==0,error(‘对角元出现0’),end%把增广矩阵消元成为上三角forp=k+1:nAug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k);endend%解上三角方程组A=Aug(:,1:n);b=Aug(:,n+1);x(n)=b(n)/A(n,n);fork=n-1:-1:1x(k)=b(k);forp=n:-1:k+1x(k)=x(k)-A(k,p)*x(p);endx(k)=x(k)/A(k,k);enddetA=det(A)2、运行结果1）LU分解解线性方程组L=1.0e+006*0.0000000-0.00000.0000000.0000-2.50000.000000.0000-2.40000.00000.0000U=1.0e+007*0.0000-0.000000.00000-0.00000.00000.0000001.50000.57500000.0000x=-0.0000-1.00001.00001.0000detA=-762.00012）列主元消去法detA=762.0001ans=0.0000-1.00001.00001.00003、问题结果1）LU分解解线性方程组L=100031000.52500000100.224000000.961U=107010062.300150000055749998.4990005.079998908x=（，，，）detA=-762.0012)列主元消去法x=（，，，）detA=762.001第六章1、程序代码（1）Jacobi迭代Clear;clc;n=6;%也可取n=8,10H=hilb(n);b=H*ones(n,1);e=0.00001;fori=1:nifH(i,i)==0'对角元为零，不能求解'returnendendx=zeros(n,1);k=0;kend=10000;r=1;whilek=kend&rex0=x;fori=1:ns=0;forj=1:i-1s=s+H(i,j)*x0(j);endforj=i+1:ns=s+H(i,j)*x0(j);endx(i)=b(i)/H(i,i)-s/H(i,i);endr=norm(x-x0,inf);k=k+1;endifkkend'迭代不收敛，失败'else'求解成功'xkend（2）SOR迭代1）程序代码functions=SOR(n,w);H=hilb(n);b=H*ones(n,1);e=0.00001;fori=1:nifH(i,i)==0‘对角线为零，不能求解’returnendendx=zeros(n,1);k=0;kend=10000;r=1;whilek=kend&rex0=x;fori=1:ns=0;forj=1:i-1s=s+H(i,j)*x(j);endforj=i+1:ns=s+H(i,j)*x0(j);endx(i)=(1-w)*x0(i)+w/H(i,i)*(b(i)-s);endr=norm(x-x0,inf);k=k+1;endifkkend'迭代不收敛，失败'else'求解成功'xend2）从命令窗口中分别输入：SOR（6,1）SOR（8,1）SOR（10,1）SOR（6,1.5）SOR（8,1.5）SOR（10，1.5）2、运行结果Jacobi迭代：ans=迭代不收敛，失败SOR迭代：第七章1、程序代码(1)不动点迭代法1）建立函数文件：g