带阵列的欧拉隐式方法

Question

我一直在为一个学校项目编写隐式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()

Answer 1

这是您看到的错误：

In [1]: 3.2 * [1,2,3]
                                     -
TypeError                                 Traceback (most recent call last)
<ipython-input-5-73b36d584f00> in <module>()
  > 1 3.2 * [1,2,3]

TypeError: can't multiply sequence by non-int of type 'float'

原因是列表上*操作符的语义。你知道吗

In [2]: 3 * [1,2,3]
Out[2]: [1, 2, 3, 1, 2, 3, 1, 2, 3]

你可能想用一个向量乘以一个标量

In [3]: 3 * numpy.array([1,2,3])
Out[3]: array([3, 6, 9])

In [4]: 3.2 * numpy.array([1,2,3])
Out[4]: array([ 3.2,  6.4,  9.6])

因此pend应该返回numpy.array(dydt)

Answer 2

您的pend()返回一个列表（return [omega, -b*omega - c*np.sin(theta)]通过变量dydt）。这个列表然后用在术语h*f(x,t,b,c)中。h是float，所以您尝试将float和list相乘。这是不可能的。将int与list相乘是可能的，并导致列表的重复：

5 * [ 42, 23 ]

结果

[ 42, 23, 42, 23, 42, 23, 42, 23, 42, 23 ]

但对于浮动，这还没有定义。请提供更多关于你期望发生什么的信息。只给出（这种）代码是很难猜测的。你知道吗

编辑：

根据后面添加的内容，我现在猜您希望将数组中的每个元素乘以float，并期望得到一个数组作为结果。这就是numpy数组所做的，因此我建议立即返回numpy数组：

dydt = np.array([omega, -b*omega - c*np.sin(theta)])
return dydt