非线性方程求解

关于非线性方程求解课堂练习某蒸馏釜的操作压力位106.7kpa，其中溶液含苯摩尔分数为0.2，甲苯为0.8，求此溶液的泡点和平衡的气相组成。苯甲苯溶液可以看做理想溶液，组分蒸汽压为：a为苯，b为甲苯第2页,共48页，星期六，2024年，5月数值计算的误差误差的分类模型误差观测误差截断误差舍入误差计算机对所储存的数据位数有限制计算时只截取有限项第3页,共48页，星期六，2024年，5月误差的分类绝对误差相对误差定义：设X*是数X的近似值，如果X*的绝对误差限时他的某一位数位的半个单位，并且从X*左起第一个非零数字到该数位共有N为，则称这个N个数字为X*的有效数字，也称X*近似X时有N为有效数字有效数字第4页,共48页，星期六，2024年，5月算法设计的原则避免两个数量级相差很大的数运算算法中应尽量减少运算次数避免两个相近的数相减递推运算中，防止误差的积累避免用绝对值对过小的数作除数第5页,共48页，星期六，2024年，5月非线性方程求解非线性方程数值求解基本原理Matlab求解非线性方程函数－roots，fzero，fsolve第6页,共48页，星期六，2024年，5月引言在945.36kPa（9.33atm）、300.2K时，容器中充以2mol氮气，试求容器体积。已知此状态下氮气的P-V-T关系符合范德华方程，其范德华常数为a＝4.17atm•L/mol2，b＝0.0371L/mol数学模型：范德华方程变形可得关于V的非线性方程第7页,共48页，星期六，2024年，5月非线性方程非线性方程包括：高次代数方程、超越方程及其它们的组合与线性方程相比，非线性方程求解问题无论从理论上还是从计算公式上都要复杂得多对于高次代数方程，当次数>4时，则没有通解公式可用，对于超越方程既不知有几个根，也没有同样的求解方式。实际上，对于n≥3代数方程以及超越方程都采用数值方法求近似根。第8页,共48页，星期六，2024年，5月非线性方程数值求解原理第9页,共48页，星期六，2024年，5月逐步扫描法逐步扫描法效率较低，常用于求根的初始近似值开始结束输入初值A，步长hX1=A，Y1=f(X1)X2=X1+hY2=f(X2)打印A，h，X1Y1*Y2>0X1=X2Y1=Y2是否第10页,共48页，星期六，2024年，5月寻找根区间functionmasearch(fun,a,b,h)n=(b-a)/h;x1=zeros(1,n);x2=zeros(1,n);a1=a;b1=a1+h;k=1;while(b1<b)iffeval(fun,a1)*feval(fun,b1)<0x1(k)=a1;x2(k)=b1;elsea1=b1;b1=a1+h;continue;enda1=b1;b1=a1+h;k=k+1;endfori=1:kifx1(i)-x2(i)~=0[x1(i),x2(i)]endend第11页,共48页，星期六，2024年，5月求方程根的精确解非线性方程(组)的求解一般采用迭代法进行。迭代法是一种重要的逐次逼近方法。这种方法用某个固定公式反复校正根的近似值，使之逐步精确化，最后得到满足精度要求的结果。常见的算法有二分法不动点迭代牛顿法第12页,共48页，星期六，2024年，5月二分法若函数f(x)在区间[a,b]内单调连续，且f(a)f(b)<0，则在闭区间[a,b]内必然存在方程f(x)=0的根x*functionx=mabisec(fun,a,b,ep)x=(a+b)/2.0;k=0;whileabs(feval(fun,x))>ep|(b-a>ep)iffeval(fun,x)*feval(fun,a)<0b=x;elsea=x;endx=(a+b)/2.0;k=k+1;end二分法的图形解释

二分法的MATLAB程序

二分法是一种可靠的算法，但计算速度较慢第13页,共48页，星期六，2024年，5月二分法计算示例计算方程的正数解第14页,共48页，星期六，2024年，5月不动点迭代法

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

第16页,共48页，星期六，2024年，5月例题：正确的收敛的迭代格式第17页,共48页，星期六，2024年，5月例题：发散的迭代格式第18页,共48页，星期六，2024年，5月不动点定义：函数g(x)的一个不动点（fixedpoint）是指一个实数x*，满足x*=g(x*)从图形角度分析，函数y=g(x)的不动点是y=g(x)和y=x的交点第19页,共48页，星期六，2024年，5月不动点定理

设函数P(x)在区间[a,b]上有连续的一阶导数，且满足：则函数P(x)在区间[a,b]上存在唯一的不动点X*，即有x*=P(x*)，且对任意的属于区间[a,b]的x0，都可以通过迭代得到X*。第20页,共48页，星期六，2024年，5月不动点迭代的图形解释

单调收敛

振荡收敛

第21页,共48页，星期六，2024年，5月不动点迭代的图形解释

单调发散

振荡发散

第22页,共48页，星期六，2024年，5月迭代法functionx=maiter(phi,x0,ep,N)

ifnargin<4N=500;end

ifnargin<3ep=1e-4;end

k=0;

whilek<N

x=feval(phi,x0);

ifabs(x-x0)<epbreak;

end

x0=x;k=k+1;

end

ifk==N,warning(’已达迭代次数上限’);

enddisp([’k=’,num2str(k)])第23页,共48页，星期六，2024年，5月牛顿法牛顿法也称为牛顿-拉普森法或者切线法。由于这个方法的计算结果颇佳，而计算过程也比较简单，所以被普遍采用。牛顿法的核心内容是通过泰勒级数将非线性方程式转化为线性方程式，然后用迭代法求解。第24页,共48页，星期六，2024年，5月牛顿法原理

