<jupyter_start><jupyter_text>Dormand-Prince Method > Disclaimer: the original code (written in Scilab) and example below are due to <NAME>.<jupyter_code>def dopri(x0, xf, y0, f, rtol, atol):
    r"""
    Dormand-Prince method for solving first-order ODE.
    
    This function implements the Dormand-Prince5(4) method
    with adaptive step through the estimation of the error
    tolerance (with norm sup) from the approximations of
    global error O(h**5) and O(h**4).
    
    It solves the following IVP:
    
    y'=f(x,y), x > x0,
    y(x0) = y0.
    
    The initial step is computed from the algorithm proposed in
    Hairer,E; Nørsett,S.P.; Wanner,G. "Solving Ordinary
    Differential Equations I. Nonstiff Problems" 2nd rev.ed. (1993),
    pg. 169.
    
    The estimatives for the computation of the step follow the ideas
    exposed in this very same reference in pages 167 and 168.
    
    Parameters
    ----------
    x0 : float
        The initial value for the independent variable.
    xf : float
        The final value for the independent variable.
    y0 : array_like
        1-D array of initial values for the dependente variable evaluated at `x_0`.
    f : callable
        The function that represents the first derivative of y with respect to x.
        It must accept two arguments: the point x at which it will be evaluated and
        the value of y at this point.
        
    Returns
    ------
    x : numpy.ndarray
        The point at which the function was evaluated.
    y : numpy.ndarray
        The approximation of the solution.
    hc : numpy.ndarray
        Steps used along the computations.
    rp : list
        Iterations where rejection of the step size occurred.
    """
    
    xi = x0
    yi = y0
    fac = 0.8
    facmax = 5
    facmin = 0.2
    
    a = np.zeros((7, 7))
    a[1, 0] = 1/5
    a[2, 0:2] = np.array([3/40, 9/40])
    a[3, 0:3] = np.array([44/45, -56/15, 32/9])
    a[4, 0:4] = np.array([19372/6561, -25360/2187, 64448/6561, -212/729])
    a[5, 0:5] = np.array([9017/3168, -355/33, 46732/5247, 49/176, -5103/18656])
    a[6, 0:6] = np.array([35/384, 0.0, 500/1113, 125/192, -2187/6784, 11/84])
    
    b = np.array([5179/57600, 0.0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40])
    
    c = np.array([0.0, 1/5, 3/10, 4/5, 8/9, 1.0, 1.0])

    i = 0
    r = 0
    
    x = []
    y = []
    hc = []
    rp = []
    
    done = False

    d0 = np.max(np.abs(y0))
    d1 = np.max(np.abs(f(x0, y0)))
    
    h0 = 1.0e-6 if np.min([d0, d1]) < 1.0e-5 else (0.01 * (d0/d1))
    
    d2 = np.max(np.abs(f(x0 + h0, y0 + h0 * f(x0, y0)) - f(x0, y0)))/h0
    
    h1 = np.max(1.0e-6, h0/1000) if np.max([d1, d2]) < 1.0e-15 else (0.01/np.max([d1, d2]))**(1/5)
    
    h = np.min([100 * h0, h1])
    
    while not done:
        k1 = f(xi, yi)
        k2 = f(xi + h * c[1], yi + h * (a[1, 0] * k1))
        k3 = f(xi + h * c[2], yi + h * (a[2, 0] * k1 + a[2, 1] * k2))
        k4 = f(xi + h * c[3], yi + h * (a[3, 0] * k1 + a[3, 1] * k2 + a[3, 2] * k3))
        k5 = f(xi + h * c[4], yi + h * (a[4, 0] * k1 + a[4, 1] * k2 + a[4, 2] * k3 + a[4, 3] * k4))
        k6 = f(xi + h * c[5], yi + h * (a[5, 0] * k1 + a[5, 1] * k2 + a[5, 2] * k3 + a[5, 3] * k4 + a[5, 4] * k5))
        
        y4 = yi + h * (a[6, 0] * k1 + a[6, 2] * k3 + a[6, 3] * k4 + a[6, 4] * k5 + a[6, 5] * k6)
        
        k7 = f(xi + h * c[6], y4)
        
        y5 = yi + h * (b[0] * k1 + b[2] * k3 + b[3] * k4 + b[4] * k5 + b[5] * k6 + b[6] * k7)
        
        sc = atol + np.max([np.abs(y5), np.abs(yi)]) * rtol
        err = np.max(np.abs(y5 - y4))/sc
        hnew = h * np.min([facmax, np.max([facmin, fac * (1/err)**(1/5)])])
        
        if err <= 1:
            if (xi + h) <= xf:
                i += 1 
                xi += h
                yi = y5
                x.append(xi)
                y.append(yi)
                hc.append(h)
                h = hnew
            else:
                h = xf - xi
                r += 1
                rp.append(i)
        else:
            h = hnew
            r += 1
            rp.append(i)

        done = xi >= xf
        
    return (np.array(x), np.array(y), np.array(hc), rp)<jupyter_output><empty_output><jupyter_text>Example 1<jupyter_code>def F(t, z):
    return np.array([z[1], z[0]**(-13) - z[0]**(-7)])

