当前位置:首页 >> 电力/水利 >>

迭代法与高斯消元法


迭代法与高斯消元法
计算地球化学

主要内容
1. 一般迭代法 2. 牛顿迭代法 3. 任意函数求根 4. 高斯消元法

1. 一般迭代法
将 程 f (x) = 0转 为 价 程 x = ?(x), 方 化 等 方
在隔根区 间 任 一 x0 , 按递推公式 内 取 点

x

n =?(xn?1) (n =1, 2,L )
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
迭代收敛 , 1.32472 为计算精度范围内的所求根 .
机动 目录 上页 下页 返回 结束

迭代法的敛散性与迭代函数的特性有关. 可以证明 下述定理:

区 满 定理. 程 定理. 方 x = ?(x) 在 间[a, b]上 足:
1 x = ?(x)连 , a ≤ ?(x) ≤ b ) 续且
2 ?′(x) 存 , ?′(x) ≤ L <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)内 唯 的 根ξ . 程 有 一 实
有如下四种情况: 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.
如此继续下去, 可得求近似根的迭代公式 : 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 ?ξ)
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]内
机动 目录 上页 下页 返回 结束

牛顿迭代法 牛顿迭代法
优点
至少二阶局部收敛,收敛速度较快, 至少二阶局部收敛,收敛速度较快,特别是当迭 代点充分靠近精确解时。 代点充分靠近精确解时。

缺点
对重根收敛速度较慢(线性收敛) 对重根收敛速度较慢(线性收敛) 对初值的选取很敏感, 对初值的选取很敏感,要求初值相当接近真解 在实际计算中, 在实际计算中,可以先用其它方法获得真解的一个粗 糙近似,然后再用牛顿法求解。 糙近似,然后再用牛顿法求解。

例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] 一 根 间
由 在[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

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

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. 任意函数求根
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)

例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 /)

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

例 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 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. 高斯消元法
矩阵的加、减、乘法 如: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乘以方程组中某个方程

乘以方程组中某个方程后加到另一个方程上去。 (3)用数 乘以方程组中某个方程后加到另一个方程上去。 )用数k乘以方程组中某个方程后加到另一个方程上去

定理1 定理

线性方程组经过初等变换后得到的新方程组 与原方程组同解。 与原方程组同解。
? 2x1 + 2x2 ? x3 = 6 ? ? x1 ? 2x2 + 4x3 = 3 ?5x + 7x + x = 28 2 3 ? 1

例1 解线性方程组

下面通过例子来说明高斯消元法的做法。 下面通过例子来说明高斯消元法的做法。

二 方 与 三 方 分 减 第 个 程 解: 第 个 程 第 个 程 别 去 一 方 的 1 5 ? 倍 与 倍 得 ? 2x + 2x ? x = 6 , 1 2 3 2 2

? ? ? ? ? ? ?

9 ?3x2 + x3 = 0 2 7 2x2 + x3 =13 2

2 第三个方程加上第二个方程的 倍, 得 3 ?
? 2x1 + 2x2 ? x3 = 6 ? 9 ? ? ?3x2 + x3 = 0 2 ? 13 ? x3 =13 ? ? 2 (4.4)

由 理知 程 (4.4)与 方 组 同 方 组 定 1 方 组 原 程 是 解 程 . 由 (4.4)易 x3 = 2, 代 第 个 程 x2 = 3,再 代 知 入 二 方 得 回 第 一 方 得 1 =1.所 原 程 的 为 = (1,3,2)T . 个 程 x 以 方 组 解 x 形 (4.4)的 程 称 阶梯形方程组. 如 方 组 为阶梯形方程组 阶梯形方程组.
从上述的解题过程可以看出,对方程组作消元变换时 实际 从上述的解题过程可以看出 对方程组作消元变换时,实际 对方程组作消元变换时 上是对原方程组施行一系列的初等变换将其化为阶梯形方程 然后通过回代求出原方程组的解。 组,然后通过回代求出原方程组的解。同时可以发现 对方程组 然后通过回代求出原方程组的解 同时可以发现,对方程组 作初等变换消元时,只是对未知量的系数和常数项进行运算 只是对未知量的系数和常数项进行运算。 作初等变换消元时 只是对未知量的系数和常数项进行运算。

