3524.sagews

# TO DO
#-Figure out plausible units for initial values and constants
#-Determine correct coefficient for arbitrary order RK (i.e. RK4 -> 1/6, RK5 -> ?)
#-Try to find a way to make RK1 = Euler method and RK2 = Midpoint method
#-Dynamic number of steps (# of oscillations, etc.)
#-International System units!

@interact
def solver(method_choice = selector(['Euler', 'Midpoint', 'Runge-Kutta'], default='Runge-Kutta', nrows=1, label='Method'),
           step_size = slider(vmin=0.0000001, vmax=0.1, step_size=0.0000001, default=1, label='Step size'),
           rk_order = input_box(default="4", label='Order of Runge-Kutta'),
           init_radius = input_box(default="1", label='Initial radius'),
           init_rad_vel = input_box(default="0", label='Initial radial velocity'),
           init_rad_acc = input_box(default="0", label='Initial radial acceleration'),
           f_bubble_pressure = input_box(default="1", label='Initial gas pressure inside bubble', width=10),
           f_v_pressure = input_box(default="1", label='Vapor pressure (p_v(T_inf))', width=10),
           f_i_pressure = input_box(default="1", label='Fluid pressure far away from cavitation', width=10),
           f_density = input_box(default="1", label='Liquid density', width=10),        #f_var indicates fluid variable
           f_viscosity = input_box(default="1", label='Liquid viscosity', width=10),
           f_surface_tension = input_box(default="1", label='Bubble surface tension', width=10),
           f_constant_k = input_box(default="1", label='Constant K', width=10)):
        
        h = step_size
        rk_order = int(rk_order)
        
        R_0 = int(init_radius)
        RV_0 = int(init_rad_vel)
        RA_0 = int(init_rad_acc)
        PB = int(f_bubble_pressure)
        PV = int(f_v_pressure)
        PI = int(f_i_pressure)
        D = int(f_density)
        V = int(f_viscosity)
        S = int(f_surface_tension)
        K = int(f_constant_k)
        
        #Function for linearized Rayleigh-Plesset
        def equation(r, rv):
            return sin(r)#return rv*((3*RV_0/r) + (4*V/r^2)) + (R_0^K)/(r^(3*K+1)) + (2*S/(3*r^3)) + ((PB-PV-PI)/D) - (3*RV_0^2)/(2*r)
        
        #Create some arrays to store values
        R_array = [[0,R_0]]
        RV_array = [[0,RV_0]]
        RA_array = [[0,equation(R_0,RV_0)]]
        slope_array = []
        extrema_array = []
        trend = True #True = ascending, False = descending
        osc_goal = 10
        osc_num = 0
        
        #and variables to use in computations
        r = R_0
        rv = RV_0
        ra = RA_0
        
        i = 1
        
        #BEGIN MAIN LOOP
        while i < 1000 or osc_goal < osc_num:
            if method_choice == 'Euler':
                #calculate the radius based on previous radius and velocity
                r = r + rv*h
                R_array.append([i*h,r])
                #calculate the velocity based on acceleration 'ra' and previous velocity
                rv = rv + ra*h
                RV_array.append([i*h,rv])
                #find acceleration based on linearized equation function
                ra = equation(r,rv)
                RA_array.append([i*h,ra])
            elif method_choice == 'Midpoint':
                r = r + rv*h
                rv = rv + ra*h
                ra = equation(r,rv)
                ra = equation(r+h/2, rv+ra/2)
                #append to arrays
                RA_array.append([i*h,ra])
                RV_array.append([i*h,rv])
                R_array.append([i*h,r])
            elif method_choice == 'Runge-Kutta':
                r = r + rv*h
                R_array.append([i*h,r])
                rv = rv + ra*h
                RV_array.append([i*h,rv])
                
                #calculation of 'ra' for arbitrary RK order with delta values stored in 'k_array'
                k_array = []
                k_array.append(h*equation(r,rv))
                for a in range(rk_order-2):
                    k_array.append(2*h*equation(r+h/2, rv+k_array[a]/2))
                k_array.append(h*equation(r+h,rv+k_array[rk_order-2]))
                
                ra = ra + sum(k_array)/6
                RA_array.append([i*h,ra])
            
            #calculate slope and detect extrema
            slope = (R_array[i][1] - R_array[i-1][1])/h
            slope_array.append(slope.n())  
            if len(extrema_array) == 0:
                trend = R_array[1][1] > R_array[0][1]    #on first iteration, set 'trend'
            elif R_array[i][1] > R_array[i-1][1] and trend == False:    #detect the new slope; if it changes, declare an extrema and switch 'trend'
                trend = True
                extrema_array.append(R_array[i-1])
            elif R_array[i][1] < R_array[i-1][1] and trend == True:    #same as above but for the opposite direction
                trend = False
                extrema_array.append(R_array[i-1])
            
            #check for low amounts of activity (overflow safeguard #1)
            avg_change = 0
            tot_change = 0    #not really the 'total' change - takes sign changes into account, unlike avg_change
            l = len(slope_array)
            if l > 5:
                for a in range(5):
                    avg_change += abs(slope_array[l-a-1] - slope_array[l-a-2])
                    tot_change += slope_array[l-a-1] - slope_array[l-a-2]
                avg_change /= 5
                tot_change /= 5
                if avg_change < 0.000005:
                    print "Computation interrupted due to lack of activity\n\t(Average change = %s)" % avg_change
                    break
            if avg_change > 100 and avg_change == abs(tot_change):
                print "Computation interrupted due to exponential behavior"
                break
            
            i += 1
        #END MAIN LOOP
        
        #Create a line that connects all of the points in r_array
        L = line([(R_array[j][0], R_array[j][1]) for j in range(len(R_array))], color=(1,0,0), legend_label="Radial position")
        show(L)
        
        #Same thing, RV_array
        M = line([(RV_array[j][0], RV_array[j][1]) for j in range(len(RV_array))], color=(0,1,0), legend_label="Radial velocity")
        show(M)
        
        #ctrl-c ctrl-v RA_array
        N = line([(RA_array[j][0], RA_array[j][1]) for j in range(len(RA_array))], color=(0,0,1), legend_label="Radial accel")
        show(N)
        show(L+M+N)

