用Python数值求解ODE

Question

我正在用Python数值求解谐振子的ODE。当我添加一个驱动力时，它没有区别，所以我猜代码有问题。有人看到问题了吗？(h/m)*f0*np.cos(wd*i)部分是驱动力。

import numpy as np
import matplotlib.pyplot as plt

# This code solves the ODE mx'' + bx' + kx = F0*cos(Wd*t)
# m is the mass of the object in kg, b is the damping constant in Ns/m
# k is the spring constant in N/m, F0 is the driving force in N,
# Wd is the frequency of the driving force and x is the position 

# Setting up

timeFinal= 16.0   # This is how far the graph will go in seconds
steps = 10000     # Number of steps
dT = timeFinal/steps      # Step length 
time = np.linspace(0, timeFinal, steps+1)   
# Creates an array with steps+1 values from 0 to timeFinal

# Allocating arrays for velocity and position
vel = np.zeros(steps+1)
pos = np.zeros(steps+1)

# Setting constants and initial values for vel. and pos.
k = 0.1
m = 0.01
vel0 = 0.05
pos0 = 0.01
freqNatural = 10.0**0.5
b = 0.0
F0 = 0.01
Wd = 7.0
vel[0] = vel0    #Sets the initial velocity
pos[0] = pos0    #Sets the initial position



# Numerical solution using Euler's
# Splitting the ODE into two first order ones
# v'(t) = -(k/m)*x(t) - (b/m)*v(t) + (F0/m)*cos(Wd*t)
# x'(t) = v(t)
# Using the definition of the derivative we get
# (v(t+dT) - v(t))/dT on the left side of the first equation
# (x(t+dT) - x(t))/dT on the left side of the second 
# In the for loop t and dT will be replaced by i and 1

for i in range(0, steps):
    vel[i+1] = (-k/m)*dT*pos[i] + vel[i]*(1-dT*b/m) + (dT/m)*F0*np.cos(Wd*i)
    pos[i+1] = dT*vel[i] + pos[i]

# Ploting
#----------------
# With no damping
plt.plot(time, pos, 'g-', label='Undampened')

# Damping set to 10% of critical damping
b = (freqNatural/50)*0.1

# Using Euler's again to compute new values for new damping
for i in range(0, steps):
    vel[i+1] = (-k/m)*dT*pos[i] + vel[i]*(1-(dT*(b/m))) + (F0*dT/m)*np.cos(Wd*i)
    pos[i+1] = dT*vel[i] + pos[i]

plt.plot(time, pos, 'b-', label = '10% of crit. damping')
plt.plot(time, 0*time, 'k-')      # This plots the x-axis
plt.legend(loc = 'upper right')

#---------------
plt.show()

Answer 1

注意，您的ODE表现良好，并且有一个解析解。所以你可以利用同情心来替代：

import sympy as sy    
sy.init_printing()  # Pretty printer for IPython

t,k,m,b,F0,Wd = sy.symbols('t,k,m,b,F0,Wd', real=True)  # constants

consts = {k:  0.1, # values
          m:  0.01,
          b:  0.0,
          F0: 0.01,
          Wd: 7.0}

x = sy.Function('x')(t)  # declare variables
dx = sy.Derivative(x, t)
d2x = sy.Derivative(x, t, 2)

# the ODE:
ode1 = sy.Eq(m*d2x + b*dx + k*x, F0*sy.cos(Wd*t))

sl1 = sy.dsolve(ode1, x) # solve ODE
xs1 = sy.simplify(sl1.subs(consts)).rhs # substitute constants


# Examining the solution, we note C3 and C4 are superfluous
xs2 = xs1.subs({'C3':0, 'C4':0})
dxs2 = xs2.diff(t)

print("Solution x(t) = ")
print(xs2)
print("Solution x'(t) = ")
print(dxs2)

给予

Solution x(t) = 
C1*sin(3.16227766016838*t) + C2*cos(3.16227766016838*t) - 0.0256410256410256*cos(7.0*t)
Solution x'(t) = 
3.16227766016838*C1*cos(3.16227766016838*t) - 3.16227766016838*C2*sin(3.16227766016838*t) + 0.179487179487179*sin(7.0*t)

常数C1,C2可以通过计算初始条件的x(0),x'(0)来确定。

Answer 2

这里的问题是术语np.cos(Wd*i)。应该是np.cos(Wd*i*dT)，也就是说，自从t = i*dT以来，dT已经被添加到正确的方程中。

如果进行了此校正，则模拟看起来是合理的。这是一个带有F0=0.001的版本。注意，在阻尼条件下的持续振荡中，驱动力是明显的。

原始方程的问题是np.cos(Wd*i)只是在圆周围随机跳跃，而不是在圆周围平滑移动，最终不会产生净效应。这可以通过直接绘制来最好地看到，但是最简单的事情是使用非常大的F0运行原始表单。下面是F0 = 10（即，1000x是正确方程式中使用的值），但使用了不正确的方程式形式，很明显，这里的驱动力只是在它随机绕圆移动时添加噪声。