非线性方程可以有两种形式：

求根问题：给定一个 $R^n$ 到 $R^n$ 的映射 (mapping) $f$ ，找到一个 $n$ 维向量 $x$ 作为 $f$ 的根，满足

$$ f(x) = 0. $$

不动点问题：给定一个 $R^n$ 到 $R^n$ 的映射 (mapping) $g$，计算一个 $n$ 维向量 $x$ 作为 $g$ 的不动点，满足

$$ x = g(x). $$

在经济学问题中，我们经常要对 $x$ 做出限制，比如限价（农产品最低价格、最低工资）、限量（配额），这时我们要处理互补问题 (complementarity problem)，一般表示为：给定两个 $n$ 维向量 $a$ 和 $b$ ($a < b$) 和函数 $f: R^n \mapsto R^n$, 求一个 $n$ 维向量 $x \in [a, b]$ 满足

$$ x_i > a_i \implies f_i(x) \ge 0 \quad \forall i = 1, \ldots, n,\newline x_, < b_i \implies f_i(x) \le 0 \quad \forall i = 1, \ldots, n. $$

Bisection Method 链接到标题

对于非线性方程组求解，我们经常用的数值方法是将非线性问题转化为一系列线性问题。其中最简单也是最可靠的方法是二分法，其原理是首先找到一个区间，函数 $f$ 在区间的两端取相反符号，那么根据中值定理，这个区间内一定有解。下一步将这个区间二等分，其中一个小区间的两端一定可以令函数的符号相反，然后在这个小区间内重复之前的步骤，直到区间范围缩小到容许范围内。

下面我给出该算法的 julia 实现：

function bisect(f, a, b, tol=1.5e-8)
    s = sign(f(a))
    x = (a + b) / 2  # the middle point
    d = (b-a) / 2    # range
    while d > tol    # when range > tolerance, iterate
        d = d / 2
        if s == sign(f(x))  # sign of middle equals sign of beginning
            x = x + d       # then the root is in the latter half
        else
            x = x - d       # otherwise the root is in the former half
        end
    end
    return x   # when range <= tolerance, return the middle point
end

Function Iteration 链接到标题

对于不动点问题的求解，我们有类似的迭代方法——函数迭代。根据这个方法，我们首先猜测一个 $x^{(0)}$ 作为迭代的起始，带入函数 $g$ 则有

$$ x^{(k+1)} \gets g(x^{(k)}). $$

如果 $g$ 是连续函数，且迭代收敛，那么迭代一定会收敛到 $g$ 的不动点。同样这里我也给出一个 julia 的实现：

function fixpoint(g, x0, tol=1.5e-8)
    x = x0
    x_new = g(x)
    while abs(x - x_new) > tol
        x = x_new
        x_new = g(x)
    end
    return x_new
end

Newton’s Method 链接到标题

在实践中，非线性方程组一般通过牛顿法求根，这一方法通过泰勒展开把难以处理的非线性方程转化为一系列线性方程。根据这一方法，首先我们猜测一个解作为迭代的起点，然后在该点我们对函数 $f$ 进行一阶泰勒展开，这个一阶泰勒展开式的解会更接近真实解，以此为基础重复前一步骤，不断迭代直到泰勒展开式的解收敛。

因此我们要不断计算一阶泰勒展开式的解：

$$ f(x) \approx f(x^{(k)}) + f’(x^{(k)}) (x - x^{(k)}) = 0, $$

由此我们可以得出迭代规则：

$$ x^{(k+1)} \gets x^{(k)} - {[f’(x^{(k)})]}^{-1} f(x^{(k)}). $$

下面是 julia 的实现：

using LinearAlgebra: norm
function newton(f, x, maxit = 1000, tol = 1.5e-8)
    for _ in 1:maxit
        fval, fjac = f(x)
        x = x - fjac \ fval
        if norm(fval) < tol
            return x
        end
    end
end

其中函数 f(x) 应该返回函数 $f$ 在 x 处的值 fval 和 Jacobian 值 fjac。

Quasi-Newton Methods 链接到标题

牛顿方法的一个缺陷是需要我们自己提供 Jacobian 值，类牛顿方法通过把 Jacobian $f’$ 换为一个通用的式子，来降低我们求解出错的可能性，当然这一方法会降低迭代收敛的速度。

下面我来介绍两种 quasi-newton methods: Secant Method 和 Broyden’s Method。

其中，Secant Method 把 $f$ 的倒数表示为近似值

$$ f’(x^{(k)}) \approx \frac{f(x^{(k)}) - f(x^{(k-1)})} {x^{(k)} - x^{(k-1)}}. $$

这一方法需要提供两个初始值来计算 $f’$。

而 Broyden’s Method 把 Secant Method 从一元函数扩展到多元函数。使用这个方法的时候，我们需要提供初始值 $x^{(0)}$ 和 Jacobian 矩阵 $A^{(0)}$。给定 $x^{(k)}$ 和 $A^{(k)}$ ，我们能解出函数 $f$ 一阶泰勒展开的根

$$ f(x) \approx f(x^{(k)}) + A^{(k)} (x - x^{(k)}) = 0, $$

因此得到迭代规则

$$ x^{(k+1)} \gets x^{(k)} - {(A^{(k)})}^{-1} f(x^{(k)}). $$

类似于导数的定义，我们可以得到 $A^{(k+1)}$：

$$ f(x^{(k+1)}) - f(x^{(k)} ) = A^{(k+1)} (x^{(k+1)} - x^{(k)} ) \newline \implies A^{(k+1)} \gets A^{(k)} + \left( f(x^{(k+1)}) - f(x^{(k)}) - A^{(k)} d^{(k)} \right) \frac{d^{(k)\top}}{d^{(k)\top} d^{(k)}}, $$

其中 $d^{(k)} = x^{(k+1)} - x^{(k)}$。

写到这里，我比较累了，对于 Broyden 方法，我就不给出 julia 的实现方式了，对于用户来说，我们直接使用 NLsolve.jl 包中的 nlsolve 函数就行。