Passing matrices as input in scipy.integrate.solve_ivp python

Question

I was solving a 2DOF spring-mass-damper system given below:

These are the 2 governing Equations

I have solved it in the following way:

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
m1 = 3
m2 = 5

k1 = 7
k2 = 9

c1 = 1
c2 = 2

f1 = 40
f2 = 4

start_time = 0
end_time = 60

initial_position_m_1 = 6
initial_velocity_m_1 = 0

initial_position_m_2 = 9
initial_velocity_m_2 = 4

delta_t = 0.1
def F(t, y):
    arr = np.array([
        y[1],
        (1/m1)*(f1*np.cos(3*t) - ((c1 + c1)*y[1] + (k1 + k2)*y[0]) + c2*y[3] + k2*y[2]),
        y[3],
        (1/m2)*(f2*np.sin(t**2) - c2*y[3] - k2*y[2] + c2*y[1] + k2*y[0])
    ])
    return arr

time_interval = np.array([start_time, end_time])
initial_conditions = np.array([initial_position_m_1, initial_velocity_m_1, initial_position_m_2, initial_velocity_m_2])

#######     solving the system of equations    ####
sol = solve_ivp(F, time_interval, initial_conditions, max_step = delta_t)

T = sol.t
Y = sol.y

Now, this is done by converting the 2 governing equations into 4 equations like this:

The problem with this is that I have to write each and every equation separately (as the function F)

Matlab has a way of solving it just with matrices using Ode45 function i.e. you don't have to write all the equations seperately in the function F in Matlab. You can enter the mass, stiffness and damping coefficients as matrices in it. Like this:

I am trying to solve a problem involving 30x30 matrices and if I do it in the above way, I will have to write 60 separate equations for the function F whereas in Matlab, I can pass the previously calculated 30x30 matrices directly into function. Is there any way of doing the same with solve_ivp in python or any such functions?

Thank you.

that doesn't look like the MATLAB code I used for ODEs years ago. That `F(t,y)` could be written as a matrix product of a 2d array and the `y` array. — hpaulj, Oct 04 '21 at 20:38
You can rewrite the (y1, y2, y3, y4) vector in order, grouping the first order derivatives first, followed by the second order derivatives. Then split the problem in two: 1) rewrite the last set of matrix equations by replacing the second derivative with the "ys" `y[0,1] = f(....)`(2D vector) 2) create the equations for the second derivative (`dy[2,3]/dt = y[0,1]`), also 2D vector. Then append both. I'll try to code it if I can find some time. But perhaps you get the idea. — Tarifazo, Oct 04 '21 at 20:48
my point being, the problem is probably that you are alternating the first and second order derivatives in your definition of `y` vector, so you need to write everything by hand. — Tarifazo, Oct 04 '21 at 20:51
I was referring to the last image that proported to be MATLAB. — hpaulj, Oct 04 '21 at 20:58

hpaulj · Accepted Answer · 2021-10-04T23:09:42.023

arr = np.array([
        y[1],
        (1/m1)*(f1*np.cos(3*t) - ((c1 + c1)*y[1] + (k1 + k2)*y[0]) + c2*y[3] + k2*y[2]),
        y[3],
        (1/m2)*(f2*np.sin(t**2) - c2*y[3] - k2*y[2] + c2*y[1] + k2*y[0])
    ])

can be written as (roughly):

f = np.array([0, f1*np.cos(3*t),0,f2*np.sin(t**2)])
M = np.array([
        [0, 1, 0, 0],
        [(k1+K2), (c1+c1), k2, c2],
        [0,0,0,1],
        [k2, c2, ....]])
arr = f[:,None] + M.dot(y)

That M array could be passed via args=(M,) (it's independent of t and y). Or just be a global to the function.

Passing matrices as input in scipy.integrate.solve_ivp python

1 Answers1