Python Scipy Odeint

In this Python Scipy tutorial, we will learn about the “Python Scipy Odeint” to solve the problem related to differential equations and how to implement it in Python Scipy with the following topics.

• What is the order of the differential equations?
• Python Scipy Odeint
• How different is Odeint from the Solve_ivp method
• How to integrate 2nd order equations
• Python Scipy Odeint Complex
• Python Scipy Odeint Rk45

What is the order of the differential equations?

One or more variables (unknowns) and some of their derivatives make up a differential equation. In other words, the derivatives of the variables are defined by the differential equation.

Differential equations are categorized according to their order. The highest derivative sometimes referred to as the differential coefficient, that is present in a differential equation determines its order.

In this tutorial, we will solve the system of ordinary differential equations using the Python Scipy method.

Python Scipy Odeint

Python Scipy has a method odeint() in a module scipy.integrate that solves a set of ordinary differential equations by integrating them.

The syntax is given below.

scipy.integrate.odeint(func full_output=0, ml=None, mu=None, , y0, t, args=(), Dfun=None, col_deriv=0, rtol=None, atol=None, tcrit=None, h0=0.0, mxhnil=0, mxordn=12, mxords=5, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, printmessg=0, tfirst=False)

Where parameters are:

• func(callable): Calculating y’s derivative at time t. The tfirst argument must be set to True if the signature is callable(t, y,…)
• y0(array): The initial state of y (can be a vector).
• t(array): A list of time intervals for which y needs to be solved. The first component of this sequence should be the beginning value point. Repeated values are permitted; this sequence must be monotonically growing or decreasing.
• args(tuple): Additional parameters to pass to the function.
• dfun: Gradient of function (Jacobian). The parameter tfirst must be set to True if the signature is callable(t, y,…)
• col_deriv(boolean): True if Dfun defines derivatives across rows instead of along columns (which takes longer).
• full_output(boolean): True if the second output should be a dictionary of optional outputs.
• printmessg(boolean): If the convergence message should be printed.
• tfirst(boolean): If this is the case, func’s first two arguments (and Dfun, if provided) must be t, y rather than the normal y, t.

The method odeint() returns y(array with the initial value of y0 in the first row and the value of y for each chosen time in t) of type array.

Let’s take an example by following the below steps:

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

The vector [theta, omega] shall be y. We use Python to implement this system.

def pend_fun(y_, t_, b_, c_):
    theta_, omega_ = y_
    dydt_ = [omega_, -b_*omega_ - c_*np.sin(theta_)]
    return dydt_

We presume that b = 0.20 and c = 4.5 are constants.

b_ = 0.20
c_ = 4.5

The pendulum is initially at rest with omega(0) = 0, and the initial circumstances are that it is approximately vertical with theta(0) = pi – 0.1. The initial conditions vector is then.

y0_ = [np.pi - 0.1, 0.0]

In the range 0 = t = 10, we shall produce a solution at 101 evenly spaced samples. So, here is our range of times.

t_ = np.linspace(0, 10, 101)

To generate the solution, call odeint. We use the args argument to supply the arguments b and c to odeint so that they may be passed to pend.

sol_ =integrate.odeint(pend_fun, y0_, t_, args=(b_, c_))

The answer is a shape-preserving array (101, 2). Theta(t) is the first column, and Omega is the second (t). The code below plots both elements.

