理工系数学と計算科学の基礎

https://www.amazon.co.jp/%E7%90%86%E5%B7%A5%E7%B3%BB%E6%95%B0%E5%AD%A6%E3%81%A8%E8%A8%88%E7%AE%97%E7%A7%91%E5%AD%A6%E3%81%AE%E5%9F%BA%E7%A4%8E-%E6%9D%BE%E6%9C%AC-%E8%81%A1/dp/4485301214
P2

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# 1. サンプルデータの生成
np.random.seed(0)  # 再現性のためのシード値
mu = 50   # 母平均
sigma = 5 # 母標準偏差
sample_size = 30  # 標本サイズ
samples = np.random.normal(mu, sigma, sample_size)  # 正規分布に従う標本データを生成

# 2. 標本平均と標本分散の計算
sample_mean = np.mean(samples)
sample_variance = np.var(samples, ddof=1)  # 不偏分散 (ddof=1 は分母が n-1)

# 3. 母平均と標本平均の差を t 検定で p 値を計算
t_statistic, p_value = stats.ttest_1samp(samples, mu)

# 4. 外れ値の検出
# 外れ値は標準正規分布の3σ範囲外 (±3σ) と定義
outliers = samples[np.abs(samples - sample_mean) > 3 * np.std(samples)]

# 5. 標準不確かさの計算
# 標準不確かさ u は、標本標準偏差 / √(標本サイズ)
standard_uncertainty = np.std(samples, ddof=1) / np.sqrt(sample_size)

# 6. 結果の表示
print(f"サンプルデータ: {samples}")
print(f"標本平均: {sample_mean:.2f}")
print(f"標本分散: {sample_variance:.2f}")
print(f"p値: {p_value:.5f}")
print(f"外れ値: {outliers}")
print(f"標準不確かさ (u): {standard_uncertainty:.2f}")

# 7. データのヒストグラムと標準正規分布のプロット
plt.figure(figsize=(10, 6))
count, bins, ignored = plt.hist(samples, bins=10, density=True, alpha=0.6, color='g', label='Sample Data')

# 標準正規分布のプロット
x = np.linspace(min(bins), max(bins), 100)
pdf = stats.norm.pdf(x, mu, sigma)
plt.plot(x, pdf, 'k', linewidth=2, label='Normal Distribution')

# 外れ値のプロット
if len(outliers) > 0:
    plt.scatter(outliers, np.zeros(len(outliers)), color='red', s=100, zorder=5, label='Outliers')

plt.title('Sample Data vs. Normal Distribution')
plt.xlabel('Value')
plt.ylabel('Density')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# サンプルデータを生成
np.random.seed(42)  # 再現性のためにシードを設定
n = 30  # サンプル数
data = np.random.normal(loc=50, scale=10, size=n)  # 平均50、標準偏差10の正規分布

# 平均値の計算
mean = np.mean(data)

# 標本分散の計算（標本分散 u^2 の計算）
sample_variance = np.var(data, ddof=1)  # ddof=1は標本分散のための設定

# 母分散の計算（sigma^2 の計算）
population_variance = np.var(data, ddof=0)  # ddof=0は母分散のための設定

# 標準偏差の計算
sample_std_dev = np.sqrt(sample_variance)
population_std_dev = np.sqrt(population_variance)

# 結果を表示
print(f"サンプル平均 (mean): {mean:.2f}")
print(f"標本分散 (sample variance): {sample_variance:.2f}")
print(f"母分散 (population variance): {population_variance:.2f}")
print(f"標本標準偏差 (sample standard deviation): {sample_std_dev:.2f}")
print(f"母標準偏差 (population standard deviation): {population_std_dev:.2f}")

# ヒストグラムをプロット
plt.hist(data, bins=10, color='lightblue', edgecolor='black', alpha=0.7)
plt.axvline(mean, color='red', linestyle='dashed', linewidth=2, label=f'Mean = {mean:.2f}')
plt.title('Sample Data Distribution')
plt.xlabel('Data Values')
plt.ylabel('Frequency')
plt.legend()
plt.grid(True)
plt.show()

# データの分布をプロット
plt.figure(figsize=(10, 6))
x = np.linspace(min(data), max(data), 100)
y = (1 / (population_std_dev * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((x - mean) / population_std_dev) ** 2)
plt.plot(x, y, label='Normal Distribution (Population)', color='orange')
plt.hist(data, bins=10, density=True, alpha=0.6, color='skyblue', edgecolor='black')
plt.axvline(mean, color='red', linestyle='dashed', linewidth=2, label=f'Mean = {mean:.2f}')
plt.title('Sample Data Distribution and Normal Distribution')
plt.xlabel('Data Values')
plt.ylabel('Density')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# x 軸の値を設定 (-5 から 5 までの範囲)
x = np.linspace(-5, 5, 500)

# 確率密度関数 (PDF) を計算
pdf = norm.pdf(x, 0, 1)  # 平均 0, 標準偏差 1 の標準正規分布
# 累積分布関数 (CDF) を計算
cdf = norm.cdf(x, 0, 1)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x, pdf, label='Probability Density Function (PDF)', color='blue')
plt.plot(x, cdf, label='Cumulative Distribution Function (CDF)', color='orange')
plt.axvline(x=1, color='gray', linestyle='--', linewidth=1)
plt.axvline(x=2, color='gray', linestyle='--', linewidth=1)
plt.axvline(x=3, color='gray', linestyle='--', linewidth=1)
plt.axvline(x=-1, color='gray', linestyle='--', linewidth=1)
plt.axvline(x=-2, color='gray', linestyle='--', linewidth=1)
plt.axvline(x=-3, color='gray', linestyle='--', linewidth=1)
plt.text(1, 0.1, r'$\pm1\sigma\ (68.2\%)$', fontsize=12, verticalalignment='bottom')
plt.text(2, 0.05, r'$\pm2\sigma\ (95.4\%)$', fontsize=12, verticalalignment='bottom')
plt.text(3, 0.02, r'$\pm3\sigma\ (99.7\%)$', fontsize=12, verticalalignment='bottom')
plt.xlabel('z')
plt.ylabel('Density / Probability')
plt.title('Standard Normal Distribution and Cumulative Distribution Function')
plt.legend()
plt.grid(True)
plt.show()

P7

import numpy as np
from scipy.optimize import minimize

# Define the function A with respect to x and t
def A(x, t):
    # Example function, modify this based on the actual problem
    return x[0] * np.exp(-x[1] * t) + x[2]

# Define the partial derivatives of A with respect to each x
def partial_derivative_A(x, t):
    dA_dx = np.zeros(len(x))
    dA_dx[0] = np.exp(-x[1] * t)  # dA/dx1
    dA_dx[1] = -x[0] * t * np.exp(-x[1] * t)  # dA/dx2
    dA_dx[2] = 1  # dA/dx3
    return dA_dx

# Define the function to minimize, which is the sum of squared differences S
def sum_of_squares(x, t_values, measured_A):
    S = 0
    for i, t in enumerate(t_values):
        delta_A_i = measured_A[i] - A(x, t)
        partials = partial_derivative_A(x, t)
        # Calculate the sum of squared differences
        S += (delta_A_i - np.dot(partials, x))**2
    return S

# Example data: measured values of A and time points t
t_values = np.linspace(0, 10, 50)  # Time points from 0 to 10
measured_A = 3.0 * np.exp(-0.5 * t_values) + 2.0  # Example measurements

