# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt


def dydx(x,y):
    y_derivs = np.zeros(2)
    
    y_derivs[0] = y[1]
    y_derivs[1] = -1*y[0]
    
    return y_derivs


# ## Core

def rk4_mv_core(dydx,xi,yi,nv,h):
    k1 = np.zeros(nv)
    k2 = np.zeros(nv)
    k3 = np.zeros(nv)
    k4 = np.zeros(nv)
    
    x_ipoh = xi +0.5*h
    
    x_ipo = xi + h
    
    y_temp = np.eros(nv)
    
    y_derivs = dydx(xi,yi)
    k1[:] = h*y_derivs[:]
    
    y_temp[:] = yi[:] + 0.5*k1[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k2[:] =h*y_derivs[:]
    
    y_temp[:] = yi[:] + 0.5*k2[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k3[:] =h*y_derivs[:]
    
    y_temp[:] = yi[:] + 0.5*k3[:]
    y_derivs = dydx(x_ipoh,y_temp)
    k4[:] =h*y_derivs[:]
    
    yipo = yi + (k1 + 2*k2 + 2*k3 + k4)/6.
    
    return yipo


# ##Driver 

def rk4_mv_ad(dydx,x_i,y_i,nv,h,tol):
    SAFETY = 0.9
    H_NEW_FAC = 2.0
    
    imax = 10000
    
    i = 0
    
    Delta = np.full(nv,2*tol)
    
    h_step = h
    
    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)
        
        Delta = np.fabs(y_2 - y_11)
        
        if(delta.max()/tol >1.0):
            h_step *= SAFETY * (Delta.max()/tol)**(-0.25)
            
        if(i>=imax):
            print("too many iterations in rk4_mv_ad()")
            raise StopIteration("Ending after i = ",i)
        
        i+=1
    h_new = np.fmin(h_step * (Delta()/tol)**(-0.9), h_step*H_NEW_FAC)
    
    return y_2, h_new, h_step


# ## Wrapper

def_rk4_mv(dydx,a,b,y_a,tol):
    
    xi = a
    yi = y_a.copy()
    
    h = 1.0e-4 * (b-a)
    
    imax = 10000
    i = 0
    nv = len(y_a)
    x = pn.full(1,a)
    y = np.full((1,nv),y_a)
    flag = 1
    while(flag):
        yi_new, h_new, h_step = rk4_mv_ad(dydx,xi,yi,nv,h,tol)
        h = n_new
        if(ii+h_step>b):
            h = b-xi
            yi_new, h_new, h_step = rk4_mv_ad(dydx,xi,yi,nv,h,tol)
            flag = 0
        

StarcoderdataJupyterScriptsDedupFiltered