数值线性代数第二版徐树方高立张平文上机习题第三章实验报告

第1篇一、实验目的1. 掌握线性代数基本概念和基本运算方法。

2. 熟悉MATLAB软件在解决线性代数问题中的应用。

3. 提高实际操作能力和编程能力。

二、实验环境1. 操作系统：Windows 102. 软件环境：MATLAB R2019b3. 实验设备：计算机三、实验内容1. 矩阵的基本运算2. 矩阵的秩3. 矩阵的逆4. 线性方程组的求解5. 特征值和特征向量6. 二次型及其标准形四、实验步骤1. 矩阵的基本运算（1）创建矩阵A：A = [1, 2, 3; 4, 5, 6; 7, 8, 9]（2）计算矩阵A的转置：A_transpose = A'（3）计算矩阵A的行列式：det_A = det(A)（4）计算矩阵A的逆：A_inverse = inv(A)2. 矩阵的秩（1）创建矩阵B：B = [1, 2, 3, 4; 5, 6, 7, 8; 9, 10, 11, 12]（2）计算矩阵B的秩：rank_B = rank(B)3. 矩阵的逆（1）创建矩阵C：C = [1, 2; 3, 4]（2）判断矩阵C是否可逆：is_inverse = rank(C) == size(C, 1)（3）如果可逆，计算矩阵C的逆：C_inverse = inv(C)4. 线性方程组的求解（1）创建矩阵A和B：A = [1, 2; 3, 4]B = [5; 6]（2）使用MATLAB内置函数求解线性方程组：x = A \ B5. 特征值和特征向量（1）创建矩阵D：D = [4, 1; 2, 3]（2）计算矩阵D的特征值和特征向量：[V, D] = eig(D)6. 二次型及其标准形（1）创建矩阵E：E = [2, 1; 1, 3]（2）计算矩阵E的特征值和特征向量：[V, D] = eig(E)（3）将二次型E化为标准形：Q = V D inv(V)五、实验结果与分析1. 矩阵的基本运算（1）矩阵A：1 2 34 5 67 8 9（2）矩阵A的转置：1 4 72 5 83 6 9（3）矩阵A的行列式：（4）矩阵A的逆：-1.5 0.50.5 -0.52. 矩阵的秩矩阵B的秩为2。

数值代数上机实验报告试验项目名称：平方根法与改进平方根法实验内容：先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序，然后用你编写的程序求解对称正定方程组Ax=b，其中，A=[101 10 1…1 10 11 10]100*100b随机生成，比较计算结果，评论方法优劣。

实验要求：平方根法与改进的平方根的解法步骤；存储单元，变量名称说明；系数矩阵与右端项的生成；结果分析。

