【计算方法】数值微积分

MorleyOlsen

已于 2023-10-15 17:01:25 修改

阅读量1k

点赞数 2

分类专栏：计算方法文章标签：算法 python 学习方法经验分享

于 2022-10-29 11:22:35 首次发布

本文链接：https://blog.csdn.net/m0_65787507/article/details/127584509

版权

计算方法专栏收录该内容

7 篇文章

订阅专栏

0.WARNINGS

本文章主要适用于计算方法代码的实现参考，由于本人是python究极小白，为了实验课速成了一些内容，因此会包含较多的暴力解法orz，如有代码错误、可优化的地方欢迎各位大佬指出，感激不尽。

此外，本人制作本文的目的，主要是懒的将代码保存在本地中，且博客方便个人的复习。（写这篇的时候已经知道下周要小测了，哇的一声哭出来！）

数值微积分主要包含的内容有：求积公式&代数精度；梯形公式；牛顿-科特斯公式；复化求积公式（Newton、Cotes、Gauss、Romberg、梯形等）；龙贝格公式；高斯型求积公式（切比雪夫、勒让德）；数值微分（向前、向后、中间）。

1.用辛普森公式求解定积分

import math

inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])

def fun(x):
    y=x*x*((math.e)**x)
    return y

def simpson(a,b):
    ret=1.0/6 * (b-a) * (fun(a) + 4*fun(a/2+b/2) + fun(b))
    return ret

print("%.5f" % simpson(a,b))

2. 用科特斯公式计算定积分

import math

inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])

def fun(x):
    #x=0的时候分母没有意义！！！！！！！！！！！
    if x==0:
        y=1
    else:
        y=math.sin(x)/x
    return y

def cotes(a,b):
    temp=float((b-a)/4)
    ret=(b-a) * (fun(a)*(7/90) + fun(a+temp)*(16/45) + fun(a+temp*2)*(2/15) + fun(a+temp*3)*(16/45) +fun(b)*(7/90))
    return ret

result = cotes(a,b)

print("%.5f" % result)

3.比较梯形辛普森柯特斯三种不同的积分公式的求积结果

import math

inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])

def fun(x):
    y=x * math.sin(x) * math.e**x
    return y

def simpson(a,b):
    ret1=1.0/6 * (b-a) * (fun(a) + 4*fun(a/2+b/2) + fun(b))
    return ret1

def cotes(a,b):
    temp=float((b-a)/4)
    ret2=(b-a) * (fun(a)*(7/90) + fun(a+temp)*(16/45) + fun(a+temp*2)*(2/15) + fun(a+temp*3)*(16/45) +fun(b)*(7/90))
    return ret2

def ti(a,b):
    ret3=1.0/2 * (fun(a) + fun(b)) * (b-a)
    return ret3

myti=ti(a,b)
mysimpson=simpson(a,b)
mycotes=cotes(a,b)
print("%.5f %.5f %.5f" % (myti,mysimpson ,mycotes))

4. 用复化梯形公式计算式子的近似值

import math

inp=input().split(' ')
a=float(inp[0])
b=float(inp[1])
n=int(inp[2])

def fun(x):
    y=math.log(x)
    return y

def ti(a,b,n):
    h=(b-a)/n
    T=0
    for i in range(1,n):
        T += fun(a+i*h)
    T = h/2 * (fun(a)+2*T+fun(b))
    return T


ret=ti(a,b,n)
print("%.5f" % ret)

5. 用复化辛卜生公式计算积分的近似值

import math

inp=input().split(' ')
a=float(inp[0])
b=float(inp[1])
n=int(inp[2])

def fun(fai):
    y=math.sqrt(4-math.sin(fai)*math.sin(fai))
    return y

def simpson(a,b,n):
    h=(b-a)/n
    S=0
    for i in range(1,n):
        S += 2*fun(a+i*h)
    j=a+h/2
    while j<b:
        S += 4*fun(j)
        j += h
    S = h/6 * (fun(a)+S+fun(b))
    return S

