Imitate ode45() Function in Python

Imitate ode45() Function in Python

Ordinary differential equations are used in MatLab to solve many scientific problems. The ode45() is used in MatLab to solve differential equations.

This article will show how we can imitate the ode45() function in Python.

Imitate the ode45() Function in Python

To imitate the ode45() function in python, we can use the solve_ivp() method defined in the scipy module. The solve_ivp() method integrates a system of ordinary differential equations (ODEs).

  • The solve_ivp() method takes a function as its first input argument. The function given in the input argument must return an array containing the coefficients of the differential equation.
  • In the second input argument, the solve_ivp() method takes a tuple or list containing two numerical values. The values represent the interval of integration, where the first value in the tuple represents the start of the interval and the second value of the tuple represents the highest value in the interval.
  • In the third input argument, the solve_ivp() method takes an array representing the initial values.
  • After execution, the solve_ivp() method returns a bunch object with various attributes.
    1. The t attribute contains a numpy array containing time points.
    2. The y attribute contains a numpy array with values and time points in t.
    3. The sol attribute contains an Odesolution object containing the solution of the differential equation. If the dense_output argument is set to false in the solve_ivp() method, the sol attribute contains None.

To understand this better, see the following example.

from scipy.integrate import solve_ivp

def function(t, y):
    return 2 * y

interval = [0, 10]
initial_values = [10, 15, 25]
solution = solve_ivp(function, interval, initial_values)
print("Time:", solution.t)
print("Y:", solution.y)

Output:

Time: [ 0.          0.07578687  0.56581063  1.18741382  1.85887096  2.55035821
  3.25007544  3.95320486  4.65775424  5.36289544  6.06828346  6.77377445
  7.47930839  8.18486026  8.89041961  9.59598208 10.        ]
Y: [[1.00000000e+01 1.16366412e+01 3.10073783e+01 1.07492109e+02
  4.11689241e+02 1.64114780e+03 6.65071446e+03 2.71362627e+04
  1.11036049e+05 4.54874443e+05 1.86437495e+06 7.64300835e+06
  3.13352156e+07 1.28474398e+08 5.26752964e+08 2.15973314e+09
  4.84541488e+09]
 [1.50000000e+01 1.74549617e+01 4.65110674e+01 1.61238163e+02
  6.17533861e+02 2.46172171e+03 9.97607169e+03 4.07043941e+04
  1.66554074e+05 6.82311665e+05 2.79656243e+06 1.14645125e+07
  4.70028233e+07 1.92711598e+08 7.90129446e+08 3.23959970e+09
  7.26812231e+09]
 [2.50000000e+01 2.90916029e+01 7.75184457e+01 2.68730272e+02
  1.02922310e+03 4.10286951e+03 1.66267862e+04 6.78406569e+04
  2.77590123e+05 1.13718611e+06 4.66093739e+06 1.91075209e+07
  7.83380389e+07 3.21185996e+08 1.31688241e+09 5.39933284e+09
  1.21135372e+10]]

In the above example, we first defined a function named function that takes t and y as its input argument and returns a value based on y.

Then, we have defined an interval and initial values for the ODE using the variables interval and initial_values, respectively. We pass function, interval, and initial_values as input arguments to the solve_ivp() function, and lastly, we get the output in the variable solution.

In the output, you can observe that the time values are spread throughout the interval 0 to 10. Similarly, the output contains a y value corresponding to each time value.

We can also explicitly specify the time points in the attribute t of the solution. For this, we need to pass an array containing the desired time values for which we need the y values to the t_eval argument of the solve_ivp() method, as shown below.

from scipy.integrate import solve_ivp

def function(t, y):
    return 2 * y

interval = [0, 10]
initial_values = [10, 15, 25]
time_values = [1, 2, 3, 6, 7, 8]
solution = solve_ivp(function, interval, initial_values,t_eval=time_values)
print("Time:", solution.t)
print("Y:", solution.y)

Output:

Time: [1 2 3 6 7 8]
Y: [[7.38683416e+01 5.46053271e+02 4.03089733e+03 1.62618365e+06
  1.20160156e+07 8.87210156e+07]
 [1.10802512e+02 8.19079906e+02 6.04634600e+03 2.43927547e+06
  1.80240234e+07 1.33081523e+08]
 [1.84670854e+02 1.36513318e+03 1.00772433e+04 4.06545912e+06
  3.00400390e+07 2.21802539e+08]]

You can see that the time values only contain those values that are passed as input arguments to the t_eval parameter. Similarly, the attribute y contains values for only the specified t values.

This approach can help you obtain values for certain points in the interval.

Related Article - Python Function

  • Exit a Function in Python
  • Optional Arguments in Python
  • Fit a Step Function in Python
  • Built-In Identity Function in Python
  • Arguments in the main() Function in Python
  • Python Functools Partial Function