matlab非线性方程求解解读

1、非线性方程的解法1引言数学物理中的许多问题归结为解函数方程的问题，即，f(x)=O(1.1)这里，f (x)可以是代数多项式，也可以是超越函数。若有数x为方程f (x) = 0的根，或称函数f (x)的零点。设函数f (x)在a,b内连续，且f (a)f (b) 0。根据连续函数的性质知道，方程f(x)=O在区间a,b内至少有一个实根；我们又知道，方程 f(x)=O的根，除了极少简单方程的根可以 用解析式表达外，一般方程的根很难用一个式子表达。即使能表示成解析式的，往往也很复 杂，不便计算。所以，具体求根时，一般先寻求根的某一个初始近似值，然后再将初始近似 值逐步加工成满足精度要求为止。如何寻

2、求根的初始值呢？简单述之，为了明确起见，不妨设 f (x)在区间a,b内有一个实的单根，且f(a) prog1021计算结果如下。err =0.4589k =3err =0.1052k =4err =0.0292k =5err =0.0078k =6err =0.0021k =7err =5.7408e-004k =8err =1.5525e-004k =9err =4.1975e-005k =10err =1.1350e-005k =11err =3.0688e-006P = Columns 1 through 60.50000.95890.85370.88290.87510.8772Col

3、umns 7 through 110.87660.87680.87670.87670.8767ans =0.87673二分法3.1二分法原理二分法是方程求解最直观、最简单的方法。二分法以连续函数的介值定理为基础的。由介值定理知道，若函数f(x)区间a,b上连续，且f(a)* f (b) ：0，即f(a)和f(b)负号相反,23则f(x)在a,b内一定有实根。二分法的基本思想是：用对分区间的方法根据分点处函数f(x)的符号逐步将有限区间缩小，使在足够小的区间内，方程有且仅有一根。下面简述其基本步 骤。首先记aa.bo =b。用中点x。二也直将区间a,bo等分成2个小区间： 凶和2xo,bo。然后

4、分析可能存在的三种情况：如果f(a)* f(X。)=0，则x。是零点，也就是方程的根。如果f(a)* f(xo) ：0，则区间ao,x。内存在零点。如果f(xo)* f (b) : 0，则区间xo,bo内存在零点。对有根的新区间施行同样的操作，于是得到一系列有空的区间:a,bo二gb二 a二二ak，bk(3.1)其中每1个区间的长度都是前一区间长度的一半，最后 1个区间的长度为: bak ba2k如果取最后1个区间ak,bk的中点:(3.2)bk akx“ 2(3.3)作为f(x)=0根的近似值，则有误差估计式:. bk 理.b -a对于所给精度；，若取k使得b a则有，(3.4)(3.5)(

5、3.6)3.2程序与实例用二分法求解方程f(x)=0在有根区间a,b内的一个根，其中f(x)在a,b只有一个根的情形*function c,err,yc = bisect(f1031,a,b,delta)% f1031 是所要求解的函数。% a和b分别为有根区间的左右限。% delta 是允许的误差界。% c 为所求近似解。% yc 为函数 f 在 c 上的值。% err 是 c 的误差估计。ya = feval(f1031,a);yb = feval(f,b);if yb = 0, c = b； return end if ya*yb0,disp(a,b)不是有根区间); returnend

