Pythonで理解する名門の森


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/de859efd-50bc-bff7-49f5-301bb88fab8a.png)
解説の方は動画でみるとよいかもしれません。ここではPythonコードでシミュレーションプログラムを作っていきます。コピペして使ってください。



https://studyphys.com/lm-of-ua/

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/f642be07-52e8-5363-60ec-210d8a5d2999.png)
```
import numpy as np
import matplotlib.pyplot as plt

# 関数定義
def x(t, initial_velocity, acceleration):
    return initial_velocity * t + 0.5 * acceleration * (t**2)

# 時間の範囲を設定
t_values = np.linspace(0, 10, 100)  # 0から10の範囲を100分割

# 初速度と加速度を指定
initial_velocity = 5
acceleration = 9.8

# x(t)の計算
position = x(t_values, initial_velocity, acceleration)

# x(t)のプロット
plt.figure(figsize=(10, 5))
plt.plot(t_values, position, label='x(t)')
plt.xlabel('Time (t)')
plt.ylabel('Position (x)')
plt.title('Position vs Time')
plt.legend()
plt.grid(True)
plt.show()

# x(t)をtで微分した結果をプロット
velocity = initial_velocity + acceleration * t_values
plt.figure(figsize=(10, 5))
plt.plot(t_values, velocity, label='Velocity')
plt.xlabel('Time (t)')
plt.ylabel('Velocity')
plt.title('Velocity vs Time')
plt.legend()
plt.grid(True)
plt.show()

# x(t)をtで2回微分した結果をプロット
acceleration_values = np.full_like(t_values, acceleration)
plt.figure(figsize=(10, 5))
plt.plot(t_values, acceleration_values, label='Acceleration')
plt.xlabel('Time (t)')
plt.ylabel('Acceleration')
plt.title('Acceleration vs Time')
plt.legend()
plt.grid(True)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/b0ef22a6-35dd-705e-3714-e037423d2415.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/c9b22eba-4aba-18e8-9d7c-b5a92319fc54.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/7d8cefb7-74c4-775a-fd01-d4d36bd9435d.png)




微分方程式の解説
ラプラス変換と伝達関数を使って微分方程式の解を求める

https://www.youtube.com/watch?v=hT2CGxDOy8c&t=918s




https://controlabo.com/laplace-transform/

https://www.momoyama-usagi.com/entry/math-seigyo06

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/a2b11432-c1da-fa02-9d4b-7acd14678810.png)
F(s)をラプラス変換で求めて、X(s)を逆ラプラス変換してX(t)を求める
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/dc4eb23e-b56f-0fed-ee22-ba731969a151.png)





![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/52a41075-aa4d-142c-f1cf-46f562f5ff08.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/587dd0b8-2f09-b456-d0be-fe831df48493.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/4a562f51-6938-0026-0f83-8bac8cb7decf.png)

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

# Parameter settings
A = 1  # Value of A
omega = 2  # Value of ω
C = 3  # Value of constant C
K = 1  # Value of K

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

# x(t)
x = A * np.cos(omega * t) + C / K

# x'(t)
x_prime = -A * omega * np.sin(omega * t)

# x''(t)
x_double_prime = -A * omega**2 * np.cos(omega * t)

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

plt.plot(t, x, label='x(t)', color='blue')
plt.plot(t, x_prime, label="x'(t)", color='red')
plt.plot(t, x_double_prime, label="x''(t)", color='green')

plt.title('Plot of x(t), x\'(t), x\'\'(t)')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.grid(True)
plt.show()


```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/2cadcb34-247c-e16d-c5ad-c48f84a1732e.png)


https://w3e.kanazawa-it.ac.jp/math/physics/category/mechanics/masspoint_mechanics/damped_harmonic_motion/henkan-tex.cgi?target=/math/physics/category/mechanics/masspoint_mechanics/damped_harmonic_motion/damped_harmonic_motion.html&pcview=2

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

# 初期条件
A = 50
d = 2
omega = 1  # 適切な値に置き換える必要があります

