C中arcsin的近似值

Question

我有一个程序可以根据泰勒级数计算反正弦值的近似值。

我和我的朋友想出了一个算法，它能够返回几乎“正确”的值，但我认为我们做得不是很清楚。看看：

double my_asin(double x)
{
    double a = 0;
    int i = 0;
    double sum = 0;
    a = x;
    for(i = 1; i < 23500; i++)
    {
        sum += a;
        a = next(a, x, i);
    }
}

double next(double a, double x, int i)
{
    return a*((my_pow(2*i-1, 2)) / ((2*i)*(2*i+1)*my_pow(x, 2)));
}

我检查了 my_pow 是否正常工作，所以我也不需要在这里发布它。基本上，一旦当前项和下一项之间的差异大于或等于我的 EPSILON (0.00001)，我希望循环结束，这是我在计算平方根时使用的精度。

---

这就是我希望它的工作方式：

while(my_abs(prev_term - next_term) >= EPSILON)

但是函数 double next 依赖于 i，所以我想我也必须在 while 语句中增加它。有什么想法我应该怎么做？

---

-1 的示例输出：

$ -1.5675516116e+00

代替：

$ -1.5707963268e+00

---

非常感谢大家。

Answer 1

您的代码和问题的问题包括：

1. 显示 arcsin 泰勒级数的图像文件有两个错误：x⁵ 项上有一个减号，而不是一个加号，以及x 显示为 xⁿ 但应为 x²ⁿ⁺¹.
2. 泰勒级数中的 x 因子在每一项中增加 x² ，但您的公式 a*((my_pow(2*i-1, 2)) / ((2*i)*(2*i+1)*my_pow(x, 2))) 在每个术语中 除 除以 x²。这对于您询问的特定值 -1 无关紧要，但对于除 1 之外的其他值，它会产生错误的结果。
3. 您询问如何在术语差异“大于或等于”您的 epsilon 时结束循环，但是对于 x 的大多数值，您实际上希望小于（或者相反，您希望继续，而不是结束，而差异大于或等于，如您在代码中所示）。
4. Taylor 级数是评估函数的一种糟糕方法，因为它的误差会随着您远离级数的中心点而增加。大多数此类函数的数学库实现都使用极小极大数列或与之相关的东西。
5. 从低阶项到高阶项评估序列会导致您先添加较大的值，然后再添加较小的值。由于浮点运算的性质，这意味着较小项的准确性会丢失，因为它被较大的值“推出”了浮点格式的宽度。这种影响会限制任何结果的准确性。
6. 最后，为了直接解决您的问题，您构建代码的方式，您直接更新a，因此您永远不会同时拥有上一期和下一期。相反，请创建另一个 double b，以便为上一个术语创建一个对象 b，为当前术语创建一个对象 a，如下所示。

例子：

double a = x, b, sum = a;
int i = 0;
do
{
    b = a;
    a = next(a, x, ++i);
    sum += a;
} while (abs(b-a) > threshold);

Answer 2

对arcsin 使用 Taylor 级数 是极其不精确的，因为这些东西的收敛性非常差，并且对于有限数量的 therms，与真实的东西会有相对较大的差异。同样使用带有整数指数的pow 也不是很精确和有效。

不过使用arctan 就可以了

arcsin(x) = arctan(x/sqrt(1-(x*x)));

当它的泰勒级数在<0.0,0.8> 范围内收敛时，范围的所有其他部分都可以通过它计算（使用三角恒等式）。所以这里是我的 C++ 实现（来自我的算术模板）：

T atan    (const T &x)                                              // = atan(x)
    {
    bool _shift=false;
    bool _invert=false;
    bool _negative=false;
    T z,dz,x1,x2,a,b; int i;
    x1=x; if (x1<0.0) { _negative=true; x1=-x1; }
    if (x1>1.0) { _invert=true; x1=1.0/x1; }
    if (x1>0.7) { _shift=true; b=::sqrt(3.0)/3.0; x1=(x1-b)/(1.0+(x1*b)); }
    x2=x1*x1;
    for (z=x1,a=x1,b=1,i=1;i<1000;i++)  // if x1>0.8 convergence is slow
        {
        a*=x2; b+=2; dz=a/b; z-=dz;
        a*=x2; b+=2; dz=a/b; z+=dz;
        if (::abs(dz)<zero) break;
        }
    if (_shift) z+=pi/6.0;
    if (_invert) z=0.5*pi-z;
    if (_negative) z=-z;
    return z;
    }
T asin    (const T &x)                                              // = asin(x)
    {
    if (x<=-1.0) return -0.5*pi;
    if (x>=+1.0) return +0.5*pi;
    return ::atan(x/::sqrt(1.0-(x*x)));
    }

其中T 是任何浮点类型（float,double，...）。如您所见，您需要实现sqrt(x)、pi=3.141592653589793238462643383279502884197169399375105、zero=1e-20 和+,-,*,/ 操作。 zero 常量是目标精度。

所以只需将T 替换为float/double 并忽略:: ...

Answer 3

所以我想我也必须在 while 语句中增加它

是的，这可能是一种方法。是什么阻止了你？

int i=0;
while(condition){
   //do something
   i++;
}

另一种方法是使用 for 条件：

for(i = 1; i < 23500 && my_abs(prev_term - next_term) >= EPSILON; i++)

Answer 4

你的公式是错误的。这是正确的公式：http://scipp.ucsc.edu/~haber/ph116A/taylor11.pdf。

附：另请注意，您的公式和您的系列不对应。

---

你可以像这样使用while：

while( std::abs(sum_prev - sum) < 1e-15 )
    {
        sum_prev = sum;
        sum += a;
        a = next(a, x, i);
    }