import sympy as sp
# シンボルの定義
x, y, u = sp.symbols('x y u')
p = sp.Derivative(u, x)
q = sp.Derivative(u, y)
# 例としてラグランジュの方程式を設定
F = sp.Function('F')(x, y, u, p, q)
F_eq = sp.Eq(F, 0)
# 完全微分方程式の例
M, N = sp.symbols('M N')
M = sp.Function('M')(x, y)
N = sp.Function('N')(x, y)
dx = sp.Symbol('dx')
dy = sp.Symbol('dy')
# 完全微分方程式を解く
eq = M*dx + N*dy
is_complete = sp.diff(M, y) == sp.diff(N, x)
if is_complete:
print("完全微分方程式です。")
else:
print("完全微分方程式ではありません。")
# シャルピーの解法に関するサンプル
# 非線形1階偏微分方程式の例
u = sp.Function('u')(x, y)
F = sp.Eq(u.diff(x)**2 + u.diff(y)**2 - 1, 0)
solution = sp.dsolve(F, u)
print(solution)
import numpy as np
import matplotlib.pyplot as plt
# Function to define the characteristic curves
def characteristic_curves(x, y):
X, Y = np.meshgrid(x, y)
U = np.sin(X) * np.cos(Y)
return X, Y, U
x = np.linspace(0, 2 * np.pi, 100)
y = np.linspace(0, 2 * np.pi, 100)
X, Y, U = characteristic_curves(x, y)
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, U, levels=50, cmap='viridis')
plt.colorbar(label='u(x, y)')
plt.title('Characteristic Curves of a First-Order PDE')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# Function to define the solution for Lagrange's PDE
def lagrange_characteristic_curves(x, y):
X, Y = np.meshgrid(x, y)
U = np.exp(-X**2 - Y**2)
return X, Y, U
X, Y, U = lagrange_characteristic_curves(x, y)
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, U, levels=50, cmap='plasma')
plt.colorbar(label='u(x, y)')
plt.title('Characteristic Curves of Lagrange\'s PDE')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Define constant coefficients
A = 1
B = 0.5
C = 1
D = 0.1
E = 0.2
F = 0.3
G = 1 # Non-homogeneous term
# Compute the discriminant
def discriminant(x, y):
return B**2 - A * C
# Compute the non-homogeneous term
def non_homogeneous_term(x, y):
return G * np.ones_like(x) # Simple constant function for illustration
# Create a grid
x = np.linspace(0, 2 * np.pi, 100)
y = np.linspace(0, 2 * np.pi, 100)
X, Y = np.meshgrid(x, y)
# Compute discriminant values
D = discriminant(X, Y)
# Compute non-homogeneous term values
G_term = non_homogeneous_term(X, Y)
# Plotting different types of PDE classification
plt.figure(figsize=(12, 6))
# Elliptic, Parabolic, and Hyperbolic regions
plt.subplot(1, 2, 1)
contour_elliptic = plt.contourf(X, Y, D, levels=[-np.inf, 0], colors='lightblue', alpha=0.6)
contour_hyperbolic = plt.contourf(X, Y, D, levels=[0, np.inf], colors='lightcoral', alpha=0.6)
plt.contour(X, Y, D, levels=[0], colors='black', linestyles='--')
plt.title('PDE Classification')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar(contour_elliptic, label='Elliptic')
plt.colorbar(contour_hyperbolic, label='Hyperbolic')
# Plotting non-homogeneous term
plt.subplot(1, 2, 2)
plt.contourf(X, Y, G_term, levels=50, cmap='plasma')
plt.colorbar(label='G(x, y)')
plt.title('Non-Homogeneous Term')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
from scipy.fft import fft, fftfreq
# Function to compute the Fourier transform
def fourier_transform(f, x):
N = len(x)
T = x[1] - x[0] # Sampling interval
yf = fft(f)
xf = fftfreq(N, T)[:N // 2]
return xf, np.abs(yf[:N // 2])
# Create x values and function f(x)
x = np.linspace(-10, 10, 1000)
f = np.sin(x) # Example function
# Compute Fourier transform
xf, yf = fourier_transform(f, x)
plt.figure(figsize=(8, 6))
plt.plot(xf, yf, label='Fourier Transform')
plt.title('Fourier Transform')
plt.xlabel('Frequency (k)')
plt.ylabel('Magnitude')
plt.legend()
plt.grid(True)
plt.show()
from scipy.special import legendre
import numpy as np
import matplotlib.pyplot as plt
# Define x values
x = np.linspace(-1, 1, 500)
# Compute Legendre polynomials
P0 = legendre(0)(x)
P1 = legendre(1)(x)
P2 = legendre(2)(x)
P3 = legendre(3)(x)
plt.figure(figsize=(8, 6))
plt.plot(x, P0, label='P0(x)')
plt.plot(x, P1, label='P1(x)')
plt.plot(x, P2, label='P2(x)')
plt.plot(x, P3, label='P3(x)')
plt.title('Legendre Polynomials')
plt.xlabel('x')
plt.ylabel('P_n(x)')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# Parameters
L = 10.0 # Length of the domain
Nx = 100 # Number of spatial points
Nt = 200 # Number of time steps
alpha = 0.01 # Diffusion coefficient
# Spatial and temporal grids
x = np.linspace(0, L, Nx)
dx = x[1] - x[0]
dt = 0.1 # Time step size
t_span = (0, dt * Nt)
# Initial condition
u0 = np.sin(np.pi * x)
def heat_eq(t, u_flat):
u = u_flat.reshape((Nx,))
u_xx = np.roll(u, -1) - 2*u + np.roll(u, 1)
u_xx[0] = u_xx[-1] = 0 # Boundary conditions (u_x = 0)
du_dt = alpha * u_xx / dx**2
return du_dt.flatten()
# Solve PDE
u0_flat = u0.flatten()
sol = solve_ivp(heat_eq, t_span, u0_flat, t_eval=np.linspace(0, t_span[1], Nt))
# Reshape solution
u_sol = sol.y.reshape((Nx, Nt))
# Plot results
plt.figure(figsize=(10, 6))
for i in range(0, Nt, Nt//10):
plt.plot(x, u_sol[:, i], label=f't = {i*dt:.2f}')
plt.title('Heat Equation Solution Over Time')
plt.xlabel('x')
plt.ylabel('u(x,t)')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Define parameters
pi = np.pi
x = np.linspace(0, 1, 100) # Spatial domain
t = np.linspace(0, 0.1, 50) # Temporal domain
# Create meshgrid for X and T
X, T = np.meshgrid(x, t)
# Compute the analytical solution
U = np.exp(-pi**2 * T) * np.sin(pi * X)
# Plot the solution
fig = plt.figure(figsize=(12, 8))
# Surface plot of the analytical solution
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.plot_surface(X, T, U, cmap='viridis')
ax1.set_title('Analytical Solution of Heat Conduction Equation')
ax1.set_xlabel('Distance x')
ax1.set_ylabel('Time t')
ax1.set_zlabel('Solution u(x,t)')
# Profile plot at different times
ax2 = fig.add_subplot(1, 2, 2)
for i in range(0, len(t), len(t) // 5):
ax2.plot(x, U[i, :], label=f't={t[i]:.2f}')
ax2.set_xlabel('Distance x')
ax2.set_ylabel('Solution u(x,t)')
ax2.set_title('Solution Profiles at Different Times')
ax2.legend()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Define parameters
xMin = 0
xMax = 1
tMin = 0
tMax = 0.1
nX = 10
nT = 40
deltaT = (tMax - tMin) / (nT - 1)
deltaX = (xMax - xMin) / (nX - 1)
r = deltaT / (deltaX**2)
# Create spatial and temporal grids
x = np.linspace(xMin, xMax, nX)
t = np.linspace(tMin, tMax, nT)
# Initialize the solution matrix
sol = np.zeros((nX, nT))
# Set initial condition
sol[:, 0] = np.sin(np.pi * x)
# Construct the matrix for the finite difference method
pdMat = np.eye(nX-2) - 2*np.eye(nX-2) + np.eye(nX-2)
# Time-stepping loop
for it in range(1, nT):
sol[1:-1, it] = sol[1:-1, it-1] + r * np.dot(pdMat, sol[1:-1, it-1])
# Plot the results
X, T = np.meshgrid(x, t)
fig = plt.figure(figsize=(12, 8))
# 3D surface plot
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.plot_surface(X, T, sol.T, cmap='viridis')
ax1.set_title('Numerical Solution of Heat Conduction Equation')
ax1.set_xlabel('Distance x')
ax1.set_ylabel('Time t')
ax1.set_zlabel('Solution u(x,t)')
# 2D profile plots at different times
ax2 = fig.add_subplot(1, 2, 2)
for i in range(0, nT, nT // 5):
ax2.plot(x, sol[:, i], label=f't={t[i]:.2f}')
ax2.set_xlabel('Distance x')
ax2.set_ylabel('Solution u(x,t)')
ax2.set_title('Solution Profiles at Different Times')
ax2.legend()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Define parameters
xMin = 0
xMax = 1
tMin = 0
tMax = 0.1
nX = 10
nT = 40
deltaT = (tMax - tMin) / (nT - 1)
deltaX = (xMax - xMin) / (nX - 1)
r = deltaT / (deltaX**2)
# Create spatial and temporal grids
x = np.linspace(xMin, xMax, nX)
t = np.linspace(tMin, tMax, nT)
# Initialize the solution matrix
sol = np.zeros((nX, nT))
# Set initial condition
sol[:, 0] = np.sin(np.pi * x)
# Construct the matrix for the implicit method (Crank-Nicolson scheme)
main_diag = (1 + 2 * r) * np.ones(nX - 2)
off_diag = -r * np.ones(nX - 3)
pdMat = np.diag(main_diag) + np.diag(off_diag, -1) + np.diag(off_diag, 1)
# Time-stepping loop
for it in range(1, nT):
# Solve the system of linear equations using np.linalg.solve
sol[1:-1, it] = np.linalg.solve(pdMat, sol[1:-1, it-1])
# Plot the results
X, T = np.meshgrid(x, t)
fig = plt.figure(figsize=(12, 8))
# 3D surface plot
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.plot_surface(X, T, sol.T, cmap='viridis')
ax1.set_title('Numerical Solution of Heat Conduction Equation (Implicit Method)')
ax1.set_xlabel('Distance x')
ax1.set_ylabel('Time t')
ax1.set_zlabel('Solution u(x,t)')
# 2D profile plots at different times
ax2 = fig.add_subplot(1, 2, 2)
for i in range(0, nT, nT // 5):
ax2.plot(x, sol[:, i], label=f't={t[i]:.2f}')
ax2.set_xlabel('Distance x')
ax2.set_ylabel('Solution u(x,t)')
ax2.set_title('Solution Profiles at Different Times')
ax2.legend()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters
length = 1.0 # Length of the rod
n_points = 100 # Number of spatial points
x = np.linspace(0, length, n_points) # Spatial grid
alpha = 1.0 # Thermal diffusivity (lambda/c)
time_steps = np.linspace(0, 0.1, 10) # Time steps for sweeping
# Initial condition
u_initial = np.sin(np.pi * x) # Example initial condition
# Plotting
plt.figure(figsize=(10, 6))
for time in time_steps:
# Heat equation solution (analytical solution)
u = np.exp(-alpha * np.pi**2 * time) * np.sin(np.pi * x)
# Plot each time step
plt.plot(x, u, label=f't={time:.2f}')
plt.xlabel('Position x')
plt.ylabel('Temperature u(x, t)')
plt.title('Heat Equation Solution Over Time')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# Parameters
length = 1.0 # Length of the string
n_points = 100 # Number of spatial points
x = np.linspace(0, length, n_points) # Spatial grid
c = 1.0 # Wave speed
dx = x[1] - x[0] # Spatial step size
dt = 0.01 # Time step size
n_steps = 200 # Number of time steps
# Initialize variables
u = np.zeros(n_points) # Current displacement
u_prev = np.zeros(n_points) # Displacement at previous time step
u_next = np.zeros(n_points) # Displacement at next time step
# Initial condition: Gaussian-shaped pulse
u = np.exp(-100 * (x - length/2)**2)
# Plot setup
fig, ax = plt.subplots(figsize=(8, 6))
line, = ax.plot(x, u, label='Wave')
ax.set_xlabel('Position x')
ax.set_ylabel('Displacement u(x, t)')
ax.set_title('Wave Equation Solution')
ax.set_ylim(-1.5, 1.5)
ax.legend()
# Animation function
def animate(t):
global u, u_prev, u_next
u_next[1:-1] = 2 * u[1:-1] - u_prev[1:-1] + (c**2 * dt**2 / dx**2) * (u[2:] - 2 * u[1:-1] + u[:-2])
u_prev[:] = u
u[:] = u_next
line.set_ydata(u)
return line,
# Create animation
ani = FuncAnimation(fig, animate, frames=n_steps, interval=dt*1000, blit=True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters
L = 1.0 # Length of the domain
n_points = 100 # Number of spatial points
x = np.linspace(0, L, n_points) # Spatial grid
# Pe numbers to evaluate
Pe_values = [0.1, 1, 10] # Example Pe numbers
# Function to compute solution for given Pe
def solve_convection_diffusion(Pe):
# General solution coefficients
A = -1 / (1 - np.exp(-Pe))
B = 1 / (1 - np.exp(-Pe))
# Solution
u = A + B * np.exp(Pe * x)
return u
# Plotting
plt.figure(figsize=(10, 6))
for Pe in Pe_values:
u = solve_convection_diffusion(Pe)
plt.plot(x, u, label=f'Pe={Pe}')
plt.xlabel('Position x')
plt.ylabel('Temperature u(x)')
plt.title('Steady-State Convection-Diffusion Equation')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters
length = 1.0 # Length of the domain
n_points = 50 # Number of spatial points
n_time_steps = 100 # Number of time steps
alpha = 1.0 # Diffusion coefficient (lambda)
dx = length / (n_points - 1) # Spatial step size
dt = 0.01 # Time step size
lambda_ = alpha * dt / dx**2 # Stability parameter
# Create spatial and time grids
x = np.linspace(0, length, n_points)
t = np.linspace(0, dt * (n_time_steps - 1), n_time_steps)
# Initialize solution matrix
u = np.zeros((n_time_steps, n_points))
# Set initial condition
u[0, :] = 0
# Set boundary conditions
u[:, 0] = 1
u[:, -1] = 1
# FTCS scheme
for n in range(0, n_time_steps - 1):
for i in range(1, n_points - 1):
u[n+1, i] = u[n, i] + lambda_ * (u[n, i+1] - 2 * u[n, i] + u[n, i-1])
# Plotting
plt.figure(figsize=(10, 6))
for n in range(0, n_time_steps, int(n_time_steps / 10)):
plt.plot(x, u[n, :], label=f't={t[n]:.2f}')
plt.xlabel('Position x')
plt.ylabel('Temperature u(x, t)')
plt.title('Unsteady Diffusion Equation Solution Over Time')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Parameters
length = 1.0 # Length of the domain
n_points = 100 # Number of spatial points
n_time_steps = 100 # Number of time steps
c = 1.0 # Wave speed
dx = length / (n_points - 1) # Spatial step size
dt = 0.01 # Time step size
alpha = (c * dt / dx)**2 # Stability parameter
# Create spatial and time grids
x = np.linspace(0, length, n_points)
t = np.linspace(0, dt * (n_time_steps - 1), n_time_steps)
# Initialize solution matrix
u = np.zeros((n_time_steps, n_points))
# Set initial condition
u[0, :] = np.sin(np.pi * x)
u[1, :] = u[0, :] # Assuming initial velocity is zero
# Time-stepping using the finite difference method
for n in range(1, n_time_steps - 1):
for i in range(1, n_points - 1):
u[n + 1, i] = 2 * (1 - alpha) * u[n, i] - u[n - 1, i] + alpha * (u[n, i + 1] + u[n, i - 1])
# Create an animation
fig, ax = plt.subplots()
line, = ax.plot(x, u[0, :], lw=2)
def animate(i):
line.set_ydata(u[i, :])
return line,
ax.set_xlim(0, length)
ax.set_ylim(-1, 1)
ax.set_xlabel('x')
ax.set_ylabel('u(x, t)')
ax.set_title('Wave Equation')
ani = animation.FuncAnimation(fig, animate, frames=n_time_steps, interval=50, blit=True)
# Save animation as a gif
ani.save('wave_animation.gif', writer='imagemagick')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Parameters
Lx = 2 * np.pi # Length of the domain
N = 256 # Number of spatial points
dx = Lx / N # Spatial step size
dt = 1.0e-2 # Time step size
nu = 2.0e-3 # Diffusion coefficient
a = 10.0 # Parameter in the initial condition
t_max = 100.0 # Maximum time
t_out = 1.0 # Output interval
# Create spatial grid
x = np.linspace(0, Lx, N, endpoint=False)
x = np.concatenate(([x[-1] + dx], x, [x[0] - dx])) # Periodic boundary conditions
u = np.zeros(N + 2) # Solution array
# Initial condition
u[1:-1] = np.exp(-a * (x[1:-1] - Lx / 2.0)**2)
# Initialize variables
t = 0.0
t_step = int(t_max / dt)
j_out = int(t_out / dt)
# Create an array to store the solution for animation
u_history = []
# Time-stepping loop
for n in range(t_step):
du = np.zeros(N + 2)
# Update internal points
du[1:-1] = nu * (u[2:] - 2 * u[1:-1] + u[:-2]) / dx**2
# Update solution
u[1:-1] += du[1:-1] * dt
# Apply periodic boundary conditions
u[0] = u[-2]
u[-1] = u[1]
# Store the result for animation
if n % j_out == 0:
u_history.append(u.copy())
t += dt
# Create an animation
fig, ax = plt.subplots()
line, = ax.plot(x, u_history[0], lw=2)
def animate(i):
line.set_ydata(u_history[i])
return line,
ax.set_xlim(0, Lx)
ax.set_ylim(np.min(u_history), np.max(u_history))
ax.set_xlabel('x')
ax.set_ylabel('u(x, t)')
ax.set_title('Diffusion Equation')
ani = animation.FuncAnimation(fig, animate, frames=len(u_history), interval=50, blit=True)
# Save animation as a gif
ani.save('diffusion_animation.gif', writer='imagemagick')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Parameters
Lx = 2 * np.pi # Length of the domain
N = 256 # Number of spatial points
dx = Lx / N # Spatial step size
dt = 1.0e-2 # Time step size
nu = 2.0e-3 # Diffusion coefficient
a = 10.0 # Parameter in the initial condition
t_max = 100.0 # Maximum time
t_out = 1.0 # Output interval
# Create spatial grid
x = np.linspace(0, Lx, N, endpoint=False)
x = np.concatenate(([x[-1] + dx], x, [x[0] - dx])) # Periodic boundary conditions
u = np.zeros(N + 2) # Solution array
# Initial condition
u[1:-1] = np.exp(-a * (x[1:-1] - Lx / 2.0)**2)
# Initialize variables
t = 0.0
t_step = int(t_max / dt)
j_out = int(t_out / dt)
# Array to store the numerical and analytical solutions for animation
u_history = []
tu_history = []
# Analytical solution calculation function
def analytical_solution(x, t, a, nu):
t0 = 1.0 / (4.0 * nu * a)
return np.sqrt(a / (4.0 * np.pi * nu * (t + t0))) * np.exp(-((x - Lx / 2.0)**2) / (4.0 * nu * (t + t0)))
# Time-stepping loop
for n in range(t_step):
du = np.zeros(N + 2)
# Update internal points
du[1:-1] = nu * (u[2:] - 2 * u[1:-1] + u[:-2]) / dx**2
# Update solution
u[1:-1] += du[1:-1] * dt
# Apply periodic boundary conditions
u[0] = u[-2]
u[-1] = u[1]
# Store the result for animation
if n % j_out == 0:
u_history.append(u.copy())
tu_history.append(analytical_solution(x, t, a, nu))
t += dt
# Create an animation
fig, ax = plt.subplots()
line_num, = ax.plot(x, u_history[0], lw=2, label='Numerical')
line_theory, = ax.plot(x, tu_history[0], lw=2, linestyle='--', label='Analytical')
def animate(i):
line_num.set_ydata(u_history[i])
line_theory.set_ydata(tu_history[i])
return line_num, line_theory
ax.set_xlim(0, Lx)
ax.set_ylim(np.min(u_history), np.max(u_history))
ax.set_xlabel('x')
ax.set_ylabel('u(x, t)')
ax.set_title('Diffusion Equation')
ax.legend()
ani = animation.FuncAnimation(fig, animate, frames=len(u_history), interval=50, blit=True)
# Save animation as a gif
ani.save('diffusion_comparison_animation.gif', writer='imagemagick')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Parameters
Lx = 2 * np.pi # Length of the domain
N = 256 # Number of spatial points
dx = Lx / N # Spatial step size
dt = 1.0e-2 # Time step size
c = 0.1 # Advection speed
t_max = 100.0 # Maximum time
t_out = 1.0 # Output interval
# Create spatial grid
x = np.linspace(0, Lx, N, endpoint=False)
x = np.concatenate(([x[-1] + dx], x, [x[0] - dx])) # Periodic boundary conditions
u = np.zeros(N + 2) # Solution array
# Initial condition
u[1:-1] = np.sin(4 * np.pi * x[1:-1] / Lx)
# Initialize variables
t = 0.0
t_step = int(t_max / dt)
j_out = int(t_out / dt)
# Array to store the numerical solution for animation
u_history = []
# Time-stepping loop
for n in range(t_step):
du = np.zeros(N + 2)
# Central difference scheme for advection
du[1:-1] = (u[2:] - u[:-2]) * c * dt / (2.0 * dx)
# Update solution
u[1:-1] -= du[1:-1]
# Apply periodic boundary conditions
u[0] = u[-2]
u[-1] = u[1]
# Store the result for animation
if n % j_out == 0:
u_history.append(u.copy())
t += dt
# Create an animation
fig, ax = plt.subplots()
line, = ax.plot(x, u_history[0], lw=2)
def animate(i):
line.set_ydata(u_history[i])
return line,
ax.set_xlim(0, Lx)
ax.set_ylim(np.min(u_history), np.max(u_history))
ax.set_xlabel('x')
ax.set_ylabel('u(x, t)')
ax.set_title('Advection Equation')
ani = animation.FuncAnimation(fig, animate, frames=len(u_history), interval=50, blit=True)
# Save animation as a gif
ani.save('advection_animation.gif', writer='imagemagick')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters
L = 10.0 # Length of the domain
T = 1.0 # Total time
nx = 50 # Number of spatial points
nt = 100 # Number of time steps
alpha = 0.01 # Diffusion coefficient
dx = L / (nx - 1) # Spatial grid size
dt = T / nt # Time step size
# Stability condition check
if alpha * dt / dx**2 > 0.5:
raise ValueError("Stability condition violated")
# Initial condition
u = np.zeros(nx)
u[int(nx/2)] = 1 / dx
# Initialize matrices for Crank-Nicolson
r = alpha * dt / (2 * dx**2)
A = np.diag((1 + r) * np.ones(nx - 2)) + np.diag(-r * np.ones(nx - 3), 1) + np.diag(-r * np.ones(nx - 3), -1)
B = np.diag((1 - r) * np.ones(nx - 2)) + np.diag(r * np.ones(nx - 3), 1) + np.diag(r * np.ones(nx - 3), -1)
# Time-stepping loop
for n in range(nt):
b = np.dot(B, u[1:-1])
u[1:-1] = np.linalg.solve(A, b)
# Plotting results
x = np.linspace(0, L, nx)
plt.plot(x, u, label='t = {:.2f}'.format(T))
plt.xlabel('x')
plt.ylabel('u')
plt.title('Diffusion Equation - Crank-Nicolson Method')
plt.legend()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters
L = 10.0 # Length of the string
T = 2.0 # Total time
nx = 100 # Number of spatial points
nt = 200 # Number of time steps
c = 1.0 # Wave speed
dx = L / (nx - 1) # Spatial grid size
dt = T / nt # Time step size
# Stability condition check
if c * dt / dx > 1:
raise ValueError("Stability condition violated")
# Initial condition
x = np.linspace(0, L, nx)
u = np.zeros(nx)
u[int(nx/4):int(3*nx/4)] = 1 # Initial displacement
# Initialize arrays for old and new time steps
u_prev = np.copy(u)
u_next = np.zeros(nx)
# Time-stepping loop
for n in range(nt):
# Update the next time step using the finite difference method
u_next[1:-1] = (2 * u[1:-1] - u_prev[1:-1] +
(c * dt / dx)**2 * (u[2:] - 2 * u[1:-1] + u[:-2]))
# Update arrays
u_prev[:] = u
u[:] = u_next
# Plot at specific intervals
if n % 20 == 0:
plt.plot(x, u, label=f't = {n * dt:.2f}')
# Plot results
plt.xlabel('x')
plt.ylabel('u')
plt.title('String Vibration - Finite Difference Method')
plt.legend()
plt.show()