Pythonプログラミング演習で学ぶ高専応用数学

Last updated at 2025-09-23Posted at 2025-09-23

ベクトル
- スカラーとベクトルの違い
- 内積と外積の定義と幾何学的意味
- ベクトル方程式（直線・平面・球面）
勾配，発散，回転
- 勾配（Gradient）の意味
- 発散（Divergence）
- 回転（Curl）
線積分と面積分
- 線積分（仕事・電位差の計算）
- 面積分（流量・磁束の計算）
- パラメトリック曲線と曲面積分
ガウスの発散定理とストークスの定理
- 発散定理（体積積分と面積分の関係）
- ストークスの定理（線積分と面積分の関係）

第2章複素関数論

複素数
- 複素数の表現（直交形式・極形式）
- 複素平面とオイラーの公式
- 複素数の演算（加減乗除・共役・絶対値）
複素関数
- 複素関数の定義域と値域
- 正則関数とコーシー・リーマンの関係式
- 基本的な複素関数（exp, sin, cos, log, 1/z）の性質
複素関数の積分
- コーシーの積分定理
- コーシーの積分公式
- 複素積分と実積分の関係
ローラン展開と留数定理
- ローラン級数展開
- 留数の計算と留数定理

第3章ラプラス変換

ラプラス変換
- 定義と基本性質（線形性・微分積分の変換）
- 基本関数のラプラス変換（指数・正弦波・ステップ関数）
- 部分分数分解と逆ラプラス変換
デルタ関数と線形システム
- デルタ関数の定義と直観的意味
- 線形システムの畳み込み積分

第4章フーリエ級数とフーリエ変換

フーリエ級数
- 周期関数のフーリエ展開
- 偶関数・奇関数の展開公式
- ギブス現象と収束性
フーリエ変換
- 定義と直観的意味（時間領域と周波数領域）
- 基本関数のフーリエ変換（ガウス関数・デルタ関数）
- 周波数特性と畳み込み定理
FFT（高速フーリエ変換）
- 離散フーリエ変換（DFT）の定義
- FFTアルゴリズムの基本原理
- FFTの応用例（音声解析・画像処理・通信工学）
- サンプリングとナイキスト周波数の注意点

第5章微分方程式

微分方程式の基礎
- 微分方程式の定義
- 常微分方程式と偏微分方程式の違い
- 階数と線形性の分類
1階常微分方程式
- 変数分離形の解法
- 同次形の解法
- 1階線形微分方程式の解法
2階常微分方程式
- 2階同次微分方程式（特性方程式による解法）
- 2階非同次微分方程式（特解の求め方）
高階常微分方程式
- n階線形微分方程式
- 特性方程式と一般解
連立微分方程式
- 1階連立微分方程式
- 行列を用いた表現と解法
変数係数・非線形微分方程式
- 変数係数方程式の例（オイラー方程式）
- 非線形微分方程式の特徴
偏微分方程式（入門）
- 偏微分方程式の基礎
- 代表的な偏微分方程式（波動・熱伝導・ラプラス方程式）

第1章ベクトル解析 — Python課題

1. ベクトル

【問題1】
ベクトル a=(2,3)、b=(1,-1) の和 a+b を Python で計算せよ。

【問題2】
ベクトル a=(3,4) の大きさ（ノルム）を計算せよ。

【問題3】
ベクトル a=(1,2,3)、b=(4,5,6) の内積を計算せよ。

【問題4】
ベクトル a=(1,0,0)、b=(0,1,0) の外積を Python で計算せよ。

【問題5】
ベクトル a=(1,2,3)、b=(2,1,0) の成す角を求めよ。

【問題6】
直線 x=1+t, y=2t, z=3t をベクトル方程式で表し、Pythonで t=0〜5 の点を計算せよ。

【問題7】
平面 2x+3y+z=6 の法線ベクトルを Python で表現せよ。

【問題8】
球面 x²+y²+z²=9 を半径 r=3 としてPythonで可視化せよ（matplotlib の 3Dプロット）。

2. 勾配，発散，回転

【問題9】
f(x,y)=x²+y² の勾配を Python でシンボリック計算せよ（SymPy 使用）。

【問題10】
f(x,y,z)=xyz の勾配を Python で求めよ。

【問題11】
ベクトル場 F=(x,y,z) の発散を Python で計算せよ。

【問題12】
ベクトル場 F=(x², y², z²) の発散を求めよ。

【問題13】
ベクトル場 F=(y, -x, 0) の回転（Curl）を求めよ。

【問題14】
ベクトル場 F=(yz, xz, xy) の回転を Python で計算せよ。

【問題15】
ベクトル場 F=(x,y,z) の発散と回転を同時に求めよ。

【問題16】
f(x,y)=sin(x)cos(y) の勾配ベクトルを Python で描画せよ（quiver プロット）。

3. 線積分と面積分

【問題17】
ベクトル場 F=(y,0) の線積分 ∫C F・dr を C: (0,0)→(1,1) で計算せよ。

【問題18】
F=(x,y) を C: x=t, y=t², t∈[0,1] に沿って線積分せよ。

【問題19】
F=(y,-x) の線積分を C: 単位円上で計算せよ。

【問題20】
F=(x,0,0) の面積分を平面 S: y=0, z=0, 0≤x≤1 で計算せよ。

【問題21】
F=(x,y,z) の面積分を単位立方体の表面で計算せよ。

【問題22】
パラメトリック曲線 C: x=cos t, y=sin t, 0≤t≤2π に沿って F=(x,y) を線積分せよ。

【問題23】
パラメトリック曲面 S: x=u cos v, y=u sin v, z=u², (0≤u≤1, 0≤v≤2π) の面積を Python で計算せよ。

【問題24】
流体の速度場 F=(x,y) において、正方形領域 [0,1]×[0,1] の流量を Python で数値積分せよ。

4. ガウスの発散定理とストークスの定理

【問題25】
F=(x,y,z) の発散を計算し、単位立方体での体積積分を Python で求めよ。

【問題26】
同じF=(x,y,z) の面積分を立方体表面で計算し、発散定理を確認せよ。

【問題27】
F=(y,-x,0) の回転を計算し、半径1の円盤上の面積分を求めよ。

【問題28】
F=(y,-x,0) の線積分を単位円上で求め、ストークスの定理を確認せよ。

【問題29】
F=(z,0,0) の発散を求め、球面 r=1 上で発散定理を確認せよ。

【問題30】
F=(0,z,-y) の回転を求め、半球面でストークスの定理をPythonで確認せよ。

答え


import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# Problem 1: Vector addition
# --------------------------
a1 = np.array([2, 3])
b1 = np.array([1, -1])
print("Problem 1:", a1 + b1)

# --------------------------
# Problem 2: Vector norm
# --------------------------
a2 = np.array([3, 4])
print("Problem 2:", np.linalg.norm(a2))

# --------------------------
# Problem 3: Dot product
# --------------------------
a3 = np.array([1, 2, 3])
b3 = np.array([4, 5, 6])
print("Problem 3:", np.dot(a3, b3))

# --------------------------
# Problem 4: Cross product
# --------------------------
a4 = np.array([1, 0, 0])
b4 = np.array([0, 1, 0])
print("Problem 4:", np.cross(a4, b4))

# --------------------------
# Problem 5: Angle between vectors
# --------------------------
a5 = np.array([1, 2, 3])
b5 = np.array([2, 1, 0])
cos_theta = np.dot(a5, b5) / (np.linalg.norm(a5) * np.linalg.norm(b5))
theta = np.arccos(cos_theta)
print("Problem 5 (radians):", theta)
print("Problem 5 (degrees):", np.degrees(theta))

# --------------------------
# Problem 6: Parametric line
# --------------------------
t = np.linspace(0, 5, 20)
x = 1 + t
y = 2 * t
z = 3 * t
points_line = np.vstack((x, y, z)).T
print("Problem 6 (first 5 points):\n", points_line[:5])

# --------------------------
# Problem 7: Plane normal vector
# --------------------------
normal = np.array([2, 3, 1])
print("Problem 7:", normal)

# --------------------------
# Problem 8: Sphere visualization
# --------------------------
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection='3d')
u = np.linspace(0, 2*np.pi, 50)
v = np.linspace(0, np.pi, 50)
r = 3
X = r * np.outer(np.cos(u), np.sin(v))
Y = r * np.outer(np.sin(u), np.sin(v))
Z = r * np.outer(np.ones_like(u), np.cos(v))
ax.plot_surface(X, Y, Z, color='lightblue', edgecolor='k', alpha=0.5)
ax.set_title("Problem 8: Sphere r=3")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
plt.show()


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

# --------------------------
# Problem 9: Gradient of f(x,y)=x^2+y^2
# --------------------------
x, y, z = sp.symbols('x y z')
f1 = x**2 + y**2
grad_f1 = [sp.diff(f1, x), sp.diff(f1, y)]
print("Problem 9:", grad_f1)

# --------------------------
# Problem 10: Gradient of f(x,y,z)=xyz
# --------------------------
f2 = x*y*z
grad_f2 = [sp.diff(f2, var) for var in (x,y,z)]
print("Problem 10:", grad_f2)

# --------------------------
# Problem 11: Divergence of F=(x,y,z)
# --------------------------
F1 = sp.Matrix([x, y, z])
div_F1 = sp.diff(F1[0], x) + sp.diff(F1[1], y) + sp.diff(F1[2], z)
print("Problem 11:", div_F1)

# --------------------------
# Problem 12: Divergence of F=(x^2,y^2,z^2)
# --------------------------
F2 = sp.Matrix([x**2, y**2, z**2])
div_F2 = sp.diff(F2[0], x) + sp.diff(F2[1], y) + sp.diff(F2[2], z)
print("Problem 12:", div_F2)

# --------------------------
# Problem 13: Curl of F=(y,-x,0)
# --------------------------
F3 = sp.Matrix([y, -x, 0])
curl_F3 = sp.Matrix([
    sp.diff(F3[2], y) - sp.diff(F3[1], z),
    sp.diff(F3[0], z) - sp.diff(F3[2], x),
    sp.diff(F3[1], x) - sp.diff(F3[0], y)
])
print("Problem 13:", curl_F3)