# Initial guess for the variables x (e.g., coefficients)
initial_guess = np.array([1.0, 0.1, 1.0])

# Minimize the sum of squares
result = minimize(sum_of_squares, initial_guess, args=(t_values, measured_A))

# Display the optimized parameters
print("Optimized parameters:", result.x)
print("Sum of squares:", result.fun)

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

# 1. 記述統計量の計算
# サンプルデータを生成
x_data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])  # 入力データ (例)
y_data = np.array([2.1, 2.9, 3.7, 4.1, 5.2, 5.9, 7.1, 8.1, 8.8, 10.2])  # 出力データ (例)

# 平均値の計算
mean_x = np.mean(x_data)
mean_y = np.mean(y_data)

# 標本分散の計算
variance_x = np.var(x_data, ddof=1)
variance_y = np.var(y_data, ddof=1)

# 標本標準偏差の計算
std_dev_x = np.sqrt(variance_x)
std_dev_y = np.sqrt(variance_y)

# 2. 最小二乗法の適用
# 線形回帰モデルの定義
def linear_model(x, a, b):
    return a * x + b

# 最小二乗法によるフィッティング
params, covariance = curve_fit(linear_model, x_data, y_data)
a, b = params

# 近似値の計算
y_fit = linear_model(x_data, a, b)

# 3. 相関係数の計算
correlation_coefficient = np.corrcoef(x_data, y_data)[0, 1]

# 4. 結果の表示
print(f"平均 (Mean) of x: {mean_x:.2f}, y: {mean_y:.2f}")
print(f"標本分散 (Sample Variance) of x: {variance_x:.2f}, y: {variance_y:.2f}")
print(f"標本標準偏差 (Sample Standard Deviation) of x: {std_dev_x:.2f}, y: {std_dev_y:.2f}")
print(f"線形回帰の係数: a = {a:.2f}, b = {b:.2f}")
print(f"相関係数 (Correlation Coefficient): {correlation_coefficient:.2f}")

# 5. プロット
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, color='blue', label='Data Points')
plt.plot(x_data, y_fit, color='red', label=f'Fitted Line: y = {a:.2f}x + {b:.2f}')
plt.axhline(mean_y, color='orange', linestyle='--', label=f'Mean y = {mean_y:.2f}')
plt.axvline(mean_x, color='green', linestyle='--', label=f'Mean x = {mean_x:.2f}')
plt.title('Least Squares Fit and Correlation')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

def calculate_measurement_error(full_scale, error_percentage, measured_value):
    """
    測定誤差を計算する関数
    
    Parameters:
        full_scale (float): 計測器のフルスケール値 (例: 100V)
        error_percentage (float): フルスケールの誤差率 (%) (例: 1%)
        measured_value (float): 実際の測定値 (例: 40V)
    
    Returns:
        max_absolute_error (float): 最大絶対誤差 (例: 1V)
        relative_error (float): 測定値に対する相対誤差 (%) (例: 2.5%)
    """
    # フルスケールの最大絶対誤差の計算
    max_absolute_error = full_scale * (error_percentage / 100)
    
    # 測定値に対する相対誤差の計算
    relative_error = (max_absolute_error / measured_value) * 100
    
    return max_absolute_error, relative_error

# 例題のパラメータを設定
full_scale = 100          # フルスケール値 (V)
error_percentage = 1      # フルスケール誤差率 (%)
measured_value = 40       # 測定値 (V)

# 誤差の計算
max_error, relative_error = calculate_measurement_error(full_scale, error_percentage, measured_value)

# 結果の表示
print(f"フルスケールの最大絶対誤差: {max_error:.2f} V")
print(f"測定値に対する相対誤差: {relative_error:.2f} %")

P9

import numpy as np

# Define the function g with respect to variables x
def g(x):
    # Example function g(x1, x2, ..., xm)
    return x[0]**2 + 2*x[1] + 0.5*x[2]  # Modify based on your specific function

# Define the partial derivatives of g with respect to each x
def partial_derivatives_g(x):
    dgd_x = np.zeros(len(x))
    dgd_x[0] = 2 * x[0]  # Partial derivative of g with respect to x1
    dgd_x[1] = 2  # Partial derivative of g with respect to x2
    dgd_x[2] = 0.5  # Partial derivative of g with respect to x3
    return dgd_x

# Example values of x and their measurement errors
x_values = np.array([1.0, 2.0, 3.0])  # x1, x2, x3
delta_x = np.array([0.05, 0.02, 0.01])  # Measurement errors for each x
E_values = np.array([0.1, 0.05, 0.02])  # Precision values for each x

# Calculate δA using the formula (1.9.2)
dA = np.dot(partial_derivatives_g(x_values), delta_x)

# Calculate ε using the root sum square formula (1.9.3)
epsilon = np.sqrt(np.sum((partial_derivatives_g(x_values) * delta_x)**2))

# Calculate E using formula (1.9.4)
E = np.sqrt(np.sum((partial_derivatives_g(x_values) * E_values)**2))

# Display results
print(f"Total Error δA: {dA}")
print(f"Cumulative Error ε: {epsilon}")
print(f"Overall Precision E: {E}")

P13

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
eta = 1  # 尺度パラメータ
t = np.linspace(0, 5, 500)  # 時間 t の範囲

# 形状パラメータ m のリスト
m_values = [0.5, 1, 2, 5]

# 累積故障率 F(t) の計算とプロット
plt.figure(figsize=(10, 6))
for m in m_values:
    # 累積故障率 F(t) の計算
    F_t = 1 - np.exp(-(t / eta)**m)
    
    # グラフのプロット
    plt.plot(t, F_t, label=f'm = {m}')

# グラフの装飾
plt.title('Weibull Distribution: Cumulative Failure Rate F(t)')
plt.xlabel('Time t')
plt.ylabel('Cumulative Failure Rate F(t)')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import curve_fit

# 乱数生成の設定
np.random.seed(0)  # 乱数のシード値を設定 (再現性のため)
sample_size = 1000  # 生成するデータ数

# 正規分布に従うデータを生成
mu, sigma = 0, 1  # 平均値と標準偏差
normal_data = np.random.normal(mu, sigma, sample_size)

# ワイブル分布に従うデータを生成
shape_param, scale_param = 2, 1  # ワイブル分布の形状パラメータと尺度パラメータ
weibull_data = np.random.weibull(shape_param, sample_size) * scale_param

# ワイブル分布の累積分布関数 (CDF)
def weibull_cdf(x, shape, scale):
    return 1 - np.exp(-(x / scale) ** shape)

# パラメータ推定のためにフィッティングを行う
params, _ = curve_fit(weibull_cdf, np.sort(weibull_data), np.linspace(0, 1, sample_size))

# パラメータの推定結果
estimated_shape, estimated_scale = params
print(f"推定された形状パラメータ (m): {estimated_shape:.2f}")
print(f"推定された尺度パラメータ (eta): {estimated_scale:.2f}")

# ワイブル分布と正規分布のヒストグラムの表示
plt.figure(figsize=(12, 6))

# 正規分布のプロット
plt.subplot(1, 2, 1)
plt.hist(normal_data, bins=30, density=True, alpha=0.6, color='g')
plt.title('Normal Distribution')
plt.xlabel('Value')
plt.ylabel('Frequency')

