matlab非线性方程求解

龙纯鹏

1引 言

数学物理中的许多问题归结为解函数方程的问题，即，

f(x)0 (1.1) 这里，f(x)可以是代数多项式，也可以是超越函数。若有数x*为方程f(x)0的根，或称函数f(x)的零点。

设函数f(x)在[a,b]内连续，且f(a)f(b)0。根据连续函数的性质知道，方程f(x)0在区间[a,b]内至少有一个实根；我们又知道，方程f(x)0的根，除了极少简单方程的根可以用解析式表达外，一般方程的根很难用一个式子表达。即使能表示成解析式的，往往也很复杂，不便计算。所以，具体求根时，一般先寻求根的某一个初始近似值，然后再将初始近似值逐步加工成满足精度要求为止。

如何寻求根的初始值呢？简单述之，为了明确起见，不妨设f(x)在区间[a,b]内有一个实的单根，且f(a)0,f(b)0。我们从左端出点x0a出发，按某一预定的步长h一步一步地向右跨，每跨一步进行一次根的“搜索”，即检查每一步的起点xk和xk1(即，xkh)的函数值是否同号。若有：

f(xk)*f(xkh)0 (1.2) 那么所求的根必在(xk,xkh)内，这时可取xk或xkh作为根的初始近似值。这种方法通常称为“定步长搜索法”。另外，还是图解法、近似方程法和解析法。

2 迭代法

2.1 迭代法的一般概念

迭代法是数值计算中一类典型方法，不仅用于方程求根，而且用于方程组求解，矩阵求

1

特征值等方面。迭代法的基本思想是一种逐次逼近的方法。首先取一个精糙的近似值，然后用同一个递推公式，反复校正这个初值，直到满足预先给定的精度要求为止。

对于迭代法，一般需要讨论的基本问题是：迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。这里，主要看看解方程迭代式的构造。 对方程(1.1)，在区间[a,b]内，可改写成为：

x(x) (2.1) 取x0[a,b]，用递推公式：

xk1(xk), k0,1,2, (2.2) 可得到序列：

x0,x1,x2,xk,{xk}k0 (2.3)

~~当k时，序列{xk}k0有极限x，且(x)在x附近连续，则在式(2.2)两边极限，得，

~x(~x)

即，~x为方程(2.1)的根。由于方式(1.1)和方程(2.1)等价，所以， x*~x 即，

limxkx*

k式(2.2)称为迭代式，也称为迭代公式；(x)可称为迭代函数。称求得的序列{xk}k0为迭代序列。

2.2 程序和实例

下面是基于MATLAB的迭代法程序，用迭代格式pn1g(xn)，求解方程xg(x)，其中初始值为p0。

**************************************************************************

function[p,k,err,P]=fixpt(f1021,p0,tol,max1) % f1021是给定的迭代函数。

2

% p0是给定的初始值。 % tol是给定的误差界。

% max1是所允许的最大迭代次数。 % k是所进行的迭代次数加1。 % p是不动点的近似值。 % err是误差。 % P = {p1,p2,…,pn}

P(1) = p0;

for k = 2:max1

P(k) = feval('f1021', P(k-1)); k, err = abs(P(k) - P(k-1))

