<jupyter_start><jupyter_text>Self-ConvergenceSolution to Problem Set 4, Problem 2 ~ <NAME>To run: **jupyter notebook SC_4_2.ipynb**  Import Libraries<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams; rcParams["figure.dpi"] = 300
from matplotlib.ticker import (AutoMinorLocator)
#rcParams['text.usetex'] = True
plt.rc('font', family='serif')
plt.rc('xtick', labelsize='x-small')
plt.rc('ytick', labelsize='x-small')<jupyter_output><empty_output><jupyter_text>Functions to define the RK4 Step and the ODE Solver  RK3 (Self-Derived)The third-order Runge-Kutta formula is:$$ \begin{array}{l} k_1 = h f(y_n,t_n)\\ k_2 = h f(y_n+\frac{h}{2}, t_n+\frac{h}{2}k_1)\\ k_3 = h f(y_n+h, t_n-hk_1+2hk_2)\\ y_{n+1} = y_n + \frac{h}{6}(k_1 + 4k_2 + k_3)\label{RK3}\end{array} $$  RK4The classic fourth-order Runge-Kutta formula is:$$ \begin{array}{l} k_1 = h f(y_n,t_n)\\ k_2 = h f(y_n+\frac{k_1}{2}, t_n+\frac{h}{2})\\ k_3 = h f(y_n+\frac{k_2}{2}, t_n+\frac{h}{2})\\ k_4 = h f(y_n+k_3, t_n+h)\\ y_{n+1} = y_n + \frac{k_1}{6}+ \frac{k_2}{3}+ \frac{k_3}{3} + \frac{k_4}{6} \label{RK4}\end{array} $$<jupyter_code>def RK3_step(t, y, h, f):
    """
    Implements a single step of a third-order, explicit Runge-Kutta scheme
    """
    k1 = h * f(t, y)
    k2 = h * f(t + 0.5*h*k1, y + 0.5*h)
    k3 = h * f(t - h*k1 + 2*h*k2, y + h)
    return y + (h/6) * (k1 + 4*k2 + k3)
def RK4_step(t, y, h, g, *P):
    """
    Implements a single step of a fourth-order, explicit Runge-Kutta scheme
    """
    thalf = t + 0.5*h
    k1 = h * g(t, y, *P)
    k2 = h * g(thalf, y + 0.5*k1, *P)
    k3 = h * g(thalf, y + 0.5*k2, *P)
    k4 = h * g(t + h, y + k3, *P)
    return y + (k1 + 2*k2 + 2*k3 + k4)/6
def odeSolve(t0, y0, tmax, h, RHS, method, *P):
    """
    ODE driver with constant step-size, allowing systems of ODE's
    """
    # make array of times and find length of array
    t = np.arange(t0,tmax+h,h)
    ntimes,  = t.shape

    # find out if we are solving a scalar ODE or a system of ODEs, and allocate space accordingly
    
    if type(y0) in [int, float]:  # check if primitive type -- means only one eqn
        neqn = 1
        y = np.zeros( ntimes )
    else:                         # otherwise assume a numpy array -- a system of more than one eqn
        neqn, = y0.shape
        y = np.zeros( (ntimes, neqn) )

    # set first element of solution to initial conditions (possibly a vector) 
    y[0] = y0

    # march on...
    for i in range(0,ntimes-1):
        y[i+1] = method(t[i], y[i], h, RHS, *P)

    return t,y<jupyter_output><empty_output><jupyter_text>Questions from Problem Set * Solve the following initial value problem $$ \frac{dy}{dt} = -2ty, 0 \leq t \leq 3,\;\; y(0)=1 $$  using the RK3 and RK4 implementations developed in class.  Function to describe RHS for the given ODE<jupyter_code>def RHS(t, y):
    """
    Implements the RHS (y'(x)) of the DE
    """
    return -2*t*y
# initial conditions
h = 0.01
t0 = 0.0
y0 = 1.0
tmax = 3.0

# solve the ODE
t, y = odeSolve(t0, y0, tmax, h, RHS, RK4_step)
T, Y = odeSolve(t0, y0, tmax, h, RHS, RK3_step)

# plot
f,a = plt.subplots()
a.plot(t,y,'b', label='RK4')
a.plot(T,Y,'r', label='RK3')
a.set_xlabel('Time')
a.legend()
a.set_title("RK3 v/s RK4 for dy/dt = -2ty", fontweight='bold')
a.xaxis.set_minor_locator(AutoMinorLocator()) 
a.yaxis.set_minor_locator(AutoMinorLocator())
a.tick_params(which='minor', length=2.5, color='k')
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>* Now set $h_2$ to the step size of the solution corresponding to N = 800, and $h_3$ to the step size corresponding to N = 1600. Then consider the step size $h$ corresponding to N = 100,200,300,400,500,600, and on a plot show the curves plot $\frac{y_{num}(h)−y_{num}(h3)}{y_{num}(h2)−y_{num}(h3)}$ vs $h$ and $\frac{(h/h3)^n−1}{2^n−1}$ vs $h$. Use this plot to argue that your code is self-convergent  Function to vary h<jupyter_code>def h(t0,tmax,N):
    return (tmax-t0)/N