# ワイブル分布のプロット
plt.subplot(1, 2, 2)
plt.hist(weibull_data, bins=30, density=True, alpha=0.6, color='b')
x = np.linspace(0, np.max(weibull_data), 100)
plt.plot(x, weibull_cdf(x, estimated_shape, estimated_scale), 'r--', label='Estimated Weibull CDF')
plt.title('Weibull Distribution')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.legend()

plt.tight_layout()
plt.show()

P19

# Function to propagate errors for addition
def error_addition(delta_a, delta_b):
    return delta_a + delta_b

# Function to propagate errors for multiplication
def error_multiplication(a, delta_a, b, delta_b):
    return abs(delta_a / a) + abs(delta_b / b)

# Function to propagate errors for power rule
def error_power(a, delta_a, n):
    return n * abs(delta_a / a)

# Example usage
a = 10
delta_a = 0.5
b = 20
delta_b = 1

# Error in addition
delta_c_add = error_addition(delta_a, delta_b)
print(f"Error in addition: {delta_c_add}")

# Error in multiplication
delta_c_mult = error_multiplication(a, delta_a, b, delta_b)
print(f"Relative error in multiplication: {delta_c_mult}")

# Error in power
n = 2
delta_c_power = error_power(a, delta_a, n)
print(f"Relative error in power: {delta_c_power}")

import numpy as np
import matplotlib.pyplot as plt

# Parameters for the Weibull distribution
m = 1.5  # shape parameter
n = 2.0  # shape parameter
V_L = 0.5  # location parameter
V_e = 1.0  # scale parameter
t_0 = 1.0  # reference time

# Function for Weibull cumulative distribution
def weibull_cdf(V, t, m, n, V_L, V_e, t_0):
    return 1 - np.exp(-((V - V_L) / V_e)**m * (t / t_0)**n)

# Values for V and t
V_values = np.linspace(0, 2, 100)  # V values from 0 to 2
t_values = np.linspace(0, 2, 100)  # t values from 0 to 2

# Create a meshgrid for V and t values
V, t = np.meshgrid(V_values, t_values)

# Calculate Weibull CDF
F = weibull_cdf(V, t, m, n, V_L, V_e, t_0)

# Plotting
plt.figure(figsize=(8, 6))
cp = plt.contourf(V, t, F, cmap='viridis')
plt.colorbar(cp)
plt.title("Weibull Cumulative Distribution Function")
plt.xlabel("V")
plt.ylabel("t")
plt.show()

P28

import numpy as np
import matplotlib.pyplot as plt

# Constants
k = 1.38e-23  # Boltzmann constant (J/K)
T = 300  # Temperature in Kelvin
m = 4.002602e-3 / 6.022e23  # Mass of helium in kg (for example)

# Maxwell-Boltzmann distribution function
def maxwell_boltzmann(v, m, T):
    factor = 4 / np.sqrt(np.pi) * (m / (2 * k * T))**(3/2)
    return factor * v**2 * np.exp(-m * v**2 / (2 * k * T))

# Velocity range
v = np.linspace(0, 2500, 500)

# Compute distribution for given mass and temperature
f_v = maxwell_boltzmann(v, m, T)

# Plot the distribution
plt.figure(figsize=(8, 6))
plt.plot(v, f_v, label='Maxwell-Boltzmann Distribution')
plt.title('Maxwell-Boltzmann Velocity Distribution')
plt.xlabel('Velocity (m/s)')
plt.ylabel('Probability Density')
plt.legend()
plt.grid(True)
plt.show()

P39

import cmath

# パラメータを設定（例として任意の値を使用）
R = 10  # 単位長さあたりの抵抗 (Ω/m)
L = 0.001  # 単位長さあたりのインダクタンス (H/m)
G = 0.0001  # 単位長さあたりのコンダクタンス (S/m)
C = 0.00001  # 単位長さあたりのキャパシタンス (F/m)
omega = 2 * cmath.pi * 50  # 角周波数 (rad/s), ここでは50 Hzを仮定

# 単位長さあたりのインピーダンスとアドミタンスの定義
z = R + 1j * omega * L
y = G + 1j * omega * C

# 伝搬定数の計算
gamma = cmath.sqrt(z * y)

# 特性インピーダンスの計算
Z0 = cmath.sqrt(z / y)

# 計算結果の表示
gamma_real = gamma.real
gamma_imag = gamma.imag
Z0_real = Z0.real
Z0_imag = Z0.imag

gamma_real, gamma_imag, Z0_real, Z0_imag

P40

import cmath

# 伝送線のパラメータを設定
Z0 = 50  # 特性インピーダンス (Ω), 例として50Ωを使用
Z = 10 + 1j * 2  # インピーダンス (Ω), 任意の値
Y = 0.01 + 1j * 0.005  # アドミタンス (S), 任意の値
l = 100  # 伝送線の長さ (m)

# 伝搬定数の計算
gamma = cmath.sqrt(Z * Y)

# 四端子パラメータの計算
A = cmath.cosh(gamma * l)
B = Z0 * cmath.sinh(gamma * l)
C = (1 / Z0) * cmath.sinh(gamma * l)
D = A  # 対称回路の場合

# 結果の表示
A_real, A_imag = A.real, A.imag
B_real, B_imag = B.real, B.imag
C_real, C_imag = C.real, C.imag
D_real, D_imag = D.real, D.imag

# 結果を表示
A_real, A_imag, B_real, B_imag, C_real, C_imag, D_real, D_imag

P41

import numpy as np

# 相電圧を定義 (例として任意の値を設定)
E_a = 100 + 0j   # A相電圧 (複素数)
E_b = -50 + 86.6j  # B相電圧 (複素数)
E_c = -50 - 86.6j  # C相電圧 (複素数)

# ベクトル演算子 a の定義
a = np.exp(1j * 2 * np.pi / 3)  # a = e^(j*2π/3)

# 電圧ベクトルの定義
E_abc = np.array([E_a, E_b, E_c])

# 変換行列 A の定義
A = (1/3) * np.array([[1, 1, 1],
                      [1, a, a**2],
                      [1, a**2, a]])

# 逆相、零相、正相電圧の計算
E_012 = np.dot(A, E_abc)

# 結果の表示
E_0, E_1, E_2 = E_012
E_0, E_1, E_2

P107
LDU分解 LU分解 LU分解 ユニタリ行列 随伴行列

import numpy as np
from scipy.linalg import lu, lu_factor, lu_solve

# 1. LU分解 (LU Decomposition)
def lu_decomposition(A):
    P, L, U = lu(A)  # P: permutation matrix, L: lower triangular, U: upper triangular
    return P, L, U

# 2. LDU分解 (LDU Decomposition)
def ldu_decomposition(A):
    # LU分解をまず行う
    P, L, U = lu(A)
    
    # 上三角行列 U から対角成分 D を抽出
    D = np.diag(np.diag(U))
    
    # 対角成分 D を用いて U を正規化
    U_normalized = np.dot(np.linalg.inv(D), U)
    
    return P, L, D, U_normalized

# 3. ユニタリ行列の作成 (Unitary Matrix)
def create_unitary_matrix(n):
    # ランダムなユニタリ行列を作成 (ユニタリ行列の定義を満たす行列を生成)
    Q, R = np.linalg.qr(np.random.randn(n, n) + 1j*np.random.randn(n, n))  # QR分解を利用してユニタリ行列を生成
    return Q

# 4. 随伴行列 (Adjoint Matrix)
def adjoint_matrix(A):
    # 随伴行列は複素共役転置
    return np.conjugate(A.T)

