数值线性代数实验

合集下载

线性代数实验报告

线性代数实验报告

2.输入:for
n=20:80 p1(n)=prod(365-n+1:365)/365^n; p(n)=1-p1(n); end plot(p)
输出:

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)” ,

输出:
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版)数值线性代数第二版徐树方高立张平文上机习题第四章实验报告

(完整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实验

算法
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)绝对经典

数值代数实验报告(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. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 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条件数在一定程度上刻画了扰动对方程组解的影响程度。

相关文档
最新文档