Analysis script for numerical solution of differential equation

#import libraries

import numpy as np

import scipy as sp

from ipywidgets import interact

from scipy.integrate import odeint

import matplotlib.pyplot as plt

%matplotlib inline

Fit data

# Solving the Newton's law model and also equations 17a and 17b from Lorenceau paper:

# Newton's law model

def DZ_dt_Newton(Z, t, args):

    cc = 0

    h = args[0]

    g = args[1]

    b = args[2]

    return [Z[1], -Z[1]**2 / (Z[0] + cc) -g * Z[0] / (Z[0] + cc)  + g * h / (Z[0] + cc) - b * Z[1] / (Z[0] + cc) - cc * g / (Z[0] + cc)]

def plot_osc(h = 11.2, g = 9.8e2, b = 23):

    # prepare data for plotting:

    z_data1 = -z_data[:] + h # change the overall level so that bottom of straw is z=0

    time_axis1 = time_data[:] # only include data for positive times (after cap is released)

    params = (h, g, b)    

    # solve Newton model:

    t_soln = time_axis1

    Z_soln_Newton = sp.integrate.odeint(DZ_dt_Newton, [0.02, 0], t_soln, args=(params,))   

    z_soln_Newton = Z_soln_Newton[:, 0]      # fluid height

    plt.clf()

    plt.plot(time_axis1,z_data1, 'b.', label = 'Data') 

    plt.plot(t_soln, z_soln_Newton, 'r', label = 'Newtonian model')

    plt.xlabel('time (sec)', fontsize = 15) 

    plt.ylabel('fluid level (cm)', fontsize = 15)

    plt.rc('xtick', labelsize = 10)

    plt.rc('ytick', labelsize = 10)

    plt.title('Fluid level oscillations with both models\nh=%2.2f, b=%2.2e\nfilename = %s'%(h,b,filename), fontsize = 15)

    plt.legend(frameon = False, loc = 1)

    plt.xlim([-0.2, 3])

    plt.grid()

    plt.show()

interact(plot_osc, h=(0.0, 20.0), g=(5.0e2, 15.0e2), b = (0.0, 30.0));