非线性方程

二分法（dichotomy.m） %非线性方程的二分法

function [x,n]=dichotomy(fun,a,b,eps) if nargin<4 eps=1e-3; end

if feval(fun,a)*feval(fun,b)<0 n=1;c=(a+b)/2; while abs(c)>eps

if feval(fun,a)*feval(fun,c)<0 b=c;c=(a+b)/2;

elseif feval(fun,c)*feval(fun,b)<0 a=c;c=(a+b)/2; else x=c;break; end n=n+1; end x=c;

elseif feval(fun,a)==0 x=a;

elseif feval(fun,b)==0 x=b;

else disp('无法获得方程根!') end

不动点法（unmovepoint.m） %非线性方程的不动点法

function [x,n]=unmovepoint(fun,x0,eps) if nargin<3 eps=1e-3; end

x1=feval(fun,x0); n=1;

while(abs(x1-x0))>=eps

x0=x1;x1=feval(fun,x0); n=n+1;

if n>100000

disp('无法收敛!'); break end end x=x1;

Newton迭代法（newton.m） %非线性方程的Newton迭代法

function [x,n]=newton(fun1,fun2,x0,eps) if nargin<4 eps=1e-3; end

x1=x0-feval(fun1,x0)/feval(fun2,x0); n=1;

while abs(x1-x0)>=eps

x0=x1;x1=x0-feval(fun1,x0)/feval(fun2,x0); n=n+1;

if n>100000

disp('无法收敛!'); break end end x=x1;

割线法（secant.m） %非线性方程的割线法

function [x,n]=secant(fun,x0,eps) if nargin<3 eps=1e-3; end

x1=x0+rand;

x2=x1-feval(fun,x1)*(x1-x0)/(feval(fun,x1)-feval(fun,x0)); n=1;

while abs(x2-x1)>=eps

x0=x1;x1=x2;x2=x1-feval(fun,x1)*(x1-x0)/(feval(fun,x1)-feval(fun,x0)); n=n+1; end x=x2;

用以上四种方法求解下列非线性方程：

fxx2x1

函数：

%非线性方程函数(适用于二分法与割线法) function f=nonlinerequ1(x) f=x^2-x-1;

%非线性方程函数(适用于不动点法) function f=nonlinerequ2(x) f=sqrt(x+1);

%非线性方程函数(适用于Newton迭代法) function f=nonlinerequ3(x) f=x^2-x-1;

%非线性方程函数导数(适用于Newton迭代法) function f=nonlinerequ4(x) f=2*x-1;

命令：

[x,n]=dichotomy(@nonlinerequ1,0,2) [x,n]=dichotomy(@nonlinerequ1,-2,0) [x,n]=unmovepoint(@nonlinerequ2,2)

[x,n]=newton(@nonlinerequ3,@nonlinerequ4,2) [x,n]=newton(@nonlinerequ3,@nonlinerequ4,-2) [x,n]=secant(@nonlinerequ1,2) [x,n]=secant(@nonlinerequ1,-2)

计算结果：（eps=0.001）

迭代方法 解析解

二分法 不动点法 Newton迭代法

割线法

迭代次数n x1

1.61803398874989

50 1.61803398874989

6 1.61835033659881

4 1.61803398874999

5 1.61803400029839

迭代次数n x2

-0.61803398874989 -0.61803398874989 50

迭代出现复数根

-0.61803398874999 5 -0.61803480416329 5

注：不动点法中负根的求解需将函数程序改为-sqrt(x+1)，在计算过程中出现复数根，

且迭代次数多精度低。

苯—甲苯溶液的初浓度为0.4，共100kmol，相对挥发度为2.45，在常压下简单蒸馏至残液量为60kmol，则釜内残液的浓度为多少？

lnx1x21x2W11ln1ln1xW21x21x11 

%非线性方程函数(简单蒸馏) function f=sampledist(x)

w1=100;w2=60;alfa=2.45;x1=0.4;x2=x;

f=log(w1/w2)-1/(alfa-1)*log(x1*(1-x2)/x2/(1-x1))-log((1-x2)/(1-x1));

fzero(@sampledist,0.2)

secant(@sampledist,0.5)……

计算结果：x20.2890

装有3m高填料逆流吸收塔，其气体处理量为0.0237kmol/m2.s，气体入塔浓度为0.05（摩尔分率，下同），要求气体出塔浓度不高于0.005，气相总体积传质系数为0.0545kmol/m3.s，相平衡常数为1.5，现由于吸收剂入塔浓度从零上升为0.001，为了保持吸收率不变，则吸收剂用量应提高为多少？

%非线性方程函数(填料塔) function f=packedtower(x)

h0=3;yb=0.05;ya=0.005;m=1.5;g=0.0237;l=x;kya=0.0545;xa=0.001; s=m*g/l;hog=g/kya;nog=h0/hog;

f=nog-log((1-s)*(yb-m*xa)/(ya-m*xa)+s)/(1-s);

fzero(@packedtower,0.05)

计算结果：0.0426

某三元理想物系，由氯丙烯、二氯丙烷、二氯丙烯组成，其液体混合物的组成及各组分的安托因参数如下：

组分 氯丙烯 二氯丙烷 二氯丙烯

xi 0.0145 0.3090 0.6765 Ai 15.9610 16.0415 18.1025 Bi 2569.00 2985.61 4329.18 Ci -42.2 -52.2 0.0

安托因方程为：lnpi0AiBi，T：K，p：kPa；试求：

CiT （1）在总压为1atm时，该物系的泡点温度；

（2）若该物系在1atm总压下进行闪蒸，并且汽化率为50%，那么需要将其加热至多高的温度？其汽化产生的汽液相组成如何？ 解：

（1）泡点计算函数（satruation.m）

fTyii1Kxiii10，Kipi0p

%非线性方程函数(泡点的计算) function f=saturation(t)

a=[15.9610,16.0415,18.1025]; b=[2569.00,2985.61,4329.18]; c=[-42.2,-52.2,0.0];

x=[0.0145,0.3090,0.6765]; p=101.3;

f=sum(x.*exp(a-b./(c+t))/p)-1; 命令：

secant(@saturation,300) 计算结果： 316.1962K

（2）闪蒸计算函数（flash.m）

物料衡算：FLV，FxFLxVy，平衡关系：yiKixi，Kipi0p 联立得：xi

xii1，

yii1

xFiV1Ki1F，yiKixFiV1Ki1F

%非线性方程函数(闪蒸) function f=flash(t)

a=[15.9610,16.0415,18.1025]; b=[2569.00,2985.61,4329.18]; c=[-42.2,-52.2,0.0];

x=[0.0145,0.3090,0.6765]; p=101.3;rate=0.5; k=exp(a-b./(c+t))/p; x0=x./(1+rate.*(k-1)) y0=k.*x0

f=sum(x0)-sum(y0); 命令：

fzero('flash',300) secant('flash',300) 计算结果： 318.1509K

气相组成：0.0256 液相组成：0.0034

0.3394 0.6350 0.2786 0.7180