# --------------------------
# Problem 14: Curl of F=(yz,xz,xy)
# --------------------------
F4 = sp.Matrix([y*z, x*z, x*y])
curl_F4 = sp.Matrix([
    sp.diff(F4[2], y) - sp.diff(F4[1], z),
    sp.diff(F4[0], z) - sp.diff(F4[2], x),
    sp.diff(F4[1], x) - sp.diff(F4[0], y)
])
print("Problem 14:", curl_F4)

# --------------------------
# Problem 15: Divergence and Curl of F=(x,y,z)
# --------------------------
div_F5 = sp.diff(x, x) + sp.diff(y, y) + sp.diff(z, z)
curl_F5 = sp.Matrix([
    sp.diff(z, y) - sp.diff(y, z),
    sp.diff(x, z) - sp.diff(z, x),
    sp.diff(y, x) - sp.diff(x, y)
])
print("Problem 15 (div):", div_F5)
print("Problem 15 (curl):", curl_F5)

# --------------------------
# Problem 16: Gradient vector field visualization f(x,y)=sin(x)cos(y)
# --------------------------
f3 = sp.sin(x)*sp.cos(y)
grad_f3 = [sp.diff(f3, x), sp.diff(f3, y)]

# Numerical evaluation on grid
X, Y = np.meshgrid(np.linspace(-2*np.pi, 2*np.pi, 20),
                   np.linspace(-2*np.pi, 2*np.pi, 20))
grad_fx = sp.lambdify((x,y), grad_f3[0], "numpy")
grad_fy = sp.lambdify((x,y), grad_f3[1], "numpy")
U = grad_fx(X,Y)
V = grad_fy(X,Y)

plt.figure(figsize=(6,6))
plt.quiver(X, Y, U, V, color="blue")
plt.title("Problem 16: Gradient field of f(x,y)=sin(x)cos(y)")
plt.xlabel("x")
plt.ylabel("y")
plt.axis("equal")
plt.show()


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

x, y, z, t, u, v, r, theta, phi = sp.symbols('x y z t u v r theta phi')

# ======================================================
# Problem 17
F = sp.Matrix([y, 0])
dx, dy = 1, 1
dr = sp.Matrix([dx, dy])
res17 = F.subs({x:1,y:1}).dot(dr)
print("Problem 17:", res17)

# ======================================================
# Problem 18
x_t, y_t = t, t**2
dxdt, dydt = sp.diff(x_t,t), sp.diff(y_t,t)
F = sp.Matrix([x, y])
dr = sp.Matrix([dxdt,dydt])
res18 = sp.integrate(F.subs({x:x_t,y:y_t}).dot(dr),(t,0,1))
print("Problem 18:", res18)

# ======================================================
# Problem 19
x_t, y_t = sp.cos(t), sp.sin(t)
dxdt, dydt = sp.diff(x_t,t), sp.diff(y_t,t)
F = sp.Matrix([y,-x])
dr = sp.Matrix([dxdt,dydt])
res19 = sp.integrate(F.subs({x:x_t,y:y_t}).dot(dr),(t,0,2*sp.pi))
print("Problem 19:", res19)

# ======================================================
# Problem 20
F = sp.Matrix([x,0,0])
res20 = sp.integrate(F[0].subs({y:0,z:0}), (x,0,1))
print("Problem 20:", res20)

# ======================================================
# Problem 21
F = sp.Matrix([x,y,z])
flux = 0
flux += sp.integrate(F[0].subs({x:1}),(y,0,1),(z,0,1)) - sp.integrate(F[0].subs({x:0}),(y,0,1),(z,0,1))
flux += sp.integrate(F[1].subs({y:1}),(x,0,1),(z,0,1)) - sp.integrate(F[1].subs({y:0}),(x,0,1),(z,0,1))
flux += sp.integrate(F[2].subs({z:1}),(x,0,1),(y,0,1)) - sp.integrate(F[2].subs({z:0}),(x,0,1),(y,0,1))
print("Problem 21:", flux)

# ======================================================
# Problem 22
x_t, y_t = sp.cos(t), sp.sin(t)
dxdt, dydt = sp.diff(x_t,t), sp.diff(y_t,t)
F = sp.Matrix([x,y])
dr = sp.Matrix([dxdt,dydt])
res22 = sp.integrate(F.subs({x:x_t,y:y_t}).dot(dr),(t,0,2*sp.pi))
print("Problem 22:", res22)

# ======================================================
# Problem 23
X = sp.Matrix([u*sp.cos(v), u*sp.sin(v), u**2])
Xu = X.diff(u)
Xv = X.diff(v)
dS = Xu.cross(Xv).norm()
res23 = sp.integrate(dS,(u,0,1),(v,0,2*sp.pi))
print("Problem 23:", res23)

# ======================================================
# Problem 24
F = sp.Matrix([x,y])
res24 = sp.integrate(F[0],(x,0,1),(y,0,1)) + sp.integrate(F[1],(x,0,1),(y,0,1))
print("Problem 24:", res24)

# ======================================================
# Problem 25
F = sp.Matrix([x,y,z])
divF = sp.diff(F[0],x)+sp.diff(F[1],y)+sp.diff(F[2],z)
res25 = sp.integrate(divF,(x,0,1),(y,0,1),(z,0,1))
print("Problem 25:", res25)

# ======================================================
# Problem 26
flux = 0
flux += sp.integrate(F[0].subs({x:1}),(y,0,1),(z,0,1)) - sp.integrate(F[0].subs({x:0}),(y,0,1),(z,0,1))
flux += sp.integrate(F[1].subs({y:1}),(x,0,1),(z,0,1)) - sp.integrate(F[1].subs({y:0}),(x,0,1),(z,0,1))
flux += sp.integrate(F[2].subs({z:1}),(x,0,1),(y,0,1)) - sp.integrate(F[2].subs({z:0}),(x,0,1),(y,0,1))
print("Problem 26:", flux)

# ======================================================
# Problem 27
F = sp.Matrix([y,-x,0])
curlF = sp.Matrix([
    sp.diff(F[2],y)-sp.diff(F[1],z),
    sp.diff(F[0],z)-sp.diff(F[2],x),
    sp.diff(F[1],x)-sp.diff(F[0],y)
])
expr = curlF[2].subs({x:r*sp.cos(theta),y:r*sp.sin(theta)})
res27 = sp.integrate(expr*r,(r,0,1),(theta,0,2*sp.pi))
print("Problem 27:", res27)

# ======================================================
# Problem 28
x_t, y_t = sp.cos(t), sp.sin(t)
dxdt, dydt = sp.diff(x_t,t), sp.diff(y_t,t)
F = sp.Matrix([y,-x,0])
dr = sp.Matrix([dxdt,dydt,0])
res28 = sp.integrate(F.subs({x:x_t,y:y_t}).dot(dr),(t,0,2*sp.pi))
print("Problem 28:", res28)

# ======================================================
# Problem 29
F = sp.Matrix([z,0,0])
divF = sp.diff(F[0],x)+sp.diff(F[1],y)+sp.diff(F[2],z)
expr = divF.subs({x:r*sp.sin(theta)*sp.cos(phi),y:r*sp.sin(theta)*sp.sin(phi),z:r*sp.cos(theta)})
jacobian = r**2*sp.sin(theta)
res29 = sp.integrate(expr*jacobian,(r,0,1),(theta,0,sp.pi),(phi,0,2*sp.pi))
print("Problem 29:", res29)

# ======================================================
# Problem 30
F = sp.Matrix([0,z,-y])
curlF = sp.Matrix([
    sp.diff(F[2],y)-sp.diff(F[1],z),
    sp.diff(F[0],z)-sp.diff(F[2],x),
    sp.diff(F[1],x)-sp.diff(F[0],y)
])
print("Problem 30 curl:", curlF)

# ======================================================
# Visualization: Vector field example for Problem 19
X,Y = np.meshgrid(np.linspace(-1,1,10),np.linspace(-1,1,10))
U,V = Y,-X
plt.figure(figsize=(5,5))
plt.quiver(X,Y,U,V)
plt.gca().set_aspect('equal')
plt.title("Vector field F=(y,-x) (Problem 19)")
plt.xlabel("x"); plt.ylabel("y")
plt.show()

第2章複素関数論 — Python課題

1. 複素数

【問題1】
複素数 z=3+4j の絶対値を Python で計算せよ。

【問題2】
複素数 z=1-√3 j の偏角（角度）を Python で計算せよ。

【問題3】
z=2+3j と w=1-4j の和と積を Python で求めよ。

【問題4】
複素数 z=exp(1j * π/3) を Python で計算し、直交形式で表示せよ。

【問題5】
オイラーの公式 exp(iθ)=cosθ+i sinθ を Python で確認せよ（θ=π/4 の場合）。

【問題6】
複素数 z=5∠60° を Python で直交形式に変換せよ。

【問題7】
z=3+4j の共役複素数を Python で求めよ。

【問題8】
交流回路のインピーダンス Z=R+jωL を R=10Ω, L=0.1H, ω=100 rad/s で Python で計算せよ。

2. 複素関数

【問題9】
f(z)=z² の値を z=1+2j で計算せよ。

【問題10】
f(z)=1/z の値を z=2+3j で Python で計算せよ。

【問題11】
f(z)=exp(z) を z=1+πj で Python で計算せよ。

【問題12】
f(z)=log(z) を z=1+1j で Python で計算せよ。

【問題13】
sin(z) を z=2+3j で Python で計算せよ。

【問題14】
cos(z) を z=2+3j で Python で計算せよ。

【問題15】
コーシー・リーマン方程式を f(z)=z² の場合に Python (SymPy) で確認せよ。

【問題16】
f(z)=u(x,y)+i v(x,y), u=x²−y², v=2xy としてコーシー・リーマンの条件を Python で確認せよ。

3. 複素関数の積分

【問題17】
f(z)=1/z の積分を単位円上（C:|z|=1）で計算し、2πi になることをPythonで確認せよ。

【問題18】
f(z)=z を単位円上で積分せよ（結果は0）。

【問題19】
f(z)=1/(z−1) の積分を半径2の円上で計算せよ。