ret=simpson(a,b,n)
print("%.5f" % ret)

6.变步长梯形积分求解定积分

import math

inp=input().split(' ')
a=float(inp[0])
b=float(inp[1])
e=float(inp[2])

def fun(x):
    if x==0:
        return 1
    else:
        y=math.sin(x)/x
        return y

h=b-a
T1 =h/2*(fun(b)+fun(a))

while 1:
    s=0
    x=a+h/2
    while(x<b):
        s+=fun(x)
        x+=h
    T2=T1/2+h/2*s
    if(math.fabs(T2-T1)<=e):
        break
    else:
        T1=T2
        h/=2

ret=T2
print("%.6f" % ret)

7. 用龙贝格积分法计算积分

import math
import numpy as np


inp=input().split(' ')
a=float(inp[0])
b=float(inp[1])

def fun(x):
    y=4/(1+x**2)
    return y

def romberg(a,b,n):
    cnt=0
    while n>1:
        n/=2
        cnt+=1
    sum=cnt+1 + 3
    array=np.zeros(sum)

    # T1 is array[0]
    h=b-a
    array[0]=h/2 * (fun(a)+fun(b))
    # print(array[0])
    for i in range(1,sum):
        s=0
        x=a+h/2
        while x<b:
            s+=fun(x)
            x+=h
        T2=array[i-1]/2+h/2*s
        h/=2
        array[i]=T2
        # print(array[i])

    # Sn
    count=sum-1
    # print(count)
    while count>0:
        array[count]=array[count]*4/3-array[count-1]*1/3
        # print(array[count])
        count -= 1

    # Cn
    count=sum-1
    while count>1:
        array[count]=array[count]*16/15-array[count-1]*1/15
        # print(array[count])
        count -= 1

    # Rn
    count=sum-1
    while count>2:
        array[count] = array[count] * 64 / 63 - array[count - 1] * 1 / 63
        # print(array[count])
        count -= 1

    # output R4
    print("%.2f" % array[sum-1])
"""
    # from R1 ~ Rn
    for i in range(3,sum):
        print(array[i])
"""

romberg(a,b,4)

8.用n=2高斯-勒让德求积公式计算下列积分

import math

inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])

def fun(x):
    y=math.sqrt(x+1.5)
    return y

def gauss(a,b):
    x0=-math.sqrt(15)/5
    x1=0
    x2=math.sqrt(15)/5
    myx0=(b-a)/2 * x0 + (a+b)/2
    myx1 = (b - a) / 2 * x1 + (a + b) / 2
    myx2 = (b - a) / 2 * x2 + (a + b) / 2
    I=5/9 * fun(myx0)+ 8/9*fun(myx1 )+5/9*fun(myx2)
    I*=(b-a)/2
    return I

ret=gauss(a,b)
print("%.5f" % ret)

9. 比较4种高精度的求积方法

import math
import numpy as np

inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])

def fun(x):
    if x==0:
        return 0
    else:
        return 1/x

def romberg(a,b,n):
    cnt=0
    while n>1:
        n/=2
        cnt+=1
    sum=cnt+1 + 3
    array=np.zeros(sum)

    # T1 is array[0]
    h=b-a
    array[0]=h/2 * (fun(a)+fun(b))
    # print(array[0])
    for i in range(1,sum):
        s=0
        x=a+h/2
        while x<b:
            s+=fun(x)
            x+=h
        T2=array[i-1]/2+h/2*s
        h/=2
        array[i]=T2
        # print(array[i])

    # Sn
    count=sum-1
    # print(count)
    while count>0:
        array[count]=array[count]*4/3-array[count-1]*1/3
        # print(array[count])
        count -= 1

    # Cn
    count=sum-1
    while count>1:
        array[count]=array[count]*16/15-array[count-1]*1/15
        # print(array[count])
        count -= 1

    # Rn
    count=sum-1
    while count>2:
        array[count] = array[count] * 64 / 63 - array[count - 1] * 1 / 63
        # print(array[count])
        count -= 1

    # return R4
    return array[sum-1]

