典数值方法求解两个非线性方程的根+p218+1

/*

 * ======================================================
 * 本程序使用四种经典数值方法求解两个非线性方程的根：
 *   1. 不动点迭代法 + Steffensen 加速 (Aitken's Δ² process)
 *   2. 牛顿 - 拉夫逊法 (Newton-Raphson Method)
 *   3. 二分法 (Bisection Method)
 *   4. 割线法 (Secant Method)
 *
 * 目标方程:
 *   (1) f(x) = x^2 - 3x + 2 - e^x = 0
 *       预期根：~0.44285
 *   (2) g(x) = x^3 + 2x^2 + 10x - 20 = 0
 *       预期根：~1.36881
 *
 * 输出：所有详细迭代过程将写入 'output.txt' 文件。
 */

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

#define TOL     1e-8        /* 收敛容差：当 |x_k - x_{k-1}| < TOL 时认为收敛 */
#define MAXITER 10000       /* 最大迭代次数：防止算法不收敛时陷入死循环 */

/* ─────────────────────────────────────────────
   函数定义 (原函数及其导数)
   ───────────────────────────────────────────── */

/* 方程 1: f(x) = x^2 - 3x + 2 - e^x */
double f1(double x) { return x*x - 3*x + 2 - exp(x); }
/* 方程 1 的导数: f'(x) = 2x - 3 - e^x (用于牛顿法) */
double df1(double x){ return 2*x - 3 - exp(x); }

/* 方程 2: g(x) = x^3 + 2x^2 + 10x - 20 */
double f2(double x) { return x*x*x + 2*x*x + 10*x - 20; }
/* 方程 2 的导数: g'(x) = 3x^2 + 4x + 10 (用于牛顿法) */
double df2(double x){ return 3*x*x + 4*x + 10; }

/* ─────────────────────────────────────────────
   §1  不动点迭代法 + Steffensen 加速
   核心思想：
     普通不动点迭代 x_{k+1} = g(x_k) 通常只有线性收敛速度。
     Steffensen 方法利用 Aitken Δ² 加速技术，无需导数即可将收敛阶提升至二次。

     加速公式推导：
     令 y = g(x), z = g(y)
     则加速后的新值 xn = x - (y-x)^2 / (z - 2y + x)
   ───────────────────────────────────────────── */

/*
 * 方程 1 的不动点形式重构：
 * 原方程：x^2 - 3x + 2 - e^x = 0  =>  3x = x^2 + 2 - e^x
 * 构造迭代函数：g1(x) = (x^2 + 2 - e^x) / 3
 * 收敛性检查：|g1'(x)| = |(2x - e^x)/3|。在根 ~0.44 处，|g1'| < 1，满足收敛条件。
 */
double g1(double x){ return (x*x + 2 - exp(x)) / 3.0; }

/*
 * 方程 2 的不动点形式重构：
 * 原方程：x^3 + 2x^2 + 10x - 20 = 0  =>  10x = 20 - x^3 - 2x^2
 * 构造迭代函数：g2(x) = (20 - x^3 - 2x^2) / 10
 * 收敛性检查：|g2'(x)| = |(-3x^2 - 4x)/10|。在根 ~1.37 处，|g2'| ≈ 0.96 < 1，收敛较慢，非常适合用 Steffensen 加速。
 */
double g2(double x){ return (20.0 - x*x*x - 2*x*x) / 10.0; }

/**
 * @brief 执行 Steffensen 迭代法
 * @param phi 不动点迭代函数 g(x)
 * @param x0 初始猜测值
 * @param label 方程标签，用于输出标识
 * @param fp 文件指针，用于写入详细日志
 */