【問題20】
コーシーの積分公式 f(a)=1/(2πi)∮ f(z)/(z−a)dz を f(z)=z², a=0 の場合で Python により確認せよ。

【問題21】
f(z)=exp(z)/(z−1) を単位円（中心0, 半径2）で積分し、結果を Python で計算せよ。

【問題22】
実積分 ∫₀^∞ cos(x)/(x²+1) dx を Python で数値積分せよ（複素積分により π/2·e^(−1) となることを確認）。

【問題23】
実積分 ∫₀^∞ sin(x)/x dx を Python で近似計算せよ（ディリクレ積分）。

【問題24】
矩形経路上で f(z)=z² を積分し、結果が経路に依存しないことを Python で確認せよ。

4. ローラン展開と留数定理

【問題25】
f(z)=1/(z−1) を z=0 のまわりでローラン展開せよ。

【問題26】
f(z)=1/(z(1−z)) を |z|<1 でローラン展開せよ。

【問題27】
f(z)=1/z² の z=0 における留数を求めよ。

【問題28】
f(z)=exp(z)/z の z=0 における留数を求めよ。

【問題29】
f(z)=1/(z²+1) の z=i における留数を Python で求めよ。

【問題30】
実積分 ∫₋∞^∞ 1/(x²+1) dx を留数定理を用いて求めよ（結果 π）。

答え


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

j = sp.I   # 複素単位
pi = sp.pi

# ======================================================
# 1. 複素数
# ======================================================

# Problem 1
z = 3 + 4*j
print("Problem 1:", abs(z))

# Problem 2
z = 1 - sp.sqrt(3)*j
print("Problem 2 (arg):", sp.arg(z))

# Problem 3
z, w = 2+3*j, 1-4*j
print("Problem 3 sum:", z+w, " product:", z*w)

# Problem 4
z = sp.exp(j*pi/3)
print("Problem 4:", sp.re(z), "+", sp.im(z), "j")

# Problem 5
theta = pi/4
lhs = sp.exp(j*theta)
rhs = sp.cos(theta) + j*sp.sin(theta)
print("Problem 5 Euler check:", lhs, " vs ", rhs)

# Problem 6
r, deg = 5, 60
z = r*(sp.cos(sp.rad(deg)) + j*sp.sin(sp.rad(deg)))
print("Problem 6:", z.simplify())

# Problem 7
z = 3+4*j
print("Problem 7:", sp.conjugate(z))

# Problem 8
R, L, w = 10, 0.1, 100
Z = R + j*w*L
print("Problem 8:", Z)

# ======================================================
# 2. 複素関数
# ======================================================

# Problem 9
z = 1+2*j
print("Problem 9:", z**2)

# Problem 10
z = 2+3*j
print("Problem 10:", 1/z)

# Problem 11
z = 1+pi*j
print("Problem 11:", sp.exp(z))

# Problem 12
z = 1+1*j
print("Problem 12:", sp.log(z))

# Problem 13
z = 2+3*j
print("Problem 13 sin(z):", sp.sin(z))

# Problem 14
print("Problem 14 cos(z):", sp.cos(z))

# Problem 15 (Cauchy-Riemann check for f=z^2)
x, y = sp.symbols('x y', real=True)
f = (x + j*y)**2
u, v = sp.re(f), sp.im(f)
print("Problem 15 CR eqns:",
      sp.diff(u,x), "=?",
      sp.diff(v,y),
      ", ",
      sp.diff(u,y), "=? -", sp.diff(v,x))

# Problem 16 (u=x^2-y^2, v=2xy)
u = x**2 - y**2
v = 2*x*y
print("Problem 16 CR eqns:",
      sp.diff(u,x), "=?",
      sp.diff(v,y),
      ", ",
      sp.diff(u,y), "=? -", sp.diff(v,x))

# ======================================================
# Visualization example (Argand diagram for Problem 6)
# ======================================================
z_val = complex(sp.N(z))
plt.figure(figsize=(5,5))
plt.quiver(0,0, z_val.real, z_val.imag, angles='xy', scale_units='xy', scale=1)
plt.xlim(-6,6); plt.ylim(-6,6)
plt.axhline(0,color='k'); plt.axvline(0,color='k')
plt.title("Problem 6: Complex number in Argand diagram")
plt.xlabel("Real"); plt.ylabel("Imag")
plt.grid(True); plt.show()


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

j = sp.I
pi = sp.pi
z, t = sp.symbols('z t', real=False)

# ======================================================
# 問題17: f(z)=1/z の積分 (単位円)
z_t = sp.exp(j*t)
dzdt = sp.diff(z_t, t)
res17 = sp.integrate((1/z_t)*dzdt, (t,0,2*pi))
print("P17:", res17)   # → 2πi

# ======================================================
# 問題18: f(z)=z の積分 (単位円)
res18 = sp.integrate(z_t*dzdt, (t,0,2*pi))
print("P18:", res18)   # → 0

# ======================================================
# 問題19: f(z)=1/(z-1), 半径2の円
z_t = 2*sp.exp(j*t)
dzdt = sp.diff(z_t, t)
res19 = sp.integrate((1/(z_t-1))*dzdt, (t,0,2*pi))
print("P19:", sp.simplify(res19))   # → 2πi

# ======================================================
# 問題20: コーシーの積分公式 f(z)=z^2, a=0
f = z**2
z_t = sp.exp(j*t)
dzdt = sp.diff(z_t, t)
res20 = sp.integrate((f/z_t)*dzdt, (t,0,2*pi))
lhs = 1/(2*pi*j)*res20
print("P20:", lhs, " vs f(0)=", f.subs(z,0))  # → 0

# ======================================================
# 問題21: f(z)=exp(z)/(z-1), 半径2の円
z_t = 2*sp.exp(j*t)
dzdt = sp.diff(z_t, t)
f = sp.exp(z)/(z-1)
res21 = sp.integrate(f.subs(z,z_t)*dzdt, (t,0,2*pi))
print("P21:", sp.simplify(res21))   # → 2πi*e

# ======================================================
# 問題22: ∫0∞ cos(x)/(x²+1) dx
res22 = mp.quad(lambda xx: mp.cos(xx)/(xx**2+1), [0, mp.inf])
print("P22:", res22, " exact:", (pi/2)*mp.e**(-1))

# ======================================================
# 問題23: ∫0∞ sin(x)/x dx (Dirichlet積分)
res23 = mp.quad(lambda xx: mp.sin(xx)/xx, [0, mp.inf])
print("P23:", res23, " exact:", pi/2)

# ======================================================
# 問題24: f(z)=z² の矩形経路 (0→1→1+i→i→0)
# 経路分割して積分
z1 = sp.integrate((z**2).subs(z, t+0*j)*1, (t,0,1))      # 0→1
z2 = sp.integrate((z**2).subs(z, 1+t*j)*j, (t,0,1))      # 1→1+i
z3 = sp.integrate((z**2).subs(z, 1-t+1*j)*(-1), (t,0,1)) # 1+i→i
z4 = sp.integrate((z**2).subs(z, 0+(1-t)*j)*(-j), (t,0,1)) # i→0
res24 = z1+z2+z3+z4
print("P24:", res24)   # → 0

# ======================================================
# 可視化: 経路の描画 (問題17, 19, 24)
theta = np.linspace(0, 2*np.pi, 400)
circle1_x, circle1_y = np.cos(theta), np.sin(theta)         # 単位円
circle2_x, circle2_y = 2*np.cos(theta), 2*np.sin(theta)     # 半径2円
rect_x = [0,1,1,0,0]
rect_y = [0,0,1,1,0]

plt.figure(figsize=(6,6))
plt.plot(circle1_x, circle1_y, 'b', label="Unit circle (P17,P18)")
plt.plot(circle2_x, circle2_y, 'g', label="Radius=2 circle (P19,P21)")
plt.plot(rect_x, rect_y, 'r--', label="Rectangle path (P24)")
plt.gca().set_aspect('equal')
plt.title("Integration paths for Section 3 (P17–P24)")
plt.xlabel("Re(z)")
plt.ylabel("Im(z)")
plt.legend()
plt.grid()
plt.show()


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

z = sp.symbols('z')
j = sp.I
pi = sp.pi

# ======================================================
# 問題25: f(z)=1/(z-1), z=0 まわりでローラン展開
f25 = 1/(z-1)
laurent25 = sp.series(f25, z, 0, 5)
print("P25 Laurent:", laurent25)

# ======================================================
# 問題26: f(z)=1/(z(1-z)), |z|<1 でローラン展開
f26 = 1/(z*(1-z))
laurent26 = sp.series(f26, z, 0, 5)
print("P26 Laurent:", laurent26)

# ======================================================
# 問題27: f(z)=1/z^2 の z=0 における留数
f27 = 1/z**2
res27 = sp.residue(f27, z, 0)
print("P27 Residue at 0:", res27)

# ======================================================
# 問題28: f(z)=exp(z)/z の z=0 における留数
f28 = sp.exp(z)/z
res28 = sp.residue(f28, z, 0)
print("P28 Residue at 0:", res28)

# ======================================================
# 問題29: f(z)=1/(z^2+1) の z=i における留数
f29 = 1/(z**2+1)
res29 = sp.residue(f29, z, sp.I)
print("P29 Residue at i:", res29)

# ======================================================
# 問題30: ∫-∞^∞ 1/(x^2+1) dx
res30 = mp.quad(lambda xx: 1/(xx**2+1), [-mp.inf, mp.inf])
print("P30 Integral:", res30, " exact:", pi)

# ======================================================
# 可視化: 複素平面の極をプロット
poles = [1, 0, sp.I, -sp.I]
labels = ["z=1 (P25)", "z=0 (P26-28)", "z=i (P29)", "z=-i (P30 contour)"]

plt.figure(figsize=(6,6))
for pole, label in zip(poles, labels):
    plt.plot(sp.re(pole), sp.im(pole), 'ro')
    plt.text(sp.re(pole)+0.05, sp.im(pole)+0.05, label, fontsize=9)

# 単位円と大円を描く
theta = np.linspace(0,2*np.pi,400)
plt.plot(np.cos(theta), np.sin(theta), 'b--', label="Unit circle")
plt.plot(2*np.cos(theta), 2*np.sin(theta), 'g--', label="Radius=2 circle")