# === デモ用の行列を定義 ===
A = np.array([[4, 3], [6, 3]])

# 1. LU分解の実行
P, L, U = lu_decomposition(A)
print("LU Decomposition:")
print("Permutation matrix P:\n", P)
print("Lower triangular matrix L:\n", L)
print("Upper triangular matrix U:\n", U)

# 2. LDU分解の実行
P, L, D, U_normalized = ldu_decomposition(A)
print("\nLDU Decomposition:")
print("Permutation matrix P:\n", P)
print("Lower triangular matrix L:\n", L)
print("Diagonal matrix D:\n", D)
print("Upper triangular matrix U (normalized):\n", U_normalized)

# 3. ユニタリ行列の作成
unitary_matrix = create_unitary_matrix(2)
print("\nUnitary Matrix:\n", unitary_matrix)

# 4. 随伴行列の作成
adjoint = adjoint_matrix(unitary_matrix)
print("\nAdjoint Matrix:\n", adjoint)

P60

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

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

# 第1種ベッセル関数 J_n(x) の例
n_values = [0, 1, 2]  # 次数 n の設定
for n in n_values:
    plt.plot(x, jn(n, x), label=f'$J_{n}(x)$')

# 第2種ベッセル関数 Y_n(x) の例
for n in n_values:
    plt.plot(x, yn(n, x), '--', label=f'$Y_{n}(x)$')

plt.xlabel('$x$')
plt.ylabel('$J_n(x), Y_n(x)$')
plt.title('Bessel Functions of the First and Second Kind')
plt.legend()
plt.grid()
plt.show()

from scipy.special import gamma

# xの範囲を設定
x = np.linspace(0.1, 5, 500)  # ガンマ関数は x=0 で特異性を持つため、0.1 からスタート

# ガンマ関数のプロット
plt.plot(x, gamma(x), label='$\Gamma(x)$')
plt.xlabel('$x$')
plt.ylabel('$\Gamma(x)$')
plt.title('Gamma Function')
plt.legend()
plt.grid()
plt.show()

P121

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

# 方程式の定義
def equations(variables):
    x, y = variables
    eq1 = x**2 - 2*x + y**2
    eq2 = np.cos(x) - y
    return [eq1, eq2]

# 初期値設定
initial_guess = [1.0, 1.0]

# fsolveを用いて非線形方程式を解く
solution = fsolve(equations, initial_guess)

# 結果の表示
x_sol, y_sol = solution
print(f"解: x = {x_sol}, y = {y_sol}")

# グラフ描画
x = np.linspace(-3, 3, 400)
y1 = np.sqrt(2*x - x**2)  # 式1の y に対する解
y2 = -np.sqrt(2*x - x**2)  # 式1の y に対する別解
y3 = np.cos(x)  # 式2の解

# 式1のプロット
plt.plot(x, y1, label=r"$x^2 - 2x + y^2 = 0$")
plt.plot(x, y2, label=r"$x^2 - 2x + y^2 = 0$")
# 式2のプロット
plt.plot(x, y3, label=r"$\cos(x) - y = 0$")

# 求めた解をプロット
plt.scatter(x_sol, y_sol, color='red', zorder=5, label=f"Solution: ({x_sol:.2f}, {y_sol:.2f})")

plt.xlabel("x")
plt.ylabel("y")
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.legend()
plt.title("Solving Nonlinear Equations using Newton-Raphson Method")
plt.grid()
plt.show()

P130

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

# 1. チェビシェフ補間 (Chebyshev Interpolation)
def chebyshev_nodes(n, a=-1, b=1):
    """
    チェビシェフ補間点を生成する関数。
    n: 補間点の数
    a: 範囲の下限
    b: 範囲の上限
    """
    i = np.arange(n)
    x_chebyshev = np.cos((2*i + 1) / (2*n) * np.pi)  # [-1, 1] の範囲のチェビシェフ点
    # [-1, 1] のチェビシェフ点を [a, b] の範囲に変換
    return 0.5 * (a + b) + 0.5 * (b - a) * x_chebyshev

# チェビシェフ補間用のサンプル関数
def sample_function(x):
    return 1 / (1 + 25 * x**2)  # ランゲ関数 (Runge's function)

# 補間点の数
n = 10

# チェビシェフ点を生成し、対応する関数値を求める
x_cheb = chebyshev_nodes(n, a=-1, b=1)
y_cheb = sample_function(x_cheb)

# 補間用の多項式を生成 (numpyの多項式フィッティング関数を使用)
coefficients = np.polyfit(x_cheb, y_cheb, deg=n-1)  # チェビシェフ点での多項式フィッティング
polynomial_cheb = np.poly1d(coefficients)  # 多項式を生成

# 連続プロット用のx点を生成
x_dense = np.linspace(-1, 1, 1000)
y_dense = sample_function(x_dense)
y_cheb_interp = polynomial_cheb(x_dense)  # チェビシェフ補間の結果

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x_dense, y_dense, 'k-', label='Original Function (Runge)')
plt.plot(x_cheb, y_cheb, 'ro', label='Chebyshev Nodes')
plt.plot(x_dense, y_cheb_interp, 'b--', label='Chebyshev Interpolation')
plt.legend()
plt.title('Chebyshev Interpolation')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.show()

# 2. スプライン補間 (Spline Interpolation)
# 補間点を均等に生成
x_spline = np.linspace(-1, 1, n)
y_spline = sample_function(x_spline)

# スプライン補間を実行
cubic_spline = CubicSpline(x_spline, y_spline)

# スプライン補間結果を生成
y_spline_interp = cubic_spline(x_dense)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x_dense, y_dense, 'k-', label='Original Function (Runge)')
plt.plot(x_spline, y_spline, 'go', label='Spline Nodes')
plt.plot(x_dense, y_spline_interp, 'm--', label='Cubic Spline Interpolation')
plt.legend()
plt.title('Cubic Spline Interpolation')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.show()

P180

import numpy as np

# Laplacian operator for 2D grid
def laplacian(u, dx, dy):
    return (np.roll(u, 1, axis=0) + np.roll(u, -1, axis=0) - 2*u) / dx**2 + \
           (np.roll(u, 1, axis=1) + np.roll(u, -1, axis=1) - 2*u) / dy**2

# Poisson 方程式: ∇²u(r) = p(r)
def poisson_equation(u, p, dx, dy):
    return laplacian(u, dx, dy) - p

# Helmholtz 方程式: ∇²u + λu = 0
def helmholtz_equation(u, lambd, dx, dy):
    return laplacian(u, dx, dy) + lambd * u

# Schrödinger 方程式: ∇²u + (E - V)u = 0
def schrodinger_equation(u, E, V, dx, dy):
    return laplacian(u, dx, dy) + (E - V) * u

# 振動波: ∂²u/∂t² = c²∇²u
def oscillation_wave(u, c, dx, dy, dt):
    laplacian_u = laplacian(u, dx, dy)
    return c**2 * laplacian_u

# 減衰振動波: ∂²u/∂t² = c²∇²u - h∂u/∂t
def damped_oscillation_wave(u, v, c, h, dx, dy, dt):
    laplacian_u = laplacian(u, dx, dy)
    return c**2 * laplacian_u - h * v

# 電信方程式: ∂²u/∂t² = c²∇²u - ku
def telegraph_equation(u, v, c, k, dx, dy, dt):
    laplacian_u = laplacian(u, dx, dy)
    return c**2 * laplacian_u - k * u

