chapr-6 方程和方程组的迭代数值解法.ppt

1、第6章 方程和方程组的迭代法,现代科技领域或工程技术的许多实际问题，常常可以归结为求解函数方程： 上述方程可能是代数方程,也可能是超越方程。 当f(x) 为代数方程（多项式）时，理论上已经证明， 大于五次的多项式一般没有代数解法。 当f(x) 为超越方程时，一般不能用代数方法求其根。 所以，对于一般的方程（6-1），只能用数值方法求解。 本章主要介绍二分法、切线法、弦切法、迭代法。,(6-1),6.1 方程求根数值法,6.1.1 二分法 二分法是方程求根最常用而且也是最保险的方法之一。 一、基本思想 将区间对分，保留有根的区间，舍去无根的区间。如此往复，以逐步逼近方程的根。 二、算法分析 第一

2、步：将a, b对分，即a, x0和x0, b，x0=(a+b)/2。若 f(a)f(x0)0，则根x*(a, x0)，并令 a1=a, b1=x0；若 f(x0)f(b)0， 则根x*( x0,b)，并令 a1 =x0, b1 =b 。 第二步：将a1, b1再对分，即a1, x0和x0, b1，x1=(a1 + b1)/2。若 f(a1)f(x1)0，则根x*(a1, x1)，并令 a2 = a1, b2 =x1；若 f(x0)f(b1)0，则根x*( x1, b1)并令a2 =x1, b1= b1 。直到满足精度要求为止，这样便得到一系列的对分结果系列：,6.1 方程求根数值法,三、程序框