def change(a,b,x):
    x=(b-a)/2 * x + (a+b)/2
    return x

def twolegendre(a,b):
    I=0
    x0=-1/math.sqrt(3)
    x1=-x0
    x0=change(a,b,x0)
    x1=change(a,b,x1)
    I=fun(x0)+fun(x1)
    I*=(b-a)/2
    return I

def threelegendre(a,b):
    I=0
    x0=-math.sqrt(15)/5
    x1=0
    x2=-x0
    x0=change(a,b,x0)
    x1=change(a,b,x1)
    x2 = change(a, b, x2)
    I=5/9*fun(x0)+8/9*fun(x1)+5/9*fun(x2)
    I*=(b-a)/2
    return I

def fivelegendre(a,b):
    I=0
    x0=-0.9061798
    x1=-0.5384693
    x2=0
    x3=-x1
    x4=-x0
    x0 = change(a, b, x0)
    x1 = change(a, b, x1)
    x2 = change(a, b, x2)
    x3 = change(a, b, x3)
    x4 = change(a, b, x4)
    A0=0.2369269
    A1=0.4786287
    A2=0.5688889
    I=A0*(fun(x0)+fun(x4))+A1*(fun(x1)+fun(x3))+A2*fun(x2)
    I *= (b - a) / 2
    return I

def complexgauss(a,b):
    I=0
    h=(b-a)/4
    for i in range(0,4):
        newa=a+i*h
        newb=a+(i+1)*h
        I+=twolegendre(newa,newb)
    return I

ret1=romberg(a,b,4)
ret2=threelegendre(a,b)
ret3=fivelegendre(a,b)
ret4=complexgauss(a,b)
print("%.5f" % ret1)
print("%.5f" % ret2)
print("%.5f" % ret3)
print("%.5f" % ret4)

10.用高斯切比雪夫求积公式求积分

import math
"""
inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])
"""
def fun(x):
    y=math.e**x
    return y

def chebyshev(n):
    A=math.pi/(n+1)
    I=0
    for i in range(0,n+1):
        x=math.cos((2*i+1)/(2*n+2)*math.pi)
        I+=fun(x)
    I*=A
    return I


n=int(input())
# x0,x1,...xn个高斯点
ret=chebyshev(n)
print("%.6f" % ret)

11.用向前差商公式、向后差商公式和中心差商公式计算定在sinx/x点的导数

import math
"""
inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])
"""
def fun(x):
    if x==0:
        return 1
    else:
        y = math.sin(x) / x
        return y

def forward(x):
    h=1
    temp=0
    while 1:
        f=1/h * (fun(x+h)-fun(x))
        if math.fabs(f-temp)<0.001:
            break
        else:
            temp=f
            h/=2
    return temp # 输出前点的导数值

def backward(x):
    h=1
    temp=0
    while 1:
        f=1/h * (fun(x)-fun(x-h))
        if math.fabs(f-temp)<0.001:
            break
        else:
            temp=f
            h/=2
    return temp

def mid(x):
    h = 1
    temp = 0
    while 1:
        f=1/2/h * (fun(x+h)-fun(x-h))
        if math.fabs(f-temp)<0.001:
            break
        else:
            temp=f
            h/=2
    return temp

myx=float(input())
print("%.3f %.3f %.3f" % (forward(myx),backward(myx),mid(myx)))

12.利用数值求导的三点公式计算人口增长率

import math
import numpy as np

year=[float(e) for e in input().split()]

population=[float(t) for t in input().split()]

length=len(year)
# print(length)

h=year[1]-year[0]
# print(h)
rate=np.zeros(length)

A=1/(2*h)
for i in range(0,length):
    if i==0:
        rate[i]=A*(-3*population[0]+4*population[1]-population[2])
    elif i==length-1:
        rate[i]=A*(population[length-3]-4*population[length-2]+3*population[length-1])
    else:
        rate[i]=A*(population[i+1]-population[i-1])
    print("%.3f" % rate[i])