偏微分方程式

```
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()


```
http://dsl4.eee.u-ryukyu.ac.jp/DOCS/NumA/p15.pdf

```
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()


```
http://nkl.cc.u-tokyo.ac.jp/14n/PDE.pdf

```

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()


```

https://ocw.kyoto-u.ac.jp/wp-content/uploads/2000/04/1999_introduction-to-simulation_08.pdf
```

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()

```


```


```


```


```


```

```


```


```


```


```