plt.plot(t_, sol_[:, 0], 'b', label='theta(t)')
plt.plot(t_, sol_[:, 1], 'g', label='omega(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()

This is how to integrate the differential equation using the method odeint() of Python Scipy.

Read Python Scipy Freqz

Python Scipy Odeint Vs Solve_ivp

The scipy.integrate.ode class and the function scipy.integrate.solve_ivp employ the system definition function, which by default requires the first two parameters of func to be in the opposite order of those arguments. It is necessary to set the argument tfirst to True in order to utilize a function with the signature func(t, y,…).

We have already learned about “Python Scipy Odeint”, in this section we will learn about the method solve_ivp() of Python Scipy.

The method solve_ivp() fix the initial value issue with an ODE system. Given an initial value, this function numerically integrates a set of ordinary differential equations.

dy / dt = f(t, y)
y(t0) = y0

Here, y(t) is an N-D vector-valued function (state), t is a 1-D independent variable (time), and f(t, y) is an N-D vector-valued function that establishes the differential equations. Given an initial value of y(t0)=y0, the aim is to find y(t) roughly satisfying the differential equations.

The right-hand side of stiff ODE solvers must be complex-differentiable, however, some solvers offer integration in the complex domain (satisfy Cauchy-Riemann equations [11]).

Pass y0 with a complex data type in order to solve a complex problem. Rewriting your problem individually for the real and imagined sections is another alternative that is always available.

The syntax is given below.

scipy.integrate.solve_ivp(fun, dense_output=False, events=None,  t_span, y0, method='RK45', t_eval=None,vectorized=False, args=None, **options)

Where parameters are:

• fun(callable): The system is on the right side. Fun with the calling signature (t, y). Here, t is a scalar, and the ndarray y has two possibilities: Alternatively, if it has shape (n,), fun must return an array-like object with shape (n,). It can also have a form (n, k), in which case fun must return an array whose shape (n, k) is such that each column is equivalent to a single column in y. Vectorized argument controls which of the two alternatives is selected. The vectorized implementation enables a quicker Jacobian through finite differences approximation (required for stiff solvers).
• _span(two tuple floats): It represents the interval of integration from t0 to t1.
• y0(array_data like shape(n)): Used for specifying the initial state.
• method(string or odsolver): It is used to specify which integration method to use like RK45, Radau, RK23, DOP853, BDF, and LSODA.
• t_eval(array_data): It is used to specify the times at which to store the calculated solution.
• dense_output(boolean): It is used to specify whether we want to calculate the continuous solution or not.
• events(callable): It is used for tracking events.
• Vectorized(boolean): It is used to specify whether we want to implement the function in a vectorized way or not.
• args(tuple): To pass extra arguments.
• **options: Options given to a selected solver Below is a list of every option that is available for solvers that have previously been implemented.
• first_step(None, float): First-time step size Since None is the default, the algorithm should make the decision.
• max_step(float): Maximum step size permitted. The step size is not bounded and is totally up to the solver when np.inf is used as the default.
• rtol, atol( array_data, float): tolerances, both relative and absolute. The solution maintains local error estimates that are less than atol plus rtol times abs (y). Here, the absolute accuracy is controlled by atol, while the relative accuracy (number of accurate digits) is controlled by rtol (number of correct decimal places). Set atol to be smaller than the least value that can be predicted from rtol * abs(y) so that rtol dominates the permissible error in order to obtain the desired rtol. There is no guarantee that there will be as many accurate numbers if atol is more than rtol * abs(y). On the other hand, to get the required atol, adjust rtol such that it is never more than atol and rtol * abs(y). It may be advantageous to pass array-like with shape (n,) for atol if the scales of the components of y are different from one another.
• jac(array_data, sparse matrix): The right-hand side of the system’s Jacobian matrix with respect to y is what the Radau, BDF, and LSODA methods require. The (n, n)-dimensional Jacobian matrix has an element I j) that is equal to d f_i / d y_j. The Jacobian can be described in three different ways.
  1. The Jacobian is taken to be constant if the matrix or array type is array-like. Not backed by “LSODA.”
  2. The Jacobian will be called as jac(t, y), if necessary, if it is callable, in which case it is presumed to depend on both t and y. The return value for the Radau and BDF algorithms could be a sparse matrix.
  3. If None (the default), finite differences will be used to approximate the Jacobian.
• jac_sparsity(array_data, sparse matrix): outlines the Jacobian matrix’s sparsity structure for a finite-difference approximation. Its form has to be (n, n). If jac is not None, this argument is not considered. The computations will be significantly sped up if the sparsity structure is provided and the Jacobian has only a few non-zero elements in each row [10]. A zero entry indicates that the Jacobian’s associated element is always zero. If None (the default), it is presumed that the Jacobian is dense. Not supported by “LSODA,” instead use lband and uband.
• lband, uband(int): The ‘LSODA’ technique parameters, i.e., jac[i, j]!= 0 only for i – lband <= j <= i + uband, define the Jacobian’s bandwidth. No default is used. In order to set these, your jac method must return the Jacobian in the packed format, which requires that the returned array have n columns and uband + lband + 1 rows with Jacobian diagonals written in them. Jac packed[uband + I – j, j] specifically equals Jac[i, j]. Scipy.linalg.solve banded employs the same format (check for an illustration). In order to lower the number of Jacobian elements estimated via finite differences, these parameters can alternatively be utilised with jac=None.
• min_step(float): The smallest step size permitted by the “LSODA” algorithm. The default value for the min step is 0

