Monte Carlo Simulation(MCS)的使用

MCS估计二维联合概率密度离散值

考虑如下非线性动力系统：

• 参数值为：
  $m,c,k1,k3,D=1,0.5,1.5,2.5,1/2$

Python 代码如下：

代码中不需要导入以下所有依赖库，因方便，于是复制粘贴所有常用库。

import time
import sympy
import numpy as np
import pandas as pd
import seaborn as sns
import scipy
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
from pysindy.feature_library import PolynomialLibrary,FourierLibrary
import torch

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,c,k1,k3):
    dxdt = y[1]
    dydt = -(c*y[1]+k1*y[0]+k3*np.power(y[0],3)) + Wt
    return np.array([dxdt,dydt])

def RK45(m, c, k1, k3, D):
    delta_t, N = 0.05, 1000000
    t = np.arange(0,N*delta_t,delta_t)
#     rng = np.random.Generator(np.random.MT19937(24))
    rng = np.random.default_rng(24)
    GWN = np.sqrt(2*D/delta_t)*rng.standard_normal(N)
#     np.random.seed(13)
#     GWN = np.sqrt(2*D/delta_t)*np.random.randn(N)
#     path = r'E:\matlb\bin\05_2021YanWork\20230724-SystemReduction\RamdonNumber.mat'
#     GWN = scipy.io.loadmat(path)['Wt']
#     GWN = GWN.reshape(GWN.shape[1])

    y = np.zeros((2,N+1))
    y[:,0] = np.array([0.01,0.01])
    for i in range(N):
        rk1 = MotionEquation(
            t[i],y[:,i],GWN[i],c,k1,k3
        )
        rk2 = MotionEquation(
            t[i]+delta_t/2,y[:,i]+rk1*delta_t/2,GWN[i],c,k1,k3
        )
        rk3 = MotionEquation(
            t[i]+delta_t/2,y[:,i]+rk2*delta_t/2,GWN[i],c,k1,k3
        )
        rk4 = MotionEquation(
            t[i]+delta_t,  y[:,i]+rk3*delta_t,  GWN[i],c,k1,k3
        )
        y[:,i+1] = y[:,i] + (rk1+2*rk2+2*rk3+rk4)*delta_t/6
    return t,y

def BasisFunc(x,xdot,m, c, k1, k3, D):
#         m, c, k1, k3, D  = 1, 0.5, 1.5, 2.5, 0.5
    dictionary = np.vstack([
        np.power(x,0), np.power(x,2), np.power(x,4), x*xdot,
        m*k3/D*np.power(x,3)*xdot,m*c*k1*k3/np.power(D,2)*np.power(x,5)*xdot,
        m*c/D*np.power(xdot,2), m*c*k3/(k1*D)*np.power(x,2)*np.power(xdot,2)
    ])
    return dictionary.T

m, c, k1, k3, D = 1, 0.5, 1.5, 2.5, 1/2
t,y = RK45(m, c, k1, k3, D)
print(t.shape,y.shape)

x1,x1dot = y[0,1:],y[1,1:]
x1_x1dot = np.vstack([x1,x1dot])

delta_x,delta_xdot = 0.05, 0.05
MCS_x1, MCS_x2 = np.arange(-4,4,delta_x), np.arange(-4,4,delta_xdot)

frequency = np.zeros((MCS_x1.shape[0],MCS_x2.shape[0]))
for i1 in range(MCS_x1.shape[0]):
    for j1 in range(MCS_x2.shape[0]):
        i,j = MCS_x1[i1], MCS_x2[j1]
        idx = np.array(
            [(x1_x1dot[0,:]>i-delta_x), (x1_x1dot[0,:]<i+delta_x), 
             (x1_x1dot[1,:]>j-delta_xdot), (x1_x1dot[1,:]<j+delta_xdot)]
        )
        frequency[i1,j1] = np.count_nonzero(np.all(idx,axis=0))
        
pdf_x1_x1dot = frequency/(x1.shape[0]*delta_x*2*delta_xdot*2)
MCS_x1_, MCS_x2_ = np.meshgrid(MCS_x1, MCS_x2)
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(
    MCS_x1_,MCS_x2_,pdf_x1_x1dot,
    cmap='Spectral_r'
)
plt.show()
plt.close()

运行结果与论文对比

[1] Tian, YP et al. Stationary response probability density of nonlinear random
vibrating systems: a data-driven method. Nonlinear Dynamics (2020) 100:2337–2352