割线法求解sinc类超越方程部分测试用例精度异常的排查修复

问题原因 #

1. 单精度浮点数精度不足：代码中使用float类型，其有效位数仅6-7位，无法支撑部分测试用例对1e-7~1e-9精度的要求，计算过程中累积误差会导致结果偏离目标值。
2. 迭代终止后返回错误值：割线法迭代停止时返回的是alph_0（上一次的迭代值），而非最新计算的alph_1，丢失了最后一次迭代的精度提升。
3. 初始点选择逻辑适配性差：对于c接近1的场景（如0.99），根的取值很小，但初始alph_1选择的是较大的角度，导致初始步长过大，迭代收敛速度慢，甚至提前触发停止条件。
4. 停止条件单一：仅依赖步长变化判断收敛，未结合函数值的误差（|Phi(alpha)-c|），可能在函数值未达到要求精度时就停止迭代。

修复步骤与代码修改 #

1. 替换浮点类型为double #

将所有float替换为double，提升计算精度，确保累积误差可控。

2. 修正迭代返回值 #

停止迭代后返回最新的alph_1，或选择函数值更接近c的点。

3. 优化初始点选择 #

针对c接近1的场景，计算更接近根的初始点，加速收敛。

4. 改进停止条件 #

结合函数值误差和步长变化，确保收敛到足够精度。

修改后的完整代码 #

#include <iostream>
#include <conio.h>
using namespace std;

const double Pi = 3.14159265358979323846;
double abs(double);
double fact(int);
double pot_sgn(int);
double potent(double, int);
double rad_to_deg(double);

double sqrt(double);
double arcsin(double, int iter = 10);
double sin(double, int iter = 10);

double Phi(double, double);
double Solve_Sinc(double, double, double c = 0.7);

int main(void){
        
    double a_size = 0.1;
    
    cout << "__________" << endl;
    cout << Solve_Sinc(a_size, 0.001, 0.35) << endl;  //正确值应为0.00713124
    cout << Solve_Sinc(a_size, 0.01, 0.35) << endl;  //正确值应为0.0713724
    cout << Solve_Sinc(a_size, 0.015, 0.35) << endl;  //正确值应为0.107173
    cout << Solve_Sinc(a_size, 0.04, 0.35) << endl;  //正确值应为0.289264
    cout << Solve_Sinc(a_size, 0.05, 0.35) << endl;  //正确值应为0.364582
    cout << Solve_Sinc(a_size, 0.06, 0.35) << endl;
    cout << Solve_Sinc(a_size, 0.08, 0.35) << endl;
    cout << Solve_Sinc(a_size, 0.1, 0.35) << endl;
    cout << Solve_Sinc(a_size, 0.11, 0.35) << endl;
    cout << Solve_Sinc(a_size, 0.3, 0.35) << endl;
    cout << Solve_Sinc(a_size, 0.5, 0.35) << endl;
    cout << "__________" << endl;
    cout << Solve_Sinc(a_size, 0.001) << endl;
    cout << Solve_Sinc(a_size, 0.01) << endl;
    cout << Solve_Sinc(a_size, 0.015) << endl;
    cout << Solve_Sinc(a_size, 0.04) << endl;  //正确值应为0.180529
    cout << Solve_Sinc(a_size, 0.05) << endl;  //正确值应为0.226366
    cout << Solve_Sinc(a_size, 0.06) << endl;
    cout << Solve_Sinc(a_size, 0.08) << endl;
    cout << Solve_Sinc(a_size, 0.1) << endl;
    cout << Solve_Sinc(a_size, 0.11) << endl;
    cout << Solve_Sinc(a_size, 0.3) << endl;
    cout << Solve_Sinc(a_size, 0.5) << endl;
    cout << "__________" << endl;
    cout << Solve_Sinc(a_size, 0.001, 0.99) << endl;  //正确值应为0.000780871
    cout << Solve_Sinc(a_size, 0.01, 0.99) << endl;  //正确值应为0.00780879
    cout << Solve_Sinc(a_size, 0.015, 0.99) << endl;  //正确值应为0.0117133
    cout << Solve_Sinc(a_size, 0.04, 0.99) << endl;  //正确值应为0.0312399
    cout << Solve_Sinc(a_size, 0.05, 0.99) << endl;  //正确值应为0.0390535
    cout << Solve_Sinc(a_size, 0.06, 0.99) << endl;  //正确值应为0.0468694
    cout << Solve_Sinc(a_size, 0.08, 0.99) << endl;  //正确值应为0.0625104
    cout << Solve_Sinc(a_size, 0.1, 0.99) << endl;  //正确值应为0.0781667
    cout << Solve_Sinc(a_size, 0.11, 0.99) << endl;  //正确值应为0.0860018
    cout << Solve_Sinc(a_size, 0.3, 0.99) << endl;  //正确值应为0.236459
    cout << Solve_Sinc(a_size, 0.5, 0.99) << endl;
    cout << "__________" << endl;
        
    getch();
    return 0;
}