The method solve_ivp() returns a bunch of objects.

Let’s take an example and understand how solve_ivp() works by following the below steps:

Import the required libraries or methods using the below python code.

from scipy import integrate

Simple exponential decay with randomly selected time points.

def exponential_decay_fun(t_, y_): return -0.4 * y_
sol = integrate.solve_ivp(exponential_decay_fun, [0, 10], [1, 3, 7])
print(sol.t)

From the above output, we can conclude the difference between odeint and solve_ivp of Python Scipy.

Read Python Scipy Minimize

Python Scipy Odeint 2nd order

A differential equation for a variable of interest frequently relies not just on the first derivative but also on higher derivatives.

We are going to use the method odeint() of Python Scipy that is covered in the above subsection.

Import the required libraries or methods using the below python code.

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt

Define the second-order function and find its derivative using the below code.

def fun(u_,x_):
  return (u_[1],-2*u_[1]-2*u_[0]+np.cos(2*x_))
y0_ =[0,0]
xs_ = np.linspace(1,10,300)
us_ = integrate.odeint(fun,y0_,xs_)
ys_ = us_[:,0]

Plot the above derivative using the below code.

plt.plot(xs_,ys_,'-')
plt.plot(xs_,ys_,'r*')
plt.xlabel('x-values')
plt.ylabel('y-values')
plt.show()

This is how to integrate the 2nd order equation using the method odeint() of Python Scipy.

Read Python Scipy Distance Matrix

Python Scipy Odeint Complex

The Python Scipy has a method complex_ode() in a module scipy.integrate, which is Ode wrapped around complex systems. Similar to ode, but before utilizing integrators, this remaps a complex-valued equation system to a real-valued one.

The syntax is given below.

scipy.integrate.complex_ode(f, jac=None)

Where parameters are:

• f(callable): The equation’s Rhs. Y.shape ==(n,) and t is a scalar. Set f_args by using set f params(*args).
• jac(callable): The rhs’ Jacobian “d f[i] / d y[j] equals jac[i,j]. Calling “set f params(*args)” sets the variable “jac args.” “.

Python Scipy Odeint Rk45

The Runge-Kutta-Fehlberg method, sometimes known as the Fehlberg method, is a numerical analytic approach used to solve ordinary differential equations. Based on the broad class of Runge-Kutta procedures, it was created by the German mathematician Erwin Fehlberg.

Fehlberg’s technique is unique because it is a Runge-Kutta family embedded method, which means it combines many evaluations of the same function to produce methods with different orders and similar error constants.

The approach introduced in Fehlberg’s 1969 publication, known as the RKF45 method, is an order 4 method with an order 5 error estimator. The higher-order embedded technique, which enables the automatic determination of an adaptive stepsize, can be used to estimate and regulate the solution error by conducting an additional calculation.

In this section, we will use the above-described method Rk45 with the method solve_ivp() to solve the differential equations because the method odeint() doesn’t support any kind of parameter method for Rk45.

We will use the same example that we have used in the above subsection by following the below steps:

Import the required libraries or methods using the below python code.

from scipy import integrate

Simple exponential decay with randomly selected time points and with the method RK45 using the below code.

def exponential_decay_fun(t_, y_): return -0.4 * y_
sol = integrate.solve_ivp(exponential_decay_fun, [0, 10], [1, 3, 7], method = 'RK45')
print(sol.t)

This is how to integrate the function using the method Rk45 with the method solve_ivp() of Python Scipy.

We have covered how to find the integration of the equation with higher order and also learned about “what is Solve_ivp method” etc, with the following topics.

• What is the order of the differential equations?
• Python Scipy Odeint
• How different is Odeint from the Solve_ivp method
• How to integrate 2nd order equations
• Python Scipy Odeint Complex
• Python Scipy Odeint Rk45

You may like the following Python Scipy tutorials:

• Python Scipy Confidence Interval
• Python Scipy Exponential
• Python Scipy Leastsq
• Python Scipy Chi-Square Test
• Scipy Find Peaks
• Python Scipy FFT