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

第6章 方程和方程组的迭代法 现代科技领域或工程技术的许多实际问题，常常可以归结为求解函数方程： 上述方程可能是代数方程,也可能是超越方程。 当f(x) 为代数方程（多项式）时，理论上已经证明， 大于五次的多项式一般没有代数解法。 当f(x) 为超越方程时，一般不能用代数方法求其根。 所以，对于一般的方程（6-1），只能用数值方法求解。 本章主要介绍二分法、切线法、弦切法、迭代法。 (6-1) 6.1 方程求根数值法 6.1.1 二分法 二分法是方程求根最常用而且也是最保险的方法之一。 一、基本思想 将区间对分，保留有根的区间，舍去无根的区间。如此往复，以逐步逼近 方程的根。 二、算法分析 第一步：将a, b对分，即a, x0和x0, b，x0=(a+b)/2。若 f(a)f(x0)E) then 60 if (f(a1)*f(x0)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)。又设 为方程的解， 根据微分中值定理： 当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(x1) 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 方程求根数值法 则由此迭代递推计算公式可以求得一系列 ，可以证明，在 所选的初值 满足条件： 时，迭代公式(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 subroutine 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）式 进行迭代，可得到根的一系列近似值 。这个系列的极限 就是方程式（6-1）式的根, 即 。 (6-7) 6.1 方程求根数值法 图6-7 一般弦割法的几何意义 图6-8 快速弦割法的几何意义 一般弦割法的几何意义如图6-7所示，过两点 ， ，作直 线 与x轴相交，则交点的横坐标即为 ，因为直线 的方程为： 6.1 方程求根数值法 该方程的解为： 上式就是当k=1时的一般弦割法的迭代计算公式。再过 作垂线 交曲线y=f(x)于 点，联接并延长与x轴相交于 点，重复上述过程， 可以获得一系列直线 ，此系列的极限位置就是P点 ，P点的横坐标就是方程的根。但因每次所作的直线都过点 ，收敛速度 虽然比二分法快，但比牛顿切线法要慢得多。 二、程序框图与程序语言设计 请读者自行设计 6.1 方程求根数值法 6.1.6 快速牛顿弦割法 一、算法分析 快速弦线法是用下面的差商代替牛顿切线法中的导线 进行求根，即： 把上式代入牛顿切线法的迭代公式（6-6）式得到快速牛顿弦割 法迭代格式为 若给定两个初始近似值x0和x1，则反复使用迭代公式（6-8）式 进行迭代，可得到根的一系列近似值 。这个系列的极限 就是方程式（6-1）式的根, 即 。 (6-8) 6.1 方程求根数值法 快速牛顿弦割法的几何意义如图6-8所示，过两点 ， ，作 直线 与x轴相交，则交点的横坐标即为 ，因为直线 的方程为 ： 该方程的解为： 再过 作垂线 交曲线y=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 subroutine 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,那么每个分量序列 就必收敛于 ， 就必然是方程组的解。这种方法就 是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) 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(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时，它前面的其它未知数的本次迭代近似值已经计算 出来了。一般来说，新值总比旧值更接近真值，因为应该优 先使用它们。这样改进所得的方法就是高斯-赛德尔迭代法, 其迭代格式为 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=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) 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)，称为逐个 超松弛迭代法或SOR法。显然 时为Gauss-seidel迭代法 。 (6-11) 6.3 非线性方程组求解数值法 在现代工程技术或科研过程中，常常会遇到非线性代数方程组的求解问题 。本节介绍如下方程组的迭代数值解法： 解非线性方程组的方法通常有两大类：一类属于线性化方法，即用一线性 代数方程组来近似逼近非线性代数方程组，由此构造一组递推公式，用于 逐次逼近所求的根，这类方法有牛顿-拉夫逊方法及其各种改进：另一类 方法是把方程组的求解问题转化为求多元函数的极小值的等效问题来解决 ，这类方法有最速下降法及其各种改进。 (6-12) 6.3 非线性方程组求解数值法 6.3.1 非线性方程组Gauss-yacobi迭代法 一、算法分析 仿照方程求根的简单迭代法，把方程组(6-12) 表示成如下的 等价方程： 作成迭代格式 (6-13) 6.3 非线性方程组求解数值法 选取一组初始向量 ，并令 ，可 以得到一组向量序列 ，如果方程组(6-13)或(6-14)只有 唯一解 ，且 收敛，则得逐次收敛于 的近似值。这样 求解方程组(6-13)的方法称简单迭代法。 (6-14) 6.3 非线性方程组求解数值法 例6-4 用简单迭代法解方程组 解 作迭代格式 取初始值 ，并令 ，得 故 6.3 非线性方程组求解数值法 一般迭代格式写成向量形式 记矩阵 可以证明 时迭代收敛。 (6-15) (6-16) 二、程序框图与通用程序设计 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 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 非线性方程组求解数值法 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 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(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.3 非线性方程组求解数值法 可以看出，梯度法实际上就是利用差商代替牛顿法中的偏导数。 3. 具体计算步骤 （1）从给定的不全为零的初值 出发 ，设已经计算到第k,得 。 （2）计算目标函数的值 。 （3）如果|F|E) then x(i) =x(i)+cx s=s+df(i)*2 x(i) =x(i)-cx sfl =f0/s call sub compf endif 6.3 非线性方程组求解数值法 10 subroutine gradmt(n,E,x) 20 dimension