拟牛顿算法c语言代码实现

下面是一个简单的拟牛顿算法的C语言代码实现，供参考：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 2  // 变量个数
#define MAX_ITER 1000  // 最大迭代次数
#define TOL 1e-5  // 精度

// 二次函数
double fun(double x[N])
{
    return 100 * pow(x[1] - pow(x[0], 2), 2) + pow(1 - x[0], 2);
}

// 二次函数梯度
void gradient(double x[N], double g[N])
{
    g[0] = -400 * x[0] * (x[1] - pow(x[0], 2)) - 2 * (1 - x[0]);
    g[1] = 200 * (x[1] - pow(x[0], 2));
}

int main()
{
    double x[N] = {0.0};  // 初始点
    double g[N] = {0.0};  // 梯度
    double h[N][N] = {{1, 0}, {0, 1}};  // 初始Hessian矩阵
    double s[N] = {0.0};  // 搜索方向
    double y[N] = {0.0};  // 梯度差
    double a = 1.0;  // 步长
    double f = 0.0;  // 函数值
    double fold = 0.0;  // 上一次迭代的函数值
    double alpha = 1.0;  // 黄金分割法中的参数
    int iter = 0;  // 迭代次数
    int i, j;

    // 计算初始点的函数值和梯度
    f = fun(x);
    gradient(x, g);

    while (iter < MAX_ITER && fabs(f - fold) > TOL) {
        iter++;
        fold = f;

        // 计算搜索方向s
        for (i = 0; i < N; i++) {
            s[i] = 0.0;
            for (j = 0; j < N; j++) {
                s[i] -= h[i][j] * g[j];
            }
        }

        // 一维搜索
        a = 1.0;
        while (1) {
            for (i = 0; i < N; i++) {
                y[i] = fun(x) - fun(x + a * s);
            }
            if (fun(x + a * s) <= fun(x) - 0.5 * a * dot(g, s, N)) {
                break;
            }
            a *= 0.5;
        }

        // 更新x和Hessian矩阵
        for (i = 0; i < N; i++) {
            y[i] = g[i] - y[i] / a;
            for (j = 0; j < N; j++) {
                h[i][j] += y[i] * s[j] / dot(y, s, N) - h[i][j] * dot(y, s, N) * dot(y, s, N);
            }
            x[i] += a * s[i];
            g[i] = y[i];
        }

        f = fun(x);  // 更新函数值
    }

    printf("Minimum found at (%.4f, %.4f)\n", x[0], x[1]);
    printf("Minimum value is: %.4f\n", f);

    return 0;
}

// 计算向量a和b的点积
double dot(double a[], double b[], int n)
{
    double res = 0.0;
    int i;
    for (i = 0; i < n; i++) {
        res += a[i] * b[i];
    }
    return res;
}

这段代码实现了一个简单的拟牛顿算法，用于求解二次函数的最小值。其中，`fun()`函数计算函数值，`gradient()`函数计算梯度，`dot()`函数计算向量的点积。在主函数中，首先初始化各个变量，然后进行迭代，直到满足精度要求或达到最大迭代次数为止。每次迭代中，计算搜索方向并进行一维搜索，然后更新变量和Hessian矩阵，最后计算新的函数值。