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

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 + − + + =

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)中点非线性方程的基本解法

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 ′′ ≠ ⇒ ′ = + + − ′ = ϕ ϕ ϕ

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

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)的最大值点

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 初值解