一般n*n线性方程组 一般n*n线性方程组
?a11x1 + a12 x2 +L+ a1n xn = b1 ?a x + a x +L+ a x = b ? 21 1 22 2 2n n 2 ? L L ? L ?an1x1 + an2 x2 +L+ ann xn = bn ?
AX = b

?a11 a12 L ? ?a21 a22 L A= ? L L ? ?a ? n1 an2 L

a1n ? ? a2n ? L ? ? ann ? ?

? x1 ? ? ? ? x2 ? X =? ? M ? ? ?x ? ? n?

?b1 ? ? ? ?b2 ? b=? ? M ? ? ?b ? ? n?

三角形方程组
=b ?a11x1 1 ?a x + a x = b2 ? 21 1 22 2 ? ?L L L ?an1x1 + an2 x2 +L+ annxn = bn ?

?a11x1 + a12 x2 +L+ a1n xn = b1 ? a22 x2 +L+ a2n xn = b2 ? ? L ? ? ann xn = bn ?

(高斯消元法)

消元公式

( ? (k ) akjk?1) j = k, k +1,L, n +1 ?akj = (k?1) akk ? ? (k ) ( ( ( aij = aijk?1) ? aikk?1) ? akjk ) ? ? ? j = k +1,L, n +1 ?i = k +1,Ln ?

回代公式

( ) ?xn = annn+1 , ? n ? ( ) ( xk = akkn+1 ? ∑akjk ) xj k = n ?1 K,1 , , ? j =k +1 ?

高斯消元法的子程序
! Subroutine to solve a linear equation group Subroutine AGAUS(A,B,N,X) DIMENSION A(N,N), X(N), B(N), JS(N) DOUBLE PRECISION A, B, X, T, D L=1 DO 50 K=1,N-1 D = 0.0 DO 210 I = K,N DO 210 J = K,N IF(dABS(A(I,J)).GT.D)THEN D=dABS(A(I,J)) JS(K)=J IS = I ENDIF

210 CONTINUE IF(D<=1d-80)THEN L =0 ELSE IF(JS(K).NE.K)THEN DO 220 I=1,N T=A(I,K) A(I,K)=A(I,JS(K)) A(I,JS(K))=T 220 CONTINUE ENDIF IF(IS.NE.K)THEN DO 230 J=K,N T=A(K,J) A(K,J)=A(IS,J) A(IS,J)=T

230 CONTINUE T=B(K) B(K)=B(IS) B(IS)=T ENDIF ENDIF IF(L.EQ.0)THEN WRITE(*,100) RETURN ENDIF DO 10 J=K+1,N A(K,J)=A(K,J)/A(K,K) 10 CONTINUE B(K)=B(K)/A(K,K) DO 30 I = K+1, N DO 20 J=K+1,N A(I,J)=A(I,J)-A(I,K)*A(K,J)

20

CONTINUE B(I)=B(I)-A(I,K)*B(K) 30 CONTINUE 50 CONTINUE IF(dABS(A(N,N))<=1d-80)THEN L=0 WRITE(*,100) RETURN ENDIF X(N)=B(N)/A(N,N) DO 70 I=N-1,1, -1 T= 0.0 DO 60 J=I+1,N T= T+A(I,J)*X(J) 60 CONTINUE

X(I)=B(I)-T 70 CONTINUE 100 FORMAT(1X,'FAIL') JS(N)=N DO 150 K=N,1,-1 IF(JS(K).NE. K)THEN T=X(K) X(K)=X(JS(K)) X(JS(K))=T ENDIF 150 CONTINUE RETURN END Subroutine AGAUS

例 5:高斯消元法
program main Implicit none double precision :: A(3,3) = (/ 1,2,3,& 1,5,6,& 2,8,9 /) double precision :: B(3) = (/ 4, 15, 18 /), C(3) ! 求解 x+y+2z=4 ! 2x+5y+8z=15 ! 3x+6y+9z=18 Call AGAUS(A,B,3,C) write(*,*) C stop end program

解法2 解法2:逆矩阵法
program main use IMSL real :: A(3,3) = (/ 1,2,3,& 1,5,6,& 2,8,9 /) real :: B(3) = (/ 4, 15, 18 /), C(3) ! 求解 x+y+2z=4 ! 2x+5y+8z=15 ! 3x+6y+9z=18 C = A .ix. B ! invert(A) * B write(*,*) C stop end program

