C中arcsin的逼近
我有一个程序可以根据Taylor的序列计算arcsin值的近似值。
我的朋友和我已经提出了一个能够返回几乎“正确”值的算法,但我认为我们并没有非常清晰地完成它。 看一看:
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
非常感谢你们。
您的代码和问题的问题包括:
- 显示arcsin的泰勒级数的图像文件有两个错误: x 5项上有一个减号而不是加号, x的幂显示为x n,但应为x 2 n +1 。
- arcsin的泰勒级数中的x因子在每个项中增加 x 2 ,但是你的公式
a*((my_pow(2*i-1, 2)) / ((2*i)*(2*i+1)*my_pow(x, 2)))
在每个项中除以 x 2 。 这与您询问的特定值-1无关,但除了1之外,它将为其他值产生错误的结果。 - 一旦术语的差异“大于或等于”你的epsilon,你就会问如何结束循环,但是,对于x的大多数值,你实际上想要小于(或者相反,你想要继续 ,而不是结束 ,而差异大于或等于,如您在代码中所示)。
- 泰勒级数是评估函数的一种不好的方法,因为当距离系列居中的位置越远时,它的误差就越大。 像这样的函数的大多数数学库实现使用minimax系列或与之相关的东西。
- 从低阶项到高阶项评估系列会导致您先添加更大的值,然后再添加更小的值。 由于浮点运算的本质,这意味着较小项的精度会丢失,因为它会被较大的值“推出”浮点格式的宽度。 此效果将限制任何结果的准确程度。
- 最后,为了直接解决您的问题,您构建代码的方式,您直接更新
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(ba) > threshold);
所以我想我也必须在while语句中增加它
是的,这可能是一种方式。 是什么阻止了你?
int i=0; while(condition){ //do something i++; }
另一种方法是使用for条件:
for(i = 1; i < 23500 && my_abs(prev_term - next_term) >= EPSILON; i++)
你的公式错了。 这是正确的公式: http : //scipp.ucsc.edu/~haber/ph116A/taylor11.pdf 。
PS还注意到您的公式和系列不相互对应。
您可以像这样使用:
while( std::abs(sum_prev - sum) < 1e-15 ) { sum_prev = sum; sum += a; a = next(a, x, i); }
使用泰勒序列进行arcsin
是非常不精确的,因为这些东西收敛得非常严重,并且对于有限数量的热量,真实的东西会有相对较大的差异。 使用具有整数指数的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)=+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
常数是目标精度。
所以只需用float/double
替换T
并忽略::
…