0

我需要使用泰勒级数在没有 math.h 库的情况下编写自己的 asin() 函数。它适用于 <-0.98;0.98> 之间的数字，但是当我接近极限时，它会停止 1604 次迭代，因此不准确。

我不知道如何使它更准确。任何建议都非常感谢！

代码如下：

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

#define EPS 0.000000000001

double my_arcsin(double x)
{
    long double a, an, b, bn;
    a = an = 1.0;
    b = bn = 2.0;
    long double n = 3.0;
    double xn;
    double xs = x;
    double xp = x;

    int iterace = 0;

    xn = xs + (a/b) * (my_pow(xp,n) / n);

    while (my_abs(xn - xs) >= EPS)
    {
        n += 2.0;
        an += 2.0;
        bn += 2.0;
        a = a * an;
        b = b * bn;

        xs = xn;
        xn = xs + (a/b) * (my_pow(xp,n) / n);
        iterace++;
    }

    //printf("%d\n", iterace);

    return xn;
}

int main(int argc, char* argv[])
{

    double x = 0.0;

    if (argc > 2)
        x = strtod(argv[2], NULL);
    if (strcmp(argv[1], "--asin") == 0)
    {
           if (x < -1 || x > 1)
               printf("nan\n");
           else
           {
               printf("%.10e\n", my_arcsin(x));
               //printf("%.10e\n", asin(x));
           }

        return 0;
    }
}

还有我的价值观和预期价值观的简短列表：

My values              Expected values        my_asin(x)
5.2359877560e-01       5.2359877560e-01       0.5
1.5567132089e+00       1.5707963268e+00       1      //problem
1.4292568534e+00       1.4292568535e+00       0.99   //problem
1.1197695150e+00       1.1197695150e+00       0.9
1.2532358975e+00       1.2532358975e+00       0.95
4

2 回答 2

5

即使您使用的级数展开式的收敛半径是 1，因此级数最终会在 -1 < x < 1 时收敛，但在接近此区间的极限时收敛确实非常缓慢。解决方案是以某种方式避免间隔的这些部分。

我建议你

  • 对 |x| 使用您的原始算法 <= 1/sqrt(2),
  • 对 1/sqrt(2) < x <= 1.0 使用恒等式 arcsin(x) = pi/2 - arcsin(sqrt(1-x^2))，
  • 对 -1.0 <= x < -1/sqrt(2) 使用恒等式 arcsin(x) = -pi/2 + arcsin(sqrt(1-x^2))。

这样，您可以将输入 x 转换为 [-1/sqrt(2),1/sqrt(2)]，其中收敛速度相对较快。

于 2013-11-25T15:26:00.167 回答
0

请注意：在这种情况下，我强烈推荐@Bence 的方法，因为您不能指望具有低数据精度的缓慢收敛方法获得任意精度。

但是，我愿意向您展示如何使用您当前的算法改进结果。

主要问题是它增长太快并且很快就变成了a（仅在大约 150 次迭代之后）。另一个类似的问题是在增长时快速增长，但是在这种情况下这并不重要，因为我们可以假设输入数据在.binfmy_pow(xp,n)n[-1, 1]

所以我刚刚a/b通过引入更改了您处理的方法ab_ratio，请参阅我编辑的代码：

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

#define EPS 0.000000000001

#include <math.h>
#define my_pow powl
#define my_abs fabsl

double my_arcsin(double x)
{
    #if 0
    long double a, an, b, bn;
    a = an = 1.0;
    b = bn = 2.0;
    #endif
    unsigned long _n = 0;
    long double ab_ratio = 0.5;
    long double n = 3.0;
    long double xn;
    long double xs = x;
    long double xp = x;

    int iterace = 0;

    xn = xs + ab_ratio * (my_pow(xp,n) / n);
    long double step = EPS;

    #if 0
    while (my_abs(step) >= EPS)
    #else
    while (1) /* manually stop it */
    #endif
    {
        n += 2.0;
        #if 0
        an += 2.0;
        bn += 2.0;
        a = a * an;
        b = b * bn;
        #endif
        _n += 1;
        ab_ratio *= (1.0 + 2.0 * _n) / (2.0 + 2.0 * _n);

        xs = xn;
        step = ab_ratio * (my_pow(xp,n) / n);
        xn = xs + step;
        iterace++;
        if (_n % 10000000 == 0)
            printf("%lu %.10g %g %g %g %g\n", _n, (double)xn, (double)ab_ratio, (double)step, (double)xn, (double)my_pow(xp, n));
    }

    //printf("%d\n", iterace);

    return xn;
}

int main(int argc, char* argv[])
{

    double x = 0.0;

    if (argc > 2)
        x = strtod(argv[2], NULL);
    if (strcmp(argv[1], "--asin") == 0)
    {
           if (x < -1 || x > 1)
               printf("nan\n");
           else
           {
               printf("%.10e\n", my_arcsin(x));
               //printf("%.10e\n", asin(x));
           }

        return 0;
    }
}

对于0.99（甚至0.9999999），它很快就会给出超过 10 个有效数字的正确结果。但是接近时会变慢1。
实际上这个过程在我的笔记本电脑上计算已经运行了将近 12 分钟--asin 1，当前的结果是1.570786871经过3560000000迭代的。

更新：现在已经 1 小时 51 分钟，结果1.570792915和迭代次数为27340000000.

于 2013-11-25T16:27:39.153 回答