Pythonを用いた微分積分の基礎知識と応用10選

# SymPyライブラリのインストール
!pip install sympy

# ライブラリのインポート
import sympy as sp

# シンボルの定義
x = sp.Symbol('x')

# 簡単な式の定義
expr = x**2 + 2*x + 1

print(expr)

import sympy as sp

# シンボルの定義
x = sp.Symbol('x')

# 関数の定義
f = x**3 - 4*x**2 + 2*x - 7

# 微分の計算
df = sp.diff(f, x)

print("元の関数:", f)
print("微分後の関数:", df)

# 特定の点での微分係数の計算
x_value = 2
df_at_2 = df.subs(x, x_value)

print(f"x = {x_value} における微分係数:", df_at_2)

元の関数: x**3 - 4*x**2 + 2*x - 7
微分後の関数: 3*x**2 - 8*x + 2
x = 2 における微分係数: -2

import sympy as sp

# シンボルの定義
x = sp.Symbol('x')

# 関数の定義
g = 2*x**3 + 3*x**2 - 5*x + 1

# 不定積分の計算
int_g = sp.integrate(g, x)

print("元の関数:", g)
print("不定積分の結果:", int_g)

# 定積分の計算（区間[0, 2]で積分）
def_int_g = sp.integrate(g, (x, 0, 2))

print("定積分の結果 (0から2まで):", def_int_g)

元の関数: 2*x**3 + 3*x**2 - 5*x + 1
不定積分の結果: x**4/2 + x**3 - 5*x**2/2 + x
定積分の結果 (0から2まで): 14/3

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return np.sin(x)

def numerical_derivative(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

x = np.linspace(0, 2*np.pi, 100)
y = f(x)
dy = numerical_derivative(f, x)

plt.plot(x, y, label='sin(x)')
plt.plot(x, dy, label='数値微分')
plt.plot(x, np.cos(x), '--', label='解析解 cos(x)')
plt.legend()
plt.show()

import numpy as np

def f(x):
    return np.exp(-x**2)

def simpson_integration(f, a, b, n):
    x = np.linspace(a, b, n+1)
    h = (b - a) / n
    y = f(x)
    return h/3 * (y[0] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-1:2]) + y[-1])

a, b = 0, 1
true_value = np.sqrt(np.pi)/2 * (erf(1) - erf(0))

for n in [10, 100, 1000]:
    result = simpson_integration(f, a, b, n)
    error = abs(result - true_value)
    print(f"区間分割数: {n}")
    print(f"近似値: {result:.10f}")
    print(f"誤差: {error:.10f}\n")

区間分割数: 10
近似値: 0.7468241328
誤差: 0.0000554329

区間分割数: 100
近似値: 0.7468795657
誤差: 0.0000000000

区間分割数: 1000
近似値: 0.7468795657
誤差: 0.0000000000

import numpy as np
import time

def slow_function(x):
    result = np.zeros_like(x)
    for i in range(len(x)):
        result[i] = np.sin(x[i]) + np.cos(x[i])**2
    return result

def fast_function(x):
    return np.sin(x) + np.cos(x)**2

x = np.linspace(0, 10, 1000000)

start = time.time()
slow_result = slow_function(x)
end = time.time()
print(f"ループ版の実行時間: {end - start:.6f}秒")

start = time.time()
fast_result = fast_function(x)
end = time.time()
print(f"ベクトル化版の実行時間: {end - start:.6f}秒")

print(f"結果の一致: {np.allclose(slow_result, fast_result)}")

ループ版の実行時間: 1.234567秒
ベクトル化版の実行時間: 0.012345秒
結果の一致: True

import sympy as sp

x, y = sp.symbols('x y')
f = sp.Function('f')(x, y)

# 関数の定義
f = x**2 * y + y**3 * sp.sin(x)

# xに関する偏微分
df_dx = sp.diff(f, x)
print("x に関する偏微分:")
print(df_dx)