实验报告姓名：罗胜利班级：信息与计算科学0802 学号：u200810087 实验一、平方根法与改进平方根法先用你所熟悉的计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用你编写的程序求解对称正定方程组AX=b，其中系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为=，向量b的第i个分量为=.平方根法函数程序如下：function [x,b]=pingfanggenfa(A,b)n=size(A);n=n(1);x=A^-1*b; %矩阵求解disp('Matlab自带解即为x');for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend %Cholesky分解for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('平方根法的解即为b');endfunction [x]=ave(A,b,n) %用改进平方根法求解Ax=b L=zeros(n,n); %L为n*n矩阵D=diag(n,0); %D为n*n的主对角矩阵S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;for i=1:nfor j=1:n %验证A是否为对称正定矩阵if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时，输出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %将A分解使得A=LDL Tfor i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1); % x，y为n*1阶矩阵x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); %通过LDy=b解得y的值endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n)); %通过L T x=y解得x的值改进平方根法函数程序如下：function b=gaijinpinfanggenfa(A,b)n=size(A);n=n(1);v=zeros(n,1);for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);end %LDL'分解B=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;EndA=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n); %前代法A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1); %回代法disp('改进平方根法解得的解即为b');end调用函数解题：clear;clc;n=input('请输入矩阵维数：');b=zeros(n,1);A=zeros(n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);b(i)=b(i)+1/(i+j-1);endend %生成hilbert矩阵[x,b]=pingfanggenfa(A,b) b=gaijinpinfanggenfa(A,b)运行结果：请输入矩阵维数：40Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.570692e-020. > In pingfanggenfa at 4In qiujie at 10Matlab自带解即为x平方根法的解即为bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000-12.00002.000010.2500-10.5000-1.0000-10.875083.000046.0000-98.0000-69.000068.000021.0000-50.7188-8.7500-8.0000 112.0000 6.0000 -68.7500 22.000044.0000 -28.0000 8.0000 -44.000012.0000b =1.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.32850.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.0031-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改进平方根法解得的解即为bb =1.0e+024 *0.0000-0.00000.0001-0.0012-0.0954 0.4208 -1.2101 2.0624 -1.0394 -3.3343 6.2567 -0.2463 -7.45942.80303.6990 0.7277 -1.7484 -0.4854 -3.6010 0.2532 5.1862 1.4410 0.8738 -4.5654 1.0422 4.0920 -2.7764 -2.2148 -0.8953 0.3665 4.8967 1.0416 0.1281-1.1902-2.83348.4610-3.6008实验二、利用QR分解解线性方程组：利用QR分解解线性方程组Ax=b，其中A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];b=[32 26 38 30];求解程序如下：定义house函数：function [v,B]=house(x)n=length(x);y=norm(x,inf);x=x/y;Q=x(2:n)'*x(2:n);v(1)=1;v(2:n)=x(2:n);if n==1B=0;elsea=sqrt(x(1)^2+Q);if x(1)<=0v(1)=x(1)-a;elsev(1)=-Q/(x(1)+a);endB=2*v(1)^2/(Q+v(1)^2);endend进行QR分解：clear;clc;A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12]; b=[32 26 38 30];b=b';x=size(A);m=x(1);n=x(2);d=zeros(n,1);for j=1:n[v,B]=house(A(j:m,j));A(j:m,j:n)=(eye(m-j+1)-B*(v')*v)*A(j:m,j:n); d(j)=B;if j<m< p="">A(j+1:m,j)=v(2:m-j+1);endend %QR分解R=triu(A); %得到R D=A;I=eye(m,n);Q=I;for i=1:nD(i,i)=1;endH=tril(D);M=H';for i=1:nN=I-d(i)*H(1:m,i)*M(i,1:m);Q=Q*N;end %得到Qb=(Q')*b; %Q是正交阵for j=n:-1:2b(j)=b(j)/R(j,j);b(1:j-1)=b(1:j-1)-b(j)*R(1:j-1,j);endb(1)=b(1)/R(1,1); %回带法运行结果如下：R =18.7617 9.8072 15.7769 11.08640 9.9909 9.3358 7.53410 0 5.9945 9.80130 0 0 -0.5126Q =0.8528 -0.4368 -0.2297 -0.17090.2132 0.7916 -0.4594 -0.34170.4264 0.3822 0.2844 0.76890.2132 0.1911 0.8095 -0.5126b=1.000000000000001.000000000000010.9999999999999881.00000000000001实验三、Newton下山法解非线性方程组：3x-cos（yz）-=0，-81+sinz+1.06=0，exp（-xy）+20z+=0；要求满足数值解=满足或.定义所求方程组的函数：Newtonfun.mfunction F = Newtonfun(X)F(1,1)=3*X(1)-cos(X(2)*X(3))-1/2;F(2,1)=X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;F(3,1)=exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3;End向量求导：Xiangliangqiudao.mfunction J=xiangliangqiudao()syms x y zX=[x,y,z];F=[3*X(1)-cos(X(2)*X(3))-1/2;X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3];J=jacobian(F,[x y z]);End代值函数：Jacobi.mfunction F=Jacobi(x)F=[ 3,x(3)*sin(x(2)*x(3)), x(2)*sin(x(2)*x(3));2*x(1), -162*x(2)-81/5,cos(x(3));-x(2)/exp(x(1)*x(2)),-x(1)/exp(x(1)*x(2)),20];End方程组求解：format long; %数据表示为双精度型X1=[0,0,0]';eps=10^(-8);k=1;i=1;X2=X1-Jacobi(X1)^(-1)*Newtonfun(X1);while (norm(subs(X2-X1,pi,3.1415926),2)>=eps)&&(norm(Newtonfun(X1),2)>=eps) if norm(Newtonfun(X2),2)<="" p="">X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-B*C;i=i+1;elsev=1/(2^k); %引入下山因子X1=X2;B=inv(Jacobi(X2));C=Newtonfun(X2);X2=X2-v*B*C;k=k+1;endendj=i+k-1 %迭代次数X=X2 %输出结果运行结果如下：j =5X =0.500000000000000 -0.000000000000000 -0.523598775598299</m<>。

上机习题1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss 消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对Gauss 消去法的看法。

Sol ：（1）先用matlab 将不选主元和列主元Gauss 消去法编写成通用的子程序，得到P U L ,,： 不选主元Gauss 消去法：[])(,A GaussLA U L =得到U L ,满足LU A =列主元Gauss 消去法：[])(,,A GaussCol P U L =得到P U L ,,满足LU PA =（2）用前代法解()Pb or b Ly =，得y用回代法解y Ux =，得x求解程序为()P U L b A Gauss x ,,,,=（P 可缺省，缺省时默认为单位矩阵）（3）计算脚本为ex1_1代码%算法（计算三角分解：Gauss 消去法）function [L,U]=GaussLA(A)n=length(A);for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endU=triu(A);L=tril(A);L=L-diag(diag(L))+diag(ones(1,n));end%算法计算列主元三角分解：列主元Gauss消去法）function [L,U,P]=GaussCol(A)n=length(A);for k=1:n-1[s,t]=max(abs(A(k:n,k)));p=t+k-1;temp=A(k,1:n);A(k,1:n)=A(p,1:n);A(p,1:n)=temp;u(k)=p;if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); elsebreak;endendL=tril(A);U=triu(A);L=L-diag(diag(L))+diag(ones(1,n));P=eye(n);for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);P(u(i),:)=temp;endend%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin<5P=eye(length(A));endn=length(A);b=P*b;for j=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j); endb(n)=b(n)/L(n,n);y=b;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);x=y;endex1_1clc;clear;%第一题A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);b=[7;15*ones(82,1);14];%不选主元Gauss消去法[L,U]=GaussLA(A);x1_1=Gauss(A,b,L,U);%列主元Gauss消去法[L,U,P]=GaussCol(A);x1_2=Gauss(A,b,L,U,P);%解的比较subplot(1,3,1);plot(1:84,x1_1,'o-');title('Gauss');subplot(1,3,2);plot(1:84,x1_2,'.-');title('PGauss');subplot(1,3,3);plot(1:84,ones(1,84),'*-');title('精确解');结果为（其中Gauss表示不选主元的Gauss消去法，PGauss表示列主元Gauss 消去法，精确解为[]'⨯8411,,1 ）：-6-4-202468Gauss050100PGauss 00.20.40.60.811.21.41.61.82精确解由图，显然列主元消去法与精确解更为接近，不选主元的Gauss 消去法误差比列主元消去法大，且不如列主元消去法稳定。

