python 实现GM21

# -*- coding=utf-8 -*-
import numpy as np
from sympy import *
from sympy.abc import x, y
import re
import math

init_printing()
# 定义符号常量x 与 f(x) g(x)。这里的f g还可以用其他字母替换，用于表示函数
f, g = symbols(‘f g‘, cls=Function)


def solving_equation(x1, equation_parameter, n):
    # 二阶齐次微分方程的根分为两个不同的实根、两个相同的实根以及两个虚根
    # 不同情况方程的形式不同，根据predict函数中求得的带参数的方程来决定使用的策略
    # 下面以两个不同的实根为例
    parameter = solve(
    [x * math.exp(equation_parameter[0] * 0) + y * math.exp(equation_parameter[1] * 0) + equation_parameter[2] - x1[0],
    x * math.exp(equation_parameter[0] * (len(x1)-1)) + y * math.exp(equation_parameter[1] * (len(x1)-1)) + equation_parameter[2] - x1[len(x1) - 1]])
    print("parameter", parameter)
    # 返回X1的预测值
    return [parameter[x] * math.exp(equation_parameter[0] * i) + parameter[y] * math.exp(equation_parameter[1] * i) + equation_parameter[2] for i in range(len(x1) + n)]


# n代表需要预测的个数，默认为0
def predict(data,n=0):
    x1 = data.cumsum()
    a_x0 = np.ediff1d(data).T
    z = (x1[:len(x1) - 1] + x1[1:]) / 2.0
    B = np.array([-data[1:], -z, np.ones([len(data) - 1])]).T
    Y = a_x0
    u = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)
    a1, a2, b = u[0], u[1], u[2]
    # 用diffeq代表微分方程
    diffeq = Eq(f(x).diff(x, x) + a1 * f(x).diff(x) + a2 * f(x), b)
    # 调用dsolve函数,返回一个Eq对象,并提取带参数方程
    differential_equation = str(dsolve(diffeq, f(x)).args[1])
    # 使用正则表达式提取齐次微分方程的根与非齐次微分方程的特解
    equation_parameter = re.findall("-?\d+.?\d+", differential_equation.replace(‘ ‘, ‘‘))
    # str转为float
    for i in range(len(equation_parameter)):
        equation_parameter[i] = float(equation_parameter[i])
    # 利用边界条件，取X1中第一个数和最后一个数,构造方程组，求参数C1和C2，并返回预测值
    return solving_equation(x1, equation_parameter, n)


if __name__ == ‘__main__‘:
    data = np.array([41, 49, 61, 78, 96, 104])
    predict_data = predict(data)  # 预测x1
    result = np.ediff1d(predict_data)  # 递减，求出最终结果
    print(‘原数据:‘, data[1:])
    print(‘预测结果:‘, result)