Method  Step size  Order of Runge-Kutta  Initial radius  Initial radial velocity  Initial radial acceleration  Initial gas pressure inside bubble  Vapor pressure (p_v(T_inf))  Fluid pressure far away from cavitation  Liquid density  Liquid viscosity  Bubble surface tension  Constant K  [removed] [removed] [removed] [removed]

#Implement this code into the main bubble simulation
#-Convert main simulations to 'while' loops
#-Add in the overflow/nonactivity detector
#-Create function to evaluate an extreme point given a slice of the array
#-HAVE THE RANGES FOR OVERFLOW DETECTION SCALE WITH STEP SIZE

#Finish this code
#-Account for a triple peak with middle peak being lowest
#-Implement run-time extrema evaluation (see below)

array = [[0,0]]
slope_array = []
extrema_array = []
n = 0.1 #step size
trend = True #True = ascending, False = descending
osc_goal = 10
osc_num = 0

i = 1

while osc_num < osc_goal:
    #calculate new x, y, and slope; add them to respective arrays
    x = i*n
    y = sin(x) + 0.1*sin(10*x)# + 0.2*sin(20*x)# + 0.3*sin(30*x)
    array.append([x,y])
    slope = (array[i][1] - array[i-1][1])/(4*pi/n)
    slope_array.append(slope.n())
    
    #detect extrema
    if i == 0:
        trend = y > array[0][1]    #on first iteration, set 'trend'
    elif array[i][1] > array[i-1][1] and trend == False:    #detect the new slope; if it changes, declare an extrema and switch 'trend'
        trend = True
        extrema_array.append(array[i-1])
    elif array[i][1] < array[i-1][1] and trend == True:    #same as above but for the opposite direction
        trend = False
        extrema_array.append(array[i-1])
    
    #guess at a typical length of oscillation
    #possible things to look for: distance traveled in y direction 
    sorted(extrema_array)
    
    #check for low amounts of activity or exponential behavior
    avg_change = 0
    tot_change = 0    #not really the 'total' change - takes sign changes into account, unlike avg_change
    l = len(slope_array)
    if l > 5:
        for a in range(5):
            avg_change += abs(slope_array[l-a-1] - slope_array[l-a-2])
            tot_change += slope_array[l-a-1] - slope_array[l-a-2]
        avg_change /= 5
        tot_change /= 5
        if avg_change < 0.000005:
            print "Computation interrupted due to lack of activity\n\t(Average change = %s)" % avg_change
            break
        if avg_change > 100 and avg_change == abs(tot_change):
            print "Computation interrupted due to exponential behavior"
            break
     
    i += 1