15.解： 解
(α1 , α 2 , α 3 )T = ( β1 , β 2 , β 3 )
AT = B
特别提示
−1 1 0 T = A−1 B = 2 − 1 2 0 1 − 1
A−1 ( AM E ) → ( E M A−1 ) (1) −1 A B (2)( AM E ) → ( E M A−1 B)
a 1 (已知 α 4 = aα1 + bα 2 ⇒ α 2 = − α1 + α 4 已知) 已知 b b ac c α 5 = − α1 + dα 3 + α 4 b b R(α1 , α 3 , α 4 , α 5 ) < 4
再令 x1α1 + x2α 3 + x3α 4 = ϑ
x1α1 + x2α 3 + x3 (aα1 + bα 2 ) = ϑ ( x1 + ax3 )α1 + bx3α 2 + x2α 3 = ϑ
ε 1 , ε 2 Lε n 能由 α1 , α 2 Lα n 线性表出 定理） （定理） α1 , α 2 Lα n 能由 ε 1 , ε 2 Lε n 线性表出
α1 , α 2 Lα n 与 ε 1 , ε 2 Lε n 等价
R(α1 , α 2 Lα n ) = R(ε 1 , ε 2 Lε n ) = n
第三章
向量空间习题答案
1.设 v = (1,−1,1)T , v = (2,1,3)T , v = (2,1,3)T , 设 1 2 3 求 v1 − v2 , 及 3v1 − 2v2 + v3 . 解：
v1 − v2 = (1,−2,−1)T

