Update week4/week4_2.py, week4/README.MD, week4/week4_4.py, week4/week4_3.py

week4/README.MD

+# Week 4 - Ordinary Differential Equations
+
+From 6.2, error at different step sizes for Euler integration (both are diverging) and Runge-Kutta with much smaller error.
+
+![](img/6_2_euler_v_runge.png)
+
+From 6.3, adaptive Runge-Kutta stepper's step size increasing to meet a desired error target. Could have used a more aggressive minimum error term to see the steps oscillate with the derivative instead of just meeting the requirements of the peak derivative.
+
+![](img/adaptive_runge_step_v_targetError.png)

week4/week4_2.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+'''
+6.2)
+For a simple harmonic oscillator y'' + y = 0, with initial conditions y(0) = 1, y'(0) = 0,
+find y(t) from t = 0 to 100π. Use an Euler method and a fixed-step fourth-order Runge-Kutta method.
+For each method check how the average local error, and the error in the final value and slope, depend on the step size
+'''
+
+# Euler method
+def euler_integrate(step):
+    t = np.arange(0, 100*np.pi, step)
+    increments = np.arange(1, len(t), 1) # we will need something indexed from 1 to final time_step to not interfere with the initial conditions
+    y = np.zeros(len(t))
+    y[0] = 1
+    dy = np.zeros(len(t))
+    dy[0] = 0
+    for time in increments:
+        y[time] = y[time-1] + step*dy[time-1]
+        dy[time] = dy[time-1] - step*(y[time-1])
+        # dy[time] = dy[time-1] - step*(y[time]) # reverse euler? Just use the expected new y 
+    return y, dy, t
+
+# Runge-Kutta method
+def rk4_integrate(step):
+    t = np.arange(0, 100*np.pi, step)
+    increments = np.arange(1, len(t), 1) # we will need something indexed from 1 to final time_step to not interfere with the initial conditions
+    y = np.zeros(len(t))
+    y[0] = 1
+    dy = np.zeros(len(t))
+    dy[0] = 0
+    for time in increments:
+        k1, k2, k3, k4 = [0, 0], [0, 0], [0, 0], [0, 0] # make position and velocity vectors
+        
+        k1[0] = step * dy[time-1] # normal euler step
+        k1[1] = - step * y[time-1]
+        k2[0] = step * (dy[time-1] + k1[1]*0.5) # half euler step in the direction of k1
+        k2[1] = -step * (y[time-1] + k1[0]*0.5)
+        k3[0] = step * (dy[time-1] + k2[1]*0.5) # half euler step in the direction of k2
+        k3[1] = -step * (y[time-1] + k2[0]*0.5)
+        k4[0] = step * (dy[time-1] + k3[1]) # full euler step in the direction of k3
+        k4[1] = -step * (y[time-1] + k3[0])
+
+        y[time] = y[time-1] + (k1[0] + 2* k2[0] + 2 * k3[0] + k4[0])/6 # divide by 6 (not 4) because we assigned 2x value to k2 and k3
+        dy[time] = dy[time-1] + (k1[1] + 2* k2[1] + 2 * k3[1] + k4[1])/6
+    return y, dy, t
+
+step1 = 0.001
+step2 = 0.01
+
+# y1, dy1, t1 = euler_integrate(step1)
+# y2, dy2, t2 = euler_integrate(step2)
+
+y1, dy1, t1 = rk4_integrate(step1)
+y2, dy2, t2 = rk4_integrate(step2)
+
+fig, ax = plt.subplots(2)
+ax[0].plot(t1, y1)
+ax[0].title.set_text('Times step = ' + str(step1) + '\n Final y: ' + str(round(y1[-1],3)) + '\n Final dy: ' + str(round(dy1[-1],3)))
+ax[1].plot(t2, y2)
+ax[1].title.set_text('Times step = ' + str(step2) + '\n Final y: ' + str(round(y2[-1],3)) + '\n Final dy: ' + str(round(dy2[-1],3)))
+plt.tight_layout()
+plt.show()
\ No newline at end of file