double Solve_Sinc(double a, double l, double c){
    if(a <= 0.0 || l <= 0.0 || c < 0.0){
        cout << "活塞尺寸与波长必须为正数；判定值必须非负！" << endl;
        return 0.0;
    }
    
    double alph_old, alph_0, alph_1, k;
    k = a/l*Pi;
    alph_old = 0.0;
    alph_0 = 0.0;

    // 优化初始点选择：针对c接近1的场景，计算更接近根的初始值
    if(c > 0.9){
        // 利用sinc(x) ≈ 1 - x²/6，求解x ≈ sqrt(6*(1-c))
        double x = sqrt(6*(1 - c));
        alph_1 = arcsin(x * l / (a * Pi));
        // 确保初始点在定义域内
        if(alph_1 > Pi/2) alph_1 = Pi/2;
    } else if(l/a > 1.0){
        if(Phi(k, Pi/2.0) - c > 0.0){return Pi/2.0;}
        alph_1 = Pi/2.0 ;
    } else if(l/a > 0.5){
        alph_1 = arcsin(l/a);
    } else {
        alph_1 = arcsin(2.0*l/a);
    }
    
    bool rep = true;
    int eq = 0;
    const double func_tol = 1e-10; // 函数值误差阈值
    const double step_tol = 1e-12; // 步长阈值
    while(rep){
        alph_old = alph_0;
        alph_0 = alph_1;
        double f0 = Phi(k, alph_0);
        double f_old = Phi(k, alph_old);
        double denom = f0 - f_old;
        
        // 避免分母为零，若接近零则改用小步长迭代
        if(abs(denom) < 1e-15){
            alph_1 = alph_0 + (c > f0 ? 1e-8 : -1e-8);
        } else {
            alph_1 = alph_0 - (f0 - c) * (alph_0 - alph_old) / denom;
        }
        
        // 检查函数值是否已满足精度要求
        if(abs(Phi(k, alph_1) - c) < func_tol){
            rep = false;
            break;
        }
        
        // Garwick步长判断
        double step_curr = abs(alph_1 - alph_0);
        double step_prev = abs(alph_0 - alph_old);
        if(step_curr > step_prev){
            rep = false;
        } else if(step_prev - step_curr <= step_tol){
            eq += 1;
            if(eq == 5){ // 减少连续相等次数，避免不必要迭代
                rep = false;
            }
        } else {
            eq = 0;
        }
    }
    
    // 返回函数值更接近c的点
    double f0 = abs(Phi(k, alph_0) - c);
    double f1 = abs(Phi(k, alph_1) - c);
    double result = (f1 < f0) ? alph_1 : alph_0;
    
    cout << Phi(k, result) - c << " : ";
    return result;
}

double Phi(double k, double x){
    if(abs(x) <= 1e-12){return 1.0;}
    double sin_x = sin(x);
    double arg = k * sin_x;
    if(abs(arg) <= 1e-12){return 1.0;}
    return sin(arg)/arg;
}

double sin(double x, int iter){
    if(x < 0.0){return -sin(-x);}
    
    while (x >= 2.0*Pi){
        x -= 2.0*Pi;
    }
    
    if(x > Pi){return -sin(2.0*Pi - x);}
    
    if(x > Pi/2.0){return sin(Pi-x);}
    
    double z = 0.0;
    for(int i=iter;i>=0;i--){
        z += pot_sgn(i)*(potent(x,2*i+1)/fact(2*i+1));
    }
    return z;
}

double arcsin(double x, int iter){
    if(x < 0.0){return -arcsin(-x);}
    
    if(x - 1.0 >= 1e-12){
        cout << x << ":" << endl;
        cout << "arcsin()的参数绝对值不能大于1！" << endl;
        return 0.0;
    }
    
    if(x >= 0.7071067811865476){return Pi/2.0 - arcsin(sqrt(1.0 - x*x));}
    
    if(x > 0.5 && x < 0.7071067811865476){iter *= 2;}
    double z = 0.0;
    for(int i=iter;i>=0;i--){
        z += (potent(x,2*i+1)*fact(2*i))/(fact(i)*fact(i)*potent(4.0,i)*(2*i+1));
    }
    return z;
}

double sqrt(double x){
    double prec = 1e-12;
    if(abs(x) < prec*10.0){return 0.0;}
    if(x < 0.0){
        cout << "sqrt()的参数必须非负！" << endl;
        return 0.0;
    }
    double y = x;
    double delta;
    double cnt = x;
    if(x >= 1.0){
        do{
            cnt *= 0.1;
            prec *= 10.0;
        }
        while(cnt >= 1.0);
    }
    
    do{
        y = (y + x/y) / 2.0;
        delta = x - y*y;
    } while(delta < -prec || delta > prec);
    
    return y;
}

double rad_to_deg(double x){
    return (x/Pi)*180.0;
}

double potent(double x,int a){
    double y = 1.0;
    for(int i=1;i<=a;i++){
        y *= x;
    }
    return y;
}

double pot_sgn(int n){
    return (n%2 == 0) ? 1.0 : -1.0;
}

double fact(int n){
    double f = 1.0;
    while(n>1){
        f *= n;
        n -= 1;
    }
    
    return f;
}

double abs(double x){
    return (x < 0.0) ? -x : x;
}

修复说明 #

1. 精度提升：所有浮点类型改为double，Pi的取值更精确，辅助函数的迭代次数增加（如sin的iter默认从4改为10），确保三角函数计算精度足够。
2. 初始点优化：针对c>0.9的场景，利用sinc函数的近似公式计算更接近根的初始点，避免大步长迭代导致的收敛问题。
3. 迭代逻辑修正：添加分母为零的处理，避免计算错误；停止条件结合函数值误差和步长变化，确保结果满足精度要求。
4. 返回值优化：比较最后两次迭代的函数值，选择更接近目标c的点作为结果，进一步降低误差。

内容的提问来源于stack exchange，提问作者Cpp Nosavvier