Matlab非线性方程求根PPT课件

1、 在工程和科学技术中许多问题常常归结为求解非线性方程式问在工程和科学技术中许多问题常常归结为求解非线性方程式问题，例如在控制系统的设计领域，研究人口增长率等。题，例如在控制系统的设计领域，研究人口增长率等。 例例 关于真实气体的状态方程（关于真实气体的状态方程（van der waals方程）为方程）为 其中，其中，p是气体压力，是气体压力，v是气体体积，是气体体积，t是绝对温度，是绝对温度，r是气体常数。是气体常数。 如果已知某气体的温度如果已知某气体的温度t及压力及压力p，那么求体积，那么求体积v的方程为：的方程为： rtbvvap)(20)()(2rtbvvapxf1、引言、引言)(2v

2、gbvaprtv定理定理1（介值定理）设函数（介值定理）设函数f(x)在区间在区间a,b连续，且连续，且f(a)f(b)0，则方程则方程f(x)=0在区间在区间a,b内至少有一个根。内至少有一个根。)(21ababkkk每次二分后，取有根区间的中点每次二分后，取有根区间的中点xk= (ak+bk) /2作作为根的近似值，则可得一近似根序列：为根的近似值，则可得一近似根序列： x0, x1, x2, 该序列必以根该序列必以根x*为极限。为极限。实际计算中，若给定充分小的正数实际计算中，若给定充分小的正数 0和允许误差和允许误差限限 1，当，当|f(xn)| 0或或bn- an0) disp(两端

3、点函数值乘积大于两端点函数值乘积大于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)eps) n=

