《数值分析》计算实习题姓名：学号：班级：第二章1、程序代码Clear;clc;x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=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.1643 q3 =- 1.3393*x^3 + 2.4107*x^2 - 1.6929*x + 1.4171 q4 =- 1.3393*x^3 + 3.2143*x^2 - 2.8179*x + 1.86293、问题结果4次牛顿差值多项式4()P x = 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三次样条差值多项式()Q x0.10.20.30.40.50.60.70.80.910.40.50.60.70.80.911.1323232321.33930.803570.40714 1.04,[0.2,0.4]1.3393 1.60710.88929 1.1643,[0.4,0.6]1.3393 2.4107 1.6929 1.4171,[0.6,0.8]1.3393 3.21432.8179 1.8629,[0.8,1.0]x x x x x x x x x x x x x x x x ⎧-+-+∈⎪-+-+∈⎪⎨-+-+∈⎪⎪-+-+∈⎩第三章1、程序代码Clear;clc;x=[0 0.1 0.2 0.3 0.5 0.8 1];y=[1 0.41 0.5 0.61 0.91 2.02 2.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.6221 12.8147 -4.6591 0.9266 p2 =2.8853 -12.3348 16.2747 -5.2987 0.9427 p3 =3.1316 -1.2400 0.73563、问题结果三次多项式拟合P1=32-6.622112.8147 4.65910.9266x x x +-+四次多项式拟合P2=4322.885312.334816.2747 5.29870.9427x x x x -+-+ 二次多项式拟合P3=23.1316 1.24000.7356x x -+第四章1、程序代码1)建立函数文件f.m: function y=fun(x); y=sqrt(x)*log(x); 2)编写程序：a. 利用复化梯形公式及复化辛普森公式求解：Clear;clc;h=0.001;%h 为步长,可分别令h=1,0.1,0.01,0.001 n=1/h;t=0;s1=0;s2=0; for i=1:n-1 t=t+f(i*h); endT=h/2*(0+2*t+f(1));T=vpa(T,7) %梯形公式0.10.20.30.40.50.60.70.80.9100.511.522.53for i=0:n-1s1=s1+f(h/2+i*h);endfor i=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;for i=1:mh=h/2;n=1/h;t=0;for j=1:n-1t=t+f(j*h);endT(i)=h/2*(0+2*t+f(1));%梯形公式endfor i=1:m-1for j=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、问题结果b. 利用龙贝格公式求解：通过15次二分，得到结果：-0.4444414第五章1、程序代码（1）LU分解解线性方程组：Clear;clc;A=[10 -7 0 1-3 2.099999 6 25 -1 5 -12 1 0 2];b=[8;5.900001;5;1];[m,n]=size(A);L=eye(n);U=zeros(n);flag='ok';for i=1:nU(1,i)=A(1,i);endfor r=2:nL(r,1)=A(r,1)/U(1,1);endfor i=2:nfor j=i:nz=0;for r=1:i-1z=z+L(i,r)*U(r,j);endU(i,j)=A(i,j)-z;endif abs(U(i,i))<epsflag='failure'return;endfor k=i+1:nm=0;for q=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）列主元消去法：function x = gauss(A,b);A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];b=[8;5.900001;5;1];[n,n] = size(A);x = zeros(n,1);Aug = [A,b]; %增广矩阵for k = 1:n-1[piv,r] = max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r r = r + k - 1; % 列主元所在大矩阵的行if r>ktemp=Aug(k,:);Aug(k,:)=Aug(r,:);Aug(r,:)=temp;endif Aug(k,k)==0, error(‘对角元出现0’), end% 把增广矩阵消元成为上三角for p = 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);for k = n-1:-1:1x(k)=b(k);for p=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.0000 0 0 0-0.0000 0.0000 0 00.0000 -2.5000 0.0000 00.0000 -2.4000 0.0000 0.0000U =1.0e+007 *0.0000 -0.0000 0 0.00000 -0.0000 0.0000 0.00000 0 1.5000 0.57500 0 0 0.0000x =-0.0000-1.00001.00001.0000detA =-762.00012）列主元消去法detA =762.0001ans =0.0000-1.00001.00001.00003、问题结果1）LU分解解线性方程组L=10003100 0.5250000010 0.224000000.961⎛⎫ ⎪-⎪ ⎪-⎪-⎝⎭U=10701006 2.3 00150000055749998.499 000 5.079998908 -⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭x=（0.0000，−1.0000，1.0000，1.0000）T detA=-762.0012)列主元消去法x=（0.0000，−1.0000，1.0000，1.0000）T detA=762.001第六章1、程序代码（1）Jacobi迭代Clear;clc;n = 6; %也可取n=8,10H = hilb(n);b = H * ones(n, 1);e = 0.00001;for i = 1:nif H(i, i)==0 '对角元为零，不能求解'returnendendx = zeros(n, 1);k = 0;kend = 10000;r = 1;while k<=kend & r>ex0 = x;for i = 1:ns = 0;for j = 1:i - 1s = s + H(i, j) * x0(j);endfor j = 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;endif k>kend '迭代不收敛，失败'else'求解成功'xkend（2）SOR迭代1）程序代码function s = SOR(n, w);H = hilb(n);b = H*ones(n, 1);e = 0.00001;for i = 1:nif H(i,i)==0 ‘对角线为零，不能求解’returnendendx = zeros(n, 1);k = 0;kend = 10000;r = 1;while k<=kend & r>ex0 = x;for i = 1:ns = 0;for j = 1:i - 1s = s + H(i, j) * x(j);endfor j = 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;endif k>kend '迭代不收敛，失败'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.mfunction f=g(x)f(1)=20/(x^2+2*x+10);2)建立函数文件：bdd.mfunction [y, n] = bdd(x, eps)if nargin==1eps=1.0e-8;elseif nargin<1errorreturnendx1 = g(x);n = 1;while (norm(x1-x)>=eps)&&(n<=10000) x = x1;x1 = g(x);n = n + 1;endy = x;n3)从命令窗口输入：bdd(0)（2）牛顿迭代clear;clc;format long;m=8; %m为迭代次数，可分别令m=2,4,6,8,10x=sym('x');f=sym('x^3+2*x^2+10*x-20');df=diff(f,x);FX=x-f/df; %牛顿迭代公式Fx=inline(FX);disp('x=');x1=0.5;disp(x1);Eps=1E-8;k=0;while 1x0=x1;k=k+1;x1=feval(Fx,x1); %将x1代入牛顿迭代公式替代x1 disp(x1); %在屏幕上显示x1if k==mbreak;endendk,x12、运行结果（1）不动点迭代法>> bdd(0)n =25ans =1.3688（2）牛顿迭代x=21.4666666666666671.3715120138059211.3688102226338951.3688081078226671.3688081078213731.3688081078213731.368808107821373k =8x1 =1.3688081078213733、问题结果(1)不动点迭代法x=1.3688 n=25 收敛太慢（2）牛顿迭代初值取0迭代次数k=8时，。