p = P(k); if(errend

if k == max1

disp('maximum number of iterations exceeded');

end end P=P;

**************************************************************************** 例2.1 用上述程序求方程sinxx20的一个近似解，给定初始值x00.5，误差界为

105。

解：先用m文件先定义一个名为f1021.m的函数文件。 function y = f1021(x)

y = sin(x)/x;

建立一个主程序prog1021.m

clc

clear all

fixpt('f1021',0.5,10^(-5),20)

然后在MATLAB命令窗口运行上述主程序，即：

>> prog1021 计算结果如下。

k = 2

3

err = 0.4589 k = 3 err = 0.1052 k = 4 err = 0.0292 k = 5 err = 0.0078 k = 6 err = 0.0021 k = 7 err = 5.7408e-004 k = 8 err = 1.5525e-004 k = 9 err = 4.1975e-005 k = 10 err = 1.1350e-005 k = 11 err = 3.0688e-006 P = Columns 1 through 6

0.5000 0.9589 0.8537 0.8829 0.8751 0.8772 Columns 7 through 11

0.8766 0.8768 0.8767 0.8767 0.8767 ans = 0.8767

3 二分法

3.1 二分法原理

二分法是方程求解最直观、最简单的方法。二分法以连续函数的介值定理为基础的。由介值定理知道，若函数f(x)区间[a,b]上连续，且f(a)*f(b)0，即f(a)和f(b)负号相反，

4

则f(x)在[a,b]内一定有实根。二分法的基本思想是：用对分区间的方法根据分点处函数f(x)的符号逐步将有限区间缩小，使在足够小的区间内，方程有且仅有一根。下面简述其基本步骤。

首先记a0a,b0b。用中点x0a0b0将区间[a0,b0]等分成2个小区间：[a0,x0]和2[x0,b0]。然后分析可能存在的三种情况：

如果f(a)*f(x0)0，则x0是零点，也就是方程的根。 如果f(a)*f(x0)0，则区间[a0,x0]内存在零点。

如果f(x0)*f(b)0，则区间[x0,b0]内存在零点。

对有根的新区间施行同样的操作，于是得到一系列有空的区间：

[a0,b0][a1,b1][a2,b2][ak,bk] (3.1) 其中每1个区间的长度都是前一区间长度的一半，最后1个区间的长度为： bkakba (3.2) k2如果取最后1个区间[ak,bk]的中点： xkbkak (3.3) 2作为f(x)0根的近似值，则有误差估计式： x*xkbkakbak1 (3.4) 22对于所给精度，若取k使得 则有，

x*xk (3.6)

ba (3.5) 2k13.2 程序与实例

用二分法求解方程f(x)0在有根区间[a,b]内的一个根，其中f(x)在[a,b]只有一个根的情形。

5

***************************************************************************

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*yb>0,

disp('(a,b)不是有根区间'); return end

max1 = 1 + round((log(b-a) - log(delta))/log(2)); for k = 1:max1 c = (a + b)/2; yc = feval('f',c); if yc == 0

a=c;

b=c; return

elseif yb*yc>0 b = c; yb = yc; else a = c; ya = yc; end

if (b-a) < delta,return,end

end k;

c = (a + b)/2; err = abs(b-a); yc = feval('f',c);

*****************************************************************************

6

例3.1 用上述程序求方程x3x10在区间[1.0,1.5]内的一个近似解，要求准确到小数点后第2位。

解：先用m文件先定义一个名为f1031.m的函数文件。 function y = f1031(x)

y = x^3 - x - 1;

建立一个主程序prog1031.m

clc

clear all

bisect('f1031', 1.0, 1.5, 0.0005) 然后在MATLAB命令窗口运行上述主程序，即：

>> prog1031 计算结果如下。 ans =

1.3247

4 牛顿法 4.1 牛顿法原理

从前面迭代法，我们知道，迭代函数(x)构造的好坏，不仅影响收敛速度，而且迭代格式有可能发散。怎样选择一个迭代函数才能保证迭代序列一定收敛呢？

构代迭代函数的一条重要途径是用近似方程来代替原方程去求根。因此，如果能将非线性方程(1.1)用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。牛顿(Newton)法就是一种将非线性方程线化的一种方法。

设xk是方程(1.1)的一个近似根，把如果f(x)在xk处作一阶Taylor展开，即： f(x)f(xk)f'(xk)(xxk) (4.1) 于是我们得到如下近似方程：

f(xk)f'(xk)(xxk)0 (4.2) 设f'(xk)0，则方程(10.2.1)的解为：

7

xk xf(xk) (4.3) f'(xk)取~x作为原方程(1.1)的新近似根xk1，即令： xk1xkf(xk)， k0,1,2, (4.4) f'(xk)上式称为牛顿迭代格式。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法。 牛顿法具有明显的几何意义。方程：

yf(xk)f'(xk)(xxk) (10.4.5) 是曲线yf(x)上点(xk,f(xk))处的切线方程。迭代格式(4.4)就是用切线式(4.5)的零点来代替曲线(1.1)的零点。正因为如此，牛顿法也称为切线法。

牛顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。一般来说，牛顿法对初值x0的要求较高，初值足够靠近x*时才能保证收敛。若要保证初值在较大范围内收敛，则需对f(x)加一些条件。如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛顿法作一些改时，即可以采用下面的迭代格式：

xk1xkf(xk)， k0,1,2, (4.6) f'(xk)式中，01，称为下山因子。因此，用这种方法求方程的根，也称为牛顿下山法。 牛顿法对单根收敛速度快，但每迭代一次，除需计算f(xk)之外，还要计算f'(xk)的值。如果f(x)比较复杂，计算f'(xk)的工作量就可能比较大。为了避免计算导数值，我们可用差商来代替导数。通常用如下几种方法： (1) 割线法。如果用

f(xk)f(xk1)

xkxk1代替f'(xk)，则得到割线法的迭代格式为： xk1xkxkxk1f(xk) (4.7)

f(xk)f(xk1) (2) 拟牛顿法。如果用

8

f(xk)f(xkf(xk1))

f(xk)代替f'(xk)，则得到拟牛顿法的迭代格式为：

f2(xk) xk1xk (4.8)

f(xk)f(xkf(xk1)) (3) Steffenson法。如果用

f(xkf(xk))f(xk)

f(xk)代替f'(xk)，则得到拟牛顿法的迭代格式为：

f2(xk) xk1xk (4.9)

f(xkf(xk))f(xk)

4.2 程序与实例

1. 牛顿法的程序

给定初值p0，用牛顿法格式pkpk1

**************************************************************************

function [p1,err,k,y] = newton(f1041,df1041,p0,delta,max1) % f1041是非线性函数。 % df1041是f1041的微商。 % p0是初始值。

% delta是给定允许误差。 % max1是迭代的最大次数。

% p1是牛顿法求得的方程的近似解。 % err是p0的误差估计。 % k是迭代次数。 % y = f(p1)

p0, feval('f1041',p0)

9

f(pk1)，k1,2,，求解非线性方程f(x)0。

f'(pk1)for k = 1:max1

p1 = p0 - feval('f1041', p0)/feval('df1041', p0); err = abs(p1-p0); p0 = p1;

p1, err, k, y = feval('f1041', p1) if (err < delta) | (y == 0),

break, end

p1, err, k, y = feval('f1041', p1) end

****************************************************************************

例4.1 用上述程序求方程x33x20的一个近似解，给定初值为1.2，误差界为106。

解：先用m文件先定义二个名为f1041.m和df1041.m的函数文件。 function y = f1041(x)

y = x^3 – 3*x + 2;

function y=df1041(x)

y=3*x^2-3;

建立一个主程序prog1041.m

clear

newton('f1041','df1041',1.2, 10^(-6), 18) 然后在MATLAB命令窗口运行上述主程序，即：

>> prog1041 计算结果如下。

p0 = 1.2000 ans = 0.1280 p1 = 1.1030 err = 0.0970 k = 1 y = 0.0329 p1 = 1.1030

10

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

11

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

12

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

13

y = 8.1170e-009 p1 = 1.0000 err = 2.6008e-005 k = 13 y = 2.0293e-009 p1 = 1.0000 err = 2.6008e-005 k = 13 y = 2.0293e-009 p1 = 1.0000 err = 1.3004e-005 k = 14 y = 5.0732e-010 p1 = 1.0000 err = 1.3004e-005 k = 14 y = 5.0732e-010 p1 = 1.0000 err = 6.5020e-006 k = 15 y = 1.2683e-010 p1 = 1.0000 err = 6.5020e-006 k = 15 y = 1.2683e-010 p1 = 1.0000 err = 3.2510e-006 k = 16 y = 3.1708e-011 p1 = 1.0000

14

err = 3.2510e-006 k = 16 y = 3.1708e-011 p1 = 1.0000 err = 1.6255e-006 k = 17 y = 7.9270e-012 p1 = 1.0000 err = 1.6255e-006 k = 17 y = 7.9270e-012 p1 = 1.0000 err = 8.1277e-007 k = 18 y = 1.9817e-012 ans = 1.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(p1)

p0, p1, feval('f1042', p0), feval('f1042',p1), k = 0; for k = 1:max1

15

p2 = p1 - feval('f1042',p1)*(p1-p0)/(feval('f1042',p1) - feval('f1042',p0)); err = abs(p2-p1); p0 = p1; p1 = p2;

p1, err, k, y = feval('f1042', p1)； if (err < delta) | (y == 0),

break,

end

end

******************************************************************************

例4.2 用上述程序求方程x3x20的一个近似解，给定初值为-1.5，-1.52, 误差界为

106。

解：先用m文件先定义一个名为f1042.m的函数文件。 function y = f1042(x)

y = x^3 – x + 2; 建立一个主程序prog1042.m

clc

clear

secant('f1042',-1.5,-1.52,10^(-6),11) 然后在MATLAB命令窗口运行上述主程序，即：

>> prog1042 计算结果如下。

p0 = -1.5000 p1 = -1.5200 ans = 0.1250 ans = 0.0082 p1 = -1.5214 err = 0.0014 k = 1 p1 = -1.5214

16

err = 2.2961e-005 k = 2 p1 = -1.5214 err = 2.4318e-008 k = 3 ans = -1.5214

这就表明，经过3次迭代得到满足精度要求的的近似解-1.5214。

17