matlab非线性方程求解

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

2、，每跨一步进行一次根的“搜索”，即检查每一步的起点和(即，)的函数值是否同号。若有： (1.2)那么所求的根必在内，这时可取或作为根的初始近似值。这种方法通常称为“定步长搜索法”。另外，还是图解法、近似方程法和解析法。2 迭代法2.1 迭代法的一般概念迭代法是数值计算中一类典型方法，不仅用于方程求根，而且用于方程组求解，矩阵求特征值等方面。迭代法的基本思想是一种逐次逼近的方法。首先取一个精糙的近似值，然后用同一个递推公式，反复校正这个初值，直到满足预先给定的精度要求为止。 对于迭代法，一般需要讨论的基本问题是：迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。这里，主要看看解方程迭代式的构

3、造。 对方程(1.1)，在区间内，可改写成为： (2.1)取，用递推公式： , (2.2)可得到序列： (2.3)当时，序列有极限，且在附近连续，则在式(2.2)两边极限，得， 即，为方程(2.1)的根。由于方式(1.1)和方程(2.1)等价，所以， 即， 式(2.2)称为迭代式，也称为迭代公式；可称为迭代函数。称求得的序列为迭代序列。 2.2 程序和实例下面是基于MATLAB的迭代法程序，用迭代格式，求解方程，其中初始值为。 *functionp,k,err,P=fixpt(f1021,p0,tol,max1)% f1021是给定的迭代函数。% p0是给定的初始值。% tol是给定的误差界。

4、% max1是所允许的最大迭代次数。% k是所进行的迭代次数加1。% p是不动点的近似值。% err是误差。% P = p1,p2,pnP(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 prog1021计算结果如下。k = 2err = 0.4589k = 3err = 0.1052k = 4err = 0.0292k = 5err = 0.0078k = 6err = 0.0021k = 7err = 5.7408e-004k = 8err = 1.

5、5525e-004k = 9err = 4.1975e-005k = 10err = 1.1350e-005k = 11err = 3.0688e-006P = 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.8767ans = 0.87673 二分法3.1 二分法原理 二分法是方程求解最直观、最简单的方法。二分法以连续函数的介值定理为基础的。由介值定理知道，若函数区间上连续，且，即和负号相反，则在内一定有实根。二分法

6、的基本思想是：用对分区间的方法根据分点处函数的符号逐步将有限区间缩小，使在足够小的区间内，方程有且仅有一根。下面简述其基本步骤。 首先记。用中点将区间等分成2个小区间：和。然后分析可能存在的三种情况： 如果，则是零点，也就是方程的根。 如果，则区间内存在零点。如果，则区间内存在零点。 对有根的新区间施行同样的操作，于是得到一系列有空的区间： (3.1)其中每1个区间的长度都是前一区间长度的一半，最后1个区间的长度为： (3.2)如果取最后1个区间的中点： (3.3)作为根的近似值，则有误差估计式： (3.4)对于所给精度，若取使得 (3.5)则有， (3.6)3.2 程序与实例用二分法求解方程

7、在有根区间内的一个根，其中在只有一个根的情形。*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； returnendif ya*yb0, disp(a,b)不是有根区间); returnendmax1 = 1 + round(log(b-a) - log(delta)/log(2);f

8、or k = 1:max1 c = (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; end if (b-a) prog1031计算结果如下。 ans = 1.32474 牛顿法4.1 牛顿法原理从前面迭代法，我们知道，迭代函数构造的好坏，不仅影响收敛速度，而且迭代格式有可能发散。怎样选择一个迭代函数才能保证迭代序列一定收敛呢？ 构代迭代函数的一条重要途径是用近似方程来代替原方程去求根。因此，如果能将非线性方程(1.1)用线性方程去代

9、替，那么，求近似根问题就很容易解决，而且十分方便。牛顿(Newton)法就是一种将非线性方程线化的一种方法。 设是方程(1.1)的一个近似根，把如果在处作一阶Taylor展开，即： (4.1)于是我们得到如下近似方程： (4.2)设，则方程(10.2.1)的解为： (4.3)取作为原方程(1.1)的新近似根，即令： ， (4.4)上式称为牛顿迭代格式。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法。 牛顿法具有明显的几何意义。方程： (10.4.5)是曲线上点处的切线方程。迭代格式(4.4)就是用切线式(4.5)的零点来代替曲线(1.1)的零点。正因为如此，牛顿法也称为切线法。 牛

10、顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。一般来说，牛顿法对初值的要求较高，初值足够靠近时才能保证收敛。若要保证初值在较大范围内收敛，则需对加一些条件。如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛顿法作一些改时，即可以采用下面的迭代格式：， (4.6)式中，称为下山因子。因此，用这种方法求方程的根，也称为牛顿下山法。 牛顿法对单根收敛速度快，但每迭代一次，除需计算之外，还要计算的值。如果比较复杂，计算的工作量就可能比较大。为了避免计算导数值，我们可用差商来代替导数。通常用如下几种方法： (1) 割线法。如果用 代替，则得到割线法的迭代格式为： (4.7) (2)

11、拟牛顿法。如果用 代替，则得到拟牛顿法的迭代格式为： (4.8) (3) Steffenson法。如果用 代替，则得到拟牛顿法的迭代格式为： (4.9)4.2 程序与实例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, f

12、eval(f1041,p0)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 prog1041计算结果如下。p0 = 1.2000ans = 0.1280p1 = 1.1030err = 0.0970k = 1y = 0.0329p1 = 1.1030err = 0.0970k = 1y = 0.0329p1 = 1.0524err = 0.0507k = 2y = 0.0084p1 =

13、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.0066err = 0.0066k = 5y = 1.3270e-004p1 = 1.0033err = 0.

14、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 = 9y = 5.1943e-007p1

15、= 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.000

16、1err = 5.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 =

17、 6.5020e-006k = 15y = 1.2683e-010p1 = 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-012ans = 1.0000 这说明，经过18次

18、迭代得到满足精度要求的值。2. 割线法的MATLAB实现 *function p1, err, k, y = secant(f1042, p0, p1, delta, max1)% f1042是给定的非线性函数。% p0,p1为初始值。% delta为给定误差界。% max1是迭代次数的上限。% p1为所求得的方程的近似解。% err为p1-p0的绝对值。% k为所需的迭代次数。% yf(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);