非线性方程组计算

在科学与工程计算中，经常遇到求解非线性方程组的问题；非线性方程组在收敛速度及收敛性比线性方程组要差，特别对于非凸的非线性方程组，其求解更是困难。下面简要介绍非线性方程组的三种解法——牛顿法、拟牛顿法、同伦算法，分析三种解法的适用性，并附Matlab 原程序。

（一）、牛顿迭代法

迭代公式为：x k+1=x k-f(x k)/f'(x k);牛顿迭代法是解非线性方程组比较经典的方法，在局部收敛点附近是平方收敛的；但其解依赖于初始解，且迭代每一步都要计算f'(x k)，不仅计算量大而且有时会发生计算困难。

（二）、拟牛顿迭代法

拟牛顿法是为了解决求Jacobi矩阵时带来的困难，现已成为解决非线性方程组和最优化问题的最有效方法之一。其迭代格式为：

x(k+1)=x(k)-A k-1F(x(k))

A k+1=A k+[(y k-A k s k)(y k-A k s k)T]/[(y k-A k s k)T s k]

在一定条件下，计算H的序列是超收敛的，但稳定性较差，有时迭代效果不理想。

（三）、同伦算法

同伦算法基本思想是从容易求解的方程组开始，逐步过渡到原方程组的求解，从而得到问题的解。非线性方程组为：F(x)=0，其解为X*。

构造泛函 G：[0，1]XR n->R n

G定义为：G（λ，x）=λ F(x)+(1-λ)[F(x)-F(x(0))]=F(x)+(λ-1)F(x(0))

(其中：x（0）为任意给的初值，假定为λ函数（λ=0）)

对于λ的方程G（λ，x）=0，

当λ=0时，0=G（0，x）=F(x)-F(x(0))；x（0）是方程的解；

当λ=1时，0=G（1，x）=F(x)；x*是方程的解，即x(1)=x*

基于这个思想我们最后可以得到如下关系式：

x'(λ)=-[J(x(λ))]-1F(x(0)) ( 0<=λ<=1,对初始值x(0) )

J为雅可比矩阵，由上面的式子，对λ在[0，1]上积分，就可得到x*=x(1)

上面的非线性方程组问题就转化为数值积分问题。

附：同伦算法的Matlab程序

%% N 为积分部数，h 为积分步长，其值视计算精度而定%% x0 为初始值，可任意给定

%% Equfun(x)为非线性方程组表达式（列向量），

%% Jacobfun(x)为雅可比矩阵

N=100;

h=1/N;

x0=[0 0 0]';

x=x0;

%format long

f=Equfun(x);

b=-h*f;

for i=1:N

A=Jacobfun(x);

k1=inv(A)*b;

A=Jacobfun(x+0.5*k1);

k2=inv(A)*b;

A=Jacobfun(x+0.5*k2);

k3=inv(A)*b;

A=Jacobfun(x+0.5*k3);

k4=inv(A)*b;

x=x+(k1+2*k2+2*k3+k4)/6;

end

disp('The Solution is:')

disp('x=');disp(x);

%%%%%%%%%%%%%%%%%%%%%%%%%

%%%下面是非线性方程组及其雅可比行列式%%%%%%%%%%%%%%%%%%%%%%%%%

function y=Equfun(x)

%nonliner functions

y=[3*x(1)-cos(x(2)*x(3))-1/2;

x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;

exp(-x(1)*x(2))+20*x(3)+(10*pi-3)/3]; %%%%%%%%%%%%%%%%%%%%%%%%%

function y=Jacobfun(x)

%Jacobi function

%J=zeros(length(x));

y(1,1)=3;

y(1,2)=x(3)*sin(x(2)*x(3));

y(1,3)=x(2)*sin(x(2)*x(3));

y(2,1)=2*x(1);

y(2,2)=-162*(x(2)+0.1);

y(2,3)=cos(x(3));

y(3,1)=-x(2)*exp(-x(1)*x(2));

y(3,2)=-x(1)*exp(-x(1)*x(2));

y(3,3)=20;