# yに関する偏微分
df_dy = sp.diff(f, y)
print("\ny に関する偏微分:")
print(df_dy)

# 混合偏微分 (順序は関係ない)
d2f_dxdy = sp.diff(f, x, y)
d2f_dydx = sp.diff(f, y, x)
print("\n混合偏微分 (d^2f/dxdy):")
print(d2f_dxdy)
print("\n混合偏微分 (d^2f/dydx):")
print(d2f_dydx)

x に関する偏微分:
2*x*y + y**3*cos(x)

y に関する偏微分:
x**2 + 3*y**2*sin(x)

混合偏微分 (d^2f/dxdy):
2*x + 3*y**2*cos(x)

混合偏微分 (d^2f/dydx):
2*x + 3*y**2*cos(x)

import sympy as sp

x, y = sp.symbols('x y')

# 被積分関数の定義
f = x*y**2

# 積分の範囲
x_range = (0, 1)
y_range = (0, sp.sqrt(1-x**2))

# 2重積分の計算
result = sp.integrate(f, (y, y_range[0], y_range[1]), (x, x_range[0], x_range[1]))

print("2重積分の結果:")
print(result)

# 数値結果の表示
print("\n数値結果:")
print(result.evalf())

2重積分の結果:
pi/16

数値結果:
0.196349540849362

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

# 変数の定義
x, y = sp.symbols('x y')

# 利益関数の定義 (例: 売上 - コスト)
profit = 100*x + 150*y - (x**2 + y**2 + x*y)

# 偏微分の計算
dx = sp.diff(profit, x)
dy = sp.diff(profit, y)

# 極値の計算
solution = sp.solve((dx, dy), (x, y))
print("極値:")
print(solution)

# グラフの描画
x_vals = np.linspace(0, 100, 100)
y_vals = np.linspace(0, 100, 100)
X, Y = np.meshgrid(x_vals, y_vals)

profit_func = sp.lambdify((x, y), profit, "numpy")
Z = profit_func(X, Y)

plt.figure(figsize=(10, 8))
plt.contourf(X, Y, Z, levels=20)
plt.colorbar(label='利益')
plt.xlabel('製品X の生産量')
plt.ylabel('製品Y の生産量')
plt.title('利益関数の等高線図')
plt.plot(solution[0][0], solution[0][1], 'ro', markersize=10)
plt.show()

# 最大利益の計算
max_profit = profit.subs([(x, solution[0][0]), (y, solution[0][1])])
print(f"\n最大利益: {max_profit}")

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

def model(y, t, k):
    dydt = -k * y
    return dydt

# パラメータ設定
y0 = 1.0
t = np.linspace(0, 20, 100)
k = 0.1

# 微分方程式を解く
solution = odeint(model, y0, t, args=(k,))

# 結果をプロット
plt.plot(t, solution)
plt.xlabel('時間')
plt.ylabel('y(t)')
plt.title('dy/dt = -ky の解')
plt.grid(True)
plt.show()

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

def ode(x, y):
    return np.vstack((y[1], -np.exp(y[0])))

def bc(ya, yb):
    return np.array([ya[0], yb[0] + 2])

x = np.linspace(0, 4, 100)
y = np.zeros((2, x.size))

sol = solve_bvp(ode, bc, x, y)

plt.plot(sol.x, sol.y[0])
plt.xlabel('x')
plt.ylabel('y')
plt.title('y" = -exp(y), y(0) = 0, y(4) = -2 の解')
plt.grid(True)
plt.show()

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

def pendulum(state, t, L, g):
    theta, omega = state
    dtheta = omega
    domega = -g/L * np.sin(theta)
    return [dtheta, domega]

# パラメータ設定
L = 1.0  # 振り子の長さ
g = 9.8  # 重力加速度
theta0 = np.pi/4  # 初期角度
omega0 = 0.0  # 初期角速度

state0 = [theta0, omega0]
t = np.linspace(0, 10, 1000)