6、max1 = 1 + round(log(b-a) - log(delta)/log(2);for k = 1:max1c = (a + b)/2;yc = feval(f,c); if yc = 0 a=c; b=c;return elseif yb*yc0 b = c;yb = yc; else a = c; ya = yc;endif (b-a) prog1031计算结果如下。ans =1.32474牛顿法4.1牛顿法原理从前面迭代法，我们知道，迭代函数(x)构造的好坏，不仅影响收敛速度，而且迭代格式有可能发散。怎样选择一个迭代函数才能保证迭代序列一定收敛呢？构代迭代函数的一条重要途径是

7、用近似方程来代替原方程去求根。因此，如果能将非线 性方程(1.1)用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。牛顿 (Newt on)法就是一种将非线性方程线化的一种方法。设Xk是方程(1.1)的一个近似根，把如果f(x)在Xk处作一阶Taylor展开，即：f (x) ： f (xk) f(Xk)(x - Xk)(4.1)于是我们得到如下近似方程：f(xj f(Xk)(x-Xk) =0(4.2)设f(Xk)=0，则方程(10.2.1)的解为:(4.3)Xk.1 =Xk -空必,k =0,12(4.4)取作为原方程(1.1)的新近似根X-，即令:f(Xk)上式称为牛顿迭代格式

8、。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法牛顿法具有明显的几何意义。方程：(1045)y = f (兀)f(Xk)(x -Xk)是曲线y = f(X)上点(Xk，f(Xk)处的切线方程。迭代格式(4.4)就是用切线式(4.5)的零点来代替 曲线(1.1)的零点。正因为如此，牛顿法也称为切线法。牛顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。一般来说， 牛顿法对初值Xo的要求较高，初值足够靠近X*时才能保证收敛。若要保证初值在较大范围内 收敛，则需对f (X)加一些条件。如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛 顿法作一些改时，即可以采用下面的迭代

9、格式：Xk 1 = Xk， k = 0,1,2，(4.6)f (Xk)式中，0 ： ：1，称为下山因子。因此，用这种方法求方程的根，也称为牛顿下山法。牛顿法对单根收敛速度快，但每迭代一次，除需计算f(xj之外，还要计算f(Xk)的值。如果f(X)比较复杂，计算f(Xk)的工作量就可能比较大。为了避免计算导数值，我们可用差 商来代替导数。通常用如下几种方法：(1)割线法。如果用f (Xk) - f (Xk 4)Xk _ Xk 4代替f(Xk)，贝U得到割线法的迭代格式为:(4.7)拟牛顿法。如果用f ( xk ) - f ( xk f ( Xk) f (Xk)代替f(Xk)，贝U得到拟牛顿法的迭

10、代格式为:xk 1f2(Xk)f (Xk) - f 做 - f (Xk)(4.8)Steffe nson法。如果用f (Xkf (Xk) - f (Xk)f(Xk)代替f(Xk)，贝U得到拟牛顿法的迭代格式为:f 2(Xk)f(Xkf (Xk) - f(Xk)(4.9)4.2程序与实例1.牛顿法的程序f (x) =0。给定初值po，用牛顿法格式Pk -Pw12山，k =1,2/，求解非线性方程 f(P*fun ctio n p1,err,k,y = newto n(f1041,df1041,p0,delta,max1)% f1041是非线性函数。% df1041 是 f1041 的微商。% p

11、0是初始值。% delta是给定允许误差。% max1是迭代的最大次数。% p1是牛顿法求得的方程的近似解。% err是p0的误差估计。% k是迭代次数。% y = f(p1)p0, feval(f1041,p0)for k = 1:max1pl = p0 - feval(f1041, p0)/feval(df1041, pO);err = abs(pl-pO);pO = p1;p1, err, k, y = feval(f1041, p1)if (err prog1041计算结果如下。p0 =1.2000ans =0.1280p1 =1.1030err =0.0970k =1y =0.032

12、9p1 =1.1030err =0.0970k =1y =0.0329p1 =1.0524err =0.0507k =2y =0.0084p1 =1.0524err =0.0507k =2y =0.0084p1 =1.0264err =0.0260k =3y =0.0021p1 =1.0264err =0.0260k =3y =0.0021p1 =1.0133err =0.0131k =4y =5.2963e-004p1 =1.0133err =0.0131k =4y =5.2963e-004p1 =1.0066err =0.0066k =5y =1.3270e-004p1 =1.0066er

13、r =0.0066k =5y =1.3270e-004p1 =1.0033err =0.0033k =6y =3.3211e-005p1 =1.0033err =0.0033k =6y =3.3211e-005p1 =1.0017err =0.0017k =7y =8.3074e-006p1 =1.0017err =0.0017k =7y =8.3074e-006p1 =1.0008err =8.3157e-004k =8y =2.0774e-006p1 =1.0008err =8.3157e-004k =8y =2.0774e-006p1 =1.0004err =4.1596e-004k =

14、9y =5.1943e-007p1 =1.0004err =4.1596e-004k =9y =5.1943e-007p1 =1.0002err =2.0802e-004k =10y =1.2987e-007p1 =1.0002err =2.0802e-004k =10y =1.2987e-007p1 =1.0001err =1.0402e-004k =11y =3.2468e-008p1 =1.0001err =1.0402e-004k =11y =3.2468e-008p1 =1.0001err =5.2014e-005k =12y =8.1170e-009p1 =1.0001err =5

15、.2014e-005k =12y =8.1170e-009p1 =1.0000err =2.6008e-005k =13y =2.0293e-009p1 =1.0000err =2.6008e-005k =13y =2.0293e-009p1 =1.0000err =1.3004e-005k =14y =5.0732e-010p1 =1.0000err =1.3004e-005k =14y =5.0732e-010p1 =1.0000err =6.5020e-006k =15y =1.2683e-010p1 =1.0000err =6.5020e-006k =15y =1.2683e-010p

16、1 =1.0000err =3.2510e-006k =16y =3.1708e-011p1 =1.0000err = 3.2510e-006k =16y =3.1708e-011p1 =1.0000err =1.6255e-006k =17y =7.9270e-012p1 =1.0000err =1.6255e-006k =17y =7.9270e-012p1 =1.0000err =8.1277e-007k =18y =1.9817e-012ans1.0000这说明，经过 18 次迭代得到满足精度要求的值。2. 割线法的 MATLAB 实现*function p1, err, k, y = secant(f1042, p0, p1, delta, max1) % f1042 是给定的非线性函数。% p0,p1 为初始值。% delta 为给定误差界。% max1 是迭代次数的上限。% p1 为所求得的方程的近似解。% err 为 p1-p0 的绝对值。% k 为所需的迭代次数。% y=f(pi)p0, p1, feval(f1042, p0), feval(f1042,p1), k = 0;for k = 1:max1p2 = pl - feval(f1042,p1)*(p1-p0)/(feval(f1042,p1) - feval(f1042,p0); err = a