(精品)Matlab非线性方程求根.ppt

《(精品)Matlab非线性方程求根.ppt》由会员分享，可在线阅读，更多相关《(精品)Matlab非线性方程求根.ppt（66页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、科学计算与科学计算与MATLABMATLAB中南大学材料科学与工程学院中南大学材料科学与工程学院中南大学材料科学与工程学院中南大学材料科学与工程学院 2010.102010.102010.102010.10第八讲非线性方程求根第八讲非线性方程求根内容提要内容提要n引言引言n二分法二分法n迭代法迭代法nNewton迭代法迭代法nMATLAB的非线性方程求根函数的非线性方程求根函数n小结小结 在工程和科学技术中许多问题常常归结为求解非线性方程式问在工程和科学技术中许多问题常常归结为求解非线性方程式问题，例如在控制系统的设计领域，研究人口增长率等。题，例如在控制系统的设计领域，研究人口增长率等。例例

2、 关于真实气体的状态方程（关于真实气体的状态方程（Van der waals方程）为方程）为 其中，其中，P是气体压力，是气体压力，V是气体体积，是气体体积，T是绝对温度，是绝对温度，R是气体常数。是气体常数。如果已知某气体的温度如果已知某气体的温度T及压力及压力P，那么求体积，那么求体积V的方程为：的方程为：1、引言、引言非线性方程的一般形式：非线性方程的一般形式：f(x)=0代数方程：代数方程：f(x)=a0+a1x+anxn（an 0）超越方程超越方程：f(x)中含三角函数、指数函数、或其他超中含三角函数、指数函数、或其他超越函数。越函数。用数值方法求解非线性方程的步骤：用数值方法求解非

3、线性方程的步骤：（1）找出有根区间；（只含一个实根的区间称隔根）找出有根区间；（只含一个实根的区间称隔根区间）区间）（2）近似根的精确化。从隔根区间内的一个或多个）近似根的精确化。从隔根区间内的一个或多个点出发，逐次逼近，寻求满足精度的根的近似值。点出发，逐次逼近，寻求满足精度的根的近似值。2、方程求根的二分法、方程求根的二分法2.1 二分法的基本思想：二分法的基本思想：假定假定f(x)=0在在a,b内有唯一单实根内有唯一单实根x*，考察有根区间，考察有根区间a,b，取中点，取中点x0=(a+b)/2，若，若f(x0)=0，则，则x*=x0，否则，否则，（1）若）若f(x0)f(a)0，则，则

4、x*在在x0右侧，令右侧，令a1=x0，b1=b;（2）若）若f(x0)f(a)0，则，则x*在在x0左侧，令左侧，令a1=a，b1=x0。定理定理1（介值定理）设函数（介值定理）设函数f(x)在区间在区间a,b连续，且连续，且f(a)f(b)0，则方程，则方程f(x)=0在区间在区间a,b内至少有一个根。内至少有一个根。以以a1,b1为新的隔根区间，且仅为为新的隔根区间，且仅为a,b的一半，对的一半，对a1,b1重复前过程，得新的隔根区间重复前过程，得新的隔根区间a2,b2，如此二分下去，得，如此二分下去，得一系列隔根区间：一系列隔根区间：a,b a1,b1 a2,b2 ak,bk 其中每个

5、区间都是前一区间的一半，故其中每个区间都是前一区间的一半，故ak,bk 的长度：的长度：当当k趋于无穷时趋于趋于无穷时趋于0。即若二分过程无限继续下去，这。即若二分过程无限继续下去，这些区间最后必收敛于一点些区间最后必收敛于一点x*，即方程的根。，即方程的根。每次二分后，取有根区间的中点每次二分后，取有根区间的中点xk=(ak+bk)/2作作为根的近似值，则可得一近似根序列：为根的近似值，则可得一近似根序列：x0,x1,x2,该序列必以根该序列必以根x*为极限。为极限。实际计算中，若给定充分小的正数实际计算中，若给定充分小的正数 0和允许误差和允许误差限限 1，当，当|f(xn)|0或或bn-

