python - Python ODE45 IndexError: list assignment out of range

Question

I am trying to duplicate an ODE script I have running in Matlab to Python. Here is the Matlab script:

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:

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

The Matlab code works perfectly fine. When I try to replicate this in Python, I receive the following error: IndexError: list assignment index out of range

Here is my attempted solution in 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()

Error message:

---------------------------------------------------------------------------
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不同odeint：

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