<jupyter_start><jupyter_text>coupled derivatives<jupyter_code>def dydx(x,y):
    #where dydx = z and dzdx = -y
    #set y = y[0] and z = y[1]
    
    #declare array
    y_derivs = np.zeros(2)
    
    #dydx= z
    y_derivs[0] = y[1]
    
    #dzdx = -y
    y_derivs[1] = -1*y[0]
    
    return y_derivs<jupyter_output><empty_output><jupyter_text>4th order RK method<jupyter_code>def rk4_mv_core(dydx,xi,yi,nv,h):
    
    #declare k4 arrays
    #nv = number of variables
    #h = width
    
    k1 = np.zeros(nv)
    k2 = np.zeros(nv)
    k3 = np.zeros(nv)
    k4 = np.zeros(nv)
    k5 = np.zeros(nv)
    k6 = np.zeros(nv)
    
    #define x at half step
    x_ipoh = xi + 0.5*h
    
    #define x at 1 step
    x_ipo = xi + h
    
    #declare a temporary y array
    y_temp = np.zeros(nv)
    
    #k1 values
    y_derivs = dydx(xi,yi)
    k1[:] = h*y_derivs[:]
    
    #k2
    y_temp[:] = yi[:] + (1/5)*k1[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k2[:] = h*y_derivs[:]
    
    #k3
    y_temp[:] = yi[:] + (9/40)*k2[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k3[:] = h*y_derivs[:]
    
    #k4
    y_temp[:] = yi[:] + (6/5)*k3[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k4[:] = h*y_derivs[:]
    
    #k5
    y_temp[:] = yi[:] + (35/27)*k4[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k5[:] = h*y_derivs[:]
    
    #k6
    y_temp[:] = yi[:] + (253/4096)*k5[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k6[:] = h*y_derivs[:]
    
    yipo = yi + h*((2825/27648)*k1 + (18575/48384)*k3 + (13525/55296)*k4 + (277/14336)*k5 + (1/4)*k6)
    
    #yout(i)=y(i)+h*(C1*dydx(i)+C3*ak3(i)+C4*ak4(i)+C6*ak6(i))
    
    return yipo<jupyter_output><empty_output><jupyter_text>adaptive step driver<jupyter_code>def rk4_mv_ad(dydx,x_i,y_i,nv,h,tol):
    
    #safety scale
    SAFETY    = 0.9
    H_NEW_FAC = 2.0
    
    imax = 10000
    
    i = 0
    
    #create error
    Delta = np.full(nv,2*tol)
    
    h_step = h
    
    #adjust step
    while(Delta.max()/tol > 1.0):
        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_i,nv,0.5*h_step)
        
        #error
        Delta = np.fabs(y_2 - y_11)
        
        #if error is too lardge, take smaller step
        if(Delta.max()/tol > 1.0):
            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
        
    #for bigger steps in the future
    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>wrapper function<jupyter_code>def rk4_mv(dydx,a,b,y_a,tol):
    #y_a are boundary conditions
    
    #starting step
    xi = a
    yi = y_a.copy()
    
    #an initial step size (small)
    h = 1.0e-4 * (b-a)
    
    imax = 10000
    
    i = 0
    
    nv = len(y_a)
    
    #initial conditions
    x = np.full(1,a)
    y = np.full((1,nv),y_a)
    
    flag = 1
    
    while(flag):
        
        #y_i+1
        yi_new, h_new, h_step = rk4_mv_ad(dydx,xi,yi,nv,h,tol)
        
        #update step
        h = h_new
        
        if(xi+h_step>b):
            h = b-xi
            
            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 step to 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
        
        if(i>=imax):
            print("Maximum iterations reached.")
            raise StopIteration("Iteration number = ",i)
            
        #iterate
        i += 1
        
        #info outputs
        s = "i = %3d\tx = %9.8f\th = %9.8f\tb = %9.8f" % (i,xi,h_step,b)
        print(s)
        
        #break if the new xi is == b
        if(xi==b):
            flag = 0
            
    return x,y<jupyter_output><empty_output><jupyter_text>perform 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

x,y = rk4_mv(dydx,a,b,y_0,tolerance)
plt.plot(x,y[:,0],'.',label='y(x)')
plt.plot(x,y[:,1],'.',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)
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.167833607907743, 'fraction_numerical': 0.04242174629324547, 'mean_word_length': 1.9944512946979038, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 26, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '101297', 'license': 'MIT', 'n_tokens_mistral': 2365, 'n_tokens_neox': 2221, 'n_words': 356, 'path': 'hw-5.ipynb', 'repo_name': 'nncastil/astr-119-hw-5', 'size': '4856'}