Modeling hanging chain using numerical methods in python

Question

I am trying to model a hanging chain using numerical methods (this pde ) by following the example here for a wave on a string. I am having trouble getting a good approximation from my differencing schemes and I am not sure why. In particular, the true solution and the model do not agree at the first time step, which they should.

Here is my code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import j0, jn_zeros #Bessel Functions

#Variables we change:
mode = 1 #Which natural mode 
tend = 1000 #Final time 
dt = 0.1 #time step

#Defining our constants:
g = 9.81 #gravitational acceleration (m/s^2)
L = 10 #Length of the chain (m)
A = 0.05 #Maximum amplitude (m)
x_end = 200
x = np.linspace(0,L,x_end) #The vertical axis (m)
bessel_scale = 2 * np.sqrt(x/g) #Scaling the vertical axis aka inside bessel equation
w = jn_zeros(0, mode)[mode-1] * np.sqrt(g/L) / 2 #Omega? - calculating the first n zeroes
    #of the zeroeth order Bessel function

dx = L/x_end
pi = np.pi

#-----------------------------Approximation Methods----------------------------------#
c1 = g * x * dt**2 / dx**2
c2 = g * dt**2 / dx

nT = 200#1200 # Number of time points
N = np.size(x)#201 # Number of position points
y_fwd = np.zeros([N, nT])
y_ctr = np.zeros([N, nT])
y_bck = np.zeros([N, nT])

y_rk4 = np.zeros([N,nT])

# Loop over initial condition
for i in range(N-1): #initial condition based on the position point x(within bessel scale
    y_fwd[i,0] = A * j0(w*bessel_scale[i])
    y_ctr[i,0] = A * j0(w*bessel_scale[i])
    y_bck[i,0] = A * j0(w*bessel_scale[i])
    y_rk4[i,0] = A * j0(w*bessel_scale[i])


#Hanging Chain Function
for t in range(tend):
    hc = A * j0(w*bessel_scale) * np.cos(w*t)#Bessel_scale depends on x
    
    for i in range(N-1): #position
        y_fwd[i,int(t+1)] = -y_fwd[i,t-1] + y_fwd[i+1,t]*(c1[i]+c2) - y_fwd[i,t]*(2*c1[i]+2+c2) + y_fwd[i-1,t]*(c1[i]) #Forward Method
        y_ctr[i,int(t+1)] = -y_ctr[i,t-1] + y_ctr[i+1,t]*(c1[i]+.5*c2) - y_ctr[i,t]*(2*c1[i]+2) + y_ctr[i-1,t]*(c1[i]-.5*c2) #Centered Method
        y_bck[i,int(t+1)] = -y_bck[i,t-1] + y_bck[i+1,t]*(c1[i]) - y_bck[i,t]*(2*c1[i]+2-c2) + y_bck[i-1,t]*(c1[i] -c2) #Forward Method

    plt.xlabel('Horizontal Position (x)')
    plt.xlim((-2*A, 2*A))
    plt.ylabel('Vertical Position (y)')
    plt.ylim((0,L))
    plt.title('Hanging Chain: mode = %i' %mode)

    plt.plot(y_fwd[:,int(t-2)],x, label='Forward Method')
    plt.plot(y_ctr[:,int(t-2)],x, label='Centered Method')
    plt.plot(y_bck[:,int(t-2)],x, label='Backward Method')
    plt.plot(hc,x, 'k', label='True Solution')

    plt.legend()

    plt.pause(0.1)
    
    plt.clf()
plt.show()

I expected that the solutions match better and that they at least start in the same position. I am not sure how to proceed.

Modeling hanging chain using numerical methods in python

0 Answers0