C语言实现Ujevic法求非线性方程根返回-nan问题咨询

问题现象 #

• 实现目标：用Ujevic迭代法求解非线性方程 $f(x)=x^3-3x+2$ 的实根，测试用例理论根值为1，预期迭代次数75次
• 测试输入：区间端点a=0.5、b=1.2，初始迭代值x0=0.5，松弛参数Alpha=0.1
• 异常表现：程序仅迭代43次就输出根结果为-nan，与预期结果完全不符

原始问题代码 #

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

#define f(x) ((x*x*x)-(3*x)+2)
#define g(x) ((3*x*x)-3)

void main() {
    float a, b, x0, xk, z0, xt, y;
    int i = 0;
    printf("\nEnter The Value Of a: ");
    scanf("%f", &a);
    printf("\nEnter The Value Of b: ");
    scanf("%f", &b);
    printf("\nEnter The Value Of x0: ");
    scanf("%f", &x0);
    printf("\nEnter The Value Of Alpha: ");
    scanf("%f", &y);
    do {
        xt = x0;
        z0 = x0 - y * (f(x0) / g(x0));
        xk = x0 + 4 * (z0 - x0) * (f(x0) / ((3 * f(x0)) - (2 * f(z0))));
        x0 = xk;
        i++;
    } while (fabs(xk - xt >= 0.000000000000000001));
    printf("The value of root is: %f", xk);
    printf("\n the no of iterations is %d", i);
}

实际运行输出 #

Enter The Value Of a: 0.5
Enter The Value Of b: 1.2
Enter The Value Of x0: 0.5
Enter The Value Of Alpha: 0.1
The value of root is: -nan
the no of iterations is 43

问题核心诱因 #

1. 循环判断逻辑完全写错
   收敛判断的括号位置错误，当前写法fabs(xk - xt >= 0.000000000000000001)的执行逻辑是：先判断xk-xt >= 1e-18得到布尔值0（假）或1（真），再对这个0/1值取绝对值，根本没有计算两次迭代的差值大小，循环不会在精度达标时正常退出。直到迭代中出现NaN，而NaN参与所有比较运算都返回假，此时布尔判断结果为0，循环才异常终止，刚好对应43次迭代就退出的现象。
2. 无除零防护触发非法运算
   迭代公式中有两处除法运算：f(x0)/g(x0)和f(x0)/(3*f(x0)-2*f(z0))，代码完全没有做分母趋近于0的判断。本次测试的函数在x=1处是二重根，导数g(1)=0，迭代靠近根时分母会无限趋近于0，直接做除法会得到无穷大值，后续迭代运算就会产生NaN。
3. 数据精度与阈值不匹配
   所有迭代变量用的是float单精度类型，有效位数仅6~7位，代码设置的1e-18收敛阈值远小于单精度能表示的最小精度，阈值完全无效。迭代到单精度精度上限后数值会无规则震荡，进一步放大除零、产生NaN的概率。另外代码中读入的a、b区间参数完全没被使用，没有做初始区间的根存在性校验。
4. 语法不规范
   标准C语言要求main函数返回值为int类型，写void main属于未定义行为，不同编译器下可能出现异常。

修正方案 #

• 修正收敛判断的括号位置，将fabs的作用范围限定在两次迭代的差值上，收敛阈值匹配所用数据类型的实际精度范围。
• 两处除法运算前增加分母判0逻辑，分母绝对值小于极小值（如1e-12）时终止迭代，避免非法运算。
• 将迭代相关变量替换为double双精度类型，提升计算精度，可支持更小的收敛阈值，匹配预期75次迭代的测试要求。
• 增加最大迭代次数限制，避免迭代不收敛时出现死循环；补全区间校验逻辑，确认输入区间内存在实根再启动迭代。
• 修正宏定义的参数括号、main函数返回值等语法问题。

修正后参考代码 #

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

// 宏定义参数加括号避免运算优先级问题
#define f(x) ((x)*(x)*(x) - 3*(x) + 2)
#define g(x) (3*(x)*(x) - 3)
#define EPS 1e-12
#define CONV_THRESHOLD 1e-10
#define MAX_ITER 1000

int main() {
    double a, b, x0, xk, z0, xt, alpha;
    int i = 0;
    printf("\nEnter The Value Of a: ");
    scanf("%lf", &a);
    printf("\nEnter The Value Of b: ");
    scanf("%lf", &b);
    printf("\nEnter The Value Of x0: ");
    scanf("%lf", &x0);
    printf("\nEnter The Value Of Alpha: ");
    scanf("%lf", &alpha);

    // 初始区间根校验
    if (f(a)*f(b) > 0) {
        printf("输入区间端点函数值同号，无法确认区间内存实根\n");
        return 1;
    }

    do {
        xt = x0;
        double gx0 = g(x0);
        // 导数位置除零防护
        if (fabs(gx0) < EPS) {
            printf("迭代到导数趋近于0的重根位置，终止迭代\n");
            break;
        }
        z0 = x0 - alpha * (f(x0) / gx0);
        double denom = 3*f(x0) - 2*f(z0);
        // 迭代分母除零防护
        if (fabs(denom) < EPS) {
            printf("迭代过程分母趋近于0，终止迭代\n");
            break;
        }
        xk = x0 + 4 * (z0 - x0) * (f(x0) / denom);
        x0 = xk;
        i++;
        // 最大迭代次数防护
        if (i >= MAX_ITER) {
            printf("达到最大迭代次数，未收敛\n");
            break;
        }
    } while (fabs(xk - xt) >= CONV_THRESHOLD); // 修正fabs括号位置

    printf("The value of root is: %lf", xk);
    printf("\n the no of iterations is %d", i);
    return 0;
}

内容的提问来源于stack exchange，提问作者Rajarshi Chanda