数值分析-常微分方程数值解

数值积分上机课写的 Python 代码，实现常微分方程数值求解。

import numpy as np


def improved_euler(f, y0, x_range, h):
    
    
    def predict(x, y):
        return y + h * f(x, y)
    
    
    def correction(x, y):
        return y + h / 2 * (f(x, y) + f(x + h, predict(x, y)))
    
    
    X = np.append(np.arange(x_range[0], x_range[1], h), x_range[1])
    Y = [y0]
    for i in range(len(X) - 1):
        Y.append(correction(X[i], Y[i]))
        
    return np.array(Y)


def fourth_order_classical_RK(f, y0, x_range, h):
    
    
    def get_K(x, y):
        K_1 = f(x, y)
        K_2 = f(x + h / 2, y + h / 2 * K_1)
        K_3 = f(x + h / 2, y + h / 2 * K_2)
        K_4 = f(x + h, y + h * K_3)
        
        return K_1, K_2, K_3, K_4
    
    
    X = np.append(np.arange(x_range[0], x_range[1], h), x_range[1])
    Y = [y0]
    for i in range(len(X) - 1):
        K = get_K(X[i], Y[i])
        Y.append(Y[i] + h / 6 * (K[0] + 2 * K[1] + 2 * K[2] + K[3]))
    
    return np.array(Y)