<jupyter_start><jupyter_text>A notebook to perform Cash-Karp 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>Define our coupled derivativs to integrate<jupyter_code>def dydx(x,y):
        
    # set the derivatives
        
    # our equation is d^2x/dx^2 = -y
        
    # so we can write
    # dydx = z
    # dzdx = -y
        
    # we will set y = y[0]   (the function y)
    # we will set z = y[1]   (the function x)
        
    # declare an array
    y_derivs = np.zeros(2)  # array of functions
        
    # set dydx = z
    y_derivs[0] = y[1]      # dydx we put in y_derivs[1]
        
    # set dydx = -y
    y_derivs[1] = -1*y[0]
        
    # here we have to return an array
    return y_derivs<jupyter_output><empty_output><jupyter_text>The 5th order Cash-Karp Runge-Kutta method<jupyter_code># now we'll operate on all the elements at the same time 

def rk5_mv_core(dydx,xi,yi,nv,h):   #dydx is function we just wrote, 
                                        # xi is value of x at step i
                                        # xi changes by value of h
                                        # nv is number of variables
                                        # yi is the array 

    # coefficient values from adaptive step-size control
    # c:               
    ci = np.array([0, 1/5, 3/10, 3/5, 1, 7/8])                    
    # b:
    bi = np.array([37/378, 0, 250/621, 125/594, 0, 512/1771])   
    # dbi:
    dbi = np.array([2825/27648, 0, 18575/48384, 13525/55296, 277/14336, 1/4])   
   
    ai = np.array([[0,      0,       0,         0,            0,        0],
               [1/5,        0,       0,         0,            0,        0],
               [3/40,       9/40,    0,         0,            0,        0],
               [3/10,      -9/10,    6/5,       0,            0,        0],
               [-11/54,     5/2,    -70/27,     35/27,        0,        0],
               [1631/55296, 175/512, 575/13824, 44275/110592, 253/4096, 0]]) 
    
    # declare k? arrays
    k1 = np.zeros(nv)               # each k is array of two emelents one that
    k2 = np.zeros(nv)               # corresponds to y and one that corresponds to z
    k3 = np.zeros(nv)
    k4 = np.zeros(nv)
    k5 = np.zeros(nv)
    k6 = np.zeros(nv)
    
    #declare a temp y array
    y_temp = np.zeros(nv)
    
    # get k1 values
    # k1 = hf(xn, yn)
    y_derivs = dydx(xi,yi)
    k1[:] = h*y_derivs[:]
    
    
    # get k2 values
    #     k2 = h*f (  xn +  c2*h, yn +  a21*k1)
    y_temp[:] = yi[:] + ai[1,0]*k1[:]
    y_derivs = dydx(xi + ci[1]*h, y_temp[:])  
    k2[:] = h*y_derivs[:]  
 
    
    # get k3 values
    #     k3 = h*f (   xn + c3*h, yn +     a31*k1     +   a32*k2)
    y_temp[:] = yi[:] + ai[2,0]*k1[:] + ai[2,1]*k2[:] 
    y_derivs = dydx(xi + ci[2]*h, y_temp[:] )  
    k3[:] = h*y_derivs[:]  

    
    # get k4 values
    #     k4 = h*f (   xn + c4*h, yn +   a41*k1       +   a42*k2       +   a43*k3)
    y_temp[:] = yi[:] + ai[3,0]*k1[:] + ai[3,1]*k2[:] + ai[3,2]*k3[:]
    y_derivs = dydx(xi + ci[3]*h, y_temp[:] ) 
    k4[:] = h*y_derivs[:]

    
    # get k5 values
    #     k5 = h*f (   xn + c5*h, yn +   a51*k1       +    a52*k2      +     a53*k3     +   a54*k4)
    y_temp[:] = yi[:] + ai[4,0]*k1[:] + ai[4,1]*k2[:] + ai[4,2]*k3[:] +ai[4,3]*k4[:]
    y_derivs = dydx(xi + ci[4]*h, y_temp[:]) 
    k5[:] = h*y_derivs[:]

    
    # get k6 values
    #     k6 = h*f (   xn + c6*h, yn +   a61*k1       +    a62*k2      +    a63*k3      +   a64*k4        +   a65*k5)
    y_temp[:] = yi[:] + ai[5,0]*k1[:] + ai[5,1]*k2[:] + ai[5,2]*k3[:] + ai[5,3]*k4[:] + ai[5,4]*k5[:]
    y_derivs = dydx(xi + ci[5]*h, y_temp[:] ) 
    k6[:] = h*y_derivs[:]

    
    # yn+1 = yn +   b1*k1  +    b2*k2    +   b3*k3  +  b4*k4   +  b5*k5   +   b6*k6 
    yipo =   yi + bi[0]*k1 + bi[1]*k2 + bi[2]*k3 + bi[3]*k4 + bi[4]*k5 + bi[5]*k6
    
    
    # embedded fourth-order formula is
    dyipo =   yi[:] + dbi[0]*k1[:] + dbi[1]*k2[:] + dbi[2]*k3[:] + dbi[3]*k4[:] + dbi[4]*k5[:] + dbi[5]*k6[:]
   
    
    Delta = np.fabs(yipo - dyipo)
    
    return yipo, Delta<jupyter_output><empty_output><jupyter_text>Adaptive step size driver for Cash-Karp Runge-Kutta method<jupyter_code># dydx functions that take drivatives