void steffensen(double (*phi)(double), double x0,
                const char *label, FILE *fp)
{
    fprintf(fp, "\n--- Steffensen  (%s) ---\n", label);
    fprintf(fp, "  k  |       x_k            |   |x_k - x_{k-1}|\n");
    fprintf(fp, "-----|----------------------|------------------\n");

    double x = x0;
    int k;
    fprintf(fp, " %3d | %20.15f |\n", 0, x);

    for (k = 1; k <= MAXITER; k++) {
        // 计算两步普通不动点迭代值：y = g(x), z = g(g(x))
        double y  = phi(x);
        double z  = phi(y);

        // 计算分母 d = z - 2y + x (即 Aitken 过程中的二阶差分 Δ²x)
        double d  = z - 2*y + x;

        // 【工程鲁棒性】安全性检查：防止分母接近零导致数值溢出或 NaN
        if (fabs(d) < 1e-30) {
            fprintf(fp, "  denominator ≈ 0, switching plain iterate\n");
            x = y; // 退化为一普通不动点迭代，避免崩溃
        } else {
            // Steffensen 加速公式：x_new = x - (Δx)^2 / Δ²x
            double xn = x - (y - x)*(y - x) / d;
            double err = fabs(xn - x);

            fprintf(fp, " %3d | %20.15f | %e\n", k, xn, err);
            x = xn;

            // 收敛判断：误差小于容差
            if (err < TOL) break;
        }
    }
    fprintf(fp, "  Converged in %d iteration(s).  Root ≈ %.15f\n", k, x);
}

/* ─────────────────────────────────────────────
   §2  牛顿法 (Newton's Method)
   核心思想：
     利用泰勒展开的一阶近似，通过切线寻找零点。
     公式：x_{k+1} = x_k - f(x_k) / f'(x_k)
   特点：
     - 二阶收敛（Quadratic convergence），速度极快。
     - 需要计算导数 f'(x)。
     - 对初值敏感，若 f'(x)≈0 或初值离根太远可能发散。
   ───────────────────────────────────────────── */

/**
 * @brief 执行牛顿迭代法
 * @param f 原函数
 * @param df 导函数
 * @param x0 初始猜测值
 * @param label 方程标签
 * @param fp 文件指针
 */
void newton(double (*f)(double), double (*df)(double),
            double x0, const char *label, FILE *fp)
{
    fprintf(fp, "\n--- Newton  (%s) ---\n", label);
    fprintf(fp, "  Initial x0 = %.6f\n", x0);
    fprintf(fp, "  k  |       x_k            |   |x_k - x_{k-1}|\n");
    fprintf(fp, "-----|----------------------|------------------\n");

    double x = x0;
    int k;
    fprintf(fp, " %3d | %20.15f |\n", 0, x);

    for (k = 1; k <= MAXITER; k++) {
        double fx  = f(x);
        double dfx = df(x);

        // 【工程鲁棒性】安全性检查：防止导数接近零导致除零错误
        if (fabs(dfx) < 1e-30) {
            fprintf(fp, "  Derivative ≈ 0! Method fails.\n");
            break;
        }

        // 牛顿迭代公式
        double xn  = x - fx / dfx;
        double err = fabs(xn - x);

        fprintf(fp, " %3d | %20.15f | %e\n", k, xn, err);
        x = xn;

        if (err < TOL) break;
    }
    fprintf(fp, "  Converged in %d iteration(s).  Root ≈ %.15f\n", k, x);
}

/* ─────────────────────────────────────────────
   §3  二分法 (Bisection Method)
   核心思想：
     基于介值定理。若 f(a)*f(b) < 0，则 [a,b] 间必有根。
     每次迭代将区间减半。
   特点：
     - 线性收敛（Linear convergence），速度较慢。
     - 只要初始区间包含根且函数连续，保证收敛（Robust）。
     - 不需要导数。
   ───────────────────────────────────────────── */

/**
 * @brief 执行二分法
 * @param f 原函数
 * @param a 区间左端点
 * @param b 区间右端点
 * @param label 方程标签
 * @param fp 文件指针
 */
void bisection(double (*f)(double), double a, double b,
               const char *label, FILE *fp)
{
    fprintf(fp, "\n--- Bisection  (%s) ---\n", label);
    fprintf(fp, "  k  |       a              |       b              |       mid            | f(mid)\n");
    fprintf(fp, "-----|----------------------|----------------------|----------------------|----------\n");

    // 检查初始区间是否有效（异号）
    if (f(a)*f(b) > 0) {
        fprintf(fp, "  Error: f(a) and f(b) have the same sign! No guarantee of root.\n");
        return;
    }

    double mid = a;
    int k;
    for (k = 1; k <= MAXITER; k++) {
        mid = (a + b) / 2.0;
        double fm = f(mid);

        fprintf(fp, " %3d | %20.15f | %20.15f | %20.15f | %+.3e\n",
                k, a, b, mid, fm);

        // 收敛判断：区间长度小于容差
        if (fabs(b - a) < TOL) break;

        // 根据符号更新区间：根位于函数值异号的子区间内
        if (f(a)*fm < 0) b = mid; // 根在 [a, mid]
        else             a = mid; // 根在 [mid, b]
    }
    fprintf(fp, "  Converged in %d iteration(s).  Root ≈ %.15f\n", k, mid);
}