6、an0)disp(两端点函数值乘积大于两端点函数值乘积大于0!);return;else root=FindRoots(f,a,b,eps);endfunction r=FindRoots(f,a,b,eps)f_1=subs(sym(f),findsym(sym(f),a);f_2=subs(sym(f),findsym(sym(f),b);mf=subs(sym(f),findsym(sym(f),(a+b)/2);if(f_1*mf0)t=(a+b)/2;r=FindRoots(f,t,b,eps);else if(f_1*mf=0)r=(a+b)/2;else if(abs(b-a)e

7、ps)n=n+1;r1=root;root=subs(sym(f),findsym(sym(f),r1)+r1;tol=abs(root-r1);end实例实例例：求方程例：求方程 x3-x-1=0 在在 x=1.5 附近的一个根。附近的一个根。r,n=StablePoint(x3-x-1,1.5)r=NaNn=9实例实例例：求方程例：求方程 1/sqrt(x)+x-2=0 在在 x=1.5 附近的一个根。附近的一个根。r,n=StablePoint(1/sqrt(x)+x-2,1.5)r=Infn=1028 r,n=StablePoint(1/sqrt(x)+x-2,0.5)r=0.3820

8、n=4假定假定(x)改变不大，近似取某个近似值改变不大，近似取某个近似值L，则有，则有3.3 迭代收敛的加速方法迭代收敛的加速方法1、Aitken加速收敛方法：加速收敛方法：由微分中值定理，有由微分中值定理，有同理同理两式相比，得两式相比，得类推可得类推可得故故上式即为上式即为Aitken加速收敛方法的迭代格式。加速收敛方法的迭代格式。2、Steffensen迭代法：迭代法：将将Aitken加速技巧与不动点结合可得加速技巧与不动点结合可得或将其写为或将其写为3.4 MATLAB实现实现例：用不动点迭代法、例：用不动点迭代法、Aitken加速收敛方法和加速收敛方法和Steffensen迭代法迭代

9、法分别求方程，并比较这三种方法。分别求方程，并比较这三种方法。r,n=StablePoint(1/sqrt(x)+x-2,0.5)r=0.3820n=4 r,n=AtkenStablePoint(1/sqrt(x)+x-2,0.5)r=0.3820n=4 r,n=AtkenStablePoint(1/sqrt(x)+x-2,0.999)r=1.0000n=4 r,n=StablePoint(1/sqrt(x)+x-2,0.999)r=0.3820n=21 r,n=StevenStablePoint(1/sqrt(x)+x-2,0.999)r=1.0000n=2 r,n=StevenStable

10、Point(1/sqrt(x)+x-2,0.5)r=0.3820n=4表明即使不动点迭代法表明即使不动点迭代法不收敛，用不收敛，用steffensen迭代法仍可能收敛。迭代法仍可能收敛。说明说明steffensen迭代法迭代法的收敛速度比不动点迭的收敛速度比不动点迭代快得多。代快得多。谢谢！4、Newton迭代法迭代法4.1 基本思想：将非线性方程逐步归结为某种基本思想：将非线性方程逐步归结为某种线性方程求解。线性方程求解。设方程设方程f(x)=0有近似根有近似根xk（f(xk)0），将），将f(x)在在xk展开：展开：（在在x和和xk之间）之间）可设可设记该线性方程的根为记该线性方程的根为x

11、k+1，则，则故故f(x)=0可近似表示为可近似表示为即为即为Newton法迭代格式法迭代格式。（k=0,1,）4.2 Newton迭代法迭代法（切线法）（切线法）xyx*xk+1xkPky=f(x)切线方程切线方程故故切线法切线法的缺陷的缺陷1.被零除错误被零除错误2.程序死循环程序死循环y=arctan x方程方程:f(x)=x3 3x+2=0在重根在重根x*=1附近，附近，f(x)近似近似为零。为零。对对 f(x)=arctan x存在存在 x0，Newton迭代迭代法陷入死循环。法陷入死循环。计算输入计算计算 输出输出stop是是否否程序设计程序设计function root=Newt

12、onRoot(f,a,b,eps)if(nargin=3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(f1=0)root=a;endif(f2=0)root=b;endif(f1*f20)disp(两端点函数值乘积大于两端点函数值乘积大于0!);return;else tol=1;fun=diff(sym(f);fa=subs(sym(f),findsym(sym(f),a);fb=subs(sym(f),findsym(sym(f),b);dfa=subs(sym(fun

13、),findsym(sym(fun),a);dfb=subs(sym(fun),findsym(sym(fun),b);if(dfadfb)root=a-fa/dfa;else root=b-fb/dfb;end while(toleps)r1=root;fx=subs(sym(f),findsym(sym(f),r1);dfx=subs(sym(fun),findsym(sym(fun),r1);root=r1-fx/dfx;tol=abs(root-r1);endend例：采用切线法求方程在区间例：采用切线法求方程在区间0.5，2上的一个根。上的一个根。r=NewtonRoot(sqrt(

14、x)-x3+2,0.5,2)r=1.4759若若|(x)|=|1-cf(x)|1，即取，即取0cf(x)0)disp(两端点函数值乘积大于两端点函数值乘积大于0!);return;else tol=1;fun=diff(sym(f);fa=subs(sym(f),findsym(sym(f),a);fb=subs(sym(f),findsym(sym(f),b);dfa=subs(sym(fun),findsym(sym(fun),a);dfb=subs(sym(fun),findsym(sym(fun),b);if(dfadfb)df0=1/dfa;root=a-df0*fa;else df

15、0=1/dfb;root=b-df0*fb;end while(toleps)r1=root;fx=subs(sym(f),findsym(sym(f),r1);root=r1-df0*fx;tol=abs(root-r1);endend例：采用平行弦法求方程在区间例：采用平行弦法求方程在区间1.2，2上的一个根。上的一个根。r=SimpleNewton(sqrt(x)-x3+2,1.2,2)r=1.4759在在Newton迭代格式中，用差商近似导数，迭代格式中，用差商近似导数，4.4 弦截法（割线法）弦截法（割线法）称弦截法。称弦截法。得得弦截法的几何意义：弦截法的几何意义：xyx*xk+1

16、xk-1Pk-1y=f(x)xkPk弦线弦线PkPk-1的方程：的方程：当当y0时，时，程序设计程序设计function root=Secant(f,a,b,eps)if(nargin=3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(f1=0)root=a;endif(f2=0)root=b;endif(f1*f20)disp(两端点函数值乘积大于两端点函数值乘积大于0!);return;else tol=1;fa=subs(sym(f),findsym(sym(f),a)

17、;fb=subs(sym(f),findsym(sym(f),b);root=a-(b-a)*fa/(fb-fa);while(toleps)r1=root;fx=subs(sym(f),findsym(sym(f),r1);s=fx*fa;if(s=0)root=r1;else if(s0)root=b-(r1-b)*fb/(fx-fb);else root=a-(r1-a)*fa/(fx-fa);end end tol=abs(root-r1);endend例：采用弦截法求方程在区间例：采用弦截法求方程在区间1，4上的一个根。上的一个根。r=Secant(log(x)+sqrt(x)-2,

18、1,4)r=1.8773例例 用简化的用简化的Newton迭代法和弦截法计算方程迭代法和弦截法计算方程x3-3x+1=0的根。的根。解：设解：设f(x)=x3-3x+1，则，则f(x)=3x2-3由简化的由简化的Newton法，得法，得由弦截法，得由弦截法，得x0=0.5x1=0.3333333333x2=0.3497942387x3=0.3468683325x4=0.3473702799x5=0.3472836048x6=0.3472985550 x7=0.3472959759x8=0.3472964208x9=0.3472963440 x10=0.3472963572x11=0.34729

19、63553x0=0.5;x1=0.4;x2=0.3430962343x3=0.3473897274x4=0.3472965093x5=0.3472963553x6=0.3472963553简化简化Newton法法弦截法弦截法要达到精度要达到精度10-8，简化，简化Newton法法迭代迭代11次，弦截法迭代次，弦截法迭代5次，次，Newton迭代法迭代迭代法迭代4次次。无论前面哪种迭代法：无论前面哪种迭代法：（Newton迭代法、迭代法、简化简化Newton法、法、弦截法）弦截法）Newton迭代法迭代法x0=2x1=-3.54x2=13.95x3=-279.34x4=122017是否收敛均与初

20、值的位置有关。是否收敛均与初值的位置有关。如如x0=1x1=-0.5708x2=0.1169x3=-0.0011x4=7.9631e-010 x5=0收敛收敛发散发散为防止为防止Newton法发散，可增加一个条件法发散，可增加一个条件：|f(xk+1)|f(xk)|，满足该条件的算法称下山法。，满足该条件的算法称下山法。可用下山法保证收敛，可用下山法保证收敛，Newton法加快速度。法加快速度。4.5 Newton下山法下山法称称Newton下山法。下山法。（01,下山下山因子）因子）记记即即 的选取：从的选取：从=1开始，逐次减半计算，开始，逐次减半计算，即按即按的顺序，直到使下降条件的顺序

21、，直到使下降条件|f(xk+1)|r=NewtonDown(sqrt(x)-x3+2,1.2,2)r=1.4759例：求解方程例：求解方程要求达到精度要求达到精度|xn-xn-1|10-5，取，取x0=-0.99。解：先用解：先用Newton迭代法：迭代法：f(x)=x2-1x2=21.69118 x3=15.15689x4=9.70724 x5=6.54091x6=4.46497 x7=3.13384x8=2.32607 x9=1.90230 x10=1.75248 x11=1.73240 x12=1.73205 x13=1.73205需迭代需迭代13次才达次才达到精度要求到精度要求用用Ne

22、wton下山法，结果如下：下山法，结果如下：k=0 x0=-0.99 f(x0)=0.666567k=1 x1=32.505829 f(x)=11416.4 =0.5 x1=15.757915 f(x)=1288.5 =0.25 x1=7.383958 f(x)=126.8 =0.125 x1=3.196979 f(x)=7.69 =0.0625 x1=1.103489 f(x)=-0.655k=2 x2=4.115071 f(x)=19.1 =0.5 x2=2.60928 f(x)=3.31 =0.25 x2=1.85638 f(x)=0.27k=3 x3=1.74352 f(x)=0.02

23、3k=4 x4=1.73216 f(x)=0.00024k=5 x5=1.73205 f(x)=0.00000k=6 x6=1.73205 f(x)=0.000000k 下山因子下山因子xk f(xk)Matlab中中zeroin的算法实现是的算法实现是fzero.x=fzero(FUN,x0)%x0可以是数，或区间可以是数，或区间x=fzero(FUN,x0,options)x,fval=fzero(FUN,x0,options)x,fval,exitflag=fzero(FUN,x0,options)5、MATLAB的非线性方程求根函数的非线性方程求根函数例：例：求方程在求方程在x=0.5

24、附近的根。附近的根。x=fzero(x2+x-1,0.5)x=0.6180例：采用牛顿下山法求方程在区间例：采用牛顿下山法求方程在区间1.2，2上的一个上的一个根。根。x=fzero(sqrt(x)-x3+2,0 2)x=1.4759 r=NewtonDown(sqrt(x)-x3+2,1.2,2)r=1.4759实例：实例：Ti-15-3合金热变形显微组织的分形合金热变形显微组织的分形 分形理论是非线性科学的三大方向之一，是一种描述分形理论是非线性科学的三大方向之一，是一种描述和处理具有自相似性非规则图形及客观对象有效的分析方和处理具有自相似性非规则图形及客观对象有效的分析方法和建模工具，其

25、根本特征在于实现了测度观的转变。因法和建模工具，其根本特征在于实现了测度观的转变。因而，分形的出现为分析解决一些复杂、非线性系统问题提而，分形的出现为分析解决一些复杂、非线性系统问题提供了新的科学方法和认识方法，具有广阔和深远的实际应供了新的科学方法和认识方法，具有广阔和深远的实际应用前景。用前景。大量研究表明，金属显微组织具有分形特征，且可以大量研究表明，金属显微组织具有分形特征，且可以采用尺度不变的分形维数对其进行客观、全面地表征。一采用尺度不变的分形维数对其进行客观、全面地表征。一些研究工作者以金属显微组织中离散的点状、线状分布的些研究工作者以金属显微组织中离散的点状、线状分布的典型夹杂

26、物为研究对象，通过对不同等级夹杂物的分维测，典型夹杂物为研究对象，通过对不同等级夹杂物的分维测，验证了夹杂物等级与分形维数之间的关系。验证了夹杂物等级与分形维数之间的关系。首先要对显微组织的首先要对显微组织的TEM像进行像进行二值化处理。二值化处理。灰度图像的二值化方法主要有全局灰度图像的二值化方法主要有全局阈值法和局部阈值法。在目标和背阈值法和局部阈值法。在目标和背景的灰度差别较明显时，全局阈值景的灰度差别较明显时，全局阈值法效果比较好且计算速度较快。而法效果比较好且计算速度较快。而采用全局阈值法对图像进行二值化采用全局阈值法对图像进行二值化处理的关键是阈值的选取，其准确处理的关键是阈值的选

27、取，其准确性直接影响对图像中特征信息进行性直接影响对图像中特征信息进行描述分析的正确性。描述分析的正确性。1显微组织图像的二值化处理显微组织图像的二值化处理2显微组织分形维数的计算显微组织分形维数的计算分形维数是定量刻划分形体复杂结构的特征参数，表征了分形维数是定量刻划分形体复杂结构的特征参数，表征了分形体的不规则程度。材料科学中常用的测定分形维数的分形体的不规则程度。材料科学中常用的测定分形维数的方法有盒维数法、码尺法和小岛法。其中，盒维数法无论方法有盒维数法、码尺法和小岛法。其中，盒维数法无论是对曲线还是对曲线围成的面都适用，与图像物理含义关是对曲线还是对曲线围成的面都适用，与图像物理含义

28、关系不大，因而广泛应用于二维图像分形维数的计算。系不大，因而广泛应用于二维图像分形维数的计算。对于二维欧式空间的一系列点，用边长为对于二维欧式空间的一系列点，用边长为的正方形盒子覆盖的正方形盒子覆盖这些点，统计出至少包含一个点的方格数目这些点，统计出至少包含一个点的方格数目N()，在双对数，在双对数坐标系中用最小二乘法直线拟合数据点坐标系中用最小二乘法直线拟合数据点(ln,lnN()，若二者，若二者具有很好的线性关系，表明图像具有分形特性，且所得直线具有很好的线性关系，表明图像具有分形特性，且所得直线的斜率就是该图像分形维数的近似值。近似计算盒维数的斜率就是该图像分形维数的近似值。近似计算盒维

29、数DB为：为：在当前的分形研究中，正方形在当前的分形研究中，正方形边长为边长为=2n(n0)，即正方形边，即正方形边长的最小尺寸为一个像素。当长的最小尺寸为一个像素。当图像边界处不能保证有完整的图像边界处不能保证有完整的正方形盒子时，包含目标像素正方形盒子时，包含目标像素点的正方形盒子数目按下式计点的正方形盒子数目按下式计算：算：式中式中N0()为包含目标像素的完为包含目标像素的完整盒子的数目；整盒子的数目；A0()为覆盖图为覆盖图像的完整盒子的面积；像的完整盒子的面积；A为整个为整个图像的面积。图像的面积。变形条件对显微组织分形维数的影变形条件对显微组织分形维数的影响响当其它变形条件相同时，

30、变形程度、变当其它变形条件相同时，变形程度、变形温度和变形速率对显微组织形温度和变形速率对显微组织(放大倍数放大倍数为为2.8万倍万倍)分形维数的影响如图分形维数的影响如图4所示。所示。由图由图4可以看出，变形温度和变形速率一可以看出，变形温度和变形速率一定时，随着变形程度的增大，显微组织定时，随着变形程度的增大，显微组织的分形维数增大。变形程度较小时，合的分形维数增大。变形程度较小时，合金晶体中的位错密度也较小，且位错线金晶体中的位错密度也较小，且位错线大多长而直；随着变形程度的增加，位大多长而直；随着变形程度的增加，位错之间相互交割加，具有更高密度的位错之间相互交割加，具有更高密度的位错缠

31、结，位错线长度因缠结而逐步缩短，错缠结，位错线长度因缠结而逐步缩短，显微组织形貌更加复杂，因而分形维数显微组织形貌更加复杂，因而分形维数增大。增大。实例：采用人工神经网络预测轧后厚比的实例：采用人工神经网络预测轧后厚比的变化规律变化规律轧制复合中影响厚比变化的因素非常多，作用也很复杂，而且轧制复合中影响厚比变化的因素非常多，作用也很复杂，而且一般影响关系都是非线性的。因而，很难用简单的数学表达式一般影响关系都是非线性的。因而，很难用简单的数学表达式来描述轧后厚比的影响因素及变化规律。来描述轧后厚比的影响因素及变化规律。有人根据轧制复合实验数据进行了回归分析，得出了简单的回有人根据轧制复合实验数

32、据进行了回归分析，得出了简单的回归经验计算式。但是运用这种方法，计算工作量很大，而且预归经验计算式。但是运用这种方法，计算工作量很大，而且预测精度也不够高。回归分析也存在着局限性：一方面，回归分测精度也不够高。回归分析也存在着局限性：一方面，回归分析不能建立考虑所有影响因素的回归方程；另一方面，当加入析不能建立考虑所有影响因素的回归方程；另一方面，当加入新的实验数据后，又必须对模型中的常数重新进行回归求解。新的实验数据后，又必须对模型中的常数重新进行回归求解。因此回归常数使用范围很窄。因此回归常数使用范围很窄。有鉴于此，准备采用人工神经网络来处理实验数据。利用人工有鉴于此，准备采用人工神经网络

33、来处理实验数据。利用人工神经网络处理数据时，不需要给定特定的数学模型，而且其学神经网络处理数据时，不需要给定特定的数学模型，而且其学习的自适应性很强，即在原有网络基础上加入新得数据时，可习的自适应性很强，即在原有网络基础上加入新得数据时，可以调整网络状态，从而适应新的实验数据。以调整网络状态，从而适应新的实验数据。1人工神经网络概论人工神经网络概论人工神经网络是采用物理可实现的系统人工神经网络是采用物理可实现的系统来模仿人脑神经细胞的结构和功能。它来模仿人脑神经细胞的结构和功能。它由很多处理单元有机地联结起来，并行由很多处理单元有机地联结起来，并行地进行工作。它的处理单元很简单，但地进行工作。

34、它的处理单元很简单，但大量处理单元构成的网络系统却可以实大量处理单元构成的网络系统却可以实现极其丰富多彩的功能。现极其丰富多彩的功能。从神经网络结构上来看，多层感知器从神经网络结构上来看，多层感知器（Multi-layer perception orMLP）是具）是具有实际应用价值，而且目前应用最为广有实际应用价值，而且目前应用最为广泛的神经网络。它克服了单层感知器只泛的神经网络。它克服了单层感知器只能满足线性分类的局限性，在输入层和能满足线性分类的局限性，在输入层和输出层之间引入了若干个隐层，如图输出层之间引入了若干个隐层，如图2-5所示。因此，多层感知器可以很好地所示。因此，多层感知器可以

35、很好地完成函数的非线性估计。完成函数的非线性估计。3轧制复合轧后厚比的人工神经网络预测轧制复合轧后厚比的人工神经网络预测采用采用MATLAB神经网络工具箱中的神经网络工具箱中的trainB-Px函数的动量函数的动量算法和学习率自适应两项先进技术，可以提高学习速度，算法和学习率自适应两项先进技术，可以提高学习速度，并增加算法的可靠性。也可以采用并增加算法的可靠性。也可以采用trainlm函数中的函数中的Levenberg-Marquardt优化算法等。优化算法等。采用采用MATLAB语言，实现语言，实现B-P算法，程序流程如图算法，程序流程如图2-7所示。所示。根据多次试验和总结经验，本文采用的根据多次试验和总结经验，本文采用的B-P网络结构为网络结构为3-9-10-1型，如图型，如图2-8所示。网络的输入模式为所示。网络的输入模式为X=(i,Ti,)，输，输出模式为出模式为Y=(f)。学习样本从前面。学习样本从前面45组实验数据中选择组实验数据中选择27组组共共729个数据样本进行神经网络学习，其余个数据样本进行神经网络学习，其余18组数据用于对组数据用于对预测结果进行验证。预测结果进行验证。小小 结结非线性方程求根方法的几何意义非线性方程求根方法的几何意义非线性方程求根函数非线性方程求根函数的理解与应用的理解与应用谢 谢！