非线性方程组求解及实现

非线性方程组求解及实现第1页，课件共51页，创作于2023年2月复习与练习按以下要求编写一个函数计算

的值，其中x>0时，y=;x<0时，y=2/x;x=0时，返回错误信息（xcann’tbezero)。

要求：1)主函数名称为excer1，x作为输如变量，A作为输出变量；2)主函数中包括一个子函数myfun用于计算y的值。第2页，课件共51页，创作于2023年2月引言在945.36kPa（9.33atm）、300.2K时，容器中充以2mol氮气，试求容器体积。已知此状态下氮气的P-V-T关系符合范德华方程，其范德华常数为a＝4.17atm•L/mol2，b＝0.0371L/mol数学模型：范德华方程变形可得关于V的非线性方程第3页，课件共51页，创作于2023年2月非线性方程（组）在化学计算中的作用多组分混合溶液的沸点、饱和蒸气压计算流体在管道中阻力计算多组分多平衡级分离操作模拟计算平衡常数法求解化学平衡问题定态操作的全混流反应器的操作分析第4页，课件共51页，创作于2023年2月非线性方程非线性方程包括：高次代数方程、超越方程及其它们的组合与线性方程相比，非线性方程求解问题无论从理论上还是从计算公式上都要复杂得多对于高次代数方程，当次数>4时，则没有通解公式可用，对于超越方程既不知有几个根，也没有同样的求解方式。实际上，对于n≥3代数方程以及超越方程都采用数值方法求近似根。第5页，课件共51页，创作于2023年2月非线性方程数值求解原理第6页，课件共51页，创作于2023年2月逐步扫描法逐步扫描法效率较低，常用于求根的初始近似值第7页，课件共51页，创作于2023年2月逐步扫描法计算示例-方程x2-2=0的正数解计算方程的正数解第8页，课件共51页，创作于2023年2月二分法若函数f(x)在区间[a,b]内单调连续，且f(a)f(b)<0，则在闭区间[a,b]内必然存在方程f(x)=0的根x*k=0;whileabs(b-a)>eps x=(a+b)/2;

ifsign(f(x))==sign(f(b)) b=x;

else a=x;

end k=k+1;end二分法的图形解释

二分法的MATLAB程序

二分法是一种可靠的算法，但计算速度较慢第9页，课件共51页，创作于2023年2月二分法计算示例-方程x2-2=0的正数解计算方程的正数解第10页，课件共51页，创作于2023年2月求方程根的精确解非线性方程(组)的求解一般采用迭代法进行。迭代法是一种重要的逐次逼近方法。这种方法用某个固定公式反复校正根的近似值，使之逐步精确化，最后得到满足精度要求的结果。常见的迭代算法有不动点迭代牛顿法弦截法抛物线法威格斯坦法（Wegstein）第11页，课件共51页，创作于2023年2月不动点迭代法

我们可以通过多种方法将方程式转化为例如方程可以转化为以下不同形式(1)(2)(3)第12页，课件共51页，创作于2023年2月不动点迭代法从给定的初值x0，按上式可以得到一个数列：{x0,x1,x2,…,xk,…}如果这个数列有极限，则迭代格式是收敛的。这时数列{xk}的极限就是方程的根上述求非线性代数方程式数值解的方法称为直接迭代法（或称为不动点迭代法）。这个方法虽然简单，但根本问题在于当k->∞时，xk是否收敛于x*，也就是必须找出收敛的充分条件

第13页，课件共51页，创作于2023年2月例题：正确的收敛的迭代格式第14页，课件共51页，创作于2023年2月例题：发散的迭代格式第15页，课件共51页，创作于2023年2月例题：错误的收敛迭代格式第16页，课件共51页，创作于2023年2月不动点定义：函数g(x)的一个不动点（fixedpoint）是指一个实数P，满足P=g(P)从图形角度分析，函数y=g(x)的不动点是y=g(x)和y=x的交点第17页，课件共51页，创作于2023年2月不动点定理

