我正在尝试评估指数积分$Ei(x)=-\int^{\infty}_{-x}\frac{e^{-t}}{t}dt$为了$x>0$（解释为柯西主值）通过使用有理切比雪夫近似，可以在 Cody & Thacher 的论文“Chebyshev Approximations for the Exponential Integral$Ei(x)$”。该论文可在ams.org在线获得。我正在尝试使用第 292 页上的等式作为间隔$0<x\le6$：

在哪里$x_0\approx{}0.37$是零$Ei(x)$,$p_j$和$q_j$是论文表 II 中的系数（我正在使用$n=9$） 和$T_j^{*}(x)=T_j(2x-1)$是移位的切比雪夫多项式。所以我在一个简单的 C 程序中实现了这个方程

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

#define N_1 10

double chebyshev(int n, double x){
    switch(n){
        case 0:
            return 1.0;
            break;
        case 1:
            return x;
            break;
        default:
            return 2.0*x*chebyshev(n-1, x)-chebyshev(n-2, x);
            break;
    }
}

double shifted_chebyshev(int n, double x){
    return chebyshev(n, 2.0*x-1.0);
}

double interval1(double x){
    static const double pj_1[N_1] = {
                            -4.1658081333604994241879E11,
                            +1.2177698136199594677580E10,
                            -2.5301823984599019348858E10,
                            +3.1984354235237738511048E8,
                            -3.5377809694431133484800E8,
                            -3.1398660864247265862050E5,
                            -1.4299841572091610380064E6,
                            -1.4287072500197005777376E4,
                            -1.2831220659262000678155E3,
                            -1.2963702602474830028590E1
    };
    static const double qj_1[N_1] = {
                            -1.7934749837151009723371E11,
                            +9.8900934262481749439886E10,
                            -2.8986272696554495342658E10,
                            +5.4229617984472955011862E9,
                            -7.0108568774215954065376E8,
                            +6.4698830956576428587653E7,
                            -4.2648434812177161405483E6,
                            +1.9418469440759880361415E5,
                            -5.5648470543369082846819E3,
                            +7.6886718750000000000000E1
    };
    static const double x0 = 0.372507410781366634461991866580;

    int j;
    double result, numerator, denominator, sum;

    result = log(x/x0);
    numerator = 0.0;
    denominator = 0.0;
    sum = 0.0;

    for(j=0; j<N_1; j++){
        numerator += pj_1[j] * shifted_chebyshev(j, x/6.0);
        denominator += qj_1[j] * shifted_chebyshev(j, x/6.0);
    }
    sum = numerator / denominator;
    sum *= (x-x0);
    result += sum;

    return result;
}

double ei(double x){
    if( x < 0.0 ){
        printf("Argument of ei(x) must be positive\n");
        exit(1);
    }

    if( x <= 6.0 ){
        return interval1(x);
    }

    printf("Argument range not yet implemented for ei(x).\n");
    exit(1);
}

int main(void){
    double x = 5.0;
    printf("ei(%g)=%g\n", x, ei(x));

    return 0;
}

但是对于有问题的间隔，结果完全不正确。例如我想评估$Ei(5)\approx{}40.2$，但我得到$\approx{}19.0654$. 该论文声称最大相对误差为$8\cdot10^{-19}$，所以我想我一定做错了什么。我仔细检查了系数和方程式的实现，但我没有看到任何错误。

当我绘制 C 计算值与实际/预期值（由 Mathematica 的 ExpIntegralEi 函数返回）相比的相对误差时，我可以看到只有第一个间隔（$0<x\le6$) 错误高得不合理：

绘制的是$y(x)=\left|\frac{C\left(x\right)-M\left(x\right)}{M\left(x\right)}\right|$， 在哪里$C(x)$是我的 C 程序返回的值，而 M(x) 与 Mathematica 返回的值相同。

任何帮助将非常感激。

参考

1. Cody、WJ 和 Henry C. Thacher。“指数积分 𝐸𝑖 (𝑥) 的切比雪夫近似。” 计算数学 23.106（1969）：289-303。