matlab非线性方程求解

1 非线性方程的解法

1引 言

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

0)(xf (1.1)

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

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

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

0)(*)(hxfxfkk (1.2)

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

2 迭代法

2.1 迭代法的一般概念

迭代法是数值计算中一类典型方法，不仅用于方程求根，而且用于方程组求解，矩阵求特征值等方面。迭代法的基本思想是一种逐次逼近的方法。首先取一个精糙的近似值，然后 2 用同一个递推公式，反复校正这个初值，直到满足预先给定的精度要求为止。

对于迭代法，一般需要讨论的基本问题是：迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。这里，主要看看解方程迭代式的构造。

对方程(1.1)，在区间],[ba内，可改写成为：

)(xx (2.1)

取],[0bax，用递推公式：

)(1kkxx, ,2,1,0k (2.2)

可得到序列：

0210}{,,,,kkkxxxxx (2.3)

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

)~(~xx

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

xx~*

即，

*limxxkk

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

2.2 程序和实例

下面是基于MATLAB的迭代法程序，用迭代格式)(1nnxgp，求解方程)(xgx，其中初始值为0p。

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

function[p,k,err,P]=fixpt(f1021,p0,tol,max1)

% f1021是给定的迭代函数。

% p0是给定的初始值。 3 % 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(err

break;

end

if k == max1

disp('maximum number of iterations exceeded');

end

end

P=P;

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

例2.1 用上述程序求方程0sin2xx的一个近似解，给定初始值5.00x，误差界为510。

解：先用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

err = 0.4589 4 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 二分法原理

二分法是方程求解最直观、最简单的方法。二分法以连续函数的介值定理为基础的。由介值定理知道，若函数)(xf区间],[ba上连续，且0)(*)(bfaf，即)(af和)(bf负号相反， 5 则)(xf在],[ba内一定有实根。二分法的基本思想是：用对分区间的方法根据分点处函数)(xf的符号逐步将有限区间缩小，使在足够小的区间内，方程有且仅有一根。下面简述其基本步骤。

首先记bbaa00,。用中点2000bax将区间],[00ba等分成2个小区间：],[00xa和],[00bx。然后分析可能存在的三种情况：

如果0)(*)(0xfaf，则0x是零点，也就是方程的根。

如果0)(*)(0xfaf，则区间],[00xa内存在零点。

如果0)(*)(0bfxf，则区间],[00bx内存在零点。

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

],[],[],[],[221100kkbabababa (3.1)

其中每1个区间的长度都是前一区间长度的一半，最后1个区间的长度为：

kkkabab2 (3.2)

如果取最后1个区间],[kkba的中点：

2kkkabx (3.3)

作为0)(xf根的近似值，则有误差估计式：

1*22kkkkababxx (3.4)

对于所给精度，若取k使得

12kab (3.5)

则有，

kxx* (3.6)

3.2 程序与实例

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

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

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);

***************************************************************************** 7