用 Python 数值求解 ODE [英] Solving ODE numerically with Python

问题描述

我正在用 Python 数值求解谐波振荡器的常微分方程.当我添加驱动力时，它没有任何区别，所以我猜代码有问题.任何人都可以看到问题吗?(h/m)*f0*np.cos(wd*i)部分是驱动力.

I am solving an ODE for an harmonic oscillator numerically with Python. When I add a driving force it makes no difference, so I'm guessing something is wrong with the code. Can anyone see the problem? The (h/m)*f0*np.cos(wd*i) part is the driving force.

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()

推荐答案

这里的问题在于术语 np.cos(Wd*i).它应该是 np.cos(Wd*i*dT)，注意 dT 已被添加到正确的方程中，因为 t = i*dT.

The problem here is with the term np.cos(Wd*i). It should be np.cos(Wd*i*dT), that is note that dT has been added into the correct equation, since t = i*dT.

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

If this correction is made, the simulation looks reasonable. Here's a version with F0=0.001. Note that the driving force is clear in the continued oscillations in the damped condition.

原方程的问题在于np.cos(Wd*i)只是绕着圆随机跳跃，而不是平滑地绕着圆移动，最终没有产生任何效果.通过直接绘制可以最好地看出这一点，但最简单的方法是运行带有 F0 非常大的原始表单.下面是F0 = 10(即正确方程中使用的值的10000倍)，但使用了错误的方程形式，很明显，这里的驱动力只是在随机移动时增加了噪音圆圈.

The problem with the original equation is that np.cos(Wd*i) just jumps randomly around the circle, rather than smoothly moving around the circle, causing no net effect in the end. This can be best seen by plotting it directly, but the easiest thing to do is run the original form with F0 very large. Below is F0 = 10 (ie, 10000x the value used in the correct equation), but using the incorrect form of the equation, and it's clear that the driving force here just adds noise as it randomly moves around the circle.

这篇关于用 Python 数值求解 ODE的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！