# データポイントの準備
t = np.array([0, 2*np.pi/omega, 4*np.pi/omega, 6*np.pi/omega])
x = np.array([A, A-4*d, A-8*d, A-12*d])

# y=d および y=-d をプロットするためのデータ
y_d = np.array([d, d, d, d])
y_minus_d = np.array([-d, -d, -d, -d])

# プロット
plt.plot(t, x, marker='o', label='Data Points')
plt.plot(t, x, linestyle='-', color='blue')  # データポイントを結ぶ線
plt.plot(t, y_d, linestyle='--', color='red', label='y=d')
plt.plot(t, y_minus_d, linestyle='--', color='green', label='y=-d')

# ラベルと凡例の追加
plt.xlabel('Time')
plt.ylabel('Position')
plt.title('Position vs Time')
plt.legend()

# グリッドを表示
plt.grid(True)

# プロットを表示
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/c6b29cf4-e26a-069d-2fe3-7b24ec9d4388.png)

https://controlabo.com/first-order-system-impulse-step-response/

https://domeblog.net/first-order-system/



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/b675b824-359a-2b20-fb0a-910d09d90e07.png)






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

# パラメータの値を設定
F = 1.0  # 力の値
D = 1.0  # 減衰の値
T = 1.0  # 時定数

# 微分方程式の右辺を定義
def diff_eq1(t):
    return (F/D) * (1 - np.exp(-t / T))

def diff_eq2(t):
    return (F/D) * np.exp(-t / T)

# 時間の範囲を設定
t = np.linspace(0, 10, 100)  # 0から10までの時間を100分割

# 微分方程式の解析解を計算
analytical_solution1 = (F/D) * (t - T * (1 - np.exp(-t / T)))
analytical_solution2 = (F/D) * (T * (np.exp(-t / T) - 1))

# 解析解と数値解のプロット
plt.plot(t, diff_eq1(t), label="X'(t) = (F/D)(1 - e^(-t/T))")
plt.plot(t, diff_eq2(t), label="X'(t) = (F/D)(e^(-t/T))")
plt.xlabel('Time')
plt.ylabel('X\'(t)')
plt.title('Analytical Solutions of the Differential Equations')
plt.legend()
plt.grid(True)
plt.show()

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/9c5ae2f2-3f55-c2a5-1b3f-ce003d511908.png)

3番

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

# 定数
e = 0.5

# θの値を生成
theta = np.linspace(0, 2*np.pi, 100)

# 関数を定義
def y1(theta):
    return np.exp(1) * np.cos(theta)

def y2(theta):
    return np.exp(2) * np.cos(theta)

def y3(theta):
    return np.exp(3) * np.cos(theta)

# プロット
plt.figure(figsize=(8, 6))
plt.plot(theta, y1(theta), label='y1 = e^1 * v(cosθ)')
plt.plot(theta, y2(theta), label='y2 = e^2 * v(cosθ)')
plt.plot(theta, y3(theta), label='y3 = e^3 * v(cosθ)')
plt.xlabel('θ')
plt.ylabel('y')
plt.title('Functions Plot')
plt.legend()
plt.grid(True)
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/aeef72c9-ec47-ef2b-d08a-b4b9dd57db5d.png)

18番

```
import matplotlib.pyplot as plt

# パラメータ
m = 1
M = 2
l = 3
θ = 30  # 角度を度数法で指定

# 点の座標計算
x1, y1 = 0, 0
x2, y2 = m * l * (1 - abs(θ)/90) / (m + M), 0
x3, y3 = l, 0

# プロット
plt.plot([x1, x2, x3], [y1, y2, y3], 'ro')  # 赤い点でプロット
plt.text(x1, y1, '(0, 0)', fontsize=12, ha='right')
plt.text(x2, y2, f'({x2}, 0)', fontsize=12, ha='right')
plt.text(x3, y3, f'({x3}, 0)', fontsize=12, ha='right')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Plot of points in two dimensions')
plt.grid(True)
plt.axis('equal')  # アスペクト比を等しくする
plt.show()

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/fc9f7d8a-85bc-0231-3599-09204d0829b8.png)

電磁気編

```
from sympy import symbols, sqrt, solve

