数值线性代数实验

合集下载

数值分析课程实验设计——数值线性代数实习题

数值分析课程实验设计——数值线性代数实习题

数值分析课程实验设计——数值线性代数
实习题
1. 实验目的
本实验的主要目的是进一步加深对数值线性代数的理解,熟悉
常见矩阵分解方法,并在此基础上解决实际问题。

2. 实验内容
本次实验将任务分为两个部分,分别是矩阵分解与求解线性方
程组。

2.1 矩阵分解
首先,我们需要熟悉三种常见的矩阵分解:QR分解、LU分解
和奇异值分解。

我们需要通过Python语言实现这三种分解方法,
并利用这些方法解决实际问题。

2.2 求解线性方程组
其次,我们需要学会用矩阵分解的方法来求解线性方程组。


们将通过两个例子来进行说明,并利用Python语言实现这些方法。

3. 实验要求
本次实验要求熟悉矩阵分解的基本方法,在此基础上解决实际问题;能够运用多种方法来求解线性方程组,并分析比较它们的优缺点。

4. 实验总结
本次实验通过矩阵分解和求解线性方程组两个部分的学习,巩固了我们对于数值线性代数的知识,并在实际问题的解决中得到了应用。

感谢老师的指导,我们会在今后的学习中持续探索数值分析方面的知识。

数值代数实验报告

数值代数实验报告

数值代数实验报告数值代数实验报告引言:数值代数是一门研究数值计算方法和算法的学科,它在科学计算和工程应用中起着重要的作用。

本实验报告旨在通过实际的数值计算问题,探讨数值代数的应用和效果。

实验一:线性方程组求解线性方程组求解是数值代数中的一个重要问题。

在实验中,我们使用了高斯消元法和LU分解法两种求解线性方程组的方法,并对比了它们的效果。

首先,我们考虑一个3×3的线性方程组:2x + 3y - z = 54x - 2y + 2z = 1x + y + z = 3通过高斯消元法,我们将该方程组转化为上三角形式,并得到解x=1, y=2, z=0。

而通过LU分解法,我们将该方程组分解为LU两个矩阵的乘积,并得到相同的解。

接下来,我们考虑一个更大的线性方程组,例如10×10的方程组。

通过比较高斯消元法和LU分解法的运行时间,我们可以发现LU分解法在处理大规模方程组时更加高效。

实验二:特征值与特征向量计算特征值与特征向量计算是数值代数中的另一个重要问题。

在实验中,我们使用了幂法和QR方法两种求解特征值与特征向量的方法,并对比了它们的效果。

首先,我们考虑一个3×3的矩阵:1 2 34 5 67 8 9通过幂法,我们可以得到该矩阵的最大特征值为15.372,对应的特征向量为[0.384, 0.707, 0.577]。

而通过QR方法,我们也可以得到相同的结果。

接下来,我们考虑一个更大的矩阵,例如10×10的矩阵。

通过比较幂法和QR 方法的运行时间,我们可以发现QR方法在处理大规模矩阵时更加高效。

实验三:奇异值分解奇异值分解是数值代数中的一种重要技术,它可以将一个矩阵分解为三个矩阵的乘积,从而实现数据降维和信息提取的目的。

在实验中,我们使用了奇异值分解方法,并通过实际的数据集进行了验证。

我们选取了一个包含1000个样本和20个特征的数据集,通过奇异值分解,我们将该数据集分解为三个矩阵U、S和V的乘积。

数值代数(第13周)实验报告()