plt.axhline(0, color='k', linewidth=0.5)
plt.axvline(0, color='k', linewidth=0.5)
plt.gca().set_aspect('equal')
plt.title("Poles and Contours (Problems 25–30)")
plt.xlabel("Re(z)")
plt.ylabel("Im(z)")
plt.legend()
plt.grid()
plt.show()

第3章ラプラス変換 — Python課題

1. ラプラス変換

【問題1】
f(t) = 1 のラプラス変換を Python (SymPy) で計算せよ。

【問題2】
f(t) = e^(-2t) のラプラス変換を計算せよ。

【問題3】
f(t) = sin(3t) のラプラス変換を計算せよ。

【問題4】
f(t) = cos(5t) のラプラス変換を計算せよ。

【問題5】
f(t) = t のラプラス変換を Python で求めよ。

【問題6】
ラプラス変換の線形性を確認せよ。f1(t)=e^(-t), f2(t)=cos(2t) の和を変換して確かめよ。

【問題7】
微分のラプラス変換性質を確認せよ。f(t)=e^(-2t) の場合、f'(t) を変換して比較せよ。

【問題8】
積分のラプラス変換性質を確認せよ。f(t)=sin(t) の積分をラプラス変換で求めよ。

【問題9】
部分分数分解を使って F(s) = (s+2)/(s²+3s+2) を分解せよ。

【問題10】
F(s) = (s+2)/(s²+3s+2) の逆ラプラス変換を Python で求めよ。

【問題11】
F(s) = 1/(s(s+1)) の逆ラプラス変換を求めよ。

【問題12】
常微分方程式 y' + y = 1, y(0)=0 をラプラス変換で解き、Pythonで確認せよ。

2. デルタ関数と線形システム

【問題13】
δ(t) のラプラス変換を Python で確認せよ。

【問題14】
δ(t-a) のラプラス変換を計算せよ。

【問題15】
システム h(t)=e^(-t)u(t) に入力 x(t)=δ(t) を加えたときの出力 y(t) を Python で計算せよ。

【問題16】
入力 x(t)=u(t), インパルス応答 h(t)=e^(-2t)u(t) の畳み込みを Python で求めよ。

【問題17】
入力 x(t)=e^(-t)u(t), h(t)=u(t) の畳み込みを Python で求めよ。

【問題18】
H(s)=1/(s+1) の伝達関数に対して、ステップ入力を与えたときの応答を Python (scipy.signal) で描け。

【問題19】
H(s)=1/(s²+2s+2) の伝達関数に対して、インパルス応答を Python で求め、時間波形を描け。

【問題20】
H(s)=10/(s²+2s+10) のボード線図を Python で描け。

答え


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

t, s = sp.symbols('t s')
y = sp.Function('y')

# ======================================================
# Problem 1
f1 = 1
res1 = sp.laplace_transform(f1, t, s)
print("P1:", res1)

# ======================================================
# Problem 2
f2 = sp.exp(-2*t)
res2 = sp.laplace_transform(f2, t, s)
print("P2:", res2)

# ======================================================
# Problem 3
f3 = sp.sin(3*t)
res3 = sp.laplace_transform(f3, t, s)
print("P3:", res3)

# ======================================================
# Problem 4
f4 = sp.cos(5*t)
res4 = sp.laplace_transform(f4, t, s)
print("P4:", res4)

# ======================================================
# Problem 5
f5 = t
res5 = sp.laplace_transform(f5, t, s)
print("P5:", res5)

# ======================================================
# Problem 6
f6 = sp.exp(-t) + sp.cos(2*t)
res6 = sp.laplace_transform(f6, t, s)
print("P6 (direct):", res6)
res6_sep = sp.laplace_transform(sp.exp(-t), t, s)[0] + sp.laplace_transform(sp.cos(2*t), t, s)[0]
print("P6 (separately):", res6_sep)

# ======================================================
# Problem 7
f7 = sp.exp(-2*t)
f7p = sp.diff(f7, t)
res7_f = sp.laplace_transform(f7, t, s)
res7_fp = sp.laplace_transform(f7p, t, s)
print("P7 f:", res7_f)
print("P7 f':", res7_fp)

# ======================================================
# Problem 8
f8 = sp.sin(t)
F8 = sp.laplace_transform(f8, t, s)[0]
res8 = F8/s
print("P8:", res8)

# ======================================================
# Problem 9
F9 = (s+2)/(s**2+3*s+2)
res9 = sp.apart(F9, s)
print("P9:", res9)

# ======================================================
# Problem 10
res10 = sp.inverse_laplace_transform(F9, s, t)
print("P10:", res10)

# ======================================================
# Problem 11
F11 = 1/(s*(s+1))
res11 = sp.inverse_laplace_transform(F11, s, t)
print("P11:", res11)

# ======================================================
# Problem 12
F12Y = 1/(s*(s+1))  # ODE Laplace form
res12 = sp.inverse_laplace_transform(F12Y, s, t)
print("P12 solution y(t):", res12)

# ======================================================
# Visualization with matplotlib
# Example: visualize Problem 12 solution
f_num = sp.lambdify(t, res12, 'numpy')
T = np.linspace(0, 10, 200)
Y = f_num(T)

plt.figure(figsize=(6,4))
plt.plot(T, Y, label="y(t) from ODE y'+y=1")
plt.xlabel("t")
plt.ylabel("y(t)")
plt.title("Solution of ODE via Laplace Transform (Problem 12)")
plt.legend()
plt.grid(True)
plt.show()


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

t, tau, s, a = sp.symbols('t tau s a', real=True, positive=True)
u = sp.Heaviside

# ======================================================
# Problem 13: δ(t) Laplace
delta = sp.DiracDelta(t)
res13 = sp.laplace_transform(delta, t, s)
print("P13 δ(t):", res13)

# ======================================================
# Problem 14: δ(t-a) Laplace
delta_a = sp.DiracDelta(t-a)
res14 = sp.laplace_transform(delta_a, t, s)
print("P14 δ(t-a):", res14)

# ======================================================
# Problem 15: h(t)=e^-t u(t), x(t)=δ(t)
h15 = sp.exp(-t)*u(t)
x15 = sp.DiracDelta(t)
# 畳み込み: y(t) = ∫ x(τ) h(t-τ) dτ
y15 = sp.integrate(x15.subs(t, tau)*h15.subs(t, t-tau), (tau, -sp.oo, sp.oo))
print("P15 y(t):", sp.simplify(y15))

# ======================================================
# Problem 16: x=u(t), h=e^-2t u(t)
x16 = u(tau)
h16 = sp.exp(-2*(t-tau))*u(t-tau)
y16 = sp.integrate(x16*h16, (tau, 0, t))
print("P16 y(t):", sp.simplify(y16))

# ======================================================
# Problem 17: x=e^-t u(t), h=u(t)
x17 = sp.exp(-tau)*u(tau)
h17 = u(t-tau)
y17 = sp.integrate(x17*h17, (tau, 0, t))
print("P17 y(t):", sp.simplify(y17))

# ======================================================
# Problem 18: H(s)=1/(s+1), step input
num18, den18 = [1], [1, 1]
system18 = signal.TransferFunction(num18, den18)
T = np.linspace(0, 10, 200)
T18, Y18 = signal.step(system18, T=T)

plt.figure()
plt.plot(T18, Y18)
plt.title("Problem 18: Step Response H(s)=1/(s+1)")
plt.xlabel("t"); plt.ylabel("y(t)")
plt.grid(True)
plt.show()

# ======================================================
# Problem 19: H(s)=1/(s^2+2s+2), impulse response
num19, den19 = [1], [1, 2, 2]
system19 = signal.TransferFunction(num19, den19)
T = np.linspace(0, 10, 200)
T19, Y19 = signal.impulse(system19, T=T)

plt.figure()
plt.plot(T19, Y19)
plt.title("Problem 19: Impulse Response H(s)=1/(s^2+2s+2)")
plt.xlabel("t"); plt.ylabel("h(t)")
plt.grid(True)
plt.show()

# ======================================================
# Problem 20: H(s)=10/(s^2+2s+10), Bode plot
num20, den20 = [10], [1, 2, 10]
system20 = signal.TransferFunction(num20, den20)

w, mag, phase = signal.bode(system20)

plt.figure()
plt.semilogx(w, mag)
plt.title("Problem 20: Bode Magnitude")
plt.xlabel("Frequency [rad/s]"); plt.ylabel("Magnitude [dB]")
plt.grid(True)

plt.figure()
plt.semilogx(w, phase)
plt.title("Problem 20: Bode Phase")
plt.xlabel("Frequency [rad/s]"); plt.ylabel("Phase [deg]")
plt.grid(True)
plt.show()

第4章フーリエ級数とフーリエ変換 — Python課題

1. フーリエ級数

【問題1】
周期関数 f(t)=1 （周期 2π）のフーリエ級数展開係数を Python で計算せよ。

【問題2】
f(t)=t （-π < t < π, 周期 2π）のフーリエ級数展開を求めよ。

【問題3】
f(t)=|t| （-π < t < π, 周期 2π）のフーリエ級数展開を求めよ。

【問題4】
偶関数 f(t)=cos(t) のフーリエ級数係数を Python で確認せよ。

【問題5】
奇関数 f(t)=sin(t) のフーリエ級数係数を Python で確認せよ。

【問題6】
周期矩形波をフーリエ級数展開し、部分和を Python で描いてギブス現象を観察せよ。

【問題7】
三角波をフーリエ級数展開し、10項までの近似波形を描け。

【問題8】
のこぎり波をフーリエ級数展開し、Pythonで再構成せよ。

【問題9】
矩形波の基本周波数と高調波成分のエネルギー比を計算せよ。

【問題10】
フーリエ級数展開による Parseval の等式を矩形波で確認せよ。

2. フーリエ変換

【問題11】
f(t)=e^(-t²) （ガウス関数）のフーリエ変換を Python で求めよ。

【問題12】
δ(t) のフーリエ変換を計算し、結果を確認せよ。