设有(i)g,g’∈C[a,b],(ii)K是一个正常数，(iii)p0∈(a,b)，(iv)对所有x∈[a,b]，有g(x)∈[a,b]如果对于所有x∈[a,b],有|g’(x)|≤K<1，则迭代pn=g(pn-1)将收敛到惟一的不动点P∈[a,b],。在这种情况下，P称为吸引（attractive）不动点。对于所有x∈[a,b],有|g’(x)|>1，则迭代pn=g(pn-1)将不会收敛到P点。在这种情况下，P称为排斥（repelling）不动点，而且迭代显示出局部发散性第18页，课件共51页，创作于2023年2月不动点迭代的图形解释

单调收敛

振荡收敛

第19页，课件共51页，创作于2023年2月不动点迭代的图形解释

单调发散

振荡发散

第20页，课件共51页，创作于2023年2月牛顿法牛顿法也称为牛顿-拉普森法或者切线法。由于这个方法的计算结果颇佳，而计算过程也比较简单，所以被普遍采用。牛顿法的核心内容是通过泰勒级数将非线性方程式转化为线性方程式，然后用迭代法求解。第21页，课件共51页，创作于2023年2月牛顿法原理

设方程式的近似根为则对的泰勒级数展开式为第22页，课件共51页，创作于2023年2月牛顿法的几何意义YOX切线方程

第23页，课件共51页，创作于2023年2月例:牛顿法计算x^2-25=0的解f(x)=x2-25，则f’(x)=2x可构造迭代公式如下：取x0=2代入上式，得x1=7.25，继续递推，依次得5.35、5.0114、5.000001、5.0000000001…第24页，课件共51页，创作于2023年2月牛顿法注意事项在有根区间[a,b]上，连续且不变号，则只要选取的初始近似根x0满足，切线法必定收敛。在单根附近，牛顿公式恒收敛，而且收敛速度很快。但是需要注意如果初始值不在根的附近，牛顿公式不一定收敛在实际使用中，牛顿法最好与逐步扫描法结合起来，先通过逐步扫描法求出根的近似值，然后用牛顿公式求其精确值，以发挥牛顿法收敛速度快的优点牛顿迭代法收敛速度快，但它要求计算函数导数的值第25页，课件共51页，创作于2023年2月弦截法牛顿迭代法收敛速度快，但它要求计算函数导数的值。在科学与工程计算中，常会碰到函数导数不易计算或者算式复杂而不便计算的情况弦截法的基本思想与牛顿法相似，即将非线性函数线性化后求解。两者的差别在于弦截法实现函数线性化的手段采用的是两点间的弦线（用差商代替导数），而不是某点的切线弦截法示意图

第26页，课件共51页，创作于2023年2月弦截法注意事项与牛顿法只需给出一个初值不同，弦截法需要给出两个迭代初值。如果与逐步扫描法结合起来，则最后搜索的区间的两个端点值常可作为初值弦截法虽比牛顿法收敛速度稍慢，但在每次迭代中只需计算一次函数值，又不必求函数的导数，且对初值要求不甚苛刻，是工程计算中常用的有效计算方法之一弦截法虽比牛顿法收敛速度稍慢，但计算量小第27页，课件共51页，创作于2023年2月逆二次插值(IQI)若已知三个点a,b,c，及其函数值f(a),f(b),f(c)，可以将这三点插值为关于y的二次函数。此抛物型一定与x轴有交点，在交点处y=0，对应点x=P(0)为下一步迭代解。IQI法在迭代终点时收敛速度很快，但整个过程中速度不稳定第28页，课件共51页，创作于2023年2月松弛迭代法有些非线性方程用前面的不动点迭代法求解时，迭代过程是发散的。这时可以引入松弛因子，利用松弛迭代法。通过选择合适的松弛因子，就可以使迭代过程收敛迭代法是计算数学的一种重要方法，用途很广，求解线性方程组和矩阵特征值时也要用到这种方法第29页，课件共51页，创作于2023年2月松弛法注意事项由上式可知，当松弛因子ω=1时，松弛迭代法变为不动点迭代法；当松弛因子ω>1时，松弛法使迭代步长加大，可加速迭代，但有可能使原理收敛的迭代变为发散；当0<ω<1时，松弛法使迭代步长减小，这适合于迭代发散或振荡收敛的情况，可使振荡收敛过程加速；当ω<0时，将使迭代反方向进行，可使一些迭代发散过程收敛

