微分方程式

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 方程式の定義
def model(y, x):
    return x * y

# 初期条件
y0 = 1
x = np.linspace(0, 5, 100)

# 微分方程式の解
solution = odeint(model, y0, x)

# プロット
plt.plot(x, solution)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of dy/dx = x * y')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 方程式の定義
def model(y, x):
    return y / x

# 初期条件
y0 = 1
x = np.linspace(0.1, 5, 100)  # xが0に近づくと問題が発生するため、0.1から開始

# 微分方程式の解
solution = odeint(model, y0, x)

# プロット
plt.plot(x, solution)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of dy/dx = y / x')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 方程式の定義
def model(y, x):
    p = 1
    q = np.sin(x)
    return q - p * y

# 初期条件
y0 = 0
x = np.linspace(0, 10, 100)

# 微分方程式の解
solution = odeint(model, y0, x)

# プロット
plt.plot(x, solution)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of dy/dx + y = sin(x)')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 方程式の定義
def model(Y, x):
    y, dydx = Y
    return [dydx, -2 * dydx - y]

# 初期条件
y0 = [1, 0]  # 初期値 [y(0), dy/dx(0)]
x = np.linspace(0, 10, 100)

# 微分方程式の解
solution = odeint(model, y0, x)

# プロット
plt.plot(x, solution[:, 0])
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of y\'\' - 2y\' - y = 0')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 方程式の定義
def model(Y, x):
    y, dydx = Y
    g = np.sin(x)
    return [dydx, -2 * dydx - y + g]

# 初期条件
y0 = [0, 0]  # 初期値 [y(0), dy/dx(0)]
x = np.linspace(0, 10, 100)

# 微分方程式の解
solution = odeint(model, y0, x)

# プロット
plt.plot(x, solution[:, 0])
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of y\'\' - 2y\' - y = sin(x)')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 方程式の定義
def model(Y, x):
    y1, y2 = Y
    dydx1 = y2
    dydx2 = -y1
    return [dydx1, dydx2]

# 初期条件
Y0 = [1, 0]  # 初期値 [y1(0), y2(0)]
x = np.linspace(0, 10, 100)

# 微分方程式の解
solution = odeint(model, Y0, x)

# プロット
plt.plot(x, solution[:, 0], label='y1')
plt.plot(x, solution[:, 1], label='y2')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of y1\'\' = -y1 and y2\' = y1\'')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 微分方程式の定義
def M(x, y):
    return x + y

def N(x, y):
    return x - y

# 完全微分方程式のチェック
def is_exact(M, N, x, y):
    dM_dy = np.gradient(M(x, y), y, edge_order=2)
    dN_dx = np.gradient(N(x, y), x, edge_order=2)
    return np.allclose(dM_dy, dN_dx)

# xとyの範囲
x = np.linspace(0, 5, 100)
y = np.linspace(0, 5, 100)

# 完全微分形であるかの確認
exact = is_exact(M, N, x, y)
print(f"Is the differential equation exact? {exact}")

# 解法の例（簡単な数値的アプローチ）
# ここではプロットを示すのみ
X, Y = np.meshgrid(x, y)
plt.streamplot(X, Y, M(X, Y), N(X, Y), color='r', linewidth=1, density=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Vector Field of M(x, y) dx + N(x, y) dy = 0')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# クレーローの微分方程式の定義
def cauchy_ode(y, x):
    return x - y

# ラグランジュの微分方程式の定義
def lagrange_ode(y, x):
    P = x
    Q = -y
    return P / Q

# 初期条件
y0 = 1
x = np.linspace(0, 5, 100)

# 微分方程式の解
solution_cauchy = odeint(cauchy_ode, y0, x)
solution_lagrange = odeint(lagrange_ode, y0, x)

# プロット
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(x, solution_cauchy, label='Cauchy\'s ODE')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of Cauchy\'s Differential Equation')
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(x, solution_lagrange, label='Lagrange\'s ODE', color='orange')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of Lagrange\'s Differential Equation')
plt.grid(True)

plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# クレーローの微分方程式の定義
def cauchy_ode(y, x):
    return x - y

# ラグランジュの微分方程式の定義
def lagrange_ode(y, x):
    P = x
    Q = -y
    return P / Q

# 初期条件
y0 = 1
x = np.linspace(0, 5, 100)

# 微分方程式の解
solution_cauchy = odeint(cauchy_ode, y0, x)
solution_lagrange = odeint(lagrange_ode, y0, x)

# プロット
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(x, solution_cauchy, label='Cauchy\'s ODE')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of Cauchy\'s Differential Equation')
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(x, solution_lagrange, label='Lagrange\'s ODE', color='orange')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of Lagrange\'s Differential Equation')
plt.grid(True)

plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 微分方程式の定義
def f(x, y):
    return x - y

# 等傾斜法の実装
def euler_method(f, y0, x0, xn, n):
    h = (xn - x0) / n
    x = np.linspace(x0, xn, n+1)
    y = np.zeros(n+1)
    y[0] = y0
    for i in range(n):
        y[i+1] = y[i] + h * f(x[i], y[i])
    return x, y

# 初期条件と設定
x0 = 0
xn = 5
y0 = 1
n = 100

# 等傾斜法での解
x, y = euler_method(f, y0, x0, xn, n)

# プロット
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Euler\'s Method')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 定数係数行列
A = np.array([[0, 1],
              [-2, -3]])

# 連立微分方程式の定義
def system(t, y):
    return np.dot(A, y)

# 初期条件と時間範囲
y0 = [1, 0]
t_span = (0, 10)
t_eval = np.linspace(*t_span, 100)

# 数値解法
sol = solve_ivp(system, t_span, y0, t_eval=t_eval)

# プロット
plt.figure(figsize=(10, 5))
plt.plot(sol.t, sol.y[0], label='y1(t)')
plt.plot(sol.t, sol.y[1], label='y2(t)')
plt.xlabel('Time t')
plt.ylabel('Solutions y1, y2')
plt.title('Constant Coefficient Linear ODE System')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 3階微分方程式の定義 (例: y''' - 3y'' + 3y' - y = 0)
def system(t, y):
    return [y[1], y[2], 3*y[2] - 3*y[1] + y[0]]

# 初期条件
y0 = [1, 0, 0]
t_span = (0, 10)
t_eval = np.linspace(*t_span, 100)

# 数値解法
sol = solve_ivp(system, t_span, y0, t_eval=t_eval)

# プロット
plt.figure(figsize=(10, 5))
plt.plot(sol.t, sol.y[0], label='y(t)')
plt.plot(sol.t, sol.y[1], label="y'(t)")
plt.plot(sol.t, sol.y[2], label="y''(t)")
plt.xlabel('Time t')
plt.ylabel('Solutions y, y\' and y\'\'')
plt.title('3rd Order ODE Converted to a System of First Order ODEs')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 3階微分方程式の定義 (例: y''' - 3y'' + 3y' - y = 0)
def system(t, y):
    return [y[1], y[2], 3*y[2] - 3*y[1] + y[0]]

# 初期条件
y0 = [1, 0, 0]
t_span = (0, 10)
t_eval = np.linspace(*t_span, 100)

# 数値解法
sol = solve_ivp(system, t_span, y0, t_eval=t_eval)

# プロット
plt.figure(figsize=(10, 5))
plt.plot(sol.t, sol.y[0], label='y(t)')
plt.plot(sol.t, sol.y[1], label="y'(t)")
plt.plot(sol.t, sol.y[2], label="y''(t)")
plt.xlabel('Time t')
plt.ylabel('Solutions y, y\' and y\'\'')
plt.title('3rd Order ODE Converted to a System of First Order ODEs')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 定数係数行列と定数項
A = np.array([[0, 1],
              [-2, -3]])
B = np.array([0, 1])

# 連立微分方程式の定義
def system(t, y):
    return np.dot(A, y) + B

# 初期条件と時間範囲
y0 = [1, 0]
t_span = (0, 10)
t_eval = np.linspace(*t_span, 100)

# 数値解法
sol = solve_ivp(system, t_span, y0, t_eval=t_eval)

# プロット
plt.figure(figsize=(10, 5))
plt.plot(sol.t, sol.y[0], label='y1(t)')
plt.plot(sol.t, sol.y[1], label='y2(t)')
plt.xlabel('Time t')
plt.ylabel('Solutions y1, y2')
plt.title('Linear ODE System with Constant Term')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 2次元自励系の定数係数行列（例: 安定性を調べる）
A = np.array([[0, 1],
              [-1, -1]])

# 連立微分方程式の定義
def system(t, y):
    return np.dot(A, y)

# 初期条件と時間範囲
y0 = [1, 0]
t_span = (0, 20)
t_eval = np.linspace(*t_span, 200)

# 数値解法
sol = solve_ivp(system, t_span, y0, t_eval=t_eval)

# プロット
plt.figure(figsize=(10, 5))
plt.plot(sol.t, sol.y[0], label='y1(t)')
plt.plot(sol.t, sol.y[1], label='y2(t)')
plt.xlabel('Time t')
plt.ylabel('Solutions y1, y2')
plt.title('2D Autonomous System')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 保存系の定数係数行列（例: ハミルトニアン系）
A = np.array([[0, 1],
              [-1, 0]])

# 連立微分方程式の定義
def system(t, y):
    return np.dot(A, y)

# 初期条件と時間範囲
y0 = [1, 0]
t_span = (0, 20)
t_eval = np.linspace(*t_span, 200)

# 数値解法
sol = solve_ivp(system, t_span, y0, t_eval=t_eval)

# プロット
plt.figure(figsize=(10, 5))
plt.plot(sol.t, sol.y[0], label='y1(t)')
plt.plot(sol.t, sol.y[1], label='y2(t)')
plt.xlabel('Time t')
plt.ylabel('Solutions y1, y2')
plt.title('Conservative System')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 偏微分方程式の定義
def pde_system(t, u):
    x, y = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
    return np.zeros_like(x)

# 解の範囲と初期条件
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)