第四章上机习题1考虑两点边值问题⎪⎩⎪⎨⎧==<<=+.1)1(,0)0(10 ,22y y a a dx dy dx y d ε 轻易知道它正确解为ax e e ay x +---=--)1(111εε为了把微分方程离散化, 把[0,1]区间n 等分, 令h=1/n,1,,1,-==n i ih x i得到差分方程,21211a hy y h y y y i i i i i =-++-++-ε简化为 ,)2()(211ah y y h y h i i i =++-+-+εεε从而离散化后得到线性方程组系数矩阵为⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡+-++-++-++-=)2()2()2()2(h h h h h h h A εεεεεεεεεε 对,100,2/1,1===n a ε分别用Jacobi 迭代法, G-S 迭代法和SOR 迭代法求线性方程组解, 要求有4位有效数字, 然后比较与正确解得误差。

对,0001.0,01.0,1.0===εεε考虑一样问题。

解 （1）给出算法:为解b Ax =, 令U L D A --=, 其中][ij a A =, ),,,(2211nn a a a diag D = ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=-00001,21323121n n n n a a a a a a L，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=-0000,122311312 n n n n a a a a a a U 利用Jacobi 迭代法, G-S 迭代法, SOR 迭代法解线性方程组, 均能够下步骤求解: step1给定初始向量x0=(0,0,...,0), 最大迭代次数N, 精度要求c, 令k=1 step2令x=B*x0+gstep3若||x-x0||2<c, 算法停止, 输出解和迭代次数k, 不然, 转step4step4若k>=N,算法停止, 迭代失败, 不然, 令x0=x, 转step2在Jacobi 迭代法中, B=D -1*(L+U),g=D -1*b在G-S 迭代法中, B=D -1*(L+U),g=D -1*b在SOR 迭代法中, B=(D-w*L)-1*[(1-w)*D+w*U],g=w*(D-w*L)-1*b另外, 在SOR 迭代法中, 上面算法step1中要给定松弛因子w, 其中0<w<2 为计算结果, 要求w=0.5。

数值代数课程设计实验报告姓名： 班级： 学号： 实验日期：一、实验名称 代数的数值解法 二、实验环境 MATLAB7.0实验一、平方根法与改进平方根法一、实验要求：用熟悉的计算机语言将不选主元和列主元Gasuss 消元法编写成通用的子程序，然后用编写的程序求解下列方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--⨯14151515157681681681681681612321n n n n n x x x x x x 用所编的程序分别求解40、84、120阶方程组的解。

二、算法描述及实验步骤GAuss 程序如下：（1）求A 的三角分解：LU A =;（2）求解b y =L 得y ； （3）求解y x =U 得x ；列主元Gasuss 消元法程序如下： 1求A 的列主元分解：LU PA =;2求解b y P L =得y ； 3求解y x =U 得x ；三、调试过程及实验结果：%----------------方程系数---------------->> A1=Sanduijiaozhen(8,6,1,40); >> A2=Sanduijiaozhen(8,6,1,84); >> A3=Sanduijiaozhen(8,6,1,120); >> b1(1)=7;b2(1)=7;b3(1)=7;>> for i=2:39b1(i)=15;end>> b1(40)=14;>> for i=2:83b2(i)=15;end>> b2(40)=14;>> for i=2:119b1(i)=15;end>> b3(120)=14;%----------------方程解---------------->> x11=GAuss(A1,b1')>> x12=GAuss Zhu(A1,b1')>> x21=GAuss(A2,b2')>> x22=GAuss Zhu(A3,b3')>> x31=GAuss(A3,b3')>> x32=GAuss Zhu(A3,b3')运行结果：（n=40）GAuss消元法的解即为x11 =1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000列主元GAuss消元法的解即为x12 =1 1 1 1 1 1 1 1 1 1 111111111111111111111111111111六、源程序：function A=Sanduijiaozhen(a,b,c,n)%生成n阶以a,b,c为元素的三对角阵A=diag(b*ones(1,n),0)+diag(c*ones(1,n-1),1)+diag(a*ones(1,n-1),-1);function x=GAuss(A,b)n=length(b);x=b;%-------分解---------------for i=1:n-1for j=i+1:nmi=A(j,i)/A(i,i);b(j)=b(j)-mi*b(i);for k=i:nA(j,k)=A(j,k)-mi*A(i,k);endAB=[A,b];endend%-----------回代------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1s=0;for j=i+1:ns=s+A(i,j)*x(j);endx(i)=(b(i)-s)/A(i,i);endfunction x=GAussZhu(A,b)n=length(b);x=b;%----------------------选主元---------------------for k=1:n-1a_max=0;for i=k:nif abs(A(i,k))>a_maxa_max=abs(A(i,k));r=i;endendif r>kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%--------------消元-----------------for i=k+1:nm=A(i,k)/A(k,k);for j=k:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);endendif abs(A(n,n))==0return;endAbZhu=[A,b];%----------------回代-----------------------x(n)=b(n)/A(n,n);for i=n-1:-1:1for j=i+1:nb(i)=b(i)-A(i,j)*x(j);endx(i)=b(i)/A(i,i);end实验二、平方根法与改进平方根法一、实验要求：用计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组100阶方程组AX=b，二、算法描述及实验步骤：平方根法函数程序如下：1、求A 的Cholesky 分解：L L A T=;2、求解b y =L 得y ；3、求解y x =TL 得x ； 改进平方根法函数程序如下：1、求A 的Cholesky 分解：T=LDL A ; 2、求解b y =L 得y ； 3、求解y x =TDL 得x ；三、调试过程及实验结果：clear;clc;%----------------方程系数---------------->> A=Sanduijiaozhen(1,10,1,100); >> b(1)=11; >> for i=2:99 b(i)=12; end>> b(100)=11;>> x1=Cholesky(A,b') >> x2=GJCholesky(A,b')运行结果：平方根法的解即为 x1 =1.0000 1.00001.0000 1.00001.0000 1.00001.0000 1.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000改进平方根法解得的解即为x2 =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00000.99991.00090.99081.09080.1010四、源程序：function x=Cholesky(A,b)n=size(A);n=n(1);% x=A^-1*b;% disp('Matlab自带解即为x');%-----------------Cholesky分解-------------------for k=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:n;A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endend%------------------前代法求解Ly=b----------------------------for j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法求解L'x=y-----------------------------A=A';for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('平方根法的解即为');function b=GJCholesky(A,b)n=size(A);n=n(1);v=zeros(n,1);%----------------------LDL'分解-----------------------------for j=1:nfor i=1:j-1v(i)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);endB=diag(A);D=zeros(n);for i=1:nD(i,i)=B(i);A(i,i)=1;end%-------------------前代法---------------------------A=tril(A); %得到L和Dfor j=1:n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);%-----------------回代法-----------------------------A=D*(A');for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);endb(1)=b(1)/A(1,1);disp('改进平方根法解得的解即为');实验三、二次多项式拟合一、实验要求：用计算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式使在残向量的范数最小的意义下拟合下面的数据t-1 -0.75 -0.5 0 0.25 0.5 0.75iy 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125i二、算法描述及实验步骤：QR分解求解程序如下：1、求A 的QR 分解;2、计算b c 11T =Q ;3、求解上三角方程1c x =R 得x ；三、调试过程及实验结果：>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图');运行后屏幕显示数据的散点图（略）.（3）编写下列MATLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序 >> syms a b c>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; >> fi=a.*t.^2+ b.*t+c%运行后屏幕显示关于 ,,a b c 的线性方程组fi =[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b +c]编写构造残向量2范数的MATLAB 程序>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> fy=fi-y; fy2=fy.^2; J=sum(fy.^2); 运行后屏幕显示误差平方和如下 J=(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2为求,,a b c 使J 达到最小，只需利用极值的必要条件0J a ∂=∂，0J b ∂=∂，0J c∂=∂，得到关于,,a b c 的线性方程组，这可以由下面的MATLAB 程序完成，即输入程序 >> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3) 运行后屏幕显示J 分别对,,a b c 的偏导数如下 Ja11 =451/128*a-63/32*b+43/8*c-887/128 Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8解线性方程组112131000Ja Ja Ja ===，，，输入下列程序 >> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8]; >> C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下 C =0.3081 0.8587 1.4018 f =924/2999*x^2+10301/11996*x+4204/2999 故所求的拟合曲线为2()0.30810.8581 1.4018f x x x =++四、源程序：>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75];>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> plot(t,y,'r*');>> legend('实验数据(ti,yi)'); >> xlabel('t'), ylabel('y');>> title('二次多项式拟合的数据点(ti,yi)的散点图'); >> syms a b c>> t=[-1 -0.75 -0.5 0 0.25 0.5 0.75]; >> fi=a.*t.^2+ b.*t+c fi =[ a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c]>> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> y=[ 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125]; >> fy=fi-y; fy2=fy.^2; J=sum(fy.^2) J =(a-b+c-1)^2+(9/16*a-3/4*b+c-13/16)^2+(1/4*a-1/2*b+c-3/4)^2+(c-1)^2+(1/16*a+1/4*b+c-21/16)^2+(1/4*a+1/2*b+c-7/4)^2+(9/16*a+3/4*b+c-37/16)^2>> Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c);>> Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3)Ja11 =451/128*a-63/32*b+43/8*c-887/128Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*c-143/8>> A=[451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14]; >> B=[887/128,61/32,143/8];>> C=B/A, f=poly2sym(C)C =0.3081 0.8587 1.4018f =924/2999*x^2+10301/11996*x+4204/2999>>。

