1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@mrrclb48z

「微分方程式の右辺を眺める：dx/dt = 0, 1, x, x^2, x^3 の時間発展」様を、 sympyでAI先生に教えてもらいました。

1
Last updated at Posted at 2026-05-25

・やり直しました。numpy　を使わないで下さい。が必要でした。
・やり直します。numpy　を使わないで下さい。が必要でした。申し訳ありません。

オリジナル

@sakaimaging(Yusuke Sakai)@ 様へ。いつもお世話になっております。

(やり直し)ChatGPT先生へ。

・同じ初期値で時間発展を比べる

以下を、sympyのソースコードに書き換えて下さい。
(ソースコードを省略)

import sympy as sp
import matplotlib.pyplot as plt


# 記号
t = sp.symbols('t', real=True)
x0 = sp.symbols('x0', real=True)


# ==========================================================
# 各微分方程式の解
# ==========================================================

def sol_zero(t, x0):
    return x0


def sol_one(t, x0):
    return x0 + t


def sol_x(t, x0):
    return x0 * sp.exp(t)


def sol_x2(t, x0):
    """
    dx/dt = x^2

    解:
        x(t) = x0 / (1 - x0*t)

    x0 > 0 の場合、
        t = 1/x0
    で発散する。
    """
    return x0 / (1 - x0*t)


def sol_x3(t, x0):
    """
    dx/dt = x^3

    解:
        x(t) = x0 / sqrt(1 - 2*x0^2*t)

    t = 1/(2*x0^2)
    で発散する。
    """
    return x0 / sp.sqrt(1 - 2*x0**2 * t)


# ==========================================================
# モデル一覧
# ==========================================================

models = {
    r"$dx/dt=0$": sol_zero,
    r"$dx/dt=1$": sol_one,
    r"$dx/dt=x$": sol_x,
    r"$dx/dt=x^2$": sol_x2,
    r"$dx/dt=x^3$": sol_x3,
}


# ==========================================================
# 数値代入
# ==========================================================

x0_val = 0.2

# t の値を SymPy で生成
t_values = [5*i/999 for i in range(1000)]


# ==========================================================
# 描画
# ==========================================================

fig, ax = plt.subplots(figsize=(8, 5))

for label, func in models.items():

    expr = func(t, x0)

    y_values = []

    for tv in t_values:

        try:
            y = expr.subs({
                t: tv,
                x0: x0_val
            }).evalf()

            y = float(y)

            # 発散対策
            if abs(y) > 10:
                y = float('nan')

        except Exception:
            y = float('nan')

        y_values.append(y)

    ax.plot(t_values, y_values, label=label)


# x軸
ax.axhline(0, color="black", linewidth=0.8)

# ラベル
ax.set_xlabel("t")
ax.set_ylabel(r"$x(t)$")

# タイトル
ax.set_title(fr"Time evolution from $x_0={x0_val}$")

# 表示範囲
ax.set_ylim(-0.5, 8)

# グリッド
ax.grid(alpha=0.3)

# 凡例
ax.legend()

plt.tight_layout()
plt.show()

0Figure_1.png

(やり直し続けて)ChatGPT先生へ。

初期値を変えてみる

sympy のソースコードに書き換えて下さい。
(ソースコードを省略)

import sympy as sp
import matplotlib.pyplot as plt


# ==========================================================
# 記号
# ==========================================================

t = sp.symbols('t', real=True)
x0 = sp.symbols('x0', real=True)


# ==========================================================
# 微分方程式の解
# ==========================================================

def sol_zero(t, x0):
    return x0


def sol_one(t, x0):
    return x0 + t


def sol_x(t, x0):
    return x0 * sp.exp(t)


def sol_x2(t, x0):
    """
    dx/dt = x^2

    解:
        x(t) = x0 / (1 - x0*t)

    x0 > 0 の場合、
        t = 1/x0
    で発散
    """
    return x0 / (1 - x0*t)


def sol_x3(t, x0):
    """
    dx/dt = x^3

    解:
        x(t) = x0 / sqrt(1 - 2*x0^2*t)

    t = 1/(2*x0^2)
    で発散
    """
    return x0 / sp.sqrt(1 - 2*x0**2 * t)


# ==========================================================
# モデル
# ==========================================================

models = {
    r"$dx/dt=0$": sol_zero,
    r"$dx/dt=1$": sol_one,
    r"$dx/dt=x$": sol_x,
    r"$dx/dt=x^2$": sol_x2,
    r"$dx/dt=x^3$": sol_x3,
}


# ==========================================================
# 初期値
# ==========================================================

