matlab非线性方程求解
非线性方程的解法
1引 言
数学物理中的许多问题归结为解函数方程的问题,即, 这里,数
f(x)
f(x)0
(1.1)
f(x)0
f(x)
可以是代数多项式,也可以是超越函数。若有数x*为方程的根,或称函
的零点。
f(x)
设函数在[a,b]内连续,且f(a)f(b)0。根据连续函数的性质知道,方程f(x)0
f(x)0
在
区间[a,b]内至少有一个实根;我们又知道,方程的根,除了极少简单方程的根可以
用解析式表达外,一般方程的根很难用一个式子表达。即使能表示成解析式的,往往也很复杂,不便计算。所以,具体求根时,一般先寻求根的某一个初始近似值,然后再将初始近似值逐步加工成满足精度要求为止。
如何寻求根的初始值呢?简单述之,为了明确起见,不妨设的单根,且
f(a)0,f(b)0。我们从左端出点x0a
f(x)
在区间[a,b]内有一个实
出发,按某一预定的步长h一步一步地
h
向右跨,每跨一步进行一次根的“搜索”,即检查每一步的起点xk和xk1(即,xk是否同号。若有:
f(xk)*f(xkh)0
h)
)的函数值
(1.2) 内,这时可取xk或xk
h
那么所求的根必在(xk,xk
作为根的初始近似值。这种方法通常称
为“定步长搜索法”。另外,还是图解法、近似方程法和解析法。
2 迭代法
2.1 迭代法的一般概念
迭代法是数值计算中一类典型方法,不仅用于方程求根,而且用于方程组求解,矩阵求特征值等方面。迭代法的基本思想是一种逐次逼近的方法。首先取一个精糙的近似值,然后
用同一个递推公式,反复校正这个初值,直到满足预先给定的精度要求为止。
对于迭代法,一般需要讨论的基本问题是:迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。这里,主要看看解方程迭代式的构造。 对方程(1.1),在区间[a,b]内,可改写成为:
x(x)
(2.1)
取x0[a,b],用递推公式:
xk1(xk)
,
k0,1,2,
(2.2)
可得到序列: 当k
x0,x1,x2,xk,{xk}k0
(2.3)
x,且(x)在~x附近连续,则在式(2.2)两边极限,得, 时,序列{xk}有极限~k0
~x(~x)
x为方程(2.1)的根。由于方式(1.1)和方程(2.1)等价,所以, 即,~
即,
*
x~x
limxkx
k
*
式(2.2)称为迭代式,也称为迭代公式;(x)可称为迭代函数。称求得的序列{xk}为迭代序k0列。
2.2 程序和实例
下面是基于MATLAB的迭代法程序,用迭代格式pn1
初始值为p0。
**************************************************************************
function[p,k,err,P]=fixpt(f1021,p0,tol,max1) % f1021是给定的迭代函数。 % p0是给定的初始值。
g(xn),求解方程xg(x)
,其中
% 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
end
if k == max1
disp('maximum number of iterations exceeded');
end end P=P;
**************************************************************************** 例2.1 用上述程序求方程sin
10
5
xx0
2
的一个近似解,给定初始值x0
0.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 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)负号相反,
则
f(x)
在[a,b]内一定有实根。二分法的基本思想是:用对分区间的方法根据分点处函数
f(x)
的符号逐步将有限区间缩小,使在足够小的区间内,方程有且仅有一根。下面简述其基本步骤。 首先记a0
a,b0b
。用中点x0
a0b0
2
将区间[a0,b0]等分成2个小区间:[a0,x0]和
[x0,b0]。然后分析可能存在的三种情况:
如果 如果
如果
f(a)*f(x0)0f(a)*f(x0)0f(x0)*f(b)0
,则x0是零点,也就是方程的根。 ,则区间[a0,x0]内存在零点。 ,则区间[x0,b0]内存在零点。
对有根的新区间施行同样的操作,于是得到一系列有空的区间:
[a0,b0][a1,b1][a2,b2][ak,bk]
(3.1)
其中每1个区间的长度都是前一区间长度的一半,最后1个区间的长度为:
bkak
ba2
k
(3.2)
如果取最后1个区间[ak,bk]的中点: 作为
xk
bkak
2
(3.3)
f(x)0
*
根的近似值,则有误差估计式:
bkak
2
ba2
k1
xxk
(3.4)
对于所给精度,若取k使得 则有,
xxk
*
ba2
k1
(3.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*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)
end k;
c = (a + b)/2; err = abs(b-a); yc = feval('f',c);
*****************************************************************************
例3.1 用上述程序求方程x3点后第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
x10
在区间[1.0,1.5]内的一个近似解,要求准确到小数
4 牛顿法 4.1 牛顿法原理
从前面迭代法,我们知道,迭代函数(x)构造的好坏,不仅影响收敛速度,而且迭代格式有可能发散。怎样选择一个迭代函数才能保证迭代序列一定收敛呢?
构代迭代函数的一条重要途径是用近似方程来代替原方程去求根。因此,如果能将非线性方程(1.1)用线性方程去代替,那么,求近似根问题就很容易解决,而且十分方便。牛顿(Newton)法就是一种将非线性方程线化的一种方法。 设xk是方程(1.1)的一个近似根,把如果
f(x)f(xk)f'(xk)(xxk)
f(x)
在xk处作一阶Taylor展开,即:
(4.1)
于是我们得到如下近似方程: 设
f(xk)f'(xk)(xxk)0
(4.2)
f'(xk)0,则方程(10.2.1)的解为:
xkx
f(xk)f'(xk)
(4.3)
x作为原方程(1.1)的新近似根xk1,即令: 取~
xk1xk
f(xk)f'(xk)
,
k0,1,2,
(4.4)
上式称为牛顿迭代格式。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法,简称牛顿法。 牛顿法具有明显的几何意义。方程: 是曲线y
yf(xk)f'(xk)(xxk)f(x)上点(xk,f(xk))
(10.4.5)
处的切线方程。迭代格式(4.4)就是用切线式(4.5)的零点来代替
曲线(1.1)的零点。正因为如此,牛顿法也称为切线法。
牛顿迭代法对单根至少是二阶局部收敛的,而对于重根是一阶局部收敛的。一般来说,牛顿法对初值x0的要求较高,初值足够靠近x*时才能保证收敛。若要保证初值在较大范围内收敛,则需对
f(x)
加一些条件。如果所加的条件不满足,而导致牛顿法不收敛时,则需对牛
顿法作一些改时,即可以采用下面的迭代格式:
xk1xk
f(xk)f'(xk)
,
k0,1,2,
(4.6)
式中,01,称为下山因子。因此,用这种方法求方程的根,也称为牛顿下山法。
牛顿法对单根收敛速度快,但每迭代一次,除需计算如果
f(x)
f(xk)之外,还要计算f'(xk)的值。
比较复杂,计算
f'(xk)的工作量就可能比较大。为了避免计算导数值,我们可用差
商来代替导数。通常用如下几种方法: (1) 割线法。如果用 代替
f(xk)f(xk1)
xkxk1
f'(xk),则得到割线法的迭代格式为:
xkxk1f(xk)f(xk1)
xk1xkf(xk)
(4.7)
(2) 拟牛顿法。如果用
代替
f(xk)f(xkf(xk1))
f(xk)
f'(xk),则得到拟牛顿法的迭代格式为:
f(xk)
f(xk)f(xkf(xk1))
2
xk1xk
(4.8)
(3) Steffenson法。如果用 代替
f(xkf(xk))f(xk)
f(xk)
f'(xk),则得到拟牛顿法的迭代格式为:
f(xk)
f(xkf(xk))f(xk)
2
xk1xk
(4.9)
4.2 程序与实例
1. 牛顿法的程序
给定初值p0,用牛顿法格式pk
**************************************************************************
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)
pk1
f(pk1)f'(pk1)
,k
1,2,
,求解非线性方程f(x)0。
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
break, end
p1, err, k, y = feval('f1041', p1) end
****************************************************************************
例4.1 用上述程序求方程x3
3x20的一个近似解,给定初值为
1.2,误差界为106。
解:先用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
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
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
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
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
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
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
break,
end
end
******************************************************************************
例4.2 用上述程序求方程x3
106x20的一个近似解,给定初值为-1.5,-1.52, 误差界为。
解:先用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
err = 2.2961e-005
k = 2
p1 = -1.5214
err = 2.4318e-008
k = 3
ans = -1.5214
这就表明,经过3次迭代得到满足精度要求的的近似解-1.5214。