!conda install -y matplotlib numpy

import matplotlib.pyplot as plt
import numpy as np

# Set the random seed to keep results reproducible
np.random.seed(0)

k = 1000

# Initialize the starting point of walks
start = np.zeros(k)

# Generate random steps for each walk
steps = np.random.choice([-1, 1], size=(1000,k))

# Broadcast the start position and add up the steps for each walk
walks = np.broadcast_to(start, (1000, k)) + np.cumsum(steps, axis=0)

# Plot each walk
for i in range(k):
    plt.plot(walks[:, i])

plt.title('%s Random Walks'%k)
plt.xlabel('Steps')
plt.ylabel('Position')
plt.sow()

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.linalg import solve

# Constants
alpha = 0.01  # Diffusion coefficient
L = 1         # Length of domain
T = 3         # Time to solve until
Nx = 50       # Number of spatial points in grid
Nt = 300      # Number of time steps
dx = L / (Nx - 1)  # Spatial discretization size
dt = T / (Nt - 1)  # Time discretization size

# Ensure the method is stable
if alpha*dt/dx**2 > 0.5:
    raise ValueError("Solution will be unstable")

# Initialize grid
grid = np.zeros(shape=(Nt, Nx))

# Initial condition
grid[0, :] = np.sin(np.pi*np.linspace(0, L, Nx))
grid[0, -1] = 0

# Finite differences matrix
k = alpha * dt / dx**2
elements = np.array([-k, 1+2*k, -k])
diagonals = [-1, 0, 1]
M = diags(elements, diagonals, shape=(Nx-2, Nx-2)).toarray()

# Time-stepping loop
for t in range(Nt - 1):
    b = grid[t, 1:-1].copy()
    b[0] += k*grid[t, 0]
    b[-1] += k*grid[t, -1]
    grid[t+1, 1:-1] = solve(M, b)

# Plot solution
for t in range(Nt):
    if t % 30 == 0:  # Only plot every 30th time step
        plt.plot(np.linspace(0, L, Nx), grid[t, :], lw=2)

plt.xlabel('x')
plt.ylabel('u')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Set the parameters
c = 1.0  # wave speed
L = 10.0  # size of domain
T = 2.0  # time interval
Nx = 100  # number of spatial steps
Nt = 100  # number of time steps

# Set up the grid
x = np.linspace(0, L, Nx+1)  # space
t = np.linspace(0, T, Nt+1)  # time

# Initialize the function
u = np.zeros((Nt+1, Nx+1))

# Set up the initial condition: Gaussian pulse in the center
u[0, :] = np.exp(-0.5*((x-L/2)**2))

# Solve the 2D wave equation with finite differences
for n in range(0, Nt):
    for i in range(1, Nx):
        u[n+1, i] = 2*u[n,i] - u[n-1,i] + c**2*(t[n+1]-t[n])**2/(x[i+1]-x[i])**2 *\
                    (u[n,i+1] - 2*u[n,i] + u[n,i-1])  
    # Enforce boundary conditions
    u[n+1, 0] = u[n+1, Nx] = 0

# Plotting
plt.imshow(u, extent=(0,L,T,0), aspect='auto')
plt.colorbar(label='Amplitude')
plt.xlabel('Position')
plt.ylabel('Time')
plt.title('2D-wave equation')
plt.sow()