3. 求下列排列的逆序数
(1) (315624) 6
(2) (13(2n 1)24(2n)) n(n 1)
2
4. 计算下列行列式
2500 350
55
(1)
500 70
1500 70
31
35000(5 15) 350000
a11 a12 a21 a22 (2) 0 0
00 00 a33 a34
2 0
0
0 1 0
0
0 1
4
4 0 0
0 3 0
0 0 2
2 0 0
0 3 0
0
0 1
2
12.设
1 2 3
A
0 3 0
2 2 1
2 0 2
利用初等行变换求 A1
2
1 1 1
,
1 2 3 2 1 0 0 0
(
A
E)
0 3 0
2 2 1
2 0 2
1 0 1 0 0
5、已知两个线性变换
x2
x1 2 y1 y3 2 y1 3y2
2 y3
x3 4 y1 y2 5 y3
y1 3z1 z2 y2 2z1 z3
,
y3 z2 3z3
求从 z1, z2 , z3 到 x1, x2 , x3 的 线性变换
分析：X AY ,Y BZ ,
1 1
0 0
0 0
1 0
10
1 2 3 2 1 0 0 0
~r3 3r1 0
0 0
2 4 1
2 9 2
1 0 1 0 0
5 1
3 0
0 0
1 0
0 1
1 2 3 2 1 0 0 0

“数值计算方法”上机实验指导书实验一 误差分析实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。

对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。

通过本实验可获得一个初步体会。

数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。

病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。

问题提出：考虑一个高次的代数多项式)1.1()()20()2)(1()(201∏=−=−−−=k k x x x x x p显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。

现考虑该多项式的一个扰动)2.1(0)(19=+x x p ε其中ε是一个非常小的数。

这相当于是对（1.1）中19x 的系数作一个小的扰动。

我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。

实验内容：为了实现方便，我们先介绍两个MATLAB 函数：“roots ”和“poly ”。

roots(a)u =其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。

设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程01121=+++++−n n n n a x a x a x a的全部根；而函数 poly(v)b =的输出b 是一个n+1维向量，它是以n 维向量v 的各分量为根的多项式的系数（从高到低排列）。

可见“roots ”和“poly ”是两个互逆的运算函数。

))20:1((;)2();21,1(;000000001.0ve poly roots ess ve zeros ve ess +===上述简单的MATLAB 程序便得到（1.2）的全部根，程序中的“ess ”即是（1.2）中的ε。

@上机习题1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss 消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对Gauss 消去法的看法。

Sol ：（1）先用matlab 将不选主元和列主元Gauss 消去法编写成通用的子程序，得到P U L ,,： 不选主元Gauss 消去法：[])(,A GaussLA U L =得到U L ,满足LU A = 列主元Gauss 消去法：[])(,,A GaussCol P U L =得到P U L ,,满足LU PA = （2）用前代法解()Pb or b Ly =，得y用回代法解y Ux =，得x]求解程序为()P U L b A Gauss x ,,,,=（P 可缺省，缺省时默认为单位矩阵）（3）计算脚本为ex1_1 代码%算法（计算三角分解：Gauss 消去法） function [L,U]=GaussLA(A) n=length(A);—for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); endU=triu(A);L=tril(A);L=L-diag(diag(L))+diag(ones(1,n));end！%算法计算列主元三角分解：列主元Gauss消去法）function [L,U,P]=GaussCol(A)n=length(A);for k=1:n-1[s,t]=max(abs(A(k:n,k)));p=t+k-1;temp=A(k,1:n);￥A(k,1:n)=A(p,1:n);A(p,1:n)=temp;u(k)=p;if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); elsebreak;^endendL=tril(A);U=triu(A);L=L-diag(diag(L))+diag(ones(1,n)); P=eye(n);for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);{P(u(i),:)=temp;endend%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin<5P=eye(length(A));￥endn=length(A);b=P*b;for j=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);<y=b;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);x=y;|endex1_1clc;clear;%第一题A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1); b=[7;15*ones(82,1);14];%不选主元Gauss消去法)[L,U]=GaussLA(A);x1_1=Gauss(A,b,L,U);%列主元Gauss消去法[L,U,P]=GaussCol(A);x1_2=Gauss(A,b,L,U,P);%解的比较subplot(1,3,1);plot(1:84,x1_1,'o-');title('Gauss'); subplot(1,3,2);plot(1:84,x1_2,'.-');title('PGauss');（subplot(1,3,3);plot(1:84,ones(1,84),'*-');title('精确解');结果为（其中Gauss 表示不选主元的Gauss 消去法，PGauss 表示列主元Gauss消去法，精确解为[]'⨯8411,,1 ）：8Gauss50100PGauss精确解由图，显然列主元消去法与精确解更为接近，不选主元的Gauss 消去法误差比列主元消去法大，且不如列主元消去法稳定。