# x_i value of x at step i
# y_i values of items in array at i
# h is step size
# tol tolerance
# this function  
#def rk5_mv_ad(dydx,x_i,y_i,h,tol):
def rk5_mv_ad(dydx, x_i, y_i, nv, h, tol):
    
    # accuracy check
    #yscal = np.abs(y_i[:]) + np.abs(dydx[:]*h) + 1e-3
    
    #define safety scale
    SAFETY    = 0.9
    H_NEW_FAC = 2.0
    
    # set a maximum number of iterations
    imax = 10000
    
    # set an iteration variable 
    i = 0
    
    #y_derivs = dydx(x_i,y_i)
    
    # step size, remember the step
    h_step = h
    
    # create an error
    y_2, Delta = rk5_mv_core(dydx, x_i, y_i, nv, h_step)
    
    # adjust steps until error is in our tolerance
    while( (Delta.max()/tol) > 1.0 ):
         
        #y_new, yerror = rk5_mv_core(x_i, y_i, y_derivs, h_step, nv)

        # error too big, reduce step size
        h_step *= SAFETY * (Delta.max()/tol)**(-0.25)
        
        y_2, Delta = rk5_mv_core(dydx, x_i, y_i, nv, h_step)
        # if step too small: 
        #if(np.abs(h_step_new) < 0.1*np.abs(h_step)):
        #    h_step_new = 0.1 * h_step
        #h_step = h_step_new
        
        if( i>=imax):
            print("Too many iterations in rk5_mv_ad()")
            raise StopIteration("Ending i =",i)
        i += 1
        
    # larger step next time
    h_new = np.fmin(h_step * (Delta.max()/tol)**(-0.9), h_step * H_NEW_FAC)
            
        
    return y_2, h_new, h_step<jupyter_output><empty_output><jupyter_text>Define a wrapper for RK5<jupyter_code># this function, you pass in the derivatives you want to evolve
# the starting place, the ending place, and tolerance
# it will call the two functions we just wrote 
# we don't know how many steps we'll take so we have initial conditions
# we just provide some tolerance and want for this wrapper to try to stay
# within the tolerance

