数值积分求解卫星轨道长(Python)

目录

1、概述

2、代码方法一

3、代码方法二

 4、运行结果

                                    

---

1、概述

        微积分在科学和工程中有着广泛的应用，例如，人造地球卫星的轨道可视为平面的椭圆，我国第一颗人造地球卫星近地球距离表面439km，远地点距离地球表面2384km，地球半径为6371km，求该卫星的轨道长度。

       本文题可用椭圆参数方程来描述人造地球卫星轨道，取a=8755km,b=6810km分别为椭圆的长短半轴，该轨道的长度就是参数方程的弧长积分。

2、代码方法一

import numpy as np
def fun(x):
    a,b=8755,6810
    return (a**2*np.sin(x)**2+b**2*np.cos(x)**2)**(1/2)
    return x
#Compound trapezoid
def tx(a,b,n):
    h=(b-a)/n
    x=a
    s=fun(x)+fun(b)
    for i in range(1,n):
        x=x+h
        s=s+2*fun(x)
    result=(h/2)*s
    return result
#Composite Simpson
def xps(a,b,n):
    h=(b-a)/n
    x=a
    s=fun(x)+fun(b)
    for i in range(1,n):
        x=x+h/2
        s=s+4*fun(x)
        x=x+h/2
        s=s+2*fun(x)
    result=(h/6)*s
    return result
a,b=0,np.pi/2
n=20000
t=tx(a,b,n)
p=xps(a,b,n)
print(t,p)

   3、代码方法二

import numpy as np
#def funcval(x):
  #  return np.exp**(1/X)
def funcvals(x1):
    a,b=8755,6810
    return (a**2*np.sin(x)**2+b**2*np.cos(x)**2)**(1/2)
#梯形公式
def trapezoid(a,b,n,f):
    sum=0
    sum=(f(a)+f(b))*(b-a)/2
    return sum
#中矩形公式
def retangle(a,b,n,f):
    sum=0
    sum=f((a+b)/2)*(b-a)
    return sum
#simpson
def simpson(a,b,n,f):
    sum=0
    sum=4*f((a+b)/2)+f(a)+f(b)
    return sum*(b-a)/6
#ftrapezoid
def trapezoid(a,b,n,f):
    h=(b-a)/(n-1)
    x=np.zeros(n,1)
    for i in range(n):
        x[i]=a+i*h
        sum=0
        sum=f(x[0])+f(x[-1])
        sum=sum+2*np.sum(f(x[1:-1]))
        sum=sum*h/2
        return sum
#fsimpson
def fsimpson(a,b,n,funcval1):
    h=(b-a)/(n-1)
    x=np.zeros(n,1)
    for i in range(n):
        x[i]=a+i*h
        h=h*2
        sum=f(x[0])+f(x[-1])
        if x.shape[0]>2:
            #griad
            fg=f(x[2:-1:2])
            fh=f(x[1:-1:2])
            #half grid

        sum=sum+2*np.sum(fg)+4*np.sum(fh)
        sum=sum*h/6
        return sum

def main():
    a,b=1,2
    n=20
    tr=trapezoid(a,b,n,funcval)
    sip=simpson(a,b,n,funcval)
    ft=ftrapezoid(a,b,nfuncval)
    fs=fsimpson(a,b,n,funcval)

    a,b=0,np.pi/2
    fs=fsimpon(a,b,n,funcval1)
    print(fs*4)
if __name__=='__main__':
    main

 4、运行结果

12272.491317172338 12272.032906443934

Process finished with exit code 0