0

I am trying to analyse a wave on a string by solving the wave equation with Python. Here are my requirements for the solution.

1) I model reflective ends by using much larger masses on first and last point on the string -> Large inertia

2)No spring on edges. Then k[0] and k[-1] will be ZERO.

I have problem with indices. In my 2nd loop I get y[i,j-1], k[i-1], y[i-1,j]. In the first iteration the loop will then use y[0,-1], k[-1], y[-1,0]. These are my last points and not my first points. How can I avoid this problem?

N.J Abel
  • 51
  • 2
  • 8

1 Answers1

2

What you need is initiating your mass array with one additional element. I mean

...
m = [mi]*N # mass array   [!!!] instead of (N-1)  [!!!]
...

Idem for your springs

...
k = [ki]*N
...

Consequently, you can keep k[0] equal to 10. since you model reflective ends. You may thus want to comment or drop this line

...
##k[0] = 0
...


For aesthetic considerations you may want to fill the gap at the end of the x-axis. In this case, simply do
N = 201 # Number of mass points

Your code thus becomes

from numpy import *
from  matplotlib.pyplot import *

# Variables
N = 201 # Number of mass points
nT = 1200 # Number of time points
mi = 0.02 # mass in kg
m = [mi]*N # mass array
m[-1] = 100 # Large last mass reflective edges
m[0] = 100 # Large first mass  reflective edges
ki = 10.#spring 
k = [ki]*N
k[-1] = 0
dx = 0.2
kappa = ki*dx
my = mi/dx
c = sqrt(kappa/my) # velocity 

dt = 0.04

#  3 vectors
x = arange( N )*dx # x points
t = arange( N )*dt # t points


y = zeros( [N, nT ] )# 2D array
# Loop over initial condition
for i in range(N-1):
    y[i,0] = sin((7.*pi*i)/(N-1)) # Initial condition dependent on mass point


# Iterating over time and position to find next position of wave
for j in range(nT-1):

    for i in range(N-1):
        y[i,j+1] = 2*y[i,j] - y[i,j-1] + (dt**2/m[i])*(k[i-1]*y[i+1,j] -2*k[i-1]*y[i,j] + k[i]*y[i-1,j] )

    #check values of edges
    print y[:2,j+1],y[-2:,j+1]

    # Creates an animation
    cla()
    ylabel("Amplitude")
    xlabel("x")    
    ylim(-10,10)
    plot(x,y[:,j-2])
    pause(0.001)
close()

which produces

enter image description here


Following your comment, I think that if you want to see the wave traveling along the string before reflection, you should not initiate conditions everywhere (in space). I mean, e.g., doing
...
# Loop over initial condition
for i in range(N-1):
    ci_i = sin(7.*pi*i/(N-1)) # Initial condition dependent on mass point
    if np.sign(ci_i*y[i-1,0])<0:
        break
    else:
        y[i,0] = ci_i
...

produces

enter image description here

New attempt after answers:

from numpy import *
from  matplotlib.pyplot import *

N = 201
nT = 1200
mi = 0.02
m = [mi]*(N)
m[-1] = 1000
m[0] = 1000
ki = 10.
k = [ki]*N
dx = 0.2
kappa = ki*dx
my = mi/dx
c = sqrt(kappa/my) 
dt = 0.04


x = arange( N )*dx
t = arange( N )*dt


y = zeros( [N, nT ] )

for i in range(N-1):
    y[i,0] = sin((7.*pi*i)/(N-1))



for j in range(nT-1):

    for i in range(N-1):
        if j == 0: # if j = 0 then ...  y[i,j-1]=y[i,j]
            y[i,j+1] = 2*y[i,j] - y[i,j] + (dt**2/m[i])*(k[i-1]*y[i+1,j] -2*k[i-1]*y[i,j] + k[i]*y[i-1,j] )
        else:
            y[i,j+1] = 2*y[i,j] - y[i,j-1] + (dt**2/m[i])*( k[i-1]*y[i+1,j] -2*k[i-1]*y[i,j] + k[i]*y[i-1,j] )

    cla()

    ylim(-1,1)
    plot(x,y[:,j-2])
    pause(0.0001)

    ylabel("Amplitude")
    xlabel("x")

    print len(x), len(t), N,nT

Here is a plot of the new attempt at solution with |amplitude| of anti node equal 1.0. Will this do anything with further solving the issue with indices?

enter image description here

keepAlive
  • 6,369
  • 5
  • 24
  • 39
  • I got the same plot with an earlier attempt, but I ignored the results. I think the amplitude is too large + we don't see the wave traveling along the string before reflection. – N.J Abel May 07 '17 at 11:20
  • If we set m[0] = 0.02 and m[-1] = 0.02 we get the same result as with 100 kg at the ends – N.J Abel May 07 '17 at 11:24
  • @N.JAbel The amplitude ? Did you notice that I changed `ylim` parameters ? – keepAlive May 07 '17 at 12:40
  • Thank you. I mean the |amplitude| of anti nodes being 10. Isn't that too large compared to our initial sin wave with an amplitude of 1. – N.J Abel May 07 '17 at 13:27
  • @N.JAbel. Is this process considered **theoretically** as continuous or discrete ? – keepAlive May 07 '17 at 15:45
  • I think it is discrete – N.J Abel May 07 '17 at 17:25
  • @N.JAbel. I asked you this because I think the problem comes from the discretization of the 1d space, which is configured via the number of mass points. Set it be 10 times less dense, i.e. `N = 20`, and you will see what I mean. – keepAlive May 07 '17 at 17:27
  • N = 20 is too few points to represent a wave. I tried with using j == 0 – N.J Abel May 07 '17 at 18:48