实验六-非线性方程求解

实验六非线性方程求解一、实验问题1.给定4种物质对应的参数ai,bi,ci和交互作用矩阵Q如下。在压强p=760mmHg下，为了形成均相共沸混合物，温度和组分分别是多少？请尽量找出所有可能解。ki1234a18.60715.84120.44319.293b2643.312755.644628.964117.07c239.73219.16252.64227.44𝑄=[1.00.1920.3161.02.1691.61610.4770.5240.3770.5240.2820.2821.02.0650.2961.0]2.假设商品在t时期的市场价格为p(t)，需求函数为D(p(t))=c-dp(t)(c,d0).而生产方的期望价格为q(t),供应函数为S(q(t)).当供销平衡时S(q(t))=D(p(t)).若期望价格与市场价格不符,商品市场不均衡，生产方t+1时期的期望价格将会调整，方式为q(t+1)-q(t)=r[p(t)q(t)](0r1)，以p(t)=[c-D(p(t))]/d=[c-S(q(t))]/d代入，得到关于q(t)的递推方程.设S(x)=arctan(μx),μ=4.8,d=0.25,r=0.3,以c为可变参数，讨论期望价格q(t)的变化规律，是否有混沌现象出现，并找出前几个分岔点，观察分岔点的极限趋势是否符合Feigenbaum常数揭示的规律.二、实验目的1.掌握用MATLAB软件求解非线性方程和方程组的基本用法，并对结果作初步分析。2.练习用非线性方程和方程组建立实际问题的模型并进行求解。三、实验假设线性代数方程组是一类最简单的方程组,在实际问题中，非线性代数方程和方程组大量存在，研究其有效求解方法是数值计算的基本任务之一.方程中的未知数也称为变量或变元，只含一个未知数的方程(即一元方程或单变量方程)可以记作f(x)=0.(1)该方程的解也称为方程的根(或函数f(x)的零点).当f(x)为k次多项式时，(1)式称为k次代数方程.对于3次、4次方程，虽然可以从数学手册中或在互联网上查到求根公式，但是太复杂，一般不会有人用它;至于5次以上的方程就没有现成的求根公式了.求高次代数方程的根是一个基本的、古老的数学问题，比如在线性代数中求k阶矩阵的特征根时就要解k次代数方程.根据代数基本定理,k次代数方程一定有k个根(包括复根，且重根应按重数计算根的个数)，但是不一定能很快算出它的根.一些由实际问题列出的方程,或由其他数学问题归结得到的方程中还常常包含三角函数、指数函数、对数函数等超越函数,如sinx,ex,Inx等，叫做超越方程,它与k(≥2)次代数方程都是非线性方程.求解超越方程不仅没有一般的公式，很难求出解析解，而且若只依据方程本身，那么连有没有解、有几个实根解，也难以判断.包含n个未知数的m个方程称为方程组，可以记作F(x)=0，(2)其中x=(x1,x2....xn)T是一个向量,F(x)=(f1(x),f2(x)...fm(x))T是一个向量值函数.当f1(x),f2(x)...fm(x)中至少有一个非线性函数时，方程组(2)称为非线性方程组。多数情况下，方程组中包含的方程的个数等于未知数的个数(即m=n).求解非线性方程(组)的一般方法是迭代法，其选代公式实际上可看作一种非线性差分方程.采用非线性迭代法时，会出现类非常有趣的现象——分岔和混沌现象。1.非线性方程(组)的基本解法1.1图形法与二分法解方程f(x)=0的第一歩通常是确定根的近似位置或大致范围.利用MATLAB的函数图形功能就能帮助我们快速判断方程有没有根,并且确定根的近似位置.例如:X6-2x4-6x3-13x2+8x+12=O,(14)X2-2x-4lnx=0.(15)用MATLAB作出方程(14),方程(15)左端函数f(x)的图形(x的范围需经适当调整确定)。从曲线f(x)与x轴的交点可以看出.在图形所表示的x的范围内,方程(14)有4个根;方程(15)有两个根。确定一般的非线性方程根的近似位置或大致范围的另一种方法.是基于连续函数f(x)的如下性质:设f(x)是连续函数.若对于ab,有f(a)*f(b)0,则在(a,b)内f(x)至少有一零点.即f(x)=0至少有一个根.通过试探，确定区间(a,b)后,可以用简单的二分法将区间缩小,具体歩骤如下:取(a,b)的中点x0=(a+b)/2若f(x0)=0,则x0即是根.否则,如f(a)*f(b)0,令a1=a.b1=x0;如f(x0)*f(b)0,令a1=x0.b1=b.在(a1,b1)内至少有一个根.且(a,b)⊃(a1,b1).再取(a1,b1)的中点x1=(a1+b1)/2,如此进行下去.包含根的区同(an,bn)的长度每次缩小一半(n=1,2...),足够大时即可达到满意的精度。采用二分法时，区间中点序列{xn}理论，上将收敛到根的真值，但收敛速度较慢，所以通常用作下面介绍的迭代方法的初步近似(初值).1.2迭代法将原方程f(x)=0改写成等价形式x=p(x)，选择适当的初值x0，按照迭代公式Xk+1=p(xk)，k=0,1...(16)计算,若迭代序列{xn}收敛到x*，则x*满足x*=p(x*),x*称为迭代函数p的不动点，即为原方程f(x)=0的根.问题的关键在于如何构造迭代函数p,使迭代序列{xn}以较快速度的收敛。下面先看一个例子.例1求方程f(x)=x2+x-14的一个正根.解:因为f(3)=-2,f(4)=6，所以有一个正根x*∈(3,4).实际上，用二次方程求根公式容易计算正根x*=(570.5-1)/2≈3.2749.下面我们尝试构造不同的迭代函数p(x)，再以x0=3为初值,用迭代法求解.x=p1(x)=14-x2，迭代公式:xk+1=14-x2k;x=p2(x)=14/(x+1),迭代公式:xk+1=14/(xk+1);x=p3(x)=x-(x2+x-14)/(2x+1)，迭代公式:Xk+1=xk-(x2k+xk+14)/(2xk+1).将这3个迭代公式前5步的计算结果列入表,可以看出:p1产生的{xn}根本不收敛;p2产生的{xn}虽呈现收敛趋势，但很慢;p3产生的{xn}收敛很快。x0x1x2x3x4x5p135-11-107p233.53.11113.40543.17793.351p333.28573.27493.27493.27493.2749为了研究序列{xn}收敛或发散的本质，通过几何图形观察按照迭代公式(16)进行迭代的过程。通过作图可以发现，在x*附近，若曲线y=p(x)的斜率(取绝对值)小于直线y=x的斜率(为1)，则序列{xn}将收敛;反之，当y=p(x)的斜率(取绝对值)大于1时，序列{xn}不收敛.关于迭代法的收敛性,理论上有如下的局部收敛性定理.设p(x)在x*的一个邻域内连续可微.且|p|1.则対于该邻域内的任意初值x0,序列{xn}收敛于x*.下面介绍收敛速度的概念.对Rn中的序列{xk},记ek=xk-x*,若lim𝑘→∞||𝒆𝑘+1||||𝒆𝑘||𝑞=𝑐0,𝑞0(17)其中.||||表示某种范数(对实数可以认为就是绝对值)，则称序列{xk}为q阶收敛.特别地,1阶收敛称为线性收敛,2阶收敛称为平方收敛;若q=1,c=0时(17)式成立,通常称{xk}为超线性收敛.显然,p越大收敛越快.利用p(xk)在x*的泰勒展开:𝑝(𝑥𝑘)=𝑝(𝑥∗)+𝑝′(𝑥∗)(𝑥𝑘−𝑥∗)+⋯+𝑝(𝑞)(𝑥∗)𝑞!(𝑥𝑘−𝑥∗)𝑞+⋯(18)并注意到xk+1=p(xk),x*=p'(x*)和ek=xk-x*,由(18)式可得𝑒𝑘+1=𝑝′(𝑥∗)𝑒𝑘+⋯+𝑝(𝑞)(𝑥∗)𝑞!𝑒𝑘𝑞+⋯(19)于是根据收敛阶的定义(17)式,若p'(x*)≠0,则{xk}为1阶收敛(线性收敛);若p'(x*)=…=p(q-1)(x*)=0,p(q)(x*)≠0,则{xk}为q阶收敛.1.3牛顿法对于方程f(x)=0,将f(x)在xk作泰勒展开，去掉2阶及2阶以上项(即线性化)后得f(x)=f(xk)+f'(z)(x-xk)(20)设f'(x)≠0，令上面的f(x)=0,用xk+1代替右端的x，就得到迭代公式𝑥𝑘+1=𝑥𝑘−𝑓(𝑥𝑘)𝑓′(𝑥𝑘)(21)即迭代函数为𝑝(𝑥)=𝑥−𝑓(𝑥)𝑓′(𝑥)(22)MN是曲线y=f(x)过(xk,f(xk))点的切线，它与x轴的交点即为xk+1.这种方法称为牛顿切线法(简称牛顿法或切线法)，它是线性化与迭代法的结合。由于𝑝′(𝑥∗)=𝑓(𝑥∗)𝑓(𝑥∗)𝑓′(𝑥)2,𝑝(𝑥∗)=𝑓(𝑥∗)𝑓′(𝑥)(23)若x*是f(x)=0的单根，即f(x*)=0,f'(x*)≠0.一般地，f(x*)≠0，则p'(x*)=0，g(x*)≠0，这时牛顿切线法产生的{xn}为2阶收敛.进一步的研究发现，当x*是f(x)=0的重根时，p'(x*)≠0，牛顿切线法只是1阶收敛,并且重数越高收敛越慢.为了避免切线法计算导数的麻烦，可以考虑用差商𝑓(𝑥𝑘)−𝑓(𝑥𝑘−1)𝑥𝑘−𝑥𝑘−1代替f'(xk)，迭代公式变为𝑥𝑘+1=𝑥𝑘−𝑓(𝑥𝑘)(𝑥𝑘−𝑥𝑘−1)𝑓(𝑥𝑘)−𝑓(𝑥𝑘−1)用割线PQ代替了原来的切线，称为割线法(或弦截法)，它的收敛速度比切线法稍慢(可以证明，对于单根其收敛阶数是1.618),并且需要两个初值x0，x1开始迭代.1.4非线性方程组的牛顿法求解单变量非线性方程的牛顿法可以推广到多变量方程组的情形.方程组(2):F(x)=0，其中x=(x1,x2....xn)T,F(x)=(f1(x),f2(x)...fm(x))T，设x(k)=(x(k)1,x(k)2....x(k)n)T是方程组(2)的第k布近似解,与单变量非线性方程的牛顿法类似.在x(k)作泰勒展升,线性化后用x(k+1)代替x可得𝑓𝑖(𝒙(𝑘+1))=𝑓𝑖(𝒙(𝑘))+𝜕𝑓𝑖(𝒙(𝑘))𝜕𝑥1(𝑥1(𝑘+1)−𝑥1(𝑘))+⋯+𝜕𝑓𝑖(𝒙(𝑘))𝜕𝑥𝑛(𝑥𝑛(𝑘+1)−𝑥𝑛(𝑘)),𝑖∈𝑁∗(25)记F的雅可比矩阵为𝑭′(𝒙)=[𝜕𝑓1𝜕𝑥1𝜕𝑓1𝜕𝑥2𝜕𝑓2𝜕𝑥1𝜕𝑓2𝜕𝑥2⋯𝜕𝑓1𝜕𝑥𝑛𝜕𝑓1𝜕𝑥𝑛⋮⋱⋮𝜕𝑓𝑛𝜕𝑥1𝜕𝑓𝑛𝜕𝑥2⋯𝜕𝑓𝑛𝜕𝑥𝑛](26)则(25)式可写作𝑭(𝒙(𝑘+1))=𝑭(𝒙(𝑘))+𝑭′(𝒙(𝑘))(𝒙(𝑘+1)−𝒙(𝑘))(27)若雅可比矩阵可逆，由(27)可得求解线性方程组性方程(25)的牛顿迭代公式𝒙(𝑘+1)=𝒙(𝑘)−[𝑭′(𝒙(𝑘))]−1𝑭(𝒙(𝑘))(28)在实际计算中，在计算过程的第k步，通常是先计算,再解线性方程组𝑭′(𝒙(𝑘))∆𝒙(𝑘)=−𝑭(𝒙(𝑘))(29)得到∆𝒙(𝑘)后,令𝒙(𝑘+1)=𝒙(𝑘)+∆𝒙(𝑘)(30)即可.有定理表明，(28)是超线性收敛的(即收敛阶不小于1)，稍加条件就至少是平方收敛的。四、实验内容1.均相共沸混合物在化学上，确定均相共沸混合物(homogeneousazeotrope)的物质成分(简称为组分)是一项重要而困难的工作，所谓共沸混合物，是指由两种或两种以上物质组成的液体混合物，当在某种压力下被蒸馏或局部汽化时，在气体状态下和在液体状态下保持相同的组分.设该混合物由n个可能的组分组成,组分i所占的比例为xi(i=1,2,...,n),则∑𝑥𝑖𝑛i=1=1,𝑥𝑖≥0(8)均相共沸混合物应该满足稳定条件，即共沸混合物的每个组分在气体状态下和在液体状态下具有相同的化学势能.在压强P不大的情况下，这个条件可以表示为:𝑃=𝛾