Solving Systems of ODEs with the FEniCS Library: A Practical Guide

Ordinary differential equations (ODEs) are fundamental in mathematical modelling across various scientific disciplines. Solving these equations analytically can be challenging, especially for complex systems. Fortunately, numerical methods and libraries like FEniCS provide powerful tools for solving such systems. In this article, we will walk through solving a system of ODEs derived from a functional using FEniCS, a rather unconventional method for solving ODEs!

Image generated with DALL-E

The Functional and Derivation of ODEs

Consider a potential functional Π dependent on the solutions (y_1) and (y_2), and their time derivatives (y_1') and (y_2'):

To find the ODEs, we minimize Π with respect to (y_1') and (y_2'):

These equations represent a system of coupled ODEs. Alternatively, they can be written in a more conventional form as:

We will solve these equations numerically using the FEniCS library, starting with the following initial conditions: (y_1(0) = 0) and (y_2(0) = 1).

Finite Element Formulation in FEniCS

For finite element formulation in the FEniCS library, we need to define the weak form of the ODEs. The weak form is obtained by multiplying the ODEs by a test function and integrating over the domain. The test functions are denoted by (v_1) and (v_2).

The weak form of the ODEs is then transformed into the question, find (y_1) and (y_2) such that:

for all test functions (v_1) and (v_2) in the function space (V) with domain (t in interval [0, T]), where T is the final time of the simulation.

Setting Up the FEniCS Environment

First, ensure you have FEniCS installed. You can find installation instructions on the FEniCS Project website.

Next, let’s dive into the code. We’ll set up the mesh, define the function space, and establish the boundary conditions:

from dolfin import *
import matplotlib.pyplot as plt

# Create mesh and function space
ncells = 100
mesh = IntervalMesh(ncells, 0, 2*pi*3)

# Define mixed element
degree_y1 = 2
degree_y2 = 1
Welm = MixedElement([FiniteElement('Lagrange', interval, degree_y1),
                     FiniteElement('Lagrange', interval, degree_y2)])

# Function Space
W = FunctionSpace(mesh, Welm)

# Define boundary conditions
bcsys = [DirichletBC(W.sub(0), Constant(0.0), 'near(x[0], 0)'),
         DirichletBC(W.sub(1), Constant(1.0), 'near(x[0], 0)')]

Defining the Problem

We define the functions for solving the system and the potential functional Π:

# Define functions for Solving the System
up = Function(W)
y1, y2 = split(up)
v1, v2 = split(TestFunction(W))

# Define variables
dy1 = y1.dx(0)
dy2 = y2.dx(0)
y1 = variable(y1)
y2 = variable(y2)
dy1 = variable(dy1)
dy2 = variable(dy2)

# Define Potential Energy, Pi
source = Expression('2.0*sin(x[0])', degree=2)
Pi_func = 0.5*dy1**2 + 0.5*dy2**2 - y2*dy1 - y1*dy2 + source*dy2

Deriving the ODEs and Solving the System

From ΠΠ, we derive the ODEs and set up the weak form:

# Obtain ODEs
F_1 = diff(Pi_func, dy1)
F_2 = diff(Pi_func, dy2)

# Define weak form and Jacobian
weak_form = F_1*v1*dx + F_2*v2*dx
Jac = derivative(weak_form, up, TrialFunction(W))

# Solve the system
solve(weak_form == 0, up, J=Jac, bcs=bcsys)

# Extract solutions
y1, y2 = up.split()

# Get the coordinates for plotting
t_values = mesh.coordinates()

# Evaluate the solutions at the mesh points
y1_values = y1.compute_vertex_values(mesh)
y2_values = y2.compute_vertex_values(mesh)

Visualizing the Results

Finally, we visualize the solutions:

# Plotting
plt.figure(figsize=(10, 6))  
plt.plot(t_values, y1_values, label=r'$y_1(t)$')
plt.plot(t_values, y2_values, label=r'$y_2(t)$')
plt.legend()
plt.title('Solutions $y_1(t)$ and $y_2(t)$')
plt.xlabel('t')
plt.ylabel('Solution values')
plt.show()

The numerical solutions look sinusoidal as shown below!

Conclusion

In this article, we’ve demonstrated how to solve a system of ODEs derived from a functional using the FEniCS library. By defining the potential (\Pi) and minimizing it, we derived the necessary ODEs, which were then solved numerically. This method is powerful for handling complex systems where analytical solutions are not feasible. FEniCS, with its flexible and robust framework, proves to be an invaluable tool for such tasks.

If you found this tutorial helpful, explore more about FEniCS on their official documentation.