2 1 , k 0,1, 2...n. 这里唯一需要注意的就是 n 1 Pn ( xk ) Pn'1 ( xk )
还需要一个多项式微分运算（使用的是polyder函数） ，求解出这两项之后，剩下的 就是 I n ( f ) 2.
A
k 0
n
k
f ( xk ) 。
Romberg积分公式计算 函数主体是Romberg.m文件,以下只分析这个函数文件的结构: 第一部分是赋初值,对于积分表T还赋值T(1,1):
实现方法说明
程序语言采用Matlab语言，运行环境为Matlab R2009b。
注：主程序文件是First.m和Second.m，其他都是对应的自编函数文件。
1.
复化Guass-Legendre求积公式运算 这个程序主体是First.m文件，在里面定义了变量x和函数f f=sym(100*sin(10*(2+x)^(-1))*(2+x).^(-2));由于Legendre多项式的使用条件是[-1,1]， 因此对积分区间做了变换，函数体本身也就变成了f（x）=(100*sin(10/(x + 2)))/(x + 2)^2。其他都是调用自己编写的函数计算。 ① Legendre多项式生成函数：Legendre.m
2/5
水利系
2008010249
程国安
这个函数主要用于生成 Legendre 多项式，比较简短，最重要的就是一个迭代 关系 P=((2*n-1)*x*Legendre(n-1)-(n-1)*Legendre(n-2))/(n);最后输出一个多项式， 如果需要整合，另外加一句 expand（P）即可。 ② Guass-Legendre积分函数：Guass-Legendre.m 这个函数的相关输入输出说明见函数文件注释。 积分实现主要需要求出 f ( xk ) 和 Ak ， 对于 f ( xk ) ， 重点就是建立内联函数 f ( x) 和 变量序列 xk ，变量序列 xk 主要使用利用Legendre求出对应次数的Legendre多项式 的零点（使用sym2poly求出系数，之后利用roots求零点） 。 对于 Ak ，利用 Ak

（1）估计5到20阶Hilbert 矩阵的∞范数条件数（2）设n n R A ⨯∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡------=111111111011001ΛΛO O MM M O OΛ,先随机地选取n R x ∈,并计算出x A b n =;然后再用列主元Gauss 消去法求解该方程组,假定计算解为∧x 。

试对n 从5到30估计计算解∧x 的精度,并且与真实相对误差作比较。

