i’m working on resolving this differential equation :
z” + w0²z=w0²zeq + a(t)

for this i’m now using odeint or solve_ivp from module scipy.integrate in Python and it seems to work well if i give an explicit formula for a : for example : a(t) = A0 * sin(w*t)

but the fact is that i would like to use the values from experiences i realised. For instance, i used a chronophotography software so as to get a [time,a(time)] array

Does someone know how to use the values of my array ?

my first ideas was : odeint chooses some values for its t variable, knowing this t variable i should be able to get from my [t,a(t)] array the a(t) value that’s nearest the value i need.

So as to test my idea, knowing a matematical solution if a(t)=Asin(wt), i first created a time array, then created the a(t) array linked to my time array: which lead me to this code

import matplotlib.pyplot as plt
import numpy as np
from math import sqrt,sin
from scipy.integrate import solve_ivp,odeint
def get_nearest_ind(l,v,previ):
    """ get the index of the maximum value which is lower to v in l """
    for i in range(previ,len(l)):
        if(l[i]>=v):
            return i-1
    print(l[0],"->",l[-1])
    #raise NameError(f"La valeur {v} n'est pas comprise dans la liste")
    print(f"The value {v} can't be in the list")
    return -1




def a(t,a_tab,t_tab,previ):
    """ return the value of a at the date t"""
    i = get_nearest_ind(t_tab,t,previ[0])
    previ[0]=i
    assert(i<len(a_tab))
    return a_tab[i]





def f(z,t,w0,ta,a_tab,zeq,l):
    #print(a(t,a_tab,ta))
    thea = a(t,a_tab,ta,l)
    ret = [z[1],-(w0**2)*z[0]+(w0**2) * zeq - thea]
    return ret


def create_tab2(t0,tf,nbr_points,k1,k2,lo1,lo2,m,ta_tab,a_tab,h):
    zeq = calcule_zeq(k1,k2,lo1,lo2,m,h)
    w0 = calcule_w0(k1,k2,m)
    t = np.linspace(t0,tf,nbr_points)
    sol=odeint(f,[zeq,0.],t,args=(w0,ta_tab,a_tab,zeq,[0]))
    return t,sol[:,0],sol[:,1]

##my values
g = 9.81 #m/s²
mn = 11.48e-3
mn3 = 3*mn
mb = 304.87e-3
A = 0.0664524209486166
w=15.93063002857565 #rad/s



## cellule 2 : à modifier si on change les conditions

H=30e-2
k1 = 9.81  #kg/s²
lo1= 4e-2  #m
k2 = 40.87 #kg/s²
lo2 = 5.5e-2 #m
m = mn #
nbr_points = 10000
t = np.linspace(0.,5.,nbr_points)




w0 = sqrt(k1/m) #rad/s
zeq = - m*g/k1 + H - lo1

C = A/(w**2-w0**2)
B=-C*w/w0
#print(f"B={B} ; C={C} ; w0 = {w0} ; w = {w} ; zeq={zeq}")
sol1 = B*np.sin(w0 * t) + zeq + C * np.sin(w*t)
solp1 = B*w0*np.cos(w0*t)+C*w*np.sin(w*t)



def fode(z,t,l,j):
    global A,w0,zeq,w
    l.append(-A*sin(w*t))

    return [z[1],-(w0**2)*z[0]+(w0**2) * zeq -A*sin(w*t)]
l1=[]
solode = odeint(fode,[zeq,0.],t,args=(l1,1))



def test():
    t = np.linspace(0.,5.,10000)
    a_tab = A*np.sin(2*np.pi*t/T)
    t,z,vz = create_tab2(0.,5.,10000,k1,0,lo1,0,m,t,a_tab,H)
    plt.plot(t,z,"r",label = "create_tab solution")
    plt.plot(t,sol1,label="mathematical solution")
    plt.plot(t,solode[:,0],"b",label="odeint solution")
    plt.legend()
    plt.show()

this is what i get…

enter image description here