62 lines
1.8 KiB
Python
62 lines
1.8 KiB
Python
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import time, sys
|
|
|
|
"""
|
|
Simulation comments :
|
|
As we increase the number of grid points, the solution tends to a hat that
|
|
have move to the right.
|
|
For nx<50, we have a bell not fully developed.
|
|
For nx=50, we have something ike a perfect bell.
|
|
For nx>50, we go back to the hat function.
|
|
|
|
For big dt, ~0.25,we have divergence of data.
|
|
Neverthe less, if we increase the timestep x2, we obtain a better solution
|
|
even if the number of grid points stays lower.
|
|
|
|
There seems to be some numerical diffusion problem, for high nmber of grid points
|
|
the solution seems more like a transitory phenomena that diverges greatly !
|
|
"""
|
|
|
|
nx = 80 # Number of grid points.
|
|
dx = 2 / (nx-1) # Distance between any pair of adjacent grid points.
|
|
nt = 25 # Number of timesteps we want to calculate.
|
|
dt = 0.025 # Amount of time each timesteps covers (delta t)
|
|
c= 1 # Assume wavespeed of c = 1.
|
|
|
|
# Boundary conditions.
|
|
# Starting all values as 1 m/s.
|
|
u = np.ones(nx)
|
|
# Setting u = 2 between 0.5 and 1 as per out initial conditions.
|
|
u[int(0.5 / dx):int(1 / dx + 1)] = 2
|
|
|
|
# Need still to write it down to discretize the space.
|
|
# Calculating the analytical solution.
|
|
u_theo = u.copy()
|
|
for i in range(1,nx+1)
|
|
u_theo[i+1] = u_theo[i] - c*(nt*dt)
|
|
|
|
|
|
# PLot initial conditions.
|
|
fig, ax = plt.subplots(nrows=1, ncols =3, sharey=True)
|
|
|
|
ax[0].plot(np.linspace(0,2,nx),u)
|
|
|
|
# Temporary array for the solutions.
|
|
un = np.ones(nx)
|
|
|
|
for n in range(nt): # Loop for values of n from 0 to nt.
|
|
un = u.copy() # Copy the existing values of u into un.
|
|
|
|
for i in range(1, nx): # Loop for values in u from 1 to nx.
|
|
u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
|
|
|
|
# PLot the solution at t = t_final.
|
|
ax[1].plot(np.linspace(0,2,nx), u)
|
|
|
|
# Analytical solution
|
|
ax[2].plot(np.linspace(0,2,nx),u_theo)
|
|
|
|
plt.show()
|
|
|