插值计算椭圆周长并与计算误差

题目

$$I(b)={\int }_0^{\pi /2}\sqrt{1+(b^2-1){\cos }^{2}t}dt$$
上述函数的导数为
$$I^{\prime }(b)={\int }_0^{\pi /2}\frac{b{\cos }^{2}t}{\sqrt{1+(b^2-1){\cos }^{2}t}}dt$$
采用三种方法中最好的方法计算这一积分

(1)利用数值积分的方法给出 $I(b)$ 在 $b=0,0.1,0.2,\dots ,1$ (可以直接计算精确值的，用精确值),用Newton插值方法得到5个椭圆的周长（五个椭圆具体见上一篇图片文件的初步处理与椭圆的最小二乘拟合）。
(2)利用数值积分的方法给出$I^{\prime }(b)$, 在$b=0,0.1,0.2,\dots ,1$ (可以直接计算精确值的，用精确值)，用Hermite插值方法得到5个椭圆的周长。

算法介绍

Netown插值

根据均差的定义，把$x$看成$[a,b]$上的一点，可以得到
$$f(x)=f(x_0)+f\left[x,x_0\right](x-x_0)$$

$$f(x,x_0)=f(x_0,x_1)+f[x,x_0,x_1](x-x_1)$$

将后一个式子带入前一个式子，得到
$$f(x)=f(x_0)+f(x_0)+f(x-x_0)+f[x_0,x_0,x_1](x-x_0)(x-x_0)$$

再由$$f\left[x_1,x_0,x_1\right]=f\left[x_0,x_1,x_2\right]+f\left[x,x_0,x_1,x_2\right](x-x_2)$$

得到
$$f(x)=f(x_0)+f(x_0)+f(x_0-x_1)+f[x_0-x_1)(x-x_0)(x-x_0)(x-x_1)$$

$$+f[x,x_0,x_1,x_2](x-x_0)(x-x_0)(x-x_1)$$

Hermite插值

设在结点 $a\le x_0<x_1<\cdots <x_n\le b$上$f(x_j)=y_j,m_j=f^{\prime }(x_j)$,$(j=0,1,\cdots ,n)$,求多项式$H(X)$,满足条件
$$H(x_j)=y_j,H^{\prime }(x_j)=m_j(j=0,1,\cdots ,n)$$

一共给出了$2n+2$个已知条件，可以唯一确定一个次数不超过$2n+1$的多项式
$$H_{2n+1}(x)=a_0+a_{1}x+\cdots +a_{2n+1}{x}^{2n+1}$$

低阶分段插值

为了避免插值高阶冗余现象导致逼近效果差的现象，我们可以采用分段低次插值得到不同区间内的分段函数。

代码模块

环境：python=3.9.0 sympy=1.12 numpy=1.23.5

导入模块

import sympy as sp
import numpy as np
import time

Newton插值

def Newton_Interger(data_x,formula,t):

    start=time.time()
    # 计算对应点的函数值
    n=len(data_x)
    y=[0]*n
    for i in range(n):
        y[i]=formula.evalf(subs={t:data_x[i]})
    # 创建差商表
    Z=[[0] *n for _ in range(n)]

    for i in range(n):
        Z[i][0]=y[i]
    # 计算差商值
    for i in range(n):
        for j in range(1,i+1):
            Z[i][j]=(Z[i][j-1]-Z[i-1][j-1])/(data_x[i]-data_x[i-j])
    # 计算插值多项式
    res=0 
    for  i in range(n):
        C=Z[i][i]
        for j in range(i):
            C*=(t-data_x[j])
        res+=C
    end=time.time()
    print("time:",end-start)
    
    return res.expand()

Netown法计算椭圆周长

def Compute_N(a,b):

    # 构建函数
    T=sp.Symbol('t')
    formula=(a**2*sp.cos(T)**2+b**2*sp.sin(T)**2)**(1/2)
    data_t=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,\
    1.2,1.3,1.4,1.5]
    formula0=Newton_Interger(formula=formula,data_x=\
        [data_t[0],data_t[-1]],t=T)
    R0=4*Romberg(formula0,T,0,sp.pi/2,1e-6)
    print(R0.evalf())
    print(formula0)
    formula1=Newton_Interger(formula=formula,data_x=\
        [data_t[0],data_t[7],data_t[-1]],t=T)
    R1=4*Romberg(formula1,T,0,sp.pi/2,1e-6)
    print(R1.evalf())
    print(formula1)
    formula2=Newton_Interger(formula=formula,data_x=\
        [data_t[0],data_t[5],data_t[10],data_t[-1]],t=T)
    R2=4*Romberg(formula2,T,0,sp.pi/2,1e-6)
    print(R2.evalf())
    print(formula2)
    formula3=Newton_Interger(formula=formula,data_x=\
        [data_t[0],data_t[4],data_t[8],data_t[12],data_t[-1]],t=T)
    R3=4*Romberg(formula3,T,0,sp.pi/2,1e-6)
    print(R3.evalf())
    print(formula3)
    formula4=Newton_Interger(formula=formula,data_x=\
        [data_t[0],data_t[3],data_t[6],data_t[9],data_t[12],\
            data_t[-1]],t=T)
    R4=4*Romberg(formula4,T,0,sp.pi/2,1e-6)
    print(R4.evalf())
    print(formula4)

计算结果

def Compute(a,b):

    print("Compute_N:")
    Compute_N(a,b)
    #print("Compute_H:")
    #Compute_H(a,b)
    #print("Separte:")
    #Separte(a,b)  

Compute(a=42.8755430086675,b=40.8633914591405)

Compute(a=41.8822717617791,b=39.0938947820857)

Compute(a=40.8633914591405,b=36.5409190520717)

Compute(a=40.3106783148759,b=34.9340215298776)

Compute(a=37.2328566002455,b=30.8883836902119)

误差衡量

对于下述的椭圆形式：
$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1(a>0,b>0).$$

可以通过拉马努金椭圆计算公式快速计算出结果
$$C\approx \pi (a+b)\left(1+\frac{3{\lambda }^2}{10+\sqrt{4-3{\lambda }^2}}\right),\lambda =\frac{a-b}{a+b}.$$

详细代码

这里只是以Newton为例，详细代码见github代码

参考资料

[1].图片文件的初步处理与椭圆的最小二乘拟合

[2].拉马努金椭圆计算公式