# Define the symbols
x, λ, m = symbols('x λ m')

# Define the equation
equation = sqrt((3 * λ)**2 + x**2) - x - (m + 1/2) * λ

# Solve for x
solution = solve(equation, x)
print(solution)
```


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

def F(m):
    return (9 - (m + 0.5)**2) / (2 * m + 1)

m_values = np.linspace(-5, 5, 400)  # Get 400 points from -5 to 5

plt.figure(figsize=(8, 6))  # Set the size of the plot
plt.plot(m_values, F(m_values), label='$F(m)$')  # Plot the values of F(m)
plt.xlabel('$m$')  # x-axis label
plt.ylabel('$F(m)$')  # y-axis label
plt.title('Graph of $F(m)$')  # Title of the graph
plt.grid(True)  # Show grid
plt.legend()  # Show legend
plt.show()  # Display the graph


```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/00832a78-69a4-e139-de9d-e724f70d610e.png)





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

# Function definition
def W(x, epsilon, C, V, I):
    return -(epsilon - 1) * C * V**2 * x / (2 * (I + (epsilon - 1) * x))

# Parameters
epsilon = 1.5  # Relative permittivity
C = 1         # Capacitance (in farads)
V = 10        # Voltage (in volts)
I = 5         # Current (in amperes)

# Generate x values
x_values = np.linspace(0, 10, 100)

# Calculate corresponding W values
W_values = W(x_values, epsilon, C, V, I)

# Plot
plt.plot(x_values, W_values)
plt.xlabel('x')
plt.ylabel('W(x)')
plt.title('Plot of W(x)')
plt.grid(True)
plt.show()


```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/aed6719a-d94c-5833-3ba9-9e39444d8483.png)




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

# Define the function
def Q(x, a, CV):
    return -((2 * x**2) / (a**2 - x**2)) * CV

# Define the range of x values
x_values = np.linspace(-5, 5, 400)  # Adjust the range as needed

# Choose a value for 'a' and 'CV'
a_value = 3
CV_value = 1

# Calculate corresponding y values
y_values = Q(x_values, a_value, CV_value)

# Plot the function
plt.plot(x_values, y_values)
plt.title('Plot of Q(x)')
plt.xlabel('x')
plt.ylabel('Q(x)')
plt.grid(True)
plt.show()



```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/6c8f156d-ad54-de78-dae6-4a9c53266778.png)





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

# Define the function
def x(n, V):
    return 2*V - V / (2**(n-1))

# Define the range of n values
n_values = np.arange(1, 10)  # Adjust the range as needed

# Choose a value for 'V'
V_value = 10

# Calculate corresponding x values
x_values = x(n_values, V_value)

# Plot the function
plt.plot(n_values, x_values, marker='o', linestyle='-')
plt.title('Plot of x(n)')
plt.xlabel('n')
plt.ylabel('x(n)')
plt.grid(True)
plt.show()



```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/e9479286-4bdd-222f-7f13-1168bc6bfe58.png)


3番
https://www.ne.jp/asahi/tokyo/nkgw/www_2/gakusyu/gensi/kousi/kansyou_kaisetu_kaisetu/kansyou_kaisetui_kaisetu_2.html

```


import numpy as np
import matplotlib.pyplot as plt

# 関数の定義
def f(theta, a, wavelength):
    numerator = np.sin((np.pi * a * np.sin(theta)) / wavelength)
    denominator = (np.pi * a * np.sin(theta)) / wavelength
    return (numerator / denominator) ** 2

# パラメータの設定
a = 1.0  # 例として a を1に設定
wavelength = 0.5  # 例として波長を0.5に設定
theta = np.linspace(-np.pi/2, np.pi/2, 1000)  # -π/2 から π/2 までの角度を1000個生成

# プロット
plt.plot(theta, f(theta, a, wavelength))
plt.xlabel('θ')
plt.ylabel('f(θ)')
plt.title('Plot of f(θ)')
plt.grid(True)
plt.show()


```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3760770/f5ab890c-404e-5ffc-daf7-bddfc842f927.png)