# 初期条件としての簡単な関数
U = np.sin(X) * np.cos(Y)

# プロット
plt.figure(figsize=(10, 6))
cp = plt.contourf(X, Y, U, levels=20, cmap='viridis')
plt.colorbar(cp)
plt.xlabel('x')
plt.ylabel('y')
plt.title('1st Order PDE: $\partial u / \partial x + \partial u / \partial y = 0$')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ラグランジュの方程式の定義
def lagrange_pde(x, y):
    return np.ones_like(x)

# 解の範囲と初期条件
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)

# 初期条件としての簡単な関数
U = X + Y

# プロット
plt.figure(figsize=(10, 6))
cp = plt.contourf(X, Y, U, levels=20, cmap='viridis')
plt.colorbar(cp)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Lagrange PDE: $\partial u / \partial x + \partial u / \partial y = 1$')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 2階線形偏微分方程式の定義
def linear_pde(x, y):
    return np.sin(x) * np.cos(y)

# 解の範囲と初期条件
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)

# 初期条件としての簡単な関数
U = np.sin(X) * np.sinh(Y)

# プロット
plt.figure(figsize=(10, 6))
cp = plt.contourf(X, Y, U, levels=20, cmap='viridis')
plt.colorbar(cp)
plt.xlabel('x')
plt.ylabel('y')
plt.title('2nd Order Linear PDE: $\partial^2 u / \partial x^2 - \partial^2 u / \partial y^2 = 0$')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 定数係数2階同次偏微分方程式の定義
def homogeneous_pde(x, y):
    return np.sin(x) * np.cos(y)

# 解の範囲と初期条件
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)

# 初期条件としての簡単な関数
U = np.sin(X) * np.sin(Y)

# プロット
plt.figure(figsize=(10, 6))
cp = plt.contourf(X, Y, U, levels=20, cmap='viridis')
plt.colorbar(cp)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Constant Coefficient 2nd Order Homogeneous PDE: $\partial^2 u / \partial x^2 + \partial^2 u / \partial y^2 = 0$')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定数係数2階非同次偏微分方程式の右辺
def nonhomogeneous_pde(x, y):
    return np.sin(x) * np.cos(y)

# 解の範囲と初期条件
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)

# 初期条件としての簡単な関数
U = np.sin(X) * np.sin(Y) + np.sin(X) * np.cos(Y)

# プロット
plt.figure(figsize=(10, 6))
cp = plt.contourf(X, Y, U, levels=20, cmap='viridis')
plt.colorbar(cp)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Constant Coefficient 2nd Order Nonhomogeneous PDE: $\partial^2 u / \partial x^2 + \partial^2 u / \partial y^2 = \sin(x) \cos(y)$')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh_tridiagonal

# 固有値問題の設定
def eigenvalue_problem(N):
    # トリディアゴナル行列の作成
    main_diag = -4 * np.ones(N)
    off_diag = np.ones(N-1)
    
    # 固有値と固有ベクトルを求める
    eigenvalues, eigenvectors = eigh_tridiagonal(main_diag, off_diag)
    
    return eigenvalues, eigenvectors

# 計算
N = 50
eigenvalues, eigenvectors = eigenvalue_problem(N)

# 固有値のプロット
plt.figure(figsize=(10, 6))
plt.plot(eigenvalues, 'bo', markersize=8)
plt.xlabel('Index')
plt.ylabel('Eigenvalue')
plt.title('Eigenvalues of the Discretized Laplacian Operator')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 固有関数の定義
def eigenfunction(x, y, n, m):
    return np.sin(n * np.pi * x) * np.sin(m * np.pi * y)

# グリッド設定
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)

# 固有関数の計算
n, m = 1, 2
U = eigenfunction(X, Y, n, m)

# プロット
plt.figure(figsize=(10, 6))
cp = plt.contourf(X, Y, U, levels=20, cmap='viridis')
plt.colorbar(cp)
plt.xlabel('x')
plt.ylabel('y')
plt.title(f'Eigenfunction for n={n}, m={m}')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# 変分法の目的関数
def objective(u, x, dx):
    du_dx = np.gradient(u, dx)
    return np.sum((du_dx**2 - 2 * u**2)) * dx

# 設定
L = 1
N = 100
x = np.linspace(0, L, N)
dx = x[1] - x[0]
initial_guess = np.sin(np.pi * x)  # 初期推定

# 最適化
result = minimize(objective, initial_guess, args=(x, dx), method='BFGS')
u_optimal = result.x

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x, u_optimal, label='Optimal Function')
plt.xlabel('x')
plt.ylabel('u(x)')
plt.title('Optimal Function for Variational Problem')
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# グリーン関数の定義
def green_function(x, x0, L):
    return np.where(x < x0, x * (1 - x0), x0 * (1 - x))

# 設定
L = 1
x = np.linspace(0, L, 100)
x0 = 0.5

# グリーン関数の計算
G = green_function(x, x0, L)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x, G, label='Green Function')
plt.xlabel('x')
plt.ylabel('G(x)')
plt.title(f'Green Function for x0={x0}')
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# パラメータ
a, b, c = 1, 0, -1  # 例: y'' - y = g(x)
def g(x):
    return np.sin(x)  # 例: 非同次項

# 微分方程式
def model(y, x):
    dydx = [y[1], -b/a * y[1] - c/a * y[0] + g(x)/a]
    return dydx

# 初期条件
y0 = [0, 0]  # 初期位置と初期速度
x = np.linspace(0, 10, 100)
solution = odeint(model, y0, x)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x, solution[:, 0], label='y(x)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of Non-Homogeneous ODE')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# パラメータ
a, b, c = 1, -2, 1  # 例: y'' - 2y' + y = 0

# 微分方程式
def model(y, x):
    dydx = [y[1], -b/a * y[1] - c/a * y[0]]
    return dydx

# 初期条件
y0 = [1, 0]  # 初期位置と初期速度
x = np.linspace(0, 10, 100)
solution = odeint(model, y0, x)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x, solution[:, 0], label='y(x)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of Constant Coefficient 2nd Order Linear ODE')
plt.legend()
plt.grid(True)
plt.show()

import sympy as sp

# 記号変数の定義
x = sp.Symbol('x')
y = sp.Function('y')(x)

# 微分方程式の定義
deqn = sp.Eq(y.diff(x, x, x) - y, 0)

# 解の級数展開
sol = sp.series(sp.dsolve(deqn, y).rhs, x, n=10)
print(f'Series Solution: {sol}')

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import legendre

# ルジャンドル関数の計算
def plot_legendre(n):
    x = np.linspace(-1, 1, 400)
    Pn = legendre(n)
    y = Pn(x)

    plt.figure(figsize=(10, 6))
    plt.plot(x, y, label=f'Legendre Polynomial P_{n}(x)')
    plt.xlabel('x')
    plt.ylabel('P_n(x)')
    plt.title(f'Legendre Polynomial of Degree {n}')
    plt.legend()
    plt.grid(True)
    plt.show()

# プロット
plot_legendre(3)  # 例: 3次ルジャンドル多項式

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn

# ベッセル関数の計算
def plot_bessel(n):
    x = np.linspace(0, 20, 400)
    y = jn(n, x)

    plt.figure(figsize=(10, 6))
    plt.plot(x, y, label=f'Bessel Function J_{n}(x)')
    plt.xlabel('x')
    plt.ylabel(f'J_{n}(x)')
    plt.title(f'Bessel Function of Order {n}')
    plt.legend()
    plt.grid(True)
    plt.show()

# プロット
plot_bessel(0)  # 例: ベッセル関数の0次

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq

# パラメータ
T = 1.0  # サンプル間隔
N = 1000  # サンプル数
x = np.linspace(0.0, N*T, N, endpoint=False)
f = np.sin(2*np.pi*50.0*x) + 0.5*np.sin(2*np.pi*80.0*x)