4、n+1; r1=root; root=subs(sym(f),findsym(sym(f),r1)+r1; tol=abs(root-r1);end r,n=stablepoint(x3-x-1,1.5)r = nann = 9 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.3820n = 4*)(*)()(*001xxxxxx*)(*01xxlxx*)(*12xxlxx*1021xxxxxxxx1、aitken加速收敛方法：加速收敛方法：由微分中值定理，有由微

5、分中值定理，有同理同理两式相比，得两式相比，得0012201001221202)(2*xxxxxxxxxxxxxx21212)(kkkkkkkxxxxxxx故故上式即为上式即为aitken加速收敛方法的迭代格式。加速收敛方法的迭代格式。)(kkxy)(kkyzkkkkkkkxyzxyxx2)(21或将其写为或将其写为)(1kkxxxxxxxxx)(2)()()(2例：用不动点迭代法、例：用不动点迭代法、 aitken加速收敛方法和加速收敛方法和steffensen迭代迭代法分别求方程，并比较这三种方法。法分别求方程，并比较这三种方法。120 xx r,n=stablepoint(1/sqrt(

6、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=stevenstablepoint(1/sqrt(x)+x-2,0.5)r = 0.3820

7、n = 4表明即使不动点迭代表明即使不动点迭代法不收敛，用法不收敛，用steffensen迭代法仍可迭代法仍可能收敛。能收敛。说明说明steffensen迭代法迭代法的收敛速度比不动点的收敛速度比不动点迭代快得多。迭代快得多。4.1 基本思想：将非线性方程逐步归结为某种基本思想：将非线性方程逐步归结为某种线性方程求解。线性方程求解。2)(! 2)()()()(kkkkxxfxxxfxfxf )()()(kkkxxxfxfxf可设可设0)()(kkkxxxfxf)()(1kkkkxfxfxx故故f(x)=0可近似表示为可近似表示为即为即为newton法迭代格式法迭代格式。（k=0,1,）)()(

8、1kkkkxfxfxxxyx*xk+1xkpky=f(x)切线方程切线方程)()(kkkxxxfxfy故故切线法切线法的缺陷的缺陷1.被零除错误被零除错误2.程序死循环程序死循环-6-4-20246-200-1000100200 xx3-3 x+2y = arctan x方程方程: f(x)=x3 3x + 2 = 0在重根在重根x*=1附近，附近，f(x)近似近似为零。为零。对对 f(x) = arctan x存在存在 x0，newton迭代迭代法陷入死循环。法陷入死循环。 计算输入计算计算 输出输出stop,00nx0, 2 , 1nk)(0 xf)(0 xf0)(0 xf1kx)()(0

9、00 xfxfxx0 xxxx 020i10i00ikxfx),(,0illl是是否否0nk function root=newtonroot(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(sy

10、m(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) 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;

11、tol=abs(root-r1); endend例：采用切线法求方程在区间例：采用切线法求方程在区间0.5，2上的一个根。上的一个根。220 xx r=newtonroot(sqrt(x)-x3+2,0.5,2)r = 1.4759)(1kkkxcfxx)()(xcfxx迭代公式：迭代公式：（c 0,k=0,1,）迭代函数：迭代函数：function root=simplenewton(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

12、(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),findsym(sym(fun),a); dfb=subs(sym(fun),findsym(sym(fun),b); if(dfadfb) df0=1/dfa; root=a-df0*fa; else

13、 df0=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上的一个根。上的一个根。220 xx r=simplenewton(sqrt(x)-x3+2,1.2,2)r = 1.475911)()()(kkkkkxxxfxfxf)()()()(111kkkkkkkxxxfxfxfxx称弦截法。称弦截法。得得)()()()(11kkkk

14、kkxxxxxfxfxfyxyx*xk+1xk-1pk-1y=f(x)xkpk弦线弦线pkpk-1的方程：的方程：当当y0时，时，)()()()(111kkkkkkkxxxfxfxfxxfunction 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!); retu

15、rn;else tol=1; fa=subs(sym(f),findsym(sym(f),a); 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例：采用弦截法

16、求方程在区间例：采用弦截法求方程在区间1，4上的一个根。上的一个根。lg2xx r=secant(log(x)+sqrt(x)-2,1,4)r = 1.8773313)()()()(21123111kkkkkkkkkkkkkkxxxxxxxxxxfxfxfxx解：设解：设f(x)=x3-3x+1，则，则f (x)=3x2-3由简化的由简化的newton法，得法，得)()(01xfxfxxkkk3313203xxxxkkk由弦截法，得由弦截法，得x0=0.5x1= 0.3333333333x2 = 0.3497942387x3 = 0.3468683325x4 = 0.3473702799x5

17、= 0.3472836048x6 = 0.3472985550 x7 = 0.3472959759x8 = 0.3472964208x9 = 0.3472963440 x10 = 0.3472963572x11 = 0.3472963553x0=0.5;x1=0.4;x2 = 0.3430962343x3 = 0.3473897274x4 = 0.3472965093x5 = 0.3472963553x6 = 0.3472963553简化简化newton法法弦截法弦截法要达到精度要达到精度10-8，简化，简化newton法法迭代迭代11次，弦截法迭代次，弦截法迭代5次，次， newton迭代法

18、迭代迭代法迭代4次。次。无论前面哪种迭代法：无论前面哪种迭代法：（newton迭代法、迭代法、 简化简化newton法、法、弦截法）弦截法）0)arctan()(xxf0 x精确解为newton迭代法迭代法)1(arctan21kkkkxxxx10 x取初值x0 = 2x1 = -3.54x2 = 13.95x3 = -279.34x4 = 122017是否收敛均与初值的位置有关。是否收敛均与初值的位置有关。如如20 x取初值x0 =1x1 = -0.5708x2 = 0.1169x3 = -0.0011x4 = 7.9631e-010 x5 = 0收敛收敛发散发散)()(1kkkkxfxfx

19、x称称newton下山法。下山法。)()(1kkkkxfxfxxkkkxxx)1 (11（01,下山因下山因子）子）记记即即,21,21,21, 132的顺序，直到使下降条件的顺序，直到使下降条件|f(xk+1)| r=newtondown(sqrt(x)-x3+2,1.2,2)r = 1.4759033 xx要求达到精度要求达到精度|xn-xn-1|10-5，取，取x0=-0.99。解：先用解：先用newton迭代法：迭代法：f (x)=x2-1)()(1kkkkxfxfxx)1(3323kkkkxxxx)1(332003001xxxxx505829.32x2=21.69118 x3=15.

20、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次才达次才达到精度要求到精度要求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.196

21、979 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.023k = 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的算法实现是的算法实现

22、是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附近的根。附近的根。210 xx x=fzero(x2+x-1,0.5)x = 0.6180例：采用牛顿下山法求方程在区间例：采用牛顿下山法求方程在区间1.2，2上的一个上的一个根。根。220 xx x=fzero(sqrt(x)-

23、x3+2,0 2)x = 1.4759 r=newtondown(sqrt(x)-x3+2,1.2,2)r = 1.4759首先要对显微组织的首先要对显微组织的tem像进行像进行二值化处理。二值化处理。灰度图像的二值化方法主要有全局灰度图像的二值化方法主要有全局阈值法和局部阈值法。在目标和背阈值法和局部阈值法。在目标和背景的灰度差别较明显时，全局阈值景的灰度差别较明显时，全局阈值法效果比较好且计算速度较快。而法效果比较好且计算速度较快。而采用全局阈值法对图像进行二值化采用全局阈值法对图像进行二值化处理的关键是阈值的选取，其准确处理的关键是阈值的选取，其准确性直接影响对图像中特征信息进行性直接影

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

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

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

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

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

29、且预测精度也不够高。回归分析也存在着局限性：一方面，回归分测精度也不够高。回归分析也存在着局限性：一方面，回归分析不能建立考虑所有影响因素的回归方程；另一方面，当加入析不能建立考虑所有影响因素的回归方程；另一方面，当加入新的实验数据后，又必须对模型中的常数重新进行回归求解。新的实验数据后，又必须对模型中的常数重新进行回归求解。因此回归常数使用范围很窄。因此回归常数使用范围很窄。有鉴于此，准备采用人工神经网络来处理实验数据。利用人工有鉴于此，准备采用人工神经网络来处理实验数据。利用人工神经网络处理数据时，不需要给定特定的数学模型，而且其学神经网络处理数据时，不需要给定特定的数学模型，而且其学习的自适应性很强，即在原有网络基础上加入新得数据时，可习的自适应性很强，即在原有网络基础上加入新得数据时，可以调整网络状态，从而适