如何使用odeint解决ODE时传递源项

Question

强迫谐振子的微分方程给出为Mx''+ Lx'+（w ^ 2）x = F（t）。这里F（t）是源项。为了解决这个问题，我写了一个代码，在函数'diff'中定义微分方程。我写了另一个函数'generate_pulse'给出了F（t）。

然后我使用'odeint'，它通过调用'diff'函数和其他参数来解决微分方程。如果我在'diff'函数中放置F = 0（即忽略任何F（t）项），我现在没有得到任何错误消息。请查看'diff'功能：

F=0         #No error detected if I put F=0 here. Comment out this line to see the error

一旦我保留了F（t），我收到一条错误消息'ValueError：设置一个带序列的数组元素'。

如何解决问题？

码：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
import itertools



def diff(y, t,param):
    M=param[0]
    L=param[1]
    w2=param[2]
    F=param[3]
    F=0         #No error detected if I put F=0 here. Comment out this line to see the error
    x,v = y
    dydt = [v, (F-L*v - w2*x)/M]
    return dydt


def generate_pulse(t,Amp,init_delay,rtime,PW):


        L=len(t)
        t0=t[0:math.ceil(init_delay*L/100)]  
        t1=t[len(t0):len(t0)+math.ceil(rtime*L/100)]
        t2=t[len(t0)+len(t1):len(t0)+len(t1)+math.ceil(PW*L/100)]
        t3=t[len(t0)+len(t1)+len(t2):len(t0)+len(t1)+len(t2)+2*math.ceil(rtime*L/100)]
        t4=t[len(t0)+len(t1)+len(t2)+len(t3):len(t0)+len(t1)+len(t2)+len(t3)+math.ceil(PW*L/100)]
        t5=t[len(t0)+len(t1)+len(t2)+len(t3)+len(t4):len(t0)+len(t1)+len(t2)+len(t3)+len(t4)+math.ceil(rtime*L/100)]
        t6=t[len(t0)+len(t1)+len(t2)+len(t3)+len(t4)+len(t5):]

        s0=0*t0

        s1=(Amp/(t1[-1]-t1[0]))*(t1-t1[0])
        s2=np.full((1,len(t2)),Amp)
        s2=list(itertools.chain.from_iterable(s2))  #The 'tuple' is converted into array

        s3=-Amp/(t1[-1]-t1[0])*(t3-t3[0])+Amp

        s4=np.full((1,len(t4)),-Amp)
        s4=list(itertools.chain.from_iterable(s4))  #The 'tuple' is converted into array

        s5=(Amp/(t5[-1]-t5[0]))*(t5-t5[0])-Amp
        s6=np.full((1,len(t6)),0)
        s6=list(itertools.chain.from_iterable(s6))  #The 'tuple' is converted into array

        s=[s0,s1,s2,s3,s4,s5,s6]
        s=list(itertools.chain.from_iterable(s))
        return s

###############################################################################
#                      Main code from here

t = np.linspace(0, 30, 200)
y0 = [- 10, 0.0]


M=5
L = 1
w2 = 15.0

Amp=5
init_delay=10
rtime=10
PW=10

F=generate_pulse(t,Amp,init_delay,rtime,PW)


Param=([M,L,w2,F],)   #Making the 'Param' a tuple. Because args of odeint takes tuple as argument. 


sol = odeint(diff, y0, t, args=Param)


plt.plot(t, sol[:, 0], 'b', label='x(t)')
plt.plot(t,F,'g',label='Force(t)')
plt.legend(loc='best')
plt.show()

Answer 1

您收到错误，因为您传递的F的值是您生成的数组。

使用numpy或scipy的插值函数从数组中生成实际函数。然后在时间t评估该函数。或者直接将强制项实现为标量t的分段定义函数。

另请注意，您在t中向odeint提供的采样时间列表（几乎）与odeint调用ODE函数diff的时间无关。如果你想控制你必须实现自己的固定步骤方法，可能不是Ru​​nge-Kutta而是一些多步骤方法。