# initial conditions
t0 = 0.0
y0 = 1.0
tmax = 3.0

# define N
N2 = 800
N3 = 1600
N = [100, 200, 300, 400, 500, 600]

#define h
h2 = h(t0,tmax,N2)
h3 = h(t0,tmax,N3)
H = np.zeros(len(N))
for i in range(len(N)):
    H[i] = h(t0,tmax,N[i])

# solve the ODE
t2, y2 = odeSolve(t0, y0, tmax, h2, RHS, RK4_step)
t3, y3 = odeSolve(t0, y0, tmax, h3, RHS, RK4_step)
T1, Y1 = odeSolve(t0, y0, tmax, H[0], RHS, RK4_step)
T2, Y2 = odeSolve(t0, y0, tmax, H[1], RHS, RK4_step)
T3, Y3 = odeSolve(t0, y0, tmax, H[2], RHS, RK4_step)
T4, Y4 = odeSolve(t0, y0, tmax, H[3], RHS, RK4_step)
T5, Y5 = odeSolve(t0, y0, tmax, H[4], RHS, RK4_step)
T6, Y6 = odeSolve(t0, y0, tmax, H[5], RHS, RK4_step)

# define LHS and RHS form equation 8 on the problem set
LHS = [(Y1[-1]-y3[-1])/(y2[-1]-y3[-1]),
       (Y2[-1]-y3[-1])/(y2[-1]-y3[-1]),
       (Y3[-1]-y3[-1])/(y2[-1]-y3[-1]),
       (Y4[-1]-y3[-1])/(y2[-1]-y3[-1]),
       (Y5[-1]-y3[-1])/(y2[-1]-y3[-1]),
       (Y6[-1]-y3[-1])/(y2[-1]-y3[-1])]
RHS = []
for i in range(6):
    RHS.append(((H[i]/h3)**N[i] - 1)/(2**N[i] - 1))
    
#plot
f,a = plt.subplots()
a.plot(H,LHS,'r', label=r'$\frac{y_{num}(h)−y_{num}(h3)}{y_{num}(h2)−y_{num}(h3)}$')
a.set_xlabel('$h$')
a.set_title(r"$\frac{y_{num}(h)−y_{num}(h3)}{y_{num}(h2)−y_{num}(h3)}$ for dy/dt = -2ty using RK4", fontweight='bold')
a.xaxis.set_minor_locator(AutoMinorLocator()) 
a.yaxis.set_minor_locator(AutoMinorLocator())
a.tick_params(which='minor', length=2.5, color='k')
a.legend()
plt.tight_layout()
plt.show()

f,a = plt.subplots()
a.plot(H,RHS,'b', label=r'$\frac{(h/h3)^n−1}{2^n−1}$')
a.set_xlabel('$h$')
a.set_title(r"$\frac{(h/h3)^n−1}{2^n−1}$ for dy/dt = -2ty using RK4", fontweight='bold')
a.xaxis.set_minor_locator(AutoMinorLocator()) 
a.yaxis.set_minor_locator(AutoMinorLocator())
a.tick_params(which='minor', length=2.5, color='k')
a.legend()
plt.tight_layout()
plt.show()
# final plot
f,a = plt.subplots()
a.plot(H,LHS,'r', label=r'$\frac{y_{num}(h)−y_{num}(h3)}{y_{num}(h2)−y_{num}(h3)}$')
a.plot(H,RHS,'b--', label=r'$\frac{(h/h3)^n−1}{2^n−1}$')
a.set_xlabel('$h$')
a.set_title("Self-Convergence for dy/dt = -2ty using RK4", fontweight='bold')
a.xaxis.set_minor_locator(AutoMinorLocator()) 
a.yaxis.set_minor_locator(AutoMinorLocator())
a.tick_params(which='minor', length=2.5, color='k')
a.legend()
plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20719016083254493, 'fraction_numerical': 0.053768527278461056, 'mean_word_length': 3.9361867704280153, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '632917', 'license': 'MIT', 'n_tokens_mistral': 3086, 'n_tokens_neox': 2917, 'n_words': 655, 'path': 'Problem Sets/Problem Set 4/SC_4_2.ipynb', 'repo_name': 'astroarshn2000/PHYS305S20', 'size': '6352'}