Skip to content
Home » Calculations in Chemical Engineering Blog » How to solve a system of differential-algebraic equations in Python

How to solve a system of differential-algebraic equations in Python

A differential-algebraic equation (DAE) is an equation involving an unknown function and its derivatives. A (first order) DAE in its most general form is given by:

How to a solve system of differential-algebraic equations in Python

DAE definition

where x(t) , the unknown function, and F=F(t, u, v) have N components, denoted by xiand Fi, i = 1,2,…,N, respectively.

Example

Let’s consider the following problem:

System of differential-algebraic equations

There are several libraries which can be used for solving such problems, but there are sometimes problems with installation. This time, I will show how to solve such problem using fsolve from SciPy.

The method involves discretization of time derivatives. Instead of dUi/dt we use ΔUiti:

Transformation of DAE

Value of ΔUi is equal to difference between two consecutive time values. Then we get:

Discretization of DAE - finite differences

This way, the system of differential-algebraic equations was transformed to a system of algebraic equations. Now, we need to move non-zero elements to left-hand-side of the equation, as follows:

Final form of DAE

Finally, at every time step we need to solve system of algebraic equations f0, f1,…f7.

Solution of a system of differential-algebraic equations

The necessary imports:

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

First, we will define function describing the algebraic equations:

def dae_eqs(U, *args):

    # U denotes vector of variables at time t + time_step. U(0) means U1, U(1) is U2 and so on.

    # *args contains Ut and time_step
        # Ut denotes vector of variables at time t.
        # value of time step is necessary for calculating discretized derivatives.

    # Unpack arguments
    Ut, time_step = args
    # set some arbitray values of parameters
    c1 = 1.0
    c2 = 2.0
    c3 = 3.0
    c4 = 4.0
    R2 = 2.0
    R3 = 3.0
    R4 = 4.0
    B1 = 1.0
    B2 = 2.0
    E = 5.0

    # Let's assign values of U array to U1, U2, ... for ease.
    U1 = U[0]
    U2 = U[1]
    U3 = U[2]
    U4 = U[3]
    U5 = U[4]
    U6 = U[5]
    U7 = U[6]
    U8 = U[7]
    # The same for array Ut.
    Ut1 = Ut[0]
    Ut2 = Ut[1]
    Ut3 = Ut[2]
    Ut4 = Ut[3]
    Ut5 = Ut[4]
    Ut6 = Ut[5]
    Ut7 = Ut[6]
    Ut8 = Ut[7]

    # Derivatives of U valuses. Some of them are not used, but let's leave them for consistency.
    U1_deriv = (U1-Ut1)/time_step
    U2_deriv = (U2-Ut2)/time_step
    U3_deriv = (U3-Ut3)/time_step
    U4_deriv = (U4-Ut4)/time_step 
    U5_deriv = (U5-Ut5)/time_step
    U6_deriv = (U6-Ut6)/time_step
    U7_deriv = (U7-Ut7)/time_step   
    U8_deriv = (U8-Ut8)/time_step   

    # calculate values of algebraic equations
    f0 = c1 * U3_deriv + c3 * U7_deriv - E
    f1 = c1 * U3_deriv - U4/R3
    f2 = c2 * U5_deriv + c4 * U8_deriv - E
    f3 = c2 * U5_deriv - U6/R4
    f4 = B1 + U7 - U4 - U3
    f5 = U6 + U5 - U8 - B2
    f6 = U7 + U8 - E*R2 - U1
    f7 = B1 - U8 - E*U2 - U1 - U3 - U4

    # Gather all fs into one array
    f = np.array([f0, f1, f2, f3, f4, f5, f6, f7])
    return f

According to the explanation in the first part of the article, the main program looks as follows:

# set initial conditions for t = 0
number_of_unknowns = 8
U = np.zeros(number_of_unknowns)

# If initial conditions are different from zero, they should be set separately, as below
U[0] = 1.0
U[1] = 1.0
U[2] = 1.0
U[3] = 1.0
U[4] = 1.0
U[5] = 1.0
U[6] = 1.0    
U[7] = 1.0   

# set time step, end_time and time
time_step = 0.1
end_time = 5
time = 0

# Initial condition must satisfy the system. Therefore it is required to find proper solution at time 0.
u_0 = fsolve(dae_eqs, U, args=(U,time_step))

# Create arrays for storing solution
time_solution = np.linspace(0, end_time, (int)(end_time/time_step)+1)
u_solution = np.zeros((time_solution.size, number_of_unknowns))
i = 0
u_solution[0] = u_0 # assign initial condition
while time < end_time - time_step:
    args = (U, time_step) # arguments for fsolve
    u_next =  fsolve(dae_eqs, U, args=args) # calculate u at time t + time_step
    U = u_next # in next iteration U is equal U from previous iteration 
    u_solution[i+1] = u_next
    time += time_step # increase time
    i += 1

We can plot some results now:

fig, ax = plt.subplots(2)
ax[0].plot(time_solution, u_solution[:,0])

ax[0].set(ylabel='U1')
ax[0].grid()

ax[1].plot(time_solution, u_solution[:,1])

ax[1].set(xlabel='time (s)', ylabel='U2')
ax[1].grid()

#Plot the solution of differential-algebraic equations
plt.show()
Result of DAE

Summary

Above, we have seen how to solve differential-algebraic equations in Python. Such problems appear in many fields of engineering, including biochemical, electrical and process engineering.

The above method is not perfect. It uses the simplest discretization and it will not always work, for example, if one deals with stiff differential equations. There are several Python libraries which can be used. Some of them are listed below:

https://pypi.org/project/diffeqpy/

https://lcvmwww.epfl.ch/software/daepy/

You can lear more about DAEs on Scholarpedia.

If you want to learn more about solving engineering problems using Python, see the offer for the Course: Python in engineering and science

Contact Us

Contact us to get professional help with your simulations and mathematical modelling problems.


Read related content

PID controller in Python

Proportional-Integral-Differential controllers are commonly used in the industry. Here, you will find an implementation of a control system with a PID control.

Lecture notes in process control

A set of lectures and exercises in process control. I prepared it during several years when I teached foreign students in CUT.

Python for engineers

A list of Python libraries every process, operations engineers should know.