# フーリエ変換
yf = fft(f)
xf = fftfreq(N, T)[:N//2]

# フーリエ逆変換
f_rec = ifft(yf)

# プロット
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(x, f, label='Original Signal')
plt.title('Original Signal')
plt.grid()

plt.subplot(3, 1, 2)
plt.plot(xf, np.abs(yf[:N//2]), label='Fourier Transform')
plt.title('Fourier Transform')
plt.grid()

plt.subplot(3, 1, 3)
plt.plot(x, f_rec, label='Reconstructed Signal')
plt.title('Reconstructed Signal')
plt.grid()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import legendre

# ルジャンドル多項式の計算
def plot_legendre_expansion(x, n_terms):
    x_vals = np.linspace(-1, 1, 400)
    y_vals = np.zeros_like(x_vals)
    for n in range(n_terms):
        Pn = legendre(n)
        y_vals += Pn(x_vals)

    plt.figure(figsize=(8, 6))
    plt.plot(x_vals, y_vals, label=f'Legendre Series Expansion ({n_terms} terms)')
    plt.xlabel('x')
    plt.ylabel('P_n(x)')
    plt.title('Legendre Polynomial Expansion')
    plt.legend()
    plt.grid(True)
    plt.show()

# プロット
plot_legendre_expansion(np.linspace(-1, 1, 400), 5)  # 例: 5項の展開

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn, jn_zeros

# パラメータ
n = 0  # ベッセル関数の次数

# ベッセル関数の零点
zeros = jn_zeros(n, 10)  # 最初の10個の零点

# プロット
x = np.linspace(0, 15, 400)
y = jn(n, x)

plt.figure(figsize=(10, 6))
plt.plot(x, y, label=f'Bessel Function J_{n}(x)')
plt.scatter(zeros, jn(n, zeros), color='red', label='Zeros')
plt.xlabel('x')
plt.ylabel(f'J_{n}(x)')
plt.title('Bessel Function and its Zeros')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
L = 10.0
T = 2.0
Nx = 100
Nt = 100
dx = L / (Nx - 1)
dt = T / Nt
c = 1.0
x = np.linspace(0, L, Nx)
u = np.zeros((Nt, Nx))
u[0, :] = np.sin(np.pi * x / L)  # 初期条件

# 数値解法
for n in range(0, Nt-1):
    for i in range(1, Nx-1):
        u[n+1, i] = (2 - 2 * (c * dt / dx)**2) * u[n, i] - u[n-1, i] + (c * dt / dx)**2 * (u[n, i+1] + u[n, i-1])

# プロット
plt.figure(figsize=(10, 6))
plt.imshow(u, extent=[0, L, 0, T], aspect='auto', cmap='jet')
plt.colorbar(label='u(x, t)')
plt.xlabel('x')
plt.ylabel('t')
plt.title('1D Wave Equation Solution')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
Lx, Ly = 10.0, 10.0
T = 2.0
Nx, Ny = 100, 100
Nt = 100
dx, dy = Lx / (Nx - 1), Ly / (Ny - 1)
dt = T / Nt
c = 1.0
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y)
u = np.zeros((Nt, Nx, Ny))
u[0, :, :] = np.sin(np.pi * X / Lx) * np.sin(np.pi * Y / Ly)  # 初期条件

# 数値解法
for n in range(0, Nt-1):
    for i in range(1, Nx-1):
        for j in range(1, Ny-1):
            u[n+1, i, j] = (2 - 2 * (c * dt / dx)**2 - 2 * (c * dt / dy)**2) * u[n, i, j] - u[n-1, i, j] + (c * dt / dx)**2 * (u[n, i+1, j] + u[n, i-1, j]) + (c * dt / dy)**2 * (u[n, i, j+1] + u[n, i, j-1])

# プロット
plt.figure(figsize=(10, 6))
plt.imshow(u[int(Nt/2), :, :], extent=[0, Lx, 0, Ly], cmap='jet')
plt.colorbar(label='u(x, y)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('2D Wave Equation Solution at t = T/2')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
Lx, Ly = 10.0, 10.0
T = 2.0
Nx, Ny = 100, 100
Nt = 100
dx, dy = Lx / (Nx - 1), Ly / (Ny - 1)
dt = T / Nt
alpha = 0.01
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y)
u = np.zeros((Nt, Nx, Ny))
u[0, :, :] = np.sin(np.pi * X / Lx) * np.sin(np.pi * Y / Ly)  # 初期条件

# 数値解法
for n in range(0, Nt-1):
    u[n+1, 1:-1, 1:-1] = u[n, 1:-1, 1:-1] + alpha * dt * (
        (u[n, 2:, 1:-1] - 2 * u[n, 1:-1, 1:-1] + u[n, :-2, 1:-1]) / dx**2 +
        (u[n, 1:-1, 2:] - 2 * u[n, 1:-1, 1:-1] + u[n, 1:-1, :-2]) / dy**2
    )

# プロット
plt.figure(figsize=(10, 6))
plt.imshow(u[int(Nt

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 微分方程式の定義
def dydx(y, x):
    p = 1  # p(x) = 1
    q = np.sin(x)  # q(x) = sin(x)
    n = 2  # n = 2
    return (q * y**n - p * y) / (y**(n-1))

# 初期条件とxの範囲
y0 = 1
x = np.linspace(0, 10, 100)

# 微分方程式の数値解法
y = odeint(dydx, y0, x)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x, y, label="Solution to dy/dx + y = sin(x) y^2")
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution to Bernoulli Differential Equation')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 関数の定義
def M(x, y):
    return x + y

def N(x, y):
    return x - y

# 確認
def is_exact(M, N):
    return np.all(np.gradient(M, axis=1) == np.gradient(N, axis=0))

# 計算
x = np.linspace(0, 10, 100)
y = np.linspace(0, 10, 100)
X, Y = np.meshgrid(x, y)

M_values = M(X, Y)
N_values = N(X, Y)

# プロット
plt.figure(figsize=(10, 6))
plt.contourf(X, Y, M_values, alpha=0.5, cmap='jet')
plt.colorbar(label='M(x, y)')
plt.contourf(X, Y, N_values, alpha=0.5, cmap='viridis')
plt.colorbar(label='N(x, y)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Exact Differential Form Visualization')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 微分不等式の定義
def dydx(y, x):
    return np.sin(x)  # f(x, y) = sin(x)

# グロンウォールの不等式の右辺
def gronwall_rhs(y, x):
    return np.sin(x) * y  # f(x) = sin(x)

# 初期条件とxの範囲
y0 = 1
x = np.linspace(0, 10, 500)

# 微分方程式の数値解法
y = odeint(dydx, y0, x)

# グロンウォールの不等式の評価
y_gw = y0 * np.exp(np.cumsum(np.sin(x)) * (x[1] - x[0]))

# プロット
plt.figure(figsize=(12, 6))
plt.plot(x, y, label='Solution to dy/dx = sin(x)', color='blue')
plt.plot(x, y_gw, label='Gronwall\'s Inequality', color='red', linestyle='--')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution and Gronwall\'s Inequality')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 例: 微分方程式の定義
def model(y, x):
    y1, y2 = y
    dy1_dx = y2
    dy2_dx = -2 * y2 - y1
    return [dy1_dx, dy2_dx]

# 初期条件とxの範囲
y0 = [1, 0]  # 初期値
x = np.linspace(0, 10, 500)

# 微分方程式の数値解法
solution = odeint(model, y0, x)
y1, y2 = solution[:, 0], solution[:, 1]

# ロンスキアンの計算
def wronskian(y1, y2, dy1_dx, dy2_dx):
    return y1 * dy2_dx - y2 * dy1_dx

dy1_dx = np.gradient(y1, x)
dy2_dx = np.gradient(y2, x)
W = wronskian(y1, y2, dy1_dx, dy2_dx)

# プロット
plt.figure(figsize=(12, 6))
plt.plot(x, y1, label='y1(x)')
plt.plot(x, y2, label='y2(x)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solutions of the Differential Equation')
plt.grid(True)
plt.legend()
plt.show()

plt.figure(figsize=(12, 6))
plt.plot(x, W, label='Wronskian', color='red')
plt.xlabel('x')
plt.ylabel('W')
plt.title('Wronskian of the Solutions')
plt.grid(True)
plt.legend()
plt.show()

import sympy as sp

# シンボリックな変数と関数の定義
x = sp.Symbol('x')
y = sp.Function('y')(x)

# 定係数2階線形微分方程式の定義
eq = sp.Eq(sp.Derivative(y, x, x) - 2 * sp.Derivative(y, x) + y, 0)

# 微分方程式の解
solution = sp.dsolve(eq, y)
solution

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 定係数 n 階線形微分方程式の例として、3階微分方程式を定義
def model(y, x):
    y1, y2, y3 = y
    dy1_dx = y2
    dy2_dx = y3
    dy3_dx = -2 * y3 - y2 - y1
    return [dy1_dx, dy2_dx, dy3_dx]

# 初期条件とxの範囲
y0 = [1, 0, 0]  # 初期値
x = np.linspace(0, 10, 500)

# 微分方程式の数値解法
solution = odeint(model, y0, x)
y1, y2, y3 = solution[:, 0], solution[:, 1], solution[:, 2]

# プロット
plt.figure(figsize=(12, 6))
plt.plot(x, y1, label='y1(x)')
plt.plot(x, y2, label='y2(x)')
plt.plot(x, y3, label='y3(x)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solutions of the 3rd Order Differential Equation')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 3階線形微分方程式の定義
def model(y, x):
    y1, y2, y3 = y
    dy1_dx = y2
    dy2_dx = y3
    dy3_dx = -2 * y3 - y2 - y1
    return [dy1_dx, dy2_dx, dy3_dx]

# 初期条件とxの範囲
y0 = [1, 0, 0]  # 初期値
x = np.linspace(0, 10, 500)

# 微分方程式の数値解法
solution = odeint(model, y0, x)
y1, y2, y3 = solution[:, 0], solution[:, 1], solution[:, 2]

# プロット
plt.figure(figsize=(12, 6))
plt.plot(x, y1, label='y1(x)')
plt.plot(x, y2, label='y2(x)')
plt.plot(x, y3, label='y3(x)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solutions of the 3rd Order Differential Equation')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from matplotlib.patches import FancyArrowPatch

# 6.2 線形連立微分方程式の定義
def linear_system(y, t):
    # 状態ベクトル y = [y1, y2]
    y1, y2 = y
    # 行列 A による線形連立微分方程式
    dydt = [2*y1 - y2, y1 + 2*y2]
    return dydt

# 時間の範囲
t = np.linspace(0, 10, 500)

# 初期条件
initial_conditions = [[1, 0], [0, 1], [1, 1], [-1, -1]]

# 数値解法
plt.figure(figsize=(12, 6))
for ic in initial_conditions:
    solution = odeint(linear_system, ic, t)
    plt.plot(solution[:, 0], solution[:, 1], label=f'IC: {ic}')

plt.xlabel('y1')
plt.ylabel('y2')
plt.title('Phase Portrait of the Linear System')
plt.grid(True)
plt.legend()
plt.show()

# 6.3 解の安定性についての分析
# 固有値の計算
A = np.array([[2, -1], [1, 2]])
eigenvalues, eigenvectors = np.linalg.eig(A)

print("Eigenvalues:", eigenvalues)
print("Eigenvectors:", eigenvectors)

# 6.4 極方程式と解の挙動
# 極方程式の定義
def polar_trajectory(r, theta):
    # r の関数として θ の変化を定義
    return np.array([np.cos(theta), np.sin(theta)])

theta_vals = np.linspace(0, 2 * np.pi, 500)
r = np.ones_like(theta_vals)
polar_coords = np.array([polar_trajectory(r_val, theta) for r_val, theta in zip(r, theta_vals)])

plt.figure(figsize=(12, 6))
plt.polar(theta_vals, r, label='Polar trajectory')
plt.title('Polar Coordinates of the System')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn  # ベッセル関数

# 7.3 ベッセルの微分方程式
def bessel_function(nu, x):
    return jn(nu, x)

# x の範囲
x = np.linspace(0, 20, 500)

# ベッセル関数の階数
nu_values = [0, 1, 2, 3]

# プロット
plt.figure(figsize=(12, 8))
for nu in nu_values:
    y = bessel_function(nu, x)
    plt.plot(x, y, label=f'ν = {nu}')

plt.xlabel('x')
plt.ylabel('Jν(x)')
plt.title('Bessel Functions of the First Kind')
plt.legend()
plt.grid(True)
plt.show()

# 7.2 よく知られた微分方程式のべき級数解（例：ラプラス方程式）
def power_series_solution(x, terms=10):
    # ベキ級数展開の計算
    a = np.zeros(terms)
    a[0] = 1  # 初期条件
    for n in range(1, terms):
        a[n] = -a[n-1] * (x**2 / (2 * n))
    return a

# x の範囲と級数の計算
x_values = np.linspace(0, 2, 100)
terms = 10
series_solutions = np.zeros((len(x_values), terms))
for i, x_val in enumerate(x_values):
    series_solutions[i, :] = power_series_solution(x_val, terms)

# プロット
plt.figure(figsize=(12, 8))
for n in range(terms):
    plt.plot(x_values, series_solutions[:, n], label=f'Term {n+1}')

plt.xlabel('x')
plt.ylabel('Series Terms')
plt.title('Power Series Solution')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq

# 9.1 放物型方程式の混合問題（熱方程式の例）
def heat_equation_solution(x, t, L=10, N=100):
    """
    放物型方程式 (熱方程式) の数値解
    """
    x = np.linspace(0, L, N)
    u = np.zeros(N)
    
    for n in range(1, 10):
        lambda_n = (n * np.pi / L)**2
        u += (2 / (n * np.pi)) * np.sin(n * np.pi * x / L) * np.exp(-lambda_n * t)
        
    return u

# 時間と位置の範囲
x = np.linspace(0, 10, 100)
t = [0.1, 0.5, 1.0]

plt.figure(figsize=(12, 8))
for time in t:
    plt.plot(x, heat_equation_solution(x, time), label=f't = {time}')

plt.xlabel('x')
plt.ylabel('u(x, t)')
plt.title('Heat Equation Solution')
plt.legend()
plt.grid(True)
plt.show()

# 9.2 フーリエ級数
def fourier_series(x, N=10):
    """
    フーリエ級数展開
    """
    a0 = 0
    a_n = lambda n: 2 / (n * np.pi)
    b_n = lambda n: 0
    
    result = np.zeros_like(x)
    for n in range(1, N+1):
        result += a_n(n) * np.sin(n * np.pi * x)
    
    return result

# x の範囲
x = np.linspace(0, 2, 500)

# フーリエ級数のプロット
plt.figure(figsize=(12, 8))
plt.plot(x, fourier_series(x, N=10), label='Fourier Series (N=10)')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Fourier Series Expansion')
plt.legend()
plt.grid(True)
plt.show()

# 9.3 双曲型方程式の混合問題（波動方程式の例）
def wave_equation_solution(x, t, L=10, N=100):
    """
    双曲型方程式 (波動方程式) の数値解
    """
    x = np.linspace(0, L, N)
    u = np.zeros(N)
    
    for n in range(1, 10):
        lambda_n = n * np.pi / L
        u += (4 / (n * np.pi)) * np.sin(n * np.pi * x / L) * np.cos(lambda_n * t)
        
    return u

# 時間と位置の範囲
x = np.linspace(0, 10, 100)
t = [0.1, 0.5, 1.0]

plt.figure(figsize=(12, 8))
for time in t:
    plt.plot(x, wave_equation_solution(x, time), label=f't = {time}')

plt.xlabel('x')
plt.ylabel('u(x, t)')
plt.title('Wave Equation Solution')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 7.1 1階微分方程式の解の追跡
def first_order_ODE(t, y):
    return y - y**3

t_span = (0, 10)
y0 = [0.1, 1, -1]
t_eval = np.linspace(*t_span, 500)

plt.figure(figsize=(12, 8))
for y_init in y0:
    sol = solve_ivp(first_order_ODE, t_span, [y_init], t_eval=t_eval)
    plt.plot(sol.t, sol.y[0], label=f'Initial y={y_init}')

plt.xlabel('Time')
plt.ylabel('y')
plt.title('1st Order ODE Solution Tracking')
plt.legend()
plt.grid(True)
plt.show()

# 7.2 2次元自励系の軌道と特異点
def self_excited_system(t, state):
    x, y = state
    dxdt = x - y
    dydt = x + 0.2 * y
    return [dxdt, dydt]

t_span = (0, 20)
state0 = [0.1, 0.1]
t_eval = np.linspace(*t_span, 500)

sol = solve_ivp(self_excited_system, t_span, state0, t_eval=t_eval)

plt.figure(figsize=(12, 8))
plt.plot(sol.y[0], sol.y[1])
plt.xlabel('x')
plt.ylabel('y')
plt.title('2D Self-Excited System Trajectory')
plt.grid(True)
plt.show()

# 7.3 安定性と漸近安定性
def stability_system(t, state):
    x, y = state
    dxdt = -x + y
    dydt = -x - y
    return [dxdt, dydt]

t_span = (0, 10)
state0 = [1, 1]
t_eval = np.linspace(*t_span, 500)

sol = solve_ivp(stability_system, t_span, state0, t_eval=t_eval)

plt.figure(figsize=(12, 8))
plt.plot(sol.t, sol.y[0], label='x(t)')
plt.plot(sol.t, sol.y[1], label='y(t)')
plt.xlabel('Time')
plt.ylabel('State')
plt.title('Stability and Asymptotic Stability')
plt.legend()
plt.grid(True)
plt.show()

# 7.4 周期係数の方程式の解の挙動
def periodic_coeff_system(t, state):
    x, y = state
    dxdt = np.sin(t) - x
    dydt = np.cos(t) - y
    return [dxdt, dydt]

t_span = (0, 10)
state0 = [0, 0]
t_eval = np.linspace(*t_span, 500)

sol = solve_ivp(periodic_coeff_system, t_span, state0, t_eval=t_eval)

plt.figure(figsize=(12, 8))
plt.plot(sol.t, sol.y[0], label='x(t)')
plt.plot(sol.t, sol.y[1], label='y(t)')
plt.xlabel('Time')
plt.ylabel('State')
plt.title('Periodic Coefficient System Behavior')
plt.legend()
plt.grid(True)
plt.show()

# 7.5 2次元自励系の極限閉軌道
def limit_cycle_system(t, state):
    x, y = state
    dxdt = y - x**2
    dydt = -x - y
    return [dxdt, dydt]

t_span = (0, 20)
state0 = [1, 1]
t_eval = np.linspace(*t_span, 500)

sol = solve_ivp(limit_cycle_system, t_span, state0, t_eval=t_eval)

plt.figure(figsize=(12, 8))
plt.plot(sol.y[0], sol.y[1])
plt.xlabel('x')
plt.ylabel('y')
plt.title('2D Self-Excited System Limit Cycle')
plt.grid(True)
plt.show()

# 7.6 高次元の力学系とアトラクタ
def lorenz_system(t, state):
    x, y, z = state
    sigma = 10
    rho = 28
    beta = 8 / 3
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]

t_span = (0, 100)
state0 = [0, 1, 1.05]
t_eval = np.linspace(*t_span, 10000)

sol = solve_ivp(lorenz_system, t_span, state0, t_eval=t_eval)

plt.figure(figsize=(12, 8))
ax = plt.axes(projection='3d')
ax.plot(sol.y[0], sol.y[1], sol.y[2])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Lorenz Attractor')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Constants
T_ambient = 20  # Ambient temperature (°C)
k = 0.1  # Cooling constant (1/min)

# Differential equation
def model(T, t):
    dTdt = -k * (T - T_ambient)
    return dTdt

# Initial conditions
T0 = 90  # Initial temperature (°C)
time = np.linspace(0, 60, 100)  # Time (minutes)

# Numerical solution
T = odeint(model, T0, time)

# Plot
plt.figure(figsize=(10, 6))
plt.plot(time, T, label='Coffee Temperature')
plt.axhline(y=T_ambient, color='r', linestyle='--', label='Ambient Temperature')
plt.xlabel('Time (minutes)')
plt.ylabel('Temperature (°C)')
plt.title('Cooling of Coffee')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Constants
g = 9.81  # Gravitational acceleration (m/s^2)
k = 0.1   # Air resistance coefficient (kg/m)

# Differential equations
def model(y, t):
    v, h = y
    dvdt = -g - k * v
    dhdt = v
    return [dvdt, dhdt]

# Initial conditions
v0 = 0  # Initial velocity (m/s)
h0 = 0  # Initial height (m)
initial_conditions = [v0, h0]
time = np.linspace(0, 100, 1000)

# Numerical solution
result = odeint(model, initial_conditions, time)
velocity = result[:, 0]
height = result[:, 1]

# Plot
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(time, velocity)
plt.xlabel('Time (seconds)')
plt.ylabel('Velocity (m/s)')
plt.title('Rocket Velocity')

plt.subplot(2, 1, 2)
plt.plot(time, height)
plt.xlabel('Time (seconds)')
plt.ylabel('Height (m)')
plt.title('Rocket Height')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Constants
r = 0.5  # Growth rate
K = 1000  # Carrying capacity

# Differential equation
def model(P, t):
    dPdt = r * P * (1 - P / K)
    return dPdt

# Initial conditions
P0 = 10  # Initial quantity
time = np.linspace(0, 20, 100)

# Numerical solution
P = odeint(model, P0, time)

# Plot
plt.figure(figsize=(10, 6))
plt.plot(time, P, label='Effect of Advertising')
plt.xlabel('Time (units)')
plt.ylabel('Effect (units)')
plt.title('Logistic Growth Model')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, Function, Eq, dsolve, exp

# Define symbols
x = symbols('x')
y = Function('y')(x)
P = symbols('P', real=True)
Q = symbols('Q', real=True)

# Differential equation
P = -1/x
Q = x
diff_eq = Eq(y.diff(x) + P*y, Q)

# Integrating factor
integrating_factor = exp(P.integrate(x))
solution = dsolve(diff_eq, y)

# Display result
print(f'Solution: {solution}')

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define the system of differential equations
def model(Y, t):
    x, y = Y
    dxdt = x + 2*y
    dydt = 3*x + 4*y
    return [dxdt, dydt]

# Initial conditions
x0 = 1
y0 = 0
initial_conditions = [x0, y0]

# Time vector
time = np.linspace(0, 10, 100)

# Numerical solution
solution = odeint(model, initial_conditions, time)
x = solution[:, 0]
y = solution[:, 1]

# Plot
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(time, x, label='x(t)')
plt.xlabel('Time (t)')
plt.ylabel('x(t)')
plt.title('Solution of Coupled Differential Equations')
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(time, y, label='y(t)', color='orange')
plt.xlabel('Time (t)')
plt.ylabel('y(t)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define Bernoulli's differential equation
def bernoulli(y, t):
    dydt = -2*y + 5*y**2
    return dydt

# Define Riccati's differential equation
def riccati(y, t):
    dydt = 1 + 2*y + y**2
    return dydt

# Time vector
t = np.linspace(0, 5, 100)

# Initial condition
y0_bernoulli = 0.5
y0_riccati = 0.5

# Solve Bernoulli's equation
solution_bernoulli = odeint(bernoulli, y0_bernoulli, t)
y_bernoulli = solution_bernoulli[:, 0]

# Solve Riccati's equation
solution_riccati = odeint(riccati, y0_riccati, t)
y_riccati = solution_riccati[:, 0]

# Plot Bernoulli's equation
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(t, y_bernoulli, label='Bernoulli\'s Equation y(t)')
plt.xlabel('Time (t)')
plt.ylabel('y(t)')
plt.title('Solution of Bernoulli\'s Equation')
plt.legend()
plt.grid(True)

# Plot Riccati's equation
plt.subplot(2, 1, 2)
plt.plot(t, y_riccati, label='Riccati\'s Equation y(t)', color='orange')
plt.xlabel('Time (t)')
plt.ylabel('y(t)')
plt.title('Solution of Riccati\'s Equation')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define the differential equation
def constant_variation(y, t):
    y1, y2 = y
    dydt = [y2, 4*y2 - 4*y1 + t * np.exp(t)]
    return dydt

# Time vector
t = np.linspace(0, 10, 100)

# Initial conditions
y0 = [0, 0]  # y(0) and y'(0)

# Solve the differential equation
solution = odeint(constant_variation, y0, t)
y = solution[:, 0]

# Plot
plt.figure(figsize=(8, 6))
plt.plot(t, y, label='y(t)')
plt.xlabel('Time (t)')
plt.ylabel('y(t)')
plt.title('Solution using Method of Constant Variation')
plt.legend()
plt.grid(True)
plt.show()

import sympy as sp

# Define symbols
x = sp.Symbol('x')
y = sp.Function('y')(x)

# Define the differential equation
diff_eq = x**2 * y.diff(x, x) - x * y.diff(x) + (x**2 - 1) * y

# Identify the regular and singular points
singular_points = sp.singularities(diff_eq, x)
print(f'Singular points: {singular_points}')

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import chebyt

# Define Chebyshev polynomials
x = np.linspace(-1, 1, 400)
for n in range(4):
    Tn = chebyt(n)
    plt.plot(x, Tn(x), label=f'T_{n}(x)')

# Plot
plt.xlabel('x')
plt.ylabel('T_n(x)')
plt.title('Chebyshev Polynomials')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

# Define symbols
x = sp.Symbol('x')
y = sp.Function('y')(x)
y_prime = y.diff(x)

# Define the Lagrangian
F = (1/2) * y_prime**2 + (1/2) * y**2

# Calculate the Euler-Lagrange equation
L = sp.lambdify((x, y, y_prime), F)
dF_dy = sp.diff(F, y)
dF_dy_prime = sp.diff(F, y_prime)
d_dF_dy_prime = sp.diff(dF_dy_prime, x)
euler_lagrange_eq = dF_dy - d_dF_dy_prime
euler_lagrange_eq_simplified = sp.simplify(euler_lagrange_eq)

# Solve the differential equation
solution = sp.dsolve(euler_lagrange_eq_simplified, y)
print("Solution of the Euler-Lagrange Equation:")
print(solution)

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Define the differential equation
def geodesic(y, x):
    dy_dx, dy2_dx2 = y
    return [dy_dx, -dy_dx / (1 + dy_dx**2)**0.5]

# Initial conditions and range
y0 = [0, 1]  # Initial position and slope
x = np.linspace(0, 10, 100)

# Solve the differential equation
sol = odeint(geodesic, y0, x)
y = sol[:, 0]

# Plot the result
plt.figure(figsize=(8, 6))
plt.plot(x, y, label='Shortest Path')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Shortest Path (Geodesic) on a Plane')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# Define the functional
def functional(y, x):
    return np.trapz(y**2 + np.gradient(y)**2, x)

# Define the range and initial guess
x = np.linspace(0, 10, 100)
initial_y = np.sin(x)  # Initial guess

# Define the objective function for optimization
def objective_function(y):
    y = np.reshape(y, x.shape)
    return functional(y, x)

# Optimize
result = minimize(objective_function, initial_y, method='BFGS')
optimized_y = result.x

# Plot the result
plt.figure(figsize=(8, 6))
plt.plot(x, optimized_y, label='Optimized Path')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Minimized Functional Path')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

# Define the matrix A
A = np.array([[2, -1], [1, 1]])

# Define the initial conditions
y0 = np.array([1, 0])

# Time range
t = np.linspace(0, 10, 100)

# Solve the system
def solution(t, A, y0):
    return np.dot(expm(A*t), y0)

sol = np.array([solution(ti, A, y0) for ti in t])

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(t, sol[:, 0], label='x(t)')
plt.plot(t, sol[:, 1], label='y(t)')
plt.xlabel('Time')
plt.ylabel('Values')
plt.title('Solution of the Homogeneous System with Constant Coefficients')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import beta

# Define x and y values
x = np.linspace(0.1, 5, 100)
y = np.linspace(0.1, 5, 100)

# Compute the Beta function
X, Y = np.meshgrid(x, y)
Z = beta(X, Y)

# Plot the Beta function
plt.figure(figsize=(10, 7))
cp = plt.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar(cp)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Beta Function B(x, y)')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma

# Define x values
x = np.linspace(0.1, 5, 100)

# Compute the Gamma function
y = gamma(x)

# Plot the Gamma function
plt.figure(figsize=(10, 6))
plt.plot(x, y, label='Gamma Function')
plt.xlabel('x')
plt.ylabel('Gamma(x)')
plt.title('Gamma Function Γ(x)')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Define the vector field (dx/dt, dy/dt)
def vector_field(t, y):
    x, y = y
    dxdt = -y
    dydt = x
    return [dxdt, dydt]

# Define the range for t and initial conditions
t_span = (0, 10)
initial_conditions = [[1, 0], [0, 1], [-1, 0], [0, -1]]

# Solve the differential equation for different initial conditions
plt.figure(figsize=(8, 8))
for initial in initial_conditions:
    sol = solve_ivp(vector_field, t_span, initial, t_eval=np.linspace(0, 10, 100))
    plt.plot(sol.y[0], sol.y[1], label=f'IC: {initial}')

plt.xlabel('x')
plt.ylabel('y')
plt.title('Integral Curves of the Vector Field')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

# Define the system of differential equations (dx/dt, dy/dt)
def system_equations(variables):
    x, y = variables
    dxdt = x - y
    dydt = x + y - 2
    return [dxdt, dydt]

# Find equilibrium points by solving the system of equations
equilibrium_points = fsolve(system_equations, [0, 0])

# Define a grid for plotting
x = np.linspace(-2, 2, 20)
y = np.linspace(-2, 2, 20)
X, Y = np.meshgrid(x, y)
U, V = system_equations([X, Y])

# Plot vector field and equilibrium points
plt.figure(figsize=(8, 8))
plt.quiver(X, Y, U, V, color='b', headlength=4)
plt.scatter(equilibrium_points[0], equilibrium_points[1], color='r', label='Equilibrium Point')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Vector Field with Equilibrium Points')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
c = 1  # wave speed
L = 10  # length of the domain
T = 5   # total time
Nx = 100  # number of spatial points
Nt = 200  # number of time points

# Define the initial conditions
def f(x):
    return np.sin(np.pi * x / L)  # initial displacement

def g(x):
    return np.zeros_like(x)  # initial velocity

# Define spatial and time grids
x = np.linspace(0, L, Nx)
t = np.linspace(0, T, Nt)

# Create a meshgrid for plotting
X, T = np.meshgrid(x, t)

# Compute the solution using D'Alembert's formula
U = np.zeros_like(X)
for i in range(Nx):
    for j in range(Nt):
        xi = X[j, i] - c * T[j, i]
        xj = X[j, i] + c * T[j, i]
        integral = np.trapz(g(np.linspace(xi, xj, 100)), np.linspace(xi, xj, 100))
        if xi < 0:
            xi = 0
        if xj > L:
            xj = L
        U[j, i] = 0.5 * (f(xi) + f(xj)) + (1 / (2 * c)) * integral

# Plot the solution
plt.figure(figsize=(10, 6))
plt.contourf(X, T, U, 20, cmap='viridis')
plt.colorbar(label='u(x, t)')
plt.xlabel('x')
plt.ylabel('t')
plt.title('D\'Alembert\'s Solution of the Wave Equation')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
c = 1  # wave speed
L = 10  # length of the domain
T = 5   # total time
Nx = 100  # number of spatial points
Nt = 200  # number of time points
N_terms = 10  # number of Fourier series terms

# Define the initial conditions
def f(x):
    return np.sin(np.pi * x / L)  # initial displacement

def B_n(n):
    return 2 * np.trapz(f(np.linspace(0, L, Nx)) * np.sin(n * np.pi * np.linspace(0, L, Nx) / L), np.linspace(0, L, Nx)) / L

# Define spatial and time grids
x = np.linspace(0, L, Nx)
t = np.linspace(0, T, Nt)

# Create a meshgrid for plotting
X, T = np.meshgrid(x, t)

# Compute the solution using Fourier's method
U = np.zeros_like(X)
for n in range(1, N_terms + 1):
    B = B_n(n)
    U += B * np.sin(n * np.pi * X / L) * np.sin(n * np.pi * c * T / L)

# Plot the solution
plt.figure(figsize=(10, 6))
plt.contourf(X, T, U, 20, cmap='viridis')
plt.colorbar(label='u(x, t)')
plt.xlabel('x')
plt.ylabel('t')
plt.title('Fourier Series Solution of the Wave Equation')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
R = 1  # radius of the circular domain
Nx = 100  # number of spatial points
Ny = 100  # number of spatial points

# Define the boundary condition
def boundary_condition(theta):
    return np.sin(theta)  # example: u = sin(theta) on the boundary

# Define spatial grid
x = np.linspace(-R, R, Nx)
y = np.linspace(-R, R, Ny)
X, Y = np.meshgrid(x, y)
R_sq = X**2 + Y**2

# Compute the solution using Poisson's formula
U = np.zeros_like(X)
for i in range(Nx):
    for j in range(Ny):
        if R_sq[i, j] <= R**2:
            r = np.sqrt(R_sq[i, j])
            theta = np.arctan2(Y[i, j], X[i, j])
            integral = 0
            for t in np.linspace(0, 2 * np.pi, 100):
                x_b = R * np.cos(t)
                y_b = R * np.sin(t)
                ds = R * np.sqrt((X[i, j] - x_b)**2 + (Y[i, j] - y_b)**2)
                integral += boundary_condition(t) / ds
            U[i, j] = integral / (2 * np.pi)

# Plot the solution
plt.figure(figsize=(8, 8))
plt.contourf(X, Y, U, 20, cmap='viridis')
plt.colorbar(label='u(x, y)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution Using Poisson\'s Formula')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
L = 1.0  # length of the domain
Nx = 50  # number of spatial points
Ny = 50  # number of spatial points
tolerance = 1e-5  # convergence tolerance
max_iter = 10000  # maximum number of iterations

# Define the grid
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

# Initialize the potential field
U = np.zeros((Nx, Ny))

# Boundary conditions
U[:, 0] = 0  # u = 0 on left boundary
U[:, -1] = 0  # u = 0 on right boundary
U[0, :] = 1  # u = 1 on top boundary
U[-1, :] = 0  # u = 0 on bottom boundary

# Iterative solver for Laplace's equation
for iteration in range(max_iter):
    U_old = U.copy()
    U[1:-1, 1:-1] = 0.25 * (U[2:, 1:-1] + U[:-2, 1:-1] + U[1:-1, 2:] + U[1:-1, :-2])
    
    # Check for convergence
    if np.max(np.abs(U - U_old)) < tolerance:
        break

# Plot the solution
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, U, cmap='viridis')
plt.colorbar(label='u(x, y)')
plt.title('Solution to Laplace\'s Equation')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
L = 1.0  # domain size
Nx = 100  # number of spatial points
Ny = 100  # number of spatial points

# Define the grid
x = np.linspace(-L, L, Nx)
y = np.linspace(-L, L, Ny)
X, Y = np.meshgrid(x, y)

# Define complex potential function
def complex_potential(z):
    return np.log(z)  # example: logarithmic potential

# Compute the real and imaginary parts
Z = X + 1j * Y
U = np.real(complex_potential(Z))
V = np.imag(complex_potential(Z))

# Plot the solution
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.contourf(X, Y, U, cmap='viridis')
plt.colorbar(label='Re(u)')
plt.title('Real Part of Complex Potential')
plt.xlabel('x')
plt.ylabel('y')

plt.subplot(1, 2, 2)
plt.contourf(X, Y, V, cmap='viridis')
plt.colorbar(label='Im(u)')
plt.title('Imaginary Part of Complex Potential')
plt.xlabel('x')
plt.ylabel('y')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
L = 10.0  # Length of the domain
Nx = 256  # Number of spatial points
Nt = 500  # Number of time steps
dx = L / Nx
dt = 0.01  # Time step
alpha = 0.01  # Diffusion coefficient

# Define the grid
x = np.linspace(0, L, Nx)
u = np.sin(np.pi * x / L)  # Initial condition
u_hat = np.fft.fft(u)

# Time evolution
for _ in range(Nt):
    k = np.fft.fftfreq(Nx, dx) * 2 * np.pi
    u_hat = u_hat * np.exp(-alpha * (k**2) * dt)
    u = np.fft.ifft(u_hat).real

# Plot the solution
plt.plot(x, u)
plt.title('Heat Equation Solution')
plt.xlabel('x')
plt.ylabel('u(x, t)')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
L = 10.0  # Length of the domain
Nx = 256  # Number of spatial points
Nt = 500  # Number of time steps
dx = L / Nx
dt = 0.01  # Time step
alpha = 0.01  # Diffusion coefficient

# Define the grid
x = np.linspace(0, L, Nx)
u = np.sin(np.pi * x / L)  # Initial condition
u_hat = np.fft.fft(u)

# Time evolution
for _ in range(Nt):
    k = np.fft.fftfreq(Nx, dx) * 2 * np.pi
    u_hat = u_hat * np.exp(-alpha * (k**2) * dt)
    u = np.fft.ifft(u_hat).real

# Plot the solution
plt.plot(x, u)
plt.title('Heat Equation Solution')
plt.xlabel('x')
plt.ylabel('u(x, t)')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Parameters
epsilon = 0.1
x = np.linspace(-1, 1, 1000)
delta = (1 / np.sqrt(np.pi * epsilon)) * np.exp(-x**2 / epsilon)

# Plot the Dirac delta function approximation
plt.figure(figsize=(8, 6))
plt.plot(x, delta, label=r'$\delta(x)$ with $\epsilon=0.1$')
plt.title('Dirac Delta Function Approximation')
plt.xlabel('x')
plt.ylabel(r'$\delta(x)$')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the function to be expanded
def f(x):
    return np.sin(x)  # Example function

# Fourier series expansion
def fourier_series(x, N):
    a0 = (1 / np.pi) * np.trapz(f(np.linspace(0, 2*np.pi, 1000)), np.linspace(0, 2*np.pi, 1000))
    series_sum = a0 / 2
    for n in range(1, N + 1):
        an = (1 / np.pi) * np.trapz(f(np.linspace(0, 2*np.pi, 1000)) * np.cos(n * np.linspace(0, 2*np.pi, 1000)), np.linspace(0, 2*np.pi, 1000))
        bn = (1 / np.pi) * np.trapz(f(np.linspace(0, 2*np.pi, 1000)) * np.sin(n * np.linspace(0, 2*np.pi, 1000)), np.linspace(0, 2*np.pi, 1000))
        series_sum += an * np.cos(n * x) + bn * np.sin(n * x)
    return series_sum

# Define domain
x = np.linspace(0, 2*np.pi, 1000)

# Calculate Fourier series expansion for different N
N = [1, 5, 10]
plt.figure(figsize=(12, 8))
for n in N:
    plt.plot(x, fourier_series(x, n), label=f'N={n}')

# Plot the original function
plt.plot(x, f(x), 'k--', label='Original function')
plt.title('Fourier Series Expansion')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

# Define a function and its derivative
def f(x):
    return np.sin(x) + 0.5 * x

def weak_derivative(x, y):
    spline = UnivariateSpline(x, y, k=1, s=0)
    return spline.derivative()(x)

# Create a range of x values and compute the function and its weak derivative
x = np.linspace(0, 2 * np.pi, 1000)
y = f(x)
y_prime = weak_derivative(x, y)

# Plot the function and its weak derivative
plt.figure(figsize=(12, 6))
plt.plot(x, y, label='Function f(x)')
plt.plot(x, y_prime, label='Weak Derivative f\'(x)', linestyle='--')
plt.title('Function and Its Weak Derivative')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np

# Define a function and its weak derivative
def f(x):
    return np.sin(x) + 0.5 * x

def weak_derivative(x, y):
    spline = UnivariateSpline(x, y, k=1, s=0)
    return spline.derivative()(x)

# Compute Sobolev norm for a function and its derivative
def sobolev_norm(x, y, p=2):
    y_prime = weak_derivative(x, y)
    norm_f = np.sqrt(np.trapz(y**p, x))
    norm_f_prime = np.sqrt(np.trapz(y_prime**p, x))
    return norm_f, norm_f_prime

# Create a range of x values
x = np.linspace(0, 2 * np.pi, 1000)
y = f(x)

# Compute Sobolev norm
norm_f, norm_f_prime = sobolev_norm(x, y)
print(f"Sobolev Norm of f(x): {norm_f}")
print(f"Sobolev Norm of f'(x): {norm_f_prime}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# Define the Green's function for the 1D Laplace operator
def green_function(x, y, L):
    if x < y:
        return (y - x) * (y - L)
    else:
        return (x - y) * (x - L)

# Define the domain
L = 1
x = np.linspace(0, L, 1000)

# Compute Green's function for a fixed point y
y = 0.5
G = np.array([green_function(xi, y, L) for xi in x])

# Plot Green's function
plt.figure(figsize=(8, 6))
plt.plot(x, G, label='Green\'s Function G(x, y)')
plt.title('Green\'s Function for the 1D Laplace Operator')
plt.xlabel('x')
plt.ylabel('G(x, y)')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the eigenfunctions for the 1D Laplace operator with Dirichlet boundary conditions
def eigenfunction(n, x, L):
    return np.sqrt(2 / L) * np.sin(n * np.pi * x / L)

# Define domain and parameters
L = 1
n_values = [1, 2, 3]
x = np.linspace(0, L, 1000)

# Plot eigenfunctions
plt.figure(figsize=(10, 6))
for n in n_values:
    plt.plot(x, eigenfunction(n, x, L), label=f'Eigenfunction n={n}')

plt.title('Eigenfunctions of the 1D Laplace Operator')
plt.xlabel('x')
plt.ylabel('Eigenfunction')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

# Define the matrix A
def matrix_A(n):
    A = np.diag([2] * n) - np.diag([1] * (n-1), k=1) - np.diag([1] * (n-1), k=-1)
    return A

# Compute eigenvalues of the matrix A
def compute_eigenvalues(n):
    A = matrix_A(n)
    eigenvalues, _ = eigh(A)
    return eigenvalues

# Define matrix size
n = 10
eigenvalues = compute_eigenvalues(n)

# Plot eigenvalues
plt.figure(figsize=(8, 6))
plt.plot(range(n), eigenvalues, 'o-', label='Eigenvalues')
plt.title('Eigenvalues of the Matrix A')
plt.xlabel('Index')
plt.ylabel('Eigenvalue')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the perturbed domain function
def perturbed_domain(x, epsilon):
    return x * (1 + epsilon * np.sin(2 * np.pi * x))

# Define the eigenfunction for the perturbed domain
def eigenfunction_perturbed(n, x, epsilon):
    return np.sin(n * np.pi * x) * (1 + epsilon * np.sin(2 * np.pi * x))

# Parameters
epsilon = 0.1
n = 3
x = np.linspace(0, 1, 1000)
eigenfunc = eigenfunction_perturbed(n, x, epsilon)

# Plot eigenfunction for the perturbed domain
plt.figure(figsize=(10, 6))
plt.plot(x, eigenfunc, label=f'Eigenfunction n={n} with perturbation ε={epsilon}')
plt.title('Eigenfunction for the Perturbed Domain')
plt.xlabel('x')
plt.ylabel('Eigenfunction')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the grid size and step size
N = 50  # Number of grid points in each dimension
L = 1.0  # Length of the domain in each dimension
h = L / (N - 1)  # Grid spacing

# Create a grid
x = np.linspace(0, L, N)
y = np.linspace(0, L, N)
X, Y = np.meshgrid(x, y)

# Initialize the solution matrix
U = np.zeros((N, N))

# Set boundary conditions
U[0, :] = 0          # u(x, 0) = 0
U[-1, :] = 0         # u(x, L) = 0
U[:, 0] = np.sin(np.pi * X[:, 0])  # u(0, y) = sin(pi * y)
U[:, -1] = np.sin(np.pi * X[:, -1])  # u(L, y) = sin(pi * y)

# Iterate to solve the PDE using the finite difference method
for _ in range(10000):  # Number of iterations
    U_old = U.copy()
    U[1:-1, 1:-1] = 0.25 * (U_old[1:-1, :-2] + U_old[1:-1, 2:] + U_old[:-2, 1:-1] + U_old[2:, 1:-1])

# Plot the result
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, U, levels=50, cmap='viridis')
plt.colorbar(label='u(x, y)')
plt.title('Numerical Solution of 2D Laplace Equation')
plt.xlabel('x')
plt.ylabel('y')
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
T = 2.0   # total time
Nx = 100   # number of spatial points
Nt = 100   # number of time steps
c = 1.0    # wave speed

# Spatial and time grids
x = np.linspace(0, L, Nx)
t = np.linspace(0, T, Nt)

# Initial conditions
def initial_conditions(x):
    u0 = np.sin(np.pi * x / L)  # Initial displacement
    v0 = np.zeros_like(x)        # Initial velocity
    return u0, v0

u0, v0 = initial_conditions(x)

# Define the PDE system as a first-order system
def pde_system(t, y):
    u = y[:Nx]
    v = y[Nx:]
    
    # Compute second derivatives using central differences
    du2_dx2 = np.zeros(Nx)
    du2_dx2[1:-1] = (u[:-2] - 2 * u[1:-1] + u[2:]) / (x[1] - x[0])**2
    
    # Compute time derivatives
    du_dt = v
    dv_dt = c**2 * du2_dx2
    
    return np.concatenate([du_dt, dv_dt])

# Initial state
y0 = np.concatenate([u0, v0])

# Time span for integration
t_span = (0, T)

# Solve the PDE system
sol = solve_ivp(pde_system, t_span, y0, t_eval=t, method='RK45')

# Extract solution
u = sol.y[:Nx, :]
v = sol.y[Nx:, :]

# Plot results
plt.figure(figsize=(12, 6))

# Plot displacement over time
plt.subplot(2, 1, 1)
plt.imshow(u, extent=[0, T, 0, L], aspect='auto', origin='lower')
plt.colorbar(label='Displacement')
plt.xlabel('Time')
plt.ylabel('Position')
plt.title('Displacement over Time')

# Plot velocity over time
plt.subplot(2, 1, 2)
plt.imshow(v, extent=[0, T, 0, L], aspect='auto', origin='lower')
plt.colorbar(label='Velocity')
plt.xlabel('Time')
plt.ylabel('Position')
plt.title('Velocity over Time')

plt.tight_layout()
plt.show()

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

# Parameters
L = 10.0  # length of the domain
Nx = 100   # number of spatial points
a = 1.0    # coefficient in the differential equation
x = np.linspace(0, L, Nx)

# Define the function f(x)
def f(x):
    return np.sin(np.pi * x / L)  # Example function

# Discretize the domain
dx = L / (Nx - 1)
x_inner = x[1:-1]

# Construct the coefficient matrix A and the right-hand side vector b
main_diag = -2 * np.ones(Nx - 1)
off_diag = np.ones(Nx - 2)
A = diags([main_diag, off_diag, off_diag], [0, -1, 1], format='csr')
A /= dx**2
A += a / (2 * dx) * diags([np.ones(Nx - 2), -np.ones(Nx - 2)], [0, 1], format='csr')

b = f(x_inner)
b = b / dx**2

# Solve the linear system A * u = b
u_inner = spsolve(A, b)
u = np.concatenate(([0], u_inner, [0]))  # Adding boundary conditions

# Plot the result
plt.plot(x, u, label='Numerical Solution')
plt.xlabel('x')
plt.ylabel('u(x)')
plt.title('Solution of the Elliptic PDE')
plt.legend()
plt.grid()
plt.show()

# Check smoothness by computing the second derivative
u_2nd_derivative = np.zeros_like(x)
u_2nd_derivative[1:-1] = (u[:-2] - 2 * u[1:-1] + u[2:]) / dx**2

# Plot second derivative
plt.plot(x, u_2nd_derivative, label='Second Derivative of u(x)')
plt.xlabel('x')
plt.ylabel('u\'\'(x)')
plt.title('Second Derivative of the Solution')
plt.legend()
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 時間領域の信号
def signal(t):
    return np.sin(2 * np.pi * 5 * t)

# フーリエ積分
def fourier_integral(t, f):
    return np.trapz(signal(t) * np.exp(-1j * 2 * np.pi * f * t), t)

# 時間領域
t = np.linspace(-1, 1, 500)
f = np.linspace(-10, 10, 500)

# フーリエ積分
signal_f = np.array([fourier_integral(t, freq) for freq in f])

# プロット
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(t, signal(t))
plt.title('Time Domain Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')

plt.subplot(2, 1, 2)
plt.plot(f, np.abs(signal_f))
plt.title('Fourier Integral')
plt.xlabel('Frequency')
plt.ylabel('Magnitude')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定義された熱伝導方程式の解
def heat_conduction(x, t, alpha):
    return (1 / np.sqrt(4 * np.pi * alpha * t)) * np.exp(-x**2 / (4 * alpha * t))

# パラメータ
alpha = 1
x = np.linspace(-10, 10, 500)
t = [0.1, 1, 5]  # 時間の異なる値

# プロット
plt.figure(figsize=(10, 6))
for ti in t:
    plt.plot(x, heat_conduction(x, ti, alpha), label=f't = {ti}')

plt.xlabel('x')
plt.ylabel('Temperature')
plt.title('Heat Conduction in an Infinite Rod')
plt.legend()
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 3次元熱伝導のグリーン関数
def green_function_3d(x, y, z, t, alpha):
    r = np.sqrt(x**2 + y**2 + z**2)
    return (1 / ( (4 * np.pi * alpha * t)**(3/2) )) * np.exp(-r**2 / (4 * alpha * t))

# パラメータ
alpha = 1
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)
Z = green_function_3d(X, Y, 0, 1, alpha)

# プロット
plt.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar(label='Temperature')
plt.title('3D Heat Conduction Green Function (z=0, t=1)')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# ディリクレ・ノイマン問題の解法（円の領域）
def dirichlet_solution(x, y):
    return np.log(np.sqrt(x**2 + y**2))

# グリッドの設定
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x, y)
Z = dirichlet_solution(X, Y)

# プロット
plt.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar(label='Potential')
plt.title('Dirichlet Solution in Circular Region')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# グリーンの定理を用いた解法
def green_theorem_solution(x, y):
    return np.sin(np.pi * x) * np.sin(np.pi * y)

# グリッドの設定
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)
Z = green_theorem_solution(X, Y)

