Solve Differential Equations with Time-Varying Inputs and Coefficients in Python – Fusion of Engineering, Control, Coding, Machine Learning, and Science

In this Python tutorial, we will learn how to solve differential equations with time-varying inputs and coefficients in Python. We use the odeint() Python function. The GitHub page with all the codes is given here.

The YouTube tutorial accompanying this post is given below.

Test case – Nonlinear Pendulum

To explain the procedure for solving ordinary differential equations with time-varying inputs and coefficients, we need a test case. As a test case, we are using an example of a pendulum system shown in the Figure below.

Figure 1: Pendulum system used as a test case for solving ordinary differential equations in Python.

A ball (red color in the figure) with a mass of $m$ is attached by using a massless rod to the pivot point $Q$ . We assume that the force $F$ is acting at the ball. The length of the rod is $l$ . The force is always perpendicular to the rod. In the figure above, $g$ is the gravitational acceleration constant. We assume that the mass of the rod is significantly smaller than the mass of the ball and consequently, we can neglect it. The free-body diagram is shown in the figure below.

Figure 2: Free body diagram of the system.

In the figure above, $N$ is the normal reaction force exerted by the rod on the ball, $mg$ is the gravitational force, $F$ is the control force, $n$ is the normal unit vector in the direction of the rod, and $\tau$ is the tangent unit vector perpendicular to the rod (tangent to the circle describing the trajectory of the ball). Note that $n$ and $\tau$ define two perpendicular axes. From Newton’s law, we obtain

(1) $\begin{align*}m\vec{a}=\vec{N}+\vec{F}+m\vec{g}\end{align*}$

where $\vec{a}$ is the acceleration vector. By scalarly multiplying this equation with $\vec{\tau}$ unit vector, we obtain the projection of this equation onto the $\tau$ axis (tangent axis):

(2) $\begin{align*}ma_{\tau} = F(t)- mg\sin(\theta)\end{align*}$

where we assume that the force $F(t)$ is a function of time $t$ . That is, in this tutorial, we explain how to numerically solve this equation when $F$ is an arbitrary function of time. For generality, we assume that this force function is only defined at certain discrete-time points and values of this function will be interpolated at other time points. The approach presented in this tutorial can easily be generalized to the case when the analytical form of the function $F(t)$ is known.

On the other hand, the tangential acceleration is given by

(3) $\begin{align*}a_{\tau}=l\ddot{\theta}\end{align*}$

where $\ddot{\theta}$ is the second time derivative of $\theta$ . The variable $\ddot{\theta}$ is called the angular acceleration. By substituting this expression in the previous equation, we obtain

(4) $\begin{align*}ml\ddot{\theta} =F(t)- mg\sin(\theta)\end{align*}$

From this equation, we obtain

(5) $\begin{align*}\ddot{\theta}+\frac{g}{l}\sin(\theta) =\frac{1}{ml}F(t)\end{align*}$

Let us write the ordinary differential equation (5) in the state-space form. This is necessary since Python’s solver needs a state-space form of differential equation in order to solve it.

First, we introduce the state-space variables

(6) $\begin{align*}x_{1}=\theta \\x_{2}=\dot{\theta}\end{align*}$

By differentiating the last two equations, we obtain

(7) $\begin{align*}\dot{x}_{1}& =\dot{\theta}=x_{2} \\\dot{x}_{2}& =\ddot{\theta}=-\frac{g}{l}\sin(\theta)+\frac{1}{ml}F(t)=-\frac{g}{l}\sin(x_{1})+\frac{1}{ml}F(t)\end{align*}$

Consequently, the state-space model has the following form

(8) $\begin{align*}\begin{bmatrix} \dot{x}_{1} \\ \dot{x}_{2} \end{bmatrix} =\begin{bmatrix} x_{2} \\ -\frac{g}{l}\sin(x_{1})+\frac{1}{ml}F(t) \end{bmatrix}\end{align*}$

Python Numerical Solution

We compute the solution of the system (8), by using the odeint() function and interpolation. We are using the interpolation approach since we are assuming that the force F(t) is given at the discrete-time points. We can easily modify approach when the analytical form of the force $F(t)$ is given. The GitHub page with all the codes is given here.

First, we need to import the necessary libraries and define the constants/parameters

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

# define the constants
gConstant=9.81 
lConstant=1
mConstant=5

Next, let us test interpolation in Python. The code presented in the sequel will interpolate values of a sinusoidal function at desired points that are not defined by a given set of discrete time instants and corresponding values of the sinusoidal function.

###############################################################################
# first, let us test the interpolation

timeArray=np.linspace(0,2*np.pi,100)
sinArray=np.sin(timeArray)
plt.plot(timeArray, sinArray)

