数值分析计算实习习题

1、 数值分析计算实习题姓名： 学号： 班级： 第二章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:j c(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) %P

2、4即为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(

3、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.784 q1 = - 1.3393*x3 + 0.80357*x2 - 0.40714*x + 1.04 q2 = - 1.3393*x3 + 1.6071*x2 - 0.88929*x + 1.1643 q3 = - 1.3393*x3 + 2.4107*x2 - 1.6929*x + 1.41

4、71 q4 = - 1.3393*x3 + 3.2143*x2 - 2.8179*x + 1.86293、问题结果4次牛顿差值多项式= 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三次样条差值多项式 第三章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

5、)%三次多项式拟合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.814

6、7 -4.6591 0.9266p2 = 2.8853 -12.3348 16.2747 -5.2987 0.9427p3 =3.1316 -1.2400 0.73563、问题结果三次多项式拟合P1=四次多项式拟合P2=二次多项式拟合P3= 第四章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.001n=1/h;t=0;s1=0;s2=0;for i=1:n-1 t=t+f(i*h);endT

7、=h/2*(0+2*t+f(1);T=vpa(T,7) %梯形公式for i=0:n-1 s1=s1+f(h/2+i*h);endfor i=1:n-1 s2=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)，而其余情况结果都相同）C

8、lear;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:m h=h/2;n=1/h;t=0; for j=1:n-1 t=t+f(j*h); end T(i)=h/2*(0+2*t+f(1);%梯形公式endfor i=1:m-1 for j=m:i+1 T(j)=4i/(4i-1)*T(j)-1/(4i-1)*T(j-1); %通过不断的迭代求得T（j)，即T表的对角线上的元素。 endendvpa(T(m

9、),7)2、运行结果T = -0.4443875 S = -0.4444345 ans = -0.44444143、问题结果a. 利用复合梯形公式及复合辛普森公式求解：步长h10.10.010.001梯形求积T=01.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 -7 0 1 -3 2.099999 6 2 5

10、-1 5 -1 2 1 0 2;b=8;5.900001;5;1;m,n=size(A); L=eye(n); U=zeros(n); flag='ok' for i=1:n U(1,i)=A(1,i); end for r=2:n L(r,1)=A(r,1)/U(1,1); end for i=2:n for j=i:n z=0; for r=1:i-1 z=z+L(i,r)*U(r,j); end U(i,j)=A(i,j)-z; end if abs(U(i,i)<eps flag='failure' return; end for k=i+1:n

11、m=0; for q=1:i-1 m=m+L(k,q)*U(q,i); end L(k,i)=(A(k,i)-m)/U(i,i); end end LUy=Lb;x=UydetA=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); %找列主元所在子矩阵的行rr

12、= r + k - 1; % 列主元所在大矩阵的行if r>k temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)=0, error(对角元出现0), end % 把增广矩阵消元成为上三角for p = k+1:n Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end end % 解上三角方程组A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n

13、:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x(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 0 0.0000 -2.5000 0.0000 0 0.0000 -2.4000 0.0000 0.0000U = 1.0e+007 * 0.0000 -0.0000 0 0.0000 0 -0.0000 0.0000 0.0000 0 0 1.5000 0.5750 0 0 0 0.0000x = -0.0000 -1.0000

14、 1.0000 1.0000detA = -762.00012）列主元消去法detA = 762.0001ans = 0.0000 -1.0000 1.0000 1.00003、 问题结果1） LU分解解线性方程组L=U=x=（0.0000，-1.0000，1.0000，1.0000）TdetA=-762.0012)列主元消去法x=（0.0000，-1.0000，1.0000，1.0000）TdetA=762.001 第六章1、程序代码（1）Jacobi迭代Clear;clc;n = 6; %也可取n=8,10H = hilb(n); b = H * ones(n, 1); e = 0.000

15、01; for i = 1:n if H(i, i)=0 '对角元为零，不能求解' return endendx = zeros(n, 1); k = 0; kend = 10000; r = 1; while k<=kend & r>e x0 = x; for i = 1:n s = 0; for j = 1:i - 1 s = s + H(i, j) * x0(j); end for j = i + 1:n s = s + H(i, j) * x0(j); end x(i) = b(i) / H(i, i) - s / H(i, i); end r =

16、norm(x - x0, inf); k = k + 1; endif k>kend '迭代不收敛，失败' else '求解成功' x k end（2）SOR迭代1）程序代码function s = SOR(n, w); H = hilb(n); b = H*ones(n, 1); e = 0.00001; for i = 1:n if H(i,i)=0 对角线为零，不能求解 return endendx = zeros(n, 1); k = 0;kend = 10000; r = 1; while k<=kend & r>e x0 =

17、 x; for i = 1:n s = 0; for j = 1:i - 1 s = s + H(i, j) * x(j); end for j = i + 1:n s = s + H(i, j) * x0(j); end x(i) = (1 - w) * x0(i) + w / H(i, i) * (b(i) - s); end r = norm(x - x0, inf); k = k + 1; endif k>kend '迭代不收敛，失败' else '求解成功' xend2）从命令窗口中分别输入：SOR（6,1）SOR（8,1）SOR（10,1）SO

18、R（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/(x2+2*x+10);2)建立函数文件：bdd.mfunction y, n = bdd(x, eps) if nargin=1 eps=1.0e-8;elseif nargin<1 error returnend x1 = g(x); n = 1; while (norm(x1-x)>=eps)&&(n<=10000) x = x1; x1 = g(x); n = n + 1; end y = x; n