Solving Ordinary Differential Equations with SciPy

Welcome to the next lesson in our SciPy course, focusing on solving Ordinary Differential Equations (ODEs). Building on your previous experiences with integration, we will now explore ODEs, which are equations involving a function and its derivatives. ODEs are widely used in fields like physics, engineering, and biology to model dynamic systems, such as the motion of a pendulum or population growth.

Our tool of choice is SciPy, a powerful library in Python that makes solving these equations manageable. By the end of this lesson, you'll be able to solve a simple first-order ODE and visualize the solution using Python's scientific libraries.

Let's start by defining the ODE we aim to solve:

$\frac{dx}{dt} = 2t + 3$

• ODE: An equation involving a function x and its derivative with respect to time t.
• Initial Conditions: Essential for solving ODEs, initial conditions specify the starting value of the function. For our problem, let's use $x(0) = 0$.
• Time Span: The interval over which we want to solve the ODE. For this lesson, we'll consider $t$ from 0 to 10.

Understanding these components is crucial for setting up and solving ODEs effectively.

To solve our ODE, we will use SciPy's solve_ivp function. Let's break this process down step-by-step.

First, we define the ODE as a Python function. This function computes the derivative at any given time, $t$.

def model(t, x):
    return 2 * t + 3

• Function model takes two arguments: time t and the function value x.
• It returns the derivative $\frac{dx}{dt}$, which is $2t + 3$.

We also need to specify the starting value and the time span for the solution.

x0 = [0] 
t_span = (0, 10)

• x0: A list containing the initial condition.
• t_span: A tuple defining the start and end times.

Finally, we'll define the time points where we want to compute the solution.

import numpy as np
t_eval = np.linspace(t_span[0], t_span[1], 100)  # 100 evenly spaced points between 0 and 10

t_eval is an array of time points using NumPy's linspace to evaluate the solution at specific intervals.

With the ODE defined and initial conditions specified, we're ready to solve and visualize the solution.

Let's use solve_ivp to compute the solution.

from scipy.integrate import solve_ivp

solution = solve_ivp(model, t_span, x0, t_eval=t_eval)

• solve_ivp: This function takes our model, time span, initial condition, and evaluation points to compute the solution.
• solution: An object containing the results, including time solution.t and solution values solution.y.

We'll use Matplotlib to plot the solution over time.

import matplotlib.pyplot as plt

plt.plot(solution.t, solution.y[0], label='x(t)')
plt.xlabel('Time t')
plt.ylabel('Solution x')
plt.title('Solution of ODE: dx/dt = 2t + 3')
plt.legend()
plt.show()

• plt.plot: Plots time against the solution.
• plt.xlabel and plt.ylabel: Set the labels for the x and y axes.
• plt.title: Provides a title for the plot.
• plt.legend: Adds a legend for clarity.

Here is the solution:

As you can see, the answer is a quadratic function, as expected.

Let's consider a slightly more complex first-order ODE that models exponential decay with an external forcing function:

$\frac{dx}{dt} = -\alpha x + \sin(t)$

Here's how to implement it in Python using solve_ivp:

def exponential_decay_with_forcing(t, x, alpha):
    return -alpha * x + np.sin(t)

# Parameter for the decay rate
alpha = 0.2

# Initial condition: Starting value of x
x0 = [1.0]
t_span = (0, 20)

# Time evaluation points
t_eval = np.linspace(t_span[0], t_span[1], 300)

# Solving the ODE
solution = solve_ivp(exponential_decay_with_forcing, t_span, x0, args=(alpha,), t_eval=t_eval)

# Plotting the results
plt.plot(solution.t, solution.y[0], label='x(t)')
plt.xlabel('Time t')
plt.ylabel('Solution x')
plt.title(r'Solution of ODE: $\frac{{dx}}{{dt}} = -\alpha x + \sin(t)$')
plt.legend()
plt.show()

• exponential_decay_with_forcing: The function defining our ODE. It includes a decay term -alpha * x and an external forcing function sin(t).
• alpha: Represents the decay rate.
• args=(alpha,): Passes the decay rate as an additional argument to the model function.

Here is the result for this code snippet:

In this lesson, we explored solving an ODE using SciPy's solve_ivp function and visualizing the results with Matplotlib. Starting with defining our ODE, setting initial conditions, solving the equation, and finally plotting the solution, you now have a framework to apply to similar problems.

Practice exercises follow this lesson, which will reinforce your understanding and give you the opportunity to apply these concepts to new situations. Keep exploring, and remember that you've already built a solid foundation for working with scientific computations in Python.