用编程实现计算方法/数值分析的概念知识 （三）

数值积分与数值微分

1. 用复化梯形公式、复化辛普森公式计算积分

   $$ I=\int^1_0 \frac{4}{1+x^2}dx $$

解：

# =============================================================================
# 定义区间长度, 子区间数, 步长
l = 1           # 区间长度
n = 8           # 子区间数
h = 1/8         # 步长


# 定义被积函数
def f(x):
    return 4/(1+x**2)

# =============================================================================
# 构造字典，存放分点数据
table = {}
for i in range(0,9):
    table[i*h] = f(i*(1/8))

# =============================================================================
# 定义复化梯形公式
def T(n, h):
    res = table[0]+table[1]
    for i in range(1, n):
        res += 2*table[i*h]
    res*=(1/2 * h)
    return res


# 定义复化辛普森公式
def S(n, h):
    res = table[0]+table[1]
    for i in range(1,n):
        res += 2*table[i*h]
    for i in range(0,n):
        res += 4*table[i*h+h/2]

    res*=(1/6 * h)
    return res

# ==============================================================================
def main():
    print('数据表格如下')
    for key in table:
       print('f('+str(key)+'):'+str(table[key]))
    print('\n')
    
    
    print('复化梯形公式得:')
    print(T(8,1/8))
    print('\n')


    print('复化辛普森公式得:')
    print(S(4,1/4))

if __name__ == '__main__':
    main()

但是这个实现方法并不精美，数值积分中每项的系数是变化的，显然不能写死在代码里，所以对自己之前写的代码进行了一段优化。
优化过后：

# =============================================================================
# 定义区间长度, 子区间数, 步长
l = 1           # 区间长度
n = 8           # 子区间数
h = 1/8         # 步长


# =============================================================================
# 定义被积函数
def f_(x):
    return 4/(1+x**2)

# 构造字典，存放分点数据
table = {}

# 定义查找函数, 若有则直接返回, 若没有则先计算, 存入字典中, 再返回
def f(x):
    if x in table.keys():
        return table[x]
    else:
        res = f_(x)
        table[x] = res
        return res

# =============================================================================
# 定义复化梯形公式
def T(s, t, n):
    res = f(s)+f(t)
    len = t-s   # 区间总长度
    h = len/n   # 子区间长度
    for i in range(1, n):
        res += 2*f(i*h)
    res*=(1/2 * h)
    return res


# 定义复化辛普森公式
def S(s, t, n):
    res = f(s)+f(t)
    len = t-s   # 区间总长度
    h = len/n   # 子区间长度
    
    for i in range(1,n):
        res += 2*f(i*h)
    for i in range(0,n):
        res += 4*f(i*h+h/2)

    res*=(1/6 * h)
    return res

# 定义复化科特斯公式
def C(s, t, n):
    res = 7*f(s)+7*f(t)
    len = t-s   # 区间总长度
    h = len/n   # 子区间长度
    
    for i in range(0,n):
        res += 32*f(i*h+1/4*h)
        res += 32*f(i*h+3/4*h)
        res += 12*f(i*h+1/2*h)
    for i in range(1,n):
        res += 14*f(i*h)
    
    res*=(1/90 * h)
    return res


def main():
    print('复化梯形公式得:')
    print(T(0, 1, 8))
    print('\n')

    print('复化辛普森公式得:')
    print(S(0, 1, 4))
    print('\n')
    
    print('复化柯特斯公式得:')
    print(C(0, 1, 2))

if __name__ == '__main__':
    main()

2. 用龙贝格求积公式计算积分

   $$ I = \int^1_0\frac{4}{1+x^2}dx $$

解：

# =============================================================================
# 定义区间长度, 子区间数, 步长
l = 1           # 区间长度
n = 8           # 子区间数
h = 1/8         # 步长


# =============================================================================
# 定义被积函数
def f_(x):
    return 4/(1+x**2)

# 构造字典，存放分点数据
table = {}

# 定义查找函数, 若有则直接返回, 若没有则先计算, 存入字典中, 再返回
def f(x):
    if x in table.keys():
        return table[x]
    else:
        res = f_(x)
        table[x] = res
        return res
    
    
    
