非线性方程组求解及matlab实现.pptx

《非线性方程组求解及matlab实现.pptx》由会员分享，可在线阅读，更多相关《非线性方程组求解及matlab实现.pptx（54页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、复习与练习复习与练习按以下要求编写一个函数计算 的值，其中x0时，y=;x4时，则没有通解公式可用，对于超越方程既不知有几个根，也没有同样的求解方式。实际上，对于n3代数方程以及超越方程都采用数值方法求近似根。第4页/共54页非线性方程数值求解原理非线性方程数值求解原理第5页/共54页逐步扫描法逐步扫描法逐步扫描法效率较低，常用于求根的初始近似值第6页/共54页逐步扫描法计算示例逐步扫描法计算示例-方程x2-2=0的正数解计算方程 的正数解第7页/共54页二分法二分法n若函数f(x)在区间a,b内单调连续，且f(a)f(b)epsx=(a+b)/2;if sign(f(x)=sign(f(b)

2、b=x;elsea=x;endk=k+1;end二分法的图形解释二分法的图形解释 二分法的二分法的MATLAB程序程序 二分法是一种可靠的算法，但计算速度较慢第8页/共54页二分法计算示例二分法计算示例-方程x2-2=0的正数解计算方程 的正数解第9页/共54页求方程根的精确解求方程根的精确解非线性方程(组)的求解一般采用迭代法进行。迭代法迭代法是一种重要的逐次逼近方法。这种方法用某个固定公式反复校正根的近似值，使之逐步精确化，最后得到满足精度要求的结果。常见的迭代算法有不动点迭代牛顿法弦截法抛物线法 威格斯坦法（Wegstein）第10页/共54页不动点迭代法不动点迭代法 我们可以通过多种方

3、法将方程式转化为例如方程可以转化为以下不同形式(1)(2)(3)第11页/共54页不动点迭代法不动点迭代法从给定的初值x0，按上式可以得到一个数列：x0,x1,x2,xk,如果这个数列有极限，则迭代格式是收敛的。这时数列xk的极限 就是方程的根 上述求非线性代数方程式数值解的方法称为直接迭代法（或称为不动点迭代法）。这个方法虽然简单，但根本问题在于当k-时，xk是否收敛于x*，也就是必须找出收敛的充分条件 第12页/共54页例题例题：正确的收敛的迭代格式第13页/共54页例题例题：发散的迭代格式第14页/共54页例题例题：错误的收敛迭代格式第15页/共54页不动点不动点定义定义：函数g(x)的

4、一个不动点（fixed point）是指一个实数P，满足P=g(P)从图形角度分析，函数y=g(x)的不动点是y=g(x)和y=x的交点 第16页/共54页不动点定理不动点定理 设有(i)g,g Ca,b,(ii)K是一个正常数，(iii)p0(a,b)，(iv)对所有x a,b，有g(x)a,b如果对于所有x a,b,有|g(x)|K1，则迭代pn=g(pn-1)将不会收敛到P点。在这种情况下，P称为排斥（repelling）不动点，而且迭代显示出局部发散性 第17页/共54页不动点迭代的图形解释不动点迭代的图形解释 单调收敛单调收敛 振荡收敛振荡收敛 第18页/共54页不动点迭代的图形解释

5、不动点迭代的图形解释 单调发散单调发散 振荡发散振荡发散 第19页/共54页牛顿法牛顿法牛顿法也称为牛顿-拉普森法或者切线法。由于这个方法的计算结果颇佳，而计算过程也比较简单，所以被普遍采用。牛顿法的核心内容是通过泰勒级数将非线性方程式转化为线性方程式，然后用迭代法求解。第20页/共54页牛顿法原理牛顿法原理 设方程式 的近似根为则 对 的泰勒级数展开式为第21页/共54页牛顿法的几何意义牛顿法的几何意义YOX切线方程切线方程 第22页/共54页例例:牛顿法计算牛顿法计算x2-25=0的解的解f(x)=x2-25，则f(x)=2x可构造迭代公式如下：取x0=2代入上式，得x1=7.25，继续递

6、推，依次得5.35、5.0114、5.000001、5.0000000001 第23页/共54页牛顿法注意事项牛顿法注意事项在有根区间a,b上，连续且不变号，则只要选取的初始近似根x0满足 ，切线法必定收敛。在单根附近，牛顿公式恒收敛，而且收敛速度很快。但是需要注意如果初始值不在根的附近，牛顿公式不一定收敛在实际使用中，牛顿法最好与逐步扫描法结合起来，先通过逐步扫描法求出根的近似值，然后用牛顿公式求其精确值，以发挥牛顿法收敛速度快的优点 牛顿迭代法收敛速度快，但它要求计算函数导数的值第24页/共54页弦截法弦截法牛顿迭代法收敛速度快，但它要求计算函数导数的值。在科学与工程计算中，常会碰到函数导

7、数不易计算或者算式复杂而不便计算的情况 弦截法的基本思想与牛顿法相似，即将非线性函数线性化后求解。两者的差别在于弦截法实现函数线性化的手段采用的是两点间的弦线（用差商代替导数），而不是某点的切线 弦截法示意图弦截法示意图 第25页/共54页弦截法注意事项弦截法注意事项与牛顿法只需给出一个初值不同，弦截法需要给出两个迭代初值。如果与逐步扫描法结合起来，则最后搜索的区间的两个端点值常可作为初值弦截法虽比牛顿法收敛速度稍慢，但在每次迭代中只需计算一次函数值，又不必求函数的导数，且对初值要求不甚苛刻，是工程计算中常用的有效计算方法之一 弦截法虽比牛顿法收敛速度稍慢，但计算量小第26页/共54页逆二次插

8、值逆二次插值(IQI)若已知三个点a,b,c，及其函数值f(a),f(b),f(c)，可以将这三点插值为关于y的二次函数。此抛物型一定与x轴有交点，在交点处y=0，对应点x=P(0)为下一步迭代解。IQI法在迭代终点时收敛速度很快，但整个过程中速度不稳定第27页/共54页松弛迭代法松弛迭代法有些非线性方程用前面的不动点迭代法求解时，迭代过程是发散的。这时可以引入松弛因子，利用松弛迭代法。通过选择合适的松弛因子，就可以使迭代过程收敛 迭代法是计算数学的一种重要方法，用途很广，求解线性方程组和矩阵特征值时也要用到这种方法第28页/共54页松弛法注意事项松弛法注意事项由上式可知，当松弛因子=1时，松

9、弛迭代法变为不动点迭代法；当松弛因子1时，松弛法使迭代步长加大，可加速迭代，但有可能使原理收敛的迭代变为发散；当01时，松弛法使迭代步长减小，这适合于迭代发散或振荡收敛的情况，可使振荡收敛过程加速；当0时，将使迭代反方向进行，可使一些迭代发散过程收敛 松弛迭代法是否有效的关键因素是松弛因子的值能否正确选定。如果值选用适当，能使迭代过程加速，或者使原来不收敛的过程变成收敛；但如果值选用不合适，则效果相反，有时甚至会使原来收敛的过程变得不收敛。松弛因子的数值往往要根据经验选定，但选用较小的松弛因子，一般可以保证迭代过程的收敛第29页/共54页威格斯坦法威格斯坦法威格斯坦法在化工流程模拟中得到了广泛

10、应用威格斯坦法在化工流程模拟中得到了广泛应用 威格斯坦法是一种迭代加速方法 第30页/共54页Wegstein法注意事项法注意事项应注意，如果x1和x2两点选择不当，则连线的斜率等于1，与直线y=x无交点，从而迭代无法进行，这就是Wegstein法应当避免的陷井。引入一个量C第31页/共54页Wegstein法注意事项法注意事项令q1-C1.当q0时，Wegstein法退化为简单的不动点迭代2.当0q0时，迭代能稳定收敛，但收敛较慢3.当q0可以加速收敛，但易导致不稳定4.为了加速收敛又避免不稳定，常取-5qc=1-1 0-1;r=roots(c)r=1.4656 -0.2328+0.7926

11、i -0.2328-0.7926ipolyval(c,r(1)ans=-2.5535e-015 第37页/共54页非线性方程求解函数非线性方程求解函数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)0，程序收敛于解；Cha2demo1则得到结果function y=myfun(x)y=x3-2*sin(x);1.编写函数2.将上述文件保存为myf

12、un的m文件3.在命令窗口中键入 x=fzero(fun,1)，则得结果第42页/共54页利用函数求解利用函数求解1.编写文件x=fzero(fun,1)function y=fun(x)y=x3-2*sin(x);2.将上述文件保存为excer1的m文件3.在命令窗口键入 excer1为求解此方程，有人编写如下程序，请问能得到正确结果吗？第43页/共54页例题例题8求 的零点，以t为自变量，取值范围为-10t10，a,b为参数，本例取值分别为0.1，0.5第44页/共54页function Cha2demo2a=0.1;b=0.5;t=-10:0.01:10;Y=sin(t).2.*exp(

13、-a*t)-b*abs(t);clf,plot(t,Y,r);hold on;plot(t,zeros(size(t),k);xlabel(t);ylabel(y(t),hold offzoom onn=input(How many zero points are there?);tt,yy=ginput(n);zoom offfor i=1:n t0(i),y(i),exitflag=fzero(t)sin(t)2*exp(-a*t)-b*abs(t),tt(i);enddisp(The zero points are:)fprintf(%.4ft,t0)fprintf(n)第45页/共54

14、页例题例题9：在945.36kPa（9.33atm）、300.2K时，容器中充以2mol氮气，试求容器体积。已知此状态下氮气的P-V-T关系符合范德华方程，其范德华常数为a4.17atmL/mol2，b0.0371L/mol。数学模型：范德华方程变形可得这是关于V的三次方程，可以由roots或fzero求解 第46页/共54页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*n2,-a*n3*b;roots(Eq)Script文件函数文件function Cha2demo3P=9.33;%a

15、tmT=300.2;%Kn=2;%mola=4.17;b=0.0371;R=0.08206;V0=n*R*T/P%5.2807V,fval=fzero(PVTeq,V0,P,T,n,a,b,R)%-function f=PVTeq(V,P,T,n,a,b,R)f=(P+a*n2/V2)*(V-n*b)-n*R*T;函数文件第47页/共54页函数函数fsolve 与fzero函数只能求解单个方程的根不同，fsolve函数可求解非线性方程组的解。其算法采用的是最小二乘法。调用格式：调用格式：x,fval,exitflag,output,jacobian=fsolve(fun,x0,options,

16、p1,p2,.)输入输出变量的意义同fzero函数输出变量中的jacobian为函数fun在x处的Jacobian矩阵定义待求解方程时，fun函数表示F(X)=0的形式！第48页/共54页例题例题10：在命令窗口输入：x0=1 1 1;x=fsolve(fun,x0)function y=fun(x)y(1)=sin(x(1)+x(2)2+log(x(3)-7;y(2)=3*x(1)+2x(2)-x(3)3+1;y(3)=x(1)+x(2)+x(3)-5;x=0.5991 2.3959 2.0050第49页/共54页fsolve函数的应用函数的应用在铜管内在在铜管内在1 atm下将异丙醇加热到

17、下将异丙醇加热到400。已知铜是生产丙酮和丙。已知铜是生产丙酮和丙醛的催化剂，或许还有某些异丙醇异构化为正丙醇。这三种产物的生醛的催化剂，或许还有某些异丙醇异构化为正丙醇。这三种产物的生成可用如下三个独立反应表示：成可用如下三个独立反应表示：iC3H7OH(IP)=n C3H7OH(NP)K1=0.064iC3H7OH(IP)=(CH3)CO(AC)+H2 K2=0.076iC3H7OH(IP)=CHHO(PR)+H2 K3=0.00012后续加工步骤需要正丙醇，虽然可含丙酮，但丙醛含量不能超过后续加工步骤需要正丙醇，虽然可含丙酮，但丙醛含量不能超过0.05(mol)%。在上述反应条件下，是否

18、存在违反这种规定的可能性。在上述反应条件下，是否存在违反这种规定的可能性？数学模型：各反应的化学平衡方程如下数学模型：各反应的化学平衡方程如下第50页/共54页function Cha2demo7x0=0.05 0.2 0.01;x=fsolve(EquiC3,x0);CAC=x(3)/sum(x)if CAC0.05 disp(The AC concentration could not be over 0.05%)else disp(The AC concentration could be over 0.05%)endfunction f=EquiC3(x)f1=x(1)-0.064*(1-x(1)-x(2)-x(3);f2=x(2)*(x(2)+x(3)-0.076*(1-x(1)-x(2)-x(3)*(1+x(2)+x(3);f3=x(3)*(x(2)+x(3)-0.00012*(1-x(1)-x(2)-x(3)*(1+x(2)+x(3);f=f1 f2 f3;fsolve函数的应用第51页/共54页练习练习求解：初始值(-5,-5)1.2.第52页/共54页复习复习不动点迭代及其收敛性roots、fzero、fsolve表示待求解方程的函数编写方法与规定函数的调用方法初值的选取原则第53页/共54页感谢您的观看！第54页/共54页