非线性方程组迭代解法不动点法（unmovepoints.m）%非线性方程组的不动点法function [x,n]=unmovepoints(fun,x0,eps)if nargin<3eps=1e-3;endx1=feval(fun,x0);n=1;while(norm(x1-x0))>=epsx0=x1;x1=feval(fun,x0);n=n+1;if n>100000disp('无法收敛!');breakendendx=x1;Newton迭代法（newtons.m）%非线性方程组的Newton迭代法function [x,n]=newtons(fun1,fun2,x0,eps)if nargin<4eps=1e-3;endx1=x0-feval(fun1,x0)/feval(fun2,x0);n=1;while norm(x1-x0)>=epsx0=x1;x1=x0-feval(fun1,x0)/feval(fun2,x0);n=n+1;if n>100000disp('无法收敛!');breakendendx=x1;注：方程组的迭代与方程迭代不同之处在于收敛的判断不能用abs 而应用norm （范数，默认值为向量各元素的平方和的开方；norm(x1-x0)即为向量x1与x0对应元素差的平方和的开方。
在对应的函数程序中应注意向量的运算与数量运算的区别。
）用以上方法求解下列非线性方程组：()0cos 2.0sin 7.02111=--=x x x X f ()0sin 2.0cos 7.02122=+-=x x x X f函数：%非线性方程组函数(适用于不动点法) function f=nonlinerequs1(x)f(1)=0.7*sin(x(1))+0.2*cos(x(2)); f(2)=0.7*cos(x(1))-0.2*sin(x(2));%非线性方程组函数(适用于Newton 迭代法) function f=nonlinerequs2(x)f(1)=x(1)-0.7*sin(x(1))-0.2*cos(x(2)); f(2)=x(2)-0.7*cos(x(1))+0.2*sin(x(2));%非线性方程组函数导数(适用于Newton 迭代法) function f=nonlinerequs3(x)f=[1-0.7*cos(x(1)),0.2*sin(x(2));0.7*sin(x(1)),1+0.2*cos(x(2))];导数为⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂22122111x f x f x f xf ，对多方程则类似。
命令：fsolve(@nonlinerequs2,[0.5,0.5])[x,n]=unmovepoints(@nonlinerequs1,[0,0],1e-6)[x,n]=newtons(@nonlinerequs2,@nonlinerequs3,[0,0],1e-6)计算结果：（eps=0.000001）迭代方法X 迭代次数n解析解[0.52652262191818 0.50791971903685] - fsolve [0.52652266171295 0.50791973020932]- 不动点法[0.52652130091388 0.50792028463452] 30 Newton 迭代法[0.52652279369020 0.50791961189450] 16在某操作条件下，有如下四个独立的反应：B A ⇔ DC A +⇔DE A +⇔ CF A +⇔其平衡常数分别为：0.08，0.06，0.001，0.05；反应前只有组分A ，没有其他物质，试求反应平衡时组分A 的摩尔分率为多少？ 解：设反应前组分A 的总摩尔数为1，反应平衡后四反应过程分别消耗组分A 的摩尔数为1x 、2x 、3x 、4x ，所以反应平衡时各组分的摩尔数为：A ：43211x x x x ----B ：1xC ：42x x +D ：32x x +E ：3xF ：4x 故有：()011432111=-----=K x x x x x X f ()()()()()0112432432132422=-+++----++=K x x x x x x x x x x x X f ()()()()011343243213233=-+++----+=K x x x x x x x x x x X f()()()()011443243214244=-+++----+=K x x x x x x x x x x X freaction.m%非线性方程组函数(化学平衡,适用于fsolve) function f=reaction(x)f(1)=x(1)/(1-x(1)-x(2)-x(3)-x(4))-0.08;f(2)=(x(2)+x(4))*(x(2)+x(3))/(1-x(1)-x(2)-x(3)-x(4))/(1+x(2)+x(3)+x(4))-0.06; f(3)=x(3)*(x(2)+x(3))/(1-x(1)-x(2)-x(3)-x(4))/(1+x(2)+x(3)+x(4))-0.001; f(4)=x(4)*(x(2)+x(4))/(1-x(1)-x(2)-x(3)-x(4))/(1+x(2)+x(3)+x(4))-0.05; 命令：x=fsolve(@reaction,[0.3,0.2,0.1,0])或x=fsolve(@reaction,[0.3,0.2,0.1,0],foptions) 计算结果：x=[0.0514, 0.1621, 0.0050, 0.1392] 因此组分A 的摩尔分率为：4917.0114324321=+++----x x x x x x xreaction1.m%非线性方程组函数(化学平衡,适用于不动点法) function f=reaction1(x)f(1)=(1-x(1)-x(2)-x(3)-x(4))*0.08;f(2)=(1-x(1)-x(2)-x(3)-x(4))*(1+x(2)+x(3)+x(4))*0.06/(x(2)+x(3))-x(4); f(3)=(1-x(1)-x(2)-x(3)-x(4))*(1+x(2)+x(3)+x(4))*0.001/(x(2)+x(3)); f(4)=(1-x(1)-x(2)-x(3)-x(4))*(1+x(2)+x(3)+x(4))*0.05/(x(2)+x(4)); 其迭代格式为：()4321111x x x x K x ----=()()()43243243212211x x x x x x x x x x K x -++++----= ()()()3243243213311x x x x x x x x x K x ++++----=()()()4243243214411x x x x x x x x x K x ++++----=但该迭代格式无法收敛，所以不同的迭代格式是否收敛有很大不同；要满足迭代收敛则必须满足以下关系：迭代格式：()()()k k X X φ=+1 迭代函数的偏导数矩阵：⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂n n n n n n x x x x x x x x x φφφφφφφφφ112221212111 迭代格式收敛的充分条件：111max <∑=≤≤nj ij ni a 或111ma x<∑=≤≤ni ij nj a 或11,2<∑=nj i ij a故在迭代格式的构造上比较困难，在实际应用中若能满足以上条件是最好的；但如果无法构造满足以上条件的迭代格式，那么还可采取如下的方法。
在不动点法中引如松弛因子1<W ，将原迭代格式改为：()()()k k X X φ=+1 ⇒ ()()()()()[]k k k k X X W X X -+=+φ1unmovepoints.m%非线性方程组的不动点法function [x,n]=unmovepoints(fun,x0,eps)if nargin<3eps=1e-3;endx1=feval(fun,x0);n=1;while(norm(x1-x0))>=epsx0=x1;x1=x0+0.1*(feval(fun,x0)-x0);无法收敛时,采用该式迭代,其中0.1为松弛因子<1,可根据具体情况修改%x1=feval(fun,x0);n=n+1;if n>100000disp('无法收敛!');breakendendx=x1;命令：（松弛因子取0.1，其数值不同迭代收敛程度不同，不收敛则降低松弛因子）x=unmovepoints(@reaction1,[0.3,0.2,0.1,0],1e-6)计算结果：x=[0.0514, 0.1621, 0.0050, 0.1392]。