x0_list = [-1.0, -0.5, -0.2, 0.2, 0.5, 1.0]

# NumPy を使わず t を生成
t_values = [5 * i / 999 for i in range(1000)]


# ==========================================================
# 描画準備
# ==========================================================

fig, axes = plt.subplots(
    2,
    3,
    figsize=(15, 8),
    sharex=True,
    sharey=True
)

axes = axes.ravel()


# ==========================================================
# 描画
# ==========================================================

for ax, x0_val in zip(axes, x0_list):

    for label, func in models.items():

        expr = func(t, x0)

        y_values = []

        for tv in t_values:

            try:

                y = expr.subs({
                    t: tv,
                    x0: x0_val
                }).evalf()

                y = float(y)

                # 発散や複素数を除外
                if abs(y) > 10:
                    y = float('nan')

            except Exception:

                y = float('nan')

            y_values.append(y)

        ax.plot(
            t_values,
            y_values,
            label=label
        )

    # 基準線
    ax.axhline(
        0,
        color="black",
        linewidth=0.8
    )

    ax.axhline(
        1,
        color="gray",
        linewidth=0.8,
        linestyle="--"
    )

    ax.axhline(
        -1,
        color="gray",
        linewidth=0.8,
        linestyle="--"
    )

    # タイトル
    ax.set_title(
        fr"$x_0={x0_val}$"
    )

    # 軸ラベル
    ax.set_xlabel("t")
    ax.set_ylabel(r"$x(t)$")

    # 範囲
    ax.set_ylim(-6, 6)

    # グリッド
    ax.grid(alpha=0.3)


# 凡例
axes[-1].legend(
    loc="upper right",
    fontsize=9
)

plt.tight_layout()
plt.show()

02Figure_1.png

以下、旧です。ーーーーーーーーーーーーーーーーーーーーーーー

以下、旧です。ーーーーーーーーーーーーーーーーーーーーーーー

以下、旧です。ーーーーーーーーーーーーーーーーーーーーーーー

以下、旧です。ーーーーーーーーーーーーーーーーーーーーーーー

(本日の)ChatGPT先生へ。

・同じ初期値で時間発展を比べる

以下を、sympyのソースコードに書き換えて下さい。
(ソースコードを省略)

同じ初期値で時間発展を比べる

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


# -----------------------------
# 記号
# -----------------------------
t = sp.symbols('t', real=True)
x0 = sp.symbols('x0', real=True)


# -----------------------------
# 解
# -----------------------------
sol_zero = x0
sol_one  = x0 + t
sol_x    = x0 * sp.exp(t)
sol_x2   = x0 / (1 - x0*t)
sol_x3   = x0 / sp.sqrt(1 - 2*x0**2*t)


# -----------------------------
# 初期値
# -----------------------------
x0_val = 0.2

sols = {
    r"$dx/dt=0$": sol_zero.subs(x0, x0_val),
    r"$dx/dt=1$": sol_one.subs(x0, x0_val),
    r"$dx/dt=x$": sol_x.subs(x0, x0_val),
    r"$dx/dt=x^2$": sol_x2.subs(x0, x0_val),
    r"$dx/dt=x^3$": sol_x3.subs(x0, x0_val),
}


# -----------------------------
# lambdify
# -----------------------------
funcs = {
    label: sp.lambdify(t, expr, "numpy")
    for label, expr in sols.items()
}


# -----------------------------
# 時間
# -----------------------------
t_vals = np.linspace(0, 5, 1000)


# -----------------------------
# 描画
# -----------------------------
fig, ax = plt.subplots(figsize=(8, 5))

for label, func in funcs.items():

    with np.errstate(divide='ignore', invalid='ignore'):
        y = func(t_vals)

    # -------------------------
    # スカラー対策
    # -------------------------
    if np.isscalar(y):
        y = np.full_like(t_vals, y, dtype=float)

    y = np.asarray(y, dtype=float)

    # 発散除去
    y = np.where(np.isfinite(y), y, np.nan)
    y = np.where(np.abs(y) < 10, y, np.nan)

    ax.plot(t_vals, y, label=label)


ax.axhline(0, color='black', linewidth=0.8)

ax.set_xlabel("t")
ax.set_ylabel(r"$x(t)$")

ax.set_title(fr"Time evolution from $x_0={x0_val}$")

ax.set_ylim(-0.5, 8)

ax.grid(alpha=0.3)
ax.legend()

plt.tight_layout()
plt.show()

2*sqrt

0Figure_1.png

(続けて)ChatGPT先生へ。

