python - 使用数组的欧拉隐式方法

Question

我一直在尝试为学校项目编写隐式欧拉方法。这是关于一个受到重力和摩擦力作用的钟摆。因此方程被分成 2

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)