# interpolate the values
np.interp((timeArray[9]+timeArray[10])/2,timeArray, sinArray)

# compare the interpolated values with sinArray[9] and sinArray[10]
# the interpolated value should be between these two values
###############################################################################

Now that we know how to interpolate the function values, we can proceed further. For illustration, we generate data points corresponding to an input force function that has the following analytical form

(9) $\begin{align*}F(t)=sin(t)+cos(2t)\end{align*}$

This is achieved by using the following code

# now simulate the dynamics

startTime=0
endTime=5
timeSteps=1000

# simulation time array 
# we will obtain the solution at the time points defined by 
# the vector simulationTime
simulationTime=np.linspace(startTime,endTime,timeSteps)

# define the force input 
forceInput = np.sin(simulationTime)+np.cos(2*simulationTime)
plt.plot(simulationTime, forceInput)
plt.xlabel('time')
plt.ylabel('Force - [N]')
plt.savefig('inputSequence.png',dpi=600)
plt.show()

The force function is shown in the figure below.

Figure 3: Time-varying force applied to the system.

To simulate the differential equation, we need to define a function that evaluates the right-hand side of the model (8) for the current simulation time $t$ , the value of the force $F(t)$ , and the state $x(t)$ .

We do that by using the following Python code.



# this function defines the state-space model, that is, its right-hand side
# x is the state 
# t is the current time - internal to solver
# timePoints - time points vector necessary for interpolation
# g,l,m constants provided by the user
# forceArray - time-varying input
def stateSpaceModel(x,t,timePoints,g,l,m,forceArray):
    # interpolate input force values
    # depending on the current time
    forceApplied=np.interp(t,timePoints, forceArray)
    # NOTE THAT IF YOU KNOW THE ANALYTICAL FORM OF THE INPUT FUNCTION 
    # YOU CAN JUST WRITE THIS ANALYTICAL FORM AS A FUNCTION OF TIME 
    # for example in our case, we can also write
    # forceApplied=np.sin(t)+np.cos(2*t)
    # and you do not need to specity forceArray as an input to the function
    # HOWEVER, IF YOU DO NOT KNOW THE ANALYTICAL FORM YOU HAVE TO USE OUR APPROACH 
    # AND INTERPOLATE VALUES
    
    # right-side of the state equation
    dxdt=[x[1],-(g/l)*np.sin(x[0])+(1/(m*l))*forceApplied]
    return dxdt

This function accepts the following input parameters:

“x” is a 2 by 1 vector (or a list) – this is the state vector at the time instant $t$
“t” is the current time that is varied by the solver
“timePoints” is a vector or an array of time points that is used for interpolating the values of the input force
“g, l, m” are parameters of the model
“forceArray” is the vector of forces at the points “timePoints”

This function will interpolate the values of the force function for the current value of time “t”. Then, for such a value, and for the value of the current state it will evaluate the right-hand side of the state equation. Here it should be noted that we can easily modify this function to accept the analytic form of the function $F$ . That is, we can simply define the force as “forceApplied=np.sin(t)+np.cos(2*t)”. However, this analytical approach is only feasible if the analytical form of the force is known a priori.

We solve the problem by using the following code lines.

# define the initial state for simulation 
# position and velocity
initialState=np.array([0,0])
# generate the state-space trajectory
solutionState=odeint(stateSpaceModel,initialState,simulationTime,args=(simulationTime,gConstant,lConstant,mConstant,forceInput))


plt.plot(simulationTime, solutionState[:,0],'b',label='x1')
plt.plot(simulationTime, solutionState[:,1],'r',label='x2')
plt.xlabel('time')
plt.ylabel('States')
plt.legend()
plt.savefig('simulationResult.png',dpi=600)
plt.show()

First, we define the initial state. The initial state is defined by the vector “initialState”. We numerically solve the differential equation by using the function “odeint()”. The input arguments of this function are:

“stateSpaceModel” – previously defined function that defines the right-hand side of the state-space model (8)
“initialState” – previously defined initial state for solving the problem
“simulationTime” – array of time steps at which the solution is computed
“(simulationTime,gConstant,lConstant,mConstant,forceInput)” – this is a tuple that is passed as a set of parameters to the function. It consists of “simulationTime”, model parameters, and the array “forceInput” that together with “simulationTime” define the set of points that are used for interpolating the input force.

This code produces the graph shown below which presents the simulated state trajectories.

Test case – Nonlinear Pendulum

Python Numerical Solution

admin

You might also like

Control Systems Lecture: From Differential Equations to Transfer Functions and Response Simulation in MATLAB

Explanation and Python Implementation of On-Policy SARSA Temporal Difference Learning – Reinforcement Learning Tutorial with OpenAI Gym

Stress on an Oblique Plane Under the Action of Axial Force