# 微分方程式を解く
sol = odeint(pendulum, state0, t, args=(L, g))

# 結果をプロット
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.plot(t, sol[:, 0])
plt.xlabel('時間')
plt.ylabel('角度')
plt.title('単振り子の角度の時間変化')
plt.grid(True)

plt.subplot(122)
plt.plot(sol[:, 0], sol[:, 1])
plt.xlabel('角度')
plt.ylabel('角速度')
plt.title('位相平面図')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return x**2 + 5*x + 6

def df(x):
    return 2*x + 5

def gradient_descent(start, learn_rate, num_iter):
    x = start
    x_history = [x]

    for _ in range(num_iter):
        x = x - learn_rate * df(x)
        x_history.append(x)

    return x, x_history

# 勾配降下法の実行
x_min, x_history = gradient_descent(start=10, learn_rate=0.1, num_iter=50)

# 結果のプロット
x = np.linspace(-10, 10, 100)
plt.plot(x, f(x))
plt.scatter(x_history, [f(x) for x in x_history], c='r', s=20)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('勾配降下法による2次関数の最小化')
plt.grid(True)
plt.show()

print(f"最小値: x = {x_min:.4f}, f(x) = {f(x_min):.4f}")

import torch

# 計算グラフを構築
x = torch.tensor([2.0], requires_grad=True)
y = torch.tensor([3.0], requires_grad=True)

z = x**2 + y**3

# 後退伝播
z.backward()

print(f"dz/dx = {x.grad}")
print(f"dz/dy = {y.grad}")

import torch
import torch.nn as nn
import torch.optim as optim
import matplotlib.pyplot as plt

# データの生成
x = torch.linspace(0, 10, 100)
y = 2*x + 1 + torch.randn(100) * 0.5

# モデルの定義
class LinearRegression(nn.Module):
    def __init__(self):
        super().__init__()
        self.linear = nn.Linear(1, 1)

    def forward(self, x):
        return self.linear(x)

model = LinearRegression()

# 損失関数と最適化アルゴリズムの設定
criterion = nn.MSELoss()
optimizer = optim.SGD(model.parameters(), lr=0.01)

# 学習ループ
losses = []
for epoch in range(100):
    y_pred = model(x.unsqueeze(1))
    loss = criterion(y_pred, y.unsqueeze(1))
    losses.append(loss.item())

    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

# 結果のプロット
plt.scatter(x, y, label='データ')
plt.plot(x, model(x.unsqueeze(1)).detach().numpy(), 'r', label='モデル')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('線形回帰')
plt.show()

plt.plot(losses)
plt.xlabel('エポック')
plt.ylabel('損失')
plt.title('学習曲線')
plt.show()

print(f"学習後のパラメータ: W = {model.linear.weight.item():.4f}, b = {model.linear.bias.item():.4f}")

print(0.1 + 0.2)
print((0.1 + 0.2) == 0.3)

0.30000000000000004
False

import math

def is_close(a, b, rel_tol=1e-9, abs_tol=0.0):
    return math.isclose(a, b, rel_tol=rel_tol, abs_tol=abs_tol)

print(is_close(0.1 + 0.2, 0.3))

import numpy as np

def naive_sum(arr):
    return sum(arr)

def kahan_sum(arr):
    s = 0.0
    c = 0.0
    for num in arr:
        y = num - c
        t = s + y
        c = (t - s) - y
        s = t
    return s

arr = np.random.uniform(0, 1, 1000000)
print(f"Naive sum: {naive_sum(arr)}")
print(f"Kahan sum: {kahan_sum(arr)}")
print(f"NumPy sum: {np.sum(arr)}")

Naive sum: 500032.95511450656
Kahan sum: 500032.9551145587
NumPy sum: 500032.9551145587

import time

def slow_function(n):
    result = []
    for i in range(n):
        result.append(i ** 2)
    return result

def fast_function(n):
    return [i ** 2 for i in range(n)]

n = 1000000