# 強制振動のある波動方程式: ∂²u/∂t² = c²∇²u + f(r,t)
def forced_oscillation_wave(u, v, c, f, dx, dy, dt):
    laplacian_u = laplacian(u, dx, dy)
    return c**2 * laplacian_u + f

# 拡散方程式: ∂u/∂t = a²∇²u
def diffusion_equation(u, a, dx, dy):
    return a**2 * laplacian(u, dx, dy)

# 側面損失のある拡散方程式: ∂u/∂t = a²∇²u - hu
def diffusion_with_loss(u, a, h, dx, dy):
    return a**2 * laplacian(u, dx, dy) - h * u

# 熱源(熱損失)を伴う拡散方程式: ∂u/∂t = a²∇²u + f(r, t)
def heat_diffusion(u, a, f, dx, dy):
    return a**2 * laplacian(u, dx, dy) + f

# テスト用の初期設定
dx, dy, dt = 0.1, 0.1, 0.01
lambd, E, V, c, h, k, a = 1.0, 1.0, 0.5, 1.0, 0.1, 0.5, 1.0
p = np.ones((100, 100))  # Poisson 方程式用のソース項
u = np.zeros((100, 100))  # 初期状態
v = np.zeros((100, 100))  # 初期速度（波動方程式用）
f = np.ones((100, 100))  # 強制項や熱源項

# 方程式の計算例
poisson_result = poisson_equation(u, p, dx, dy)
helmholtz_result = helmholtz_equation(u, lambd, dx, dy)
schrodinger_result = schrodinger_equation(u, E, V, dx, dy)
oscillation_result = oscillation_wave(u, c, dx, dy, dt)
damped_result = damped_oscillation_wave(u, v, c, h, dx, dy, dt)
telegraph_result = telegraph_equation(u, v, c, k, dx, dy, dt)
forced_oscillation_result = forced_oscillation_wave(u, v, c, f, dx, dy, dt)
diffusion_result = diffusion_equation(u, a, dx, dy)
diffusion_loss_result = diffusion_with_loss(u, a, h, dx, dy)
heat_diffusion_result = heat_diffusion(u, a, f, dx, dy)

P181

import numpy as np

def forward_difference_2nd_derivative(u, h_x):
    """
    2次偏導関数の前進差分近似を求める関数
    Args:
    u (numpy array): 格子点における関数 u の値
    h_x (float): x 方向の格子幅

    Returns:
    numpy array: 2次偏導関数の前進差分近似
    """
    # 格子点の個数
    n = len(u)
    
    # 2次偏導関数の格子点における前進差分近似
    d2u_dx2 = np.zeros_like(u)
    
    for i in range(n - 2):
        # 前進差分近似の式を使用して計算
        d2u_dx2[i] = (u[i + 2] - 2 * u[i + 1] + u[i]) / h_x**2
    
    return d2u_dx2

# 格子点における関数 u の例 (例えば、u = sin(x))
x = np.linspace(0, np.pi, 10)  # x 座標を等間隔に生成
u = np.sin(x)                  # 関数値 u を計算
h_x = x[1] - x[0]              # 格子幅を計算

# 2次偏導関数の前進差分近似を計算
d2u_dx2_approx = forward_difference_2nd_derivative(u, h_x)

# 結果を表示
print("2次偏導関数の前進差分近似:")
print(d2u_dx2_approx)

P183

import numpy as np

def finite_difference_approximations(u, hx, hy):
    """
    Calculates finite difference approximations of first and second derivatives.

    Parameters:
    - u: 2D numpy array representing the grid of function values.
    - hx: step size in the x-direction.
    - hy: step size in the y-direction.

    Returns:
    - derivative_x_backward: Backward difference approximation in x-direction.
    - derivative_x_forward: Forward difference approximation in x-direction.
    - derivative_x_central: Central difference approximation in x-direction.
    - second_derivative_x: Second-order derivative in x-direction.
    - second_derivative_y: Second-order derivative in y-direction.
    """
    # Backward difference approximation (partial derivative with respect to x)
    derivative_x_backward = (u[1:, :] - u[:-1, :]) / hx

    # Forward difference approximation (partial derivative with respect to x)
    derivative_x_forward = (u[2:, :] - u[1:-1, :]) / hx

    # Central difference approximation (partial derivative with respect to x)
    derivative_x_central = (u[2:, :] - u[:-2, :]) / (2 * hx)

    # Second-order derivative approximation (with respect to x)
    second_derivative_x = (u[2:, :] - 2 * u[1:-1, :] + u[:-2, :]) / hx**2

    # Second-order derivative approximation (with respect to y)
    second_derivative_y = (u[:, 2:] - 2 * u[:, 1:-1] + u[:, :-2]) / hy**2

    return derivative_x_backward, derivative_x_forward, derivative_x_central, second_derivative_x, second_derivative_y

# Example usage:
# Create a sample grid
u = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])
hx = 1.0  # Step size in x-direction
hy = 1.0  # Step size in y-direction

# Calculate finite difference approximations
results = finite_difference_approximations(u, hx, hy)

for result in results:
    print(result)

import numpy as np
import matplotlib.pyplot as plt

# Parameters for the 2D grid
Lx, Ly = 1.4, 1.0  # Length in x and y directions
Nx, Ny = 8, 5      # Number of grid points in x and y directions (hx = 1/5 as specified)

# Grid spacing
hx = Lx / (Nx - 1)
hy = Ly / (Ny - 1)

# Create grid points
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)

# Initialize phi with zeros
phi = np.zeros((Ny, Nx))

# Apply boundary conditions
phi[:, 0] = 0                      # Left boundary: phi(0, y) = 0
phi[:, -1] = 0                     # Right boundary: phi(Lx, y) = 0
phi[0, :] = 0                      # Bottom boundary: phi(x, 0) = 0
phi[-1, int(0.2/hx):int(1.2/hx)+1] = 100  # Top boundary: phi(x, Ly) = 100 for 0.2 <= x <= 1.2

# Finite difference method parameters
tolerance = 1e-4  # Convergence criterion
max_iterations = 1000

# Solve the 2D Laplace equation using the finite difference method
for it in range(max_iterations):
    phi_old = phi.copy()
    
    # Update interior points using the average of neighboring points
    for i in range(1, Ny - 1):
        for j in range(1, Nx - 1):
            phi[i, j] = 0.25 * (phi_old[i+1, j] + phi_old[i-1, j] + phi_old[i, j+1] + phi_old[i, j-1])
    
    # Check for convergence
    if np.max(np.abs(phi - phi_old)) < tolerance:
        break

# Create a meshgrid for plotting
X, Y = np.meshgrid(x, y)

# Plot the solution
plt.figure(figsize=(8, 5))
cp = plt.contourf(X, Y, phi, cmap='viridis')
plt.colorbar(cp)
plt.title('Solution to the 2D Laplace Equation')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

P217

import numpy as np
import matplotlib.pyplot as plt

# 離散フーリエ変換 (DFT) を手動で定義
def DFT(x):
    N = len(x)
    X = np.zeros(N, dtype=complex)
    for k in range(N):
        for n in range(N):
            X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return X

# 逆離散フーリエ変換 (IDFT) を手動で定義
def IDFT(X):
    N = len(X)
    x = np.zeros(N, dtype=complex)
    for n in range(N):
        for k in range(N):
            x[n] += X[k] * np.exp(2j * np.pi * k * n / N)
    return x / N

