Issue understanding scipy.integrate.RK45 requirements

Question

I am attempting to solve a system of first order differential equations with scipy.integrate.RK45(). I have scripted out the model function that I am hoping to plot (displacement vs time), but RK45() requires that this function takes 2 arguments, namely 't' and 'y', where 't' is a scalar and 'y' is an array in my case. This is better described below:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html

My script is as follows:

import numpy as np
from scipy.integrate import RK45, RK23

def model(t, y):
    
    m = 2
    c = 10
    k = 1500

    if (t >= 0 and t < 0.1):
        F = 200 * t
    if (t >= 0.1 and t < 0.25):
        F = 20
    else:
        F = 0
    
    E_matrix = [[0, 1], [(-k / m), (-c / m)]]
    Q_matrix = [0, F / m]

    return E_matrix * X + Q_matrix

time_step = 0.01
t_upper = 0.5
t_lower = 0

initial_conditions = [0, 0]

points_to_plot = RK45(fun=model(t, y), t0=t_lower, y0=initial_conditions, t_bound=t_upper, vectorized=True)

A picture of the system I'd like to solve is seen below:

I have found very few examples of this being used as most solutions use odeint().

What are these two arguments (t, y) and how would I effectively incorporate them into my function?

Answer 1

You have used t already. Now, change def model(t, y): to def model(t, X): and you will use X as well. Note that t and y are positional arguments, you can call them however you like in your function. You have another issue, which is that you multiply Python lists! In Python, as opposed to Matlab, you need to specify that you make an array:

Change

E_matrix = [[0, 1], [(-k / m), (-c / m)]]
Q_matrix = [0, F / m]

to

E_matrix = np.array([[0, 1], [(-k / m), (-c / m)]])
Q_matrix = np.array([0, F / m])

and

return E_matrix*X + Q_matrix

to

return E_matrix @ X + Q_matrix

as @ is the matrix product in NumPy. * performs element-wise product.

EDIT: I hadn't spotted the call RK45(fun=model(t, y),. By doing that, you are passing the value of the function model at t,y. You need to give the function itself: RK45(fun=model, ...