added codes
This commit is contained in:
mhjensen
@@ -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.")
|
||||
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
Reference in New Issue
Block a user