松弛迭代法是否有效的关键因素是松弛因子的值能否正确选定。如果值选用适当，能使迭代过程加速，或者使原来不收敛的过程变成收敛；但如果值选用不合适，则效果相反，有时甚至会使原来收敛的过程变得不收敛。松弛因子的数值往往要根据经验选定，但选用较小的松弛因子，一般可以保证迭代过程的收敛第30页，课件共51页，创作于2023年2月威格斯坦法威格斯坦法在化工流程模拟中得到了广泛应用

威格斯坦法是一种迭代加速方法第31页，课件共51页，创作于2023年2月Wegstein法注意事项应注意，如果x1和x2两点选择不当，则连线的斜率等于1，与直线y=x无交点，从而迭代无法进行，这就是Wegstein法应当避免的陷井。引入一个量C第32页，课件共51页，创作于2023年2月Wegstein法注意事项令q＝1-C当q＝0时，Wegstein法退化为简单的不动点迭代当0<q<1时，则变为有阻尼的迭代法。通常q>0时，迭代能稳定收敛，但收敛较慢当q<0可以加速收敛，但易导致不稳定为了加速收敛又避免不稳定，常取-5<q<0，这是称为有界的Wegstein法第33页，课件共51页，创作于2023年2月MATLAB求解非线性方程方法第34页，课件共51页，创作于2023年2月MATLAB求解非线性方程函数非线性方程非线性方程组非线性方程多项式函数rootsfzerofsolve第35页，课件共51页，创作于2023年2月多项式求根函数roots多项式的表达式约定如下：对于多项式，用以下行向量表示：这样就把多项式问题转化为向量问题

Matlab提供了多种多项式计算函数，如多项式求根函数roots，求多项式的值，polyval；多项式乘法，conv；多项式除法，deconv；多项式微分，polyder；多项式拟合，polyfit第36页，课件共51页，创作于2023年2月函数rootsr=roots(c)，用于求解多项式的根其中，行向量c的元素是多项式的系数，按多项式次数降序排列如果c中含有n+1个元素，则多项式为n次roots可以获得多项式的所有根其算法为计算伴随矩阵的特征值第37页，课件共51页，创作于2023年2月例题6：求方程的根>>c=[1-10-1];>>r=roots(c)r=1.4656-0.2328+0.7926i-0.2328-0.7926i>>polyval(c,r(1))ans=-2.5535e-015第38页，课件共51页，创作于2023年2月非线性方程求解函数fzerofzero对于一般的单个超越方程，可以采用fzero函数求解fzero函数结合使用二分法、割线法和可逆二次内插法从两个函数值异号的点a,b开始利用a,b获得割线点c重复以下步骤直至abs(b-a)<*abs(b)或f(b)=0重新排列a,b,c使得f(a)和f(b)异号abs(f(b))<abs(f(a))C替代原来的b如果ca,采用IQI方法如果c=a,采用割线法如果IQI或割线法获得的新区间在[a;b]内，则接受如果IQI或割线法获得的新区间不在[a;b]内，采用区间中点。第39页，课件共51页，创作于2023年2月函数fzero[x,fval,exitflag,output]=fzero(fun,x0,options,p1,p2,...)

