非线性方程的数值解法省公共课一等奖全国赛课获奖课件.pptx

第第5章章 非线性方程数值解法非线性方程数值解法5.1问题背景问题背景人口增加问题人口增加问题 设N(t)表示t 时刻人口数量，而表示人口固定出生率，那么人口满足以下微分方程：(5.11)上述方程解为：(5.12)其中N0表示初始人口数量。(5.12)式表示人口指数模型只有当人口是封闭没有迁移时是有效。假如允许以固定速率 迁移，那么微分方程应改为：(5.13)第1页该微分方程解为：(5.14)假设最初人口为1000万，即N0=1000万，第一年有435万人迁入，而在当年末时人口为1564万人，为了确定人口出生率，我们必须要求解下面关于方程：(5.15)则式（5.15）改为：第2页如图5.11所表示。方程（5.15）表明要求解当年出生率0，使满足 很显然，方程(5.14)是非线性方程，其解析解求解困难，只能求其近似解。本章讨论数值计算方法就是当准确解无法经过代数方法取得时用于求该类非线性方程近似解。普通而言，讨论非线性问题数值方法，与线性问题相比，不论从理论上还是计算方法上都要复杂多。非线性方程（5.16）解通常称为方程根，或称之为函数f(x)零点。第3页1.代数方程代数方程假如f(x)是多项式，即则称方程（5.16）为代数方程。2.超越方程超越方程 假如f(x)中含有三角函数，指数函数或其它超越函数，则称方程（5.16）为超越方程超越方程。比如：(5.17)(5.18)方程（5.17）是一个二次代数方程，它两个根能够表示为有限形式，即 方程（5.18）是一个超越方程，普通而言，必须用数值方法求它解。方程根可能是实数，也可能是复数，对应称之为实根和复根。第4页3.单根单根假如对于xk有则称xk为方程f(x)=0单根单根。4.重根重根假如有 则称xk为f(x)=0m+1次重根重根。对于代数方程，由代数学基本定理可知（实根或复根）个数与其次数相同，但对于超越方程却复杂多，假如有解，其解可能有一个或几个，也可能有没有穷多个。在大多数情况下，对于高于四次代数方程及超越方程没有准确准确求根公式。其实，在实际应用中并非一定要得到准确根，只要得到满足一定精度要求根近似值就能够了，就象导弹击中飞机一样，击中飞机各个部位都能够，并非一定要击中飞机发动机。下面将介绍几个求解普通非线性方程惯用数值计算方法。第5页5.2 二分法（二分法（the bisection method）设f(x)在区间a,b上连续，且f(a)f(b)0，那么依据连续函数零点定理，方程f(x)=0在(a,b)内最少有一个实根。为简便起见，不妨设f(x)=0在(a,b)内只有一个实根p。1.二分法基本思想二分法基本思想 用二分法求实根p近似值基本思绪就是：逐步将包含有根p区间二分，经过判断函数值符号，逐步对半缩小有根区间，直到区间缩小到允许误差范围之内，然后取区间中点为根p近似值。详细方法以下：为讨论方便，不妨设 如图5.21所表示。第6页Step1:区间a,b中点，再计算函数值 假如恰好 则表示待求根就是 不然判断 Step2:（即取原来区间左半部分）；（即取原来区间右半部）；这么二分区间a2,b2是方程新有根区间，它被包含在旧有根区间a1,b1即a,b之内，而且其长度仅是a1,b1二分之一。对缩小了区间a2,b2再计算其中点 Step3:（即取区间a2,b2左半部）；这么就把有根p区间缩小到更小范围a3，b3，再计算其中点二分以后求得有根p区间ak,bk中点为 如此重复重复上述计算过程，经过k次这么第7页而且（5.21）很显然，假如无限二分下去，则有因而零点p（根）为（5.22）但在实际计算时，不可能也完全没有必要进行这么无限计算过程，只要有根区间长度满足要求精度即可停顿计算，即满足：（5.23）则pk就是方程f(x)=0一个满足给定精度近似根了。其中是预先给定一个任意小正实数，如10-6，或者满足（5.24）则表示pk是方程f(x)=0一个满足给定精度要求近似根。第8页2.二分法算法源程序（二分法算法源程序（bisection.m）clear；format longa=input(a=?)；b=input(b=?)；tol=input(tol=?)；N=input(N=?)；n=1；fa=f(a)；flag=0；while(n=N)p=(a+b)./2；fp=f(p)；if(fp=0|(b-a)./20)a=p；fa=fp；else b=p end end if(flag=1)P=，P else Method failed after N iterations，N=，N end第9页说明说明：程序中函数f(x)应预先自定义，并取函数名存盘。以方程 为例，自定义一个名cubicf.m函数，源程序以下：function y=cubicf(x)y=x.3-2*x-5经过上述函数定义以后，将bisection.m源程序中f(a)改为cubicf(a)，f(p)改为cubicf(p)，那么就可运行bisection.m程序计算方程 根了。3.总总 结结 二分法优点是方法简单，编程轻易，且对函数性质要求不高，仅仅要求它连续且在区间两端点函数值异号即可。其收敛速度与公比为1/2等比数列相同。不过，二分法只能用于求实函数实根，不能用于求复根及偶数重根。第10页例例5.1 用二分法计算方程（5.15）中人口出生率。解解：经编程计算，人口出生率为0.10099887847900。源程序以下(ex51.m)：clear;format longa=0;b=1;tol=input(tol=?);N=input(N=?);n=1;fa=popu(a);flag=0;while(n=N)p=(a+b)./2;fp=popu(p);if(fp=0|(b-a)./20)a=p;fa=fp;else b=p;endif(flag=1)P=,p else Method failed after N iterations,N=,N endend第11页其中使用自定义函数（popu.m）为：function y=popu(x)if x=0y=-129;elsey=1000*exp(x)+435*(exp(x)-1)./x-1564;end第12页5.3迭代法迭代法 迭代法是数值计算中一类主要方法，是求解方程或方程组一类基本方法。在科学计算和工程实际中经常使用，尤其是伴随计算机普遍使用，使迭代法应用愈加广泛。1.迭代法基本思绪迭代法基本思绪 迭代法是一个逐次迫近方法，其基本思想是使用某个固定公式重复校正根近似值，从而得到一个近似根序列xk，使得该序列极限就是方程f(x)=0根。(1)对于非线性方程f(x)=0，用迭代法求根详细方法是：先将方程改写为便于迭代等价形式：（5.31）式中 称为迭代函数。因为（5.31）式是隐式方程，其右端含有未知x，因而不能直接求解，但假如用根某个猜测值x0代入式（5.31）右端，即可将它转化为显式计算公式：第13页若再取x1作为新猜测值，又有如此重复计算，其计算公式为（5.32）（5.32）式被称为迭代公式。假如迭代值xk有限，则称迭代收敛迭代收敛，这时极限值 显然就是方程(5.31)亦即原方程f(x)=0根。将方程f(x)=0等价转化为隐式方程，其实质就是求y=x 交点p，即 假如p是f(x)=0根，即f(p)=0，那么有 据此有（5.33）（5.33）就是依据f(x)结构 依据，从而得到（5.31）等价迭代形式。假如p是 y=x 交点p，那末，p一定是方程f(x)=0根，反之亦然。由此可知，迭代公式（5.32）也称定点迭定点迭代公式代公式。第14页2.线性迭代函数启示线性迭代函数启示 为了使迭代有效，必须确保它收敛性。一个发散（即不收敛）迭代过程，即使进行千百次迭代，其结果也毫无可用价值。为了确保迭代过程收敛性，迭代函数 结构十分关键。以最简单线性迭代函数为例能够轻易验证 迭代才是收敛，即 是确保迭代收敛线性 迭代函数基本特征（由作图法很轻易验证这一特征）。3.压缩映像原理压缩映像原理对于迭代函数 普通情形，设p为方程 根，则由微分中值定理有式中 是p与xk之间某一点，即 由此可知，假如存在，使得对于任意 成立：第15页则有（5.34）如此重复递推，对迭代误差 即迭代收敛。需要指出是，在上述论证过程中，应该确保一切迭代值xk全落在区间a,b内，为此要求对任意 总而言之有以下压缩映像原理压缩映像原理：定理定理5.1 设 在a,b上含有连续一阶导数，且满足以下两个条件：使得对于任意 成立 （5.35）第16页则迭代过程 对任意初值 均收敛于方程 根p，且有以下误差公式：（5.36）（5.37）第17页n=n+1；p0=p；%更新p0endif flag=1p=，p%The procedure was successful.elseThe method failed after N iterations N=，Nend4.定点迭代法源程序（定点迭代法源程序（fixedp.m）clear；format long p0=input(p0=?)；Tol=input(tolerance Tol=?)；N=input(maximum number of iterations N=?)；n=1；flag=0；While n=N p=f(P0)；if abs(p-p0)tolflag=1；break；end第18页例例5.2 求方程 唯一正根。解解：因为f(1)=-1，f(2)=5，所以区间1,2含有所求根。为使用迭代，将所给方程改写为：这时迭代函数 在1,2上恒有，由定理5.1可知，迭代公式 对任意初始值 均收敛。取初值x0=1.5，迭代结果如表5.1所表示。表5.1 迭代结果k012345678xk1.5 1.357211.330861.325881.324941.324761.324731.324721.32472第19页例例5.3用迭代法求解方程 根（已知方程一个根为解解：把方程改写为以下三种等价迭代形式：对于（对于（1）：）：所以用 作迭代函数时，无法收敛。对于（对于（2）：）：所以用 作迭代函数也无法收敛。对于（对于（3）：）：所以可用其作迭代函数。以x0=3为初值，迭代结果以下：第20页经过2次迭代即可到达3位有效数字精度。若将f(x)=0改为 所以用 为迭代函数时也收敛，结果以下：由以上计算结果可知，两种迭代函数都可收敛，但收敛速度却不相同，前者很快收敛，而后者却收敛较慢，所以，结构出适当迭代函数 十分主要。5.迭代过程收敛速度迭代过程收敛速度 一个迭代法要含有实用价值，不但需要必定它是收敛，而且还要求它收敛快。所谓收敛速度收敛速度是指在靠近收敛时迭代误差下降速度。详细而言就是，假如迭代误差 时成立：则称迭代过程是p阶收敛。当p=1时称线性收敛；p=2时称平方收敛第21页定理定理5.2 设 根p邻近有连续二阶导数，且 时迭代过程 为线性收敛；而当 时为平方收敛。第22页5.4迭代过程加速收敛方法迭代过程加速收敛方法1.迭代公式加工迭代公式加工 对于收敛迭代过程，只要迭代足够屡次，就能够使迭代结果到达任意精度，但有时迭代过程收敛较慢，从而大大增加了计算量，所以，迭代过程加速是个主要研究课题。设xk是根p某个近似值，用迭代公式校正一次得在所考虑范围内改变不大，其估值为L，则由微分中值定理有(5.41)从而得这表明，假如将迭代值 与xk加权平均，可望所得到 第23页是比 更加好近似根。其加速计算过程是：迭代迭代：加工加工：它们合并为：(5.42)例例5.4用迭代法和加速迭代法求方程 附近一个根p，要求精度为10-5。解解：设迭代函数为 所以迭代方程 是收敛。(1)使用非加速迭代时，经过18次迭代即得到满足精度要求根0.567141（准确值为0.56713）第24页结果如表5.2所表示。k01234567xk0.5000000.6065310.5452390.5797030.5600650.5711720.5648630.568438k89101112131415xk0.5664090.5675600.5669070.5672770.5670670.5671860.5671190.567157k161718xk0.5671350.5671480.567141(2)使用加速迭代公式(5.42)时，故取L=-0.6，此时加速迭代公式为：计算结果以下：即只要迭代3次即可得到满足精度要求结果，显然加速效果显著。第25页2.埃特金算法埃特金算法 上述加速方案缺点是，因为其中含有导数 相关信息而不便于实际应用。一样假定 在所考虑范围内改变不大，其估值为L，将迭代值 再迭代一次，即得到 由微分中值定理得利用式(5.41)：联立消去L得从而得若以上式右端得出结果作为新改进值，那末这么结构出加速公式不再含相关于导数信息第26页不过，它需要两次迭代值 进行加工，详细算法以下：迭代迭代：再迭代再迭代：加工加工：(5.43)上述方法称之为埃特金（埃特金（Aitken）加速方法）加速方法。3.埃特金加速算法源程序埃特金加速算法源程序(aitken.m)clear；format long x0=input(x0=?)；N=input(N=?)；Tol=input(Tol=?)；n=1；flag=0；While n=N x1=(x0)；x2=(x1)；x2=x2(x2x1).2./(x2-2*x1+x0)；if abs(x2-x0)Tolflag=1；breakend n=n+1；x0=x2；end if flag=1 x2=,x,The procedure was successful else The method failed after N iterations,N=,Nend 第27页5.5 牛顿迭代法牛顿迭代法 牛顿迭代法是求解非线性方程f(x)=0一个主要和惯用迭代方法，其基本思想是将非线性函数f(x)=0逐步线性化，从而将非线性方程f(x)=0近似地转化为线性方程求解。1.牛顿迭代公式导出牛顿迭代公式导出对于方程f(x)=0设已知它近似根为xk，则函数f(x)=0在点xk处作泰勒展开：取其前两项近似代替f(x)（称为f(x)线性化），即用线性方程(5.51)来近似方程f(x)=0。很显然，方程(5.51)是线性方程，它求根是很轻易。第28页则用线性方程(5.51)根作为非线性方程f(x)=0新近似根，记为xk+1，则有牛顿迭代公式牛顿迭代公式：(5.52)这就是著名牛顿迭代公式牛顿迭代公式，对应迭代函数是：(5.53)假如p是方程f(x)=0一个单根，即f(p)=0，由此可见，在单根p附近，对于任意初值x0，由收敛性定理5.1可知，由迭代公式(5.52)得到迭代序列都收敛于根p，而且收敛速度很快（由定理5.2知）。第29页 牛顿法有显著几何意义。方程f(x)=0根p在几何上表示曲线y=f(x)与x轴交点p。当求得p近似值xk以后，过曲线y=f(x)上对应点(xk，f(xk)作f(x)切线，其切线方程为：则该切线与x轴交点横坐标恰好是 这就是牛顿迭代公式(5.52)计算结果xk+1。继续取点(xk+1，f(xk+1)，再作曲线y=f(x)切线与x轴相交，又可得xk+2，由图5.51可知，只要初值x取得充分靠近根p，序列xk就会很快收敛于p。正因为牛顿法有这一显著几何意义，所以牛顿法牛顿法也常称之为切线法切线法。第30页2.牛顿法收敛性牛顿法收敛性 对于方程f(x)=0，假如f(x)在根p邻近含有连续二阶导数，且p是f(x)=0一个单根，则在跟p附近，对于任意初始值x0，由牛顿迭代法产生序列xk收敛于p，所以牛顿法含有局部收敛性，而且可深入证实牛顿迭代法在单根p附近是平方收敛，即含有二阶收敛速度。不过，假如p是f(x)=0重根时，则可证实牛顿法仅有线性收敛速度，且重数越高收敛越慢。牛顿迭代法对初值x0选取要求比较高，即x0必须充分靠近p才能确保局部收敛。为了确保牛顿法含有非局部收敛性，下面不加证实地给出一个判断牛顿非局部收敛充分定理：第31页定理定理5.3 设函数f(x)在区间a，b上存在二阶导数，且满足以下条件：在a，b上不为零；在a，b上不变号；(4)在a，b上任意选取满足条件 初始近似值x0。则有牛顿法产生序列xk单调收敛于方程f(x)=0在a,b上唯一根3 牛顿迭代法源程序（牛顿迭代法源程序（newtoniter.m）clear；format longx0=input(x0=?)；N=input(N=?)；TOL=input(TOL=?)；k=1；flag=0；While n=Nfdot=f(x0)；第32页if fdot=0 flag=0；breakendx1=x0f(x0)./fdot；if abs(x1x0)tolflag=1；breakendif n=Nflag=2；breakendn=n+1；x0=x1；endif flag=0 This equation has no solution；else if flag=1 p=，x1else This method failed after N iterations,N=，Nend第33页例例 5.5 用牛顿法求解非线性方程 解解：牛顿迭代公式为：所以方程 在0，1内有解。取初始值x=0.5，计算结果以下：与例例5.4相比可知，牛顿迭代法收敛速度与加速迭代法一样快。第34页例例5.6用牛顿迭代法（5.52）式计算 近似值。解解 准确值为2.4494897。则方程 解就是 近似值。所以方程 在2，3内有解。牛顿迭代公式以下：取初值x0=2.5，计算结果以下：第35页 由牛顿法收敛性定理5.3可知，牛顿法对初始值选取要求是很高。普通而言，牛顿法只有局部收敛性。当初始值取得离根p太远时，迭代将发散，而一旦xk进入收敛域内，牛顿法就有平方收敛速度。为了扬长避短，扩大初始值x0选择范围，下面介绍牛顿法一个改进方法牛顿下山法牛顿下山法。4.牛顿下山法牛顿下山法将牛顿迭代公式（5.52）修改为：（5.54）其中，是一个参数，且 称为下山因子下山因子。适当选取下山因子能够使单调性条件（5.55）成立。第36页 下山因子选择是个逐步探索过程，从=1开始重复将因子值减半进行计算，一旦单调性条件（5.55）成立，则称“下山成功下山成功”；不然，假如在上述过程中找不到使条件（5.55）成立下山因子，则称“下山失败下山失败”，这时需要另选初值x0重新计算。例例5.7 已知非线性方程一个根为p=1.32472。若取初值x0=0.6，用牛顿法 计算得到第一次近似值为x1=17.9，反而比x0=0.6更偏离了根p。若改用牛顿下山法：计算，仍取x0=0.6，从=1开始逐次搜索，当搜索到下山因子由牛顿下山公式（5.54）得到计算值第37页已满足条件（5.55），即 因而牛顿下山法修正了原来x1=17.9严重偏差，使迭代收敛。计算结果如表5.3所表示：k0123411/32111xk0.61.140631.366811.326281.32472由例例5.7分析能够知道，假如初值x0选取偏离根p太远，从而使迭代计算值x1偏离根p更远，那末就可采取牛顿下山法采取牛顿下山法，逐次将下山因子减半，重新计算x1值，即第38页并判断单调性条件 是否成立。若条件仍不成立，再将减半，重新计算x1，直到单调性条件成立为止，这么即可将第一次迭代计算值x1出现偏差得到修正，使之进入根p收敛区域之内，以后迭代过程，下山因子恒定为1。第39页5.6弦截法弦截法 牛顿法突出优点是收敛速度快，但也有显著缺点，这就是需要提供导数值f(xk)。假如函数f(x)=0比较复杂，致使导数计算困难，那末用牛顿法公式是很不方便。1.为了避开导数计算，能够改用差商 来代替牛顿公式（5.52）中导数，从而得到以下离散形式：（5.61）这个公式是依据f(x)=0等价形式：（5.62）第40页建立迭代公式，迭代函数为：（5.63）迭代公式（5.61）几何解释如图5.61所表示，曲线y=f(x)上横坐标为xk点记为pk，则差商 表示弦线 斜率，由公式（5.61）计算得到xk+1，实际上就是弦线 与x轴 交点，因而这种方法称为弦截法弦截法。图5.61 弦截法示意图第41页2.弦截法收敛性弦截法收敛性为提升收敛速度，改用差商（5.52）中导数f(xk)，从而导出以下迭代公式：来代替牛顿公式（5.64）这种迭代法称为快速弦截法快速弦截法 快速弦截法也称两步法两步法，这是因为在计算xk+1时，要用到前面两步信息xk-1和xk，所以在使用快速弦截法时，在计算之前必须预先提供两初始值x0和x1。第42页例例5.8 用快速弦截法求解非线性方程解解：所以方程 取初值 用快速弦截法求得结果以下：很显然，快速弦截法收敛速度是很快，与牛顿下山法(例例5.7)和加速迭代法（例例5.4）收敛速度相当。例例5.9 编程计算以下方程（见（5.11）式）根。解解：因为 所以方程f(x)=0在0，1内有解。第43页采取快速弦截法计算，程序以下：(1)f()函数定义（population2.m）function y=popu(x)n=length(x)；for j=1:n if x(j)=0 y(j)=-129；else y(j)=1000*exp(x(j)+435*(exp(x(j)-1)./x(j)-1564；endend第44页（2）求方程f()=0解源程序（secant.m）clear；format longx0=input(x0=?)；x1=input(x1=?)；N=input(N=?)；Tol=input(Tol=?)；n=2;flag=0；y0=population2(x0)；y1=population2(x1)；while n=N x=x1-y1.*(x1-x0)./(y1-y0)；if abs(x-x1)Tol flag=1；break end n=n+1；x0=x1；x1=x；y0=y1；y1=population2(x1)；endif flag=0 The method failed after N iterations，N=，Nelse x=，xend取初值 精度Tol10-6时，经六次迭代计算，计算结果为：0.1009979。第45页