<jupyter_start><jupyter_text>Grade 5/7No implementation of Cash-Karp coefficientsImplementation of adaptive step is not set by Cash Karp algorithm  Create a notebook to perform Runge-Kutta integration for multiple coupled variables.<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>This cell and the one following it are not a requirement, it is only for looks<jupyter_code>#use colors.subclass(or command; e.g bold).colorname to print
#examples: print(colors.bold, colors.fg.blue, "this will be bold and blue")
#everything after this will have that format until the following command
#is given: print(colors.reset, "now, this text will be normal")

class colors:
    reset='\033[0m' #reset all colors with colors.reset
    bold='\033[01m'
    underline='\033[04m'
    strikethrough='\033[09m'
    reverse='\033[07m'
    class fg:  #foreground subclass
        black='\033[30m'
        red='\033[31m'
        green='\033[32m'
        orange='\033[33m'
        blue='\033[34m'
        purple='\033[35m'
        cyan='\033[36m'
        lightgrey='\033[37m'
        darkgrey='\033[90m'
        lightred='\033[91m'
        lightgreen='\033[92m'
        yellow='\033[93m'
        lightblue='\033[94m'
        pink='\033[95m'
        lightcyan='\033[96m'
    class bg: #background subclass
        black='\033[40m'
        red='\033[41m'
        green='\033[42m'
        orange='\033[43m'
        blue='\033[44m'
        purple='\033[45m'
        cyan='\033[46m'
        lightgrey='\033[47m'<jupyter_output><empty_output><jupyter_text>The above code was provided by https://www.geeksforgeeks.org/print-colors-python-terminal/  Define our coupled derivatives to integrate<jupyter_code>def dydx(x,y):
    
    #set the derivatives
    
    #so we can write 
    #dydx = z
    #dzdx = -y
    
    #we will set y = y[0]
    #we will set z = y[1]
    
    #declare an array
    y_derivs = np.zeros(2)
    
    #set dydx = z
    y_derivs[0] =    y[1]
    
    #set dzdx = -y
    y_derivs[1] = -1*y[0]
    
    #here we have to return an array with dydx
    return y_derivs<jupyter_output><empty_output><jupyter_text>Define the 4th order RK method<jupyter_code>def rk4_mv_core(dydx,xi,yi,nv,h):
    
    #declare k? arrays; (? is a wildcard, used for k1,k2,...,kn)
    k1=np.zeros(nv)
    k2=np.zeros(nv)
    k3=np.zeros(nv)
    k4=np.zeros(nv)
    
    #define x at 1/2 step
    x_ipoh = xi + 0.5*h
    
    #define x at 1 step
    x_ipo = xi + h
    
    #declare a temp y array
    y_temp = np.zeros(nv)
    
    #get k1 values
    y_derivs = dydx(xi,yi)
    k1[:] = h*y_derivs[:]
    
    #get k2 values
    y_temp[:] = yi[:] + 0.5*k1[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k2[:] = h*y_derivs[:]
    
    #get k3 values
    y_temp[:] = yi[:] + 0.5*k2[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k3[:] = h*y_derivs[:]
    
    #get k4 values
    y_temp[:] = yi[:] + k3[:]
    y_derivs = dydx(x_ipo,y_temp)
    k4[:] = h*y_derivs[:]
    
    
    #advance y by a step h
    yipo = yi + (k1 + 2*k2 + 2*k3 + k4)/6.
    
    return yipo<jupyter_output><empty_output><jupyter_text>Define an adaptive step size driver for RK4<jupyter_code>def rk4_mv_ad(dydx,x_i,y_i,nv,h,tol):
    
    #define safety scale
    SAFETY    = 0.9
    H_NEW_FAC = 2.0
    
    #set a max number of iterations
    imax = 10000
    
    #set an iteration variable
    i = 0
    
    #create an error
    Delta = np.full(nv,2*tol)
    
    #remember the step
    h_step = h
    
    #adjust step
    while(Delta.max()/tol > 1.0):
        
        #estimate our error by taking one step of size h
        #vs. two steps of size h/2
        y_2  = rk4_mv_core(dydx,x_i,y_i,nv,h_step)
        y_1  = rk4_mv_core(dydx,x_i,y_i,nv,0.5*h_step)
        y_11 = rk4_mv_core(dydx,x_i+0.5*h_step,y_1,nv,0.5*h_step)
        
        #compute an error
        Delta = np.fabs(y_2 - y_11)
        
        #if the error is too large, take a smaller step
        if(Delta.max()/tol > 1.0):
            
            #our error is too large, decrease the step
            h_step *= SAFETY * (Delta.max()/tol)**(-0.25)
            
        #check iteration
        if(i>=imax):
            print("Too many iterations in rk4_mv_ad()")
            raise StopIteration("Ending after i = ",i)
            
        #iterate
        i+=1
        
        
    #next time, try to take a bigger step
    h_new = np.fmin(h_step * (Delta.max()/tol)**(-0.9), h_step*H_NEW_FAC)
    
    #return the answer, a new step, and the step we actually took
    return y_2, h_new, h_step<jupyter_output><empty_output><jupyter_text>Define a wrapper for RK4<jupyter_code>def rk4_mv(dydx,a,b,y_a,tol):
    
    #dfdx is   derivative w.r.t. x
    #a    is   lower bound
    #b    is   upper bound
    #y_a  are  boundary conditions
    #tol  is   tolerance for integrating y
    
    #define our starting step
    xi = a
    yi = y_a.copy()
    
    #an initial step size == make very small
    h = 1.0e-4 * (b-a)
    
    #set max number of iterations
    imax = 10000
    
    #set an iteration variable
    i = 0
    
    #set the number of coupled ODEs to the size of y_a
    nv = len(y_a)
    
    
    #set the initial conditions
    x = np.full(1,a)
    y = np.full((1,nv),y_a)
    
    #set a flag
    flag = 1
    
    #loop until we reach the right side
    while(flag):
        
        #calculate y_i+1
        yi_new, h_new, h_step = rk4_mv_ad(dydx,xi,yi,nv,h,tol)
        
        #update the step
        h = h_new
        
        #prevent an overshoot
        if(xi+h_step>b):
            
            #take a smaller step
            h = b-xi
            
            #recalculate y_i+1
            yi_new, h_new, h_step = rk4_mv_ad(dydx,xi,yi,nv,h,tol)
            
            #break
            flag = 0
            
            
        #update values
        xi += h_step
        yi[:] = yi_new[:]
        
        #add the step to the arrays
        x = np.append(x,xi)
        y_new = np.zeros((len(x),nv))
        y_new[0:len(x)-1,:] = y
        y_new[-1,:] = yi[:]
        del y
        y = y_new
        
        #prevent too many iterations
        if(i>=imax):
            
            print(colors.bold, colors.fg.red, colors.bg.black, "Maximum iterations reached.", colors.reset)
            raise StopIteration("Iteration number = ",i)
            
            
        #iterate
        i += 1
        
        #output some information
        s = "i =%3d\tx = %9.8f\th = %9.8f\tb=%9.8f" % (i,xi,h_step,b)
        print(s)
        
        #break if new xi is == b
        if(xi==b):
            flag = 0
            
    #return the answer
    print(colors.bold, colors.fg.purple, "Number of iterations:", i, colors.reset)
    return x,y<jupyter_output><empty_output><jupyter_text>Perform the integration<jupyter_code>a = 0.0
b = 2.0 * np.pi

y_0 = np.zeros(2)
y_0[0] = 0.0
y_0[1] = 1.0
nv = 2

tolerance = 1.0e-6

#perform the integration
x,y = rk4_mv(dydx,a,b,y_0,tolerance)<jupyter_output><empty_output><jupyter_text>Plot the result<jupyter_code>plt.plot(x,y[:,0],'o',color='y',label='y(x)')
plt.plot(x,y[:,1],'o',color='silver',label='dydx(x)')
xx = np.linspace(0,2.0*np.pi,1000)
plt.plot(xx,np.sin(xx),color='k',label='sin(x)')
plt.plot(xx,np.cos(xx),color='purple',label='cos(x)')
plt.xlabel('x',color='chocolate',fontsize='large')
plt.ylabel('y, dy/dx',color='chocolate',fontsize='large')
plt.legend(frameon=False,fontsize='medium')<jupyter_output><empty_output><jupyter_text>Plot the error Notice the errors will exceed our "tolerance".<jupyter_code>sin_X = np.sin(x)
cos_X = np.cos(x)

y_error = (y[:,0]-sin_X)
dydx_error = (y[:,1]-cos_X)

plt.plot(x, y_error, label="y(x) Error", color='fuchsia')
plt.plot(x, dydx_error, label="dydx(x) Error", color='aqua')
plt.legend(frameon=False, fontsize='medium')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13252858958068614, 'fraction_numerical': 0.03672172808132147, 'mean_word_length': 2.2578642384105962, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 42, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '10515', 'license': 'MIT', 'n_tokens_mistral': 3284, 'n_tokens_neox': 3065, 'n_words': 676, 'path': 'rojenkin.ipynb', 'repo_name': 'rojenkin/astr-119_hw-5', 'size': '7870'}