此函数的作用求函数fun在x0附件的零值点，x0是标量x 所求解fval 函数在解x处的值exitflag 程序结束情况:>0，程序收敛于解；<0，程序没有收敛；＝0，计算达到了最大次数output 是一个结构体，提供程序运行的信息；output.iterations，迭代次数；output.functions，函数fun的计算次数；output.algorithm，使用的算法options 选项，可用optimset函数设定选项的新值fun可以是函数句柄或匿名函数。Fun函数如何编写？X0如何选取？第40页，课件共51页，创作于2023年2月例题7：计算以下方程的根1)求sinx在3附近的零点;2)求cosx在[1,2]范围内的零点；3)4)本例较简单，可直接在命令窗口输入命令求解：1)fzero(@sin,3)2)fzero(@cos,[1,2])3)fzero(@(x)x^3-2*x-5,1);roots([10-2-5])4)fzero(@(x)x^3-2*sin(x),1)第41页，课件共51页，创作于2023年2月说明1)第2小题中，如果所给区间两端方程不异号，则程序出错2)除了采用匿名函数外，当然可以采用句柄函数定义函数，例如第4小题可以采用如下程序：

functionCha2demo1x=fzero(@fun,1)functiony=fun(x)y=x^3-2*sin(x);3)初值的选择对于解有影响，不同的初值可能获得不同的解可以根据感兴趣的解的区间确定初值范围可以作出函数在一定范围内的曲线，直观的确定解的大致范围4)fzero不能获得多项式的多重根，尤其是复数根。而roots函数求解，则可获得所有根第42页，课件共51页，创作于2023年2月利用函数求解1.直接编写函数functionCha2demo1x=fzero(@fun,1)functiony=fun(x)y=x^3-2*sin(x);2.将上述文件保存为Cha2demo1的m文件3.在命令窗口键入>>Cha2demo1则得到结果functiony=myfun(x)y=x^3-2*sin(x);1.编写函数2.将上述文件保存为myfun的m文件3.在命令窗口中键入>>x=fzero(@fun,1)，则得结果第43页，课件共51页，创作于2023年2月利用函数求解1.编写文件x=fzero(@fun,1)functiony=fun(x)y=x^3-2*sin(x);2.将上述文件保存为excer1的m文件3.在命令窗口键入>>excer1为求解此方程，有人编写如下程序，请问能得到正确结果吗？第44页，课件共51页，创作于2023年2月例题8求的零点，以t为自变量，取值范围为-10<t<10，a,b为参数，本例取值分别为0.1，0.5第45页，课件共51页，创作于2023年2月functionCha2demo2a=0.1;b=0.5;t=-10:0.01:10;Y=sin(t).^2.*exp(-a*t)-b*abs(t);clf,plot(t,Y,'r');holdon;plot(t,zeros(size(t)),'k');xlabel('t');ylabel('y(t)'),holdoffzoomonn=input('Howmanyzeropointsarethere?');[tt,yy]=ginput(n);zoomofffori=1:n[t0(i),y(i),exitflag]=fzero(@(t)sin(t)^2*exp(-a*t)-b*abs(t),tt(i));enddisp('Thezeropointsare:')fprintf('%.4f\t',t0)fprintf('\n')第46页，课件共51页，创作于2023年2月例题9：在945.36kPa（9.33atm）、300.2K时，容器中充以2mol氮气，试求容器体积。已知此状态下氮气的P-V-T关系符合范德华方程，其范德华常数为a＝4.17atm•L/mol2，b＝0.0371L/mol。数学模型：范德华方程变形可得这是关于V的三次方程，可以由roots或fzero求解第47页，课件共51页，创作于2023年2月P=9.33;%atmT=300.2;%Kn=2;%mola=4.17;b=0.0371;R=0.08206;Eq=[P,-(P*n*b+n*R*T),a*n^2,-a*n^3*b];roots(Eq)Script文件函数文件functionCha2demo3P=9.33;%atmT=300.2;%Kn=2;%mola=4.17;b=0.0371;R=0.08206;V0=n*R*T/P%5.2807[V,fval]=fzero(@PVTeq,V0,[],P,T,n,a,b,R)

%-------------------------------------------functionf=PVTeq(V,P,T,n,a,b,R)f=(P+