设方程式的近似根为则对的泰勒级数展开式为第25页,共48页，星期六，2024年，5月牛顿法的几何意义YOX切线方程

第26页,共48页，星期六，2024年，5月例:牛顿法计算x^2-25=0的解f(x)=x2-25，则f’(x)=2x可构造迭代公式如下：取x0=2代入上式，得x1=7.25，继续递推，依次得5.35、5.0114、5.000001、5.0000000001…第27页,共48页，星期六，2024年，5月牛顿法functionx=maqnewt(fun,x0,x1,ep,N)ifnargin<5,N=500;endifnargin<4,ep=1e-4;endk=0;whilek<Nx=x1-(x1-x0)*feval(fun,x1)/(feval(fun,x1)-feval(fun,x0));ifabs(x-x1)<epbreak;endx0=x1;x1=x;k=k+1;endifk==N,warning(’已达迭代次数上限’);enddisp([’k=’,num2str(k)])第28页,共48页，星期六，2024年，5月牛顿法注意事项在有根区间[a,b]上，连续且不变号，则只要选取的初始近似根x0满足，切线法必定收敛。在单根附近，牛顿公式恒收敛，而且收敛速度很快。但是需要注意如果初始值不在根的附近，牛顿公式不一定收敛在实际使用中，牛顿法最好与逐步扫描法结合起来，先通过逐步扫描法求出根的近似值，然后用牛顿公式求其精确值，以发挥牛顿法收敛速度快的优点牛顿迭代法收敛速度快，但它要求计算函数导数的值第29页,共48页，星期六，2024年，5月MATLAB求解非线性方程方法第30页,共48页，星期六，2024年，5月MATLAB求解非线性方程函数非线性方程非线性方程组非线性方程多项式函数rootsfzerofsolve第31页,共48页，星期六，2024年，5月多项式求根函数roots多项式的表达式约定如下：对于多项式，用以下行向量表示：这样就把多项式问题转化为向量问题

Matlab提供了多种多项式计算函数，如多项式求根函数roots，求多项式的值，polyval；多项式乘法，conv；多项式除法，deconv；多项式微分，polyder；多项式拟合，polyfit第32页,共48页，星期六，2024年，5月函数rootsr=roots(c)，用于求解多项式的根其中，行向量c的元素是多项式的系数，按多项式次数降序排列如果c中含有n+1个元素，则多项式为n次roots可以获得多项式的所有根第33页,共48页，星期六，2024年，5月例题：求方程的根>>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第34页,共48页，星期六，2024年，5月非线性方程求解函数fzerofzero对于一般的单个超越方程，可以采用fzero函数求解fzero函数结合使用二分法、割线法和可逆二次内插法第35页,共48页，星期六，2024年，5月函数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如何选取？第36页,共48页，星期六，2024年，5月例题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)第37页,共48页，星期六，2024年，5月说明除了采用匿名函数外，当然可以采用句柄函数定义函数functionCha2demo1x=fzero(@fun,1)functiony=fun(x)y=x^3-2*sin(x);2)初值的选择对于解有影响，不同的初值可能获得不同的解可以根据感兴趣的解的区间确定初值范围可以作出函数在一定范围内的曲线，直观的确定解的大致范围3)fzero不能获得多项式的多重根，尤其是复数根。而roots函数求解，则可获得所有根第38页,共48页，星期六，2024年，5月例题：在945.36kPa（9.33atm）、300.2K时，容器中充以2mol氮气，试求容器体积。已知此状态下氮气的P-V-T关系符合范德华方程，其范德华常数为a＝4.17atm•L/mol2，b＝0.0371L/mol。数学模型：范德华方程变形可得这是关于V的三次方程，可以由roots或fzero求解第39页,共48页，星期六，2024年，5月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+a*n^2/V^2)*(V-n*b)-n*R*T;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+a*n^2/V^2)*(V-n*b)-n*R*T;函数文件第40页,共48页，星期六，2024年，5月函数fsolve

与fzero函数只能求解单个方程的根不同，fsolve函数可求解非线性方程组的解。其算法采用的是最小二乘法。调用格式：[x,fval,exitflag,output,jacobian]=fsolve(fun,x0,options,p1,p2,...)输入输出变量的意义同fzero函数输出变量中的jacobian为函数fun在x处的Jacobian矩阵定义待求解方程时，必须首先将方程组变换成F(X)=0的形式！

第41页,共48页，星期六，2024年，5月例题10：

在命令窗口输入：x0=[111];x=fsolve(@fun,x0)functiony=fun(x)y(1)=sin(x(1))+x(2)^2+log(x(3))-7;y(2)=3*x(1)+2^x(2)-x(3)^3+1;y(3)=x(1)+x(2)+x(3)-5;x=0.59912.39592.0050第42页,共48页，星期六，2024年，5月fsolve函数的应用在铜管内在1atm下将异丙醇加热到400℃。已知铜是生产丙酮和丙醛的催化剂，或许还有某些异丙醇异构化为正丙醇。这三种产物的生成可用如下三个独立反应表示：iC3H7OH(IP)=nC3H7OH(NP)K1=0.064iC3H7OH(IP)=(CH3)CO(AC)+H2K2=0.076iC3H7OH(IP)=C２H５ＣHO(PR)+H2K3=0.00012后续加工步骤需要正丙醇，虽然可含丙酮，但丙醛含量不能超过0.05(mol)%。在上述反应条件下，是否存在违反这种规定的可能性？数学模型：各反应的化学平衡方程如下第43页,共48页，星期六，2024年，5月functionCha2demo7x0=[0.050.20.01];x=fsolve(@EquiC3,x0);CAC=x(3)/sum(x)ifCAC<0.05disp('TheACconcentrationcouldnotb