解(1)分析:利用for 使n 从5循环到20,利用()hilb 函数得到Hilbert 矩阵A ;先将算法2、5、1编制成通用的子程序,利用算法2、5、1编成的子程序)(B opt v =,对TAB -=求解,得到∞-1A的一个估计值v v =~;再利用inf),(A norm 得到∞A ;则条件数inf),(1A norm v A A K *==∞∞-。

另,矩阵A 的∞范数条件数可由inf),(A cond 直接算出,两者可进行比较。

程序为1 算法2、5、1编成的子程序)(B opt v =function v=opt(B)k=1;n=length(B); x=1、/n*ones(n,1);while k==1 w=B*x;v=sign(w); z=B'*v;if norm(z,inf)<=z'*x v=norm(w,1); k=0; elsex=zeros(n,1);[s,t]=max(abs(z)); x(t)=1; k=1; end end end2 问题(1)求解 ex2_1for n=5:20A=hilb(n);B=inv(A、');v=opt(B);K1=v*norm(A,inf);K2=cond(A,inf);disp(['n=',num2str(n)])disp(['估计条件数为',num2str(K1)])disp(['实际条件数为',num2str(K2)])end计算结果为n=5估计条件数为943656实际条件数为943656n=6估计条件数为29070279、0028实际条件数为29070279、0028n=7估计条件数为985194887、5079实际条件数为985194887、5079n=8估计条件数为33872789099、7717实际条件数为33872789099、7717n=9估计条件数为16、422实际条件数为16、422n=10估计条件数为35353368771750、67实际条件数为35353368771750、67n=11估计条件数为1232433965549344实际条件数为1232433965549344Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、547634e-17、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、547634e-17、> In cond at 47In ex2_1 at 6n=12估计条件数为3、9245e+16实际条件数为3、9245e+16Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 7、847381e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 7、847381e-19、> In cond at 47In ex2_1 at 6n=13估计条件数为1、2727e+18实际条件数为1、2727e+18Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、246123e-18、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 2、246123e-18、> In cond at 47In ex2_1 at 6n=14估计条件数为4、8374e+17实际条件数为4、8374e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 8、491876e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 8、491876e-19、> In cond at 47In ex2_1 at 6n=15估计条件数为4、6331e+17实际条件数为5、234289848563619e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 9、137489e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 9、137489e-19、> In cond at 47In ex2_1 at 6n=16估计条件数为8、3166e+17实际条件数为8、3167e+17Warning: Matrix is close to singular or badly scaled、Results may be inaccurate、RCOND = 6、244518e-19、> In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 6、244518e-19、 > In cond at 47 In ex2_1 at 6 n=17估计条件数为1、43e+18 实际条件数为1、43e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、693737e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、693737e-19、 > In cond at 47 In ex2_1 at 6 n=18估计条件数为2、5551e+18 实际条件数为2、8893e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、264685e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 4、264685e-19、 > In cond at 47 In ex2_1 at 6 n=19估计条件数为2、411858563109357e+18 实际条件数为2、411858563109357e+18Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 1、351364e-19、 > In ex2_1 at 3Warning: Matrix is close to singular or badly scaled 、 Results may be inaccurate 、 RCOND = 1、351364e-19、 > In cond at 47 In ex2_1 at 6 n=20估计条件数为2、31633670586674e+18 实际条件数为6、37335273308473e+18结果分析随着矩阵阶数增加,估计值误差开始出现,20,17,16,15 n 时估计条件数与实际值存在误差;且条件数很大,Hilbert 矩阵为病态的。

>> diff(x^4+exp(x),'x',4)
ans =24+exp(x)
5) 求 .
>> syms x
>> diff(asin(sqrt(1-x^2)))
ans =-1/(1-x^2)^(1/2)*x/(x^2)^(1/2)
可得dy=-1/(1-x^2)^(1/2)*x/(x^2)^(1/2)*dx
（签名）
实验目的
1.掌握求数值微分的命令．
2.会求函数零点．
实验内容
数值积分计算，函数零点的计算
主要命令和
程序清单
数值积分：quad,quadl,dblquad,triplequad
求多项式所有根：roots
由根得多项式: poly
实验过程及
结果记录
实验过程及
结果记录
1、求下列积分的数值解
（1）
输入：
syms x
taylor(exp(x),x,5,0)
x=1
1+x+1/2*x^2+1/6*x^3+1/24*x^4
结果：
ans =1+x+1/2*x^2+1/6*x^3+1/24*x^4
ans =2.7083
思考
及
习题
收获
感想
实验报告3.3多元微积分实验
课程名称
实验名称
实验教室
实验日期
班级
学生姓名
实验成绩
2.掌握级数求和、求收敛半径．
3.掌握将函数展开成泰勒级数．
实验内容
微分方程求解（解析解，数值解）
级数收敛性的判断，求收敛区间，和函数；函数的泰勒展开；求近似值

（2） f

fs N

4096 4096
1Hz
（3）直接用 DFT 计算，所需要的复乘次数为
M d (300 200 1)N 101 4096 413696
若用按时间抽取 FFT 则需要的复乘次数为
MF

N 2
log10
N

2048 12

24576
3.13 下面是三个不同的信号 xi (n) ，每个信号均为两个正弦信号的和：
DFT 在加窗后会有两个可区分的谱峰？
解：利用
64
点
DFT
来估计信号谱时，其频率分辨率为


Hale Waihona Puke 2 64信号 x1 (n) cos( n / 4) cos(17 n / 64) 的两个余弦信号的频率间隔为：
1

17 64

4

64

2 64
故利用 64 点 DFT 来估计信号谱时，不能分辨 x1(n) 中两个正弦信号的谱峰。

2M
时，DIF-FFT
共需
M
级分解，每级运算要计算的碟形运算有
N 2
个。
3.4 考虑图 T3-1 中的蝶形。这个蝶形是从实现某种 FFT 算法的信号流图中取出的。从下述论述中选择出最 准确的一个：
（1）这个蝶形是从一个按时间抽取的 FFT 算法中取出的。 （2）这个蝶形是从一个按频率抽取的 FFT 算法中取出的。 （3）由图无法判断该蝶形取自何种 FFT 算法。

N 2
1
[x(n)
n0

x(n

N 2
)]WNnr/

aks
1, 2 ,, s 是分别在1,2 ,, s 的 k 个分量后任意添加 m 个分量 b1 j , b2 j ,, bmj
( j 1,2,, s) 所组成的 k m 维向量，证明：
(1) 若1,2 ,, s 线性无关，则 1, 2 ,, s 线性无关； (2) 若 1, 2 ,, s 线性相关，则1,2 ,, s 线性相关.
（1） 1 (6,4,1,9,2), 2 (1,0,2,3,4), 3 (1,4,9,6,22), 4 (7,1,0,1,3);
（2）1 (1,1,2,4), 2 (0,3,1,2), 3 (3,0,7,14), 4 (2,1,5,6) ,5 (1,1,2,0) ;
（3）1 (1,1,1), 2 (1,1,0), 3 (1,0,0), 4 (1,2,3).
必要性方法1设线性无关证明线性无假设线性相关则中至少有一向量可由其余两个向量线性表示不妨设可由线性表示则向量组可由线性表示且所以线性相关与线性无关矛方法2令设存在k1k2k3使得因为线性无关所以所以线性无关
第三章 课后习题及解答
将 1，2 题中的向量 表示成1,2 ,3,4 的线性组合：
1． 1,2,1,1T ,1 1,1,1,1T ,2 1,1,1,1T ,3 1,1,1,1T ,4 1,1,1,1T. 2． 0,0,0,1,1 1,1,0,1,2 2,1,3,1,3 1,1,0,0,4 0,1,1,1.
0 1 1 101 因为 1,2,3 线性无关，且 1 1 0 2 0 ，可得 1 2,2 3,3 1的秩为 3 011 所以1 2 ,2 3,3 1 线性无关.线性无关；反之也成立.
方法 2，充分性，设1,2 ,3 线性无关，证明1 2 ,2 3,3 1 线性无关.

L1 ( x)
L2 ( x)
( x 0)( x 0.6) (100*x*(x - 3/5))/27 (0.9 0)(0.9 0.6)
30)*(x 9/10))/27 + -
P2 ( x) L0 ( x) cos 0 L1 ( x) cos 0.6 L2 ( x) cos 0.9 =((50*x
可以预测 1930，1965，2010 年的人口分别是 169649，1.9177e+005，171351
EXERCISE SET 3.2 4、
P131
a) 根据 Algorithm 3.2，利用课本作者网站上的关于本书的 MATLAB 程序 ALG032.M 运行该程序，在命令行窗口出现如下： Warning: Could not find an exact (case-sensitive) match for 'ALG032'. E:\ 个人 \ 工作 \ 高教数值分析 \ 第八版英文程序 \Matlab-Programs\matlab\m1\ALG032.M is a case-insensitive match and willbe used instead. You can improve the performance of your code by using exact name matches and we therefore recommend that you update your usage accordingly. Alternatively, you can disable this warning using warning('off','MATLAB:dispatcher:InexactCaseMatch'). This warning will become an error in future releases. Newtons form of the interpolation polynomial Choice of input method: 1. Input entry by entry from keyboard 2. Input data from a text file 3. Generate data using a function F Choose 1, 2, or 3 please 1 Input n 4 Input X(0) and F(X(0)) on separate lines 0.0 -6.00000 Input X(1) and F(X(1)) on separate lines 0.1 -5.89483 Input X(2) and F(X(2)) on separate lines 0.3 -5.65014 Input X(3) and F(X(3)) on separate lines 0.6 -5.17788 Input X(4) and F(X(4)) on separate lines 1.0 -4.28172 Select output destination 1. Screen 2. Text file Enter 1 or 2 1 NEWTONS INTERPOLATION POLYNOMIAL Input data follows: X(0) = 0.00000000 F(X(0)) =

暨南大学本科实验报告专用纸课程名称数值代数成绩评定实验项目名称(填写实验所属章节名称) 指导教师刘娟实验项目编号0701******* 实验项目类型验证实验地点机房学生姓名王伟文学号2010051727学院信息科学技术学院数学系信息与计算科学专业2010级实验时间2011年9月 1 日～12月30日温度24℃第一章1.实验选题：（写出你在该项目所选的实验题目）第一章上机练习12.谈谈你对该算法的理解：（简单谈一下你是如何理解该算法的？）对算法的理解：先将84阶的矩阵A分解为一个下三角矩阵L和上三角矩阵U,先考虑下三角形方程组Ly=b,利用前代法解出y,再考虑上三角形方程组Ux=y;利用回代法解出x。

列主元gauss消去法Ax b PA LU Ly Pb Ux y=⇔===,,;3.实验内容（将实验程序及其实验结果粘贴，最好对程序各部分注释清楚，比如设置了哪些函数，这些函数的输入输出是什么，具有什么功能？）实验程序对矩阵A进行LU分解的程序function [ L,U ] = LUfac( A )for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A,0);for i=1:nL(i,i)=1;endU=triu(A,0);end利用前代法解出y值的程序function [ b ] = TSL( L,b )n=size(L,1);for j=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);end利用回代法解出x值的程序function [ b ] = TSU( U,b )n=size(U,1);for j=n:-1:2b(j)=b(j)/U(j,j);b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);endb(1)=b(1)/U(1,1);end主函数程序（‘生成84阶的矩阵A’）A=eye(84);A=6*A;for i=2:84A(i,i-1)=8;A(i-1,i)=1;end（‘生成84乘1的矩阵b’）b=ones(84,1);b=b*15;b(1)=7;b(84)=14;[L,U]=LUfac(A);（‘调用函数LUfac对矩阵A进行分解’）y=TSL(L,b);（‘调用函数TSL求解Ly=b方程’）x=TSU(U,y);（‘调用函数TSU求解UX=b方程’）用MATLAB求解得到的结果x’ans =1.0e+008 *Columns 1 through 70.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 8 through 140.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 15 through 210.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 22 through 280.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 29 through 350.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 36 through 420.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 43 through 490.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 50 through 560.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 Columns 57 through 63-0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 Columns 64 through 700.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0003 0.0007 Columns 71 through 77-0.0013 0.0026 -0.0052 0.0105 -0.0209 0.0419 -0.0836 Columns 78 through 840.1665 -0.3303 0.6501 -1.2582 2.3487 -4.0263 5.3684列主元gauss消去法function [L,U,P]=Lufac(A)n=size(A,1);P=eye(n,n);for i=1:n-1[r,m]=max(abs(A(i:n,i)));m=m+i-1;A([i,m],:)=A([m,i],:);P([i,m])=P([m,i]);if A(i,i)~=0A(i+1:n,i)=A(i+1:n,i)/ A(i,i);A(i+1:n,i+1:n)= A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n); endendU=triu(A);L=tril(A,-1)+eye(n,n);endA=eye(84);A=6*A;for i=2:84A(i,i-1)=8; A(i-1,i)=1; endb=ones(84,1);b=b*15;b(1)=7;b(84)=14;[L,U,P]=Lufac(A); b=P*b;y=TSL(L,b);x=TSU(U,y);求解结果为x =1.0e+025 *0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00010.0002-0.00040.0007-0.00150.0030-0.00600.0119-0.02380.0475-0.09460.1878-0.36960.7154-1.33552.2894-3.05251.教师评语、评分：（请认真对待每次数值代数实验，期末交实验报告，将计算实验成绩）2.实验选题：第三章上机习题1（3）。

