问题描述
我是一名工程专业的学生,我试图从scipy.integrate模块中找出如何使用odeint函数(我只在MATLAB中使用过ode45)。我正在尝试以数值方式解决一个简单的二阶质量,弹簧,阻尼系统。下面是我编写的代码(特别是我正在使用Jupyter Notebook并运行最新版本的Python 3):
import numpy as np
from scipy.integrate import odeint
from matplotlib.pyplot as plt
%matplotlib inline
# Numerical solution to mx" + bx' + kx = f(t)
# Define state vector y and its derivative
def translational(x,t,m,b,k,f):
y = [x[0],x[1]] # state vector
ydot = [x[1],f -b/m*x[1] - k/m*x[0]] # derivative of state vector
return ydot
# Parameters for the system
t = np.arange(0,10,0.01)
IC = [0,0] #[x0 v0]
m = 10 # kg
b = 2 # N*s/m
k = 5 # N/m
f = 5*np.cos(10*t)
y = odeint(translational,IC,args=(m,f))
TypeError Traceback (most recent call last)
TypeError: only size-1 arrays can be converted to Python scalars
The above exception was the direct cause of the following exception:
ValueError Traceback (most recent call last)
<ipython-input-7-423018367c52> in <module>
20 k = 5 # N/m
21 f = 5*np.cos(10*t)
---> 22 y = odeint(translational,f))
~\anaconda3\lib\site-packages\scipy\integrate\odepack.py in odeint(func,y0,args,Dfun,col_deriv,full_output,ml,mu,rtol,atol,tcrit,h0,hmax,hmin,ixpr,mxstep,mxhnil,mxordn,mxords,printmessg,tfirst)
239 t = copy(t)
240 y0 = copy(y0)
--> 241 output = _odepack.odeint(func,242 full_output,243 ixpr,ValueError: setting an array element with a sequence.
对于我的一生,我无法弄清楚出了什么问题。任何帮助深表感谢!谢谢。
解决方法
f
是一个数字数组,因此f -b/m*x[1] - k/m*x[0]
也是如此,因此函数translational
的返回值不正确。
您应该做的是使用f
中函数的表达式,而不是尝试预先计算translational
的值:
def translational(x,t,m,b,k):
y = [x[0],x[1]] # state vector
f = 5*np.cos(10*t)
ydot = [x[1],f -b/m*x[1] - k/m*x[0]] # derivative of state vector
return ydot
并从f
函数调用的args
参数中删除odeint
。