【問題13】
矩形関数 rect(t) のフーリエ変換が sinc 関数になることを Python で確認せよ。

【問題14】
三角関数 cos(2πft) のフーリエ変換を Python で求めよ。

【問題15】
畳み込み定理を確認せよ。f1(t)=rect(t), f2(t)=rect(t) の畳み込みを Python で計算し、フーリエ変換と比較せよ。

【問題16】
f(t)=sin(t)/t のフーリエ変換を Python で近似計算せよ。

【問題17】
時間シフト特性を確認せよ。f(t)=rect(t) をシフトした場合のフーリエ変換を求めよ。

【問題18】
周波数シフト特性を確認せよ。f(t)=rect(t) に cos 波を掛けた場合のスペクトルを求めよ。

【問題19】
ガウス関数の自己フーリエ性を Python で確認せよ。

【問題20】
ローパスフィルタ H(f)=rect(f/2) を定義し、入力信号のスペクトルに掛けてフィルタリングを実装せよ。

3. FFT（高速フーリエ変換）

【問題21】
離散フーリエ変換（DFT）の定義式を Python で実装せよ。

【問題22】
Python の numpy.fft を用いて DFT を計算し、手実装の結果と比較せよ。

【問題23】
サンプリング周波数 fs=100Hz、f=5Hz の正弦波をFFTで解析し、ピークを確認せよ。

【問題24】
矩形波をサンプリングし、FFTを実行して高調波スペクトルを観察せよ。

【問題25】
サンプリング周波数 fs=100Hz で f=60Hz の信号をFFTすると折り返し雑音が発生することを確認せよ。

【問題26】
サンプリングとナイキスト周波数の関係を、f=10Hz, 20Hz, 40Hz, fs=60Hz のケースで確認せよ。

【問題27】
ノイズを加えた正弦波をFFTで解析し、S/N比を評価せよ。

【問題28】
音声信号（wavファイル）を読み込み、FFTでスペクトルを描け。

【問題29】
2D FFTを用いて白黒画像の周波数成分を解析し、ローパスフィルタを適用せよ。

【問題30】
QPSK信号をPythonで生成し、FFTでスペクトルを観察せよ。

答え


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

t, n = sp.symbols('t n', real=True, integer=True)

# ======================================================
# Problem 1: f(t)=1
f1 = 1
a0 = (1/sp.pi)*sp.integrate(f1, (t, -sp.pi, sp.pi))
an = (1/sp.pi)*sp.integrate(f1*sp.cos(n*t), (t, -sp.pi, sp.pi))
bn = (1/sp.pi)*sp.integrate(f1*sp.sin(n*t), (t, -sp.pi, sp.pi))
print("P1 Fourier coefficients: a0=", a0, ", an=", sp.simplify(an), ", bn=", sp.simplify(bn))

# ======================================================
# Problem 2: f(t)=t
f2 = t
a0_2 = (1/sp.pi)*sp.integrate(f2, (t, -sp.pi, sp.pi))
an_2 = (1/sp.pi)*sp.integrate(f2*sp.cos(n*t), (t, -sp.pi, sp.pi))
bn_2 = (1/sp.pi)*sp.integrate(f2*sp.sin(n*t), (t, -sp.pi, sp.pi))
print("P2 Fourier coeffs: a0=", a0_2, "an(n)=", sp.simplify(an_2), "bn(n)=", sp.simplify(bn_2))

# ======================================================
# Problem 3: f(t)=|t|
f3 = sp.Abs(t)
a0_3 = (1/sp.pi)*sp.integrate(f3, (t, -sp.pi, sp.pi))
an_3 = (1/sp.pi)*sp.integrate(f3*sp.cos(n*t), (t, -sp.pi, sp.pi))
bn_3 = (1/sp.pi)*sp.integrate(f3*sp.sin(n*t), (t, -sp.pi, sp.pi))
print("P3 Fourier coeffs: a0=", a0_3, "an(n)=", sp.simplify(an_3), "bn(n)=", sp.simplify(bn_3))

# ======================================================
# Problem 4: f(t)=cos(t)
f4 = sp.cos(t)
a0_4 = (1/sp.pi)*sp.integrate(f4, (t, -sp.pi, sp.pi))
an_4 = (1/sp.pi)*sp.integrate(f4*sp.cos(n*t), (t, -sp.pi, sp.pi))
bn_4 = (1/sp.pi)*sp.integrate(f4*sp.sin(n*t), (t, -sp.pi, sp.pi))
print("P4 Fourier coeffs cos(t): a0=", a0_4, "an=", sp.simplify(an_4), "bn=", sp.simplify(bn_4))

# ======================================================
# Problem 5: f(t)=sin(t)
f5 = sp.sin(t)
a0_5 = (1/sp.pi)*sp.integrate(f5, (t, -sp.pi, sp.pi))
an_5 = (1/sp.pi)*sp.integrate(f5*sp.cos(n*t), (t, -sp.pi, sp.pi))
bn_5 = (1/sp.pi)*sp.integrate(f5*sp.sin(n*t), (t, -sp.pi, sp.pi))
print("P5 Fourier coeffs sin(t): a0=", a0_5, "an=", sp.simplify(an_5), "bn=", sp.simplify(bn_5))

# ======================================================
# Problem 6: Rectangular wave and Gibbs phenomenon
N = 20
T = np.linspace(-2*np.pi, 2*np.pi, 1000)
f_rect = np.sign(np.sin(T))  # 矩形波

def fourier_rect(T, N):
    s = np.zeros_like(T)
    for k in range(1, N*2, 2):  # odd harmonics
        s += (4/np.pi)*(1/k)*np.sin(k*T)
    return s

plt.figure()
plt.plot(T, f_rect, 'k', label="Original square wave")
plt.plot(T, fourier_rect(T, 5), label="N=5")
plt.plot(T, fourier_rect(T, 20), label="N=20")
plt.title("P6 Gibbs phenomenon in Square Wave")
plt.legend(); plt.grid(True)
plt.show()

