In [89]:
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,10,50)
y = np.array(x)
plt.scatter(x, y,color="grey")
plt.show()
Good examples for reference:
Steps:
The lefthand side of the equation represents the second order derivative of space, and the righthand side the second order derivative of time. Using a computational cluster, we can visualize what's going on.
[insert drawing]
Page 96 provides an example of a 1D diffuusion equation. The discretization is as follows for our equation:
$$\frac{u_j^{n+1} - 2 u^n_j + u_j^{n-1}} { \Delta t^2} = v^2 \left( \frac{u^n_{j+1} - 2u^n_j + u^n_{j-1}}{ \Delta x^2} \right) $$Next, we want to know what $u_j^{n+1}$ is. Therefore, we start with removing the fraction on the lefthand side by bringing $(\Delta t)^2$ to the righthand side.
$$u_j^{n+1} - 2 u^n_j + u_j^{n-1} = \frac{ v^2 \Delta t^2}{\Delta x^2 } \left( u^n_{j+1} - 2u^n_j + u^n_{j-1} \right) $$And bring it to the left side, creates an update rule:
$$u_j^{n+1} = 2 u^n_j - u_j^{n-1} + \frac{ v^2 \Delta t^2}{\Delta x^2 } \left( u^n_{j+1} - 2u^n_j + u^n_{j-1} \right) $$This is a second-order approximation in time and space.
The Von Neumann stability analysis helps to find out how big of a step size can you use. Equation (9.21) on page 96 of the textbook provides an example. The Van Neuman for this assignment is as follows.
The first step is to create an ansatz:
$$u_j^n = A^n e^{ijk \Delta x} $$And plug it in the equation:
$$A^{n+1}e^{ijk \Delta x} =2 A^n e^{ijk \Delta x} - A^{n-1} e^{ijk \Delta x} + \frac{ v^2 \Delta t^2}{\Delta x^2 } A^n \left( e^{i(j+1)k\Delta x} - 2 e^{ijk \Delta x} + e^{i(j-1)k \Delta x } \right) $$For more readability, bring some elements to the left side and isolate A.
$$\left( A^{n+1} - 2 A^n + A^{n-1} \right) e^{ijk \Delta x} = \frac{ v^2 \Delta t^2}{\Delta x^2 } A^n \left( e^{i(j+1)k\Delta x} - 2 e^{ijk \Delta x} + e^{i(j-1)k \Delta x } \right) $$Divide by $A^n$ from righthand side.
$$\left( A - 2 + \frac{1}{A} \right) e^{ijk \Delta x} = \frac{ v^2 \Delta t^2}{\Delta x^2 } \left( e^{i(j+1)k\Delta x} - 2 e^{ijk \Delta x} + e^{i(j-1)k \Delta x } \right) $$Divide by $e^{i(j+1)k\Delta x}$ from lefthand side.
$$A - 2 + \frac{1}{A} = \frac{ v^2 \Delta t^2}{\Delta x^2 } \left( e^{ik\Delta x} - 2 + e^{-ik \Delta x } \right) $$Next:
$$ = \frac{2 v^2 \Delta t^2}{\Delta x^2 } \left( \cos(k \Delta x) -1 \right) $$Apply object-oriented algebra by writing the equation as $\alpha$.
The product of the two solutions are in the positive and negative.
$$A_+ A_- = \left( -\alpha + \sqrt{\alpha ^2 -1} \right) \left( -\alpha - \sqrt{\alpha ^2 -1} \right) $$Cross-multiplication results in:
$$A_+ A_- = \alpha ^2 - (\alpha ^2 - 1) $$With a solution of:
$$A_+ A_- = 1 $$No modes decay or grow, because the amplitude of both solutions is 1.
Advice from Amira: Perhaps some coding here. Just some points and plot it.Eyal did this in matplotlib: http://fab.cba.mit.edu/classes/864.20/people/eyal/assets/html/week5_diffusion.html
Have the update rule ready. This is what you're going to code in numpy.
First, draw a line of dots.
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,10,50)
y = np.array(x)
plt.scatter(x, y,color="grey")
plt.show()
Next, add the constants and equations.
The difference between np.arange and np.array is that arange returns an array. np.arange(3) = array([0, 1, 2])
# Adapted from: https://fab.cba.mit.edu/classes/864.20/people/allan/week5/index.html
# Numerically solve the wave equation
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
fig = plt.figure()
ax = fig.add_subplot(111)
fig.show()
fig.canvas.draw()
# Constants
v = 1
dt = 0.1
dx = 1 # given by Neil
N = 500 # number of points given by Neil
H = 10
T = 10000
x = np.arange(0, N) # Create an array
past = np.zeros(N) # u^n-1 of j (creates an array of zeroes)
present = np.zeros(N) #u^n of j
future = np.zeros(N) #u^n+1 of j
present[N//2] = H - 2 # double division = integer division,
α = (v * dt / dx) ** 2 # plug in equation (v delta t / )
for t in range(T):
future = 2 * present - past # 2x present - past (from equation unjn+1)
for j in range(1, N-1): # go through all of the points with j, start with 1 and 1 before last because otherwise you don't have neighbors.
future[j] += α * (present[j + 1] - 2 * present[j] + present[j-1])
past, present = present, future
if t % 20 == 0: # 20 timesteps
ax.clear()
ax.plot(x, present)
fig.canvas.draw()
/var/folders/zd/7vg9k68x0fg0qk2mmd4t0zyc0000gn/T/ipykernel_868/2589672503.py:10: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. fig.show()
matrix
First, let's start with creating a 3D plot. This Python handbook offers guidance on how to do 3D plots. More explanation on 3D plots can be found in Python Numerical Models.
import numpy as np
import matplotlib.pyplot as plt
ax = plt.axes(projection='3d')
x = np.linspace(0, 4*np.pi, 1000) # This creates a line at (start, stop, number of samples to generate)
y = np.sin(x) # Adding a sinus function of X
z = np.linspace(1,2,1) #Ading z axis
ax.plot(x, y, z)
ax.set_ylabel("Temperature")
ax.set_xlabel("Time")
ax.set_zlabel("Density", rotation=60)
plt.title("Experimental 3D plot", fontsize = 14)
plt.show()
Meaning of +=
i=0
i+=1
i=i+1