/* ─────────────────────────────────────────────
   §4  割线法 (Secant Method)
   核心思想：
     用割线斜率近似代替牛顿法中的导数，避免直接求导。
     公式：x_{k+1} = x_k - f(x_k) * (x_k - x_{k-1}) / (f(x_k) - f(x_{k-1}))
   特点：
     - 超线性收敛（Superlinear, order ≈ 1.618，黄金分割比）。
     - 不需要显式计算导数。
     - 需要两个初始值。
   ───────────────────────────────────────────── */

/**
 * @brief 执行割线法
 * @param f 原函数
 * @param x0 第一个初始猜测
 * @param x1 第二个初始猜测
 * @param label 方程标签
 * @param fp 文件指针
 */
void secant(double (*f)(double), double x0, double x1,
            const char *label, FILE *fp)
{
    fprintf(fp, "\n--- Secant  (%s) ---\n", label);
    fprintf(fp, "  k  |       x_k            |   |x_k - x_{k-1}|\n");
    fprintf(fp, "-----|----------------------|------------------\n");
    fprintf(fp, " %3d | %20.15f |\n", 0, x0);
    fprintf(fp, " %3d | %20.15f |\n", 1, x1);

    double xp = x0, xc = x1; // xp = x_{k-1}, xc = x_k
    int k;
    for (k = 2; k <= MAXITER; k++) {
        double fp_ = f(xp), fc = f(xc);

        // 【工程鲁棒性】安全性检查：防止分母 (f(xc) - f(xp)) 接近零
        if (fabs(fc - fp_) < 1e-30) {
            fprintf(fp, "  Denominator ≈ 0! Method fails.\n");
            break;
        }

        // 割线法迭代公式
        double xn  = xc - fc*(xc - xp)/(fc - fp_);
        double err = fabs(xn - xc);

        fprintf(fp, " %3d | %20.15f | %e\n", k, xn, err);

        // 更新前两个点，准备下一轮迭代
        xp = xc;
        xc = xn;

        if (err < TOL) break;
    }
    // 注意：循环从 k=2 开始，实际有效迭代次数为 k-1
    fprintf(fp, "  Converged in %d iteration(s).  Root ≈ %.15f\n", k-1, xc);
}

/* ─────────────────────────────────────────────
   §5  辅助函数 (无输出的快速运行版本)
   用于生成最终的对比表格，不包含繁琐的文件打印逻辑，以提高统计效率
   ───────────────────────────────────────────── */
typedef struct { int iters; double root; } Result;

// 快速运行 Steffensen 并返回结果
Result run_steffensen(double (*phi)(double), double x0){
    double x = x0; int k;
    for(k=1;k<=MAXITER;k++){
        double y=phi(x),z=phi(y),d=z-2*y+x;
        if(fabs(d)<1e-30){x=y;continue;}
        double xn=x-(y-x)*(y-x)/d;
        if(fabs(xn-x)<TOL){x=xn;break;}
        x=xn;
    }
    return (Result){k,x};
}

// 快速运行 Newton 并返回结果
Result run_newton(double (*f)(double),double (*df)(double),double x0){
    double x=x0; int k;
    for(k=1;k<=MAXITER;k++){
        double fx=f(x),dfx=df(x);
        if(fabs(dfx)<1e-30) break;
        double xn=x-fx/dfx;
        if(fabs(xn-x)<TOL){x=xn;break;}
        x=xn;
    }
    return (Result){k,x};
}