sympy のソースコードに書き換えて下さい。
(ソースコードを省略)

初期値を変えてみる

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


# ==========================================
# 記号定義
# ==========================================
t, x0 = sp.symbols('t x0', real=True)


# ==========================================
# 微分方程式の厳密解
# ==========================================

# dx/dt = 0
sol_zero = x0

# dx/dt = 1
sol_one = x0 + t

# dx/dt = x
sol_x = x0 * sp.exp(t)

# dx/dt = x^2
sol_x2 = x0 / (1 - x0*t)

# dx/dt = x^3
sol_x3 = x0 / sp.sqrt(1 - 2*x0**2*t)


# ==========================================
# 数式表示
# ==========================================
print("\n=== Exact Solutions ===\n")

print("dx/dt = 0")
sp.pprint(sol_zero)

print("\ndx/dt = 1")
sp.pprint(sol_one)

print("\ndx/dt = x")
sp.pprint(sol_x)

print("\ndx/dt = x^2")
sp.pprint(sol_x2)

print("\ndx/dt = x^3")
sp.pprint(sol_x3)


# ==========================================
# lambdify
# ==========================================
f_zero = sp.lambdify((t, x0), sol_zero, "numpy")
f_one  = sp.lambdify((t, x0), sol_one,  "numpy")
f_x    = sp.lambdify((t, x0), sol_x,    "numpy")
f_x2   = sp.lambdify((t, x0), sol_x2,   "numpy")
f_x3   = sp.lambdify((t, x0), sol_x3,   "numpy")


models = {
    r"$dx/dt=0$": f_zero,
    r"$dx/dt=1$": f_one,
    r"$dx/dt=x$": f_x,
    r"$dx/dt=x^2$": f_x2,
    r"$dx/dt=x^3$": f_x3,
}


# ==========================================
# 初期値
# ==========================================
x0_list = [-1.0, -0.5, -0.2, 0.2, 0.5, 1.0]

# 時間
t_vals = np.linspace(0, 5, 1000)


# ==========================================
# 描画
# ==========================================
fig, axes = plt.subplots(
    2,
    3,
    figsize=(15, 8),
    sharex=True,
    sharey=True
)

axes = axes.ravel()


for ax, x0_val in zip(axes, x0_list):

    for label, func in models.items():

        with np.errstate(divide='ignore', invalid='ignore'):

            y = func(t_vals, x0_val)

        # ----------------------------------
        # スカラー対策
        # ----------------------------------
        if np.isscalar(y):
            y = np.full_like(t_vals, y, dtype=float)

        y = np.asarray(y, dtype=float)

        # ----------------------------------
        # 発散除去
        # ----------------------------------
        y = np.where(np.isfinite(y), y, np.nan)
        y = np.where(np.abs(y) < 10, y, np.nan)

        # ----------------------------------
        # x^2, x^3 の発散点除去
        # ----------------------------------

        # dx/dt = x^2
        if label == r"$dx/dt=x^2$":

            denom = 1 - x0_val * t_vals

            if x0_val > 0:
                y = np.where(denom > 0, y, np.nan)

        # dx/dt = x^3
        if label == r"$dx/dt=x^3$":

            denom = 1 - 2 * x0_val**2 * t_vals

            y = np.where(denom > 0, y, np.nan)

        ax.plot(t_vals, y, label=label)

    # --------------------------------------
    # 補助線
    # --------------------------------------
    ax.axhline(0, color="black", linewidth=0.8)

    ax.axhline(
        1,
        color="gray",
        linewidth=0.8,
        linestyle="--"
    )

    ax.axhline(
        -1,
        color="gray",
        linewidth=0.8,
        linestyle="--"
    )

    # --------------------------------------
    # 軸設定
    # --------------------------------------
    ax.set_title(fr"$x_0={x0_val}$")

    ax.set_xlabel("t")

    ax.set_ylabel(r"$x(t)$")

    ax.set_ylim(-6, 6)

    ax.grid(alpha=0.3)


# 凡例
axes[-1].legend(loc="upper right", fontsize=9)

plt.tight_layout()

plt.show()

=== Exact Solutions ===

dx/dt = 0
x₀

dx/dt = 1
t + x₀

dx/dt = x
    t
x₀⋅ℯ 

dx/dt = x^2
    x₀   
─────────
-t⋅x₀ + 1

dx/dt = x^3
        x₀        
──────────────────
   _______________
  ╱         2     
╲╱  - 2⋅t⋅x₀  + 1 
2*sqrt(5)

02Figure_1.png

いつもと違うおすすめです。

wikipedia で

sympy doc

・

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?