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)