t, sol, hi, rejp = dopri(0.0, 100.0, np.array([2.0, 0.0]), F, atol=1.0e-8, rtol=1.0e-8)

q = np.insert(sol[:, 0], 0, 2)
p = np.insert(sol[:, 1], 0, 0)

H = (1/2) * (p**2) + (1/12) * (q**(-12) - 2 * q**(-6))
fig, ax = plt.subplots(figsize=(13, 8))
plt.plot(q, p)
plt.scatter(q[0], p[0], marker='o', color='red')
plt.scatter(q[-1], p[-1], marker='x', color='black')
ax.set(title="Phase space", xlabel="$q$", ylabel="$p$");
fig, ax = plt.subplots(figsize=(13, 8))
plt.plot(np.abs(H/H[0] - 1))
ax.set(title="Hamiltonian's relative error along the appproximations", xlabel="t", ylabel="$H$");
fig, ax = plt.subplots(figsize=(13, 8))
plt.plot((q**-13 - q**-7)**2, label=r"$(dV/dt)^{2}$", color='red')
plt.plot((p**2) * (13 * q**-14 - 7 * q**-8), label=r"$d^{2}V/dt^{2}$", color='blue')
plt.legend(loc='upper right', fancybox=True, framealpha=1, shadow=True, borderpad=1, frameon=True)
ax.set(title=r"Graph of $(dV/dt)^{2}$ and $d^{2}V/dt^{2}$", xlabel="t");<jupyter_output><empty_output><jupyter_text>Question 1<jupyter_code>def F(t):
    return 1/2 + sum(5/n * math.sin(n*math.pi/10) * math.cos(n*math.sqrt(10)*t) for n in range(1, 6))


def f(t, u_k):
    return np.array([u_k[1], F(t) * math.cos(u_k[0]) - (1/10) * u_k[1] - 10 * math.sin(u_k[0])])


def mechanical_energy(theta, dtheta):
    return (1/2) * dtheta**2 - 10 * np.cos(theta)

t, sol, hi, rejp = dopri(0.0, 20.0, np.array([0.0, 0.0]), f, atol=1.0e-10, rtol=1.0e-10)
theta = sol[:, 0]
dtheta = sol[:, 1]

energy = mechanical_energy(theta, dtheta)

max_energy_index = np.argmax(energy)
t_max_energy = t[np.argmax(energy)]
max_energy = energy[np.argmax(energy)]


print("Mechanical energy in t = {:.4} is {:.4} J.".format(t_max_energy, max_energy))
fig, ax = plt.subplots(figsize=(13, 8))
plt.plot(t, theta, color='blue', label=r"$\theta$")
plt.plot(t, dtheta, color='red', label=r"$d\theta/dt$")
plt.legend(loc='upper right', fancybox=True, framealpha=1, shadow=True, borderpad=1, frameon=True)
ax.set(title="Solution of the PVI", xlabel="Time", ylabel="Angle and its derivative");
fig, ax = plt.subplots(figsize=(13, 8))
plt.plot(t, energy, color='blue')
plt.scatter(t_max_energy, max_energy, color='red', marker='x', s=100)
plt.text(0.0, max_energy, "Maximum energy in t = {:.2} is {:.2} J.".format(t_max_energy, max_energy))
ax.set(title="Mechanical energy", xlabel="Time", ylabel="Energy");<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19621008654133096, 'fraction_numerical': 0.0783348254252462, 'mean_word_length': 2.38023197175996, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '367707', 'license': 'MIT', 'n_tokens_mistral': 3129, 'n_tokens_neox': 2831, 'n_words': 716, 'path': 'notebooks/math/dormand_prince.ipynb', 'repo_name': 'kmyokoyama/machine-learning', 'size': '6710'}