diff --git a/doc/src/week43/programs/dill.py~ b/doc/src/week43/programs/dill.py~ new file mode 100644 index 000000000..c5805a983 --- /dev/null +++ b/doc/src/week43/programs/dill.py~ @@ -0,0 +1,203 @@ +import numpy as np +import matplotlib.pyplot as plt +import tensorflow as tf +from math import * +import time +import tensorflow.keras +# Different methods from Keras needed to create an RNN +# This is not necessary but it shortened function calls +# that need to be used in the code. +from tensorflow.keras import datasets, layers, models +from tensorflow.keras.layers import Input +from tensorflow.keras import regularizers +from tensorflow.keras import optimizers +from tensorflow.keras.models import Model, Sequential +#from tensorflow.keras.layers.core import Dense, Activation +from tensorflow.keras.layers import Dense, SimpleRNN, LSTM, GRU + +# Define Analytical, Euler-Cromer, and Velocity-Verlet methods of solving +def analytical(k,m,x0,v0,dt,tfinal): + t = np.arange(0,tfinal+dt,dt) + v = -x0 * np.sin(t) + v0 * np.cos(t) + x = x0 * np.cos(t) + v0 * np.sin(t) + K = 1/2 *m*v**2 + U = 1/2 *k*x**2 + return x, v, K, U, t + +def euler_cromer(k,m,x0,v0,dt,tfinal): + n = ceil(tfinal/dt) # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + x = np.zeros(n) + K = np.zeros(n) + U = np.zeros(n) + # Define Initial Conditions + x[0] = x0 + v[0] = v0 + K[0] = 1/2 *m*v0**2 + U[0] = 1/2 *k*x0**2 + + # Integrate using the Euler-Cromer Method + for i in range(n-1): + a = -x[i] + v[i+1] = v[i] + dt*a + x[i+1] = x[i] + dt*v[i+1] + K[i+1] = 1/2 *m*v[i+1]**2 + U[i+1] = 1/2 *k*x[i+1]**2 + t[i+1] = t[i] + dt + return x, v, K, U, t + +def velocity_verlet(k,m,x0,v0,dt,tfinal): + n = ceil(tfinal/dt) + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + x = np.zeros(n) + K = np.zeros(n) + U = np.zeros(n) + # Define Initial Conditions + x[0] = x0 + v[0] = v0 + K[0] = 1/2 *m*v0**2 + U[0] = 1/2 *k*x0**2 + + # Integrate using the Velocity-Verlet Method + for i in range(n-1): + a = -x[i] + x[i+1] = x[i] + dt*v[i] + dt**2 /2*a + a1 = -x[i+1] + v[i+1] = v[i] + dt/2*(a+a1) + K[i+1] = 1/2 *m*v[i+1]**2 + U[i+1] = 1/2 *k*x[i+1]**2 + t[i+1] = t[i] + dt + return x, v, K, U, t + + +# Analytical Solution to position +def ana_cos(r0,t,k=1,m=1): + w0 = sqrt(k/m) + return r0*np.cos(w0*t) + +# Trial Function for Neural Net +def trial_func(x,y,y0=1): + return x*y + y0 + +# Loss Function for Position and Velocity +def combined_right(trialv,trialy,k=1,m=1): + return -k/m*trialy, trialv + +# Loss Wrapper in order to pass the Loss Function to Neural Net +def con_loss_wrapper(input_tensor): + def con_loss_function(y,y_pred): + trialy = trial_func(input_tensor,y_pred) + trialv = trial_func(input_tensor,y_pred) + righty, rightv = combined_right(trialv,trialy) + + leftv = tf.gradients(trialv,input_tensor)[0] + lefty = tf.gradients(leftv,input_tensor)[0] + + + loss = tf.reduce_mean((tf.math.squared_difference(lefty,righty))) + return loss + return con_loss_function + +# Creates the input data for the Neural Net +def create_input_data(x0=0,xmax=1,num_batch=5,len_batch=15): + input_data = np.linspace(x0,xmax,num_batch*len_batch) + input_data = input_data.reshape(num_batch,len_batch) + + return input_data + +# Creates the Neural Net +def create_net(data,len_batch,lr,epochs,right_side,loss,n_hidden_layer=50): + + input_tensor = Input(shape=(len_batch,)) + + hidden1 = Dense(n_hidden_layer,activation='tanh', + kernel_initializer='random_uniform',bias_initializer='random_uniform')(input_tensor) + hidden2 = Dense(n_hidden_layer,activation='tanh', + kernel_initializer='random_uniform',bias_initializer='random_uniform')(hidden1) + hidden3 = Dense(n_hidden_layer,activation='tanh', + kernel_initializer='random_uniform',bias_initializer='random_uniform')(hidden2) + hidden4 = Dense(n_hidden_layer,activation='tanh', + kernel_initializer='random_uniform',bias_initializer='random_uniform')(hidden3) + output = Dense(len_batch)(hidden4) + + model = Model(input_tensor,output) + + gd = optimizers.SGD(lr=lr) # May need to change first 'lr' to 'learning_rate' depending on TF/Keras version + model.compile(loss=loss(input_tensor),optimizer=gd) + model.fit(data,np.zeros((data.shape[0])),epochs=epochs) + + res = model.predict(data) + + del model + + return res + +# Define Function for Mean Squared Error +def mean_squared_error(analytical, results): + mse = 0 + for i in range(len(analytical)): + mse += (analytical[i] - results[i])**2 + + mse = mse/len(analytical) + return mse + + +# Define Constants +dt = 0.01 +tfinal = 50 +x0 = 1 +v0 = 0 +m = 1 +k = 1 + +num_batch = 1000 +len_batch = 1 + +# Create input data +data = create_input_data(x0=0,xmax=10,num_batch=num_batch,len_batch=len_batch) + +# Create and Run the Neural Net +nn_start_time = time.time() + +velocity = create_net(data,len_batch=len_batch,lr=0.001,n_hidden_layer=50, + epochs=1000,right_side=combined_right,loss=con_loss_wrapper) + +nn_end_time = time.time() + + +# Reshape Neural Net Output for easy graphing and analysis +n = num_batch*len_batch +velocity = velocity.reshape(1,n) +t = data.reshape(1,n)[0] +results_v = trial_func(t,velocity)[0] + +# Create Comparison data from Analytical Solution +analyt = ana_cos(1,t,m=1,k=1) + + + +# Euler-Cromer +ec_start_time = time.time() +ec_x = euler_cromer(k=1,m,x0,v0,dt,tfinal)[0] +ec_end_time = time.time() + +# Velocity-Verlet +vv_start_time = time.time() +vv_x = velocity_verlet(k,m,x0,v0,dt,tfinal)[0] +vv_end_time = time.time() + +# Plot +plt.plot(t,analyt,label='analytical') +plt.plot(t,results_v,label='net') +plt.plot(t, ec_x, label = "Euler-Cromer") +plt.plot(t, vv_x, label = "Velocity-Verlet") +plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + +# Find Mean Squared Error +print("The Mean Squared Error of the Neural Net Solution is", mean_squared_error(analyt, results_v), "with a runtime of", nn_end_time-nn_start_time,"seconds.") +print("The Mean Squared Error of the Euler-Cromer is", mean_squared_error(analyt, ec_x), "with a runtime of", ec_end_time-ec_start_time,"seconds.") +print("The Mean Squared Error of the Velocity-Verlet is", mean_squared_error(analyt, vv_x), "with a runtime of", vv_end_time-vv_start_time,"seconds.") + diff --git a/doc/src/week43/programs/odetf.py b/doc/src/week43/programs/odetf.py new file mode 100644 index 000000000..f2f78c758 --- /dev/null +++ b/doc/src/week43/programs/odetf.py @@ -0,0 +1,251 @@ +import numpy as np +import matplotlib.pyplot as plt +import tensorflow as tf +from math import * +import time +# User-friendly machine learning library +# Front end for TensorFlow +import tensorflow.keras +# Different methods from Keras needed to create an RNN +# This is not necessary but it shortened function calls +# that need to be used in the code. +from tensorflow.keras import datasets, layers, models +from tensorflow.keras.layers import Input +from tensorflow.keras import regularizers +from tensorflow.keras.models import Model, Sequential +#from tensorflow.keras.layers.core import Dense, Activation +from tensorflow.keras.layers import Dense, SimpleRNN, LSTM, GRU +# For timing the code +from timeit import default_timer as timer +# For plotting +import matplotlib.pyplot as plt + + +# Define Analytical, Euler-Cromer, and Velocity-Verlet methods of solving +def analytical(k,m,x0,v0,dt,tfinal): + t = np.arange(0,tfinal+dt,dt) + v = -x0 * np.sin(t) + v0 * np.cos(t) + x = x0 * np.cos(t) + v0 * np.sin(t) + K = 1/2 *m*v**2 + U = 1/2 *k*x**2 + return x, v, K, U, t + +def euler_cromer(k,m,x0,v0,dt,tfinal): + n = np.ceil(tfinal/dt)\ + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + x = np.zeros(n) + K = np.zeros(n) + U = np.zeros(n) + # Define Initial Conditions + x[0] = x0 + v[0] = v0 + K[0] = 1/2 *m*v0**2 + U[0] = 1/2 *k*x0**2 + + # Integrate using the Euler-Cromer Method + for i in range(n-1): + a = -x[i] + v[i+1] = v[i] + dt*a + x[i+1] = x[i] + dt*v[i+1] + K[i+1] = 1/2 *m*v[i+1]**2 + U[i+1] = 1/2 *k*x[i+1]**2 + t[i+1] = t[i] + dt + return x, v, K, U, t + +def velocity_verlet(k,m,x0,v0,dt,tfinal): + n = np.ceil(tfinal/dt) + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + x = np.zeros(n) + K = np.zeros(n) + U = np.zeros(n) + # Define Initial Conditions + x[0] = x0 + v[0] = v0 + K[0] = 1/2 *m*v0**2 + U[0] = 1/2 *k*x0**2 + + # Integrate using the Velocity-Verlet Method + for i in range(n-1): + a = -x[i] + x[i+1] = x[i] + dt*v[i] + dt**2 /2*a + a1 = -x[i+1] + v[i+1] = v[i] + dt/2*(a+a1) + K[i+1] = 1/2 *m*v[i+1]**2 + U[i+1] = 1/2 *k*x[i+1]**2 + t[i+1] = t[i] + dt + return x, v, K, U, t + +# Define Constants +dt = 0.01 +tfinal = 50 +x0 = 1 +v0 = 0 +m = 1 +k = 1 + +# Call the Integration Function +ax, av, aK, aU, at = analytical(k,m,x0,v0,dt,tfinal) +ecx, ecv, ecK, ecU, ect = euler_cromer(k,m,x0,v0,dt,tfinal) +vvx, vvv, vvK, vvU,vvt = velocity_verlet(k,m,x0,v0,dt,tfinal) + +# Plots +fig, axes = plt.subplots(1,3, figsize = (15,5)) +fig.suptitle("System of an Undamped Spring", fontsize=14, y = 1.05) + +axes[0].plot(at, ax, label = "Analytical Method") +axes[0].plot(ect, ecx, label = "Euler-Cromer Method") +axes[0].plot(vvt, vvx, label = "Velocity-Verlet Method") +axes[0].set_title("Position as a Function of Dimensionless Time") +axes[0].set_xlabel("Dimensionless Time") +axes[0].set_ylabel("Position") + +axes[1].plot(at, av, label = "Analytical Method") +axes[1].plot(ect, ecv, label = "Euler-Cromer Method") +axes[1].plot(vvt, vvv, label = "Velocity-Verlet Method") +axes[1].set_title("Velocity as a Function of Dimensionless Time") +axes[1].set_xlabel("Dimensionless Time") +axes[1].set_ylabel("Velocity") + +axes[2].plot(at, aU+aK, label = "Analytical Method") +axes[2].plot(ect, ecU+ecK, label = "Euler-Cromer Method") +axes[2].plot(vvt, vvU+vvK, label = "Velocity-Verlet Method") +axes[2].set_title("Energy as a Function of Dimensionless Time") +axes[2].set_xlabel("Dimensionless Time") +axes[2].set_ylabel("Total Energy") + +plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + +plt.tight_layout() + + +def damp(gamma,m,x0,v0,DeltaT,tfinal): + n = np.ceil(tfinal/DeltaT) + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + r = np.zeros(n) + # Define Initial Conditions + r[0] = x0 + v[0] = v0 + # Integrate over using the Velocity Verlet Method + for i in range(n-1): + a = -r[i] - 2*gamma*v[i] + r[i+1] = r[i] + DeltaT*v[i] + DeltaT**2 /2 *a + a1 = -r[i+1] - 2*gamma*v[i] + v[i+1] = v[i] + DeltaT/2*(a+a1) + t[i+1] = t[i] + DeltaT + + return t,r,v + +# Define Constants +under_gamma = 0.1 +crit_gamma = 1 +over_gamma = 2 +m = 1 +x0 = 1 +v0 = 0 +dt = 0.1 +tfinal = 50 + +# Call the Integration Function +t, under_r, under_v = damp(under_gamma,m,x0,v0,dt,tfinal) +t, crit_r, crit_v = damp(crit_gamma,m,x0,v0,dt,tfinal) +t, over_r, over_v = damp(over_gamma,m,x0,v0,dt,tfinal) + + +# Plots +fig, axes = plt.subplots(1,2, figsize = (15,5)) +fig.suptitle("System of a Damped Spring", fontsize=14, y = 1.05) + +axes[0].plot(t, under_r, label = "Under Damping") +axes[0].plot(t, crit_r, label = "Critical Damping") +axes[0].plot(t, over_r, label = "Over Damping") +axes[0].set_title("Position as a Function of Dimensionless Time") +axes[0].set_xlabel("Dimensionless Time") +axes[0].set_ylabel("Position") + +axes[1].plot(t, under_v, label = "Under Damping") +axes[1].plot(t, crit_v, label = "Critical Damping") +axes[1].plot(t, over_v, label = "Over Damping") +axes[1].set_title("Velocity as a Function of Dimensionless Time") +axes[1].set_xlabel("Dimensionless Time") +axes[1].set_ylabel("Velocity") + +plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + +plt.tight_layout() + + +def forced(gamma,r0,v0,F0,omega,d,DeltaT,tfinal): + n = np.ceil(tfinal/DeltaT) + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + r = np.zeros(n) + # Define Initial Conditions + r[0] = r0 + v[0] = v0 + + # Integrate using the 4th-Order RK Method + for i in range(n-1): + t[i+1] = t[i] + DeltaT + + Force = (-r[i] - 2*gamma*v[i] - F0*np.cos(omega*t[i]-d))*m + k1x = DeltaT*v[i] + k1v = DeltaT*Force + + vv = v[i]+k1v*0.5*DeltaT + rr = r[i]+k1x*0.5*DeltaT + Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m + k2x = DeltaT*vv + k2v = DeltaT*Force + + vv = v[i]+k2v*0.5*DeltaT + rr = r[i]+k2x*0.5*DeltaT + Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m + k3x = DeltaT*vv + k3v = DeltaT*Force + + vv = v[i]+k3v*DeltaT + rr = r[i]+k3x*DeltaT + Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m + k4x = DeltaT*vv + k4v = DeltaT*Force + + r[i+1] = r[i]+(k1x+2*k2x+2*k3x+k4x)/6. + v[i+1] = v[i]+(k1v+2*k2v+2*k3v+k4v)/6. + + return t,r,v + +# Define Constants +gamma = 0.1 +r0 = 1 +v0 = 0 +F0 = 3 +omega = 3 +d = 0 +dt = 0.1 +tfinal = 100 + +# Call the Integration Function +t, r, v = forced(gamma,r0,v0,F0,omega,d,dt,tfinal) + +# Plots +fig, axes = plt.subplots(1,2, figsize = (15,5)) +fig.suptitle("System of a Damped Spring with a Driving Force", fontsize=14, y = 1.05) + +axes[0].plot(t, r) +axes[0].set_title("Position as a Function of Dimensionless Time") +axes[0].set_xlabel("Dimensionless Time") +axes[0].set_ylabel("Position") + +axes[1].plot(t, v) +axes[1].set_title("Velocity as a Function of Dimensionless Time") +axes[1].set_xlabel("Dimensionless Time") +axes[1].set_ylabel("Velocity") + +plt.tight_layout() diff --git a/doc/src/week43/programs/odetf.py~ b/doc/src/week43/programs/odetf.py~ new file mode 100644 index 000000000..d974171c6 --- /dev/null +++ b/doc/src/week43/programs/odetf.py~ @@ -0,0 +1,250 @@ +# For matrices and calculations +import numpy as np +# For machine learning (backend for keras) +import tensorflow as tf +# User-friendly machine learning library +# Front end for TensorFlow +import tensorflow.keras +# Different methods from Keras needed to create an RNN +# This is not necessary but it shortened function calls +# that need to be used in the code. +from tensorflow.keras import datasets, layers, models +from tensorflow.keras.layers import Input +from tensorflow.keras import regularizers +from tensorflow.keras.models import Model, Sequential +#from tensorflow.keras.layers.core import Dense, Activation +from tensorflow.keras.layers import Dense, SimpleRNN, LSTM, GRU +# For timing the code +from timeit import default_timer as timer +# For plotting +import matplotlib.pyplot as plt + + +# Define Analytical, Euler-Cromer, and Velocity-Verlet methods of solving +def analytical(k,m,x0,v0,dt,tfinal): + t = np.arange(0,tfinal+dt,dt) + v = -x0 * np.sin(t) + v0 * np.cos(t) + x = x0 * np.cos(t) + v0 * np.sin(t) + K = 1/2 *m*v**2 + U = 1/2 *k*x**2 + return x, v, K, U, t + +def euler_cromer(k,m,x0,v0,dt,tfinal): + n = np.ceil(tfinal/dt)\ + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + x = np.zeros(n) + K = np.zeros(n) + U = np.zeros(n) + # Define Initial Conditions + x[0] = x0 + v[0] = v0 + K[0] = 1/2 *m*v0**2 + U[0] = 1/2 *k*x0**2 + + # Integrate using the Euler-Cromer Method + for i in range(n-1): + a = -x[i] + v[i+1] = v[i] + dt*a + x[i+1] = x[i] + dt*v[i+1] + K[i+1] = 1/2 *m*v[i+1]**2 + U[i+1] = 1/2 *k*x[i+1]**2 + t[i+1] = t[i] + dt + return x, v, K, U, t + +def velocity_verlet(k,m,x0,v0,dt,tfinal): + n = np.ceil(tfinal/dt) + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + x = np.zeros(n) + K = np.zeros(n) + U = np.zeros(n) + # Define Initial Conditions + x[0] = x0 + v[0] = v0 + K[0] = 1/2 *m*v0**2 + U[0] = 1/2 *k*x0**2 + + # Integrate using the Velocity-Verlet Method + for i in range(n-1): + a = -x[i] + x[i+1] = x[i] + dt*v[i] + dt**2 /2*a + a1 = -x[i+1] + v[i+1] = v[i] + dt/2*(a+a1) + K[i+1] = 1/2 *m*v[i+1]**2 + U[i+1] = 1/2 *k*x[i+1]**2 + t[i+1] = t[i] + dt + return x, v, K, U, t + +# Define Constants +dt = 0.01 +tfinal = 50 +x0 = 1 +v0 = 0 +m = 1 +k = 1 + +# Call the Integration Function +ax, av, aK, aU, at = analytical(k,m,x0,v0,dt,tfinal) +ecx, ecv, ecK, ecU, ect = euler_cromer(k,m,x0,v0,dt,tfinal) +vvx, vvv, vvK, vvU,vvt = velocity_verlet(k,m,x0,v0,dt,tfinal) + +# Plots +fig, axes = plt.subplots(1,3, figsize = (15,5)) +fig.suptitle("System of an Undamped Spring", fontsize=14, y = 1.05) + +axes[0].plot(at, ax, label = "Analytical Method") +axes[0].plot(ect, ecx, label = "Euler-Cromer Method") +axes[0].plot(vvt, vvx, label = "Velocity-Verlet Method") +axes[0].set_title("Position as a Function of Dimensionless Time") +axes[0].set_xlabel("Dimensionless Time") +axes[0].set_ylabel("Position") + +axes[1].plot(at, av, label = "Analytical Method") +axes[1].plot(ect, ecv, label = "Euler-Cromer Method") +axes[1].plot(vvt, vvv, label = "Velocity-Verlet Method") +axes[1].set_title("Velocity as a Function of Dimensionless Time") +axes[1].set_xlabel("Dimensionless Time") +axes[1].set_ylabel("Velocity") + +axes[2].plot(at, aU+aK, label = "Analytical Method") +axes[2].plot(ect, ecU+ecK, label = "Euler-Cromer Method") +axes[2].plot(vvt, vvU+vvK, label = "Velocity-Verlet Method") +axes[2].set_title("Energy as a Function of Dimensionless Time") +axes[2].set_xlabel("Dimensionless Time") +axes[2].set_ylabel("Total Energy") + +plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + +plt.tight_layout() + + +def damp(gamma,m,x0,v0,DeltaT,tfinal): + n = np.ceil(tfinal/DeltaT) + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + r = np.zeros(n) + # Define Initial Conditions + r[0] = x0 + v[0] = v0 + # Integrate over using the Velocity Verlet Method + for i in range(n-1): + a = -r[i] - 2*gamma*v[i] + r[i+1] = r[i] + DeltaT*v[i] + DeltaT**2 /2 *a + a1 = -r[i+1] - 2*gamma*v[i] + v[i+1] = v[i] + DeltaT/2*(a+a1) + t[i+1] = t[i] + DeltaT + + return t,r,v + +# Define Constants +under_gamma = 0.1 +crit_gamma = 1 +over_gamma = 2 +m = 1 +x0 = 1 +v0 = 0 +dt = 0.1 +tfinal = 50 + +# Call the Integration Function +t, under_r, under_v = damp(under_gamma,m,x0,v0,dt,tfinal) +t, crit_r, crit_v = damp(crit_gamma,m,x0,v0,dt,tfinal) +t, over_r, over_v = damp(over_gamma,m,x0,v0,dt,tfinal) + + +# Plots +fig, axes = plt.subplots(1,2, figsize = (15,5)) +fig.suptitle("System of a Damped Spring", fontsize=14, y = 1.05) + +axes[0].plot(t, under_r, label = "Under Damping") +axes[0].plot(t, crit_r, label = "Critical Damping") +axes[0].plot(t, over_r, label = "Over Damping") +axes[0].set_title("Position as a Function of Dimensionless Time") +axes[0].set_xlabel("Dimensionless Time") +axes[0].set_ylabel("Position") + +axes[1].plot(t, under_v, label = "Under Damping") +axes[1].plot(t, crit_v, label = "Critical Damping") +axes[1].plot(t, over_v, label = "Over Damping") +axes[1].set_title("Velocity as a Function of Dimensionless Time") +axes[1].set_xlabel("Dimensionless Time") +axes[1].set_ylabel("Velocity") + +plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + +plt.tight_layout() + + +def forced(gamma,r0,v0,F0,omega,d,DeltaT,tfinal): + n = np.ceil(tfinal/DeltaT) + # Set up arrays + t = np.zeros(n) + v = np.zeros(n) + r = np.zeros(n) + # Define Initial Conditions + r[0] = r0 + v[0] = v0 + + # Integrate using the 4th-Order RK Method + for i in range(n-1): + t[i+1] = t[i] + DeltaT + + Force = (-r[i] - 2*gamma*v[i] - F0*np.cos(omega*t[i]-d))*m + k1x = DeltaT*v[i] + k1v = DeltaT*Force + + vv = v[i]+k1v*0.5*DeltaT + rr = r[i]+k1x*0.5*DeltaT + Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m + k2x = DeltaT*vv + k2v = DeltaT*Force + + vv = v[i]+k2v*0.5*DeltaT + rr = r[i]+k2x*0.5*DeltaT + Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m + k3x = DeltaT*vv + k3v = DeltaT*Force + + vv = v[i]+k3v*DeltaT + rr = r[i]+k3x*DeltaT + Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m + k4x = DeltaT*vv + k4v = DeltaT*Force + + r[i+1] = r[i]+(k1x+2*k2x+2*k3x+k4x)/6. + v[i+1] = v[i]+(k1v+2*k2v+2*k3v+k4v)/6. + + return t,r,v + +# Define Constants +gamma = 0.1 +r0 = 1 +v0 = 0 +F0 = 3 +omega = 3 +d = 0 +dt = 0.1 +tfinal = 100 + +# Call the Integration Function +t, r, v = forced(gamma,r0,v0,F0,omega,d,dt,tfinal) + +# Plots +fig, axes = plt.subplots(1,2, figsize = (15,5)) +fig.suptitle("System of a Damped Spring with a Driving Force", fontsize=14, y = 1.05) + +axes[0].plot(t, r) +axes[0].set_title("Position as a Function of Dimensionless Time") +axes[0].set_xlabel("Dimensionless Time") +axes[0].set_ylabel("Position") + +axes[1].plot(t, v) +axes[1].set_title("Velocity as a Function of Dimensionless Time") +axes[1].set_xlabel("Dimensionless Time") +axes[1].set_ylabel("Velocity") + +plt.tight_layout()