迭代法与高斯消元法

4. 高斯消元法
矩阵的加、减、乘法 如：real 如：real a(m,n), b(m,n), c(m,n) c=a+b c=a–b 如：real 如：real a(m,n), b(n,L), c(m,L) c = matmul(a,b)
高斯消元法
在中学里，我们已学过用消元法解二元、 在中学里，我们已学过用消元法解二元、三元线性 方程组，现在我们把它推广到n个方程 个方程n个未知量的一般情 方程组，现在我们把它推广到 个方程 个未知量的一般情 即高斯消元法。 形，即高斯消元法。 其基本思想是通过消元变形把方程组化为容易求解的 同解方程组。而且对求解未知元较多的方程组， 同解方程组。而且对求解未知元较多的方程组，此消元步 骤规范而又简便。 骤规范而又简便。 为叙述方便， 为叙述方便，先给出在高斯消元过程中常用到的对 方程组消元的三种变换（ 方程组消元的三种变换（统称为线性方程组的初等变 换）： （1）对换方程组中某两个方程的位置； ）对换方程组中某两个方程的位置； 乘以方程组中某个方程； （2）以非零的常数 乘以方程组中某个方程； ）以非零的常数k乘以方程组中某个方程
n→∞
①
生 数 { xn }, 若lim xn = ξ , 则 ξ 即为原方程的根 . 成 列 为 代 为 代 数 ①称为迭代格式 , ϕ(x) 称 迭 函 , x0 称 迭
初值 .
若lim xn 存 称 代 敛, 否则称为发散 . 在 迭 收
n→∞
机动 目录 上页 下页 返回 结束
例1. 用迭代法求方程 x3 − x −1 = 0在[1, 2]内 实 . 的 根 解法1 解法 将方程变形为 x = x3 −1, 迭代格式为
3 xn = xn−1 −1, 取x0 =1.5
n 0 xn 1.5
1 2 3 L 2.375 12.396 1903.779 L
xn = 3 xn−1 +1, 取x0 =1.5
发散 !
解法2 解法 将方程变形为 x = 3 x +1, 迭代格式为
n 0 7 8 1 2 L xn 1.5 1.35721 1.33086 L 1.32472 1.32472
f (xn ) Q f (ξ) = 0, ∴xn −ξ = f ′(η)
f (xn−1) xn = xn−1 − f ′(xn−1)
f (xn ) 记m = min f ′(x) > 0, 则得 xn −ξ ≤ [a,b] m
说明: 说明 用牛顿法时, 若过纵坐标与 f ′′(x)异号的端点作
. 切线 , 则切线与 x 轴焦点的横坐标未必在 [a, b]内
CALL ZANLY (F, ERRABS, ERRREL, NKNOWN, NNEW, & NGUESS, ZINIT,ITMAX, Z, INFO) write(*,"(3('('F5.2','F5.2') '))") Z stop end program complex function F (Z) complex Z F = (Z-5.0)*(Z-4.0)*(Z) return end function
例 4：求复数根
Subroutine ZANLY (F, ERRABS, ERRREL, NKNOWN, NNEW, NGUESS, ZINIT, ITMAX, Z, INFO)
program main use IMSL implicit none integer, parameter :: N = 3 integer :: INFO(N) integer, parameter :: ITMAX = 100 integer, parameter :: NKNOWN = 0 integer, parameter :: NNEW = N, NGUESS = N complex :: ZINIT(N) = (/ (1.0, 0.0), (2.0,0.0), (3.0,0.0) /) real, parameter :: ERRABS = 0.0001 real, parameter :: ERRREL = 0.0001 complex, external :: F complex :: Z(N)
call ZREAL (f, ERRABS, ERRREL, EPS, ETA, & NROOT, ITMAX, XGUESS, X, INFO) write(*,*) X stop end program real function f(x) implicit none real x f=sin(x) return end function
有如下四种情况: y y y a a b o x o bx oa bx f ′ >0 f ′ >0 f ′ <0 f ′′ < 0 f ′′ > 0 f ′′ > 0
机动 目录 上页
f ′ <0 f ′′ < 0
下页 返回 结束
牛顿迭代法的基本思想 牛顿迭代法的基本思想: 用切线近似代替曲线弧求方 基本思想: 程的近似根 . 记纵坐标与 f ′′(x) 同号的端点为
y
a ξ o
x2 x1 x0 x b
(x0 , f (x0 )), 在此点作切线 , 其方程为
y − f (x0 ) = f ′(x0 )(x − x0 ) f (x0 ) 令 y = 0 得它与 x 轴的交点 (x1 , 0),其中 x1 = x0 − f ′(x0 ) 再在点(x1 , f (x1))作切线 , 可得近似根 x2.
乘以方程组中某个方程后加到另一个方程上去。 （3）用数 乘以方程组中某个方程后加到另一个方程上去。 ）用数k乘以方程组中某个方程后加到另一个方程上去
定理1 定理
线性方程组经过初等变换后得到的新方程组 与原方程组同解。 与原方程组同解。
2x1 + 2x2 − x3 = 6 x1 − 2x2 + 4x3 = 3 5x + 7x + x = 28 2 3 1
1 方 x = ϕ(x) 在[a, b] 上 唯 解ξ， ） 程 有 一
2 x0 ∈[a, b] ） n→∞ xn = ϕ(xn−1) ξ →
(证明略)
机动
目录
上页
下页
返回
结束
2. 牛顿迭代法
f (x) 满 : 足 1) 在[a,b] 上 续 f (a) f (b) < 0 连 ,
2) 在[a,b] 上 f ′(x) 及 f ′′(x)不 号 变 方 f (x) = 0 在(a,b)内 唯 的 根ξ . 程 有 一 实
例2. 用牛顿法求方程
误差不超过 0.00001.
3
x3 − 2x2 − 4x − 7 = 0的近似解, 使
2பைடு நூலகம்
解: 设 f (x) = x − 2x − 4x − 7. 由草图可见方程有唯一的正实根 ξ ,且
y
f (3) = −10, f (4) = 9
o
34 x
因 [3， 为 隔 区 . 此 4] 一 根 间
2 第三个方程加上第二个方程的 倍, 得 3
2x1 + 2x2 − x3 = 6 9 −3x2 + x3 = 0 2 13 x3 =13 2 (4.4)
iter = 0 do call fdf( x0, f, df ) if( abs(f)<1.0d-5 ) exit x = x0 - f/df x0 = x iter = iter + 1 if( iter>60 ) then write(*,*) "Wrong x0, no root!" goto 100 end if end do
如此继续下去, 可得求近似根的迭代公式 : f (xn−1) (n =1, 2,L ) xn = xn−1 − f ′(xn−1) 称为牛顿迭代公式 牛顿迭代公式
机动 目录 上页 下页 返回 结束
牛顿法的误差估计: 牛顿法的误差估计: 由微分中值定理得
y
a ξ o
x2 x1 x0 x b
f (xn ) − f (ξ ) = f ′(η)(xn −ξ)
3. 任意函数求根
Subroutine ZREAL(F, ERRABS, ERRREL, EPS, ETA, NROOT, ITMAX, XGUESS, X, INFO) Real, external :: F Real ERRABS Real ERRREL Real EPS Real ETA Integer NROOT Integer ITMAX Real XGUESS(NROOT) Real X(NROOT) Integer INFO(NROOT)
迭代收敛 , 1.32472 为计算精度范围内的所求根 .
机动 目录 上页 下页 返回 结束
迭代法的敛散性与迭代函数的特性有关. 可以证明 下述定理:
区 满 定理. 程 定理. 方 x = ϕ(x) 在 间[a, b]上 足:
1 x = ϕ(x)连 ， a ≤ ϕ(x) ≤ b ） 续且
2 ϕ′(x) 存 ， ϕ′(x) ≤ L <1 ） 在且
迭代法与高斯消元法
计算地球化学
主要内容
1. 一般迭代法 2. 牛顿迭代法 3. 任意函数求根 4. 高斯消元法
1. 一般迭代法
将 程 f (x) = 0转 为 价 程 x = ϕ(x), 方 化 等 方
在隔根区 间 任 一 x0 , 按递推公式 内 取 点
xn =ϕ(xn−1) (n =1, 2,L )
write(*,*) write(*,'(a6,f10.5)') "x =", x write(*,'(a6,i4)') "iter =", iter write(*,*) stop end program newton subroutine fdf( x, f, df ) implicit none double precision :: x, f, df f = x**3 - 2.0*x**2 - 4.0*x - 7.0 df = 3.0*x**2 - 4.0*x - 4.0 return end subroutine fdf
由 在[3， 上 于 4]
′(x) = 3x2 − 4x − 4 = (3x + 2)(x − 2) > 0 f
f ′′(x) = 6x − 4 = 2(3x − 2) > 0
机动
目录
上页
下页
返回
结束
例2的计算程序
program newton implicit none double precision :: x0, x double precision :: f, df integer :: iter 100 write(*,*) "请输入初值 请输入初值X0:" 请输入初值 read(*,*) x0
例1 解线性方程组
下面通过例子来说明高斯消元法的做法。 下面通过例子来说明高斯消元法的做法。
二 方 与 三 方 分 减 第 个 程 解: 第 个 程 第 个 程 别 去 一 方 的 1 5 倍 与 倍 得 2x + 2x − x = 6 , 1 2 3 2 2

9 −3x2 + x3 = 0 2 7 2x2 + x3 =13 2
机动 目录 上页 下页 返回 结束
牛顿迭代法 牛顿迭代法
优点
至少二阶局部收敛，收敛速度较快， 至少二阶局部收敛，收敛速度较快，特别是当迭 代点充分靠近精确解时。 代点充分靠近精确解时。
缺点
对重根收敛速度较慢（线性收敛） 对重根收敛速度较慢（线性收敛） 对初值的选取很敏感， 对初值的选取很敏感，要求初值相当接近真解 在实际计算中， 在实际计算中，可以先用其它方法获得真解的一个粗 糙近似，然后再用牛顿法求解。 糙近似，然后再用牛顿法求解。
例3
program main use IMSL implicit none integer, parameter :: ITMAX = 100 integer, parameter :: NROOT = 2 real, parameter :: EPS = 1.0E-6 real, parameter :: ERRABS = 1.0E-6 real, parameter :: ERRREL = 1.0E-6 real, parameter :: ETA = 1.0 integer INFO(NROOT) real, external :: f real :: X(NROOT) real :: XGUESS(NROOT) = (/ 1.0, 4.0 /)