科学计算—理论、方法
及其基于MATLAB 的程序实现与分析
三、 解线性方程组的迭代法（Iteration ）
线性方程组的理论求解公式
b
A x 1
-= （１）
在应用于实际问题的计算时，通常面临两方面的问题 １、计算过程复杂， ２、不能保证算法的稳定性；
此外，当初始数据（可能）存在误差时，按公式（１）即使求出了“精确解”意义也不大，因此，对于存在初始数据误差、特别是大型的线性方程组求解，需要寻求能达到精度要求的、操作和计算过程相对简单的求解方法。下面将要介绍的迭代法就属于这类方法。
迭代法求解线性方程组的基本思想是
１） 不追求“一下子”得到方程组的解，而是在逐步逼近方程组的精
确解的迭代过程中获得满足精度要求的近似解，这一点与直接法不同；
２） 通过对问题的转化，避免（困难的）矩阵求逆运算。
用迭代法求解线性方程组，首先要把线性方程组写成等价的形式
f Mx x b Ax +=?= （２）
式（２）的右端称为迭代格式，由迭代格式（２）确定如下的迭代算法：
,2,101=?
??∈?+=+k R x f
Mx x n
k k （３）
对于给定的线性方程组，可以写成不同的（无穷多）迭代格式，有意义的（可用的）迭代格式应具有收敛性―生成的解向量序列{}x n 收敛于方程组的解；而好的迭代法应具有较高的收敛速度。
关于迭代法收敛性的两个判别条件：
a 、充分必要条件是：矩阵M 的谱半径
(){}1,,2,1max <==n i M i
i
λρ
b 、充分条件是：矩阵M 的某个算子
范数
M <1。
设x 是方程组（２）的解，{}m x 是迭代法（３）生成的任一序列，因为
f Mx x +=，f Mx x m m +=+1
所以
()()()022
1x x M
x x M
x x M x x m
m m m -
==-
=-=--- （４）
设1
1--=?=TJT
M J MT
T ，其中矩阵J 是矩阵M 的Jordan 标准型，那么容易验证1
-=T
TJ M m m
，并且
()[
]()
()()1
lim ,2,10lim lim 0
lim lim 01
0
?=-
?=+∞
→+∞
→-+∞
→+∞→+∞
→M M
n i x x T
J
T x x M
x x m
m m
i n m
m m
m m m ρρλ （５）
此外，因为
(
)0
2
2
11x x M
x x M x x M x x M x x m m m m m -≤≤-≤-≤-=---- （６）
所以
x x x x M
M m m m m m
m =?=-?=?<+∞
→+∞
→+∞
→lim 0lim
0lim
1
（７）
注：迭代格式（２）所确定的迭代法收敛与否，完全由系数矩阵M 决定，而与常数项f 无关．
常用的迭代法
１、Jacobian 迭代法:
U
L D a a a a a a a a a a a a a a a a a a A n n n nn n nn nn n n n
n --=??
?
??
?
??????----????????????--
--????????????=?????
?
???
???=--00
0000000
000011121
1
2122112
122221
11211
()()
?
??=+=?++=???
?--==----+b D f U L D M b D x U L D x U L D A b Ax m m 1
11
11 （８）
例1 解下面方程组（精确解为T x )1,1,1(*=）.
???
??=++-=+-=++.
14103,53102,
14310321
321321x x x x x x x x x
解 1） 改写成等价形式
???
?
?
?
???--=++=----=--=).314(101),325(101)325(101),314(101213213
12321x x x x x x x x x x x 2） 构造迭代公式，即为雅可比迭代公式
???
?
??
???=--=++=--=+++.,2,1,0),314(101),325(101),314(101)
(2)(1)1(3)
(2)(1)1(2
)
(3)(2)1(1 k x x x x x x x x x k k k k k k k k k 3） 取初始向量T x )0,0,0()0(=，即,0)0(3)0(2)0(1===x x x 代入上式，求出
4.110
14,5.010
5,4.110
14)
1(3)
1(2
)
1(1
==
===
x x x .
再代回公式中，求出
11
.1)4.15.0314(101)
2(1
=-?-=
x ， 2
.1)5.034.125(101)
2(2
=?+?+=x ，
11.1)2.1311.114(10
1)
2(3
=?--=
x .
依次迭代，计算结果如表4-1.
Example intera_j.m
itera_j
2、Gauss-Seidel 迭代法:
()()()()?
??-=-=?-+-=???
?--==----+b L D f U L D M b L D Ux L D x U L D A b Ax n n 1
1
1
11 （9）
根据GS 迭代法（9），可进一步得到
()()()b
D Ux Lx D
x b Ux x L D b
L D Ux L D x n n n n n n n 1
11
111
1
1)(-+-++--+++=?+=-?-+-= （10）
即
???
???????????+???
????
?????
?????????---+??????? ???????????
??
?--
-=???????
??-+++--+++b x x x a a a x x x a a a D x x x N
n
n n n n n N n n n nn n N n n n
2
1
111212111
1
1
21112
11100
00000000
（11）
式（11）表明：Gauss-Seidel 迭代法在计算第k 个迭代值k n x 1
+时，及时
地利用了在此步迭代中得到的新的迭代值：i
n x
1
+，1,,2,1-=k i ，由于第
1+n 步的迭代值通常比第n 步的迭代值更接近方程组的精确解，所以，
在Jacobian 迭代法和Gauss-Seidel 迭代法都收敛的情况下，Gauss-Seidel 迭代法的收敛速度比Jacobian 迭代法的收敛速度快。
例2 解下面方程组（与例1相同，精确解为T x )1,1,1(*=）.
???
??=++-=+-=++.
14103,53102,14310321
321321x x x x x x x x x 解 1） 原方程组改为等价方程组
???
?
??
???
--=++=--=).314(101),325(101),314(101213212
321x x x x x x x x x . 2） 构造迭代公式，即为高斯-赛德尔迭代公式
???
?
??
???=--=++=--=++++++.,2,1,0),314(101),325(101),314(10
1)
1(2)1(1)1(3
)(2)1(1)1(2
)
(3)(2)1(1
k x x x x x x x x x k k k k k k k k k . 3） 取初始向量T x )0,0,0()0(=，即,0)0(3)0(2)0(1===x x x 代入上式，求出
4
.1)00314(10
1)
1(1
=-?-=x ， 78
.0)034.125(101)
1(2=?+?+=
x ， 026
.1)78.034.114(10
1
)
1(3=?--=
x .
迭代计算下去，得表4-2.
Example itera_gs.m
itera_gs
对Gauss-Seidel 迭代法做进一步的研究，式（９）还可以写成如下的形式：
()()()()
()
n n J
n S
G n n n n n n n n n n n n n n x x L D
x x x x L D
b D x U L D x b D Lx Ux Lx Lx D x b Ux x L D b
L D Ux L D x -+=?-+++=?+-++=?+=-?-+-=+-++---+-+-++--+11
_1
11
1
111
11111
1
1)()( （12）
即Gauss-Seidel 迭代法是在Jacobian 迭代法的基础上增加了一个修正项：
()n n x x L D
-+-11
并且这个修正项起到了加速的作用，受这一事实的启发，为获得更快的收敛速度，我们对Gauss-Seidel 迭代法再做进一步的研究
()()()()()()
b x D U Lx D
x x b x D U Lx Dx Dx b
Ux x L D b
L D Ux L D x n n n n n n n n n n n n +-++=?+-++=?+=-?-+-=+-++++--+11
11111
1
1 （13）
式（13）表明Gauss-Seidel 迭代法的第1+n 步迭代值是在第n 步迭代值的基础上增加了一个修正项：
()()b x D U Lx D
n n +-++-11
（14）
并且这个修正项使第1+n 步迭代值更接近方程组的精确解，受这一事实的启发，我们希望通过用一个适当的因子乘修正项（14）的办法达到获得更快收敛速度的目的：
()()
()()()()[]()()[]()b
wL
D w x wU D
w wL
D x wb
x wU D w x wL D b x D U Lx w Dx Dx b x D U Lx wD x x n n n n n n n n n n n n 11111111
111--+++++-+-++--=?++-=-?
+-++=?+-++= （15）
从而得到
3、SOR （Successive Over Relaxation Method ） 迭代法:
()
()[]()
x D w L w D w U x w D w L b n n
+--=--++-11
1
1 （16）
其中w 称为松弛因子，由于
()()[]
()()()?
?
?
??
?
?????
?------?????
??
???
??=
+--=----nn n n n nn nn n a w wa a w wa wa a w a wa
wa a wa a wU
D
w wL D M 10
1010001122
112
111
1
1
2221111
（17） ()()()()∏∏∏====
-=-=
?n
k k
n
n
k kk
n
k kk
n M w a
a w M 1
1
111det λ （1８）
()
()
1
12
0111
>??>?
>∨<-?∏
=M
M
w w w
k n
k k n
λλ （1９） 所以，为保证SOR 迭代法的收敛性，要求2
0<<
w 。
例3 用超松驰迭代法解下面方程组，取松驰因子4.1=ω.
?????
???????=?????????????????
???????------010
121
1210012100124321x x x x .
解 方程组的精确解为：T x )8.0,6.1,4.1,2.1(=.取初值T x )1,1,1,1()0(=，用高斯-赛德尔迭代10次，得T x )798216.0,5964336.1,3955917.1,1966324.1()10(=.
利用SOR 方法，构造迭代公式
??
???
???
??
?-+
=+-++=+-+=+-+=+++++++).
2(2
),
21(2
),
2(2
),
21(2
)
(4)
1(3
)
(4
)1(4)
(4)
(3
)
1(2)
(3
)
1(3)
(3)
(2)
1(1
)
(2)1(2)
(2)
(1)
(1
)1(1k k k k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x x x x ωωωω
与高斯-赛德尔方法相同，初值为T x )1,1,1,1()0(=.迭代计算结果列于表4-3.
表4-3 SOR 迭代的数值结果
k
)
(1
k x
)
(2
k x
)
(3
k x
)
(1
k x
1 2 3 4 5
1 1
1.343 1.19545
1.203472 1
1.49 1.4753 1.40236
1.4028735 1.7 1.616 1.65095 1.6019815
1.5939059
0.79 0.8152 0.829585 0.789553 0.7999129
由表4-3可知，用超松驰迭代法只迭代了5次，结果与G-S 法迭代10次的结果大体相同，可见，SOR 方法的松驰因子起到了加速收敛的重要作用.
Example itera_sor.m
itera_sor
关于迭代法收敛性的两个判别条件：
a 、充分必要条件是：矩阵M 的谱半径
(){}1,,2,1max <==n i M i
i
λρ
b 、充分条件是：矩阵M
的某个算子范数
M <1。
特例：对于下面方程组，证明：雅可比迭代法收敛而高斯-赛德尔迭代法发散.
????
?
?????=??????????????
?
?????-11112
2
111221321X X X . 证明 1) 对于雅可比迭代矩阵
???
?
?
?????-----=+=-02
2
101
220)(1
U L D B J ， J
B 的特征方程为
02
2
11
22
)det(=-=-λ
λλ
λJ B I .
J
B 的特征多项式0)(3==λλf ，所以0=λ为J B 的特征根，显然10)(<=J B ρ,因
此，由迭代收敛基本定理可知雅可比迭代法收敛.
2) 对于高斯-赛德尔迭代法，其迭代矩阵为
???
?
?
?????---=-=--20
320
220
)(1
U L D B S
G ， S G B -的特征方程为
0)2(2
320
2
2)det(2
=+=+-+-=--λλλλλ
λS G B I .
显然，12)(>=-S G B ρ，因此高斯-赛德尔迭代法发散.
关于收敛性重要结论： 设线性方程组b
Ax
=.
（1）A为严格对角占优矩阵，则解b
Ax=的雅可比迭代过程和高斯-赛德尔迭代过程均收敛.
（2）若A是对称正定方阵，则解b
Ax=的高斯-赛德尔迭代过程收敛.
（3）A为严格对角占优矩阵，1
<ω，则解b
0≤
Ax=的SOR迭代过程收敛.
数值分析实习报告 姓名：gestepoA 学号：201******* 班级：***班
序言 随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。 由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。它具有以下优点： 1友好的工作平台和编程环境。MATLAB界面精致，人机交互性强，操作简单。 2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M 文件）后再一起运行。 3强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。 4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。
目录 1 第一题 (4) 1.1 实验目的 (4) 1.2 实验原理和方法 (4) 1.3 实验结果 (5) 1.3.1 最佳平方逼近法 (5) 1.3.2 拉格朗日插值法 (7) 1.3.3 对比 (8) 2 第二题 (9) 2.1实验目的 (9) 2.2 实验原理和方法 (10) 2.3 实验结果 (10) 2.3.1 第一问 (10) 2.3.2 第二问 (11) 2.3.3 第三问 (11) 3 第三题 (12) 3.1实验目的 (12) 3.2 实验原理和方法 (12) 3.3 实验结果 (12) 4 MATLAB程序 (14)
MATLAB数值计算 MATLAB数值计算 (1) 1创建矩阵 (3) 1.1直接输入 (3) 1.2向量 (3) 1.2.1linspace：线性分布 (3) 1.2.2冒号法 (3) 1.3函数创建 (4) 1.3.1eye:单位矩阵 (4) 1.3.2rand：随机矩阵 (4)
1.3.3zeros：全0矩阵 (4) 1.3.4ones：全1矩阵 (5) 2矩阵运算 (5) 2.1加减 (5) 2.1.1[M×N]±[M×N] (5) 2.2乘 (6) 2.2.1[M×N]*a (6) 2.2.2[M×N]*[N×M] (6) 2.3乘方 (7) 2.3.1[M×M]^a (7) 2.3.2a^[M×M] (7) 2.4特殊运算 (8) 2.4.1求逆inv (8) 2.4.2行列式det (8) 2.4.3特征值eig (8) 2.4.4转置'和.' (9) 2.4.5变形reshape (10) 2.4.6翻转rot90,fliplr,flipud (11) 2.4.7抽取diag,tril,triu (12) 2.5数组运算 (12) 2.5.1乘 (12) [M×N].*[M×N] (12) 2.5.2除 (13) [M×N]./[M×N] (14) [M×N].\[M×N] (14) 2.5.3乘方 (14) [M×N].^[M×N] (15) a.^[M×N] (15) 2.6除法 (15) 2.6.1求解线性方程组 (15) 3多项式 (16) 3.1系数表示法poly (16) 3.2求根roots (16) 3.3乘法conv (16) 3.4除法deconv (17) 3.5求值polyval (17) 3.6微分polyder (18)
科学计算—理论、方法 及其基于MATLAB 的程序实现与分析 三、 解线性方程组的迭代法（Iteration ） 线性方程组的理论求解公式 b A x 1 -= （１） 在应用于实际问题的计算时，通常面临两方面的问题 １、计算过程复杂， ２、不能保证算法的稳定性； 此外，当初始数据（可能）存在误差时，按公式（１）即使求出了“精确解”意义也不大，因此，对于存在初始数据误差、特别是大型的线性方程组求解，需要寻求能达到精度要求的、操作和计算过程相对简单的求解方法。下面将要介绍的迭代法就属于这类方法。 迭代法求解线性方程组的基本思想是 １） 不追求“一下子”得到方程组的解，而是在逐步逼近方程组的精 确解的迭代过程中获得满足精度要求的近似解，这一点与直接法不同； ２） 通过对问题的转化，避免（困难的）矩阵求逆运算。 用迭代法求解线性方程组，首先要把线性方程组写成等价的形式 f Mx x b Ax +=?= （２） 式（２）的右端称为迭代格式，由迭代格式（２）确定如下的迭代算法： ,2,101=? ??∈?+=+k R x f Mx x n k k （３）
对于给定的线性方程组，可以写成不同的（无穷多）迭代格式，有意义的（可用的）迭代格式应具有收敛性―生成的解向量序列{}x n 收敛于方程组的解；而好的迭代法应具有较高的收敛速度。 关于迭代法收敛性的两个判别条件： a 、充分必要条件是：矩阵M 的谱半径 (){}1,,2,1max <==n i M i i λρ b 、充分条件是：矩阵M 的某个算子 范数 M <1。 设x 是方程组（２）的解，{}m x 是迭代法（３）生成的任一序列，因为 f Mx x +=，f Mx x m m +=+1 所以 ()()()022 1x x M x x M x x M x x m m m m - ==- =-=--- （４） 设1 1--=?=TJT M J MT T ，其中矩阵J 是矩阵M 的Jordan 标准型，那么容易验证1 -=T TJ M m m ，并且 ()[ ]() ()()1 lim ,2,10lim lim 0 lim lim 01 0
电子科技大学电子工程学院标准实验报告（实验）课程名称MATLAB与数值分析 学生姓名：李培睿 学号：2013020904026 指导教师：程建
一、实验名称 《MATLAB与数值分析》第一次上机实验 二、实验目的 1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算 操作。（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序） 2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号 转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。（用.m文件编写进行符号因式分解和函数求反的程序） 3. 掌握Matlab函数的编写规范。 4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、 三维曲线和面的填充、三维等高线等。（用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释） 5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。 三、实验内容 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以x， y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图，允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort 函数。 四、实验数据及结果分析 题目一： ①在Editor窗口编写函数代码如下：
Matlab作业3（数值分析） 机电工程学院（院、系）专业班组 学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。 答： 2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。（2）求解在x=8时多项 式(x-1)(x-2) (x-3)(x-4)的值。 答：（1） （2）
3. y=sin(x)，x从0到2π，?x=0.02π，求y的最大值、最小值、均值和标准差。 4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。解：x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909
5.实验数据处理：已知某压力传感器的测试数据如下表 p为压力值，u为电压值，试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟 合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。解： >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267
第2章牛顿插值法实现 参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序： function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; （1）保存名为newpoly.m 的M 文件 （2）输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 *
1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组， 1231231234748212515 x x x x x x x x x -+=?? -+=-??-++=? （1）给出Jacobi 迭代法的迭代方程，并判定Jacobi 迭代法求解此方程组是否收敛。 （2）若收敛，编程求解该线性方程组。 解（1）：A=[4 -1 1;4 -8 1;-2 1 5] %线性方程组系数矩阵 A = 4 -1 1 4 -8 1 -2 1 5 >> D=diag(diag(A)) D = 4 0 0 0 -8 0 0 0 5 >> L=-tril(A,-1) % A 的下三角矩阵 L = 0 0 0 -4 0 0 2 -1 0 >> U=-triu(A,1) % A 的上三角矩阵 U = 0 1 -1 0 0 -1 0 0 0 B=inv(D)*(L+U) % B 为雅可比迭代矩阵 B = 0 0.2500 -0.2500 0.5000 0 0.1250 0.4000 -0.2000 0 >> r=eigs(B,1) %B 的谱半径
r = 0.3347 < 1 Jacobi迭代法收敛。 (2)在matlab上编写程序如下： A=[4 -1 1;4 -8 1;-2 1 5]; >> b=[7 -21 15]'; >> x0=[0 0 0]'; >> [x,k]=jacobi(A,b,x0,1e-7) x = 2.0000 4.0000 3.0000 k = 17 附jacobi迭代法的matlab程序如下： function [x,k]=jacobi(A,b,x0,eps) % 采用Jacobi迭代法求Ax=b的解 % A为系数矩阵 % b为常数向量 % x0为迭代初始向量 % eps为解的精度控制 max1= 300; %默认最多迭代300，超过300次给出警告D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; x=B*x0+f; k=1; %迭代次数 while norm(x-x0)>=eps x0=x; x=B*x0+f; k=k+1; if(k>=max1) disp('迭代超过300次，方程组可能不收敛'); return; end end
第6章 MATLAB 数值计算 例6.1 求矩阵A 的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元素。 1356 78256323578255631 01-???? -? ?=???? -??A A=[13,-56,78;25,63,-235;78,25,563;1,0,-1]; max(A,[],2) %求每行最大元素 min(A,[],2) %求每行最小元素 max(A) %求每列最大元素 min(A) %求每列最小元素 max(max(A)) %求整个矩阵的最大元素。也可使用命令：max(A(:)) min(min(A)) %求整个矩阵的最小元素。也可使用命令：min(A(:)) 例6.2 求矩阵A 的每行元素的乘积和全部元素的乘积。 A=[1,2,3,4;5,6,7,8;9,10,11,12]; S=prod(A,2) prod(S) %求A 的全部元素的乘积。也可以使用命令prod(A(:)) 例6.3 求向量X =(1！，2！，3！，…，10！)。 X=cumprod(1:10) 例6.4 对二维矩阵x ，从不同维方向求出其标准方差。 x=[4,5,6;1,4,8] %产生一个二维矩阵x y1=std(x,0,1) y2=std(x,1,1) y3=std(x,0,2) y4=std(x,1,2) 例6.5 生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。 X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X)
例6.6 对下列矩阵做各种排序。 185412613713-?? ??=?? ??-?? A A=[1,-8,5;4,12,6;13,7,-13]; sort(A) %对A 的每列按升序排序 -sort(-A,2) %对A 的每行按降序排序 [X,I]=sort(A) %对A 按列排序，并将每个元素所在行号送矩阵I 例6.7 给出概率积分 2 (d x x f x x -? e 的数据表如表6.1所示，用不同的插值方法计算f (0.472)。 x=0.46:0.01:0.49; %给出x ，f(x) f=[0.4846555,0.4937542,0.5027498,0.5116683]; format long interp1(x,f,0.472) %用默认方法，即线性插值方法计算f(x) interp1(x,f,0.472,'nearest') %用最近点插值方法计算f(x) interp1(x,f,0.472,'spline') %用3次样条插值方法计算f(x) interp1(x,f,0.472,'cubic') %用3次多项式插值方法计算f(x) format short 例6.8 某检测参数f 随时间t 的采样结果如表6.2，用数据插值法计算t =2，7，12，17，22，17，32，37，42，47，52，57时的f 值。 T=0:5:65; X=2:5:57;
列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);
end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')
Matlab实现数值分析插值及积分 摘要： 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。 分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。 其中Aitken插值计算的结果图如下： 对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下：
问题重述 问题一：已知列表函数 表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。 问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。 问题解决 问题一：插值方法 对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。 一、拉格朗日插值法： 拉格朗日插值多项式如下： 首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只 可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得：
MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解： x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*')
2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数 10 的近似根，并写出调用方式： 精度为10 解： >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8
MATLAB数值计算 (读书日记及程序编写) 第四章方程求根 (2)
第四章 方程求根 #二分法 求2的值 转化成方程02-2 =x 最慢的方法是取初值1001=x 02-21>x ,取502=x 这样得到 也可以x0=a, x1=x0+h, 进行扫描，若f(x0)*f(x1)<0, 则扫描成功，有根区间为[x0,x1],否则继续扫描，如果出现x1>b ，表面扫描失败，再缩小步长h, 再次扫描。 >> format long %让显示的值为 M=2,a=1,b=2,k=0; while b-a>eps x=(a+b)/2; if x^2>M b=x else a=x end k=k+1 end 执行后得到的值为： k = 50 b = 1.414213562373095 k = 51 b = 1.414213562373095 k = 52 最后得到的值就是Matlab 能表达的最接近的值。 #牛顿法
求解f(x)=0的牛顿法是在f(x)画一条切线，确定切线与x 轴的焦点，通过迭代 ) (x f )f(x -n n 1'=+n n x x 对于平方根的问题，牛顿法简洁有效， 换成f(x)=x^2-M, )(x f n ' =2x 这样 ??? ? ??+==+n n n n x M x M x x 212x -x -n 2n 1 该算法就是反复求x 和M/x 的平均值，Matlab 的程序为： format long %让显示的值为 xprev=2; %取的不等于初值x 的一个值，让判断能继续 x=100; %取的初值为3 while abs(x-xprev)>eps*abs(x) xprev=x; x=0.5*(x+2/x) end x = 1.833333333333333 x = 1.462121212121212 x = 1.414998429894803 x = 1.414213780047198 x = 1.414213562373112 x = 1.414213562373095 x = 1.414213562373095 可见6步很快就收敛 然而，若f(x)不具有连续的、有界的一阶、二阶导数，牛顿法的收敛将变得很慢。 #fzero 函数直接求根 求x^3-1在区间[0,10]上的根 fzero(@(x)x^3-1,[0,10]) ans = 1 fzerogui(@(x)x^3-1,[0,10]) 可以通过在图形界面上选择割点来得到
1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6) disp('牛顿法'); i=1; n0=180; p0=pi/3; tol=10^(-6); for i=1:n0 p=p0-(p0-sin(p0))/(1-cos(p0)); if abs(p-p0)<=10^(-6) disp('用牛顿法求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次牛顿迭代后无法求出方程的解') end 2、disp('Steffensen加速'); p0=pi/3; for i=1:n0 p1=0.5*p0+0.5*cos(p0); p2=0.5*p1+0.5*cos(p1); p=p0-((p1-p0).^2)./(p2-2.*p1+p0); if abs(p-p0)<=10^(-6) disp('用Steffensen加速求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次Steffensen加速后无法求出方程的解') end 1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根， %误差限为 e=10^-4 disp('二分法')
a=0.2;b=0.26; tol=0.0001; n0=10; fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1; for i=1:n0 p=(a+b)/2; fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1; if fp==0||(abs((b-a)/2) roots([1 -1 -1]) x=linspace(0,2*pi,10); y=sin(x); xi=linspace(0,2*pi,100); y1=interp1(x,y,xi); y2=interp1(x,y,xi,'spline'); y3=interp1(x,y,xi,'cublic'); plot(x,y,'o',xi,y1,xi,y2,xi,y3) x=[0 300 600 1000 1500 2000]; y=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; xi=linspace(0,2000,20); yi=1.0332*exp(-(xi+500)/7756); y1=interp1(x,y,xi,'spline'); subplot(2,1,1);plot(x,y,'o',xi,yi,xi,y1,'*') p=polyfit(x,y,2); y2=polyval(p,xi); subplot(2,1,2);plot(x,y,'o',xi,yi,xi,y2,'*') x=[0 300 600 1000 1500 2000]; y=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; xi=linspace(0,2000,20); y1=interp1(x,y,xi,'spline'); subplot(2,1,1);plot(x,y,'-o', xi,y1,'-*') p=polyfit(x,y,2); y2=polyval(p,xi); subplot(2,1,2);plot(x,y,'-o',xi,y2,'-*') 第二讲MATLAB的数值分析 2-1矩阵运算与数组运算 矩阵运算和数组运算是MATLAB数值运算的两大类型，矩阵运算是按矩阵的运算规则进行的，而数组运算则是按数组元素逐一进行的。因此，在进行某些运算（如乘、除）时，矩阵运算和数组运算有着较大的差别。在MATLAB中，可以对矩阵进行数组运算，这时是把矩阵视为数组，运算按数组的运算规则。也可以对数组进行矩阵运算，这时是把数组视为矩阵，运算按矩阵的运算规则进行。 1、矩阵加减与数组加减 矩阵加减与数组加减运算效果一致，运算符也相同，可分为两种情况： （1）若参与运算的两矩阵（数组）的维数相同，则加减运算的结果是将两矩阵的对应元素进行加减，如 A=[1 1 1;2 2 2;3 3 3]; B=A; A+B ans= 2 2 2 4 4 4 6 6 6 （2）若参与运算的两矩阵之一为标量（1*1的矩阵），则加减运算的结果是将矩阵（数组）的每一元素与该标量逐一相加减，如 A=[1 1 1;2 2 2;3 3 3]; A+2 ans= 3 3 3 4 4 4 5 5 5 2、矩阵乘与数组乘 （1）矩阵乘 矩阵乘与数组乘有着较大差别，运算结果也完全不同。矩阵乘的运算符为“*”，运算是按矩阵的乘法规则进行，即参与乘运算的两矩阵的内维必须相同。设A、B为参与乘运算的 =A m×k B k×n。因此，参与运两矩阵，C为A和B的矩阵乘的结果，则它们必须满足关系C m ×n 算的两矩阵的顺序不能任意调换，因为A*B和B*A计算结果很可能是完全不一样的。如：A=[1 1 1;2 2 2;3 3 3]; B=A; A*B ans= 6 6 6 12 12 12 18 18 18 F=ones(1,3); G=ones(3,1); F*G ans 3 G*F ans= 1 1 1 1 1 1 1 1 1 （2）数组乘 数组乘的运算符为“.*”，运算符中的点号不能遗漏，也不能随意加空格符。参加数组乘运算的两数组的大小必须相等（即同维数组）。数组乘的结果是将两同维数组（矩阵）的对应元素逐一相乘，因此，A.*B和B.*A的计算结果是完全相同的，如： A=[1 1 1 1 1;2 2 2 2 2;3 3 3 3 3]; B=A; A.*B ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 B.*A ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 由于矩阵运算和数组运算的差异，能进行数组乘运算的两矩阵，不一定能进行矩阵乘运算。如 A=ones(1,3); B=A; A.*B ans= 1 1 1 A*A ???Error using= => 数值分析MATLAB计算实验报告 姓名班级学号 一、实验名称 用MATLAB编程实现数值微分的外推法计算。 二、实验目的 1．掌握数值微分和定义和外推法的计算过程； 2．了解数值微分外推法的计算方法并且编写出与其算法对应的MATLAB程序代码；3．体会利用MATLAB软件进行数值计算。 三、实验内容 用外推法计算f(x)=x2e?x在x=0.5的导数。 四、算法描述 1.命名函数。 2.如果输入未知数少于四个，默认精度10^-3 3.描述T表矩阵坐标 4.依次赋值计算 T表第一列 5.根据数值微分计算公式求出T表矩阵的值 6.若达到精度则运算结束，若未达到循环计算 7.输出T表，得出的值就是导数值 五、实验结果 六、实验结果分析 此实验通过MATLAB实现外推法数值微分计算，得到相应的数据，方便对数据进行分析。从结果可以看出，当步长h=0.025时用中点微分公式只有3位有效数字，外推一次达到5位有效数字，外推两次达到9位有效数字。 七、附录（程序） function g=waituifa(fname,x,h,e) if nargin<4,e=1e-3; end; i=1; j=1; G(1,1)=(feval(fname,x+h)-feval(fname,x-h))/(2*h); G(i+1,1)=(feval(fname,x+h/2)-feval(fname,x-h/2))/h; G(i+1,j+1)=(4^j*G(i+1,j)-G(i,j))/(4^j-1); while abs(G(i+1,i+1)-G(i+1,i))>e i=i+1; G(i+1,1)=(feval(fname,x+h/2^i)-feval(fname,x-h/2^i))/(2*h/2^i); for j=1:i G(i+1,j+1)=((4^j)*G(i+1,j)-G(i,j))/(4^j-1); end end G g=G(i+1,i+1); 2015.1.9 上机作业题报告 JONMMX 2000 1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果 1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序 数值分析幂法与反幂法 matlab程序 随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。 要求 1）比较不同的原点位移和初值说明收敛性 2）给出迭代结果，生成DOC文件。 3）程序清单，生成M文件。 解答： >> A=rand(5) %随机产生5*5矩阵求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A矩阵和A的转置相加，得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286 B=?? ????? ???? ?? ???6286.13929.03467.15979.08044 .03929.00944 .12144.18506 .13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613 .09173 .04187.1 编写幂法、反幂法程序： function [m,u,index,k]=pow(A,u,ep,it_max) % 求矩阵最大特征值的幂法，其中 % A 为矩阵； % ep 为精度要求，缺省为1e-5； % it_max 为最大迭代次数，缺省为100； % m 为绝对值最大的特征值； % u 为对应最大特征值的特征向量； % index ，当index=1时，迭代成功，当index=0时，迭代失败 if nargin<4 it_max=100; end if nargin<3 ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0.01; % 修改移位参数，原点移位法加速收敛，为0时，即为幂法 I=eye(n) T=A-m0*I while k<=it_max v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1) 1. 使用444（x）对数据进行插值，并写出误差分析理论。 建立脚本 x1=[0.2 0.4 0.6 0.8 1.0]; y1=[0.98 0.92 0.81 0.64 0.38]; n=length(y1); c=y1(:); for j=2:n %求差商 for i=n:-1:j c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1)); end end syms x df d; df(1)=1;d(1)=y1(1); for i=2:n %求牛顿差值多项式 df(i)=df(i-1)*(x-x1(i-1)); d(i)=c(i)*df(i); end disp('4次牛顿插值多项式'); P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,保留小数点后5位数figure ezplot(P4,[0.2,1.08]); 输出结果为 P4 =- 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98 插值余项：R4(x)=f(5)( ξ)/ (5!)* (x - 0.6)*(x - 0.4)*(x - 0.8)*(x - 1)*(x-0.2) 新建一个M-file syms x l; x1=[0.2 0.4 0.6 0.8 1.0]; y1=[0.98 0.92 0.81 0.64 0.38]; n=length(x1); Ls=sym(0); for i=1:n l=sym(y1(i)); for k=1:i-1 l=l*(x-x1(k))/(x1(i)-x1(k)); end for k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k)); end Ls=Ls+l; end Ls=simplify(Ls) %为所求插值多项式Ls(x). figure ezplot(Ls,[0.2,1.08]); 输出结果为 Ls =- (25*x^4)/48 + (5*x^3)/6 - (53*x^2)/48 + (23*x)/120 + 49/50 插值余项：R4(x)=f(5)( ξ)/ (5!)* (x - 0.6)*(x - 0.4)*(x - 0.8)*(x - 1)*(x-0.2) 2. 试求3次、 合曲线，用图示数据曲线及相应的三种拟合曲线。 建立脚本 X=[0.0 0.1 0.2 0.3 0.5 0.8 1.0]; Y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46]; p1=polyfit(X,Y,3) p2=polyfit(X,Y,4)第3章 MATLAB数值计算-习题 答案
第2讲 matlab的数值分析
数值分析第四章外推法计算数值微分MATLAB计算实验报告.docx
东南大学-数值分析上机题作业-MATLAB版
数值分析幂法与反幂法-matlab程序
MATLAB解决数值分析问题
相关主题
文本预览