// 快速运行 Bisection 并返回结果
Result run_bisection(double (*f)(double),double a,double b){
    double mid=a; int k;
    // 简单假设输入区间有效，实际工程中应加检查
    for(k=1;k<=MAXITER;k++){
        mid=(a+b)/2.0;
        if(fabs(b-a)<TOL) break;
        if(f(a)*f(mid)<0) b=mid; else a=mid;
    }
    return (Result){k,mid};
}

// 快速运行 Secant 并返回结果
Result run_secant(double (*f)(double),double x0,double x1){
    double xp=x0,xc=x1; int k;
    for(k=2;k<=MAXITER;k++){
        double fp_=f(xp),fc=f(xc);
        if(fabs(fc-fp_)<1e-30) break;
        double xn=xc-fc*(xc-xp)/(fc-fp_);
        if(fabs(xn-xc)<TOL){xc=xn;break;}
        xp=xc;xc=xn;
    }
    return (Result){k-1,xc};
}

/* ════════════════════════════════════════════ */
int main(void)
{
    // 打开输出文件
    FILE *fp = fopen("output.txt", "w");
    if (!fp) { perror("fopen"); return 1; }

    // 打印报告头
    fprintf(fp, "╔══════════════════════════════════════════════════════════╗\n");
    fprintf(fp, "║   Numerical Analysis – Root Finding: Full Results        ║\n");
    fprintf(fp, "║   Tolerance: 1e-8  (|x_k – x_{k-1}| < 1e-8)              ║\n");
    fprintf(fp, "╚══════════════════════════════════════════════════════════╝\n");

    /* ── §1 Steffensen ────────────────────────────────────────────── */
    fprintf(fp, "\n══════════════════════════════════════════════\n");
    fprintf(fp, "  §1  Fixed-Point Iteration + Steffensen\n");
    fprintf(fp, "══════════════════════════════════════════════\n");
    fprintf(fp, "Eq1: g1(x) = (x^2 + 2 - e^x)/3,  x0 = 0.0\n");
    steffensen(g1, 0.0, "Eq1", fp);
    fprintf(fp, "\nEq2: g2(x) = (20 - x^3 - 2x^2)/10,  x0 = 1.5\n");
    steffensen(g2, 1.5, "Eq2", fp);

    /* ── §2 Newton ────────────────────────────────────────────────── */
    fprintf(fp, "\n══════════════════════════════════════════════\n");
    fprintf(fp, "  §2  Newton's Method\n");
    fprintf(fp, "══════════════════════════════════════════════\n");
    newton(f1, df1, 0.0, "Eq1, x0=0.0", fp);
    newton(f2, df2, 1.5, "Eq2, x0=1.5", fp);

    /* ── §3 Bisection ─────────────────────────────────────────────── */
    fprintf(fp, "\n══════════════════════════════════════════════\n");
    fprintf(fp, "  §3  Bisection Method\n");
    fprintf(fp, "══════════════════════════════════════════════\n");
    bisection(f1, 0.0, 1.0, "Eq1, [0,1]",  fp);
    bisection(f2, 1.0, 2.0, "Eq2, [1,2]",  fp);

    /* ── §4 Secant ────────────────────────────────────────────────── */
    fprintf(fp, "\n══════════════════════════════════════════════\n");
    fprintf(fp, "  §4  Secant Method\n");
    fprintf(fp, "══════════════════════════════════════════════\n");
    secant(f1, 0.0, 1.0,  "Eq1, x0=0, x1=1",     fp);
    secant(f2, 1.0, 1.5,  "Eq2, x0=1, x1=1.5",   fp);

    /* ── §5 Comparison Table ──────────────────────────────────────── */
    fprintf(fp, "\n══════════════════════════════════════════════\n");
    fprintf(fp, "  §5  Performance Comparison Summary\n");
    fprintf(fp, "══════════════════════════════════════════════\n");

    Result r; // 未使用变量，可移除
    fprintf(fp, "\n┌────────────────────────────┬──────────────┬──────────────┬────────────────────────┐\n");
    fprintf(fp, "│ Method                     │  Eq1 iters   │  Eq2 iters   │  Root (Eq1 / Eq2)      │\n");
    fprintf(fp, "├────────────────────────────┼──────────────┼──────────────┼────────────────────────┤\n");

    // 运行所有方法并收集结果以生成对比表
    Result s1 = run_steffensen(g1, 0.0);
    Result s2 = run_steffensen(g2, 1.5);
    fprintf(fp, "│ Steffensen                 │ %12d │ %12d │ %10.8f / %10.8f │\n",
            s1.iters, s2.iters, s1.root, s2.root);

    Result n1 = run_newton(f1,df1,0.0);
    Result n2 = run_newton(f2,df2,1.5);
    fprintf(fp, "│ Newton                     │ %12d │ %12d │ %10.8f / %10.8f │\n",
            n1.iters, n2.iters, n1.root, n2.root);

    Result b1 = run_bisection(f1,0.0,1.0);
    Result b2 = run_bisection(f2,1.0,2.0);
    fprintf(fp, "│ Bisection                  │ %12d │ %12d │ %10.8f / %10.8f │\n",
            b1.iters, b2.iters, b1.root, b2.root);

    Result sc1 = run_secant(f1,0.0,1.0);
    Result sc2 = run_secant(f2,1.0,1.5);
    fprintf(fp, "│ Secant                     │ %12d │ %12d │ %10.8f / %10.8f │\n",
            sc1.iters, sc2.iters, sc1.root, sc2.root);

    fprintf(fp, "└────────────────────────────┴──────────────┴──────────────┴────────────────────────┘\n");

    // 输出理论分析部分：收敛阶与特性
    fprintf(fp, "\n──────────────────────────────────────────────────────────────\n");
    fprintf(fp, "  Theoretical Convergence Order & Properties\n");
    fprintf(fp, "──────────────────────────────────────────────────────────────\n");
    fprintf(fp, "\n  Method         | Order      | Needs deriv? | Bracket? | Notes\n");
    fprintf(fp, "  ---------------|------------|--------------|----------|-----------------------------\n");
    fprintf(fp, "  Bisection      | linear(1)  | No           | Yes      | Always converges; slow\n");
    fprintf(fp, "  Fixed-Point    | linear(1)  | No           | No       | Requires |g'|<1 near root\n");
    fprintf(fp, "  Steffensen     | quadratic  | No           | No       | Aitken Δ² acceleration\n");
    fprintf(fp, "  Secant         | ~1.618     | No           | No       | Superlinear; may diverge\n");
    fprintf(fp, "  Newton         | quadratic  | Yes (f')     | No       | Fastest; needs good x0\n");

    // 输出适用场景建议
    fprintf(fp, "\n──────────────────────────────────────────────────────────────\n");
    fprintf(fp, "  Applicable Scenarios\n");
    fprintf(fp, "──────────────────────────────────────────────────────────────\n");
    fprintf(fp, "\n  Bisection:\n");
    fprintf(fp, "    Best when a bracketing interval is known and the function is\n");
    fprintf(fp, "    continuous.  Guaranteed to converge but requires ~27 iterations\n");
    fprintf(fp, "    per decimal digit (log2(1/tol) ≈ 27 for tol=1e-8).\n");
    fprintf(fp, "\n  Fixed-Point + Steffensen:\n");
    fprintf(fp, "    Useful when a convergent iteration form g(x) can be constructed\n");
    fprintf(fp, "    easily.  Steffensen's Aitken Δ² step elevates order to quadratic\n");
    fprintf(fp, "    at virtually no extra cost; no derivative required.\n");
    fprintf(fp, "\n  Newton:\n");
    fprintf(fp, "    Ideal when f'(x) is cheap to evaluate and a good initial guess\n");
    fprintf(fp, "    exists.  Quadratic convergence means very few iterations.\n");
    fprintf(fp, "    May fail if f'(x0)≈0 or x0 is far from the root.\n");
    fprintf(fp, "\n  Secant:\n");
    fprintf(fp, "    A practical alternative to Newton when the derivative is\n");
    fprintf(fp, "    expensive or unavailable.  Superlinear (order≈1.618) convergence\n");
    fprintf(fp, "    with only two function evaluations per step.\n");

    fprintf(fp, "\n══════════════════════════════════════════════\n");
    fprintf(fp, "  End of Report\n");
    fprintf(fp, "══════════════════════════════════════════════\n");

    fclose(fp);
    printf("Done – results written to output.txt\n");
    return 0;
}