数值线性代数第二版徐树方高立张平文上机习题第三章实验报告.doc
数值线性代数第二版徐树方高立张平文上机习题第二章实验报告
（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 矩阵为病态的。
数值代数上机实验报告
数值代数上机实验报告试验项目名称：平方根法与改进平方根法实验内容：先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序，然后用你编写的程序求解对称正定方程组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 消去法误差比列主元消去法大，且不如列主元消去法稳定。
线性代数简明教程 (第二版)科学出版社第三章、向量空间x习题答案
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
2021年数值线性代数第二版徐树方高立张平文上机习题实验报告2
第四章上机习题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。
《数值计算》实验报告
《数值计算》实验报告第一部分：简答题（请简要回答以下问题，每小题字数不少于200字）1、Matlab变量命名有什么要求？以下变量名是否合法？对不合法的变量名说明理由。
abcd-2xyz_33chan NaN ABCDefgh2、插值、拟合、回归这三种方法是用来解决什么问题的？面对一组数据，如何选择用什么方法？3、数值积分的主要思想是什么？常用的数值积分公式有哪几个？4、请结合自己的学习，举例说明《数值计算》课程中所学方法在解决实际问题中是如何应用的。
第二部分：基础题（请完成以下问题，要求给出程序语句及计算结果，用截图方式附在各题目下方）1、已知点(1,3.0),(2,3.7),(5,3.9),(6,4.2),(7,5.7),(8,6.6),(10,7.1),(13,6.7), (17,4.5)，绘出经过这些点的函数曲线图形，并给出曲线方程。
答：采取三次样条插值法，九个输入数据分成八段，每一段就是一个三次函数。
这八段的函数形式为y = a0 + a1*x + a2*x^2 + a3*x^3，每个分段函数的参数构成下图所示的coefs 矩阵。
2、在我国某海域测得海洋不同深度处的水温如表1所示，求水深为800m和1500m处的温度。
答：采取线性插值法求得800m和1500m处的温度3、求解方程组⎪⎪⎩⎪⎪⎨⎧=-++=--=-++=++56533332821w z y x w y x w z y x z y x ，请至少使用两种方法求解，并对这两种方法的计算结果进行说明。
高斯消元法LU 分解QR分解Jacobi迭代法使用 Jacobi 迭代法无法求出结果，表示迭代的过程中不收敛4、计算积分dx eI x ⎰-=1022，精度为10-6。
被积函数总共调用 13 次，求得积分值为 0.85565、求方程t et t f t5.0)(sin )(1.02-⋅=-在[0.5,1]内的根。
6、求解微分方程0)1(22=+'--''y y y y ，0)0(,1)0(,300='=≤≤y y x ，绘出解函数的图形。
《数值计算方法》上机实验报告
《数值计算方法》上机实验报告华北电力大学实验名称数值il•算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一.各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程*对于非线性方程，若已知根的一个近似值，将在处展开成一阶xxfx ()0, fx ()xkk泰勒公式"f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2!忽略高次项，有，fxfxfxxx 0 ()()()，，, kkk右端是直线方程，用这个直线方程来近似非线性方程。
将非线性方程的**根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkkfx 0 fx 0 0,解出fX 0 *k XX,, k' fx 0 k水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ikfx ()k 八XX, Ikk* fx()k这就是牛顿迭代公式。
,2,计算机程序框图:，见,,3,输入变量、输出变量说明:X输入变量:迭代初值，迭代精度，迭代最大次数，\0输出变量:当前迭代次数，当前迭代值xkl,4,具体算例及求解结果:2/16华北电力大学实验报吿开始读入l>k/fx()0?,0fx 0 Oxx，，01* fx ()0XX,,，?10kk, ,1，kN, ?xx, 10输出迭代输出X输出奇异标志1失败标志,3,输入变量、输出变量说明: 结束例：导出计算的牛顿迭代公式，并il •算。
（课本P39例2-16） 115cc （0）, 求解结果:10. 75000010.72383710. 72380510. 7238052、列主元素消去法求解线性方程组，1,算法原理:高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角3/16华北电力大学实验报告方程组求解。
线性代数简明教程-第二版-答案
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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第三章上机习题用你所熟悉的的计算机语言编制利用QR分解求解线性方程组和线性最小二乘问题的通用子程序，并用你编制的子程序完成下面的计算任务：（1）求解第一章上机习题中的三个线性方程组，并将所得的计算结果与前面的结果相比较，说明各方法的优劣；（2）求一个二次多项式y=at 2+bt+c ，使得在残向量的 2 范数下最小的意义下拟合表中的数据；表t i -1 0y i 1 1（3）在房产估价的线性模型y x0 a1x1 a2 x2 a11x11中， a1 ,a2 ,, a11分别表示税、浴室数目、占地面积、车库数目、房屋数目、居室数目、房龄、建筑类型、户型及壁炉数目，y 代表房屋价格。
现根据表和表给出的28 组数据，求出模型中参数的最小二乘结果。
（表和表见课本P99-100 ）解分析：（1）计算一个 Householder 变换 H：由于H I 2ww T Ivv T，则计算一个Householder 变换 H 等价于计算相应的、 v 。
其中v x ||x||2e1, 2 /( T )v v 。
在实际计算中，为避免出现两个相近的数出现的情形，当x1 0 时，令v1 - ( x22 x n2 ) ；x1 || x ||2为便于储存，将v 规格化为 v v / v1，相应的，变为2v2 /(v T v)1为防止溢出现象，用x / || x || 代替（2） QR分解：利用 Householder 变换逐步将 A m n , m n 转化为上三角矩阵H n H n 1 H 1 A ，则有RA Q，其中Q H1H 2 H n， R (1: n,:) 。
~在实际计算中，从j 1: n ，若j m ，依次计算x A(( j : m, j )) 对应的( H j)( m k 1) ( m k 1) 即对应的 v j，j，将 v j (2 : m j 1) 储存到 A( j 1: m, j) ，j储存到 d ( j) ，迭代结束后再次计算 Q ，有 H j Ij 1 0H n（ n m 时 Q H 1H 2~ ， Q H1H 2 H n-1 ）0 H j（3）求解线性方程组Ax b 或最小二乘问题的步骤为i计算 A 的QR分解；ii计算 c1Q1T b ，其中 Q1Q (:,1: n)iii利用回代法求解上三角方程组 Rx c1（4）对第一章第一个线性方程组，由于 R 的结果最后一行为零，故使用前代法时不计最后一行，而用运行结果计算 x84。
运算 matlab 程序为1 计算 Householder变换[v,belta]=house(x)function [v,belta]=house(x)n=length(x);x=x/norm(x,inf);sigma=x(2:n)'*x(2:n);v=zeros(n,1);v(2:n,1)=x(2:n);if sigma==0belta=0;elsealpha=sqrt(x(1)^2+sigma);if x(1)<=0v(1)=x(1)-alpha;elsev(1)=-sigma/(x(1)+alpha);endbelta=2*v(1)^2/(sigma+v(1)^2);v=v/v(1,1);endend2 计算A的 QR分解[Q,R]=QRfenjie(A)function [Q,R]=QRfenjie(A)[m,n]=size(A);Q=eye(m);for j=1:nif j<m[v,belta]=house(A(j:m,j));H=eye(m-j+1)-belta*v*v';A(j:m,j:n)=H*A(j:m,j:n);d(j)=belta;A(j+1:m,j)=v(2:m-j+1);endendR=triu(A(1:n,:));for j=1:nif j<mH=eye(m);temp=[1;A(j+1:m,j)];H(j:m,j:m)=H(j:m,j:m)-d(j)*temp*temp';Q=Q*H;endendend3 解下三角形方程组的前代法x=qiandaifa(L,b)function x=qiandaifa(L,b)n=length(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);x=b;end4 求解第一章上机习题中的三个线性方程组ex3_1 clear;clc;%第一题A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1); b=[7;15*ones(82,1);14];n=length(A);%QR分解[Q,R]=QRfenjie(A);c=Q'*b;x1=huidaifa(R(1:n-1,1:n-1),c(1:n-1));x1(n)=c(n)-R(n,1:n-1)*x1;%不选主元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);%解的比较figure(1);subplot(1,3,1);plot(1:n,x1);title( subplot(1,3,2);plot(1:84,x1_1);title( subplot(1,3,3);plot(1:84,x1_2);title( 'QR 分解 ' );'Gauss' );'PGauss' );%第二题第一问A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1));n=length(A);%QR分解tic;[Q,R]=QRfenjie(A);c=Q'*b;x2=huidaifa(R,c);toc;%不选主元Gauss消去法tic;[L,U]=GaussLA(A);x2_1=Gauss(A,b,L,U);toc;%列主元Gauss 消去法tic;[L,U,P]=GaussCol(A);x2_2=Gauss(A,b,L,U,P);toc;%平方根法tic;L=Cholesky(A);x2_3=Gauss(A,b,L,L');toc;%改进的平方根法tic;[L,D]=LDLt(A);x2_4=Gauss(A,b,L,D*L');toc;%解的比较figure(2);subplot(1,5,1);plot(1:n,x2);title( 'QR 分解 ' );subplot(1,5,2);plot(1:n,x2_1);title( 'Gauss' );subplot(1,5,3);plot(1:n,x2_2);title( 'PGauss' );subplot(1,5,4);plot(1:n,x2_3);title( subplot(1,5,5);plot(1:n,x2_4);title( ' 平方根法 ' );' 改进的平方根法' ); %第二题第二问A=hilb(40);b=sum(A);b=b';n=length(A);[Q,R]=QRfenjie(A);c=Q'*b;x3=huidaifa(R,c);%不选主元 Gauss消去法[L,U]=GaussLA(A);x3_1=Gauss(A,b,L,U);%列主元 Gauss 消去法[L,U,P]=GaussCol(A);x3_2=Gauss(A,b,L,U,P);%平方根法L=Cholesky(A);x3_3=Gauss(A,b,L,L');%改进的平方根法[L,D]=LDLt(A);x3_4=Gauss(A,b,L,D*L');%解的比较figure(3);subplot(1,5,1);plot(1:n,x3);title( 'QR 分解 ' );subplot(1,5,2);plot(1:n,x3_1);title( 'Gauss' );subplot(1,5,3);plot(1:n,x3_2);title( 'PGauss' );subplot(1,5,4);plot(1:n,x3_3);title( ' 平方根法 ' ); subplot(1,5,5);plot(1:n,x3_4);title( ' 改进的平方根法 ' );5 求解二次多项式ex3_2clear;clc;t=[-1 0 ];y=[1 1 ];A=ones(7,3);A(:,1)=t'.^2;A(:,2)=t';[Q,R]=QRfenjie(A);Q1=Q(:,1:3);c=Q1'*y';x=huidaifa(R,c)6 求解房产估价的线性模型ex3_3clear;clc;A=xlsread( 专业课数值代数 \' , 'A2:L29' ); y=xlsread( 专业课数值代数 \' , 'M2:M29' ); [Q,R]=QRfenjie(A);Q1=Q(:,1:12);c=Q1'*y;x=huidaifa(R,c);x=x'计算结果为（1）第一章上机习题中的三个线性方程组结果对比图依次为以第二个线性方程组为例，比较各方法的运行速度。
依次为 QR分解，不选主元的Gauss 消去法，列主元 Gauss 消去法，平方根法，改进的平方根法。
Elapsed time is seconds.Elapsed time is seconds.Elapsed time is seconds.Elapsed time is seconds.Elapsed time is seconds.（2）二次多项式的系数为x =（3）房产估价的线性模型的系数为x =Columns 1 through 6Columns 7 through 12结果分析对第一章上机习题中的第二个线性方程组利用五种求解方法求解所需时间可知，不选主元的 Gauss 消去法，列主元 Gauss消去法，改进的平方根法较快，所需时间大致在一个数量级， QR分解，平方根法，所需时间较慢，所需时间在一个数量级上。