使用 odeint 求解 ODE 时如何传递源项答案

【问题标题】：How to pass the source term when solving ODE using odeint使用 odeint 求解 ODE 时如何传递源项
【发布时间】：2017-12-18 11:04:49
【问题描述】：

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

然后我使用“odeint”，它通过调用“diff”函数和其他参数来求解微分方程。现在，如果我将 F=0 放在“diff”函数中（即忽略任何 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: setting an array element with a sequence。”

如何解决问题？

代码：

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()

【问题讨论】：

标签： python ode

【解决方案1】：

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

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

还请注意，您在t 到odeint 中给出的采样时间列表（几乎）与odeint 调用ODE 函数diff 的时间无关。如果你想控制你必须实现自己的固定步法，可能不是龙格库塔，而是一些多步法。

【讨论】：