chapr方程和方程组的迭代数值解法.pptx

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 subroutine 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很大，很小；所以在 上，可视为常数。令 （根据收敛条件 ），则 (6-4)第9页/共53页6.1 6.1 方程求根数值法用误差补偿:改写为迭代过程：(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(x1)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化为简单的线性方程来求解。其数学依据如下：设是方程式(6-1)式的近似根，则在处将f(x)展开为泰勒级数并取前两项得：从而有 作迭代递推计算格式(6-6)第13页/共53页6.1 6.1 方程求根数值法则由此迭代递推计算公式可以求得一系列 ，可以证明，在所选的初值 满足条件：时，迭代公式(6-6)式一定收敛，而且收敛速度很快。为了帮助理解牛顿切线法的原理，下面给出该法的几何解释，如图6-5所示。图6-5 牛顿法几何意义相当于过 作切线与x轴相交即得x1,过 作切线与x轴相交即得x1,依次类推，直到满足精度要求为止。第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=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）式的根,即 。(6-7)第17页/共53页6.1 6.1 方程求根数值法 图6-7 一般弦割法的几何意义 图6-8 快速弦割法的几何意义 一般弦割法的几何意义如图6-7所示，过两点 ，作直线 与x轴相交，则交点的横坐标即为 ，因为直线 的方程为：第18页/共53页6.1 6.1 方程求根数值法该方程的解为：上式就是当k=1时的一般弦割法的迭代计算公式。再过 作垂线 交曲线y=f(x)于 点，联接并延长与x轴相交于 点，重复上述过程，可以获得一系列直线 ，此系列的极限位置就是P点，P点的横坐标就是方程的根。但因每次所作的直线都过点 ，收敛速度虽然比二分法快，但比牛顿切线法要慢得多。二、程序框图与程序语言设计 请读者自行设计第19页/共53页6.1 6.1 方程求根数值法6.1.6 快速牛顿弦割法 一、算法分析 快速弦线法是用下面的差商代替牛顿切线法中的导线 进行求根，即：把上式代入牛顿切线法的迭代公式（6-6）式得到快速牛顿弦割法迭代格式为若给定两个初始近似值x0和x1，则反复使用迭代公式（6-8）式进行迭代，可得到根的一系列近似值 。这个系列的极限就是方程式（6-1）式的根,即 。(6-8)第20页/共53页6.1 6.1 方程求根数值法快速牛顿弦割法的几何意义如图6-8所示，过两点 ，作直线 与x轴相交，则交点的横坐标即为 ，因为直线 的方程为：该方程的解为：再过 作垂线 交曲线y=f(x)于 点，联接并延长与x轴相交于 点，描述直线 的方程式为：该方程的解为：第21页/共53页6.1 6.1 方程求根数值法重复上述过程，可以获得一系列直线 ，此系列的极限位置就是P点，P点的横坐标就是方程的根。但因每次所作的直线的终点都是下一次所作直线的起点，即所作的直线系列是 ，而不是 ，所以，其收敛速度比一般弦线法快得多。二、程序框图与程序语言设计subroutine 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 goto 2070 endif80 end subroutine fntsp第23页/共53页6.2 6.2 线性方程组求解数值法6.2.1 线性方程组Jacobi迭代法 一、算法分析 为求线性代数方程组(4-1)的解，仿照(6-1)方程求根的办法，可将代数方程组(4-1)改写为等价方程组作构造格式第24页/共53页6.2 6.2 线性方程组求解数值法或简写为给定初值 ，并令 ，由此可得向量序列 。显然，如果此序列收敛于x x,那么每个分量序列就必收敛于 ，就必然是方程组的解。这种方法就是yacobi迭代法。例6-2 用yacobi迭代法求解方程组(6-9)第25页/共53页6.2 6.2 线性方程组求解数值法解 用yacobi迭代格式有取初始值 ，并令 ，得 ，故 第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二 程序框图与程序设计第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 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时，它前面的其它未知数的本次迭代近似值已经计算出来了。一般来说，新值总比旧值更接近真值，因为应该优先使用它们。这样改进所得的方法就是高斯-赛德尔迭代法,其迭代格式为 第29页/共53页6.2 6.2 线性方程组求解数值法或简写为例6-3 用Gauss-seidel迭代法求解例6-2方程组解：用(6-10)式，例6-2方程组的Gauss-seidel迭代格式为(6-10)第30页/共53页6.2 6.2 线性方程组求解数值法取初始值 ，并令 ，得故 ，解毕。二、程序框图与程序设计第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 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第33页/共53页6.2 6.2 线性方程组求解数值法 为加速Gauss-seidel迭代法的收敛性，仿方程求根的松弛法，将迭代公式(6-10)改为这里 为松弛因子。按此公式迭代求解方程组(4-1)，称为逐个超松弛迭代法或SOR法。显然 时为Gauss-seidel迭代法。(6-11)第34页/共53页6.3 6.3 非线性方程组求解数值法 在现代工程技术或科研过程中，常常会遇到非线性代数方程组的求解问题。本节介绍如下方程组的迭代数值解法：解非线性方程组的方法通常有两大类：一类属于线性化方法，即用一线性代数方程组来近似逼近非线性代数方程组，由此构造一组递推公式，用于逐次逼近所求的根，这类方法有牛顿-拉夫逊方法及其各种改进：另一类方法是把方程组的求解问题转化为求多元函数的极小值的等效问题来解决，这类方法有最速下降法及其各种改进。(6-12)第35页/共53页6.3 6.3 非线性方程组求解数值法6.3.1 非线性方程组Gauss-yacobi迭代法 一、算法分析 仿照方程求根的简单迭代法，把方程组(6-12)表示成如下的等价方程：作成迭代格式(6-13)第36页/共53页6.3 6.3 非线性方程组求解数值法选取一组初始向量 ，并令 ，可以得到一组向量序列 ，如果方程组(6-13)或(6-14)只有唯一解 ，且 收敛，则得逐次收敛于 的近似值。这样求解方程组(6-13)的方法称简单迭代法。(6-14)第37页/共53页6.3 6.3 非线性方程组求解数值法例6-4 用简单迭代法解方程组解 作迭代格式取初始值 ，并令 ，得 故 第38页/共53页6.3 6.3 非线性方程组求解数值法一般迭代格式写成向量形式记矩阵可以证明 时迭代收敛。(6-15)(6-16)第39页/共53页二、程序框图与通用程序设计 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)=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 非线性方程组求解数值法6.3.2 非线性方程组Gauss-yacobi-seidel迭代法 一、算法分析 仿照线性代数方程组的Gauss-Seidel迭代方法，可作非线性方程组Gauss-yacobi-seidel迭代法 为二、程序框图与通用程序设计(6-17)第42页/共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,n70 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 非线性方程组求解数值法6.3.3 非线性方程组最速下降迭代法(Gradient Iteration Method)一、算法分析 1.由已知方程组(6.12)式构造目标函数于是，方程组(6.12)式的解就上式的零极小值点，反之亦然。2.计算差商其中 ，(i=1,2,3,n)。式中c为控制常数，一般取c=0.00001。第45页/共53页6.3 6.3 非线性方程组求解数值法可以看出，梯度法实际上就是利用差商代替牛顿法中的偏导数。3.具体计算步骤（1）从给定的不全为零的初值 出发，设已经计算到第k,得 。（2）计算目标函数的值 。（3）如果|F|E)then x(i)=x(i)+cxs=s+df(i)*2x(i)=x(i)-cxsfl=f0/scall sub compfendif第48页/共53页6.3 6.3 非线性方程组求解数值法10 subroutine gradmt(n,E,x)20 dimension x(n),df(n)30 do i=1,n40 read(*,2x,f6.2)x(i)50 end do60 call compf(x,n,f)70 f0=f80 if(abs(f)E)then90 c=0.00001;s=0100 do i=1,n 110 if(abs(x(i)=0)then 120 cx=c 130 else140 cx=c*x(i)150 endif160 x(i)=x(i)+cx 170 call compf(x,n,f)180 df(i)=(f-f0)/cx190 s=s+df(i)*2 200 x(i)=x(i)-cx 210 end do220 sfl=f0/s230 do i=1,n 240 x(i)=x(i)-sfl*df(i)250 end do260 endif270 end subroutine gradmt第49页/共53页6.3 6.3 非线性方程组求解数值法10 subroutine compf(x,n,ff)20 dimension x(n),f(n)30 f(1)=expr140 f(2)=expr250 f(3)=expr360 .70 f(n)=exprn80 ff=090 do i=1,n 100 ff=ff+f(i)*f(i)110 end do120 end subroutine compf第50页/共53页练习题练习题6-1 方程 在 附近有根，将方程作3重通解变换，可得3种迭代式：判断各种迭代式在 附近的收敛性；选择一种收敛最快的迭代式，计算 附近的根，准确到4位小数。并和用牛顿法和牛顿弦割法计算的结果作比较。第51页/共53页上机实习题6-2 利用简单迭代法、加速迭代法、牛顿法和牛顿弦割法求解习题6-1题，准确到4位有效数字。6-2 利用yacobi和yacobi-seidel法求解下列方程组第52页/共53页感谢您的欣赏！第53页/共53页