# プロット
plt.contourf(X, Y, Z, levels=50, cmap='viridis')
plt.colorbar(label='Potential')
plt.title('Solution Using Green\'s Theorem')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
N = 10  # フーリエ級数の項数
r = np.linspace(0, 1, 100)  # 半径の範囲
theta = np.linspace(0, 2 * np.pi, 100)  # 角度の範囲
R, THETA = np.meshgrid(r, theta)

# 境界条件の設定
def boundary_condition(theta):
    return np.sin(theta)  # ここでは sin(θ) を例にしています

# フーリエ級数による解の計算
def laplace_solution(r, theta, N):
    solution = np.zeros_like(r)
    for n in range(1, N + 1):
        A_n = 1 / np.pi * np.trapz(boundary_condition(theta) * np.sin(n * theta), theta)  # フーリエ係数
        solution += A_n * r**n  # フーリエ級数の和
    return solution

# 解の計算
u = laplace_solution(R, THETA, N)

# プロット
plt.figure(figsize=(10, 6))
plt.contourf(R * np.cos(THETA), R * np.sin(THETA), u, levels=50, cmap='viridis')
plt.colorbar(label='Potential')
plt.title('Laplace Equation Solution on Unit Disk (Fourier Series)')
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# グリーン関数の定義
def green_function_2d(x, y, xi, eta):
    r = np.sqrt((x - xi)**2 + (y - eta)**2)
    return -1 / (2 * np.pi) * np.log(r)

# パラメータ設定
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)
xi, eta = 0.5, 0.5  # グリーン関数の評価点

# グリーン関数の計算
G = green_function_2d(X, Y, xi, eta)

# プロット
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, G, levels=50, cmap='viridis')
plt.colorbar(label='G(x, y, xi, eta)')
plt.title('2D Green Function')
plt.xlabel('x')
plt.ylabel('y')
plt.show()