# 入力データ例 (時間領域信号)
x = np.array([0, 1, 2, 3], dtype=float)

# DFT の計算
X = DFT(x)
print("DFT の結果:", X)

# IDFT の計算
x_reconstructed = IDFT(X)
print("再構築された信号:", x_reconstructed)

# 結果のプロット
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.stem(np.abs(X), use_line_collection=True)
plt.title("DFT (Magnitude)")
plt.xlabel("Frequency index")
plt.ylabel("Magnitude")

plt.subplot(1, 2, 2)
plt.stem(x_reconstructed.real, use_line_collection=True)
plt.title("Reconstructed Signal (IDFT)")
plt.xlabel("Time index")
plt.ylabel("Amplitude")

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Sample data (you can modify this based on your specific data)
x = np.array([1, 2, 3, 4, 5, 4, 3, 2])  # Example signal
N = len(x)  # Number of points

# Discrete Fourier Transform (DFT)
def dft(x):
    N = len(x)
    X = np.zeros(N, dtype=complex)
    for k in range(N):
        for n in range(N):
            X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return X

# Inverse Discrete Fourier Transform (IDFT)
def idft(X):
    N = len(X)
    x_reconstructed = np.zeros(N, dtype=complex)
    for n in range(N):
        for k in range(N):
            x_reconstructed[n] += X[k] * np.exp(2j * np.pi * k * n / N)
        x_reconstructed[n] /= N
    return x_reconstructed

# Compute the DFT
X = dft(x)

# Compute the inverse DFT
x_reconstructed = idft(X)

# Plot the original signal and its DFT magnitude
plt.figure(figsize=(10, 6))

# Original signal
plt.subplot(2, 1, 1)
plt.stem(x, use_line_collection=True)
plt.title("Original Signal")
plt.xlabel("Sample")
plt.ylabel("Amplitude")

# DFT Magnitude
plt.subplot(2, 1, 2)
plt.stem(np.abs(X), use_line_collection=True)
plt.title("DFT Magnitude")
plt.xlabel("Frequency (k)")
plt.ylabel("Magnitude")

plt.tight_layout()
plt.show()

# Print the results
print("Original signal:", x)
print("DFT of the signal:", X)
print("Reconstructed signal using IDFT:", np.real(x_reconstructed))

import numpy as np
import matplotlib.pyplot as plt

# Given sequence (approximate values)
x = np.array([1, 5 / (4 * np.sqrt(2)), 3, 3 / (4 * np.sqrt(2)), -7 / (4 * np.sqrt(2)), -3, -5 / (4 * np.sqrt(2)), 0])

# Compute DFT using numpy's fft
X = np.fft.fft(x)

# Plotting the magnitude of the DFT
plt.figure(figsize=(10, 6))

# Plot real and imaginary components
plt.subplot(2, 1, 1)
plt.stem(np.abs(X), use_line_collection=True)
plt.title('DFT Magnitude')
plt.xlabel('Frequency (k)')
plt.ylabel('Magnitude')

# Plot phase angle of the DFT
plt.subplot(2, 1, 2)
plt.stem(np.angle(X), use_line_collection=True)
plt.title('DFT Phase')
plt.xlabel('Frequency (k)')
plt.ylabel('Phase (radians)')

plt.tight_layout()
plt.show()

# Output DFT result
print("DFT of the sequence: ", X)

import numpy as np
import matplotlib.pyplot as plt

# Define the signal (example sine wave + noise)
Fs = 1024  # Sampling frequency
T = 1 / Fs  # Sampling interval
t = np.arange(0, 1, T)  # Time vector (1 second of data)
f1 = 10  # Frequency of the first sine wave
f2 = 60  # Frequency of the second sine wave
signal = np.sin(2 * np.pi * f1 * t) + 0.5 * np.cos(2 * np.pi * f2 * t) + 0.1 * np.random.randn(len(t))

