欧拉法与四阶龙格库塔方法

欧拉法

• 考虑如下系统 [1]
  
• 欧拉迭代格式
  
  其中$B_1$为标准正态分布随机数

Python 实现

import time
import sympy
import numpy as np
import pandas as pd
import seaborn as sns
from scipy import stats
import matplotlib as mpl
from matplotlib import cm
import matplotlib.pyplot as plt
from scipy.special import erfc
from scipy.integrate import quad
from scipy.optimize import root,fsolve
from mpl_toolkits.mplot3d import Axes3D
from multiprocessing.pool import ThreadPool
from matplotlib.font_manager import FontProperties
from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error
from scipy import special

import warnings
warnings.filterwarnings('ignore')

mpl.rcParams['font.sans-serif'] = 'Times New Roman'
mpl.rcParams['axes.unicode_minus'] = 'True'
mpl.rcParams['font.size'] = 15

plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['font.sans-serif'] = 'Times New Roman'
plt.rcParams['axes.unicode_minus'] = 'True'
# 用tex公式的形式输入英文和公式，以显示Times New Roman字体
plt.rcParams['mathtext.fontset'] = 'stix'

plt.rcParams['figure.dpi'] = 600 #分辨率

'----------parameter seeting----------'
delta_t, N, M, epsilon, delta = 0.001, 10**7, 3, 1, 5
c0, c2 = -1, 0.5
mu1, mu3 = 1, -2
theta1,theta2 = 0.5, 0.5
a1, a2 = 1, 1
lambda1, zeta1 = 1, 1
lambda2, zeta2 = 0.5, 1

b1 = lambda y : y 
b2 = lambda x,y,v,i : -( 
    (c0 + c2*np.power(x,2)) * y + \
    (mu1*x + mu3*np.power(x,3)) + theta1*v + theta2*i
)
a_x = lambda x : a1 + a2*x
b3 = lambda v,y : -lambda1*v + zeta1*y
b4 = lambda i,y : -lambda2*i + zeta2*y


rng = np.random.default_rng(10)
#  Z初值生成
Z0 = rng.uniform(-3,3,4 * 10**7).reshape(4,10**7)

'Gaussian white noise production'
GWN = rng.standard_normal(N)

b_1 = Z0[1,:] * delta_t
b_2 = b2(Z0[0,:],Z0[1,:],Z0[2,:],Z0[3,:]) * delta_t + \
        np.power(delta_t,1/2) * a_x(Z0[0,:]) * GWN
b_3 = b3(Z0[2,:],Z0[1,:]) * delta_t
b_4 = b4(Z0[3,:],Z0[1,:]) * delta_t

b_ = np.array((b_1,b_2,b_3,b_4))
Z1 = Z0 + b_
Z1.shape

[1] Y.-H. Sun et al. Identification of hybrid energy harvesting systems with non-Gaussian process, Acta Mech. Sin. 39, 523154 (2023)

四阶龙格库塔方法

• 考虑如下系统 [2]
  
• 参数取值
  
• 迭代格式可于csdn.net中搜索

Python实现

import time
import sympy
import numpy as np
import pandas as pd
import seaborn as sns
from scipy import stats
import matplotlib as mpl
from matplotlib import cm
import matplotlib.pyplot as plt
from scipy.special import erfc
from scipy.integrate import quad
from scipy.optimize import root,fsolve
from mpl_toolkits.mplot3d import Axes3D
from multiprocessing.pool import ThreadPool
from matplotlib.font_manager import FontProperties
from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import LinearRegression,Ridge
from pysindy.feature_library import PolynomialLibrary,FourierLibrary

import warnings
warnings.filterwarnings('ignore')

mpl.rcParams['font.sans-serif'] = 'Times New Roman'
mpl.rcParams['axes.unicode_minus'] = 'True'
mpl.rcParams['font.size'] = 15

plt.rcParams['font.sans-serif'] = 'Times New Roman'
plt.rcParams['axes.unicode_minus'] = 'True'
# 用tex公式的形式输入英文和公式，以显示Times New Roman字体
plt.rcParams['mathtext.fontset'] = 'stix'

plt.rcParams['figure.dpi'] = 600 #分辨率

def MotionEquation(t,Y,Wt,c0, c2, c4,delta1, delta3, delta5,zeta1,zeta2,F0,Omega,omega,lambda1,mu1,lambda2,mu2):
    dy1dt = Y[1]
    dy2dt = -( 
        (c4*np.power(Y[0],4)+c2*np.power(Y[0],2)+c0)*Y[1] + \
        (delta1*Y[0]+delta3*np.power(Y[0],3)+delta5*np.power(Y[0],5))+zeta1*Y[2]+zeta2*Y[3]
    ) + F0*np.cos(Omega*t/omega) + Wt
    dy3dt = -lambda1*Y[2] + mu1*Y[1]
    dy4dt = -lambda2*Y[3] + mu2*Y[1]
    dYdt = np.array([dy1dt,dy2dt,dy3dt,dy4dt])
    return dYdt

'----------setting params----------'
c0, c2, c4 = -0.5, 0.5, -0.1
delta1, delta3, delta5 = 1, -3, 1
zeta1, zeta2 = 0.5, 0.5
F0, Omega, omega, gamma = 1, 1, 1, 0
lambda1, mu1 = 1, 1
lambda2, mu2 = 0.5, 1

t0, tn, h = 0, 100, 0.001
N = int( (tn-t0)/h )
t = np.arange(t0,tn+h,h)
Z = np.zeros( (4,N+1) )
Z[:,0] = np.array([0.01, 0.01, 0.01, 0.01])

rng = np.random.default_rng(13)
Wt = np.sqrt(gamma/h)*rng.standard_normal(N)
for i in range(N):
    k1 = MotionEquation(
        t[i],Z[:,i],Wt[i],c0, c2, c4,delta1, delta3, delta5,zeta1,zeta2,
        F0,Omega,omega,lambda1,mu1,lambda2,mu2
    )
    k2 = MotionEquation(
        t[i]+h/2,Z[:,i]+h*k1/2,Wt[i],c0, c2, c4,delta1, delta3, delta5,zeta1,zeta2,
        F0,Omega,omega,lambda1,mu1,lambda2,mu2
    )
    k3 = MotionEquation(
        t[i]+h/2,Z[:,i]+h*k2/2,Wt[i],c0, c2, c4,delta1, delta3, delta5,zeta1,zeta2,
        F0,Omega,omega,lambda1,mu1,lambda2,mu2
    )
    k4 = MotionEquation(
        t[i]+h,Z[:,i]+h*k3,Wt[i],c0, c2, c4,delta1, delta3, delta5,zeta1,zeta2,
        F0,Omega,omega,lambda1,mu1,lambda2,mu2
    )
    Z[:,i+1] = Z[:,i]+h*(k1+2*k2+2*k3+k4)/6

print(t.shape,Z.shape)

[2] Ya-Hui Sun, et al. Sparse identification method of extracting hybrid energy harvesting system from observed data. Chin. Phys. B, 2022, 31 (12): 120203.