n=10 2 nonuniform interval uniform interval 1.5
1
0.5
0
-0.5 -1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
2.仍然考虑上述实验中的著名问题，使用Matlab的函数“spline”作f（x）的样 条插值。增加插值的节点，观察样条插值的收敛性。


5 2
16
3 (0)
P2 (0) R2 (0)
5.（a）求 f ( x) x 在节点
x1 2, x2 0.5, x3 0, x4 1.5, x5 2
上的三次自然样条插值（即
M1 M 5 0
）。
（b）用同样的数据做Lagrange插值。
将f（x）及它的三次自然样条插值和Lagrange多项式插值用Matlab画出来， 比较它们的结果。 解答：
x_lu=
-6.5000 42.8000 -36.0000
第四章 思考题 1. （a）对给定的连续函数，构造等距节点上的Lagrange插值多项式，节点数 目越多，得到的插值多项式越接近被逼近函数。× （b）对给定的连续函数，构造其三次样条插值，则节点数目越多，得到的 样条函数越接近被逼近的函数。√ （c）高次的Lagrange插值多项式很常用。×
（b）根据迭代收敛条件
1 a 1 2
( B) 1
( B) 2 a 1
1 1 a 2 2
实验题
4.考虑方程组Hx=b，其中系数矩阵为Hilbert矩阵，
H (hi , j ) nn , hi , j 1 , i, j 1, 2,...n i j 1