# =============================================================================
# 定义查找梯形公式值, 若有则直接返回, 若没有则先计算, 存入字典中, 再返回
T_tab = {}
S_tab = {}
C_tab = {}
R_tab = {}
def T_(s, t, x):
    if x in T_tab.keys():
        return T_tab[x]
    else:
        res = T(s, t, x)
        T_tab[x] = res
        return res


def S_(s, t, x):
    if x in S_tab.keys():
        return S_tab[x]
    else:
        res =  ( 4*T_(s, t, 2*x) - T_(s, t, x) )/3
        S_tab[x] = res
        return res


def C_(s, t, x):
    if x in C_tab.keys():
        return C_tab[x]
    else:
        res =  ( 16*S_(s, t, 2*x) - S_(s, t, x) )/15
        C_tab[x] = res
        return res


def R_(s, t, x):
    if x in R_tab.keys():
        return R_tab[x]
    else:
        res =  ( 64*C_(s, t, 2*x) - C_(s, t, x) )/63
        R_tab[x] = res
        return res

# =============================================================================
# 定义复化梯形公式
def T(s, t, n):
    res = f(s)+f(t)
    len = t-s   # 区间总长度
    h = len/n   # 子区间长度
    for i in range(1, n):
        res += 2*f(i*h)
    res*=(1/2 * h)
    return res


def main():
    err = 1
    i = 0
    while(err>0.000001):
        err = R_(0, 1, i+1)-R_(0, 1, i)
        i+=1
    print(R_(0, 1, i))
        
if __name__ == '__main__':
    main()

常微分方程数值解法

1. 用欧拉公式求解初值问题

   $$ y'=y-\frac{2x}{y},0<x\leqslant1 $$

   $$ y(0)=1 $$

   取步长$h=0.2$。

解：参数代入公式，求解该问题的欧拉公式是：

$$ y_{i+1}=y_i+h{y_i-\frac{2x_i}{y_i}} $$

$$ y_0=1 $$

编程如下：

# =============================================================================
# 定义该问题的常量
h = 0.2         # 步长
prec_res = 6        # 结果的精度

# =============================================================================
# 定义公式
def Euler(x, y):
    return (y + h*(y-2*x/y))

def exact(x):
    return (1+2*x)**(1/2)

# =============================================================================
# 用于遍历的量
k  = 0
x  = 0.0           # 0 < x <= 1
y  = round(1,prec_res)           # y0=1
ye = (exact(x))    # 精确值
err= abs(y-ye)

print('k', '\t', 'x', '\t\t', 'y', '\t\t', 'ye', '\t\t', 'err')
print(k, '\t', x, '\t', y, '\t\t', ye, '\t\t', err)
err = 100
while x<1:
    k  += 1
    y  =  Euler(x, y)
    x  = round(x+h, 1)
    ye =  exact(x)
    err=  abs(y-ye)
    print(k, '\t', x, '\t', round(y, prec_res), '\t', round(ye, prec_res), '\t', round(err, prec_res))

自己输出的结果：

用高等数学的微分方程解法，可以解得该题的精确解是$y=\sqrt{1+2x}$，跟程序输出相差不大。

2. 用改进欧拉公式对上式求值

解：改进欧拉公式其实就是转换运算方向，多求了一个后退的欧拉公式，然后取了平均。附上代码：

# =============================================================================
# 定义该问题的常量
h = 0.2         # 步长
prec_res = 6        # 结果的精度

# =============================================================================
# 定义公式

# 欧拉冲冲冲公式
def Eulerp(x, y):
    return (y + h*(y-2*x/y))

# 欧拉且战且退公式
def Eulerc(x, yi, yp):
    return (yi + h*(yp-2*x/yp))

# 微分方程精确值
def exact(x):
    return (1+2*x)**(1/2)


# =============================================================================
# 用于遍历的量
k  = 0
x  = 0.0           # 0 < x <= 1
y  = round(1,prec_res)           # y0=1
ye = (exact(x))    # 精确值
err= abs(y-ye)

print(k, '\t', x, '\t', y, '\t', ye, '\t', err)
err = 100
while x<1:
    k  += 1
    yp = Eulerp(x, y)
    x  = round(x+h, 1)
    yc = Eulerc(x, y, yp)
    y=(yp+yc)/2
    ye =  exact(x)
    err=  abs(y-ye)
    print(k, '\t', x, '\t', round(y, prec_res), '\t', round(ye, prec_res), '\t', round(err, prec_res))

附上输出结果：