# ======================================================
# Problem 7: Triangular wave
def fourier_triangle(T, N):
    s = np.zeros_like(T)
    for k in range(1, N*2, 2):
        s += (8/(np.pi**2))*(1/(k**2))*(-1)**((k-1)//2)*np.cos(k*T)
    return s

plt.figure()
plt.plot(T, fourier_triangle(T, 10), label="Triangular wave approx N=10")
plt.title("P7 Triangular Wave Fourier Approximation")
plt.grid(True); plt.legend()
plt.show()

# ======================================================
# Problem 8: Sawtooth wave
def fourier_sawtooth(T, N):
    s = np.zeros_like(T)
    for k in range(1, N+1):
        s += (-2/np.pi)*(1/k)*np.sin(k*T)
    return s

plt.figure()
plt.plot(T, fourier_sawtooth(T, 20), label="Sawtooth wave approx N=20")
plt.title("P8 Sawtooth Wave Fourier Approximation")
plt.grid(True); plt.legend()
plt.show()

# ======================================================
# Problem 9: Energy ratio square wave
Nmax = 20
energies = [(4/(np.pi*k))**2/2 for k in range(1, Nmax*2, 2)]
fundamental = energies[0]
harmonics = sum(energies[1:])
print("P9 Energy ratio fundamental/harmonics =", fundamental, "/", harmonics)

# ======================================================
# Problem 10: Parseval's theorem for square wave
total_energy = sum(energies)
parseval = (1/(2*np.pi))*np.trapz(f_rect**2, T)
print("P10 Parseval check: Fourier sum =", total_energy, "Integral =", parseval)


import numpy as np
import matplotlib.pyplot as plt

# ========================
# Parameters
# ========================
fs = 200          # Sampling frequency [Hz]
N = 2048          # Number of points
t = np.linspace(-5, 5, N)
freq = np.fft.fftshift(np.fft.fftfreq(N, d=1/fs))

# Helper: FFT
def compute_fft(sig):
    return np.fft.fftshift(np.fft.fft(sig)), np.fft.fftshift(np.fft.fftfreq(len(sig), d=1/fs))

# ========================
# Problem 11 Gaussian
# ========================
f11 = np.exp(-t**2)
F11, f_axis = compute_fft(f11)

# ========================
# Problem 12 Delta (approx as impulse)
# ========================
f12 = np.zeros_like(t)
f12[N//2] = 1
F12, _ = compute_fft(f12)

# ========================
# Problem 13 Rect -> sinc
# ========================
f13 = np.where(np.abs(t) < 0.5, 1, 0)
F13, _ = compute_fft(f13)

# ========================
# Problem 14 Cosine
# ========================
f14 = np.cos(2*np.pi*5*t)
F14, _ = compute_fft(f14)

# ========================
# Problem 15 Convolution Rect * Rect
# ========================
conv15 = np.convolve(f13, f13, mode="same")
F15, _ = compute_fft(conv15)

# ========================
# Problem 16 sinc(t)
# ========================
f16 = np.sinc(t)
F16, _ = compute_fft(f16)

# ========================
# Problem 17 Time Shift
# ========================
f17 = np.where(np.abs(t-2) < 0.5, 1, 0)
F17, _ = compute_fft(f17)

# ========================
# Problem 18 Frequency Shift
# ========================
f18 = f13 * np.cos(2*np.pi*10*t)
F18, _ = compute_fft(f18)

# ========================
# Problem 19 Gaussian self-Fourier
# ========================
f19 = np.exp(-t**2)
F19, _ = compute_fft(f19)

# ========================
# Problem 20 Low-pass filter
# ========================
F20in, _ = compute_fft(f14)
H20 = np.where(np.abs(freq) < 5, 1, 0)
F20out = F20in * H20
f20out = np.fft.ifft(np.fft.ifftshift(F20out))

# ========================
# Plot Results
# ========================
plt.figure(figsize=(12,18))

# Problem 11 Gaussian
plt.subplot(5,2,1)
plt.plot(t, f11)
plt.title("Problem 11: Gaussian f(t)")
plt.subplot(5,2,2)
plt.plot(f_axis, np.abs(F11))
plt.title("Gaussian Spectrum")

# Problem 12 Delta
plt.subplot(5,2,3)
plt.plot(t[::100], f12[::100], "o")
plt.title("Problem 12: Delta approx")
plt.subplot(5,2,4)
plt.plot(f_axis, np.abs(F12))
plt.title("Delta Spectrum")

# Problem 13 Rect
plt.subplot(5,2,5)
plt.plot(t, f13)
plt.title("Problem 13: Rect(t)")
plt.subplot(5,2,6)
plt.plot(f_axis, np.abs(F13))
plt.title("Rect Spectrum (sinc)")

# Problem 14 Cosine
plt.subplot(5,2,7)
plt.plot(t, f14)
plt.title("Problem 14: cos(2πft)")
plt.subplot(5,2,8)
plt.plot(f_axis, np.abs(F14))
plt.title("Cosine Spectrum")

# Problem 15 Convolution
plt.subplot(5,2,9)
plt.plot(conv15)
plt.title("Problem 15: rect*rect (Triangle)")
plt.subplot(5,2,10)
plt.plot(f_axis, np.abs(F15))
plt.title("Convolution Spectrum")

plt.tight_layout()
plt.show()

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

# Problem 16 sinc
plt.subplot(5,2,1)
plt.plot(t, f16)
plt.title("Problem 16: sinc(t)")
plt.subplot(5,2,2)
plt.plot(f_axis, np.abs(F16))
plt.title("sinc Spectrum")

# Problem 17 Shift
plt.subplot(5,2,3)
plt.plot(t, f17)
plt.title("Problem 17: Shifted rect")
plt.subplot(5,2,4)
plt.plot(f_axis, np.abs(F17))
plt.title("Shifted Spectrum")

# Problem 18 Frequency Shift
plt.subplot(5,2,5)
plt.plot(t, f18)
plt.title("Problem 18: rect * cos")
plt.subplot(5,2,6)
plt.plot(f_axis, np.abs(F18))
plt.title("Frequency Shift Spectrum")

# Problem 19 Gaussian Self
plt.subplot(5,2,7)
plt.plot(t, f19)
plt.title("Problem 19: Gaussian")
plt.subplot(5,2,8)
plt.plot(f_axis, np.abs(F19))
plt.title("Gaussian Self-Spectrum")

# Problem 20 Filtered cosine
plt.subplot(5,2,9)
plt.plot(t, f14, label="Input cos")
plt.plot(t, np.real(f20out), label="Filtered")
plt.legend()
plt.title("Problem 20: Filtering")

plt.tight_layout()
plt.show()


import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile

# ========================
# Problem 21: Manual DFT
# ========================
def dft(x):
    N = len(x)
    n = np.arange(N)
    k = n.reshape((N,1))
    M = np.exp(-2j*np.pi*k*n/N)
    return np.dot(M, x)

# ========================
# Problem 22: Compare manual DFT with numpy FFT
# ========================
x22 = np.array([1,2,3,4])
X22_manual = dft(x22)
X22_np = np.fft.fft(x22)

# ========================
# Problem 23: FFT of sine
# ========================
fs = 100
N = 256
t = np.arange(N)/fs
f23 = np.sin(2*np.pi*5*t)
F23 = np.fft.fft(f23)
freq23 = np.fft.fftfreq(N, d=1/fs)

# ========================
# Problem 24: FFT of square wave
# ========================
f24 = np.sign(np.sin(2*np.pi*5*t))
F24 = np.fft.fft(f24)
freq24 = np.fft.fftfreq(N, d=1/fs)

# ========================
# Problem 25: Aliasing
# ========================
fs25 = 100
t25 = np.arange(N)/fs25
f25 = np.sin(2*np.pi*60*t25)  # fs/2=50Hz, so aliasing
F25 = np.fft.fft(f25)
freq25 = np.fft.fftfreq(N, d=1/fs25)

# ========================
# Problem 26: Nyquist cases
# ========================
fs26 = 60
freqs = [10,20,40]
signals26 = [np.sin(2*np.pi*f*t25) for f in freqs]
spectra26 = [np.fft.fft(sig) for sig in signals26]
freq26 = np.fft.fftfreq(N, d=1/fs26)

# ========================
# Problem 27: Sine + noise
# ========================
f27 = np.sin(2*np.pi*5*t) + 0.5*np.random.randn(len(t))
F27 = np.fft.fft(f27)
freq27 = np.fft.fftfreq(N, d=1/fs)

# ========================
# Problem 28: Audio FFT (dummy sine if no wav)
# ========================
try:
    rate, data = wavfile.read("test.wav")
    if data.ndim > 1:
        data = data[:,0]
    N28 = 2048
    F28 = np.fft.fft(data[:N28])
    freq28 = np.fft.fftfreq(N28, d=1/rate)
except:
    rate = 8000
    t28 = np.arange(2048)/rate
    data = np.sin(2*np.pi*440*t28)
    F28 = np.fft.fft(data)
    freq28 = np.fft.fftfreq(len(data), d=1/rate)

# ========================
# Problem 29: 2D FFT image
# ========================
img = np.zeros((64,64))
img[20:40,20:40] = 1
F29 = np.fft.fftshift(np.fft.fft2(img))
lowpass = np.zeros_like(img)
lowpass[28:36,28:36] = 1
img_lp = np.fft.ifft2(np.fft.ifftshift(F29*lowpass))

# ========================
# Problem 30: QPSK
# ========================
Nsym = 128
bits = np.random.randint(0,2,2*Nsym)
symbols = (2*bits[0::2]-1) + 1j*(2*bits[1::2]-1)
f30 = np.repeat(symbols, 16)
F30 = np.fft.fftshift(np.fft.fft(f30))
freq30 = np.fft.fftshift(np.fft.fftfreq(len(f30), d=1/fs))

# ========================
# Plotting
# ========================
plt.figure(figsize=(12,18))

# Problem 22 comparison
plt.subplot(5,2,1)
plt.stem(np.abs(X22_manual), linefmt="b-", markerfmt="bo")
plt.title("Problem 22: Manual DFT Magnitude")
plt.subplot(5,2,2)
plt.stem(np.abs(X22_np), linefmt="r-", markerfmt="ro")
plt.title("Problem 22: numpy FFT Magnitude")

# Problem 23
plt.subplot(5,2,3)
plt.plot(t, f23)
plt.title("Problem 23: Sine 5Hz")
plt.subplot(5,2,4)
plt.plot(freq23, np.abs(F23))
plt.title("FFT Spectrum")

# Problem 24
plt.subplot(5,2,5)
plt.plot(t, f24)
plt.title("Problem 24: Square Wave")
plt.subplot(5,2,6)
plt.plot(freq24, np.abs(F24))
plt.title("Square Wave Spectrum")

# Problem 25
plt.subplot(5,2,7)
plt.plot(t25, f25)
plt.title("Problem 25: Aliasing 60Hz sampled 100Hz")
plt.subplot(5,2,8)
plt.plot(freq25, np.abs(F25))
plt.title("Aliased Spectrum")

# Problem 26
plt.subplot(5,2,9)
for i, sig in enumerate(signals26):
    plt.plot(t25[:50], sig[:50], label=f"{freqs[i]}Hz")
plt.legend()
plt.title("Problem 26: Signals near Nyquist")
plt.subplot(5,2,10)
for i, spec in enumerate(spectra26):
    plt.plot(freq26, np.abs(spec), label=f"{freqs[i]}Hz")
plt.legend()
plt.title("Nyquist FFT Spectra")

plt.tight_layout()
plt.show()

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

# Problem 27
plt.subplot(5,2,1)
plt.plot(t, f27)
plt.title("Problem 27: Sine + Noise")
plt.subplot(5,2,2)
plt.plot(freq27, np.abs(F27))
plt.title("FFT with Noise")

# Problem 28
plt.subplot(5,2,3)
plt.plot(data[:500])
plt.title("Problem 28: Audio Signal (or dummy)")
plt.subplot(5,2,4)
plt.plot(freq28, np.abs(F28))
plt.title("Audio FFT Spectrum")

# Problem 29
plt.subplot(5,2,5)
plt.imshow(img, cmap="gray")
plt.title("Problem 29: Image")
plt.subplot(5,2,6)
plt.imshow(np.log(np.abs(F29)+1), cmap="inferno")
plt.title("Image FFT Spectrum")
plt.subplot(5,2,7)
plt.imshow(np.real(img_lp), cmap="gray")
plt.title("Low-pass filtered Image")

# Problem 30
plt.subplot(5,2,8)
plt.plot(np.real(f30[:200]), label="I")
plt.plot(np.imag(f30[:200]), label="Q")
plt.legend()
plt.title("Problem 30: QPSK Signal")
plt.subplot(5,2,9)
plt.plot(freq30, np.abs(F30))
plt.title("QPSK Spectrum")

plt.tight_layout()
plt.show()

第5章微分方程式 — Python課題

1. 微分方程式の基礎

【問題1】
dy/dx = 3x² を Python (SymPy) で解き、一般解を求めよ。

【問題2】
dy/dx = y を Python で解き、解曲線を 0 ≤ x ≤ 5 でプロットせよ。

【問題3】
dy/dx = sin(x) を解き、初期条件 y(0)=0 の解をプロットせよ。

【問題4】
d²y/dx² = 0 を Python で解き、一般解を確認せよ。

【問題5】
常微分方程式と偏微分方程式の違いを例として、dy/dx = y と ∂u/∂t = ∂²u/∂x² を Python で数式表示せよ。

2. 1階常微分方程式

【問題6】
dy/dx = x y を変数分離法で解き、一般解を Python で計算せよ。

【問題7】
dy/dx = (x²+y²)/(xy) を Python で同次形変換し、一般解を求めよ。

【問題8】
dy/dx + y = e^x を Python で解き、初期条件 y(0)=1 を満たす解を描け。

【問題9】
dy/dx + 2y = cos(x) を Python で解き、特解をプロットせよ。

【問題10】
dy/dx = y(1−y) を Python で解き、ロジスティック曲線を描け。

3. 2階常微分方程式

【問題11】
y'' + y = 0 を Python で解き、解が sin, cos の組み合わせになることを確認せよ。

【問題12】
y'' − 4y = 0 を解き、指数関数の組み合わせになることを確認せよ。

【問題13】
y'' + y' − 6y = 0 を解き、特性方程式の解を確認せよ。

【問題14】
y'' + y = sin(x) を解き、特解と一般解を Python で確認せよ。

【問題15】
y'' + 2y' + y = 0 を解き、重解の場合の解形式を確認せよ。

【問題16】
y'' + 4y = cos(2x) を解き、共振の特解を Python で求めよ。

【問題17】
y'' − y = e^x を解き、特解が C e^x + Ax e^x になることを確認せよ。

【問題18】
質点の単振動 y'' + ω² y = 0 を ω=2 の場合に数値的に解いてプロットせよ。

【問題19】
減衰振動 y'' + 0.5y' + 4y = 0 を数値的に解き、減衰を確認せよ。

【問題20】
強制振動 y'' + y = cos(2x) を数値解法で解き、時間応答をプロットせよ。

4. 高階常微分方程式

【問題21】
y''' = 0 を Python で解き、一般解を求めよ。

【問題22】
y''' − y = 0 を Python で解き、特性方程式の3つの解を確認せよ。

【問題23】
y^(4) + y = 0 を Python で解き、sin, cos, sinh, cosh が混じる解を確認せよ。

【問題24】
y^(4) − 16y = 0 を Python で解き、指数関数の形で表せ。

【問題25】
y''' + y'' − y' − y = 0 を Python で解き、複素解を含む一般解を確認せよ。

5. 連立微分方程式

【問題26】
dx/dt = x, dy/dt = −y を Python で数値的に解き、位相図を描け。

【問題27】
dx/dt = y, dy/dt = −x を Python で解き、円運動になることを確認せよ。

【問題28】
dx/dt = x+y, dy/dt = x−y を解き、行列を使って一般解を求めよ。

【問題29】
dx/dt = 3x + 4y, dy/dt = −4x + 3y の連立方程式を Python で解き、回転行列の性質を確認せよ。

【問題30】
人口モデル dx/dt = 0.1x − 0.01xy, dy/dt = −0.1y + 0.01xy を数値的に解き、捕食-被食モデルを描け。

6. 変数係数・非線形微分方程式

【問題31】
オイラー方程式 x² y'' + x y' − y = 0 を Python で解け。

【問題32】
x² y'' − 2x y' + 2y = 0 を解き、解形式を確認せよ。

【問題33】
非線形微分方程式 dy/dx = y² を Python で解き、爆発的解を確認せよ。

【問題34】
dy/dx = sin(y) を Python で解き、数値的にプロットせよ。

【問題35】
dy/dx = y² − x を数値的に解き、x=0～2 の範囲で解をプロットせよ。

7. 偏微分方程式（入門）

【問題36】
1次元熱伝導方程式 ∂u/∂t = ∂²u/∂x² を Python で数値的に解け（初期条件: u(x,0)=sin(πx), 区間 [0,1]）。

【問題37】
波動方程式 ∂²u/∂t² = ∂²u/∂x² を Python で数値的に解け。

【問題38】
ラプラス方程式 ∂²u/∂x² + ∂²u/∂y² = 0 を Python で格子点上で数値的に解け。

【問題39】
拡散方程式 ∂u/∂t = D ∂²u/∂x² を Python で数値シミュレーションせよ（D=0.1）。

【問題40】
波動方程式を有限差分法で解き、t=0, 0.5, 1.0 の解をプロットせよ。

答え


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

# Symbolic variables
x = sp.Symbol('x')
y = sp.Function('y')(x)

# 追加: 積分定数を明示的にシンボルとして定義
C1, C2 = sp.symbols('C1 C2')

# ========================
# Problem 1: dy/dx = 3x^2
# ========================
ode1 = sp.Eq(sp.diff(y, x), 3*x**2)
sol1 = sp.dsolve(ode1)

# Problem 2: dy/dx = y
ode2 = sp.Eq(sp.diff(y, x), y)
sol2 = sp.dsolve(ode2)

# Problem 3: dy/dx = sin(x), y(0)=0
ode3 = sp.Eq(sp.diff(y, x), sp.sin(x))
sol3 = sp.dsolve(ode3, ics={y.subs(x,0):0})

# Problem 4: y''=0
ode4 = sp.Eq(sp.diff(y, x, 2), 0)
sol4 = sp.dsolve(ode4)

# Problem 5: ODE vs PDE (symbolic display)
t = sp.Symbol('t')
u = sp.Function('u')(x,t)
ode5 = sp.Eq(sp.diff(y, x), y)
pde5 = sp.Eq(sp.diff(u,t), sp.diff(u,x,2))

# ========================
# Problem 6: dy/dx = x y
# ========================
ode6 = sp.Eq(sp.diff(y, x), x*y)
sol6 = sp.dsolve(ode6)

# Problem 7: dy/dx = (x^2+y^2)/(x y)
ode7 = sp.Eq(sp.diff(y, x), (x**2+y**2)/(x*y))
sol7 = sp.dsolve(ode7)

# Problem 8: dy/dx + y = e^x, y(0)=1
ode8 = sp.Eq(sp.diff(y, x) + y, sp.exp(x))
sol8 = sp.dsolve(ode8, ics={y.subs(x,0):1})

# Problem 9: dy/dx + 2y = cos(x)
ode9 = sp.Eq(sp.diff(y, x) + 2*y, sp.cos(x))
sol9 = sp.dsolve(ode9)

# Problem 10: dy/dx = y(1-y) (Logistic)
ode10 = sp.Eq(sp.diff(y, x), y*(1-y))
sol10 = sp.dsolve(ode10)

# ========================
# Problem 11: y''+y=0
# ========================
ode11 = sp.Eq(sp.diff(y, x, 2)+y,0)
sol11 = sp.dsolve(ode11)

# Problem 12: y''-4y=0
ode12 = sp.Eq(sp.diff(y, x, 2)-4*y,0)
sol12 = sp.dsolve(ode12)

# Problem 13: y''+y'-6y=0
ode13 = sp.Eq(sp.diff(y, x, 2)+sp.diff(y, x)-6*y,0)
sol13 = sp.dsolve(ode13)

# Problem 14: y''+y=sin(x)
ode14 = sp.Eq(sp.diff(y, x, 2)+y, sp.sin(x))
sol14 = sp.dsolve(ode14)

# Problem 15: y''+2y'+y=0
ode15 = sp.Eq(sp.diff(y, x, 2)+2*sp.diff(y, x)+y,0)
sol15 = sp.dsolve(ode15)

# Problem 16: y''+4y=cos(2x)
ode16 = sp.Eq(sp.diff(y, x, 2)+4*y, sp.cos(2*x))
sol16 = sp.dsolve(ode16)

# Problem 17: y''-y=e^x
ode17 = sp.Eq(sp.diff(y, x, 2)-y, sp.exp(x))
sol17 = sp.dsolve(ode17)

# Problem 18: y''+ω^2 y=0 with ω=2
def prob18(t,y): return [y[1], -4*y[0]]
sol18 = solve_ivp(prob18, [0,10], [1,0], t_eval=np.linspace(0,10,200))

# Problem 19: y''+0.5y'+4y=0
def prob19(t,y): return [y[1], -0.5*y[1]-4*y[0]]
sol19 = solve_ivp(prob19, [0,20], [1,0], t_eval=np.linspace(0,20,200))

# Problem 20: y''+y=cos(2x)
def prob20(t,y): return [y[1], -y[0]+np.cos(2*t)]
sol20 = solve_ivp(prob20, [0,20], [0,0], t_eval=np.linspace(0,20,200))

# ========================
# Plotting
# ========================
plt.figure(figsize=(12,18))

# Problem 2
xvals = np.linspace(0,5,100)
f2 = [float(sol2.rhs.subs({x:val,C1:1})) for val in xvals]
plt.subplot(5,2,1)
plt.plot(xvals, f2)
plt.title("Problem 2: dy/dx=y")

# Problem 3
f3 = [float(sol3.rhs.subs(x,val)) for val in xvals]
plt.subplot(5,2,2)
plt.plot(xvals, f3)
plt.title("Problem 3: dy/dx=sin(x), y(0)=0")

# Problem 8
f8 = [float(sol8.rhs.subs(x,val)) for val in xvals]
plt.subplot(5,2,3)
plt.plot(xvals, f8)
plt.title("Problem 8: dy/dx+y=e^x")

# Problem 9
f9 = [float(sol9.rhs.subs({x:val,C1:0})) for val in xvals]
plt.subplot(5,2,4)
plt.plot(xvals, f9)
plt.title("Problem 9: dy/dx+2y=cos(x)")

# Problem 10 Logistic
f10 = [float(sol10.rhs.subs({x:val,C1:1})) for val in xvals]
plt.subplot(5,2,5)
plt.plot(xvals, f10)
plt.title("Problem 10: Logistic Equation")

# Problem 11
f11 = [float(sol11.rhs.subs({x:val,C1:1,C2:0})) for val in xvals]
plt.subplot(5,2,6)
plt.plot(xvals, f11)
plt.title("Problem 11: y''+y=0")

# Problem 14
f14 = [float(sol14.rhs.subs({x:val,C1:0,C2:0})) for val in xvals]
plt.subplot(5,2,7)
plt.plot(xvals, f14)
plt.title("Problem 14: y''+y=sin(x)")

# Problem 15
f15 = [float(sol15.rhs.subs({x:val,C1:1,C2:0})) for val in xvals]
plt.subplot(5,2,8)
plt.plot(xvals, f15)
plt.title("Problem 15: y''+2y'+y=0")

# Problem 18
plt.subplot(5,2,9)
plt.plot(sol18.t, sol18.y[0])
plt.title("Problem 18: Harmonic Oscillator ω=2")

# Problem 19
plt.subplot(5,2,10)
plt.plot(sol19.t, sol19.y[0])
plt.title("Problem 19: Damped Oscillator")

plt.tight_layout()
plt.show()

plt.figure(figsize=(10,4))
plt.plot(sol20.t, sol20.y[0])
plt.title("Problem 20: Forced Oscillator y''+y=cos(2x)")
plt.show()


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

# ========================
# Symbolic definitions
# ========================
x = sp.Symbol('x')
y = sp.Function('y')(x)
C1, C2, C3, C4 = sp.symbols('C1 C2 C3 C4')

# ========================
# Problem 21: y'''=0
# ========================
ode21 = sp.Eq(sp.diff(y,x,3), 0)
sol21 = sp.dsolve(ode21)

# Problem 22: y'''-y=0
ode22 = sp.Eq(sp.diff(y,x,3)-y,0)
sol22 = sp.dsolve(ode22)

# Problem 23: y^(4)+y=0
ode23 = sp.Eq(sp.diff(y,x,4)+y,0)
sol23 = sp.dsolve(ode23)

# Problem 24: y^(4)-16y=0
ode24 = sp.Eq(sp.diff(y,x,4)-16*y,0)
sol24 = sp.dsolve(ode24)

# Problem 25: y'''+y''-y'-y=0
ode25 = sp.Eq(sp.diff(y,x,3)+sp.diff(y,x,2)-sp.diff(y,x)-y,0)
sol25 = sp.dsolve(ode25)

# ========================
# System ODEs
# ========================
# Problem 26: dx/dt = x, dy/dt = -y
def sys26(t,vars): return [vars[0], -vars[1]]
sol26 = solve_ivp(sys26, [0,5], [1,1], t_eval=np.linspace(0,5,200))

# Problem 27: dx/dt = y, dy/dt = -x (circle)
def sys27(t,vars): return [vars[1], -vars[0]]
sol27 = solve_ivp(sys27, [0,10], [1,0], t_eval=np.linspace(0,10,400))

# Problem 28: dx/dt = x+y, dy/dt = x-y
A28 = np.array([[1,1],[1,-1]])
def sys28(t,vars): return A28 @ vars
sol28 = solve_ivp(sys28, [0,5], [1,1], t_eval=np.linspace(0,5,200))

# Problem 29: dx/dt = 3x+4y, dy/dt = -4x+3y
A29 = np.array([[3,4],[-4,3]])
def sys29(t,vars): return A29 @ vars
sol29 = solve_ivp(sys29, [0,5], [1,0], t_eval=np.linspace(0,5,200))

# Problem 30: Predator-prey model
def sys30(t,vars):
    x,y = vars
    return [0.1*x - 0.01*x*y, -0.1*y + 0.01*x*y]
sol30 = solve_ivp(sys30, [0,200], [40,9], t_eval=np.linspace(0,200,1000))

# ========================
# Plotting
# ========================
plt.figure(figsize=(12,18))

# Problem 26 phase portrait
plt.subplot(3,2,1)
plt.plot(sol26.y[0], sol26.y[1])
plt.xlabel("x"); plt.ylabel("y"); plt.title("Problem 26: dx/dt=x, dy/dt=-y")

# Problem 27 circle
plt.subplot(3,2,2)
plt.plot(sol27.y[0], sol27.y[1])
plt.xlabel("x"); plt.ylabel("y"); plt.title("Problem 27: Circle Motion")

# Problem 28 trajectory
plt.subplot(3,2,3)
plt.plot(sol28.t, sol28.y[0], label="x")
plt.plot(sol28.t, sol28.y[1], label="y")
plt.legend(); plt.title("Problem 28: dx/dt=x+y, dy/dt=x-y")

# Problem 29 rotation-like
plt.subplot(3,2,4)
plt.plot(sol29.y[0], sol29.y[1])
plt.xlabel("x"); plt.ylabel("y"); plt.title("Problem 29: Rotation System")

# Problem 30 predator-prey trajectory
plt.subplot(3,2,5)
plt.plot(sol30.t, sol30.y[0], label="Prey x")
plt.plot(sol30.t, sol30.y[1], label="Predator y")
plt.legend(); plt.title("Problem 30: Predator-Prey Model")

plt.subplot(3,2,6)
plt.plot(sol30.y[0], sol30.y[1])
plt.xlabel("Prey"); plt.ylabel("Predator")
plt.title("Problem 30: Phase Plane")

plt.tight_layout()
plt.show()

# ========================
# Print symbolic results
# ========================
print("Problem 21:", sol21)
print("Problem 22:", sol22)
print("Problem 23:", sol23)
print("Problem 24:", sol24)
print("Problem 25:", sol25)


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

# ========================
# Symbolic ODEs (31-34)
# ========================
x = sp.Symbol('x')
y = sp.Function('y')(x)

# Problem 31 Euler eq: x^2 y'' + x y' - y = 0
ode31 = sp.Eq(x**2*sp.diff(y,x,2) + x*sp.diff(y,x) - y, 0)
sol31 = sp.dsolve(ode31)

# Problem 32: x^2 y'' - 2x y' + 2y = 0
ode32 = sp.Eq(x**2*sp.diff(y,x,2) - 2*x*sp.diff(y,x) + 2*y, 0)
sol32 = sp.dsolve(ode32)

# Problem 33: dy/dx = y^2
y33 = sp.Function('y')(x)
ode33 = sp.Eq(sp.diff(y33,x), y33**2)
sol33 = sp.dsolve(ode33)

# Problem 34: dy/dx = sin(y)
y34 = sp.Function('y')(x)
ode34 = sp.Eq(sp.diff(y34,x), sp.sin(y34))
sol34 = sp.dsolve(ode34)

# ========================
# Problem 35: dy/dx = y^2 - x numeric
# ========================
from scipy.integrate import solve_ivp
def f35(t,y): return [y[0]**2 - t]
sol35 = solve_ivp(f35, [0,2], [0], t_eval=np.linspace(0,2,200))

# ========================
# PDEs
# ========================

# Problem 36: Heat equation 1D u_t = u_xx, u(x,0)=sin(pi x), x in [0,1]
nx, nt = 50, 200
dx, dt = 1/nx, 0.0005
alpha = dt/dx**2
u = np.sin(np.pi*np.linspace(0,1,nx+1))
U_heat = [u.copy()]
for n in range(nt):
    u_new = u.copy()
    u_new[1:-1] = u[1:-1] + alpha*(u[2:]-2*u[1:-1]+u[:-2])
    u = u_new
    if n%50==0: U_heat.append(u.copy())

# Problem 37: Wave equation u_tt = u_xx, initial u(x,0)=sin(pi x)
nx, nt = 50, 200
dx, dt = 1/nx, 0.005
c = 1
xgrid = np.linspace(0,1,nx+1)
u_prev = np.sin(np.pi*xgrid)
u = u_prev.copy()
u_next = u.copy()
U_wave = [u_prev.copy()]
for n in range(1,nt):
    if n==1:
        u = u_prev.copy()
    else:
        u_next[1:-1] = 2*u[1:-1]-u_prev[1:-1]+(c*dt/dx)**2*(u[2:]-2*u[1:-1]+u[:-2])
        u_prev, u = u, u_next.copy()
    if n%50==0: U_wave.append(u.copy())

# Problem 38: Laplace equation u_xx+u_yy=0 on grid, Dirichlet boundaries
N = 20
u = np.zeros((N,N))
u[0,:]=1; u[-1,:]=0; u[:,0]=0; u[:,-1]=0
for it in range(2000):
    u[1:-1,1:-1] = 0.25*(u[2:,1:-1]+u[:-2,1:-1]+u[1:-1,2:]+u[1:-1,:-2])
U_laplace = u.copy()

# Problem 39: Diffusion equation u_t = D u_xx, D=0.1
D=0.1
nx, nt = 50, 200
dx, dt = 1/nx, 0.001
alpha = D*dt/dx**2
u = np.exp(-100*(np.linspace(0,1,nx+1)-0.5)**2)
U_diff=[u.copy()]
for n in range(nt):
    u_new=u.copy()
    u_new[1:-1]=u[1:-1]+alpha*(u[2:]-2*u[1:-1]+u[:-2])
    u=u_new
    if n%50==0: U_diff.append(u.copy())

# Problem 40: Wave equation finite difference snapshots
nx=50; dx=1/nx; dt=0.005; nt=200
c=1
xgrid=np.linspace(0,1,nx+1)
u_prev=np.sin(np.pi*xgrid)
u=u_prev.copy()
u_next=u.copy()
snapshots=[]
for n in range(nt):
    if n==0: snapshots.append(u.copy())
    if n==1: continue
    u_next[1:-1]=2*u[1:-1]-u_prev[1:-1]+(c*dt/dx)**2*(u[2:]-2*u[1:-1]+u[:-2])
    u_prev,u=u,u_next.copy()
    if n in [0,int(nt/4),int(nt/2)]: snapshots.append(u.copy())

# ========================
# Plotting
# ========================
plt.figure(figsize=(12,18))

# Problem 35
plt.subplot(3,2,1)
plt.plot(sol35.t, sol35.y[0])
plt.title("Problem 35: dy/dx=y^2-x")

# Problem 36 heat snapshots
plt.subplot(3,2,2)
for k,u in enumerate(U_heat):
    plt.plot(np.linspace(0,1,nx+1),u,label=f"t={k*dt*50:.2f}")
plt.legend(); plt.title("Problem 36: Heat Equation")

# Problem 37 wave snapshots
plt.subplot(3,2,3)
for k,u in enumerate(U_wave):
    plt.plot(xgrid,u,label=f"t={k*dt*50:.2f}")
plt.legend(); plt.title("Problem 37: Wave Equation")

# Problem 38 Laplace
plt.subplot(3,2,4)
plt.imshow(U_laplace,origin="lower",extent=[0,1,0,1])
plt.colorbar(); plt.title("Problem 38: Laplace Equation")

# Problem 39 diffusion
plt.subplot(3,2,5)
for k,u in enumerate(U_diff):
    plt.plot(np.linspace(0,1,nx+1),u,label=f"t={k*dt*50:.2f}")
plt.legend(); plt.title("Problem 39: Diffusion Equation")

# Problem 40 wave snapshots
plt.subplot(3,2,6)
for k,u in enumerate(snapshots):
    plt.plot(xgrid,u,label=f"t={k*dt*50:.2f}")
plt.legend(); plt.title("Problem 40: Wave snapshots")

plt.tight_layout()
plt.show()

# ========================
# Print symbolic results
# ========================
print("Problem 31:", sol31)
print("Problem 32:", sol32)
print("Problem 33:", sol33)
print("Problem 34:", sol34)

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up