数学表达式：
x
¨
+
δ
x
˙
+
α
x
+
β
x
3
=
γ
c
o
s
(
w
t
)
ddot{x}+deltadot{x}+alpha x+beta x^3 = gamma cos(wt)
x¨+δx˙+αx+βx3=γcos(wt)
x
1
=
x
˙
x_1 = dot{x}
x1=x˙
x
2
=
x
¨
x_2 = ddot{x}
x2=x¨
我们可以变化得到：
x
1
=
x
˙
x_1 = dot{x}
x1=x˙
x
2
=
−
δ
x
˙
−
α
x
−
β
x
3
+
γ
c
o
s
(
w
t
)
x_2 = -delta dot{x} -alpha x - beta x^3 + gamma cos(wt)
x2=−δx˙−αx−βx3+γcos(wt)
matlab实现
1.写一个函数
function y = duffineq(t,x) global delta alpha beta gamma omega x1 = x(2); x2 = -delta * x(2) - alpha * x(1) - beta * x(1)^3 + gamma * cos(omega * t); y = [x1;x2]; end
2.画图
global delta alpha beta gamma omega
delta = 1;%阻尼系数,(Ns/m or kg/s)
alpha = 1;%线性刚度,(N/m)
beta = 1;%非线性刚度,(N/m^3)
gamma = 1;%驱动力的大小,(N)
omega = 1;%驱动力的角速度,(rad/s)
fsize = 18;%字号大小
x0 = [0.1 0]; %初始状态 [初始位移，初始速度]
tspan = 0:0.001:100; %仿真时长
[t, y] = ode45(@duffineq, tspan, x0);
figure
plot(y(:,1),y(:,2)) %第一列y(:,1)是x的位移,第二列y(:,2)是x的速度。
xlabel('$it x$','Interpreter','latex','FontSize',fsize)
ylabel('$dot{x}$', 'Interpreter','latex','FontSize',fsize)
3.运行结果
Python实现
1.代码
import numpy as np
import math
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# function
def model(x, t):
delta = 1
alpha = 1
beta = 1
gamma = 1
omega = 1
x1 = x[1]
x2 = -delta * x[1] - alpha * x[0] - beta * x[0] ** 3 + gamma * math.cos(omega * t)
return [x1, x2]
# initial condition
x0 = [0.1, 0]
# time
t = np.linspace(0, 100, 100000) # 起始点，终止点，个数。
# solve ode
y = odeint(model, x0, t)
# plot
plt.plot(y[:, 0], y[:, 1])
plt.xlabel('x')
plt.ylabel(r"$. x $") # $. x $ 为x头上一点
plt.title(r"x vs. $. x $")
plt.show()
2.运行结果