数值代数(第13周)实验报告()
[x,k]=fastest(A,b,eps);
运行结果:
x =
-4.0000
3.0000
2.0000
k =
17
(2)第二题的第一个方程组
A=[1 0.4 0.4;0.4 1 0.8;0.4 0.8 1];
b=[1 2 3]';
eps=0.00001;
[x,k]=fastest(A,b,eps);
运行结果:
(本次实验进行总结,包括但不限于对算法的理解、编程实验中遇到的问题、实验中所取得的经验等;如程序未能通过,应分析错误原因.)
通过不断试验最速下降法和范数norm(x-x0)取值范围有关.前两个方程组当eps取得越小误差越小。上面的运算结果则是最接近值的迭代次数
x =
-0.1351
-1.0811
3.9189
k =
51
(3)第二题的第二个方程组
A=[1 2 -2;1 1 1;2 2 1];
b=[1 1 1]';
eps=0.1;
[x,k]=fastest(A,b,eps);
运行结果:x =
-3.0000
3.0000
1.0000
k =
5
4实验总结【必填栏目】
得分:
x0=x;
ifk>=m
disp('迭代次数太多,可能 Nhomakorabea收敛!');
return;
end
end
x
k
3.2实验(或测试)代码和结果
得分:
用最速下降法求解课本209页习题1、2中的三个方程组.
(1)第一题的方程组
A=[5 2 1;-1 4 2;2 -3 10];

【清华】7.0_实验5-线性代数方程组的数值解法

【清华】7.0_实验5-线性代数方程组的数值解法

0
0.5000
0.7500
0.9926
0.99979 0.99999 0.99999
0
0.5000
0.7500
0.9930
0.99981 0.99999 0.99999
0
0.5000
0.7500
0.9935
0.99983 0.99999
1
0
0.5000
0.7569
0.9943
0.99985 0.99999
B1=D\(L+U);
f1=D\b; x(:,1)=x0; x(:,2)=B1*x(:,1)+f1; k=1; while norm((x(:,k+1)-x(:,k)),inf)>m
x(:,k+2)=B1*x(:,k+1)+f1; k=k+1; end
%赋初值 %以差值向量的行范数作为误差判断标准
0.22222
0 0.083333 0.095238
0.21996
0 0.083333 0.095016
0 0.33333
0.6521
0 0.16667
0.22152
0 0.083333 0.095181
0 0.33333
0.66094
0 0.16667
0.22207
0 0.083333 0.095231
0 0.33333
0.6643
0 0.16667
0
0.7083
0.9069
0.9967
0.99999
1
1
0
0.6806
0.8960
0.9962
0.99998 0.99999

数值线性代数_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 迭代法来求解方程组,求解过程中了解两种迭代方法的基本思想与迭代过程,并分析两种迭代方法的收敛性和实用用以及他们的异同点。

数值实验 线性代数方程组的性态与条件数的估计

数值实验 线性代数方程组的性态与条件数的估计
#34; " ( )中的结果给出 1 x x 的理论估计,并将它与1) (
给出的计算结果进行比 较,分析所得结果。注 意, 如果给出了cond( A)和 A 的估计,马上就可以给 出 A1 的估计
数值实验二
邹昌文
线性代数方程组的性态与条件 数的估计
• 问题的提出:
理论上,线性代数方程 Ax b的摄动满足 组 x A b cond( A) ( ) 1 x b 1 A A A
实验内容
Matlab中提供有函数 condest 可以用来估计 " " 矩阵的条件数,它给出 的是按1 范数的条件 数。首先构造非奇异矩 A和右端,使得方 阵 程是可以精确求解的。 再人为地引进系数矩 阵和右端的摄动 A和b,使得 A 和 b 充分 小
实验要求
ˆ 1 )假设方程Ax b的解为x,求解方程( A A) x b b 以1 范数,给出 x x ˆ xx x 的计算结果
2)选择一系列维数递增 的矩阵(可以是随机生 成的), 比较函数" condest 所需机器时间的差别。 " 考虑若干逆 是已知矩阵,借助函数 eig" 很容易给出cond2 ( A)的数值, " 将它与函数 cond( A,2)"所得到的结果比较 "

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

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

