数学实验之2.非线性方程(组)求解(1)

1、重庆大学信息与计算科学系报告人: 何仁斌 非线性方程(组)求解(1) 引言及引例引言及引例非线性方程的常用求解方法非线性方程的常用求解方法使用使用MATLAB的命令进行方程求解的命令进行方程求解范例范例:波音公司最优定价策略问题波音公司最优定价策略问题实实 验目的与内容验目的与内容 求高次代数方程的根是一个基本的、古老的数学问题。0011102nnnnaxaxaxacbxax高阶方程能否求解？定义：一些由实际问题列出的方程中常常包含三角函数、指数函数等，它与n(2)代数方程一起统称为非线性方程非线性方程(组组)，记作0)(xf引言例：exsin(x)+5x3=0一元二次方程一元二次方程高阶多项

2、式方程高阶多项式方程 计算机技术和数学软件的飞速发展，为一些古老数学问题的解决提供了极大的方便。如利用MATLAB的图形功能就能帮助我们判断方程有没有根，确定根的近似位置，例如：0e50305813261282345xxxxxxx-1-0.500.511.522.53-60-40-20020406080100120140 x8 x5-12 x4-26 x3-13 x2+58 x+308x5-12x4-26x3-13x2+58x+30 00.511.522.53-4-3-2-10123x5 x-exp(x)5x-exp(x) -0.51.52.50.252.5 对于一般的非线性方程，怎样得到比较

3、精确的近似解，正是该试验要研究或解决的问题。 下面给出一个实际问题，它们都归结为求非线性方程的近似解。放射性废物的处理问题 一段时间一段时间, 美国原子能委员会是按以下方式处理浓缩放射美国原子能委员会是按以下方式处理浓缩放射性废物的性废物的. 他们将废物装入密封性能很好的圆桶中他们将废物装入密封性能很好的圆桶中, 然后扔到水然后扔到水深深300英尺的海里英尺的海里. 这种做法是否会造成放射性污染这种做法是否会造成放射性污染, 很自然地很自然地引起了生态学家及社会各界的关注引起了生态学家及社会各界的关注. 原子能委员会一再保证原子能委员会一再保证, 圆圆桶非常坚固桶非常坚固, 决不会破漏决不会破

4、漏, 这种做法是绝对安全的这种做法是绝对安全的. 然而一些工然而一些工程师们却对此表示怀疑程师们却对此表示怀疑, 他们认为圆桶在海底相撞时有可能发他们认为圆桶在海底相撞时有可能发生破裂生破裂. 由此双方展开了一场笔墨官司由此双方展开了一场笔墨官司. 究竟谁的意见正确呢究竟谁的意见正确呢? 只能让事实说话了只能让事实说话了!引例 问题的关键在于圆桶到底能承受多大速度问题的关键在于圆桶到底能承受多大速度的碰撞的碰撞? 圆桶和海底碰撞时的速度有多大圆桶和海底碰撞时的速度有多大? 工程师们进行了大量破坏性的实验工程师们进行了大量破坏性的实验, 发现发现圆桶在直线速度为圆桶在直线速度为40 ft/s 的

5、冲撞下会发生破裂的冲撞下会发生破裂, 剩下的问题就是计算圆桶沉入剩下的问题就是计算圆桶沉入300 ft 深的海底深的海底时时, 其末速度究竟有多大其末速度究竟有多大? 问题分析问题分析 1. 使用使用55加仑的圆桶加仑的圆桶; ( 1加仑加仑 = 3.7854升升 ) 2. 装满放射性废物时的圆桶重量为装满放射性废物时的圆桶重量为 W = 527.436磅磅 (1 磅磅 = 0.4526公斤公斤 ) 3. 在海水中圆桶受到的浮力在海水中圆桶受到的浮力 B = 470.327磅磅 4. 圆桶下沉时受到海水的阻力圆桶下沉时受到海水的阻力 D = C v C 为常数为常数, 经测算得经测算得: C

6、= 0.08. 5. 建立坐标系建立坐标系, 取垂直向下为坐标方向取垂直向下为坐标方向 y , 海平面为坐标原点海平面为坐标原点.y0 假设假设:根据牛顿第二定律根据牛顿第二定律, 圆桶下沉时应满足微分方程圆桶下沉时应满足微分方程:)(/(86.713)(lim)1 ()(0)0(:dd,dd22时当速度：容易计算出圆桶的极限其解：初值条件其中阻力）浮力（重力tsftCBWtveCBWtvvBCvWdtdvmvtyCvDgWmDBWtymtWCg327.470B 为了求出圆桶与海底的碰撞速度为了求出圆桶与海底的碰撞速度v(t), 需要求出需要求出圆桶下沉到海底圆桶下沉到海底300英尺时的时间英

7、尺时的时间 t, 再计算再计算v(t)，要要做到这一点是十分困难的做到这一点是十分困难的. 若将速度若将速度v 看成是海水深看成是海水深度度y 的函数的函数, 即即vyvtyyvtvtytyvtvdddddddddd)()(22由复合函数的求导法知由复合函数的求导法知WgyBWCvBWCBWCvyvWgyvCvBWvCvBWyvmvln0)0(, 0)0(dddd2积分，得：初值条件：或微分方程变为：如何从上述方程中解出函数如何从上述方程中解出函数 v = v(y)？只能借助？只能借助数值方法求出数值方法求出v(300)=45.1ft/s，显然大于，显然大于40ft/s。结论：工程师们的怀疑有

8、道理。结论：工程师们的怀疑有道理。数学建模的一般步骤数学建模的一般步骤模型准备模型准备模型假设模型假设模型构成模型构成模型检验模型检验模型分析模型分析模型求解模型求解模型应用模型应用引言及引例引言及引例非线性方程的常用求解方法非线性方程的常用求解方法迭代迭代使用使用MATLAB的命令进行方程求解的命令进行方程求解范例范例:波音公司最优定价策略问题波音公司最优定价策略问题实实 验目的与内容验目的与内容求解方法求解方法图形放大法图形放大法 简单迭代算法简单迭代算法数值求解法数值求解法非线性方程非线性方程三种形式的图形描述三种形式的图形描述 f(x) = 0软件求解软件求解 方程 f(x)=0 1）

9、建立坐标系，画曲线f(x)； 2）观察曲线f(x)与x轴相交的交点； 3）将其中一个交点进行局部放大； 4）该交点的横坐标值就是方程的根。 1图形放大法例例: : 求方程 x5 2x2 + 4 = 0 的一个根.程序（tuxfd.m）x=-6:0.01:6;y=x.5+2.*x.2+4;y1=zeros(size(x);plot(x,y,x,y1)或ezplot(f(x),a,b)-6-4-20246-8000-6000-4000-200002000400060008000将区间将区间-6，6缩小至缩小至-2，2； 显然该方程有5个根，欲寻找其中一个实根，并且达到一定的精度。逐次缩小区间，观察

10、一个根在-1.55-1.5之间。-202-50050-2-1.5-1-20-10010-2-1.5-20-10010-1.6-1.5-1.4-2-101MATLAB(tuxfd1.m)2. 2. 简单的迭代算法简单的迭代算法以例说明方法。例：以例说明方法。例： 3xex = 0 1）该方程有多少个根？如何判）该方程有多少个根？如何判 断？断？2）如何进行迭代求解？）如何进行迭代求解？将方程变形为将方程变形为： x = ex/3 0 1/3 1/3 0.4652 0.4652 0.5308序号0123456左边00.333 0.465 0.531 0.567 0.588 0.599右边 0.33

11、3 0.465 0.531 0.567 0.588 0.599 0.607序号78910左边 0.607 0.612 0.615 0.616右边 0.612 0.615 0.616 0.617方程3xex=0的迭代求解表： x = ex/300.5100.20.40.60.810.616y=xy =ex/3迭代过程如下图所示： 方程： f (x) = 0 经过简单变形：x = j (x)。 或 j (x) = f (x)+x。 x被称为不动点，该方法称为不动点法不动点法。迭代过程如下： xn+1 =j (xn)，n =0，1， x0 定义为迭代初值。迭代算法 1.1.表达式表达式x = j (

12、x)是否唯一是否唯一？2.2.迭代产生的序列是否一定会收敛？迭代产生的序列是否一定会收敛？3. 3. 迭代与初始值迭代与初始值x0是否有关？是否有关？解：解： 第一步第一步 构造迭代函数： x=j (x)(111)(1)(132232123xxxxxxxxxxxxjjj例：例：用迭代方法求解方程 x3 x2 x1 0。第二步第二步 迭代设定初值x0=1， xn+1 =j (xn)，n =0，1， 用MATLAB编程（died2.m)x=1;y=1;z=1;(初始点）for k=1:20 x=x3-x2-1; % j1 (x) y=(y2+y+1)(1/3); % j2 (y) z=1+1/z+

13、1/z2; % j3 (z)endxyz计算结果序号j2(x)j3(x)序号j2(x)j3(x)11.4422 3.000081.8175 1.813621.6537 1.444491.8385 1.855431.7532 2.1716101.8389 1.829441.7995 1.6725111.8391 1.845451.8209 1.9554121.8392 1.835561.8308 1.7730131.8392 1.841671.8354 1.8822j1(x)的迭代是失败的（迭代不收敛 ）。精确解：x=1.8393 迭代函数j2(x)和j3(x)的选取是成功的。近似解为 x=1.

14、8393。 并且选取函数j2(x)、j3(x)其收敛速度不一致，前者的速度快些！结论1、当遇到迭代不收敛时有什么解决办法？2、如何提高收敛速度？ 对于给定的方程f(x) = 0, 有多种方式将它改写成等价的形式 x = j(x)。但重要的是如何改写使得序如何改写使得序列收敛？列收敛？ 若 x= j (x) 迭代不收敛，则不直接使用j(x)迭代，而用由j (x)与x的加权平均，h(x) =l j (x) +(1l)x 进行迭代，其中l为参数。显然 x = h (x) x = j (x) 加速迭代收敛方法加速迭代收敛方法xn+1 = h (xn) xn+1 = j (xn) 关键是如何确定函数关键

15、是如何确定函数h(x) 中的参数中的参数l ？理论证明：在满足|h(x)|1的条件下， 令 h(a)=0，解出)(11ajl)(11nxjl用用xn替换替换a，得，得从而迭代过程：从而迭代过程：)(1)()()1 ()()(1nnnnnnnnxxxxxxxhxjjjllj即 l j (a) +(1l)=0 例如：当j1(x)= x3x21时，进行改进得：12312)23(1)23() 1()(2232223xxxxxxxxxxxxh12312)(2231nnnnnnxxxxxhx经实验发现，它比经实验发现，它比j2(x) ，j3(x)的收敛速度要快！的收敛速度要快！函数的迭代函数的迭代 二次函

16、数二次函数为参数其中axaxxf),1 ()(考虑a 的取值范围0, 4。 给定实数域上光滑的实值函数f(x)以及初值x0, 构造迭代格式：xk+1 = f(xk)。产生的数列xk，k = 1,2, 称为函数f(x)的一个迭代序列。 函数的迭代是数学研究中的一个非常重要的思想工具。哪怕是对一个相当简单的函数进行迭代，都可以产生异常复杂的行为，并由此衍生了一些崭新的学科分支，如分形与混沌。例如函数的迭代函数的迭代 4 , 0),1 (1axaxxkkk注意：不同的参数a，产生的序列xn可能收敛或可能发散。观察a取什么值时收敛（或发散）？显示的图形方式有三种：1）线性连接图线性连接图 用点(n,

17、xn)描述迭代点，并用直 线连接这些点所形成的折线图（横坐标表示n，纵坐标表示 xn）。图形表示如下：我们可以用图形的方式来反映迭代的过程。如：线性连接图线性连接图震荡发散情形收敛情形0102000.10.20.30.4a=0.5,x0=0.5010200.20.40.60.81a=2.5,x0=0.5010200.20.40.60.81a=3.1,x0=0.5010200.20.40.60.81a=3.5,x0=0.5a=0.5;x1=;x1(1)=0.5;for i=2:20 x1(i)=a*x1(i-1)*(1-x1(i-1);endn=1:20;subplot(2,2,1),plot(

18、n,x1),title(a=0.5,x0=0.5)MATLAB程序程序Died4.m2）蛛网图）蛛网图 在直角坐标系中，首先画出直线y= x和曲线y= f(x)，其中f(x)为迭代曲线。迭代过程如下图。 从直线y= x上点An(xn,xn)到曲线y= f(x)上点Bn(xn,xn+1)；（垂直） 从Bn(xn,xn+1)到An(xn+1,xn+1); (水平） 重复。a=2.9, x0=0.2 -0.200.20.40.60.811.2-0.200.20.40.60.81AnBnAny = xy = f(x)died8.mx1=;a=2.9;x=-0.2:0.01:1.2;y=a.*x.*(1

19、-x);plot(x,y),hold on % 画二次函数曲线ezplot(x,-0.2,1.2)% 画直线x1(1)=0.2; % 初始点for i=2:50 x1(i)=a*x1(i-1)*(1-x1(i-1); plot(x1(i-1),x1(i-1),x1(i-1),x1(i); plot(x1(i-1),x1(i),x1(i),x1(i);End % 画折线x1;died8.m00.51-0.200.20.40.60.811.2xxa=2.5, x0=0.1 00.51-0.200.20.40.60.811.2xxa=1.5, x0=0.1 00.51-0.200.20.40.60.

20、811.2xxa=3.4, x0=0.2 00.51-0.200.20.40.60.811.2xxa=3.7, x0=0.2, n=80 蛛网图蛛网图发散情形收敛情形2）分枝与混沌图）分枝与混沌图 前面两种图形中，我们观察到对不同参数参数迭代序列收敛性的变化情况。但由于参数仅取了几个离散值，导致序列收敛与发散的参数a的分界值究竟在什么位置？能否对参数连续取值观察序列收敛性的变化？axn a1a22）分枝与混沌图）分枝与混沌图 具体做法是：将区间（0，4以某个步长a离散化。对每个离散的a值做迭代，忽略前50个迭代值x1,x2,x50，而把点(a, x51), (a, x52), , (a, x1

21、00)显示在坐标平面上，最后形成的图形称为Feigenbaum图。该图主要目的是：对不同的参数a系统地观察迭代的行为。Feigenbaum图图该图形特点：观察当参数a连续变化时，函数f(x)=ax(1-x)的收敛与发散情况。a=2.9a=3.4 需要编写一个对含参变量函数 f(x,a) 进行迭代可调用的通用程序。编制一个函数文件 iter.mfunction root=iter(x,a)x=;x(1)=0.2;for i=2:100 x(i)=a*x(i-1)*(1-x(i-1);endroot=x; % 产生100个迭代序列组成的数组clf; % 清除现有的图形清除现有的图形x=;hold

22、on;for a=2:0.01:4 root=iter(x,a); plot(a.*ones(size(root(51:100),root(51:100),.)endxlabel(parameter a);ylabel(迭代序列（迭代序列（51-10051-100）); 运行程序：运行程序：died9.m混沌的特性混沌的特性 从极限分支点之后，Feigenbaum图显得很杂乱，似乎没有任何规律。实际上，对任何初始值做迭代都会得到同样的结果。这就是所谓的混沌现象。 我们以a=4为例来说明迭代序列对初值的敏感性。 观察：无论两个初值如何接近，在迭代过程中，与之产生的对应序列xn、yn将渐渐分开，称

23、迭代对初值是敏感的。这是任何混沌系统都具有的特性之一。这种特性使得混沌系统会产生似乎是随机的没有规律的现象。混沌的特性混沌的特性迭代格式： xn+1=4xn(1-xn), n = 0,1,1、取初始点x0=0.2, 记产生的序列为xn；2、取初始点x0=0.201, 记产生的序列为yn；3、作线性连接图。即 plot(n, xn-yn)；01020304050-1-0.8-0.6-0.4-0.200.20.4a=4,x0=0.5,x0=0.5101020304050-1-0.500.51a=4,x0=0.2,x0=0.201混沌的特性混沌的特性 另外，混沌不是随机的，由图形可以发现，它们具有自相