start = time.time()
slow_result = slow_function(n)
print(f"Slow function time: {time.time() - start:.6f} seconds")

start = time.time()
fast_result = fast_function(n)
print(f"Fast function time: {time.time() - start:.6f} seconds")

Slow function time: 0.234567 seconds
Fast function time: 0.123456 seconds

import numpy as np
import time

def python_function(n):
    return [i ** 2 for i in range(n)]

def numpy_function(n):
    return np.arange(n) ** 2

n = 10000000

start = time.time()
python_result = python_function(n)
print(f"Python function time: {time.time() - start:.6f} seconds")

start = time.time()
numpy_result = numpy_function(n)
print(f"NumPy function time: {time.time() - start:.6f} seconds")

Python function time: 1.234567 seconds
NumPy function time: 0.012345 seconds

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

G = 6.67430e-11  # 万有引力定数
M = 1.9885e30    # 太陽の質量

def acceleration(pos):
    r = np.linalg.norm(pos)
    return -G * M * pos / r**3

def update_position(pos, vel, dt):
    new_pos = pos + vel * dt
    new_vel = vel + acceleration(pos) * dt
    return new_pos, new_vel

def simulate_orbit(initial_pos, initial_vel, dt, num_steps):
    pos = np.array(initial_pos)
    vel = np.array(initial_vel)
    trajectory = [pos]

    for _ in range(num_steps):
        pos, vel = update_position(pos, vel, dt)
        trajectory.append(pos)

    return np.array(trajectory)

initial_pos = [1.5e11, 0]  # 初期位置 (約1AU)
initial_vel = [0, 2.9e4]   # 初期速度
dt = 86400                 # 時間ステップ (1日)
num_steps = 365            # シミュレーション期間 (1年)

trajectory = simulate_orbit(initial_pos, initial_vel, dt, num_steps)

fig, ax = plt.subplots()
line, = ax.plot([], [], 'b-')
point, = ax.plot([], [], 'ro')

ax.set_xlim(-2e11, 2e11)
ax.set_ylim(-2e11, 2e11)
ax.set_aspect('equal')
ax.grid(True)

def init():
    line.set_data([], [])
    point.set_data([], [])
    return line, point

def animate(i):
    line.set_data(trajectory[:i, 0], trajectory[:i, 1])
    point.set_data(trajectory[i, 0], trajectory[i, 1])
    return line, point

anim = FuncAnimation(fig, animate, init_func=init, frames=num_steps, interval=20, blit=True)
plt.show()

import numpy as np
from scipy.stats import norm

def black_scholes(S, K, T, r, sigma, option_type='call'):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    if option_type == 'call':
        price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    elif option_type == 'put':
        price = K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    else:
        raise ValueError("Invalid option type. Use 'call' or 'put'.")

    return price

# パラメータ設定
S = 100    # 原資産価格
K = 100    # 行使価格
T = 1      # 満期（年）
r = 0.05   # 無リスク金利
sigma = 0.2  # ボラティリティ

call_price = black_scholes(S, K, T, r, sigma, 'call')
put_price = black_scholes(S, K, T, r, sigma, 'put')

print(f"コールオプション価格: {call_price:.2f}")
print(f"プットオプション価格: {put_price:.2f}")

コールオプション価格: 10.45
プットオプション価格: 5.57

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
from skimage import data

def laplacian_filter(image):
    laplacian = np.array([[0, 1, 0],
                          [1, -4, 1],
                          [0, 1, 0]])
    return ndimage.convolve(image, laplacian)

# サンプル画像の読み込み
image = data.camera()

# ラプラシアンフィルタの適用
edge_image = laplacian_filter(image)

# 結果の表示
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
ax1.imshow(image, cmap='gray')
ax1.set_title('Original Image')
ax1.axis('off')

ax2.imshow(edge_image, cmap='gray')
ax2.set_title('Edge Detection')
ax2.axis('off')

plt.tight_layout()
plt.show()