(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 矩阵为病态的。

实验6线性代数方程组的数值解法

实验6线性代数方程组的数值解法

实验6线性代数⽅程组的数值解法实验6 线性代数⽅程组的数值解法[实验⽬的]1. 1. 学会⽤MATLAB 软件数值求解线性代数⽅程组,对迭代法的收敛性和解的稳定性作初步分析;2. 2. 通过实例学习⽤线性代数⽅程组解决简化的实际问题。

[实验内容]5-5 输电⽹络:⼀种⼤型输电⽹络可简化为图5.5(见书)所⽰电路,其中R 1,R 2,…,R n 表⽰负载电阻,r 1,r 2,…,r n 表⽰线路内阻,I 1,I 2,…,I n 表⽰负载上的电流。

设电源电压为V 。

(1)列出求各负载电阻R 1,R 2,…,R n 的⽅程;(2)设I 1=I 2=…=I n =I ,r 1=r 2=…=r n =r ,在r=1,I=0.5,V=18,n=10的情况下求R 1,R 2,…,R n 及总电阻R 0。

[问题分析、模型建⽴及求解](1) 设电源负极为电势为0,电阻R 1上对应节点电压为V 1,对于任意节点,根据KCL 定律列出⽅程:111++----=k k k k k k k k r V V r V V R V⽽k kk R V I =,可得:111111)(++++--++-=k k k k k k k k k k k k R r IR r I r I R r I Ik=2,3,……,n-1;k=1时2221211R r I R r I I +-=,为与上式形式⼀致,化为22212111111)(R r IR r I r I r V I +--=-k=m (12-≤≤n m )时 111111)(++++--+--+=m m m n m m m m m m m m R r IR r I r I R r I Ik=n 时R r I I -=--11设以上⽅程组的矩阵形式为:b AR = 则 []Tn R R R R ΛΛΛΛ21=Tn I I I r V I b ?-=ΛΛΛ3211---------=------n n nn n n nn n n n n r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I A 11 11114443333233322221222111000000ΛΛΛO O M O O O O M M O O O O O M M O M O ΛΛ5.021=====I I I I n Λ,121=====r r r r n Λ,V=18,n=10,-----=5.05.0005.015.005.015.005.015.0005.01ΛΛΛO M O O O O M M O O O O O MM O M O ΛΛΛA[]Tb 5.05.05.05.17ΛΛΛ-=在命令窗⼝输⼊MA TLAB 程序如下:clear all ;n=10; %由题⽬要求设定A11=sparse(1:n-1,1:n-1,-1,n,n); %定义A 的对⾓元素,除(n,n) A12=sparse(n,n,-0.5,n,n); %定义(n,n)A1=A11+A12; %对⾓元素A2=sparse(1:n-1,2:n,0.5,n,n); %输⼊A 的上次对⾓元素 A3=sparse(2:n,1:n-1,0.5,n,n); %输⼊A 的下次对⾓元素 A=A1+A2+A3; b1=0.5*ones(n,1); %b 的除第⼀项元素 b2=sparse(1,1,18,n,1); %b 的第⼀项元素 b=b1-b2; R=A\b输出结果如下:R =26.000017.0000 9.0000 2.0000 -4.0000 -9.0000 -13.0000 -16.0000 -18.0000 -19.0000所以各阻值为(R 1,R 2,…,R 10)=(26,17,9,2,-4,-9,-13,-16,-18,-19)总电阻R 0(即输⼊等效电阻)为00I V R =,⼜nI得到 )(6.35.010180Ω=?=R5-6 有5个反应器连接如图5.6(见书),各个Q 表⽰外部输⼊、输出及反应器间的流量(m 3/min ),各个c 表⽰外部输⼊及反应器内某物质的浓度(mg/m 3)。

数值分析 解线性代数方程组的迭代法和非线性方程的求根实验

数值分析 解线性代数方程组的迭代法和非线性方程的求根实验

华北科技学院上机报告系(部)基础部专业、班级计算B111姓名学号课程名称数值分析上机题目解线性代数方程组的迭代法和非线性方程的求根实验任课教师李慧指导教师李慧成绩(优、良、中、及格、不及格)华北科技学院基础部一. 实验目的1. 熟悉求解线性方程组的有关理论和方法。

会利用雅可比迭代和高斯—塞得尔迭代法求解线性方程组。

通过实际计算,进一步了解各算法的优缺点,选择合适的数值方法。

2.通过对二分法与牛顿迭代法作编程练习与上机运算,进一步体会二分法与牛顿迭代法的不同特点;3.编写割线迭代法的程序,求非线性方程的解,并与牛顿迭代法作比较。

二.实验环境及要求MATLAB 软件等。

三.实验学时2学时四.实验内容1. 1.解线性方程组的迭代法:雅可比迭代和高斯—塞得尔迭代法;2.非线性方程的求根:二分法、牛顿法。

五.实验题及结果1.用雅克比迭代法解方程组⎪⎩⎪⎨⎧=--=-+=++16321382825321321321x x x x x x x x x程序:Jacobi.m 文件:function solution=Jacobi(Ab,epsilon)if nargin==1disp('请输入精度epsilon');returnendrow=size(Ab,1);col=size(Ab,2);if ~ismatrix(Ab)||col-row~=1disp('矩阵大小有误,不能使用Jacobi 迭代法')returnendA=Ab(:,1:row);b=Ab(:,col);ddet=abs(det(A));ddiag=abs(det(diag(diag(A))));if ddet<eps||ddiag<epsdisp('该方程的系数矩阵行列式为0,方程组无解或有无穷多解,或系数矩阵对角线有零元,不能使用Jacobi 迭代法');returnendU=-triu(A,1);L=-tril(A,-1);Dinv=diag(diag(A).^(-1));B=Dinv*(L+U);f=Dinv*b;if max(abs(eig(B)))>=1disp('迭代法不收敛!!!');returnenderror=10;start=zeros(row,1);xk=start;xknext=B*xk+f;error=norm(xk-xknext);while error>epsilonxk=xknext;xknext=B*xk+f;error=norm(xk-xknext);endsolution=xknext;Jacobi1.m 文件:A=[5,2,1,8;2,8,-3,21;1,-3,-6,1];Jacobi(A,0.0001)结果:>> Jacobi1ans =1.00002.0000-1.0000可知结果为:X1=1.0000,X2=2.0000,X3=-1.0000;2. 用高斯-塞德尔迭代法解方程组⎪⎩⎪⎨⎧=++=-+=+-3631223163332111420322318x x x x x x x x x程序:GaussSeidel.m函数文件:function solution=GaussSeidel(Ab,epsilon)if nargin==1disp('请输入精度要求epsilon')returnendrow=size(Ab,1);col=size(Ab,2);if ~ismatrix(Ab)||col-row~=1disp('矩阵大小有误,不能使用高斯塞德尔迭代法')returnendA=Ab(:,1:row);b=Ab(:,col);ddet=abs(det(A));ddiag=abs(det(diag(diag(A))));if ddet<eps || ddiag<epsdisp('该方程的系数矩阵行列式为零,方程组无解或有无穷多解,或系数矩阵的对角线有零元,不能使用高斯塞德尔迭代')returnendU=-triu(A,1);L=-tril(A,-1);D=diag(diag(A));G=(D-L)^(-1)*U;if max(abs(eig(G)))>=1disp('迭代法不收敛!!!')returnenderror=10;n=row;start=zeros(row,1);xk=start;xknext=start;while error>epsilonxk=xknext;xknext(1)=1/A(1,1)*(b(1)-sum(A(1,2:n).*xk(2:n)'));for i=2:n-1Ssum1=sum(A(i,1:i-1).*xk(1:i-1)');Ssum2=sum(A(i,i+1:n).*xknext(i+1:n)');xknext(i)=1/A(i,i)*(b(i)-Ssum1-Ssum2);endxknext(n)=1/A(n,n)*(b(n)-sum(A(n,1:n-1).*xknext(1:n-1)'));error=norm(xk-xknext);endsolution=xknext;GaussSeidel1.m 文件:A=[8.0,-3.0,2.0,20.0; 4.0,11.0,-1.0,33.0;6.0,3.0,12.0,36.0];GaussSeidel(A,0.0001)结果:>> GaussSeidel1ans =3.00002.00001.0000其中X1=3.0000,X2=2.0000,X3=1.00003. 用二分法求方程033)(23=--+=x x x x f 在1.5附近的根。

数值线性代数大作业报告

数值线性代数大作业报告

数值线性代数实验大报告指导老师:赵国忠姓名: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_线性代数方程组的数值解法

实验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计算线性代数方程组的数值解法

数学实验报告——利用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 黄浩 43一、实验目的1.学会用MATLAB软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2.通过实例学习用线性代数方程组解决简化的实际问题。

二、For personal use only in study and research; not forcommercial use三、四、实验内容1.《数学实验》第二版(问题1)问题叙述:通过求解线性方程组,理解条件数的意义和方程组性态对解的影响,其中是n阶范德蒙矩阵,即是n阶希尔伯特矩阵,b1,b2分别是的行和。

(1)编程构造(可直接用命令产生)和b1,b2;你能预先知道方程组和的解吗?令n=5,用左除命令求解(用预先知道的解可验证程序)。

(2)令n=5,7,9,…,计算和的条件数。

为观察他们是否病态,做以下试验:b1,b2不变,和的元素,分别加扰动后求解;和不变,b1,b2的分量b1(n),b2(n)分别加扰动后求解。

分析A与b的微小扰动对解的影响。

取10^-10,10^-8,10^-6。

(3)经扰动得到的解记做,计算误差,与用条件数估计的误差相比较。

模型转换及实验过程:(1)小题.由b1,b2为,的行和,可知方程组和的精确解均为n 行全1的列向量。

在n=5的情况下,用matlab编程(程序见四.1),构造,和b1,b2,使用高斯消去法得到的解x1,x2及其相对误差e1,e2(使用excel计算而得)为:由上表可见,当n=5时,所得的解都接近真值,误差在10^-12的量级左右。

(2)小题分别取n=5,7,9,11,13,15,计算和的条件数c1和c2,(程序见四.2),结果如下:由上表可见,二者的条件数都比较大,可能是病态的。

为证实和是否为病态,先保持b不变,对做扰动,得到该情况下的高斯消元解,(程序见四.3),结果如下:(为使结果清晰简洁,在此仅列出n=5,9,13的情况,n=7,11,15略去)=10^-10时:=10^-8时:=10^-6时:由上表可见:a)对于希尔伯特阵,随着阶数的增加,微小扰动对解带来的影响越来越大,到了n=9时,已经有了6倍误差的解,到了n=13时,甚至出现了22倍误差的解元素;而随着的增加,解的偏差似乎也有增加的趋势,但仅凭上述表格无法具体判断(在下一小题中具体叙述)。