3、图与程序语言设计,subroutine bisect,a1=a,b1=b,x0=a1+b1/2,comp f(a1),f(b1),f(a1)f(b1)0?,end sub,abs(f(x0)E?,comp f(a1),f(x0),f(a1)f(x0)0?,end sub,b1=x0,a1=x0,N,Y,Y,N,N,Y,图6-1 二分法程序框图,6.1 方程求根数值法,10 subroutine bisect(a,b,E,x0) 20 a1=a;b1=b 30 if(f(a1)*f(b1)E) then 60 if (f(a1)*f(x0)0) then 70 b1=x0 80 else 90

4、a1=x0 100 endif 110 goto 40 120 endif 130 endif 140 end subroutine bisect,6.1 方程求根数值法,6.1.2 简单迭代法 一、算法分析 将方程(6-1)改写成等价形式 ，则自然希望使用相同的公式 若给定根的一个初始近似值 ，则由迭代式（6-2）式可以计算得一系列的 值： 。利用（6-2）式求根近似值的方法称为简单迭代法。 称为迭代序列， 称为迭代函数，上式称为迭代格式。显然，如果极限 存在，那么必定有 成立，即就是方程（6-1）式的根。,(6-2),6.1 方程求根数值法,例6-1 用简单迭代法求方程 在1附近的根。 解

5、：将原方程变形为如下的等价形式： 作迭代格式 取 ，迭代得,6.1 方程求根数值法,可见经过四次迭代，方程的根的近似值已经精确到小数点后第六位，第五次以后的迭代值则不在改变。所以原方程的近似根为 。 对于本题，若作迭代格式 取 ，迭代得 随着k的增大 也增大,而且不不趋向于任何极限值，这时迭代过程是发散的。,6.1 方程求根数值法,从以上2种解法来看，迭代序列是否收敛或收敛快慢的问题与选取迭代函数g(x)有关。下面研究如何相关。 先看迭代法的几何意义。,图6-2 迭代法几何意义,在直解坐标系中同时作y=x 和y=g(x)两条曲线，如图6-2所示，则这两条曲线的交点的横坐标就是方程式(6-2)式

6、的根 。也就是方程式(6.1)式的根。迭代法(6-2)由 求 ，相当于过曲线上 作水平线与直线y=x相交，过,6.1 方程求根数值法,交点作x轴的垂线，此时垂足至原点距离等于 ，故垂足横坐标为 。由图可见，曲线斜率 时迭代序列收敛，且 越小收敛越快；反之，若 ，则迭代序列发散。于是有如下条件。 设迭代函数g(x)为定义在区间a,b上的连续函数，且对任意属于该区间内的x都有g(x)也属于该区间，即 成立 也成立，那么在区间a,b内必有下式成立： 式(6-3)就是迭代函数g(x)收敛的充分必要条件。,(6-3),6.1 方程求根数值法,二、程序框图与程序语言设计,subroutine suppo,

7、i=i+1,comp g(x),f(x),abs(f(x)E?,end sub,N,Y,10 subroutine suppo(x0,E,x) 20 x=x0 30 x=g(x) 40 i=i+1 50 if(abs(f(x)E) then 60 goto 30 70 endif 80 end subroutine suppo,图6-3 迭代法程序框图,6.1 方程求根数值法,6.1.3 加速迭代法 一、算法分析 加速迭代法的基本原理类似Romberg积分中所采用的方法，用误差补偿对所求的近似值进行修正。 设方程f(x)=0，作迭代格式xk+1=g(xk)。又设 为方程的解，根据微分中值定理：

8、 当k很大， 很小；所以在 上， 可视为常数。令 （根据收敛条件 ），则,(6-4),6.1 方程求根数值法,用误差补偿: 改写为 迭代过程：,(6-5),6.1 方程求根数值法,二、程序框图与程序语言设计,subroutine spsppo,x2=(x2-q*x1)/(1-q),comp g(x1),dg(x2),f(x2),abs(f(x2)E?,end sub,N,Y,x1=x0,x2=g(x1),q=dg(x2),x1=x2,i=i+1,图6-4 加速迭代法程序框图,6.1 方程求根数值法,10 subroutine spsppo(x0,E,x2) 20 x1=x0 30 x2=g(x

9、1) 40 q=dg(x2) 50 x2=(x2-q*x1)/(1-q) 60 x1=x2 70 i=i+1 80 if(f(x2)E) then 90 goto 30 100 endif 110 end subroutine spsppo,6.1 方程求根数值法,6.1.4 牛顿切线法 一、算法分析 牛顿切线法是将复杂的方程f(x)=0化为简单的线性方程来求解。其数学依据如下： 设是方程式(6-1)式的近似根，则在处将f(x)展开为泰勒级数并取前两项得： 从而有 作迭代递推计算格式,(6-6),6.1 方程求根数值法,则由此迭代递推计算公式可以求得一系列 ，可以证明，在所选的初值 满足条件：

10、 时，迭代公式(6-6)式一定收敛，而且收敛速度很快。为了帮助理解牛顿切线法的原理，下面给出该法的几何解释，如图6-5所示。,图6-5 牛顿法几何意义,相当于过 作切线与x轴相交即得x1,过 作切线与x轴相交即得x1, 依次类推，直到满足精度要求为止。,6.1 方程求根数值法,二、程序框图与程序语言设计,subroutine ntsp,x2=x1-f/f1,comp f(x1),df(x1),f(x2),abs(f(x2)E?,end sub,N,Y,x1=x0,f=f(x1),f1=df(x1),x1=x2,i=i+1,图6-6 牛顿法程序框图,6.1 方程求根数值法,10 subrouti

11、ne ntsp(x0,E,x2) 20 x1=x0 30 f=f(x1) 40 f1=df(x1) 50 x2=x1-f/f1 60 x1=x2 70 i=i+1 80 if(f(x2)E) then 90 goto 30 100 endif 110 end subroutine ntsp,6.1 方程求根数值法,6.1.5 牛顿弦割法 一、算法分析 弦割法的基本思想是利用差商代替导数来求方程的根。一般弦割法是利用下列差商代替牛顿切线法中的导数 ，即 写成迭代格式 这就是用一般弦割法求方程式(6-1)式的根的迭代计算公式，若给定两个初始近似值x0和x1，则反复使用迭代公式（6-7）式进行迭代，

12、可得到根的一系列近似值 。这个系列的极限就是方程式（6-1）式的根, 即 。,(6-7),6.1 方程求根数值法,图6-7 一般弦割法的几何意义,图6-8 快速弦割法的几何意义,一般弦割法的几何意义如图6-7所示，过两点 ， ，作直线 与x轴相交，则交点的横坐标即为 ，因为直线 的方程为：,6.1 方程求根数值法,该方程的解为： 上式就是当k=1时的一般弦割法的迭代计算公式。再过 作垂线 交曲线y=f(x)于 点，联接并延长与x轴相交于 点，重复上述过程，可以获得一系列直线 ，此系列的极限位置就是P点，P点的横坐标就是方程的根。但因每次所作的直线都过点 ，收敛速度虽然比二分法快，但比牛顿切线法

13、要慢得多。 二、程序框图与程序语言设计 请读者自行设计,6.1 方程求根数值法,6.1.6 快速牛顿弦割法 一、算法分析 快速弦线法是用下面的差商代替牛顿切线法中的导线 进行求根，即： 把上式代入牛顿切线法的迭代公式（6-6）式得到快速牛顿弦割法迭代格式为 若给定两个初始近似值x0和x1，则反复使用迭代公式（6-8）式进行迭代，可得到根的一系列近似值 。这个系列的极限就是方程式（6-1）式的根, 即 。,(6-8),6.1 方程求根数值法,快速牛顿弦割法的几何意义如图6-8所示，过两点 ， ，作直线 与x轴相交，则交点的横坐标即为 ，因为直线 的方程为： 该方程的解为： 再过 作垂线 交曲线y

14、=f(x)于 点，联接并延长与x轴相交于 点，描述直线 的方程式为： 该方程的解为：,6.1 方程求根数值法,重复上述过程，可以获得一系列直线 ，此系列的极限位置就是P点，P点的横坐标就是方程的根。但因每次所作的直线的终点都是下一次所作直线的起点，即所作的直线系列是 ，而不是 ，所以，其收敛速度比一般弦线法快得多。 二、程序框图与程序语言设计,subroutine fntsp,comp f(x0),f(x1),x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0),abs(x2-x1)E?,N,Y,x1=x2;x0=x1,end sub,6.1 方程求根数值法,10 subrouti

15、ne fntsp(x0,x1,E,x2) 20 x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0) 30 if(abs(x2-x1)E) then 40 x0=x1 50 x1=x2 60 goto 20 70 endif 80 end subroutine fntsp,6.2 线性方程组求解数值法,6.2.1 线性方程组Jacobi迭代法 一、算法分析 为求线性代数方程组(4-1)的解，仿照(6-1)方程求根的办法，可将代数方程组(4-1)改写为等价方程组 作构造格式,6.2 线性方程组求解数值法,或简写为 给定初值 ，并令 ，由此可得向量序列 。显然，如果此序列收敛于x,那么

16、每个分量序列 就必收敛于 ， 就必然是方程组的解。这种方法就是yacobi迭代法。 例6-2 用yacobi迭代法求解方程组,(6-9),6.2 线性方程组求解数值法,解 用yacobi迭代格式有 取初始值 ，并令 ，得 ，故,subroutine yacspp,if(i.ne.j) then,end sub,i=1,n,x(i)=0,end do,d=0,do i=1,n,y(i)=b(i),end do,do j=1,n,y(i)= y(i) a(i,j)*x(j),endif,if(abs(x(i)-y(i)d) then,d=abs(x(i)-y(i),y(i)=y(i)/a(i,i)

17、,end do,x(i)= y(i),i=1,n,end do,if(dE) then,go to,endif,endif,二 程序框图与程序设计,6.2 线性方程组求解数值法,10 subroutine yacspp(a,b,n,E,x) 20 dimension a(n,n),b(n),x(n) 30 do i=1,n 40 x(i)=0 50 end do 60 d=0 70 do i=1,n 80 y(i)=b(i) 90 do j=1,n 100 if(i.ne.j) then 110 y(i)=y(i)-a(i,j)*x(j) 120 endif 130 end do 140 y(

18、i)=y(i)/a(i,i),150 if(abs(x(i)-y(i)d) then 160 d=abs(x(i)-y(i) 170 endif 180 end do 190 do i=1,n 200 x(i)= y(i) 210 end do 220 if(dE) then 230 go to 60 240 endif 250 end subroutine yacspp,6.2 线性方程组求解数值法,6.2.2 线性方程组Gauss-seidel迭代法 一、算法分析 在迭代递推计算通式(6-9)中，第(k+1)次迭代用的只是第k次迭代的近似值。可是，解的各个分量是依次计算的，显然，在计算xi

19、时，它前面的其它未知数的本次迭代近似值已经计算出来了。一般来说，新值总比旧值更接近真值，因为应该优先使用它们。这样改进所得的方法就是高斯-赛德尔迭代法,其迭代格式为,6.2 线性方程组求解数值法,或简写为 例6-3 用Gauss-seidel迭代法求解例6-2方程组 解：用(6-10)式，例6-2方程组的Gauss-seidel迭代格式为,(6-10),6.2 线性方程组求解数值法,取初始值 ，并令 ，得 故 ，解毕。 二、程序框图与程序设计,subroutine gausei,if(i.ne.j) then,end sub,i=1,n,x(i)=0,end do,d=0,do i=1,n,y

20、=b(i),end do,do j=1,n,y= ya(i,j)*x(j),endif,if(abs(x(i)-y)d) then,d=abs(x(i)-y),y=y/a(i,i),end do,x(i)= y,if(dE) then,go to,endif,endif,6.2 线性方程组求解数值法,10 subroutine gausei(a,b,n,E,x) 20 dimension a(n,n),b(n),x(n) 30 do i=1,n 40 x(i)=0 50 end do 60 d=0 70 do i=1,n 80 y=b(i) 90 do j=1,n 100 if(i.ne.j)

21、 then 110 y=y-a(i,j)*x(j) 120 endif 130 end do 140 y=y/a(i,i),150 if(abs(x(i)-y)d) then 160 d=abs(x(i)-y) 170 endif 180 x(i)= y 190 end do 200 if(dE) then 210 go to 60 220 endif 230 end subroutine gausei,6.2 线性方程组求解数值法,为加速Gauss-seidel迭代法的收敛性，仿方程求根的松弛法，将迭代公式(6-10)改为 这里 为松弛因子。按此公式迭代求解方程组(4-1)，称为逐个超松弛迭

22、代法或SOR法。显然 时为Gauss-seidel迭代法。,(6-11),6.3 非线性方程组求解数值法,在现代工程技术或科研过程中，常常会遇到非线性代数方程组的求解问题。本节介绍如下方程组的迭代数值解法： 解非线性方程组的方法通常有两大类：一类属于线性化方法，即用一线性代数方程组来近似逼近非线性代数方程组，由此构造一组递推公式，用于逐次逼近所求的根，这类方法有牛顿-拉夫逊方法及其各种改进：另一类方法是把方程组的求解问题转化为求多元函数的极小值的等效问题来解决，这类方法有最速下降法及其各种改进。,(6-12),6.3 非线性方程组求解数值法,6.3.1 非线性方程组Gauss-yacobi迭代

23、法 一、算法分析 仿照方程求根的简单迭代法，把方程组(6-12) 表示成如下的等价方程： 作成迭代格式,(6-13),6.3 非线性方程组求解数值法,选取一组初始向量 ，并令 ，可以得到一组向量序列 ，如果方程组(6-13)或(6-14)只有唯一解 ，且 收敛，则得逐次收敛于 的近似值。这样求解方程组(6-13)的方法称简单迭代法。,(6-14),6.3 非线性方程组求解数值法,例6-4 用简单迭代法解方程组 解 作迭代格式 取初始值 ，并令 ，得 故,6.3 非线性方程组求解数值法,一般迭代格式写成向量形式 记矩阵 可以证明 时迭代收敛。,(6-15),(6-16),二、程序框图与通用程序设

24、计,subroutine gauyak,do i=1,n,read x(i),end do,d=0,y(i)=g(x(i),do i=1,n,end do,if(abs(y(i)-x(i)d) then,d=abs(y(i)-x(i),endif,x(i)= y(i),i=1,n,end do,if(dE) then,go to,endif,end sub,6.3 非线性方程组求解数值法,10 subroutine gauyac(n,E,x) 20 dimension y(n),x(n) 30 do i=1,n 40 read(*,2x,f6.2) x(i) 50 end do 60 d=0

25、70 y(1)=g1(x(1),x(2),x(n) 80 y(2)=g2(x(1),x(2),x(n) 90 y(3)=g3(x(1),x(2),x(n) 100 110 y(n)=gn(x(1),x(2),x(n) 120 do i=1,n 130 if(abs(y(i)-x(i)d) then 140 d=abs(y(i)-x(i),150 endif 160 end do 170 do i=1,n 180 x(i)= y(i) 190 end do 200 if(dE) then 210 go to 60 220 endif 230 end subroutine gauyac,6.3 非

26、线性方程组求解数值法,6.3.2 非线性方程组Gauss-yacobi-seidel迭代法 一、算法分析 仿照线性代数方程组的Gauss-Seidel迭代方法，可作非线性方程组Gauss-yacobi-seidel迭代法 为 二、程序框图与通用程序设计,(6-17),subroutine gayase,do i=1,n,read x(i),end do,d=0,x(i)=g(x(i),do i=1,n,end do,if(abs(x(i)-y(i)d) then,d=abs(x(i)-y(i),endif,y(i)= x(i),i=1,n,end do,go to,endif,end sub,

27、if(dE) then,6.3 非线性方程组求解数值法,10 subroutine gayase(n,E,x) 20 dimension y(n),x(n) 30 do i=1,n 40 read(*,2x,f6.2) x(i) 50 end do 60 do i=1,n 70 y(i)= x(i) 80 end do 90 d=0 100 x(1)=g3(x(1),x(2),x(n) 110 x(2)=g3(x(1),x(2),x(n) 120 x(3)=g3(x(1),x(2),x(n) 130 140 x(n)=gn(x(1),x(2),x(n),150 do i=1,n 160 if(

28、abs(x(i)-y(i)d) then 170 d=abs(x(i)-y(i) 180 endif 190 end do 200 if(dE) then 210 go to 60 220 endif 230 end subroutine gauyac,6.3 非线性方程组求解数值法,6.3.3 非线性方程组最速下降迭代法(Gradient Iteration Method) 一、算法分析 1. 由已知方程组(6.12)式构造目标函数 于是，方程组(6.12)式的解就上式的零极小值点，反之亦然。 2.计算差商 其中 ， (i=1,2,3,n)。 式中c为控制常数，一般取c=0.00001。,6

29、.3 非线性方程组求解数值法,可以看出，梯度法实际上就是利用差商代替牛顿法中的偏导数。 3. 具体计算步骤 （1）从给定的不全为零的初值 出发，设已经计算到第k,得 。 （2）计算目标函数的值 。 （3）如果|F|E（容许误差），则认为 就是所求的一组解，否则继续做下一步。 （4）计算差商 其中 ， (i=1,2,3,n)。 式中c为控制常数，一般取c=0.00001。,6.3 非线性方程组求解数值法,（4）计算 其中 然后再从第二步开始。以上各步中的上标k表示计算次数。 二、通用程序框图与程序设计,subroutine gradmt,do i=1,n,read x(i),end do,c=0.00001;s=0,cx=c*x(i),do i=1,n,end do,if(abs(x(i)=0) then,cx=c,else,x(i)= x(i)-sfl*df(i),i=1,n,end do,endif,end sub,d