# Compute FFT
N = len(signal)  # Number of samples
fft_vals = np.fft.fft(signal)
fft_vals = fft_vals[:N // 2]  # Take the positive frequency half
frequencies = np.fft.fftfreq(N, T)[:N // 2]  # Frequency values for FFT

# Compute Power Spectral Density (PSD)
psd = (np.abs(fft_vals) ** 2) / N  # Power Spectral Density (normalized)

# Plot the signal
plt.figure(figsize=(10, 8))

plt.subplot(3, 1, 1)
plt.plot(t, signal)
plt.title('Original Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

# Plot the FFT
plt.subplot(3, 1, 2)
plt.stem(frequencies, np.abs(fft_vals), use_line_collection=True)
plt.title('FFT of the Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')

# Plot the Power Spectral Density (PSD)
plt.subplot(3, 1, 3)
plt.plot(frequencies, psd)
plt.title('Power Spectral Density (PSD)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power (dB/Hz)')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the parameters for the function
A = 1  # Amplitude constant
a = 1  # Scale constant
u_vals = np.linspace(-10, 10, 400)  # Range for u
v_vals = np.linspace(-10, 10, 400)  # Range for v

# Create meshgrid for u and v
U, V = np.meshgrid(u_vals, v_vals)

# Define the 2D Fourier Transform function
def F_uv(u, v, A, a):
    # Avoid division by zero by replacing zero values with a small number
    u = np.where(u == 0, 1e-9, u)
    v = np.where(v == 0, 1e-9, v)
    return 4 * (a ** 2) * A * (np.sin(2 * np.pi * a * u) / (2 * np.pi * a * u)) * (np.sin(2 * np.pi * a * v) / (2 * np.pi * a * v))

# Compute the function for all (u, v) values
F = F_uv(U, V, A, a)

# Plot the 2D Fourier Transform
plt.figure(figsize=(8, 6))
plt.contourf(U, V, F, cmap='jet')
plt.colorbar(label='Amplitude')
plt.title('2D Fourier Transform $F(u,v)$')
plt.xlabel('u')
plt.ylabel('v')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the grid
N_values = [20, 100]
u_values = [1]
v_values = [0, 1]
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x, y)

# Plot cos and sin for each combination of u, v, and N
for N in N_values:
    for u in u_values:
        for v in v_values:
            # Compute cos and sin
            cos_vals = np.cos(2 * np.pi * (u * X + v * Y))
            sin_vals = np.sin(2 * np.pi * (u * X + v * Y))

            # Plot cos
            plt.figure(figsize=(6, 6))
            plt.contourf(X, Y, cos_vals, cmap='viridis')
            plt.title(f'cos[2π({u}x + {v}y)] with N={N}')
            plt.colorbar()
            plt.show()

            # Plot sin
            plt.figure(figsize=(6, 6))
            plt.contourf(X, Y, sin_vals, cmap='plasma')
            plt.title(f'sin[2π({u}x + {v}y)] with N={N}')
            plt.colorbar()
            plt.show()

P218

import numpy as np
import matplotlib.pyplot as plt

サンプリング数

N = 8

サンプル点 (0 <= θ < 2π)

theta = np.linspace(0, 2 * np.pi, N, endpoint=False)

入力信号の定義

x = 2 * np.cos(2 * theta) - np.cos(3 * theta) + (1/4) * np.sin(3 * theta)

DFT の計算 (numpy を使用)

X = np.fft.fft(x)

DFT の結果表示

print("入力信号 x(n):", x)
print("DFT の結果 X(k):", X)

実部と虚部を分けて表示

print("実部: ", X.real)
print("虚部: ", X.imag)

DFT の結果をプロット

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

入力信号 x(n) のプロット

plt.subplot(1, 2, 1)
plt.stem(x, use_line_collection=True)
plt.title("Input Signal x(n)")
plt.xlabel("Sample Index")
plt.ylabel("Amplitude")

DFT の実部のプロット

plt.subplot(1, 2, 2)
plt.stem(np.abs(X), use_line_collection=True)
plt.title("DFT Magnitude |X(k)|")
plt.xlabel("Frequency Index")
plt.ylabel("Magnitude")

plt.tight_layout()
plt.show()

P268

import numpy as np
import matplotlib.pyplot as plt

サンプルデータの生成 (sin波 + 雑音)

np.random.seed(0)
t = np.linspace(0, 1, 500)
signal = np.sin(2 * np.pi * 5 * t) # 5Hzのsin波
noise = np.random.normal(0, 0.5, t.shape)
x = signal + noise

移動平均法の実装

def moving_average(x, window_size=5):
""" 単純移動平均法 """
return np.convolve(x, np.ones(window_size)/window_size, mode='valid')

移動平均法の適用

window_size = 15
smoothed_signal = moving_average(x, window_size)

結果のプロット

plt.figure(figsize=(12, 6))
plt.plot(t, x, label='Noisy Signal')
plt.plot(t[(window_size-1)//2: -(window_size//2)], smoothed_signal, label='Smoothed Signal', color='orange')
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.title("Moving Average Filter")
plt.legend()
plt.show()

FFTを用いた周波数フィルタリング

def frequency_filter(x, cutoff_freq=10, sampling_rate=500):
""" 高周波成分をカットする低域フィルタの例 """
# フーリエ変換
X = np.fft.fft(x)
frequencies = np.fft.fftfreq(len(x), 1/sampling_rate)

# フィルタの適用
X[np.abs(frequencies) > cutoff_freq] = 0

# 逆フーリエ変換
return np.fft.ifft(X)

フィルタ適用

filtered_signal = frequency_filter(x, cutoff_freq=10)

結果のプロット

plt.figure(figsize=(12, 6))
plt.plot(t, x, label='Noisy Signal')
plt.plot(t, filtered_signal.real, label='Filtered Signal', color='orange')
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.title("Frequency Domain Filtering")
plt.legend()
plt.show()

同じ信号を複数回観測し、積算平均を行う例

def accumulation_average(signals):
""" 積算平均化処理 """
return np.mean(signals, axis=0)

シミュレーション用データの生成

num_trials = 50
signals = [signal + np.random.normal(0, 0.5, t.shape) for _ in range(num_trials)]

積算平均化の適用

averaged_signal = accumulation_average(signals)

結果のプロット

plt.figure(figsize=(12, 6))
plt.plot(t, x, label='Noisy Signal')
plt.plot(t, averaged_signal, label='Accumulated Average Signal', color='orange')
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.title("Accumulation Average Filtering")
plt.legend()
plt.show()

自己相関関数の計算

def autocorrelation(x):
""" 自己相関関数の計算 """
result = np.correlate(x, x, mode='full')
return result[result.size // 2:]

相互相関関数の計算

def cross_correlation(x, y):
""" 相互相関関数の計算 """
return np.correlate(x, y, mode='full')

自己相関のプロット

auto_corr = autocorrelation(signal)
plt.figure(figsize=(12, 6))
plt.plot(auto_corr)
plt.xlabel("Lag")
plt.ylabel("Correlation")
plt.title("Autocorrelation of Signal")
plt.show()

相互相関の例（同じsin波と位相をずらしたsin波の相互相関）

y = np.sin(2 * np.pi * 5 * (t + 0.1)) # 位相を0.1ずらしたsin波
cross_corr = cross_correlation(signal, y)
plt.figure(figsize=(12, 6))
plt.plot(cross_corr)
plt.xlabel("Lag")
plt.ylabel("Correlation")
plt.title("Cross-correlation of Signal and Shifted Signal")
plt.show()

P274

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftshift
from scipy.signal import get_window

窓関数のリストとパラメータ

window_types = ['boxcar', 'hann', 'hamming', 'blackman', 'kaiser']
beta = 8.6 # Kaiser 窓のパラメータ (β)
N = 512 # 窓関数の長さ
sampling_rate = 1000 # サンプリング周波数 (Hz)

各窓関数を生成して周波数スペクトルをプロット

plt.figure(figsize=(12, 12))

for i, win_type in enumerate(window_types):
# 窓関数の生成
if win_type == 'kaiser':
window = get_window((win_type, beta), N)
else:
window = get_window(win_type, N)

# 窓関数の周波数スペクトルを計算
freq_response = fftshift(fft(window, 2048))  # ゼロパディングを追加してスペクトルを計算
freq_axis = np.linspace(-sampling_rate/2, sampling_rate/2, len(freq_response))

# 時間領域の窓関数をプロット
plt.subplot(len(window_types), 2, 2*i+1)
plt.plot(window)
plt.title(f"{win_type.capitalize()} Window (Time Domain)")
plt.xlabel("Sample Index")
plt.ylabel("Amplitude")
plt.grid()

# 周波数領域の窓関数をプロット (20 * log10(|X|) を表示)
plt.subplot(len(window_types), 2, 2*i+2)
plt.plot(freq_axis, 20 * np.log10(np.abs(freq_response) / np.max(np.abs(freq_response))))
plt.title(f"{win_type.capitalize()} Window (Frequency Domain)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.ylim(-100, 0)
plt.grid()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import pywt

サンプル信号の生成

def create_signal():
t = np.linspace(0, 1, 400, endpoint=False)
signal = np.cos(2 * np.pi * 7 * t) + np.sin(2 * np.pi * 3 * t) + np.random.normal(0, 0.5, t.shape)
return t, signal

1. 離散ウェーブレット変換（DWT）

def discrete_wavelet_transform(signal, wavelet='db4', level=3):
coeffs = pywt.wavedec(signal, wavelet, level=level)
return coeffs

2. 連続ウェーブレット変換（CWT）

def continuous_wavelet_transform(signal, wavelet='cmor'):
scales = np.arange(1, 128) # スケールの定義
coefficients, frequencies = pywt.cwt(signal, scales, wavelet)
return coefficients, frequencies

3. マザーウェーブレットの生成とプロット

def plot_mother_wavelet(wavelet='cmor'):
w = pywt.ContinuousWavelet(wavelet)
psi, x = w.wavefun(level=8)
plt.figure(figsize=(6, 3))
plt.plot(x, psi.real, label='Real part')
plt.plot(x, psi.imag, label='Imaginary part')
plt.title(f'Mother Wavelet: {wavelet}')
plt.legend()
plt.show()

4. 多重解像度解析（MRA）

def multiresolution_analysis(signal, wavelet='db4', level=3):
coeffs = discrete_wavelet_transform(signal, wavelet, level)
approx = coeffs[0] # 近似成分
details = coeffs[1:] # 詳細成分

plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.plot(approx)
plt.title('Approximation Coefficients')

for i, detail in enumerate(details):
    plt.subplot(level, 1, i + 2)
    plt.plot(detail)
    plt.title(f'Detail Coefficients Level {i+1}')
plt.tight_layout()
plt.show()
return approx, details

メイン関数

def main():
# 信号生成
t, signal = create_signal()

# 1. 離散ウェーブレット変換の実行とプロット
coeffs = discrete_wavelet_transform(signal)
plt.figure(figsize=(10, 6))
plt.subplot(3, 1, 1)
plt.plot(signal)
plt.title("Original Signal")

for i, coeff in enumerate(coeffs):
    plt.subplot(3, 1, i + 2)
    plt.plot(coeff)
    plt.title(f"Level {i} Coefficients")
plt.tight_layout()
plt.show()

# 2. 連続ウェーブレット変換の実行とプロット
cwt_coeffs, freqs = continuous_wavelet_transform(signal)
plt.figure(figsize=(10, 6))
plt.imshow(np.abs(cwt_coeffs), extent=[0, 1, 1, 128], cmap='viridis', aspect='auto')
plt.title("Continuous Wavelet Transform (CWT)")
plt.ylabel("Scales")
plt.xlabel("Time")
plt.colorbar(label='Magnitude')
plt.show()

# 3. マザーウェーブレットのプロット
plot_mother_wavelet()

# 4. 多重解像度解析の実行とプロット
multiresolution_analysis(signal)

スクリプトの実行

if name == "main":
main()

P282

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

窓関数の種類と長さを設定

window_length = 51
windows = {
"Rectangular (Boxcar)": signal.boxcar(window_length),
"Hann": signal.hann(window_length),
"Hamming": signal.hamming(window_length),
"Blackman": signal.blackman(window_length),
"Kaiser (beta=14)": signal.kaiser(window_length, beta=14),
"Gaussian (std=7)": signal.gaussian(window_length, std=7),
"Bartlett": signal.bartlett(window_length),
"Parzen": signal.parzen(window_length),
"Sine": np.sin(np.linspace(0, np.pi, window_length)),
"Welch": signal.welch(np.linspace(0, 1, window_length), nperseg=window_length)[1]
}

各窓関数のプロット

plt.figure(figsize=(12, 8))
for idx, (name, window) in enumerate(windows.items()):
plt.subplot(5, 2, idx+1)
plt.plot(window)
plt.title(name)
plt.grid()

plt.tight_layout()
plt.show()

import numpy as np

1. 拡大・縮小行列

def scaling_matrix(sx, sy, sz=1):
return np.array([[sx, 0, 0, 0],
[0, sy, 0, 0],
[0, 0, sz, 0],
[0, 0, 0, 1]])

2. 平行移動行列

def translation_matrix(tx, ty, tz=0):
return np.array([[1, 0, 0, tx],
[0, 1, 0, ty],
[0, 0, 1, tz],
[0, 0, 0, 1]])

3. 回転行列（方向余弦行列）

def rotation_matrix_z(theta):
cos_theta = np.cos(theta)
sin_theta = np.sin(theta)
return np.array([[cos_theta, -sin_theta, 0, 0],
[sin_theta, cos_theta, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])

4. 鏡映行列

def reflection_matrix(axis='x'):
if axis == 'x':
return np.array([[-1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])
elif axis == 'y':
return np.array([[1, 0, 0, 0],
[0, -1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])
elif axis == 'z':
return np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, -1, 0],
[0, 0, 0, 1]])

5. せん断（スキュー変換）行列

def shear_matrix(shx=0, shy=0, shz=0):
return np.array([[1, shx, shz, 0],
[shy, 1, shz, 0],
[shx, shy, 1, 0],
[0, 0, 0, 1]])

6. 剛体変換とアフィン変換

def rigid_body_transform(translation, rotation_angle):
translation_mat = translation_matrix(*translation)
rotation_mat = rotation_matrix_z(rotation_angle)
return np.dot(translation_mat, rotation_mat)

7. アフィン変換の合成

def affine_transform(translation, scaling, rotation_angle, shear):
translation_mat = translation_matrix(*translation)
scaling_mat = scaling_matrix(*scaling)
rotation_mat = rotation_matrix_z(rotation_angle)
shear_mat = shear_matrix(*shear)
return np.dot(translation_mat, np.dot(rotation_mat, np.dot(scaling_mat, shear_mat)))

8. ラドン変換 (ここでは簡略化のために2Dラドン変換の一例を使用)

def radon_transform(image, theta):
from skimage.transform import radon
return radon(image, theta=theta)

9. 投影定理（2Dの場合の実装）

def projection_2d(points, axis='x'):
if axis == 'x':
return points[:, 1:] # y, z 方向への投影
elif axis == 'y':
return points[:, [0, 2]] # x, z 方向への投影
elif axis == 'z':
return points[:, :2] # x, y 方向への投影

10. 円を楕円に変換する行列

def circle_to_ellipse(a, b):
return scaling_matrix(a, b)

11. 球を楕円体に変換する行列

def sphere_to_ellipsoid(a, b, c):
return scaling_matrix(a, b, c)

サンプルの実行例

if name == "main":
# 各種行列の確認例
print("Scaling Matrix (2x, 3x):")
print(scaling_matrix(2, 3))

print("\nTranslation Matrix (2, 3):")
print(translation_matrix(2, 3))

print("\nRotation Matrix (45 degrees):")
print(rotation_matrix_z(np.pi / 4))

print("\nReflection Matrix (x-axis):")
print(reflection_matrix('x'))

print("\nShear Matrix (shx=0.5, shy=0.2):")
print(shear_matrix(0.5, 0.2))

print("\nRigid Body Transform (translation=(2, 3), rotation=45 degrees):")
print(rigid_body_transform([2, 3], np.pi / 4))

print("\nAffine Transform (translation=(2, 3), scaling=(2, 3), rotation=45 degrees, shear=(0.5, 0.2)):")
print(affine_transform([2, 3], [2, 3], np.pi / 4, [0.5, 0.2]))

print("\nMatrix to transform circle to ellipse (a=2, b=3):")
print(circle_to_ellipse(2, 3))

print("\nMatrix to transform sphere to ellipsoid (a=2, b=3, c=4):")
print(sphere_to_ellipsoid(2, 3, 4))

import sympy as sp
import numpy as np
from scipy.integrate import quad

勾配の計算

スカラ関数 f(x, y, z) = x2 + y2 + z**2

x, y, z = sp.symbols('x y z')
f = x2 + y2 + z**2

スカラ関数の勾配を計算

grad_f = sp.Matrix([sp.diff(f, var) for var in (x, y, z)])
print("スカラ関数の勾配 (grad f):")
sp.pprint(grad_f)

線積分の計算例

ベクトル場 A = [y, -x, 0]

A = sp.Matrix([y, -x, 0])

線積分のルートを定義

例として、円 (x, y, z) = (cos(t), sin(t), 0) をパラメータ t で表現する

t = sp.symbols('t')
r = sp.Matrix([sp.cos(t), sp.sin(t), 0])

線積分 ∫_C A ・ dr の計算

dr = r.diff(t)
line_integral = A.dot(dr)

0 から 2π までの線積分を計算

result = sp.integrate(line_integral, (t, 0, 2 * sp.pi))
print("\n線積分の結果 ∫_C A ・ dr:")
sp.pprint(result)

流線方向微分の例

流線方向の単位ベクトルを定義 (例として [1, 1, 1] を単位ベクトル化)

e = sp.Matrix([1, 1, 1])
e_unit = e / e.norm()

勾配と単位ベクトルとの内積を計算

flow_direction_derivative = grad_f.dot(e_unit)
print("\n流線方向微分:")
sp.pprint(flow_direction_derivative)