非线性方程组数值求解

Posted by baoduge on 2017年7月1日

非线性系统的数值求解是个大坑，今天先开个头，之后慢慢填。今天先讨论不含偏微分的非线性方程组的数值解法。
假设有一个很简单的非线性方程组
$\alpha (x+1)(x-3) e^{\beta y}=0$
$\lambda (y-5)(y+7) e^{\gamma x}=0$
我们需要对其求解。其中 $\alpha$ , $\beta$ , $\lambda$ , $\gamma$ 均为参数，为了适用于一些需要调参的情况，我们依然在方程中保留这四个参数。
我们知道如果 $\alpha \neq 0$ 且 $\lambda \neq 0$ ，那么我们可以有四个解：
$(x,y)=(-1,5)$
$(x,y)=(-1,-7)$
$(x,y)=(3,5)$
$(x,y)=(3,-7)$
因此数值求解的时候，初始条件非常重要。

1.MATLAB

MATLAB中有三种方法：fsolve，fmincon和lsqnonlin。

可以将上述非线性方程组表达成一个function如下：

function F = myfun(var,para)
alpha=para(1);
beta=para(2);
lambda=para(3);
gamma=para(4);
x=var(1);
y=var(2);
F(1) = alpha*(x+1)*(x-3)*exp(beta*y);
F(2) = lambda*(y-5)*(y+7)*exp(gamma*x);
return

function F = myfun(var,para)

alpha=para(1);

beta=para(2);

lambda=para(3);

gamma=para(4);

x=var(1);

y=var(2);

F(1) = alpha*(x+1)*(x-3)*exp(beta*y);

F(2) = lambda*(y-5)*(y+7)*exp(gamma*x);

return

在此 $\alpha$ , $\beta$ , $\lambda$ , $\gamma$ 除了作为输入参数传入myfun还可以使用global命令。在此 $\alpha$ , $\beta$ , $\lambda$ , $\gamma$ 都是数值型数据，所以para是一个vector。当遇到需要传递数值型和字符型数据混合的参数时，可以使用cell封装。

1.1.fsolve

fsolve:
Solves a problem specified by
F(x) = 0
for x, where F(x) is a function that returns a vector value.
x is a vector or a matrix; see Matrix Arguments.

fsolve是MATLAB中最常用的求解非线性方程组的方法，MATLAB中提供三种算法trust-region-dogleg,trust-region和levenberg-marquardt

options = optimoptions('fsolve');
options.Algorithm = 'trust-region-dogleg';
%   'trust-region-dogleg' (default)
%   'trust-region-reflective'
%   'levenberg-marquardt'
options.MaxFunEvals = Inf;
options.MaxIter = 500;
options.Display = 'iter';
para = [1,3,1，3];
F = @(var)myfun(var,para);
var0=[10,10];
[fs,~,exitflag] = fsolve(F,var0,options);
ceq=myfun(fs,para);

options = optimoptions('fsolve');

options.Algorithm = 'trust-region-dogleg';

% 'trust-region-dogleg' (default)

% 'trust-region-reflective'

% 'levenberg-marquardt'

options.MaxFunEvals = Inf;

options.MaxIter = 500;

options.Display = 'iter';

para = [1,3,1，3];

F = @(var)myfun(var,para);

var0=[10,10];

[fs,~,exitflag] = fsolve(F,var0,options);

ceq=myfun(fs,para);

计算结果如下：

&gt;&gt;fs
fs =
    3     5
&gt;&gt; exitflag
exitflag =
    1
&gt;&gt; ceq
ceq =
   1.0e-09 *
    0.7167    0.1599

>>fs

fs =

3 5

>> exitflag

exitflag =

>> ceq

ceq =

1.0e-09 *

0.7167 0.1599

1.2.fmincon

算法有
interior-point
Trust Region Reflective Algorithm
sqp【2,3】
sqp-legacy
active-set

1.3.lsqnonlin

trust-region-reflective
levenberg-marquardt

2.Fortran

IMSL Numerical Library中的NEQNF方法。

Solves a system of nonlinear equations using a modified Powell hybrid
algorithm and a finite-difference approximation to the Jacobian.

Powell's Method
有限差分估计（finite-difference approximation）去近似雅各比矩阵（Jacobian Matrix）

3.Python

scipy.optimize中的fsolve

fsolve is a wrapper around MINPACK’s hybrd and hybrj algorithms.

MINIPACK的HYBRD

HYBRD is a modification of the Powell hybrid method. Two of its
main characteristics involve the choice of the correction as a
convex combination of the Newton and scaled gradient directions,
and the updating of the Jacobian by the rank-1 method of Broy-
den.

4.R

nleqslv package

The nleqslv package provides two algorithms for solving (dense) nonlinear systems of equations.
The methods provided are
• a Broyden Secant method where the matrix of derivatives is updated after each major iteration
using the Broyden rank 1 update.
• a full Newton method where the Jacobian matrix of derivatives is recalculated at each iteration

Broyden's method