《经济数学基础》完整课件（推荐下载）：Exp06

1 大学数学实验 Experiments in Mathematics 实验6 非线性方程求解 清华大学数学科学系 什么叫方程（组） ? 方程：包含未知数（或 /和未知函数）的等式 方程组：包含未知数（或 /和未知函数）的一组等式 不包含未知函数的方程组的一般形式： F(x)=0 x= (x 1 ,x 2 ,…,x n ) T , F(x)=(f 1 (x),f 2 (x),…,f m (x)) T (一般 m=n) 满足方程（组）的未知数的取值称为方程（组）的 解，或称为 F(x)的零点。 单变量方程（一元方程）： f(x)=0, “解 ”也称为 “根 ” 非线性方程的特点 方程根的特点： ? n次代数方程有且只有 n个根 (包括复根、重根 )； ? 5次以上的代数方程无求根公式； ? 超越方程有无根，有几个根通常难以判断。 方程分类： ? 代数方程 : a 0 x n +a 1 x n-1 +…+a n =0； ? 超越方程：包含超越函数 (如 sinx, lnx)的方程； ? 非线性方程： n(≥2)次代数方程和超越方程。 实验6的基本内容 3.实际问题中非线性方程的数值解 1.非线性方程 f(x)=0 的数值解法: ? 迭代方法的基本原理; ? 牛顿法; 拟牛顿法. 2.推广到解非线性方程组 4.分岔和混沌现象 实例1 路灯照明 道路两侧分别安装路灯，在漆黑的夜晚，当两只路灯开启时， 两只路灯连线的路面上最暗的点和最亮的点在哪里？ h 2 P 2 P 1 s h 1 ? 如果 P 2 的高度可以在 3米到 9米之间变化，如何使路面上最暗 点的亮度最大？ ? 如果两只路灯的高度均可以在 3米到 9米之间变化呢？ s=20(米 ) P 1 =2, P 2 =3(千瓦 ) h 1 =5, h 2 =6(米 ) 实例1 路灯照明 建立坐标系如图，两个光源在点 Q(x,0)的照度分别为 (k是由量纲单位决定的比 例系数，不妨记 k=1)2 2 22 2 2 1 11 1 sin , sin r P kI r P kI αα == 点 Q的照度 22 2 2 2 22 1 2 1 )(, xshrxhr ?+=+= x x α 2 α 1 O h 2 P 2 r 1 P 1 s r 2 h 1 y Q 22 1 1 1 1 1 sin xh h r h + ==α 22 2 2 2 2 2 )( sin xsh h r h ?+ ==α ()() 3 22 2 22 3 22 1 11 )( )( xsh hP xh hP xC ?+ + + = 2 实例1 路灯照明 为求最暗点和最亮点，先求 C(x)的驻点 x x α 2 α 1 O h 2 P 2 r 1 P 1 s r 2 h 1 y 令 C’ (x)=0：解析解是难以求出，需数值求解 Q ()() 3 22 2 22 3 22 1 11 )( )( xsh hP xh hP xC ?+ + + = ()() 5 22 2 22 5 22 1 11 )( )( 33)(' xsh xshP xh xhP xC ?+ ? + + ?= 实例2 均相共沸混合物的组分 均相共沸混合物（ homogeneous azeotrope） 是由两种或两种 以上物质组成的液体混合物，当在某种压力下被蒸馏或局部汽 化时，在气体状态下和在液体状态下保持相同的组分（比例） 给定几种物质，如何确定它们构成均相共沸混合物时的比例？ 设该混合物由 n个可能的物质组成，物质 i所占的比例为 x i 模型建立 0,1 1 ≥= ∑ = i n i i xx 实例2 均相共沸混合物的组分 均相共沸混合物应该满足 稳定条件 ，即共沸混合物的每个组分 在气体和液体状态下具有相同的化学势能。在压强 P不大的情 况下，这个条件可以表示为： P i 是物质 i的饱和汽相压强，与温度 T有关，可以如下确定： niPP ii ,,1, L==γ ni cT b aP i i ii ,,1,ln L= + ?= 是组分 i的液相活度系数，可以根据如下表达式确定 :i γ ni qx qx qx n j n k jkk jij n j ijji ,,1,ln1ln 1 1 1 L= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= ∑ ∑ ∑ = = = γ (q ij 表示组分 i与组分 j的交互作用参数 ，可以通过实验近似得到 ) (a i , b i ， c i 是常数 ) 实例2 均相共沸混合物的组分 ii PP γ= i i ii cT b aP + ?=ln ∑ ∑ ∑ = = = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= n j n k jkk jij n j ijji qx qx qx 1 1 1 ln1lnγ 只有当物质 i参与到该共沸混合物中时才需要满足上式，故得到 0ln1ln 1 1 1 =+?? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? + + ∑ ∑ ∑ = = = Pa qx qx qx cT b i n j n k jkk jij n j ijj i i 0ln1ln 1 1 1 = ? ? ? ? ? ? ? ? ? ? ? ? +?? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? + + ∑ ∑ ∑ = = = Pa qx qx qx cT b x i n j n k jkk jij n j ijj i i i nixx i n i i ,...,2,1,0,1 1 =≥= ∑ = 解析解是难以求 出，需数值求解 ? 根的隔离:二分法 解方程 f(x)=0第一步 ——确定根的大致范围 ? 图形法：作 f(x)图形,观察 f(x)与 x轴的交点 非线性方程的基本解法 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -50 -40 -30 -20 -10 0 10 20 30 01281362 2346 =++??? xxxxx 图形法 有 4个根分别位于 x= -1.75, -0.75, 1.00， 2.40附近 Exam0600.m 若对于 ba < , 有 0)()( <? bfaf , 则在 ),( ba 内 )(xf 至少有一个零点，即 0)( =xf 至少有一个根。 二分法 的原理 二分法 的实现 不足 收敛速度较慢 0)( <af 0)( >bf 0)( 0 <xf0)( 0 >xf LL ???? ),(),(),( 11 nn bababa 011 , xbaa == bbxa == 101 , 区间每次缩小一半， n足 够大时，可确定根的范围 a b )(xf x xa b )(xf 0 x 0 x 中点 ),(: 0 bax 非线性方程的基本解法 3 例1 014)( 2 =?+= xxxf 2 1 14)( xxx ?==? , 迭代公式： 2 1 14 kk xx ?= + )1/(14)( 2 +== xxx ? ，迭代公式： )1/(14 1 += + kk xx )12/()14()( 2 3 +?+?== xxxxxx ? , 迭代公式： )12/()14( 2 1 +?+?= + kkkkk xxxxx 6)4(,2)3( =?= ff )4,3(∈x存在根 )(xx ?=? 非线性方程迭代法的基本思想 将原方程 0)( =xf 改写成容易迭代的形式 )(xx ?= , 选 合适的初值 0 x , 进行迭代： ),2,1,0()( 1 L== + kxx kk ? x 0 x 1 x 2 x 3 x 4 x 5 1 ? 3.0000 5.0000 -11.0000 -107.0000 2 ? 3.0000 3.5000 3.1111 3.4054 3.1779 3.3510 3 ? 3.0000 3.2857 3.2749 3.2749 3.2749 3.2749 1 ? 根本不收敛； 2 ? 虽呈现收敛趋势，但很慢； 3 ? 收敛很快。 2749.3 2 )571( , 2 )41411( * 2,1 ≈ +? = ×+±? = xx 非线性方程的迭代法（例） 迭代法的几何解释 y=φ (x) x * x y y=x 0 x 1 x 0 P 0 (x 0 ,x 1 ) x 2 P 1 (x 1 ,x 2 ) x 3 P 2 P 3 x y y=x y=φ (x) 0 x * x 0 x 1 x 2 x 3 P 0 P 1 P 3 {x k }收敛于 x * {x k }不收敛于 x * 取决于曲线 φ(x)的斜率 )( ** xx ?=不动点L,1,0),( 1 == + kxx kk ? 迭代法的收敛性 设 )(xy ?= 在 bxa ≤≤ 连续，且 bya ≤≤ ，若存在 1<L 使 ,)( Lx ≤′? 则 )(xx ?= 在 bxa ≤≤ 有唯一解 * x ，且 1） 对于 ),( 0 bax ∈ ，迭代公式 )2,1,0()( 1 L== + kxx kk ? 产生 的序列 }{ k x 收敛于 * x ； 2） 01 *** 1 1 , xx L L xxxxLxx k kkk ? ? ≤??≤? + 。 局部收敛性 ： 只要 )(x? 在 * x 的一个邻域连续且 ,1)( * <′ x? 则对于该邻域内的任意初值 0 x ，序列 }{ k x 就收敛于 * x 。 L不易确定 放宽定理条件，缩小初值范围 迭代法的收敛速度（收敛阶） 记 * xxe kk ?= ，若 0 1 lim ≠= + ∞→ c e e p k k k （ p 为一正数） 称序列 }{ k x p 阶收敛。 显然， p 越大收敛越快 。 LL +?++?′+= p k p kk xx p x xxxxx )( ! )( ))(()()( *)( *** ? ??? 0)( * ≠′ x? , }{ k x 0)(,0)()( *)(*)1(* ≠===′ ? xxx pp ??? L ， }{ k x p阶收敛 LL +++′= + p k p kk e p x exe ! )( )( *)( * 1 ? ? 1阶收敛(线性收敛) 结论： φ(x)的构造决定收敛速度 )12/()14()( 2 3 +?+?= xxxxx? 迭代法的收敛速度（例） )1/(14)( 2 += xx? 014)( 2 =?+= xxxf例题 阶收敛1}{0)(, )1( 14 )( * 2 2 2 k xx x x ?≠′ + ? =′ ?? 阶收敛2}{0)( ,0)(, )12( )14(2 )( * 3 * 3 2 2 3 k xx x x xx x ?≠′′ =′ + ?+ =′ ? ?? 4 牛顿切线法 )(xf 在 k x 作 Tayl or 展开，去掉 2阶及 2阶以上项得 ))(()()( kkk xxxfxfxf ?′+= )( )( )(, )( )()( )( * * * 2* ** * xf xf x xf xfxf x ′ ′′ =′′ ′ ′′ =′ ?? 牛顿切线法 2阶收敛 ,0)( * =xf 0)(),( ** ≠′′′ xfxf 0)(,0)( ** ≠′′=′ xx ?? 若 x * 为单根 设 0)( ≠′ k xf ，用 1+k x 代替右端的 x ，由 0)( =xf 得 )( )( 1 k k kk xf xf xx ′ ?= + ， 即令 )( )( )( xf xf xx ′ ?=? x k x k+1 M N x O y=f(x) y f(x k ) 牛顿 割线法 用差商 1 1 )()( ? ? ? ? kk kk xx xfxf 代替 )( k xf ′ )()( ))(( 1 1 1 ? ? + ? ? ?= kk kkk kk xfxf xxxf xx x * 为单根时收敛阶数是 1.618 y O x k+1 x y=f(x) P Q x k-1 x k 若 x * 为重根 0)( * ≠′ x? 牛顿切线法 1阶收敛 收敛速度比牛顿切线法稍慢 0)()( ** =′= xfxf 重数越高，收敛越慢 解非线性方程组的牛顿法 T n T n xfxfxFxxxxF )](),([)(,],[,0)( 11 LL === ))(()()( 11 kkkkk xxxfxfxf ?′+= ++ )( )( 1 k k kk xf xf xx ′ ?= + 解方程 f(x)=0 的 牛顿切线法 推广到解方程组 nixx x xf xx x xf xfxf k n k n n k i kk k ik i k i LL ,2,1),( )( )( )( )()( 1 1 1 1 1 1 =? ? ? ++ ? ? ? += + ++ ))(()()( 11 kkkkk xxxFxFxF ?′+= ++ )()]([ 11 kkkk xFxFxx ?+ ′?= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ? =′ × n nn n nn j i x f x f x f x f x f xF LL LL LL 1 1 1 1 ][)( T n T n xfxfxFxxxxF )](),([)(,],[,0)( 11 LL === ))(()()( 11 kkkkk xxxfxfxf ?′+= ++ )( )( 1 k k kk xf xf xx ′ ?= + 解方程 f(x)=0 kkk xxx ?+= +1 )()( kkk xFxxF ?=?′ 解非线性方程组的牛顿法 解方程 f(x)=0 的 牛顿割线法 )()()( 11 kkkkk xxaxfxf ?+= ++ k k kk a xf xx )( 1 ?= + 1 1 )()( ? ? ? ? = kk kk k xx xfxf a 解方程组的拟牛顿法 ——用 A k 代替 F′(x k ） )()()( 11 ?? ?=? kkkkkk xFxFxxAA 满足使 矩阵 A k (n 2 个未知数 )不能由这样的 n个方程确定 kkkk AAAAA 计算有不同的构造用迭代方法 )( 11 ??+= ?? )()( 11 kkkk xFAxx ?+ ?=再求 拟牛顿法 (Quasi-Newton) MATLAB优化工具箱解非线性方程 fzero: 单变量方程 f(x)=0 求根 (变号点) 最简形式x= fzero(@f, x0 ) 可选输入 : “P1,P2,...”是传给 f.m的参数 (如果需要的话 ) ’opt’是一个结构变量，控制参数 (如精度 TolX ) opt可用 optimset设定 , 不指定或指定为 ’[ ]’时将采用缺省值 如 : opt=optimset(‘TolX’,1e-8) 输出: ’fv’是函数值; ’ef’是程序停止运行的原因(1,0,-1); ’out’是一个结构变量，包含: iterations(迭代次数), funcCount (函数调用次数), algorithm(所用算法) 一般形式[x, fv, ef, out] = fzero(@f, x0, opt, P1, P2, ... ) 必须输入 : ’f’为 f.m函数名， ’x0’是迭代初值 (或有根区间 ) 输出: ’x’是变号点的近似值(函数不连续时不一定是根) 演示: exampleFzero.m 5 fsolve: 多变量方程组F(x)=0求解 输出 ---- 与fzero类似, 但: ’out’中还输出’firstorderopt’, 即结果（x点）处梯度向量 的范数(实际上是1-范数，即分量按绝对值取最大的值); ’jac’输出x点所对应的雅可比矩阵 输入 ---- 与fzero类似, 但: ’x0’是迭代初值， ’opt’中控制参数更多(如MaxFunEvals, MaxIter等) x= fsolve(@f, x0 ) 最简形式 [x, fv, ef, out, jac ] = fsolve(@f, x0, opt, P1, P2, ... ) 一般形式 注: solve函数也可求解(符号工具箱) MATLAB优化工具箱解非线性方程组 需自行编制程序，如对切线法编写名为 newton.m的 m文件 牛顿法 当 )(xf 为多项式时可用 r=roots(c) 输入多项式的系数 c（按降幂排列） ，输出 r 为 0)( =xf 的全部根； c=poly(r) 输入 0)( =xf 的全部根 r，输出 c 为多 项 式的系数（按降幂排列） ； df=polyder(c) 输入多项式的系数 c（按降幂排列） ，输出 df为多项式的微分的系数。 多项式求根 MATLAB优化工具箱解非线性方程 function [y,z]=newton(fv,df,x0,n,tol) x(1)=x0; y=x0; b=1; k=1; while abs(b)>tol*abs(x(k)) x(k+1)=x(k)-feval(fv,x(k))/feval(df,x(k)); b=x(k+1)-x(k); y=x(k+1); k=k+1; if(k>n) error('Error: Reached maximum iteration times'); end end k=k-1; fv是 f(x)的函数句柄， df是 f ’(x)的函数句柄 求解 f(x)=0的 newton.m文件 Newton.m 1,4 2 2 2 1 2 2 2 1 =?=+ xxxx解例 ? ? ? ? ? ? ? =′ ? ? ? ? ? ? ? ? ?? ?+ = 21 21 2 2 2 1 2 2 2 1 2)(, 1 4 )( xx xx xF xx xx xF ),()( kkk xFxxF ?=?′ L,1,0, 1 =?+= + kxxx kkk T x )1,1( 0 =取 kx k 0 (1.000000, 1.000000) 1 (1.750000, 1.250000) 2 (1.589286, 1.225000) 3 (1.581160, 1.224645) 4 (1.581139, 1.224745) 5 (1.581139, 1.224745) T T x )22474487.1,58113883.1( )2/3,2/5( = =精确解 演示exam0602Newton.m; exam0602Fsolve.m; exam0602Solve.m 实例1 路灯照明 ()() 0 )( )( 33)(' 5 22 2 22 5 22 1 11 = ?+ ? + + ?= xsh xshP xh xhP xC 2065 32 21 21 === == shh PP ，， ， 0 5 10 15 20 0 0.02 0.04 0.06 0.08 0.1 C(x) 0 5 10 15 20 -4 -2 0 2 4 6 x 10 -3 dC(x) ()() 3 22 2 22 3 22 1 11 )( )( xsh hP xh hP xC ?+ + + = C (x)有 3个驻点 : (9,10)内的是最小点， 0或 20附近的是最大点 实例1 路灯照明 function y=zhaoming(x) y=2*5*x/(5^2+x^2)^(5/2)-3*6*(20-x)/(6^2+(20-x)^2)^(5/2); x0=[0,10,20]; for k=1:3 x(k)=fzero(@zhaoming,x0(k)); c(k)=2*5/(5^2+x(k)^2)^(3/2)+3*6/(6^2+(20-x(k))^2)^(3/2); end [x;c] 0.084474680.084476550.018243930.081981040.08197716 C(x) 2019.976695819.338299140.028489970 x x =9.3383是 C(x)的最小值点， x =19.9767是 C(x)的最大值点 6 实例1 路灯照明 问题： P 2 =3千瓦路灯的高度在 3~9米变化，如何使路面上最暗点 的照度最大 ? 用 fzero命令解方程，得到的结果是： x=9.5032， h 2 =7.4224， C(x, h 2 )= 0.018556(最暗点的最大照度 ) ()() 3 22 2 22 3 22 1 11 2 )( )( xsh hP xh hP hxC ?+ + + =， ()() 5 22 2 22 5 22 1 11 )( )( 33 xsh xshP xh xhP x C ?+ ? + + ?= ? ? ()() 5 22 2 2 22 3 22 2 2 2 )( 3 )( xsh hP xsh P h C ?+ ? ?+ = ? ? =0 ? )( 2 1 2 xsh ?= =0 ? () 0 )( 439 3 2 5 22 1 11 = ? ? + xs P xh xhP 实例1 路灯照明 问题：讨论两只路灯的高度均可以在 3~9米之间变化的情况 实际数据计算，得到 x=9.3253, 最暗点的照度达到最大的路灯高 度 h 1 =6.5940, h 2 =7.5482 ()() 3 22 2 22 3 22 1 11 21 )( ),,( xsh hP xh hP hhxC ?+ + + = ()() 5 22 2 22 5 22 1 11 )( )( 33 xsh xshP xh xhP x C ?+ ? + + ?= ? ? ()() 5 22 2 2 22 3 22 2 2 2 )( 3 )( xsh hP xsh P h C ?+ ? ?+ = ? ? =0 ? )( 2 1 2 xsh ?= =0 ? 0 1 = ? ? h C =0 ? xh 2 1 1 = 3 2 3 1 )( xs P x P ? = 3 2 3 1 3 1 PP sP x + = 实例1 路灯照明 讨论 1： 若 P 1 =P 2 , 则 x=0.5s (中点 )， 与直觉符合 思考： 2只以上路灯的情形（如篮球场四周安装照明灯） )( 2 1 2 xsh ?= xh 2 1 1 = 3 2 3 1 3 1 PP sP x + = x x α 2 α 1 O h 2 P 2 r 1 P 1 s r 2 h 1 y Q 讨论 2： 2/2 21 == αα tgtg 6135 21 ′== o αα （这个角度与路灯的功率和道路宽度均无关 ） 实例2 均相共沸混合物的组分 0ln1ln 1 1 1 = ? ? ? ? ? ? ? ? ? ? ? ? +?? ? ? ? ? ? ? ? ? ? ? ? ? + ? ? ? ? ? ? ? ? + + ∑ ∑ ∑ = = = Pa qx qx qx cT b x i n j n k jkk jij n j ijj i i i a 1 =16.388, a 2 =16.268, a 3 =18.607； b 1 =2787.50, b 2 =2665.54, b 3 =3643.31； c 1 =229.66, c 2 =219.73, c 3 =239.73； P=760（ mmHg） Q=[1.0 0.48 0.768 1.55 1.0 0.544 0.566 0.65 1.0 ] 1 1 = ∑ = n i i x ∑ ? = ?= 1 1 1 n i in xx 给定 n=3种物资：丙酮、乙酸甲脂、甲醇 (分别记为 1、 2、 3) n个变量： T, x i 实例2 均相共沸混合物的组分 XT = [0.2740 0.4636 54.2560] Y = 1.0e-006 * [0.4195 -0.3112 0.2083] function f=azeofun(XT,n,P,a,b,c,Q) x(n)=1; for i=1:n-1 x(i)=XT(i); x(n)=x(n)-x(i); end T=XT(n); p=log(P); for i=1:n d(i) = x * Q(i,1:n)'; dd(i)=x(i)/d(i); end for i=1:n f(i)=x(i)*(b(i)/(T+c(i)) + log(x*Q(i,1:n)') + dd*Q(1:n,i) - a(i) - 1 + p); end n=3; P=760; a=[16.388,16.268,18.607]'; b=[2787.50,2665.54,3643.31]'; c=[229.66,219.73,239.73]'; Q=[1.0 0.48 0.768 1.55 1.0 0.544 0.566 0.65 1.0]; XT0=[0.333,0.333,50]; [XT,Y]=fsolve(@azeofun,XT 0,[],n,P,a,b,c,Q) 实例2 均相共沸混合物的组分 55.67640.00000.46720.5328[0.5, 0.5, 54] 54.50400.25250.00000.7475 [0.5, 0, 54] 54.35790.32340.67660.0000[0, 0.5, 54] 54.25600.26240.46360.2740 [0.333, 0.333, 50] Tx 3 x 2 x 1 XT0 解初值 7 分岔与混沌现象 离散形式的阻滞增长模型（见实验 2） L,2,1,0,)1( 1 =?=? + kx N x rxx k k kk 问题： 种群数量总是趋于稳定吗？ (x k 是第 k代的种群数量， r是固有增长率， N是种群最大容量 ) 取 N=1， r=0.3, 1.8, 2.2, 2.5, 2.55, 2.7，初值 x0=0.1，按照迭代方 程用 MATLAB计算 x k ，观察得到结果 分 岔 与 混 沌 现 象 分岔与混沌现象 k k kk x N x rxx )1( 1 ?=? + 即 1< b <3, 0 < r <2 平衡点 y=y*=1-1/b（相应于 x*=N）稳定的条件为 | f ’( y*)| <1 1, )1( += + = rbx Nr r y kk L,2,1,0 )1()( 1 = ?== + k ybyyfy kkkk 分岔与混沌现象 隔代收敛分析 )(),(,10 * 1 * 2 * 2 * 1 * 2 ** 1 yfyyfyyyy ==<<<< 1))(( * 2,1 )2( <′yf 平衡点 y 1,2 *稳定的条件是 )21)(21()()())(( * 2 * 1 2* 2 * 1 * 2,1 )2( yybyfyfyf ??=′′=′ 449.361 ≈+<b 449.2<r L,2,1,0),())(()( )2( 12 ==== ++ kyfyffyfy kkkk 平衡点 (除 y=y*=1-1/b 外 ) )]1(1)[1()( )2( ybyybbyyfy ???== b bbb y 2 321 2 * 2,1 ??+ = m 分岔与混沌现象 类似地可以得到：迭代方程有 4个稳定平衡点的条件 544.3449.3 <<b 544.2449.2 << r 记有 2 n 个收敛子序列的 b的上限为 b n ，上面的分析给出： 进一步研究表明：当 n→∞时 b n → 3.57，若 b>3.57(即 r>2.57)，就 不再存在任何 2 n 收敛子序列，序列 x k 的趋势似乎呈现一片混 乱，这就是所谓 混沌现象（ Chaos） 。 b 0 =3， b 1 =3.449， b 2 =3.544 L6692.4lim 1 1 = ? ? + ? ∞→ nn nn n bb bb 混沌现象实际上有其内在的规律性 , 如 是普遍存在于不同混沌现象中的常数 (费根鲍姆常数 ) 分岔与混沌现象 混沌现象的一个典型特征是对初始条件的敏感性 ? “差之毫厘，失之千里 ” ? 著名的 “蝴蝶效应 ” 非线性迭代过程 ---- 混沌现象 ---- 非线性科学 8 分岔与混沌现象 function chaos(iter_fun,x0,r,n) % 该函 数没有返回值； iter_fun是迭代函数 （句柄）； x0是迭代初值； kr=0; for rr=r(1):r(3):r(2) % 输入中 [r(1),r(2)]是参数变化的范 围， r(3) 是步长 kr=kr+1; y(kr,1)=feval(iter_fun,x0,rr); for i=2:n(2) %输入中 n(2)是迭代序列的长度，但 画图时前 n(1)个迭代值被舍弃 y(kr,i)=feval(iter_fun,y(kr,i-1),rr); end end plot([r(1):r(3):r(2)],y(:,n(1)+1:n(2)),'k.'); 本例迭代函数为： function y=iter01(x,r) y=r*x*(1-x); 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 输入如下命令： chaos(@iter01,0.5,[2,4,0.01], [100,200]) 布置实验 目的 1) 用 MATLAB软件掌握求解非线性方程的 迭代法和牛顿法，并对结果作初步分析； 内容课上布置，或参见网络学堂 2) 通过实例练习用非线性方程求解的实际问题。