CODEWITHSUNDEEP
pythonscipynumerical-integrationodeint
TMoover
TMoover

Reputation: 662

scipy.integrate.odeint update initial conditions every timestep

I´m using scipy.integrate.odeint to solve the equations of motion of a given system with a script from where I selected the most relevant part to this specific problem:

# Equations of Motion function to integrate 
def solveEquationsofMotion(y0, t, nRigidBodies, nCoordinates, nConstraintsByType, dataConst, Phi, dPhidq, niu, gamma, massMatrix, gVector, alpha, beta, sda_Parameters):
    ...
    Some calculations 

    matA = numpy.array
    ...

    dydt = np.hstack((qp,qpp))
    return dydt

#Integrator results 
solution = odeint(solveEquationsofMotion, y0, time_span,args=(nRigidBodies, nCoordinates, nConstraintsByType, dataConst, Phi, dPhidq, niu, gamma, massMatrix, gVector, alpha, beta), full_output=0)

and it works fine.

However now I need to multiply part of the integration result (solution variable) by matA variable in each timestep to use again as the initial conditions for the next timestep.

I've looked in the scipy.integrate.odeint documentation but I haven't seen any relevant information. Any help would be very much appreciated.

Kind Regards

Ivo

Upvotes: 1

Views: 1242

Answers (1)

user6655984
user6655984

Reputation:

If you have to change the solution at every step, it is more logical to use the step-by-step integrator ode. It is supposed to be used in a loop anyway, so one may as well change the conditions meanwhile. Here is an example of solving y' = -sqrt(t)*y (vector valued) where y is multiplied by matA after every step.

The steps in t variable are determined by the array t. The main step is y[k, :] = r.integrate(t[k]) which gets the next value of solution, and then the initial condition is changed by r.set_initial_value(matA.dot(y[k, :]), t[k]). 

import numpy as np
from scipy.integrate import ode
def f(t, y):
    return -np.sqrt(t)*y 

matA = np.array([[0, 1], [-1, 0]])
t = np.linspace(0, 10, 20)
y_0 = [2, 3]
y = np.zeros((len(t), len(y_0)))
y[0, :] = y_0
r = ode(f)
r.set_initial_value(y[0], t[0])

for k in range(1, len(t)):
    y[k, :] = r.integrate(t[k])
    r.set_initial_value(matA.dot(y[k, :]), t[k])

The values of y thus obtained are neither monotone nor positive, as the actual solution of the ODE would be - this shows that the multiplication by matA had an effect. 

[[ 2.00000000e+00  3.00000000e+00]
 [ 1.55052494e+00  2.32578740e+00]
 [ 1.46027833e+00 -9.73518889e-01]
 [-5.32831945e-01 -7.99247918e-01]
 [-3.91483887e-01  2.60989258e-01]
 [ 1.16154133e-01  1.74231200e-01]
 [ 7.11807536e-02 -4.74538357e-02]
 [-1.79307961e-02 -2.68961942e-02]
 [-9.45453427e-03  6.30302285e-03]
 [ 2.07088441e-03  3.10632661e-03]
 [ 9.57623940e-04 -6.38415960e-04]
 [-1.85274552e-04 -2.77911827e-04]
 [-7.61389508e-05  5.07593005e-05]
 [ 1.31604315e-05  1.97406472e-05]
 [ 4.85413044e-06 -3.23608696e-06]
 [-7.56142819e-07 -1.13421423e-06]
 [-2.52269779e-07  1.68179853e-07]
 [ 3.56625306e-08  5.34937959e-08]
 [ 1.08295735e-08 -7.21971567e-09]
 [-1.39690370e-09 -2.09535555e-09]]

Upvotes: 4

Related Questions