chapr方程和方程组的迭代数值解法PPT教学课件

1、6.1 6.1 方程求根数值法三、程序框图与程序语言设计subroutine bisecta1=a,b1=bx0=a1+b1/2comp f(a1),f(b1)f(a1)f(b1)0?end subabs(f(x0)E?comp f(a1),f(x0)f(a1)f(x0)0?end subb1=x0a1=x0NYYNNY图6-1 二分法程序框图第1页/共53页6.1 6.1 方程求根数值法10 subroutine bisect(a,b,E,x0)20 a1=a;b1=b30 if(f(a1)*f(b1)E) then60 if (f(a1)*f(x0)E?end subNY10 subrou

2、tine suppo(x0,E,x)20 x=x030 x=g(x)40 i=i+150 if(abs(f(x)E) then60 goto 3070 endif80 end subroutine suppo图6-3 迭代法程序框图第8页/共53页6.1 6.1 方程求根数值法6.1.3 加速迭代法一、算法分析 加速迭代法的基本原理类似Romberg积分中所采用的方法，用误差补偿对所求的近似值进行修正。 设方程f(x)=0，作迭代格式xk+1=g(xk)。又设 为方程的解，根据微分中值定理： 当k很大， 很小；所以在 上， 可视为常数。令 （根据收敛条件 ），则 1( )() ,kkkxgxx

3、 (6-4)kx ,kx( )g( )qg( )1g11()1kkkqxxxq第9页/共53页6.1 6.1 方程求根数值法用误差补偿:改写为迭代过程： 11()1kkkqxxxq11111kkkqxxxqq001011112122223233():(,)():( ,)():(,)xg xxx xxxg xxx xxxg xxx xx(6-5)第10页/共53页6.1 6.1 方程求根数值法二、程序框图与程序语言设计subroutine spsppox2=(x2-q*x1)/(1-q)comp g(x1),dg(x2),f(x2)abs(f(x2)E?end subNYx1=x0 x2=g(x

4、1)q=dg(x2)x1=x2i=i+1图6-4 加速迭代法程序框图第11页/共53页6.1 6.1 方程求根数值法10 subroutine spsppo(x0,E,x2)20 x1=x030 x2=g(x1)40 q=dg(x2)50 x2=(x2-q*x1)/(1-q)60 x1=x270 i=i+180 if(f(x2)E) then90 goto 30100 endif110 end subroutine spsppo第12页/共53页6.1 6.1 方程求根数值法6.1.4 牛顿切线法 一、算法分析牛顿切线法是将复杂的方程f(x)=0化为简单的线性方程来求解。其数学依据如下： 设是

5、方程式(6-1)式的近似根，则在处将f(x)展开为泰勒级数并取前两项得：从而有 作迭代递推计算格式 000( )()()()0f xf xfxxx000()()f xxxfx1()()kkkkf xxxfx(6-6)第13页/共53页6.1 6.1 方程求根数值法则由此迭代递推计算公式可以求得一系列 ，可以证明，在所选的初值 满足条件： 时，迭代公式(6-6)式一定收敛，而且收敛速度很快。为了帮助理解牛顿切线法的原理，下面给出该法的几何解释，如图6-5所示。 kx0 x0000()()1()()f xfxfxfx图6-5 牛顿法几何意义xy0 x1x2x相当于过 作切线与x轴相交即得x1,过

6、作切线与x轴相交即得x1, 依次类推，直到满足精度要求为止。00(,()xf x11( ,( )xf x第14页/共53页6.1 6.1 方程求根数值法二、程序框图与程序语言设计subroutine ntspx2=x1-f/f1comp f(x1),df(x1),f(x2)abs(f(x2)E?end subNYx1=x0f=f(x1)f1=df(x1)x1=x2i=i+1图6-6 牛顿法程序框图第15页/共53页6.1 6.1 方程求根数值法10 subroutine ntsp(x0,E,x2)20 x1=x030 f=f(x1)40 f1=df(x1)50 x2=x1-f/f160 x1=

7、x270 i=i+180 if(f(x2)E) then90 goto 30100 endif110 end subroutine ntsp第16页/共53页6.1 6.1 方程求根数值法6.1.5 牛顿弦割法 一、算法分析弦割法的基本思想是利用差商代替导数来求方程的根。一般弦割法是利用下列差商代替牛顿切线法中的导数 ，即 写成迭代格式这就是用一般弦割法求方程式(6-1)式的根的迭代计算公式，若给定两个初始近似值x0和x1，则反复使用迭代公式（6-7）式进行迭代，可得到根的一系列近似值 。这个系列的极限就是方程式（6-1）式的根, 即 。 ( )fx00()()()kkkf xf xfxxx0

8、10()()()kkkkkxxxxf xf xf x(6-7)kxlimkxx第17页/共53页6.1 6.1 方程求根数值法xy0 x1x2x3x0P1P2PPx0 x1x2x3xy0P1P2PP 图6-7 一般弦割法的几何意义 图6-8 快速弦割法的几何意义 一般弦割法的几何意义如图6-7所示，过两点 ， ，作直线 与x轴相交，则交点的横坐标即为 ，因为直线 的方程为： 00(,()xf x11( ,()xf x0 1PP2x0 1PP第18页/共53页6.1 6.1 方程求根数值法该方程的解为： 上式就是当k=1时的一般弦割法的迭代计算公式。再过 作垂线 交曲线y=f(x)于 点，联接并

9、延长与x轴相交于 点，重复上述过程，可以获得一系列直线 ，此系列的极限位置就是P点，P点的横坐标就是方程的根。但因每次所作的直线都过点 ，收敛速度虽然比二分法快，但比牛顿切线法要慢得多。二、程序框图与程序语言设计 请读者自行设计101010()()()()0f xf xf xxxxx101110()()()xxxxf xf xf x2(,0)x22x P2P3x0102030,nP P P P P PP P0P第19页/共53页6.1 6.1 方程求根数值法6.1.6 快速牛顿弦割法 一、算法分析 快速弦线法是用下面的差商代替牛顿切线法中的导线 进行求根，即：把上式代入牛顿切线法的迭代公式（6

10、-6）式得到快速牛顿弦割法迭代格式为若给定两个初始近似值x0和x1，则反复使用迭代公式（6-8）式进行迭代，可得到根的一系列近似值 。这个系列的极限就是方程式（6-1）式的根, 即 。 ( )fx11()()()kkkkkf xf xfxxx111()()()kkkkkkkxxxxf xf xf x(6-8)kxlimkxx第20页/共53页6.1 6.1 方程求根数值法快速牛顿弦割法的几何意义如图6-8所示，过两点 ， ，作直线 与x轴相交，则交点的横坐标即为 ，因为直线 的方程为：该方程的解为：再过 作垂线 交曲线y=f(x)于 点，联接并延长与x轴相交于 点，描述直线 的方程式为：该方程

11、的解为： 00(,()xf x11( ,()xf x0 1PP2x0 1PP101010()()()()0f xf xf xxxxx101110()()()xxxxf xf xf x2(,0)x22x P2P3x12PP212121()()()()0f xf xf xxxxx212221()()()xxxxf xf xf x第21页/共53页6.1 6.1 方程求根数值法重复上述过程，可以获得一系列直线 ，此系列的极限位置就是P点，P点的横坐标就是方程的根。但因每次所作的直线的终点都是下一次所作直线的起点，即所作的直线系列是 ，而不是 ，所以，其收敛速度比一般弦线法快得多。二、程序框图与程序语

12、言设计0112231,nnP P PP P PP P0112231,nnP P PP P PP P0102030,nP P P P P PP Psubroutine fntspcomp f(x0),f(x1)x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0)abs(x2-x1)E?NYx1=x2;x0=x1end sub第22页/共53页6.1 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) then40 x0=x150 x1=x260

13、goto 2070 endif80 end subroutine fntsp第23页/共53页6.2 6.2 线性方程组求解数值法6.2.1 线性方程组Jacobi迭代法 一、算法分析 为求线性代数方程组(4-1)的解，仿照(6-1)方程求根的办法，可将代数方程组(4-1)改写为等价方程组作构造格式11221331111221 123322221 12211()/()/()/nnnnnnnnnnnnnxa xa xa xbaxa xa xa xbaxa xa xaxba (1)( )( )( )( )11221331111(1)( )( )( )( )221 12332211(1)( )( )

14、( )( )1 12211()/()/()/kkkkknnkkkkknnkkkkknnnnnnnnnxa xa xa xbaxa xa xa xbaxa xa xaxba 第24页/共53页6.2 6.2 线性方程组求解数值法或简写为给定初值 ，并令 ，由此可得向量序列 。显然，如果此序列收敛于x x,那么每个分量序列就必收敛于 ， 就必然是方程组的解。这种方法就是yacobi迭代法。例6-2 用yacobi迭代法求解方程组(1)( )( )( )( )1 1221111( )( )()/,1,kkkkkiiiiiiiiikkinniiixa xa xaxaxa xbain (6-9)(0)(

15、0)(0)(0)T12(,)nxxxx0,1,k ( )kx( )kixix(1, )ix in 1238118172421912xxx 第25页/共53页6.2 6.2 线性方程组求解数值法解 用yacobi迭代格式有取初始值 ，并令 ，得 ，故 ( )( )(1)231( )( )(1)132( )( )(1)123884271229kkkkkkkkkxxxxxxxxx(0)(0)(0)(0)TT123(,)(0,0,0)xxxx0,1,k (1)T(2)T(3)T(11)T(12)(1,0.571249,1.333333)(1.095238,1.095238,1.058201)(0.99

16、5370,1.040812,0.972369)(1.000000,1.000000,1.00000)xxxxx(11)xx第26页/共53页subroutine yacsppif(i.ne.j) then end subi=1,nx(i)=0end dod=0do i=1,ny(i)=b(i)end dodo j=1,ny(i)= y(i) a(i,j)*x(j)endifif(abs(x(i)-y(i)d) then d=abs(x(i)-y(i)y(i)=y(i)/a(i,i)end dox(i)= y(i)i=1,nend doif(dE) then go toendifendif二 程

17、序框图与程序设计第27页/共53页6.2 6.2 线性方程组求解数值法10 subroutine yacspp(a,b,n,E,x)20 dimension a(n,n),b(n),x(n)30 do i=1,n40 x(i)=050 end do60 d=070 do i=1,n 80 y(i)=b(i)90 do j=1,n100 if(i.ne.j) then110 y(i)=y(i)-a(i,j)*x(j) 120 endif130 end do140 y(i)=y(i)/a(i,i)150 if(abs(x(i)-y(i)d) then 160 d=abs(x(i)-y(i)170

18、endif180 end do190 do i=1,n200 x(i)= y(i)210 end do 220 if(dE) then 230 go to 60240 endif250 end subroutine yacspp第28页/共53页6.2 6.2 线性方程组求解数值法6.2.2 线性方程组Gauss-seidel迭代法 一、算法分析 在迭代递推计算通式(6-9)中，第(k+1)次迭代用的只是第k次迭代的近似值。可是，解的各个分量是依次计算的，显然，在计算xi时，它前面的其它未知数的本次迭代近似值已经计算出来了。一般来说，新值总比旧值更接近真值，因为应该优先使用它们。这样改进所得的

19、方法就是高斯-赛德尔迭代法,其迭代格式为 (1)( )( )( )( )11221331111(1)(1)( )( )( )221 12332211(1)(1)(1)(1)( )1 12211()/()/()/kkkkknnkkkkknnkkkkknnnnnnnnnxa xa xa xbaxa xa xa xbaxa xa xaxba 第29页/共53页6.2 6.2 线性方程组求解数值法或简写为例6-3 用Gauss-seidel迭代法求解例6-2方程组解：用(6-10)式，例6-2方程组的Gauss-seidel迭代格式为(1)(1)(1)(1)( )1 1221111( )( )()/,

20、1,kkkkkiiiiiiiiikkinniiixa xa xaxaxa xbain (6-10)( )( )(1)231(1)( )(1)132(1)(1)(1)123884271229kkkkkkkkkxxxxxxxxx第30页/共53页6.2 6.2 线性方程组求解数值法取初始值 ，并令 ，得故 ，解毕。二、程序框图与程序设计(0)(0)(0)(0)TT123(,)(0,0,0)xxxx0,1,k (1)T(2)T(3)T(8)T(9)(1,0.714286,1.031746)(1.039683,1.014740,0.988544)(0.997000,0.996563,1.001049)

21、(1.000000,1.000000,1.00000)xxxxx(9)xx第31页/共53页subroutine gauseiif(i.ne.j) then end subi=1,nx(i)=0end dod=0do i=1,ny=b(i)end dodo j=1,ny= ya(i,j)*x(j)endifif(abs(x(i)-y)d) then d=abs(x(i)-y)y=y/a(i,i)end dox(i)= yif(dE) then go toendifendif第32页/共53页6.2 6.2 线性方程组求解数值法10 subroutine gausei(a,b,n,E,x)20

22、dimension a(n,n),b(n),x(n)30 do i=1,n40 x(i)=050 end do60 d=070 do i=1,n 80 y=b(i)90 do j=1,n100 if(i.ne.j) then110 y=y-a(i,j)*x(j) 120 endif130 end do140 y=y/a(i,i)150 if(abs(x(i)-y)d) then 160 d=abs(x(i)-y)170 endif180 x(i)= y190 end do200 if(dE) then 210 go to 60220 endif230 end subroutine gausei

23、第33页/共53页6.2 6.2 线性方程组求解数值法 为加速Gauss-seidel迭代法的收敛性，仿方程求根的松弛法，将迭代公式(6-10)改为这里 为松弛因子。按此公式迭代求解方程组(4-1)，称为逐个超松弛迭代法或SOR法。显然 时为Gauss-seidel迭代法。(1)(1)(1)(1)( )1 1221111( )( )( )()/(1),1,kkkkkiiiiiiiiikkkinniiiixa xa xaxaxa xbaxin(6-11)1第34页/共53页6.3 6.3 非线性方程组求解数值法 在现代工程技术或科研过程中，常常会遇到非线性代数方程组的求解问题。本节介绍如下方程组

24、的迭代数值解法： 解非线性方程组的方法通常有两大类：一类属于线性化方法，即用一线性代数方程组来近似逼近非线性代数方程组，由此构造一组递推公式，用于逐次逼近所求的根，这类方法有牛顿-拉夫逊方法及其各种改进：另一类方法是把方程组的求解问题转化为求多元函数的极小值的等效问题来解决，这类方法有最速下降法及其各种改进。 11232123123123( ,)0( ,)0( ,)0( ,)0ininiinninf x x xxxfx x xxxf x x xxxfx x xxx(6-12)第35页/共53页6.3 6.3 非线性方程组求解数值法6.3.1 非线性方程组Gauss-yacobi迭代法 一、算法

25、分析 仿照方程求根的简单迭代法，把方程组(6-12) 表示成如下的等价方程：作成迭代格式 1112322123123123( ,)( ,)( ,)( ,)ininiiinnninxg x x xxxxgx x xxxxg x x xxxxgx x xxx(6-13)第36页/共53页6.3 6.3 非线性方程组求解数值法选取一组初始向量 ，并令 ，可以得到一组向量序列 ，如果方程组(6-13)或(6-14)只有唯一解 ，且 收敛，则得逐次收敛于 的近似值。这样求解方程组(6-13)的方法称简单迭代法。 (1)( )( )( )( )( )11123(1)( )( )( )( )( )22123

26、(1)( )( )( )( )( )123(1)( )( )( )123(,)(,)(,)(,kkkkkkinkkkkkkinkkkkkkiiinkkkknnixg xxxxxxgxxxxxxg xxxxxxgxxxx( )( ),)kknx(6-14)(0)(0)(0)(0)(0)(0)T123(,)inxxxxxx0,1,k ( )kxx( )kx( )kx第37页/共53页6.3 6.3 非线性方程组求解数值法例6-4 用简单迭代法解方程组解 作迭代格式取初始值 ，并令 ，得 故 112212140.114/80 xxxexxx( )1( )12( )( )2211(10.1)/4()

27、/8/4kxkkkxxexxx(0)T(0,0)x0,1,k (1)T(2)T(3)T(18)T(19)(0.225,0)(0.2186919321,0.05466796875)(0.2325557961,0.05317841549)(0.2325670051,0.05645151965)xxxxx(18)xx第38页/共53页6.3 6.3 非线性方程组求解数值法一般迭代格式写成向量形式记矩阵可以证明 时迭代收敛。(1)( )()kkgxx(6-15)111122221212( )nnnnnngggxxxgggxxxggggxxxx(6-16)( )1gx第39页/共53页二、程序框图与通用

28、程序设计 subroutine gauyakdo i=1,nread x(i)end dod=0y(i)=g(x(i)do i=1,nend doif(abs(y(i)-x(i)d) then d=abs(y(i)-x(i)endifx(i)= y(i)i=1,nend doif(dE) then go toendifend sub第40页/共53页6.3 6.3 非线性方程组求解数值法10 subroutine gauyac(n,E,x)20 dimension y(n),x(n)30 do i=1,n40 read(*,2x,f6.2) x(i)50 end do60 d=070 y(1)

29、=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,n130 if(abs(y(i)-x(i)d) then 140 d=abs(y(i)-x(i)150 endif160 end do170 do i=1,n180 x(i)= y(i)190 end do 200 if(dE) then 210 go to 60220 endif230 end subroutine gauyac第41页/共53页6.3 6.3 非线性方

30、程组求解数值法6.3.2 非线性方程组Gauss-yacobi-seidel迭代法 一、算法分析 仿照线性代数方程组的Gauss-Seidel迭代方法，可作非线性方程组Gauss-yacobi-seidel迭代法 为二、程序框图与通用程序设计(1)( )( )( )( )( )11123(1)(1)( )( )( )( )22123(1)(1)(2)(1)( )( )121(1)(1)(112(,)(,)(,)(,kkkkkkinkkkkkkinkkkkkkiiiinkkknnxg xxxxxxgxxxxxxg xxxxxxgxx)(1)(1)( )1,)kkkinnxxx(6-17)第42页

31、/共53页subroutine gayasedo i=1,nread x(i)end dod=0 x(i)=g(x(i)do i=1,nend doif(abs(x(i)-y(i)d) then d=abs(x(i)-y(i)endify(i)= x(i)i=1,nend dogo toendifend subif(dE) then 第43页/共53页6.3 6.3 非线性方程组求解数值法10 subroutine gayase(n,E,x)20 dimension y(n),x(n)30 do i=1,n40 read(*,2x,f6.2) x(i)50 end do60 do i=1,n7

32、0 y(i)= x(i)80 end do90 d=0100 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,n160 if(abs(x(i)-y(i)d) then 170 d=abs(x(i)-y(i)180 endif190 end do 200 if(dE) then 210 go to 60220 endif230 end subroutine gauyac第44页/共53页6.3 6.3 非线性

33、方程组求解数值法6.3.3 非线性方程组最速下降迭代法(Gradient Iteration Method) 一、算法分析 1. 由已知方程组(6.12)式构造目标函数于是，方程组(6.12)式的解就上式的零极小值点，反之亦然。 2.计算差商其中 ， (i=1,2,3,n)。式中c为控制常数，一般取c=0.00001。 21231231( ,)( ,)niniiniF x xxxxfx xxxx123123( ,)( ,)iiininiiFF x x xxxxF x x xxxxx00iiiicxxxcx第45页/共53页6.3 6.3 非线性方程组求解数值法可以看出，梯度法实际上就是利用差商代替牛顿法中的偏导数。 3. 具体计算步骤 （1）从给定的不全为零的初值