我一直在为一个学校项目编写隐式euler方法。 它是关于一个钟摆在重力作用下的摩擦力。 因此,方程被一分为二
def pend(y, t, b, c):
theta, omega = y
dydt = [omega, -b*omega - c*np.sin(theta)]
return dydt
b
,c
是常量。
这是欧拉方法
def Euler(f,tinit,tfinal, THinit,N,b,c):
h= (tfinal-tinit)/N #step
TH = [THinit] #initialization TH = [theta,omega]
lt = [tinit]
t=tinit
y=THinit
while t<tfinal:
t= t+h
lt.append(t) #list of time
result =scipy.optimize.fsolve(lambda x:x-y-h*f(x,t,b,c), y)
TH.append(r)
y=result
return TH,lt
将此代码用于f = pend
,我得到了错误
TypeError: can't multiply sequence by non-int of type 'float'
也许我们不能对具有数组参数和/或返回数组的函数使用fsolve
。谢谢你的帮助。你知道吗
完整代码:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
g=9.81
L=0.5
def pend(y, t, b, c):
theta, omega = y
dydt = [omega, -b*omega - c*np.sin(theta)]
return dydt
def Euler(f,tinit,tfinal, THinit,N,b,c):
h= (tfinal-tinit)/N
TH = [THinit]
lt = [tinit]
t=tinit
y=THinit
while t<tfinal:
t= t+h
lt.append(t)
result =scipy.optimize.fsolve(lambda x:x-y-h*f(x,t,b,c), y)
TH.append(r)
y=result
return TH,lt
Y,X = Euler(pend,0,10,[np.pi/8,0],100,0.25,5)
plt.plot(X,Y)
plt.show()
这是您看到的错误:
原因是列表上
*
操作符的语义。你知道吗你可能想用一个向量乘以一个标量
因此pend应该返回
numpy.array(dydt)
您的
pend()
返回一个列表(return [omega, -b*omega - c*np.sin(theta)]
通过变量dydt
)。这个列表然后用在术语h*f(x,t,b,c)
中。h
是float
,所以您尝试将float
和list
相乘。这是不可能的。将int
与list
相乘是可能的,并导致列表的重复:结果
但对于浮动,这还没有定义。请提供更多关于你期望发生什么的信息。只给出(这种)代码是很难猜测的。你知道吗
编辑:
根据后面添加的内容,我现在猜您希望将数组中的每个元素乘以float,并期望得到一个数组作为结果。这就是
numpy
数组所做的,因此我建议立即返回numpy
数组:相关问题 更多 >
编程相关推荐