acosx的泰勒展开、C语言实现以及实际应用

123qwertybobo

acosx的泰勒展开在大学书本基本就是一带而过，这个公式其实既简单又复杂，但用acosx求θ又是一项很重要的事情。为了一个下午的时间，mark一下！

acosx=π/2-(x+x^3/3+……+通项an+……)，an=(n-2)!!/(n-1)!!*(x^n)/n,这里n是奇数，其中n!!=n*(n-2)*(n-4)*…*1，具体举例如5！！=5*3*1,4！！=4*2。

先分析函数，可以拆成两个括号，acosx=(π/2-x)-(0.5*x^3/3+……+通项+……)。这样的话，后一个括号才能用循环去表示。
该循环采用for去思考。通项由两部分构成，就是(n-2)!!/(n-1)!!以及(x^n)/n。n一般就是for里的i=3;i<=XXX;i+=2这个问题不大。
唯一的难点就在于(n-2)!!/(n-1)!!的计算。将其设为k，那么k就是间隔为2的阶乘。那么转为双重循环的问题。
先考虑如果XXX=5的时候，那么就有n=3以及n=5两个式子。而k(n=3)=1/2,k(n=5)=1*3/2/4。可以知道：k(n=7)=k(n=5)*5/6，即k(n=XXX)=k(n=XXX-2)*(XXX-2)/(XXX-1)，那么k*=j/(j+1).
内层循环就是for (j=1;j<=i-2;j+=2)。当i=3只有j=1，k(i=3,j=1)=1/2，没问题。当i=5时，j可以为1或3，那么k（i=5,j=1）=1/2*1/2，k（i=5,j=3）=1/2*1/2*2/4。这就有问题了。此时必须不让j=1，所以必须用if把j<i-2的k都过滤掉。见下：
       
                for (j=1;j<=i-2;j+=2)
                {
                        if (j<i-2)
                continue;  //这步就是遇到所有j<i-2的情况一经判断就返回到for的循环，直至j=i-2的情况。
                else
                        {
                                k *= j/(j+1);
                                sum+=-k*pow(x,i)/i;
                        }       
                }
采用excel的acos函数，acos（-0.999854）=3.124504438，这个是弧度制，单位为rad。让XXX=40000，sum=3.124504说明代码的可行性。有一点千万要注意，i，j这些参数千万不能设置成int，否则一运算就是0！！！

该问题最相关的一个例子就是两个面的夹角，它源于“批量求化学结构中的特定四个原子的二面角”。如果采用awk的话，由于awk内置cos函数，可以将该问题简化。假设一个α在[0,π]之间，cosα=-0.999854，那么α=0,cosα=1。由于cosα在[0,π]单调递减，那么只要cosα>-0.999854,让α一直增加进而逼近实际值。如do {a+=i/10000;} while （cosα>-0.999854）。

附C语言：
#include<stdio.h>
#include<math.h>



int main ()
{
        double i,j;
        double x=-0.999854;
        double k =1;
        double sum = 3.141592653590/2-x;
        for (i=1;i<=40000;i+=2)
        {
                for (j=1;j<=i-2;j+=2)
                {
                        if (j<i-2)
                continue;
                else
                        {
                                k *= j/(j+1);
                                sum+=-k*pow(x,i)/i;
                        }       
                }
        }
        printf("%lf次迭代%lf",i,sum);
}