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
where x = x(t) , the unknown function, and F=F(t, u, v) have N components, denoted by xiand Fi, i = 1,2,…,N, respectively.
Let’s consider the following problem:
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 ΔUi/Δti:
Value of ΔUi is equal to difference between two consecutive time values. Then we get:
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:
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 U2 = U U3 = U U4 = U U5 = U U6 = U U7 = U U8 = U # The same for array Ut. Ut1 = Ut Ut2 = Ut Ut3 = Ut Ut4 = Ut Ut5 = Ut Ut6 = Ut Ut7 = Ut Ut8 = Ut # 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 = 1.0 U = 1.0 U = 1.0 U = 1.0 U = 1.0 U = 1.0 U = 1.0 U = 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 = 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.plot(time_solution, u_solution[:,0]) ax.set(ylabel='U1') ax.grid() ax.plot(time_solution, u_solution[:,1]) ax.set(xlabel='time (s)', ylabel='U2') ax.grid() #Plot the solution of differential-algebraic equations plt.show()
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:
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 to get professional help with your simulations and mathematical modelling problems.