《经济数学基础》课程教学资源：实验6 非线性方程求解

什么叫方程(组)? 大学数学实验 方程:包含未知数(或和未知函数)的等式 方程组:包含未知数(或和未知函数)的一组等式 Experiments in Mathematics 不包含未知函数的方程组的一般形式:F(x)=0 x=(x1x2…,x)2,F(x)=(f(x22(x),…,m(x)(一般m=n) 实验6非线性方程求解 满足方程(组)的未知数的取值称为方程(组)的 解,或称为F(x)的零点 半大总数科系 单变量方程(一元方程):f(x)=0,“解”也称为“根 非线性方程的特点 实验6的基本内容 方程分类: ·代数方程:ax+ax1+…+an=0; 1.非线性方程∫(x)=0的数值解法: ·超越方程:包含超越函数如sinx,nx)的方程 迭代方法的基本原理 非线性方程:m(22)次代数方程和超越方程 牛顿法;拟牛顿法 方程根的特点: 2.推广到解非线性方程组 n次代数方程有且只有n个根包括复根、重根); 3.实际问题中非线性方程的数值解 5次以上的代数方程无求根公式 4.分仓和混沌现象 ·超越方程有无根,有几个根通常难以判断 (学学奖 实例1路灯照明 实例1路灯照明 员购灯续线然铸置工最的点亮的总在厘灯 建立坐标系如图,两个光源在点Q(x,0)的照度分别为 如果P的高度可以在米到9米之间变化,如何使路面上最暗 ,l2= P2 sina2是由量纲单位决定的比 例系数,不妨记k 点的亮度最大? 点Q的照度C(x)=+ 如果两只路灯的高度均可以在3米到9米之间变化呢? r2=2+x2,2=+(x-x)2 P PP2,P=3(千瓦) h5,h26(米) h2√2+(-

1 大学数学实验 Experiments in Mathematics 实验6 非线性方程求解 清华大学数学科学系 什么叫方程（组）? 方程：包含未知数（或/和未知函数）的等式 方程组：包含未知数（或/和未知函数）的一组等式 不包含未知函数的方程组的一般形式： F(x)=0 x= (x1,x2,…,xn)T, F(x)=(f1(x),f2(x),…,fm(x))T (一般m=n) 满足方程（组）的未知数的取值称为方程（组）的 解，或称为F(x)的零点。 单变量方程（一元方程）：f(x)=0, “解”也称为“根” 非线性方程的特点 方程根的特点： • n次代数方程有且只有n个根(包括复根、重根)； • 5次以上的代数方程无求根公式； • 超越方程有无根，有几个根通常难以判断。 方程分类： • 代数方程: a0xn+a1xn-1+…+an=0； • 超越方程：包含超越函数(如sinx, lnx)的方程； • 非线性方程：n(≥2)次代数方程和超越方程。 实验6的基本内容 3.实际问题中非线性方程的数值解 1.非线性方程 f(x)=0 的数值解法: • 迭代方法的基本原理; • 牛顿法; 拟牛顿法. 2.推广到解非线性方程组 4.分岔和混沌现象 实例1 路灯照明 道路两侧分别安装路灯，在漆黑的夜晚，当两只路灯开启时， 两只路灯连线的路面上最暗的点和最亮的点在哪里？ h2 P P2 1 s h1 • 如果P2的高度可以在3米到9米之间变化，如何使路面上最暗 点的亮度最大？ • 如果两只路灯的高度均可以在3米到9米之间变化呢？ s=20(米) P1=2, P2=3(千瓦) h1=5, h2=6(米) 实例1 路灯照明 建立坐标系如图，两个光源在点Q(x,0)的照度分别为 (k是由量纲单位决定的比 2 例系数，不妨记k=1) 2 2 2 2 2 1 1 1 1 sin , sin r P I k r P I k α α = = 点Q的照度 2 2 2 2 2 2 2 1 2 1r = h + x , r = h + (s − x) x x α1 α2 O h2 P2 r1 P1 s r2 h1 y Q 2 2 1 1 1 1 1 sin h x h r h + α = = 2 2 2 2 2 2 2 ( ) sin h s x h r h + − α = = ( ) ( )3 2 2 2 2 2 3 2 2 1 1 1 ( ) ( ) h s x P h h x Ph C x + − + + =

实例1路灯照明 实例2均相共沸混合物的组分 为求最暗点和最亮点,先求C(x)的驻点 均相共沸混合物〔 homogeneous azeotrope)是由两种或两种 C(x)= Ph,(s-x) 以上物质组成的液体混合物,当在某种压力下被燕馆或局部汽 化时,在气体状态下和在液体状态下保持相同的组分(比例) 给定几种物质,如何确定它们构成均相共沸混合物时的比例? 令C'(x)=0:解析解是难以求出,需数值求解 模型建立 设该混合物由n个可能的物质组成,物质所占的比例为x x2=1,x1≥0 ag 实例2均相共沸混合物的组分 实例2均相共沸混合物的组分 均相共沸混合物应该足穗定条件,即共沸混合物的每个组分 在气体和液体状态下具有相同的化学勢能。在压强P不大的情 P=nP InP=a 况下,这个条件可以表示为:P=yP,1=1…,n P是物质的饱和汽相压强,与温度有关,可以如下确定 2--1-a+hP=0 i=1,…,n(anb,c是常数 Y是组分液相活度系数,可以根据如下表达式确定 只有当物质选与到该共沸混合物中时才需要满足上式,故得到 lnx=1-2x9,-∑2-=1…n (q表示组分;与组分的交互作用参数,可以通过实验近似得到) ∑x=1,x20 解析解是难以求 i=1,2,,n 出,需数值求解 (学学奖 非线性方程的基本解法 非线性方程的基本解法 解方程几x)=0第一步确定根的大致范 分法若对于a0 图形法 的实现」4x0 xx0:(a,b)中 x6-2x4-6x3-13x2+8x+12=0 f(x)0da1=a,b=x01x0)<0a1=x0,b=b (a,b)=(a,b)→…→(an,bn)→…区间每次缩小一半,呢足 有4个根分别位于 够大时,可确定根的范围 -1.75,-0.75,1.00,240附近

2 实例1 路灯照明 为求最暗点和最亮点，先求C(x)的驻点 x x α1 α2 O h2 P2 r1 P1 s r2 h1 y 令C’ (x)=0：解析解是难以求出，需数值求解 Q ( ) ( )3 2 2 2 2 2 3 2 2 1 1 1 ( ) ( ) h s x P h h x Ph C x + − + + = ( ) ( )5 2 2 2 2 2 5 2 2 1 1 1 ( ) ( ) '( ) 3 3 h s x P h s x h x Ph x C x + − − + + = − 实例2 均相共沸混合物的组分 均相共沸混合物（homogeneous azeotrope） 是由两种或两种 以上物质组成的液体混合物，当在某种压力下被蒸馏或局部汽 化时，在气体状态下和在液体状态下保持相同的组分（比例） 给定几种物质，如何确定它们构成均相共沸混合物时的比例？ 设该混合物由n个可能的物质组成，物质i所占的比例为xi 模型建立 1, 0 1 ∑ = ≥ = i n i i x x 实例2 均相共沸混合物的组分 均相共沸混合物应该满足稳定条件，即共沸混合物的每个组分 在气体和液体状态下具有相同的化学势能。在压强P不大的情 况下，这个条件可以表示为： Pi 是物质i的饱和汽相压强，与温度T有关，可以如下确定： P P i n i i = γ , = 1,L, i n T c b P a i i i i ln , = 1,L, + = − i 是组分i的液相活度系数，可以根据如下表达式确定: γ i n x q x q x q n j n k k jk j ji n j i j ij ln 1 ln , 1, , 1 1 1 = L ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ −⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − ∑ ∑ ∑ = = = γ (qij表示组分i与组分j的交互作用参数 ，可以通过实验近似得到) (ai , bi ，ci 是常数) 实例2 均相共沸混合物的组分 P iPi = γ i i i i T c b P a + ln = − ∑ ∑ ∑ = = = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ −⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − n j n k k jk j ji n j i j ij x q x q x q 1 1 1 lnγ 1 ln 只有当物质i参与到该共沸混合物中时才需要满足上式，故得到 ln 1 ln 0 1 1 1 − − + = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ +⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + ∑ ∑ ∑ = = = a P x q x q x q T c b i n j n k k jk j ji n j j ij i i ln 1 ln 0 1 1 1 = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ − − + ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ +⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + ∑ ∑ ∑ = = = a P x q x q x q T c b x i n j n k k jk j ji n j j ij i i i x x i n i n i i 1, 0, 1,2,..., 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 2 6 13 8 12 0 6 4 3 2 x − x − x − x + x + = 图形法 有4个根分别位于 x= -1.75, -0.75, 1.00，2.40附近 Exam0600.m 若对于 a 0 ( ) 0 f (x0) > 0 f x0 < (a,b) ⇒(a1,b1)⇒L⇒(an,bn)⇒L 1 1 0 a = a, b = x a = x b = b 1 0 1 , 区间每次缩小一半，n足 够大时，可确定根的范围 a b f (x) x a b x f (x) 0x 0x x0 :(a,b)中点 非线性方程的基本解法

非线性方程选代法的基本思想 非线性方程的选代法(例) 将原方程f(x)=0改写成容易迭代的形式x=g(x),选 合适的初值x,进行选代:x=9(x.)(k=0,2,…) 例1 f(x)=x2+x-14=0→x=9(x) 3.00003.50003.Ill3.40543.17793.3510 3.000032857 f(3)=-2,f(4)=6存在根x∈(3,4) x=q(x)=14-x2,迭代公式:x=14-x2 1t0+1x:=(1+4573279 x=q(x)=14/(x+1),迭代公式xa=14(x+D) x=,(x)=x-(x2+x-14)/2x+1) 根本不收敛;g2虽呈现收敛趋势,但很慢:φ,收敛很快 迭代公式:x,=x4-(x+x4-14)2x+1) 选代法的几何解释 迭代法的收敛性 x=(x4),k=0,1 不动点x=c(x) 设y=φ(x)在a≤x≤b连续,且a≤y≤b,若存在L<1使 p(x)≤L,则x=9(x)在a≤x≤b有唯一解x,且 1)对于x∈(a,b),迭代公式x=p(x)(k=0,1,2…)产生 的序列{x}收敛于x (x1x2) 2)k-x4-x-x11k-x L不易确定口放宽定理条件,缩小初值范围 局部收敛性:只要φ(x)在x的一个邻域连续且 1 收敏于x 取决于曲线的斜率不收敛于 p(x)<,则对于该邻域内的任意初值xn,序列 x}就收敛于x。 (学学奖 (学学实验) 迭代法的收敛速度(收敛阶 选代法的收敛速度(例) 记e,=x.-x,若m=c≠0(P为一正数) 例题f(x)=x2+x-14=02(x)=14/(x+1) 称序列{x}P阶收敛。显然,P越大收敛越快 92(x)= 叭4)=mxx-x)+…+gB x+1202(x)=0={xk;1阶收敛 O(P)(r) 02(x)=x-(x2+x-14)/(2x+1 ek+1=q(x)ek+…+ q(x)≠0,{x,}I阶收戴跳他 93(x) 2(x2+x-14) q(x)=…=q"(x)=0,q"(x)≠0,{x}P阶收敛 g(x)≠0→{xk}2阶收敛 站论:()的构造决定收斂遗度

3 例1 ( ) 14 0 2 f x = x + x − = 2 x =ϕ1 (x) =14− x , 迭代公式： 2 1 14 k k x = − x + ( ) 14/( 1) x =ϕ2 x = x+ ，迭代公式： 14/( 1) xk+1 = xk + ( ) ( 14) /(2 1) 2 x = ϕ 3 x = x − x + x − x + , 迭代公式： ( 14) /(2 1) 2 xk +1 = xk − xk + xk − xk + f (3) = −2, f (4) = 6 存在根 x ∈ (3,4) ⇒ x = ϕ(x) 非线性方程迭代法的基本思想 将原方程 f (x) = 0 改写成容易迭代的形式 x = ϕ(x) , 选 合适的初值 0 x , 进行迭代： ( ) ( 0,1,2, ) xk +1 = ϕ xk k = L x0 x1 x2 x3 x4 x5 ϕ 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 收敛很快。 3.2749 2 ( 1 57) , 2 ( 1 1 14 4) * 1,2 ≈ − + = − ± + × x = x 非线性方程的迭代法（例） 迭代法的几何解释 y=φ(x) x* x y y=x 0 x1 x0 P0(x0,x1) x2 P1(x1,x2) x3 P2 P3 x y y=x y=φ(x) 0 x* x0 x1 x2 x3 P0 P1 P3 {xk}收敛于x* {xk}不收敛于x* 取决于曲线 φ(x)的斜率 ( ) * * xk+1 = ϕ(xk ), k = 0,1,L 不动点x =ϕ x 迭代法的收敛性 设 y = ϕ(x) 在 a ≤ x ≤ b 连续，且 a ≤ y ≤ b ，若存在 L < 1 使 ϕ′(x) ≤ L, 则 x = ϕ(x) 在 a ≤ x ≤ b 有唯一解 * x ，且 1） 对于 ( , ) 0 x ∈ a b ，迭代公式 ( ) ( 0,1,2 ) xk +1 = ϕ xk k = L 产生 的序列{ }k x 收敛于 * x ； 2） 1 0 * * * 1 1 , x x L L x x L x x x x k k k k − − + − ≤ − − ≤ 。 局部收敛性：只要ϕ(x) 在 * x 的一个邻域连续且 ( ) 1, * ϕ′ x < 则对于该邻域内的任意初值 0 x ，序列 { }k x 就收敛于 * x 。 L不易确定 放宽定理条件，缩小初值范围 迭代法的收敛速度（收敛阶） 记 * e x x k = k − ，若 0 1 lim = ≠ + → ∞ c e e p k k k （ p 为一正数） 称序列 { }k x p 阶收敛。显然， p 越大收敛越快。 = + ′ − +L+ k − p +L p k k x x p x x x x x x ( ) ! ( ) ( ) ( ) ( )( ) ( ) * * * * ϕ ϕ ϕ ϕ ( ) 0 * ϕ ′ x ≠ ,{ }k x ( ) ( ) 0, ( ) 0 * ( 1) * ( ) * ′ = = = ≠ − x x x p p ϕ L ϕ ϕ ，{ }k x p 阶收敛 + = ′ +L+ k p +L p k k e p x e x e ! ( ) ( ) ( ) * * 1 ϕ ϕ 1阶收敛(线性收敛) 结论：φ(x)的构造决定收敛速度 ( ) ( 14)/(2 1) 2 ϕ3 x = x − x + x − x + 迭代法的收敛速度（例） ( ) 14 /( 1) ( ) 14 0 ϕ2 x = x + 2 例题 f x = x + x − = , ( ) 0 { } 1阶收敛 ( 1) 14 ( ) * 2 2 2 k x x x x ′ ≠ ⇒ + − ϕ ′ = ϕ ( ) 0 { }2阶收敛 , ( ) 0, (2 1) 2( 14) ( ) * 3 * 3 2 2 3 k x x x x x x x ′′ ≠ ⇒ ′ = + + − ′ = ϕ ϕ ϕ

牛顿切线法 若x为重根f(x’)=f(x)=0口(x≠0 f(x)在x作Ta展开,去掉2阶及2阶以上项得 f(x)=f()+f(x x-x) 牛顿切线法1阶收斂重数越高,收敛越慢 设∫'(x,)≠0,用x代替右端的x,由f(x)=0得 牛顿 f(x 即令g(x)=x-(x) 制线法用差商(x)-/()代替f(x) o(x)=x1).p(x)=x f(x4)(x4-x,-) f(x) y=若x为单根fx)=0,f(x)f(x)≠0 收敛速度比牛顿切线法稍慢 口g(x)=0,9(x)≠0 xx为单根时收敛阶数是1618 牛顿切线法2阶收敛 大学学实编) 解非线性方程组的牛顿法 解非线性方程组的牛顿法 解方程x)=0f(xk+1)=f(x4)+f(xk(xk+1-xk) F(x)=0,x=[x,…xn},F(x)=[1(x)…fn(x) 的牛顿切线法 k f(r) 9.….9解方程x)=0 推广到解方程组 F(x)=0,x=[x,…xn},F(x)=[1(x)…fn(x) n可fn f(x4)=f(x2)+ (x11-x4) F(x2+)=F(x2)+F(x)(x1-x2) +…+可(x)-x2)1=1.2,…n r-F(rf() E F()Ar=-F(r) x +=r+Ar (学学奖 大学费学买验 拟牛顿法( Quasi-Newton) MATLAB优化工具箱解非线性方程 解方程∫(x)=0 的牛顿割线法f(x+1 fzero:单变量方程f(x)=0求根变号成) 录彩或 xu fzero(afx0) f(x) a4=(x)-f(x和 必须输入:'为m函数名,'x0是选代初值(或有根区间 k k-l 一形式|x,f,ef, outl-fzero(af,x0,opt,Pl,P2,…) 解方程组的拟牛顿法—用A代誉F(x4) 可选输入:“P1,P2”是传给Lm的参数(如果需要的话) pr是一个结构变量,控制参数(如精度Io) 使满足A(x4-xk-1)=F(x)-F(x- op可用 optimise设定,不指定或指定为时将采用缺省值 矩阵A(n2个未知数)不能由这样的n个方程确定 用迭代方法A4=A4-1+A4-1(△4有不同的构造)计算A 输出:'’是腾;'e『是程序尊止觉行的圆(1,01); 是一个首物量,包 再求x+1=x-()-1F(xk iterations代次款), uncount( 宗: exampleFzero.m

4 牛顿切线法 f (x)在 k x 作Taylor 展开，去掉2阶及2阶以上项得 ( ) ( ) ( )( ) k k k f x = f x + f ′ x x− x ( ) ( ) , ( ) ( ) ( ) ( ) ( ) * * * * 2 * * * f x f x x f x f x f x x ′ ′′ ′′ = ′ ′′ ϕ′ = ϕ 牛顿切线法2阶收敛 ( ) 0, * f x = ( ), ( ) 0 * * f ′ x f ′′ x ≠ ( ) 0, ( ) 0 * * ϕ′ x = ϕ′′ x ≠ 若x*为单根 设 ′( ) ≠ 0 k f x ， 用 k + 1 x 代替右端的 x ， 由 f ( x ) = 0 得 ( ) ( ) 1 k k k k f x f x x x ′ + = − ， 即 令 ( ) ( ) ( ) f x f x x x ′ ϕ = − xk xk+1 M N x O y=f(x) y f(xk) 牛顿 割线法 用差商 1 1 ( ) ( ) − − − − k k k k x x f x f x 代替 ( ) k f ′ x ( ) ( ) ( )( ) 1 1 1 − − + − − = − k k k k k k k f x f x f x x x x x x* 为单根时收敛阶数是 1.618 y O xk+1 x y=f(x) P Q xk-1 xk 若x*为重根 ( ) 0 * ϕ′ x ≠ 牛顿切线法1阶收敛 收敛速度比牛顿切线法稍慢 ( ) ( ) 0 * * f x = f ′ x = 重数越高，收敛越慢 解非线性方程组的牛顿法 T n T n F(x) 0, x [x , x ] , F(x) [ f (x), f (x)] = = 1 L = 1 L ( ) ( ) ( )( ) k 1 k k k 1 k f x = f x + f ′ x x − x + + ( ) 1 ( ) k k k k f x f x x x ′ = − + 解方程 f(x)=0 的牛顿切线法 推广到解方程组 x x i n x f x x x x f x f x f x k n k n n k i k k k k i i k i L ( ), 1,2,L ( ) ( ) ( ) ( ) ( ) 1 1 1 1 1 1 − = ∂ ∂ + + − ∂ ∂ = + + + + ( ) ( ) ( )( ) k 1 k k k 1 k F x = F x + F′ x x − x + + [ ( )] ( ) k 1 k k 1 k x x F x F x + − = − ′ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = ∂ ∂ ′ = × n n n n n n j i x f x f x f x f x f F x LL LL LL 1 1 1 1 ( ) [ ] T n T n F(x) 0, x [x , x ] , F(x) [ f (x), f (x)] = = 1 L = 1 L ( ) ( ) ( )( ) k 1 k k k 1 k f x = f x + f ′ x x −x + + ( ) 1 ( ) k k k k f x f x x x ′ = − + 解方程 f(x)=0 k k k x = x + ∆x +1 ( ) ( ) k k k F′ x ∆x = −F x 解非线性方程组的牛顿法 解方程 f(x)=0 的牛顿割线法 ( ) ( ) ( ) k 1 k k k 1 k f x = f x + a x − x + + k k k k a f x x x 1 ( ) = − + 1 1 ( ) ( ) − − − − = k k k k k x x f x f x a 解方程组的拟牛顿法——用Ak代替F′(xk） ( ) ( ) ( ) −1 −1 − = − k k k k k k 使A 满足 A x x F x F x 矩阵Ak(n2个未知数)不能由这样的n个方程确定 k k k k 用迭代方法A A A ( A有不同的构造)计算A 1 1 = + ∆ ∆ − − ( ) ( ) k 1 k k 1 k x x A F x + − 再求 = − 拟牛顿法(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

MATLAB优化工具箱解非线性方程组 MATLAB优化工具箱解非线性方程 fsolve:●量F(x)=0求 e fsolve(af, x0) 量彩式 需自行编制程序,如对切线法编写名为 newton,m的m文件 Ix, fv, ef, out, jac I -fsolve(@f, xO, opt, Pl, P2, 多项式求根 输入 r=roots(c) 输入多项式的系数c(按降幂排列),输出 op中蛇制数冥(如 MaxFun Evals,Maer等) 输出一与 fzero类似 r为f(x)=0的全部根 ou中场着出' firstorderopt,即最(点)是弹度向量 入f(x)=0的全部根r,输出c为多项 (限上是1,P分量续隐增取最大的良 式的系数(按降幂排列) jac’着出x胤开对度的可比气阵 df= polder(c)输入多项式的系数c(按降幂排列),输出 df为多项式的微分的系数。 滤:Sol函也可求市奇子工具) 求解八x)=0的 newton. m文件 例解x+x2=4,x2-x2=1 function ly, z]-newton(fv, df, xO, n, tol) F(x)= x+xi F(x) x(IxO; y-x0; b-l: k-l: while abs(bAtol"abs(x(k)) x(k+lx(k)-feval(fv, x(k)feval(df, x(k)); F(x2)△x2=-F(x2),x4+=x4+Ax,k=0,1 x(k+l)-x(k); x(k+1); 被录exam0602 Newton.m;0|0000 error(Error: Reached maximum iteration times): exam0602 Fsolve m 2.589286,122500 kek-1; 精确解x=(√5/2,√3/2 3(s8110.1.22545 4(15811,1224745) fv是f(x)的函数句柄,d是f(x)的函数句柄 =(581183122444755819124745 (学学奖 实例1路灯照明 实例1路灯照明 function y=zhaoming(x) 2+(-x)2)h=5,h2=6,=20 y=2+5*x(5^2+x^2)(5/2)3°6°(20-x62+(20-x)2)(5/2 C"(x)= xO=0,10,20 +x2)V2+(x-x x(kizero(zhaoming, xo(k)), enc(k=2*5/(5~2+x(k)2y(3/2H+36/62+(20-x(k)2y(3 0.028489979.3382991419.9766958120 16142434so474 C(x)有3个脏点:(9,10)内的是最小点,0或20附近的是最大点 x=93383是C(x)的最小值点,x=19.9767是C(x)的最大值点

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 文件 牛顿法 当 f (x) 为多项式时可用 r=roots(c) 输入多项式的系数 c（按降幂排列），输出 r 为 f (x) = 0 的全部根； c=poly(r) 输入 f (x) = 0 的全部根 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 4, 1 2 2 2 1 2 2 2 例 解 x1 + x = x − x = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − ′ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − − + − = 1 2 1 2 2 2 2 1 2 2 2 1 , ( ) 2 1 4 ( ) x x x x F x x x x x F x ( ) ( ), k k k F′ x ∆x = −F x xk+1 = xk + ∆xk , k = 0,1,L T x (1,1) 0 取 = k xk 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 (1.58113883,1.22474487) ( 5 / 2, 3/ 2) = 精确解 = 演示 exam0602Newton.m; exam0602Fsolve.m; exam0602Solve.m 实例1 路灯照明 ( ) ( ) 0 ( ) ( ) '( ) 3 3 5 2 2 2 2 2 5 2 2 1 1 1 = + − − + + = − h s x P h s x h x Ph x C x 5 6 20 2 3 1 2 1 2 = = = = = h h s P P ， ， ， 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 2 2 2 2 2 3 2 2 1 1 1 ( ) ( ) h s x P h h x Ph C x + − + + = 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] C(x) 0.08197716 0.08198104 0.01824393 0.08447655 0.08447468 x 0 0.02848997 9.33829914 19.97669581 20 x =9.3383是C(x)的最小值点，x =19.9767是C(x)的最大值点

实例1路灯照明 实例1路灯照明 问题:P=3千瓦路灯的高度在39米变化,如何使路面上最暗点 问题:讨论两只路灯的高度均可以在3-9米之间变化的情况 C(r,h,)=-Rih h=+ +(-x =0→h2=(3-x) +(-x+(4-xF =0→h2=(-x) 9√5PA +x)√时+(-x 用ero命令解方程,得到的结果是 x=9.5032,h2-7.224,Cx,h2=0.01855最暗点的最大照度) 实际数据计算,得到r=9.3253,最暗点的照度达到最大的路灯高 度h=65940,h2=75482 实例1路灯照明 实例2均相共沸混合物的组分 h2=(-x) 1-a,+InP=0 ∑x=1口x=1-∑x n个变量:T,x 讨论1:若PP2则=0.5(中点),与直觉符合 给定m=3种物资:丙酮、乙酸甲脂、甲醇(分别记为1、2、3) 讨论2:ga1=8a2=212a1=a2=3516 (这个角度与路灯的功率和道路宽度均无关) a=16388,a2=16.268,a3=18.607:P=760(mmHg) b=2787.50,b2=2665.54 b2=3643.31 1.551.00.544 思考:2只以上路灯的情形(如篮球场四周安装照明灯) c1=229.66,c2219.73,c3=23973 0.5660.651.0 (学学奖 大学费学买验 实例2均相共沸混合物的组分 实例2均相共沸混合物的组分 388,16.268,1860 X(IFXT(O5 c=2296621973239.73T TEXT(n: Q=[1.00480.768 03300274004636 54.2560 Eog(Pk 00554000040676603234543579 X0=10.30.3350] d0)=x· Q(LIny 15040747500002525545040 dd(ixid(i5 I, n, P, a,b,c, QI 1050554_0532046720556764 i=x(°(bu(T dQ( i)-a(n lou(xQL,I:n)xr=10.2740046365425601 10.4195.311202083

6 实例1 路灯照明 问题：P2=3千瓦路灯的高度在3~9米变化，如何使路面上最暗点 的照度最大? 用fzero命令解方程，得到的结果是： x=9.5032，h2=7.4224，C(x, h2)= 0.018556(最暗点的最大照度) ( ) ( )3 2 2 2 2 2 3 2 2 1 1 1 2 ( ) ( ) h s x P h h x Ph C x h + − + + ， = ( ) ( )5 2 2 2 2 2 5 2 2 1 1 1 ( ) ( ) 3 3 h s x P h s x h x Ph x x C + − − + + = − ∂ ∂ ( ) ( )5 2 2 2 2 2 2 3 2 2 2 2 2 ( ) 3 ( ) h s x P h h s x P h C + − − + − = ∂ ∂ =0 Î ( ) 2 1 2 h = s − x =0 Î ( ) 0 ( ) 9 3 4 3 2 5 2 2 1 1 1 = − − + s x P h x P h x 实例1 路灯照明 问题：讨论两只路灯的高度均可以在3~9米之间变化的情况 实际数据计算，得到x=9.3253, 最暗点的照度达到最大的路灯高 度 h1=6.5940, h2=7.5482 ( ) ( )3 2 2 2 2 2 3 2 2 1 1 1 1 2 ( ) ( , , ) h s x P h h x Ph C x h h + − + + = ( ) ( )5 2 2 2 2 2 5 2 2 1 1 1 ( ) ( ) 3 3 h s x P h s x h x Ph x x C + − − + + = − ∂ ∂ ( ) ( )5 2 2 2 2 2 2 3 2 2 2 2 2 ( ) 3 ( ) h s x P h h s x P h C + − − + − = ∂ ∂ =0 Î ( ) 2 1 2 h = s − x =0 Î 0 1 = ∂ ∂ h C =0 Îh x 2 1 1 = 3 2 3 1 (s x) P x P − = 3 2 3 1 3 1 P P P s x + = 实例1 路灯照明 讨论1：若P1=P2, 则x=0.5s (中点)， 与直觉符合 思考：2只以上路灯的情形（如篮球场四周安装照明灯） ( ) 2 1 2 h = s − x h x 2 1 1 = 3 2 3 1 3 1 P P P s x + = x x α1 α2 O h2 P2 r1 P1 s r2 h1 y Q 讨论2： 2 / 2 tgα1 = tgα 2 = 1 = 2 = 35 16′ o α α （这个角度与路灯的功率和道路宽度均无关 ） 实例2 均相共沸混合物的组分 ln 1 ln 0 1 1 1 = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ − − + ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ +⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + ∑ ∑ ∑ = = = a P x q x q x q T c b x i n j n k k jk j ji n j j ij i i i a1=16.388, a2=16.268, a3=18.607； b1=2787.50, b2=2665.54, b3=3643.31； c1=229.66, c2=219.73, c3=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 n i x x 给定n=3种物资：丙酮、乙酸甲脂、甲醇(分别记为1、2、3) n个变量：T, xi 实例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 均相共沸混合物的组分 [0.5, 0.5, 54] 0.5328 0.4672 0.0000 55.6764 [0.5, 0, 54] 0.7475 0.0000 0.2525 54.5040 [0, 0.5, 54] 0.0000 0.6766 0.3234 54.3579 [0.333, 0.333, 50] 0.2740 0.4636 0.2624 54.2560 XT0 x1 x2 x3 T 初值 解

网y 分岔与混沌现象 高散形式的阻滞增长模型(见实验2) x4-x4=r(1-3)x4,k=012 x是第代的种群数量,是固有增长率,N是种群最大容量) 问题:种群数量总是趋于稳定吗? 分岔与混沌现象 取№1,0.3,1.8,22,25,255,2.7,初值0=0.1,按照选代方 程用 MATLAB计算x观察得到结果 1020304050 1020304050 分岔与混沌现象 分岔与混沌现象 xe-X =( -XE)X* yu+=f(A)=by, (1-y) 隔代收敢分析 y42=f(V)=f(f(y)=f(y4),k=0,2 b k=0.L2 y=f(y)=bby(1-y)1-b(1-y) 平衡点(除y-1-1/b外) b+1千√b2-2b-3 平衡点yy=1-/b(相应于x*=N)穗定的条件为|(y)357即257) 不再存在任何收做子序列,序列x的趋势似乎呈现一片混 乱,这就是所谓混 混沌现象实际上有其内在的规律性如 limb-bm1=4692 =0101 是普遍存在于不同混池现象中的常数(费根鲍姆常数) 非线性迭代过程一混池现象一非线性科学

7 分岔与混沌现象 离散形式的阻滞增长模型（见实验2） +1 − = (1− )x , k = 0,1,2,L N x x x r k k k k 问题：种群数量总是趋于稳定吗？ (xk是第k代的种群数量，r是固有增长率，N是种群最大容量) 取N=1，r=0.3, 1.8, 2.2, 2.5, 2.55, 2.7，初值x0=0.1，按照迭代方 程用MATLAB计算xk，观察得到结果 分 岔 与 混 沌 现 象 分岔与混沌现象 k k k k x N x x x r(1 ) +1 − = − 即13.57(即r>2.57)，就 不再存在任何2n收敛子序列，序列xk的趋势似乎呈现一片混 乱，这就是所谓混沌现象（Chaos）。 b0=3，b1=3.449，b2=3.544 lim 4.6692L 1 1 = − − + − →∞ n n n n n b b 混沌现象实际上有其内在的规律性 b b , 如 是普遍存在于不同混沌现象中的常数(费根鲍姆常数) 分岔与混沌现象 混沌现象的一个典型特征是对初始条件的敏感性 • “差之毫厘，失之千里” • 著名的“蝴蝶效应” 非线性迭代过程 ---- 混沌现象 ---- 非线性科学

分岔与混沌现象 布置实验 数没值代该例选代通数 目 x0是达代初值 y=rx·(l-x) 用 MATLAB软件掌握求解非线性方程的 输入如下命令: 2)是参数变化的范 迭代法和牛顿法,并对结果作初步分析 围,(3)是步长 chaos(@ iter01,0.5, 2, 4,0.011, 00200D 2)通过实例练习用非线性方程求解的实际问题 wkr, IFfeval(iter fun,xo rr) 內率课上布量,或参见酬络学堂 %输入中m(2)是达代序列的长度,但 y(kr, iFfeval( iter fun, y(kr, i-I ),rk plo(rr3)2月3y(:n(1)+1m2)k方 8

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) 通过实例练习用非线性方程求解的实际问题