解法3 lin_sol_gen法 解法3:lin_sol_gen法
program main use IMSL implicit none real :: A(3,3) = (/ 1,3,2,& 1,2,1,& 2,1,3 /) real :: B(3,1) = (/ 4,6,6 /) real :: X(3,1) call lin_sol_gen(A,B,X) ! A*X=B,解X 解 write(*,"(3F5.2)") X stop end program

注意: 注意:接受的参数必 须是二维数组,故B和 须是二维数组, 和 X声明成二维数组。 声明成二维数组。 声明成二维数组

作业
1. 用牛顿迭代法编程计算方程 sin(x)+cos(x) = 0的 0的 一个近似实数根,误差不超过0.00001。 一个近似实数根,误差不超过0.00001。 2. 用高斯消元法编程解方程组: x+y+2z+w=4 2x+5y+8z-w=15 3x+6y+9z+2w=18 4x+5y+6z+7w = 30


相关文章:
线性方程组的直接法和迭代法
2 。这样,我们利用高斯消元法 解决了这个方程组。 2. 线性方程组迭代法迭代法就是用某种极限过程去逐步逼近线性方程组精确解的方法.该方法 具有对计算机的存贮...
高斯消去法高斯塞德尔迭代法
数值计算高斯消去法和高斯-塞德尔迭代法 摘要虽然已学过加减消元法、代入消元法、矩阵变换法和 Cramer 法则等,但是 无法满足实际计算需要, 故在此讨论在计算机上...
迭代法实验
实验目的 (1)深入理解线性方程组的迭代法的设计思想,学会利用系数矩阵的性质以...平方根法的运算量和存贮量约为高斯消元法的二分之一,因此它是求解对称正定 ...
数值分析-高斯消元法
数值分析-牛顿迭代法实验报告 一、 实验内容和要求用列主元高斯消去法解线性方程组 Ax=b 3.01 6.03 1.99 方程 1: 1.27 4.16 ?1.23 0.987 ?4.81 9.34...
数值分析编程及运行结果(高斯顺序消元法)
数值分析编程及运行结果(高斯顺序消元法)_数学_自然科学_专业资料。高斯顺序消元法;高斯选列主元消元法;追赶法;高斯-塞德尔迭代格式;雅戈比迭代格式;逐次超松弛法...
LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel...
一、实验目的及题目 实验目的及题目 1.1 实验目的: 验目的: (1)学会用高斯列主元消去法,LU 分解法,Jacobi 迭代法和 Gauss-Seidel 迭代法解线性 方程组。 ...
数值分析试题及答案
。 选取初值 1、 若 A 是 n 阶非奇异矩阵, 则线性方程组 AX=b 一定可以使用高斯消元法求解。 ( × ) 2、 解非线性方程 f(x)=0 的牛顿迭代法在单根...
牛顿迭代法高斯消元法拉格朗日插值公式 复化Simpson 二...
牛顿迭代法高斯消元法拉格朗日插值公式 复化Simpson 二分法 数值分析 牛顿迭代法高斯消元法拉格朗日插值公式 复化Simpson 二分法数值分析 牛顿迭代法高斯消元法拉格朗日...
计算方法矩阵直接法与迭代法
高斯消元法,结果为 [0,…0] 即出现计算错误,弹出,系数矩阵奇异, jacobi B 的谱半径为 86.3375,Jacobi 迭代法 发散,结果为[NaN,…,NaN] Gauss-S 方法 ...
高斯迭代法解线性方程组
高斯迭代法解线性方程组_工学_高等教育_教育专区。高斯迭代法解线性方程组 一. 问题描述迭代法求解线性方程组 用 Gauss-Seidel 迭代法求解线性方程组 由 Jacobi ...
更多相关标签:
高斯消元法和迭代法 | 高斯赛德尔迭代法 | 高斯牛顿迭代法 | 高斯迭代法 | 高斯赛德尔迭代法原理 | 高斯牛顿迭代法步骤 | 高斯赛德尔迭代法例题 | 高斯牛顿迭代法matlab |