Python ODE45索引器错误：列表分配超出范围

Question

我试图复制一个ODE脚本，我在Matlab中运行的Python。以下是Matlab脚本：

t0 = 0;
tfinal = 25;
q1 = 1;
q2 = 1;
q1dot = 0;
q2dot = 0;

% ODE variables
times = [t0 tfinal];
stateVars=[q1 q1dot q2 q2dot];

% ODE45
options = odeset('reltol',1E-12,'abstol',1E-12);
[t,xx]=ode45('hw1SS',times,stateVars,options);

hw1SS功能：

function [xdot] = hw1SS(t,x)
    % Given Parameters
    m = 1;
    K = 49;
    k = 1;
    C = 2;
    c = 0.5;
    k = 1;

    % Derivative variables solved for in part b)
    xdot = zeros(4,1);
    xdot(1) = x(2);
    xdot(2) = (1 / m) * (-c * x(2) - k * x(1) + C * (x(4) - x(2)) + K * (x(3) - x(1)));
    xdot(3) = x(4);
    xdot(4) = (1 / m) * (-c * x(4) - k * x(3) + C * (x(2) - x(4)) + K * (x(1) - x(3)));
end

Matlab代码运行得非常好。当我尝试在Python中复制它时，会收到以下错误：IndexError:list assignment index out of range

下面是我尝试的Python解决方案：

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

def model(x,t):
    m, K, C, c, k = 1, 49, 2, 0.5, 1
    xdot = [[],[]]
    xdot[0] = x[1]
    xdot[1]  = (1/m)*(-c*x[1]-k*x[0]+C*(x[3]-x[1])+k*(x[2]-x[0]))
    xdot[2] = x[3]
    xdot[3]  = (1/m)*(-c*x[3]-k*x[2]+C*(x[1]-x[3])+k*(x[0]-x[2]))
    return xdot
q1, q1dot, q2, q2dot = 1, 0, 1, 0
t0, tf = 0, 25
t = np.linspace(t0, tf, 100)
y0 = [q1, q1dot, q2, q2dot]
y = odeint(model, t, y0)

plt.plot(t, y)
plt.show()

错误消息：

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-61-1038b9266c62> in <module>()
      3 y0 = [q1, q1dot, q2, q2dot]
      4 
----> 5 y = odeint(model, t, y0)
      6 
      7 plt.plot(t, y)

~\Anaconda3\lib\site-packages\scipy\integrate\odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
    231                              full_output, rtol, atol, tcrit, h0, hmax, hmin,
    232                              ixpr, mxstep, mxhnil, mxordn, mxords,
--> 233                              int(bool(tfirst)))
    234     if output[-1] < 0:
    235         warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."

<ipython-input-53-3e8f74e5a043> in model(x, t)
      5     xdot[0] = x[1]
      6     xdot[1]  = (1/m)*(-c*x[1]-k*x[0]+C*(x[3]-x[1])+k*(x[2]-x[0]))
----> 7     xdot[2] = x[3]
      8     xdot[3]  = (1/m)*(-c*x[3]-k*x[2]+C*(x[1]-x[3])+k*(x[0]-x[2]))
      9 

IndexError: list assignment index out of range

Answer 1

首先，函数ode45与函数^{}在参数的位置问题上不同：

scipy.integrate.odeint(func, y0, t, args=(), ...)

正如我们看到的，第二个参数是初始值，而不是时间范围。你知道吗

另一个错误是：

xdot = [[],[]]

这是一个列表，您应该创建一个包含4个元素的列表，其中的任何值都将被覆盖（通常这些值是0）

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


def model(x, t):
    m, K, C, c, k = 1, 49, 2, 0.5, 1
    xdot = [0, 0, 0, 0]
    xdot[0] = x[1]
    xdot[1]  = (1/m)*(-c*x[1]-k*x[0]+C*(x[3]-x[1])+k*(x[2]-x[0]))
    xdot[2] = x[3]
    xdot[3]  = (1/m)*(-c*x[3]-k*x[2]+C*(x[1]-x[3])+k*(x[0]-x[2]))
    return xdot

q1, q1dot, q2, q2dot = 1, 0, 1, 0
t0, tf = 0, 25
t = np.linspace(t0, tf, 100)
y0 = [q1, q1dot, q2, q2dot]
y = odeint(model, y0, t)
plt.plot(t, y)
plt.show()