数值线性代数实验
线性代数实验报告
2.输入:for
n=20:80 p1(n)=prod(365-n+1:365)/365^n; p(n)=1-p1(n); end plot(p)
输出:
3
3: (1) (2) 输入: R = binornd(20,0.25,3,6) 输出: R= 9 8 3 4 6 6 6 3 4 5 6 2 5 6 6 4 7 4 (3)(4) R = binopdf([0:9],20,0.45) R= 0.0000 0.0001 0.0008 0.0040 0.0746 0.1221 0.1623 0.1771
0.0139
0.0365
4:输入: 1.在单元格 A1 中输入“样本数据” ,在单元格 B4 中输入“指标名称” ,在 单元格 C4 中输入“指标数值” ,并在单元格 A2:A21 中输入样本数据。 2.在单元格 B5 中输入“样本容量” ,在单元格 C5 中输入“20” 。 3.计算样本平均行驶里程。在单元格 B6 中输入“样本均值” ,在单元格 C6 中输入公式: “=AVERAGE(A2,A21), ” 4.计算样本标准差。在单元格 B7 中输入“样本标准差” ,在单元格 C7 中 输入公式: “=STDEV(A2,A21)” ,
4
输出:
5: 输入: R = normrnd(0.5,0.015) load 0.497,0.506 0.518
0.524
0.498
0.511
0.520
0.515
0.512
histfit(0.497 0.506 0.518 0.524 0.498 0.511 0.520 0.515 0.512 ); normplot(0.497 0.506 0.518 0.524 0.498 0.511 0.520 0.515 0.512 ); 输出: R = 0.5066
数值方法迭代法解线性方程组实验报告
附录1:源程序
附录2:实验报告填写说明
1.实验项目名称:要求与实验教学大纲一致。
2.实验目的:目的要明确,要抓住重点,符合实验教学大纲要求。
3.实验原理:简要说明本实验项目所涉及的理论知识。
4.实验环境:实验用的软、硬件环境。
5.实验方案(思路、步骤和方法等):这是实验报告极其重要的内容。
概括整个实验过程。
对于验证性实验,要写明依据何种原理、操作方法进行实验,要写明需要经过哪几个步骤来实现其操作。
对于设计性和综合性实验,在上述内容基础上还应该画出流程图、设计思路和设计方法,再配以相应的文字说明。
对于创新性,还应注明其创新点、特色。
6.实验过程(实验中涉及的记录、数据、分析):写明具体实验方案的具体实施步骤,包括实验过程中的记录、数据和相应的分析。
7.实验结论(结果):根据实验过程中得到的结果,做出结论。
8.实验小结:本次实验心得体会、思考和建议。
9.指导教师评语及成绩:指导教师依据学生的实际报告内容,给出本次实验报告的评价。
数值计算线性方程直接法实验
数值计算线性方程直接法实验实验报告二一、实验目的通过本实验的学习,应掌握各种直接法,如列主元消去法,追赶法等的基本思想和原理,了解它们各自的优缺点及适用范围。
二、实验题目(1)实验一分别用列主元法和顺序高斯消去法求解下面的线性方程组,分析对结果的影响。
--?-2113.6291.52592.1111211314.59103.0164321x x x x78.461297.51 (2) LU 分解的优点实验题目:给定矩阵A 和向量b :--=n n nnn nA 123121,??=1000 b 1)求A 的LU 分解,n 值自己确定;2)利用A 的LU 分解求解下列方程组(a )Ax=b (b)A 2x=b (c) A 3x=b对方程组(c),若先求LU= A 3,再解(LU )x=b 有何缺点?(3)追赶法的优点实验题目:用追赶法分别对n=10,n=100,n=1000解方程组Ax=b ,其中------=4114114114A ,=5556 b 再用LU 分解法求解此方程组,并对二者进行比较。
三、实验原理顺序Gauss 消去法:步1 输入系数矩阵A ,右端项b ,置k:=1;步2 消元:对k=1,,1-n ,计算)()(/k kk k ik ik a a m =,o a k ik=+)1(,)()1()1()()()1(,k k ik k ik ik kj ik k ijk ijb m b b a m a a -=-=+++.n k j n k i ,,1;,,1 +=+= 步3 回代:,/)()(n nn n nn a b x =对k=,1,,1 -n 计算)(1)()(/)(k kk nk j j k kjk kk a x abx ∑+=-=列主元Gauss 消去法:步1输入系数矩阵A ,右端项b ,置k:=1;步2对k=1,,1-n ,进行如下操作(1)选列主元,确定k r ,使,m a x )()(k ik ni k k k r a a k≤≤=若)(k k r ka =0,则停止计算,否则,进行下一步;(2)若两行;的第交换k k r k r k,)b ,(A ,kk > (3)消元:对n k j n k i ,,1;,,1 +=+=,计算)()(/k kk k ik ik a a m =,o a k ik=+)1(, )()1()1()()()1(,k k ik k ik ik kj ik k ijk ijb m b b a m a a -=-=+++.步3 回代:,/)()(n nn n nn a b x =对k=,1,,1 -n 计算)(1)()(/)(k kk nk j j k kjk kk a x abx ∑+=-=LU 分解法:步1输入系数矩阵A ,右端项b ,置k:=1;步2 LU 分解.,,1,/)(;,,,u n,,2,k ;,,2,/;,,1,1111k j 111111n k i u u la l n k j u la n i u a l n j a u kk k r rk irik ik k r rj krkj i i j j +=-==-======∑∑-=-=计算对步3 用向前消去法解下三角方程组Ux=y: ∑-=-===11;,1,1,-n k ,/k j j kjk k nn n n y lb y u y x 计算对步4 用回代法解上三角方程组Ux=y :,/nn n n u y x = 对k=,1,,1 -n 计算kk nk j j kjk u x uy /)(x 1k ∑+=-=追赶法:先进行替换n k d b a d d c b a b b dd b b k k k kk k k k kk ,,2,,11111111=-=-===---- 1,2,,1,,1-=-==+n k b x c d x b d x kk k k k nn n四、实验内容(1)A=[0.3e-16, 59.14, 3, 1; 1, 2, 1, 1; 11.2, 9, 5, 2; 5.291, -6.13, -1, 2]; b=[51.97, 2, 1, 46.78]'; x1=magauss(A,b) x2=magauss2(A,b) (2)1)第一步function[A,b]=juzhen(n)A=zeros(n);b=zeros(n,1);b(n,1)=1;for i=1:nfor j=1:nA(i,i)=n;if i<j< p="">A(i,j)=n-j+i;elseA(i,j)==0;endendendAb第二步LU分解function [l,u]=lufj(A)[n,m]=size(A);n=length(b)u=zeros(n,n);l=eye(n,n);u(1,:)=A(1,:); l(2:n,1)=A(2:n,1)/u(1,1);for k=2:nu(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k); end 2)(a)[A,b]=juzhen(5);[x,l,u]=malu(A,b)(b)解下三角方程组Ly=bfunction y=lowt(l,b)n=length(b);y=zeros(n,1);y(1)=b(1);for k=2:ny(k)=b(k)-l(k,1:k-1)*y(1:k-1);end解上三角方程组Ux=yfunction x=upt(u,y)n=length(y);x=zeros(n,1);x(n)=y(n)/u(n,n);for k=n-1:-1:1x(k)=(y(k)-u(k,k+1:n)*x(k+1:n))/u(k,k); end 第三步y=lowt(l,b);z=upt(u,y);w=lowt(l,z);x2=upt(u,w)xb2=malu(A*A,b)(c)y=lowt(l,b);z=upt(u,y);w=lowt(l,z);s=upt(u,w);t=lowt(l,s);x3=upt(u,t)xb3=malu(A*A*A,b)(3)A=zeros(n,n);for(i=1:n)A(i,i)=4;endfor(i=1:n-1)A(i,i+1)=-1;endfor(i=2:n)A(i,i-1)=-1;追赶法:n=7;a=-ones(n,1); c=a;b=4*ones(n,1);d=5*ones(n,1); d(1)=6;xb1=machase(a,b,c,d);x1=xb1'LU分解求解x2=malu(A,d)五、实验结果实验1的结果:Warning: Divide by zero.> In magauss at 14In sinyan12 at 3x1 =NaNNaNNaNNaNx2 =3.78761.4680-15.062010.3384实验2 的结果:(1)l =1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 u = 5 4 32 1 0 5 43 2 0 0 54 3 0 0 05 4 0 0 0 0 5xb1 =0.01150.00960.0080-0.16000.2000xb2 =0.00160.00130.0288-0.06400.0400xb3 =0.0002-0.00450.0163-0.01920.0080实验3的结果:当n=10时 x1 =2.0981 2.3923 2.4711 2.4920 2.4970 2.4960 2.4870 2.4519 2.3205 1.8301 x2 =2.0981 2.3923 2.4711 2.4920 2.4970 2.4960 2.4870 2.45192.3205 1.8301六、实验结果分析1、实验一表明顺序高斯消元法有局限。
(完整word版)数值线性代数第二版徐树方高立张平文上机习题第四章实验报告
第四章上机习题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 ,否则,转step4 step4若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。
数值线性代数_MathCAD实验
算法
MathCAD代码
下 三 角 方
方程组
⎛ ⎜ ⎜⎜⎝
ll1211 ... ln1
l22 ... ln 2
... ...
lnn
⎞ ⎟ ⎟⎟⎠
⎛ ⎜ ⎜⎜⎝
yy12 ... yn
⎞ ⎟ ⎟⎟⎠
=
⎛ ⎜ ⎜⎜⎝
bb12 ... bn
⎞ ⎟ ⎟⎟⎠
解法:
程 (1) y1 = b1 / l11
∑ (前
代)
u x n
k=n+1 nk k
表示零。但MathCAD认为变量 k
=
n
+ 1, n
,由于U
是n
阶矩
阵, un,n+1 没有定义,计算出错。因此用MathCAD编写程序时,应特别注意循环变量的初值和终值。参
看以下代码:
Ub(U, b) := "求解上三角线性系统,下标1开始 n ← rows(U) ⊕ "正确程序"
n ← rows( A)
"先定义n-1个矩阵,存放Gauss变换 "
for i∈ 1 .. n − 1
Li ← identity(n) ⊕ "定义了n-1个单位矩阵"
for k ∈ 1 .. n − 1
for i∈ k + 1 .. n
( )Lk
i,k ←
−Ai , k Ak , k
⊕
"第k个Gauss变换阵"
⎛ a11 " a1,i-1 a1,i a1,i+1 " a1,n b1 ⎞
⎜
⎜ ⎜
0 0
%# " ai-1,i-1
数值分析实验报告-清华大学--线性代数方程组的数值解法
数值分析实验报告-清华大学--线性代数方程组的数值解法(总15页)--本页仅作为文档封面,使用时请直接删除即可----内页可以根据需求调整合适字体及大小--线性代数方程组的数值解法实验1. 主元的选取与算法的稳定性问题提出:Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss 消去法作为数值算法的稳定性呢?Gauss 消去法从理论算法到数值算法,其关键是主元的选择。
主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。
实验内容:考虑线性方程组 n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求:(1)取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ,则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元,结果如何?(2)现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。
(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。
(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。
重复上述实验,观察记录并分析实验结果。
程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A);Ab=[A b];r=input('选主元方式(0:自动;1:手动),r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元,请输入第',num2str(i),'列所选元素所处的行数:']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx运行结果(1)n=10,矩阵的条件数及自动选主元Cond(A,1) =×103Cond(A,2) = ×103Cond(A,inf) =×103程序自动选择主元(列主元)a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=0b.计算结果x=[1,1,1,1,1,1,1,1,1,1]T(2)n=10,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2…(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[, , , , , , , , , ]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k 步选择主元处于第k+1行) 最终计算得x=[1,1,1,1,1,1,1,1,1,1]T(3)n=20,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2…(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[,,,,,,,,,,,,,,,,,,,]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k步选择主元处于第k+1行)最终计算得x=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]T(4)A分别为幻方矩阵,Hilbert矩阵,pascal矩阵和随机矩阵简要分析计算(1)表明:对于同一矩阵,不同范数定义的条件数是不同的;Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。
数值代数实验报告(1)绝对经典
数值代数实验报告Numerical Linear Algebra And ItsApplications学生所在学院:理学院学生所在班级:计算数学10-1学生姓名:戈东潮指导教师:于春肖教务处2012年12月实验一实验名称: Poisson 方程边值问题的五点差分格式 实验时间: 2012年12月13日 星期四 实验成绩: 一、实验目的:通过上机利用Matlab 数学软件实现Poisson 方程的边值问题的五点差分格式的线性代数方程组来认真解读Poisson 方程边值问题的具体思想与方法,使我们掌握得更加深刻,并且做到学习与使用并存,增加学习的实际动手性,不再让学习局限于书本和纸上,而是利用计算机学习来增加我们的学习兴趣。
二、实验内容:利用Poisson 方程来解下列问题:⎪⎩⎪⎨⎧∂∈=∈=∂∂+∂∂-ΩΩ),(,),(),(,sin sin )(y x y x u y x y x y u x u 0222222πππ 其中}{的边界是ΩΩ∂<<∈Ω,1,0),(y x y x 。
边值问题的解释是y x y x u ππsin sin ),(=(1)用例题中模型问题的方法取N=5,10(也可以取N=20)列出五点差分格式的线性代数方程组三、实验过程:系数矩阵A 的Matlab 实现函数如下:%文件名:poisson.m function T=poisson(N) for i=1:N a(i)=4; end b=diag(a); for j=1:N-1b(j,j+1)=-1;endfor j=2:Nb(j,j-1)=-1;endfor i=1:NI(i)=-1;endI=diag(I);for i=1:N:N*Nfor j=1:NT(i+j-1,i:i+N-1)=b(j,:);endendfor i=1:N:N*N-Nfor j=1:NT(i+j-1,i+N:i+N-1+N)=I(j,:);endendfor i=1+N:N:N*Nfor j=1:NT(i+j-1,i-N:i-1)=I(j,:);endend求解常数项b的Matlab实现函数如下:%函数名:constant.mfunction b=constant(N)n=N+1;h=1/n;i=1;for y=1/n:1/n:N/nfor x=1/n:1/n:N/nf(i)=2*pi*pi*sin(pi*x)*sin(pi*y);i=i+1;endendb=h*h*f;b=b';四、实验结果(总结/方案)在Matlab运行窗口输入>>A=poisson(5) Enter输出结果A =Columns 1 through 114 -1 0 0 0 -1 0 0 0 0 0-1 4 -1 0 0 0 -1 0 0 0 00 -1 4 -1 0 0 0 -1 0 0 00 0 -1 4 -1 0 0 0 -1 0 00 0 0 -1 4 0 0 0 0 -1 0-1 0 0 0 0 4 -1 0 0 0 -10 -1 0 0 0 -1 4 -1 0 0 00 0 -1 0 0 0 -1 4 -1 0 00 0 0 -1 0 0 0 -1 4 -1 00 0 0 0 -1 0 0 0 -1 4 00 0 0 0 0 -1 0 0 0 0 40 0 0 0 0 0 -1 0 0 0 -10 0 0 0 0 0 0 -1 0 0 00 0 0 0 0 0 0 0 -1 0 00 0 0 0 0 0 0 0 0 -1 00 0 0 0 0 0 0 0 0 0 -10 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 Columns 12 through 220 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 00 -1 0 0 0 0 0 0 0 0 00 0 -1 0 0 0 0 0 0 0 00 0 0 -1 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 0 0 04 -1 0 0 0 -1 0 0 0 0 0 -1 4 -1 0 0 0 -1 0 0 0 00 -1 4 -1 0 0 0 -1 0 0 00 0 -1 4 0 0 0 0 -1 0 00 0 0 0 4 -1 0 0 0 -1 0 -1 0 0 0 -1 4 -1 0 0 0 -10 -1 0 0 0 -1 4 -1 0 0 00 0 -1 0 0 0 -1 4 -1 0 00 0 0 -1 0 0 0 -1 4 0 00 0 0 0 -1 0 0 0 0 4 -10 0 0 0 0 -1 0 0 0 -1 40 0 0 0 0 0 -1 0 0 0 -10 0 0 0 0 0 0 -1 0 0 00 0 0 0 0 0 0 0 -1 0 0 Columns 23 through 250 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 4 -1 0 -1 4 -10 -1 4 在运行窗口输入 >>b=constant(5) Enter 常数项的输出结果b =[ 0.1371 0.2374 0.2742 0.2374 0.1371 0.2374 0.4112 0.4749 0.4112 0.2374 0.2742 0.4749 0.5483 0.4749 0.2742 0.2374 0.4112 0.4749 0.4112 0.2374 0.1371 0.2374 0.2742 0.2374 0.1371]T所以A*u=b 的五点差分格式为j i j i j i j i j i j i f h u u u u u ,,,,,,211114=-----+-+ (i,j=1,2……)实验二实验名称: 用Jacobi 迭代法和SOR 迭代法求解方程组 实验时间: 2012年12月13日 星期四 实验成绩: 一、实验目的:通过上机利用Matlab 数学软件采用Jacobi 迭代法和SOR 迭代法来求解方程组,求解过程中了解两种迭代方法的基本思想与迭代过程,并分析两种迭代方法的收敛性和实用用以及他们的异同点。
数值分析线性方程组的实验报告包含代码解析
线性代数方程组的直接解法实验目的:线性方程组求解的直接法编程实现,实验内容:线性方程组求解的高斯消去法算法实现线性方程组求解的主元素消去法算法实现线性方程组求解的LU 分解得方法算法实现线性方程组求解追赶法算法实现实验比较:高斯消去、主元素消去、LU 分解都用实例 ⎪⎩⎪⎨⎧-=-+-=+-=++15.3181533126321321321x x x x x x x x x 这个进行比较。
知识理论解线性方程组的方法大致分为直接法和迭代法。
直接法是指假设计算过程中不产生舍入误差,经过有限次的运算可求的方程组精确解的方法。
方程(2-1)⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++nn nn n n n n n n b x a x a x a b x a x a x a b x a x a x a ...... (22112)222212111212111 记为:AX=b ;一、高斯(Gauss )消元法(1).Gauss 消元法是最基本的一种方法。
先逐次消去变量,将方程组化成同解的上三角方程组。
消元过程:先逐次消去变量,将方程组化成同解的上三角方程组 基本思想回代过程:按方程的相反顺序求解三角形方程组,得到原方程组的解。
(2) Gauss 消去法的求解思路为:若先让第一个方程组保持不变,利用它消去其余方程组中的,使之变成一个关于变元的n-1阶方程组。
按照(1)中的思路继续运算得到更为低阶的方程组。
经过n-1步的消元后,得到一个三角方程。
利用求解公式回代得到线性方程组的解。
①消元过程:第一次消元,,0)1(11≠a设 )1(11)1(11a a l i i =记(i=1,2,3……,n ).将(2-1)中第i 个方程减去第一个方程乘以1i l (i=1,2,3……,n ),完成第一次消元,⎪⎩⎪⎨⎧=++=++=+++)2()2(2)2(2)2(2)2(22)2(22)1(1)1(12)1(121)1(11...... ... n n nn n n n n n b x a x a b x a x a b x a x a x a 得等价方程组 (2-2))2()2(:b x A =简记为其中:)1(11)1()2(j i ij ij a l a a -=;)1(11)1()2(b l b b i i i -=次消元后经过1-k :⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧=+++=+++=+++=+++++=++++++++++++++++++++... ...... ......... ......)()(1)(1,)()(1)(,11)(1,1)(,1)()(1)(1,)()2(2)2(21)2(1,2)2(22)2(22)1(1)1(11)1(1,1)1(12)1(121)1(11k n n k nn k k k n k k nk k k n k n k k k k k k k k k k kn k kn k k k k k k kk n n k k k k n n k k k k b x a x a x a b x a x a x a b x a x a x a b x a x a x a x a b x a x a x a x a x a 简记为:)1()1(++=k k b x A其中:)()()1(k kjik k ij k ija l a a -=+ )()()1(k k ik k i k ib l b b -=+),,1,(n k j i +=按上述方法完成n-1次消元后得到。
数值线性代数第二版徐树方高立张平文上机习题第二章实验报告
(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 矩阵为病态的。
数值计算线性方程迭代法实验
实验报告一一、实验目的理解线性方程组直接法与迭代法思想,掌握常用算法的设计,掌握用MATLAB 实现的数值解法。
二、实验题目实验一 线性方程组迭代法实验 1、 迭代法的收敛速度用迭代法分别对n=20,n=200解方程组,b Ax =其中nn A ⨯⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛------------------=431513143151513143151513143151513143151314(1) 选取不同的初值0x 和不同的右端向量b,给定迭代误差,用两种迭代法计算,观测得到的迭代向量并分析计算结果给出结论;(2) 取定初值0x 和右端向量b,给定迭代误差,将A 的主对角元成倍放大,其余元素不变,用Jacobi 迭代法计算多次,比较收敛速度,分析计算结果并给出结论。
2、 SOR 迭代法松弛因子的选取用逐次超松弛(SOR )迭代法求解方程组,b Ax =其中 .5555551221212211212212121221121221212200199198321⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛----------=x x x x x x A (1) 给定迭代误差,选取不同的超松弛因子1>ω进行计算,观测得到的近似值向量并分析计算结果,给出你的结论;(2) 给定迭代误差,选取不同的低松弛因子1<ω进行计算,观测得到的近似值向量并分析计算结果,给出你的结论。
三、实验原理Jacobi 迭代法算法:步1 取初始点()0x ,精度要求ε,最大迭代次数N ,置0:=k ;步 2 由()n i x a b a x ni j j j ij i ii k i,,1,1,11 =⎪⎪⎭⎫⎝⎛-=∑≠=+ 或()b D x A D D x k 111)(--++-= 计算()1+k x ; 步3 若()()ε≤-∞+k k xx1,则停算,输出()1+k x作为方程组的近似解; 步4 若N k =,则停算,输出迭代失败信息;否则置1+=k k ,转布2。
数值线性代数大作业报告
数值线性代数实验大报告指导老师:赵国忠姓名:1108300001 刘帅1108300004 王敏1108300032 郭蒙一、实验名称:16题P75上机习题二、实验目的:编制通用的子程序,完成习题的计算任务三、实验内容与要求:P75上机习题先用熟悉的计算机语言将算法2.5.1编制成通用的子程序,然后再用所编制的子程序完成下面两个计算任务:(1) 估计5到20阶Hilbert 矩阵的无穷范数条件数。
(2) 设A n = 11...111................1-1 (01)-- 先随机地选取x ∈R n ,并计算出b=An x;然后再用列主元Gauss 消去法求解该方程组,假定计算解为∧x .试对n 从5到30估计计算解∧x 的精度,并且与真实的相对误差作比较。
四、 实验原理:(1)矩阵范数(martix norm )是数学上向量范数对矩阵的一个自然推广。
利用for循环和cond (a )Hilbert 求解Hilbert 矩阵的无穷范数,再利用norm(a,inf)求矩阵的无穷范数条件数。
(2)本题分为4步来求解。
先运用rand 随机选取x ∈R n,输入A n 矩阵,编制一个M 文件计算出b 。
第二步用列主元高斯消去法求解出方程的解X2。
第三步建立M 文件: soluerr.m 估计计算解∧x 的精度。
第四步, 建立M 文件: bijiao.m ,与真实相对误差作比较。
五、 实验过程:(1)程序:clearfor n=5:20for i=1:nfor j=1:na(i,j)=1/(i+j-1);endendc=cond(a);f=norm(c,inf);fprintf('n=%3.0f\nnorm(c,inf)%e\n',n,f) end运行结果:n= 5norm(c,inf)4.766073e+005n= 6norm(c,inf)1.495106e+007n= 7norm(c,inf)4.753674e+008n= 8norm(c,inf)1.525758e+010n= 9norm(c,inf)4.931542e+011n= 10norm(c,inf)1.602467e+013n= 11norm(c,inf)5.224376e+014n= 12norm(c,inf)1.698855e+016n= 13norm(c,inf)3.459404e+017n= 14norm(c,inf)4.696757e+017n= 15norm(c,inf)2.569881e+017n= 16norm(c,inf)7.356249e+017n= 17norm(c,inf)4.362844e+017n= 18norm(c,inf)1.229633e+018n= 19norm(c,inf)9.759023e+017n= 20norm(c,inf)1.644051e+018(2)程序:M文件:matrix1.mfunction [a,b,x1]=matrix1(n) format longA1=-1*ones(n,n)A2=tril(A1)for i=1:nA2(i,i)=1endA2(:,n)=1a=A2x1=rand(n,1)b=A2*x1end运行结果:>> A1 =-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1A2 =-1 0 0 0 0 -1 -1 0 0 0 -1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 -1 0 0 0 -1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0-1 -1 -1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 -1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 1 0 -1 -1 -1 -1 -1A2 =1 0 0 0 0 -1 1 0 0 0 -1 -1 1 0 0 -1 -1 -1 1 0 -1 -1 -1 -1 1A2 =1 0 0 0 1 -1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1a =-1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1x1 =0.8147236863931790.9057919370756190.1269868162935060.9133758561390190.632359246225410b =1.4470829326185890.723427496907850-0.961169560949882-0.301767337397875-2.128519049675914a =1 0 0 0 1 -1 1 0 0 1 -1 -1 1 0 1 -1 -1 -1 1 1 -1 -1 -1 -1 1b =1.4470829326185890.723427496907850-0.961169560949882-0.301767337397875-2.128519049675914x1 =0.8147236863931790.9057919370756190.1269868162935060.9133758561390190.632359246225410M文件:LZYgauss.mfunction[x2]=LZYgauss(a,b)format longn=length(a);x2=zeros(n,1);a=[a b];for k=1:n-1max=k;for i=k+1:nif a(i,k)>a(max,k)max=i;endendtemp=a(k,k:n+1);a(k,k:n+1)=a(max,k:n+1);a(max,k:n+1)=temp;for i=k+1:na(i,k)=-a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)+a(i,k)*a(k,k+1:n+1);endendx2(n,1)=a(n,n+1)/a(n,n);for i=n-1:-1:1sum=0;for j=i+1:nsum=sum+x2(j,1)*a(i,j);endx2(i,1)=(a(i,n+1)-sum)/a(i,i);end运行结果:>> LZYgauss(a,b)ans =0.8147236863931790.9057919370756190.1269868162935060.913375856139020 0.632359246225410估计计算解x的精度:M文件: soluerr.mfunction [x,error]=soluerr(a,b)format long%估计计算解的精度% 算法:列主元Gauss消去法,其中% A --- 系数矩阵% b-右端项% index --- index=0表示计算成败;index=1表示计算成功%输出结果:error--本算法给出的计算解的估计% normA--逆矩阵无穷范数估计% rnorm--计算解的残量[n,m]=size(a); nb=length(b);if n~=merror('The rows and columns of matrix a must be equal!');return;endif m~=nberror('The columns of a must be equal the dimension of b!'); return;endindex=1;%列主元矩阵三角分解[L,U,u,index_col]=Gauss_col(a);%解下三角方程组Ly=Pb[y,index_low]=Gauss_low(L,b(u));%解上三角方程组Ux=y[x,index_upp]=Gauss_upp(U,y);%输出数值解xpause(0.3)%估计矩阵逆的无穷大范数normA=normAinv(L,U,u);%估计计算解的残量rnorm=norm(b-a*x,inf);%计算右端项bnorm=norm(b,inf);%计算矩阵A的范数Anorm=norm(a,inf);%计算解的精度error=normA*Anorm*rnorm/bnorm;运行结果:>> [x,error]=soluerr(a,b)x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410error =5.215941218821592e-016X1与真实相对误差做比较M文件:bijiao.mfunction [x,x1,error,error1]=bijiao(a,b,n) [a,x1]=matrix1(n)[x,error]=soluerr(a,b)error1=abs((x-x1)/x1)运行结果:>>x =0.8147236863931790.9057919370756190.1269868162935060.9133758561390200.632359246225410x1 =0.9578935505663481.6132601689718680.629241553693582-0.799716694069468-1.770467991515150error =5.215941218821592e-016error1 =Columns 1 through 40 0 0 00 0 0 00 0 0 00 0 0 00 0 0 0Column 50.0808655478999350.3995939126190040.2836847318376260.9675930648949141.357170674226221六、结果分析:1、矩阵范数(martix norm)是数学上向量范数对矩阵的一个自然推广。
实验5_线性代数方程组的数值解法
实验5 线性代数方程组的数值解法化工系 分0班 2010011805 张亚清【实验目的】1、 学会用MATLAB 软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2、 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】 1、 题目3已知方程组Ax=b ,其中2020⨯ℜ∈A ,定义为试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求:(1)选取不同的初始向量)0(x和不同的方程组右端向量b ,给出迭代误差要求,用雅克比迭代法和高斯-赛德尔迭代法计算,观测得到的迭代向量序列是否均收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论;(2)取定右端向量b 和初始向量)0(x,将A 的主对角线元素成倍增长若干次,非主对角线元素不变,每次用雅克比迭代法计算,要求迭代误差满足5)()1(10-∞+<-k k x x ,比较收敛速度,分析现象并得出你的结论。
【问题分析】对于线性方程组Ax=b ,满足0≠ii a (i=1,2,...,n ),将A 分解为A=D-L-U ,则雅克比迭代公式等价于如下的矩阵形式,...2,1,0,)(1)(1)1(=++=--+k b D x D L D x k k或b D f U L D B f x B xJ J J k J k 11)()1(),(,--+=+=+=。
类似的,Ax=b 的高斯-赛德尔迭代公式等价于如下矩阵形式b L D f U L D B f x B x J S G S G k S G k 11)()1()(,)(,-----+-=-=+=。
【问题解答】(1)选取初始向量T x )1,...,1,1()0(=,T b )1,...,1,1(=,迭代要求为4)()1(10-∞+<-k k x x 。
将A 按A=D-L-U 分解为如下三个矩阵:①对方程组进行雅克比迭代,利用MATLAB 编程计算。
数值线性代数实验
数值线性代数实验题目:数值线性代数专业:信息与计算科学班级:班姓名:山东科技大学2013年 1 月16日实验报告说明学院:信息学院专业:信息班级10-2 姓名:一、主要参考资料:(1)《Matlab数值计算-案例分析》北京航空出版(2)《Matlab数值分析》机械工业出版二、课程设计应解决的主要问题:(1)平方根(2)QR方法(3)最小二乘法三、应用软件:(1)Matlab7.0(2)数学公式编辑器四、发出日期:课程设计完成日期:指导教师签字:系主任签字:指导教师对课程设计的评语指导教师签字:年月日一、问题描述先用你所熟悉的计算机语言将平方根和改进的平方根法编成写通用的子程序,然后用你编写的程序求解对称正定方程组b x =A ,其中 (1)b 随机的选取,系数矩阵位100阶矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1011101110111011101110(2)系数矩阵为40阶Hilbert 矩阵,即系数矩阵A 的第i 行第j 列元素为11-+=j i a ij ,向量b 的第i 个分量为∑=-+=nj i j i b 111。
二、分析与程序1. 平方根法函数程序如下: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 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=bL=zeros(n,n);D=diag(n,0);S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;endfor i=1:nfor j=1:nif (eig(A)<=0)|(A(i,j)~=A(j,i))disp('wrong');break;endendendD(1,1)=A(1,1);for 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=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); endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));end2.改进平方根法函数程序如下: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);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=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'); end3.调用函数解题: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[x,b]=pingfanggenfa(A,b)b=gaijinpinfanggenfa(A,b)4.运行结果:请输入矩阵维数: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.000012.0000 -69.000068.000021.000017.0000 -50.7188-8.7500-8.0000 112.00006.0000 -68.750022.000044.0000 -28.00008.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.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.00310.0036-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.00120.0139-0.09540.4208-1.21012.0624-1.0394-3.33436.2567-0.2463-7.45942.80303.69900.7277-1.7484-0.4854-3.60100.25325.1862-2.12991.44100.8738-4.56541.04224.0920-2.7764-2.2148-0.89530.36654.89671.04160.1281 -4.3387 -1.1902 -2.8334 8.4610 -3.6008一、问题描述先用你所熟悉的计算机语言将算法2.5.1编成写通用的子程序,然后用你编写的程序完成下面两个计算任务:(1) 估计5到20阶Hilbert 矩阵的∞范数条件数;(2)设n *n n R 11-1-1-111-1-101-1001A ∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=先随机的选取n R x ∈,并计算出x A b n =;然后再用列主元Gauss 消去法求解该方程组,假定计算解为xˆ。
大学数学实验五_线性代数方程组的数值解法
【实验目的】 1、学会用 MATLAB 软件数值求解线性代数方程组,对迭代法的收敛性和解
的稳定性作初步分析。 2、通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】
3 已知方程组 Ax=b,其中
,定义为
试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对 收敛速度的影响。实验要求: (1) 选取不同的初始向量 x(0)和不同的方程组的右端项向量 b,给定迭代误差要求,用雅
k=k+1; xj=Bj*xj+fj; 多输出了矩阵 P,矩阵 P 可视为一个行向量,其每个元素均为迭代 k 次后得到的 xk。这样以 k 为横轴,解向量为纵轴,可输出图形观察 xk 是否收敛。函数 GaussSeidel 也需作同样修改,修改后的函数在此不再赘述。
模型: 已知某年该植物的数量为 x0,记第 k 年的植物数量为 xk,那么有 xk + pxk-1 + qxk-2 = 0 (k = 2, 3, …… , n)
其中 p = -a1bc,q = -a2b(1-a1)bc。若要求 n 年后数量达到 xn,则 Ax = b
其中
,
,
7
① 用稀疏系数矩阵求解。
这个函数中,n 表示矩阵 A 的阶数,在本题中恒取 20,a 表示主对角线元素的值,b 在 本题中恒取-1/4,c 在本题中恒取-1/2。
编写用雅可比迭代法求方程解的函数 Jacobi。
function [xj,k]=Jacobi(A,X0,b,e) D=diag(diag(A)); n=length(A); L=-(tril(A)-D); U=-(triu(A)-D); fj=D\b; Bj=D\(L+U); xj=X0; k=0; while norm(A*xj-b)/norm(b)>e
数学实验报告——利用MALTAB计算线性代数方程组的数值解法
实验三线性代数方程组的数值解法一、迭代法求解方程组㈠问题描述给定方程组的矩阵A,通过迭代法求解方程组。
1、选取不同的初始向量和不同的右端项向量,给定误差要求,用两种迭代法计算;2、去顶右端项向量和初始向量,将A的主对角线元素成倍增长若干次,非主对角线元素不变,用雅克比迭代法计算。
㈡方法与公式1、雅克比迭代法2、高斯-赛德尔迭代法㈢结果与分析1、不同初始向量、不同右端项向量、不同精度要求(1)初始向量定为zeros(n,1);①b=zeros(n,1)迭代次数为0,直接得到结果。
③b = 1:n(2)初始向量定为one s(n,1)①b=zeros(n,1)事实上,迭代100次时,所得结果约为10^-32,已经可以认为是0,但是由于没有达到精度要求,故不算收敛。
②b=ones(n,1)④b = n:1(3)初始向量定为1:n①b=zeros(n,1)②b=ones(n,1)(4)初始向量定为n:1①b=zeros(n,1)②b=ones(n,1)④b = n:1(5)简要小结a.在个别情况下雅可比迭代法收敛速度极慢,但事实上没有达到收敛时其计算结果已经可以接受;b.要求的精度越高,迭代的次数越多,迭代的次数与所要求的精度的对数值近似呈线性,也就是说两者近似呈指数关系;b.高斯-赛德尔迭代法有着比雅可比更好的迭代特性;2、更改A的主对角线元素(1) b =20:1;初值 1:20(2) b =20:1;初值 20:1(3) b =1:20;初值 20:1(4) b =[3;5;2;6;8;23;5;8;32;63;23;5;2;12;0;23;1;564;2;65]; 初值 ones(20,1)(5)简要小结a.迭代的次数随着对角线元素的成倍的增长而降低,趋于一稳定值;b.右端项以及迭代初值仅当对角线元素较小时对迭代次数起有作用,对角线元素成倍数增加后,迭代次数不变。
3、总结由以上各个对比可以得出以下结论:a.使用迭代法求解方程组时时,要求的精度越高,迭代次数越大;b.高斯-赛德尔迭代法的迭代次数要比雅可比迭代法迭代次数低;c.雅可比迭代的次数随着矩阵A对角线元素的成倍的增长而降低,d.当矩阵A的对角线元素足够大时,雅可比迭代法的迭代次数趋于稳定值;㈣程序清单1、第一问中的雅可比迭代function [y,k] = jacobi(A,b,m,tol)D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);y = ones(n,1);BJ=D\(L+U);fJ=D\b;k=0;while norm(A*y-b)/norm(b)>tol && k<mk = k+1;y = BJ*y+fJ;end2、高斯-赛德尔迭代法function [y,k] = GuassSeidel(A,b,m,tol)D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);BG = (D-L)\U;FG = (D-L)\b;y = zeros(n,1);k=0;while norm(A*y-b)/norm(b)>tol && k<mk = k+1;y = BG*y+FG;end3、第二问中的雅可比迭代function [y,k] = jacobi1(A,b,m,tol) D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);yk = ones(20,1);BJ=D\(L+U);fJ=D\b;k=0;yk1 = BJ*yk+fJ;k = k+1;ttol = norm(yk1-yk,inf);while ttol>tol && k<mk = k+1;yk = yk1;yk1 = BJ*yk+fJ;ttol = norm(yk1-yk,inf);endy=yk1;4、第一问脚本n=20;A1 = sparse(1:n,1:n,3,n,n);A2 = sparse(1:n-1,2:n,-1/2,n,n);A3 = sparse(2:n,1:n-1,-1/2,n,n);A4 = sparse(1:n-2,3:n,-1/4,n,n);A5 = sparse(3:n,1:n-2,-1/4,n,n);A = A1+A2+A3+A4+A5;b =zeros(20,1);m=10000;for i = 3:10tol = 10^(-i);[xJ,k1(i)] = jacobi(A,b,m,tol); [xG,k2(i)] = GuassSeidel(A,b,m,tol); endk1k25、第二问脚本n=20;k = 1;A1 = sparse(1:n,1:n,3,n,n);A2 = sparse(1:n-1,2:n,-1/2,n,n);A3 = sparse(2:n,1:n-1,-1/2,n,n);A4 = sparse(1:n-2,3:n,-1/4,n,n);A5 = sparse(3:n,1:n-2,-1/4,n,n);b =[3;5;2;6;8;23;5;8;32;63;23;5;2;12;0;23;1;564;2;65]; m=10000;tol = 10^(-5);for k = 1:10A = k*A1+A2+A3+A4+A5;[xJ,k1(k)] = jacobi1(A,b,m,tol);endk1二、一年生植物的繁殖㈠问题描述对于5.1.2的假设,给定参数,试分别用追赶法、稀疏系数矩阵和满矩阵求解;若b 有10%误差,估计对结果的影响。
数学实验——线性代数方程组的数值解
实验5 线性代数方程组的数值解法分1 黄浩 2011011743一、实验目的1.学会用MATLAB软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2.通过实例学习用线性代数方程组解决简化的实际问题。
二、实验内容1.《数学实验》第二版(问题1)问题叙述:通过求解线性方程组A1x=b1,A2x=b2,理解条件数的意义和方程组性态对解的影响,其中A1是n阶范德蒙矩阵,即A1=[1x0x02…x0n−11x1x12…x1n−1……………1x n−1x n−12 (x)n−1n−1], x k=1+0.1k , k=0,1,…,n−1A2是n阶希尔伯特矩阵,b1,b2分别是A1,A2的行和。
(1)编程构造A1(A2可直接用命令产生)和b1,b2;你能预先知道方程组A1x=b1和A2x=b2的解吗?令n=5,用左除命令求解(用预先知道的解可验证程序)。
(2)令n=5,7,9,…,计算A1和A2的条件数。
为观察他们是否病态,做以下试验:b1,b2不变,A1和A2的元素A1(n,n),A2(n,n)分别加扰动ε后求解;A1和A2不变,b1,b2的分量b1(n),b2(n)分别加扰动ε后求解。
分析A与b的微小扰动对解的影响。
ε取10^-10,10^-8,10^-6。
(3)经扰动得到的解记做x̃,计算误差‖x−x̃‖‖x‖,与用条件数估计的误差相比较。
模型转换及实验过程:(1)小题.由b1,b2为A1,A2的行和,可知方程组A1x=b1和A2x=b2的精确解均为n 行全1的列向量。
在n=5的情况下,用matlab编程(程序见四.1),构造A1,A2和b1,b2,使用高斯消去法得到的解x1,x2及其相对误差e1,e2(使用excel计算而得)为:由上表可见,当n=5时,所得的解都接近真值,误差在10^-12的量级左右。
(2)小题分别取n=5,7,9,11,13,15,计算A1和A2的条件数c1和c2,(程序见四.2),结果如下:由上表可见,二者的条件数都比较大,可能是病态的。
数学实验 5:线性代数方程组的数值解法
大倍数 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000
小值q 0.4893 0.2447 0.1631 0.1223 0.0979 0.0816 0.0699 0.0612 0.0544 0.0489 0.0445 0.0408 0.0376 0.0350 0.0326 0.0306 0.0288 0.0272 0.0258 0.0245 21.0000 12.0000 9.0000 8.0000 8.0000 7.0000 7.0000 7.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 5.0000 5.0000
实验 5:线性代数方程组的数值解法
习题3:
已知方程组,其中,定义为: 试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方 程组系数矩阵性质对收敛速度的影响。实验要求: (1) 选取不同的初始向量x0和不同的方程组右端向量b,给定迭 代误差要求,用雅可比迭代法和高斯-赛德尔迭代法计算, 观测得到的迭代向量序列是否均收敛?若收敛,记录迭代 次数,分析计算结果并得出结论; (2) 取定右端向量b和初始向量x0,将A的主对角线元素成倍的 增长若干次,非主对角元素不变,每次用雅可比迭代法计 算,要求迭代误差满足,比较收敛速度,分析现象并得出结 论。 1、 程序设计(可直接粘贴运行) 1) Jacobi迭代法 function y=jacobi(a,b,x0,e,m) %定义jacobi函数,其中:a,b为线性方程组中的矩阵和右端向量;x0 为初始值; %e和m分别为人为设定的精度和预计迭代次数;运行结果y为迭代的结 果和所有中间值组成的 %矩阵 y=0; %对y初始化 d=diag(diag(a)); %按雅可比迭代标准形 形式取主对角元素作为矩阵D u=-triu(a,1); %取上三角矩阵u l=-tril(a,-1); %取下三角矩阵l bj=d^-1*(l+u); fj=d^-1*b; x=[x0,zeros(20,m-1)]; %初始化x,其中x1=x0,即 初始值 for k=1:m %人为规定迭代次 数,防止不收敛迭代导致死循环 x(:,k+1)=bj*x(:,k)+fj; %jacobi迭代 if norm(x(:,k+1)-x(:,k),inf)<e
数值计算方法实验报告
数值计算方法实验报告一、实验目的本实验旨在通过数值计算方法的实验操作,深入理解数值计算方法的原理与应用,掌握数值计算方法的相关技能,提高数值计算方法的实际应用能力。
二、实验内容1.数值微积分2.数值代数3.数值微分方程4.数值线性代数5.数值优化6.数值统计分析7.数值随机模拟8.数值傅立叶分析9.数值偏微分方程三、实验步骤1.数值微积分:通过不同的数值积分方法,计算给定函数的定积分值,并对不同数值积分方法的误差进行分析。
2.数值代数:通过使用线性代数方法,求解给定的线性方程组,并分析不同线性方程组求解方法的优劣。
3.数值微分方程:通过使用常微分方程数值解法,求解给定的微分方程,并比较不同求解方法的精度和稳定性。
4.数值线性代数:通过使用特征值分解方法,对给定的矩阵进行特征值分解,并分析不同特征值分解方法的优缺点。
5.数值优化:通过使用不同的优化方法,求解给定的优化问题,并比较不同的优化方法的效率和精度。
6.数值统计分析:通过使用不同的统计分析方法,对给定的数据进行统计分析,并分析不同的统计方法的优缺点。
7.数值随机模拟:通过使用随机模拟方法,模拟给定的概率分布,并分析不同随机模拟方法的效率和精度。
8.数值傅立叶分析:通过使用傅立叶分析方法,对给定的信号进行频谱分析,并分析不同的傅立叶分析方法的优缺点。
9.数值偏微分方程:通过使用偏微分方程数值解法,求解给定的偏微分方程,并比较不同求解方法的精度和稳定性。
四、实验结果与分析本实验中,通过对不同的数值计算方法的实验操作,我们可以更深入地理解数值计算方法的原理与应用,并掌握数值计算方法的相关技能,提高数值计算方法的实际应用能力。
同时,通过实验结果的分析,我们可以更好地比较不同数值计算方法的优缺点,为实际应用提供参考依据。
五、实验总结本实验旨在通过数值计算方法的实验操作,深入理解数值计算方法的原理与应用,掌握数值计算方法的相关技能,提高数值计算方法的实际应用能力。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
数值线性代数实验题目:数值线性代数专业:信息与计算科学班级:班姓名:山东科技大学2013年 1 月16日实验报告说明学院:信息学院专业:信息班级10-2 姓名:一、主要参考资料:(1)《Matlab数值计算-案例分析》北京航空出版(2)《Matlab数值分析》机械工业出版二、课程设计应解决的主要问题:(1)平方根(2)QR方法(3)最小二乘法三、应用软件:(1)Matlab7.0(2)数学公式编辑器四、发出日期:课程设计完成日期:指导教师签字:系主任签字:指导教师对课程设计的评语指导教师签字:年月日一、问题描述先用你所熟悉的计算机语言将平方根和改进的平方根法编成写通用的子程序,然后用你编写的程序求解对称正定方程组b x =A ,其中 (1)b 随机的选取,系数矩阵位100阶矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1011101110111011101110(2)系数矩阵为40阶Hilbert 矩阵,即系数矩阵A 的第i 行第j 列元素为11-+=j i a ij ,向量b 的第i 个分量为∑=-+=nj i j i b 111。
二、分析与程序1. 平方根法函数程序如下: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 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=bL=zeros(n,n);D=diag(n,0);S=L*D;for i=1:n %L的主对角元素均为1L(i,i)=1;endfor i=1:nfor j=1:nif (eig(A)<=0)|(A(i,j)~=A(j,i))disp('wrong');break;endendendD(1,1)=A(1,1);for 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=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); endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));end2.改进平方根法函数程序如下: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);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=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'); end3.调用函数解题: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[x,b]=pingfanggenfa(A,b)b=gaijinpinfanggenfa(A,b)4.运行结果:请输入矩阵维数: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.000012.0000 -69.000068.000021.000017.0000 -50.7188-8.7500-8.0000 112.00006.0000 -68.750022.000044.0000 -28.00008.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.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.00310.0036-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.00120.0139-0.09540.4208-1.21012.0624-1.0394-3.33436.2567-0.2463-7.45942.80303.69900.7277-1.7484-0.4854-3.60100.25325.1862-2.12991.44100.8738-4.56541.04224.0920-2.7764-2.2148-0.89530.36654.89671.04160.1281 -4.3387 -1.1902 -2.8334 8.4610 -3.6008一、问题描述先用你所熟悉的计算机语言将算法2.5.1编成写通用的子程序,然后用你编写的程序完成下面两个计算任务:(1) 估计5到20阶Hilbert 矩阵的∞范数条件数;(2)设n *n n R 11-1-1-111-1-101-1001A ∈⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=先随机的选取n R x ∈,并计算出x A b n =;然后再用列主元Gauss 消去法求解该方程组,假定计算解为xˆ。
试对n 从5到30估计计算解x ˆ的精度,并且与真实相对误差做比较。
二、分析与程序Hilbert 矩阵:Hilbert 矩阵的分量满足H(i,j)=1/(i+j-1) 比如3阶Hilbert 矩阵是 1/1 1/2 1/3 1/2 1/3 1/4 1/3 1/4 1/5程序1. 用MathCAD 计算方法如下:H n ()A i j,1i j +1-←j 1n ..∈for i 1n ..∈for A:=κn ()normi H n ()()normi H n ()1-()⋅:=其中()H n 表示n 阶Hilbert 矩阵,()n κ即是n 阶Hilbert 矩阵的∞范数条件数2. 比较估计精度与真实相对误差用MathCAD 计算方法如下:F n ()A i i ,1←Ai n,1←i 1n ..∈for A i j,1-←j 1i 1-..∈for i 2n ..∈for A:=X n ()x i1.19998←i 1n ..∈for x:=L n ()Li i,1←i 1n ..∈for L i j,1-←j 1i 1-..∈for i 2n ..∈for L:=b n ()F n ()X n ()⋅:= U n ()L n ()1-F n ()⋅:=y n ()L n ()1-b n ()⋅:= x n ()U n ()1-y n ()⋅:= r n ()b n ()F n ()x n ()⋅-:=P n ()p i 1,normiF i ()()maxr i ()()⋅normi F i ()()T ()1-⎡⎣⎤⎦⋅maxb i ()()←i 5n ..∈for p:=T n ()ti 1,max X i ()x i ()-()max X i ()()←i 5n..∈for t:=其中()F n 表示n 阶矩阵n A ,()X n 是任选的n x R ∈,()b n 即是n b A x =,()L n 和()U n 分别是对n A 进行LU 分解得到的下三角矩阵和上三角矩阵,()x n 是用列主元Gauss 消去法求得的解,()P n 是n 阶矩阵n A 计算解的精度,()T n 是n 阶矩阵n A 计算解的真实相对误差。
3. 估计∞范数条件数用MathCAD 处理计算的结果如下:当矩阵阶数 m 520..:= 时,()m κ即是m 阶Hilbert 矩阵的∞范数条件数。
以矩阵阶数m 为横坐标,条件范数()m κ为纵坐标,其图像如下:51015205.10171.1018κm ()m条件数在一定程度上刻画了扰动对方程组解的影响程度。