import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
= 0.2 # epsilon value
ep = 2 # max time
tmax = 20 # number of steps
N = np.linspace(0, tmax, N) # mesh used for plotting
t = t[1] - t[0]
dt
# Define function for the ODE
def f(t, Y, ep):
= Y
y, yp = -1/(1 + ep*y)**2
ypp return np.array([yp, ypp])
# define the initial condition
= [0.0, 1.0]
Y = 0
ti
# define the solution vector
for i in range(1, N):
= ti + dt # Increment time
ti = Y + f(ti, Y, ep)*dt # Euler step
Y 0], 'k.')
plt.plot(ti, Y[
# Asymptotic solutions
= -1/2*t**2 + t
y0 = -1/12*t**4 + 1/3*t**3
y1 '--')
plt.plot(t, y0, + ep*y1, '--')
plt.plot(t, y0 't');
plt.xlabel('y'); plt.ylabel(
7 Numerical solutions of IVPs
In the previous section, we used built-in ODE solvers to develop numerical solutions. It is important to gain an understanding how a simple ODE solver works. The simplest scheme is called Euler’s method, and this we now explain.
Begin from the system (Equation 6.2). We assume that the solution is represented by a discrete set of points,
Rearranging yields a very simple algorithm for solving the ODE:
This would be implemented via the following pseudocode:
Euler’s method is conceptually simple but quite inaccurate. But in this case, we see that it works fairly well in comparison to the built-in solvers.