def rk5_mv(dydx,a,b,y_a,tol):
    
    # dfdx is the derivative wrt x (copied that from slide, idk what wrt x means though)
    # a,b is the lower bound,upper bound
    # y_a are the initial conditions. (the boundary conditions)
    # tol is the tolerance for integrating y
    
    xi = a          # starting step
    yi = y_a.copy()
    
    h = 1.0e-4 * (b-a)   # starting h, our adaptive function will further change this. Start small
                            # note: this might work fine if integrating backwards, it will just be a negative h
    
    imax = 10000   # reasonably we will probably never reach this, it is here however to prevent infinite loops
    
    i = 0    # set our iteration variable
    
    nv = len(y_a)  # number of dimensions used will be the length of y_a
    
    x = np.full(1,a)        # why not np.array([a]) is it to just match syntax of y below?
    y = np.full((1,nv),y_a) # 1 x nv matrix to put our different variables into
    
    flag = 1   # arbitrary switch for the while loop. Obviously switched on to run the the loop. Note: 1 = on, 0 = off
    
    while(flag):
        yi_new, h_new, h_step = rk5_mv_ad(dydx, xi, yi, nv, h, tol)  # see function above to know what this does
        
        h = h_new  # replace our original guess above
        
        # once we reach the end, we want to prevent an overshoot
        if(xi+h_step>b):
            h = b-xi   # make the last step perfectly reach b
            
            yi_new, h_new, h_step = rk5_mv_ad(dydx, xi, yi, nv, h, tol) # recalc y_i+1
            
            flag = 0  # since we will only be here when we reach the end, turn off the while loop
            
        # if we have not reached the end, this is how we do the next step
        
        xi += h_step        # change xi,yi with the output of rk4_mv_ad
        yi[:] = yi_new[:]
        
        x = np.append(x,xi) # add the new x to the x_array
        
        # it's a little more complicated to set up adding the new i. This is because we are changing the dimensions
        y_new = np.zeros((len(x),nv))  # set up the shell of y
        y_new [0:len(x)-1,:] = y       # fill the first part of the shell with the previous y values to preserve them
        y_new[-1,:] = yi[:]            # put our newest value at the end of the shell
        del y                          # delete y in order to
        y = y_new                      # replace y with the now filled shell
        
        if(i>=imax):    # prevent inifinite loop
            print('Maximum iterations reached.')
            raise StopIteration('Max iterations set to ',i)
            
        i += 1   # end of the first step, so increment i and go through the loop again until the flag above
        
        s = 'i = %3d\tx = %9.8f\th = %9.8f\tb = %9.8f' % (i,xi,h_step,b)
        print(s)   # print each step.
        if(xi==b):
            flag = 0 # makes sure to turn off the while loop
    return x,y  # return answer<jupyter_output><empty_output><jupyter_text>Perform the integration<jupyter_code>a = 0.0
b = 2.0 * np.pi

# initial conditions
y_0 = np.zeros(2)  # array of size 2
y_0[0] = 0.0
y_0[1] = 1.0
#nv = 2

tolerance = 1.0e-6

# perform the integration 
x, y = rk5_mv(dydx,a,b,y_0,tolerance)<jupyter_output><empty_output><jupyter_text>Plot the result<jupyter_code>plt.plot(x, y[:,0], 'o', label='y(x)')
plt.plot(x, y[:,1], 'o', label='dydx(x)')

xx = np.linspace(0,2.0 * np.pi, 1000)

plt.plot(xx, np.sin(xx), label='sin(x)')
plt.plot(xx, np.cos(xx), label='cos(x)')

plt.xlabel('x')
plt.ylabel('y, dy/dx')
plt.legend(frameon=False)<jupyter_output><empty_output><jupyter_text>Plot the error<jupyter_code>sine = np.sin(x)
cosine = np.cos(x)

y_error = (y[:,0] - sine)
dydx_error = (y[:,1] - cosine)

plt.plot(x, y_error, label='y(x) Error')
plt.plot(x, dydx_error, label='dydx(x) Error')
plt.legend(frameon=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13027613412228797, 'fraction_numerical': 0.04358974358974359, 'mean_word_length': 1.9309248554913294, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '653596', 'license': 'MIT', 'n_tokens_mistral': 4187, 'n_tokens_neox': 3915, 'n_words': 1091, 'path': 'hw-5.ipynb', 'repo_name': 'Veronika131/astr-119-hw-5', 'size': '10140'}