#Convert this into a run-time evaluation of extrema
#REMEMBER: You can't have 2 consecutive maxima - must be a minimum in between!
#Think about comparing y-distance of extrema to known min and max

    extrema_array_2 = [] #this will hold all the "significant" extrema
    prev_value = ""
    for a in range(len(extrema_array)-2):
        if a < 2: continue
    
        search_point = list(extrema_array[a])
        prev_point = list(extrema_array[a-1])
        prev_ext = list(extrema_array[a-2])
        next_ext = list(extrema_array[a+2])
    
        if (search_point[1] > prev_point[1] and search_point[1] >= prev_ext[1] and search_point[1] >= next_ext[1]):
            search_point.append("max")
            if (len(extrema_array_2) != 0):
                if (extrema_array_2[len(extrema_array_2)-1][2] == "min"):
                    extrema_array_2.append(search_point)
                elif (search_point[1].n() > extrema_array_2[len(extrema_array_2)-1][1].n()):
                    extrema_array_2.pop()
                    extrema_array_2.append(search_point)
                elif (search_point[1].n() < extrema_array_2[len(extrema_array_2)-1][1].n()):
                    pass
            else:
                extrema_array_2.append(search_point)
        elif (search_point[1] < prev_point[1] and search_point[1] <= prev_ext[1] and search_point[1] <= next_ext[1]):
            search_point.append("min")
            if (len(extrema_array_2) > 0):
                if (extrema_array_2[len(extrema_array_2)-1][2] == "max"):
                    extrema_array_2.append(search_point)
                elif (search_point[1].n() < extrema_array_2[len(extrema_array_2)-1][1].n()):
                    extrema_array_2.pop()
                    extrema_array_2.append(search_point)
                elif (search_point[1].n() > extrema_array_2[len(extrema_array_2)-1][1].n()):
                    pass
            else:
                extrema_array_2.append(search_point)
    
    osc_num = len(extrema_array_2)/2

print len(extrema_array_2)
graph = line([(array[j][0].n(), array[j][1].n()) for j in range(len(array))])
points = point([(extrema_array_2[j][0].n(), extrema_array_2[j][1].n()) for j in    range(len(extrema_array_2))], rgbcolor=(1,0,0), pointsize=18)
points2 = point([(extrema_array[j][0].n(), extrema_array[j][1].n()) for j in range(len(extrema_array))], rgbcolor=(0,1,0), pointsize=18)
show(points2 + points + graph)

20

#DETERMINING 'SIGNIFICANT' EXTREMA
#-An extrema is significant if the two surrounding extrema of the same type (minimum/maximum) are both higher or both lower, respectively
#-An extrema is not significant if its neighboring extrema is of the same type and of greater magnitude
#-An extrema is *probably* not significant if it the change in Y between it and the neighboring extrema of the opposite type is significantly less than that of the furthest separated pair 
#    -Thought: Margin of acceptability defined by: (distance between absolute min and max) - (distance between furthest apart extrema)

def runge_kutta_two(a,b,N,alpha):
    #create functions we will need
    f = lambda x,y: 1.0 + (y/x) + (y/x)^2    #this is the first equation
    g = lambda x: x*tan(ln(x))               # this is the exact solution for the IVP
    #continuing with page 279
    h = (b-a)/N
    t = a
    omega = alpha
    #make a list to store values
    w = []
    w.append((0,omega, numerical_approx(abs(g(1)-omega))))    # this stores y(1) = 0 & actual error
    for i in range(1,N):
        omega = omega + h * f(t + (h/2), omega + (h/t)*f(t,omega))        # compute omega_i
        w.append( ( i, omega, numerical_approx(abs(g(1)-omega))) )
        #store values
        t = a + i * h     # compute t_i
        return w

print runge_kutta_two(1,100,100,0)

#a = start point
#b = endpoint
#N = number of steps
#alpha = 

#h = step size

if method_choice == 'Euler':
            for i in range(math.ceil(n/h)):
                #calculate the radius based on previous radius and velocity
                r = r + rv*h
                R_array.append([i*h,r])
                #calculate the velocity based on acceleration 'ra' and previous velocity
                rv = rv + ra*h
                RV_array.append([i*h,rv])
                #find acceleration based on linearized equation function
                ra = equation(r,rv)
                RA_array.append([i*h,ra])
        elif method_choice == 'Midpoint':
            for i in range(math.ceil(n/h)):
                r = r + rv*h
                rv = rv + ra*h
                ra = equation(r,rv)
                ra = equation(r+h/2, rv+ra/2)
                #append to arrays
                RA_array.append([i*h,ra])
                RV_array.append([i*h,rv])
                R_array.append([i*h,r])
        elif method_choice == 'Runge-Kutta':
            for i in range(math.ceil(n/h)):
                r = r + rv*h
                R_array.append([i*h,r])
                rv = rv + ra*h
                RV_array.append([i*h,rv])
                
                #calculation of 'ra' for arbitrary RK order with delta values stored in 'k_array'
                k_array = []
                k_array.append(h*equation(r,rv))
                for a in range(rk_order-2):
                    k_array.append(2*h*equation(r+h/2, rv+k_array[a]/2))
                k_array.append(h*equation(r+h,rv+k_array[rk_order-2]))
                
                ra = ra + sum(k_array)/6
                RA_array.append([i*h,ra])

Traceback (most recent call last): rv = rv + ra*h File "", line 1, in <module> File "/tmp/tmpL_gBWK/___code___.py", line 14 elif method_choice == 'Midpoint': ^ IndentationError: unindent does not match any outer indentation level