week4/week4_3.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+'''
+6.3)
+Write a fourth-order Runge Kutta adaptive stepper for the preceding problem,
+and check how the average step size that it finds depends on the desired local error.
+'''
+
+# Adaptive Runge-Kutta method (rewritten to be a step call)
+def rk4_step(step, y, dy):
+    k1, k2, k3, k4 = [0, 0], [0, 0], [0, 0], [0, 0] # make position and velocity vector
+    k1[0] = step * dy # normal euler step for velocity
+    k1[1] = - step * y  # normal euler step for position
+    k2[0] = step * (dy + k1[1]*0.5) # half euler step in the direction of k1
+    k2[1] = -step * (y + k1[0]*0.5)
+    k3[0] = step * (dy + k2[1]*0.5) # half euler step in the direction of k2
+    k3[1] = -step * (y + k2[0]*0.5)
+    k4[0] = step * (dy + k3[1]) # full euler step in the direction of k3
+    k4[1] = -step * (y + k3[0])
+    y_new = y + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6 # divide by 6 (not 4) because we assigned 2x value to k2 and k3
+    dy_new = dy + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])/6
+    return y_new, dy_new
+
+def rk4_adaptive(step_0, step_nudge, error_min, error_max):
+    ''' In this case we don't know how many timesteps we will take, so it's best to
+    start with lists and convert to numpy arrays at the end...
+    '''
+    t_final = 100*np.pi
+    t = [0]
+    y = [1]
+    dy = [0]
+    steps = [0.1]
+    
+    step = steps[0]
+    time = step
+    count = 1
+
+    while time < t_final:
+        y_full, dy_full = rk4_step(step, y[count-1], dy[count-1])
+        y_half1, dy_half1 = rk4_step(step*0.5, y[count-1], dy[count-1])
+        y_half2, dy_half2 = rk4_step(step*0.5, y_half1, dy_half1)
+
+        error = np.abs(y_full - y_half2) # error is how closely the full and half steps agree
+        # if we are too high or too low, nudge the step size and try again
+        if error > error_max:
+            step *= step_nudge
+        elif error < error_min:
+            step *= 1 - (step_nudge-1)
+        # otherwise go ahead and append and increment time and count
+        else:
+            y.append(y_full)
+            dy.append(dy_full)
+            steps.append(step)
+            time += step
+            t.append(time)
+            count += 1
+    return t, y, dy, steps
+
+step_0 = 0.1
+nudge = 0.99 # increase or decrease the step by this percentage to nudge the step size
+error_min =  1e-6 # easier to choose if theses are informed by something physical
+error_max = 1e-4
+
+t, y, dy, steps = rk4_adaptive(step_0, nudge, error_min, error_max)
+
+fig, ax = plt.subplots(2)
+ax[0].plot(t, y)
+ax[0].title.set_text('\n Max Error: ' + str(round(error_max, 10)) + '\n Min Error: ' + str(round(error_min, 10)))
+ax[1].plot(t, steps)
+ax[1].title.set_text('Total number of steps: ' + str(len(t)) + '\n Average step size: ' + str(round(np.mean(steps), 6)))
+plt.tight_layout()
+plt.show()

week4/week4_4.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+'''
+6.4)
+Numerically solve the differential equation found in Problem 4.3:
+l*ddθ + (g + ddz)*sinθ = 0
+Take the motion of the platform to be periodic, and interactively explore the dynamics of
+the pendulum as a function of the amplitude and frequency of the excitation.
+
+Consider a pendulum with a mass m connected by a massless rod of length l to a moveable platform at height z (Figure 6.4)
+'''
+step = 0.001
+t_final = 10*np.pi
+t = np.arange(0, t_final, step)
+increments = np.arange(1, len(t), 1)
+time = 0
+
+# Platform Parameters (periodic motion and very differentiable)
+A = 1 # amplitude of the platform's motion [m]
+platform_omega = 1 # angular frequency of the platform's motion [rad/s]
+zs = A*np.sin(platform_omega*t) # all platform positions as a function of time
+ddzs = A*platform_omega**2*np.cos(platform_omega*t) # all platform accelerations as a function of time
+
+# Pendulum Parameters
+thetas = np.zeros(len(t))
+d_thetas = np.zeros(len(t))
+dd_thetas = np.zeros(len(t))
+l = 1 # length of pendulum
+g = 9.81 # acceleration due to gravity
+theta_0 = 0 # initial angle of the pendulum [rad]
+d_theta_0 = 1 # initial angular velocity of the pendulum [rad/s]
+
+def d_theta(theta, omega, t):
+    return omega
+
+def dd_theta(theta, omega, t):
+    solution = -np.sin(theta)/l * (g - ddzs[int(t/step)])
+    print(solution)
+    return solution
+
+def rk4_step(prev_values, dif_eqs, step, time):
+    variables = np.array(prev_values)
+    k1, k2, k3, k4 = [0, 0], [0, 0], [0, 0], [0, 0] # make position and velocity vector
+    ddy = dif_eqs[1](prev_values[0], prev_values[1], time)
+    dy = dif_eqs[0](prev_values[0], prev_values[1], time)
+    print(ddy)
+    k1[0] = step * ddy # normal euler step for velocity
+    k1[1] = - step * dy  # normal euler step for position
+    k2[0] = step * (ddy + k1[1]*0.5) # half euler step in the direction of k1
+    k2[1] = -step * (dy + k1[0]*0.5)
+    k3[0] = step * (ddy + k2[1]*0.5) # half euler step in the direction of k2
+    k3[1] = -step * (dy + k2[0]*0.5)
+    k4[0] = step * (ddy + k3[1]) # full euler step in the direction of k3
+    k4[1] = -step * (dy + k3[0])
+    variables = variables + np.array(k1)/6 + np.array(k2)/3 + np.array(k3)/3 + np.array(k4)/6
+    return variables
+
+for increment in increments:
+    time = t[increment]
+    thetas[increment], d_theta[increment] = rk4_step([thetas[increment-1], d_thetas[increment-1]], [d_theta, dd_theta], step, time)
+
+# let's plot the chaotic theta over the periodic z
+plt.plot(t, zs)
+plt.plot(t, thetas)
\ No newline at end of file