数值方法求π和π的x次方

学习了数值分析这门课，老师留了一个作业是设计数值方法求π的x次方（x范围为1到10）
具体步骤如下（老师要求的，就是想让你多用用几种数值方法，让你多解几步）

1. 求π
2. 利用上一步的π求lnπ
3. 利用上一步的lnπ求π的x次方
4. 误差分析
   不难
   但是误差分析很麻烦，因为要考虑每一步的误差对下一步的影响，还需要考虑舍入误差和方法误差，误差分析在文章的末尾

求π

求解π的方法有很多，现在只讨论用计算机怎么求π。一般有下面这三种求法
1.泰勒展开
2.BBP公式
3.高斯勒让德算法（最好的一种）
这三种方法具体可以参考https://blog.csdn.net/xfxyy_sxfancy/article/details/48378121
这个博文很详细我就不再废话了，下面给出泰勒展开求π的C++代码，展开40项即可得到很精确的结果

    int h1,i,t;
    double dsqrt3=2*sqrt(3);//dsqrt为2倍的根号3，是一个常量
    pi=dsqrt3;
    h1=1; //h1为正1或负1
    for(i=1;i<40;i++)//采用泰勒展开方法求π
    {
        t=2*i+1;
        h1=-h1;
        pi=pi+h1*dsqrt3/t/pow(3,i);
    }
    cout << "π的值：";
    cout << fixed << setprecision(15)<< pi << endl;

求lnπ

可以用数值积分的方式求

注意：上式中的π用的是第一步求出的π
数值积分有很多方法，比如梯形公式，辛普森公式等等
我采用的是复合辛普森公式（可能不是最好的），原理如下

C++代码如下

    double h=(pi-1)/500;//pi即为上一步求的π
    lnpi=1+1/pi;//设定初值
    for(i=1;i<500;i++)//采用复合辛普森公式求数值积分
    {
        x2=1+i*h;
        x1=1+i*h-h/2;
        lnpi=lnpi+4/x1+2/x2;
    }
    lnpi=lnpi+4/(pi-h/2);
    lnpi=lnpi/6*h;
    cout << "lnπ的值：";
    cout << fixed << setprecision(15)<< lnpi << endl;

求π的x次方

要求用到前两步的结果，我采用求解常微分方程的方法

求解常微分方程也有很多方法，如各种欧拉法，龙格库塔，等等
我采用的是四阶龙格库塔公式，如下

C++代码如下，给了个龙格库塔的函数，该函数输入
x：几次方
n：区间分为几份，份数越多越精确
输出：pix，即结果，表示π的x次方

double RungKutta(double x,double n)//四阶龙格-库塔公式
{
    double pix,K1,K2,K3,K4,h;
    pix = 1.0;
    h=x/n;
    for(int i=0;i<n;i++)
    {
        K1 = derivative(pix);
        K2 = derivative(pix+0.5*h*K1);
        K3 = derivative(pix+0.5*h*K2);
        K4 = derivative(pix+h*K3);
        pix = pix+h*(K1+2*K2+2*K3+K4)/6;
    }
    return pix;
}
double derivative(double y)//y（x）的导数，用于微分方程求解
{
    return (lnpi*y);//lnpi即自己求出的lnπ
}

重中之重：误差分析

这一部分的分析很多，足足有四页，因此在这里就不放了，详情请见评论区我给的链接，点开就是