数学实验 5:线性代数方程组的数值解法

数学实验 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. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

教学科研实验项目表
课程名称:数值线性代数实验层次:(√本科、研究生、科研)专业方向:信息与计算科学 2005年4月 7日
,编制不选主元、列主元、全主元Gauss
通用子程序。

比较三种算法的计算结果。

求解正定方程的Chol esky 分解算法
掌握Chol esky 三角分解的编程技巧。

与第一个实验结果比较,认识求解正定方程不必要选取主元。

(1) b 随机选取,(2)系数矩阵为100阶矩阵
(3)系数矩阵为40阶Hilbert 矩阵;b 的第i 个分量
编写算法子程序。

用列主元Gauss 算本实验方程组。

比较两种算法的实验结果,评论各个算法的优劣。

(1)求前两个实验中的三个线性方程组,并根据计算结果,评价各种算法的优劣。

(2)求
使在残向量的2范数最
小的意义下拟合下面的数据
-1.00 -0.75 -0.50 0.00
0.25 0.50
0.75
1.00
0.8125
0.75
1.00
1.3125
1.75
2.3125
编制通用子程序。

分析计算结果。

考虑两点边值问题
的数值解。

将 [0,1] 区间n等分,令
,,得到差分方程用Jacobi、G-S、
超松驰迭代法编
程计算数值解,
并与精确解
+
进行比较,给出
精度分析。

4
考虑如下的Dirichlet问题
其中为正方形区域的边界。

其差分方程编制算法程序,分析计算结果。

计算50阶到100阶三对角矩阵的特征值和特征向量编制
法、二分法通用子程序。

用Jacobi
阵的全部特征值和特征向量。

用二分法求其最大和最小特征值,再利用反幂法求对应特征向量。

相关文档
最新文档