SymPy で（あえて）数値解析

# SymPy を使うときのおまじない
from sympy import *
from sympy.abc import *
# 虚数単位，円周率，ネイピア数
from sympy import I, pi, E
# グラフは SymPy Plotting Backends (SPB) で描く
from spb import *

# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']

integrate(sin(sin(x)), (x, 0, pi))

integrate(sin(sin(x)), (x, 0, pi)).evalf()

plot(besselj(0,x), 0, (x, 0, 10), 
     legend = False, # 凡例は無しに
     xlabel = "$x$", ylabel = "$J_0(x)$", 
     xlim = [0, 10], 
     title = "ベッセル関数 $J_0(x)$");

# besselj(0,x) = 0 となる解を
# x が 2.5 のあたりから始めて 15 桁の精度で求める

nsolve(besselj(0,x), x, 2.5, prec=15)

nsolve(besselj(0,x), x, [2, 4], prec=15)

# ほんとに解になってるか確認。ゼロに近ければよい。
besselj(0,_).evalf()

def rk1(func, x0, t0, t1, N):
    """ dx/dt = func(t, x) を
        初期条件 t = t0 で x = x0 とし，
        [t0, t1] を N 分割して
        4次のルンゲ・クッタ法で解く
    """
    h = (t1 - t0)/N
    T = [t0]
    X = [x0]
    t = t0
    x = x0
    for i in range(N):
        k1 = h*func(t, x)
        k2 = h*func(t + h/2, x + k1/2)
        k3 = h*func(t + h/2, x + k2/2)
        k4 = h*func(t + h, x + k3)
        x = x + (k1 + 2*k2 + 2*k3 + k4)/6
        t = t + h
        T.append(t)
        X.append(x)
    return [T, X]

# dx/dt = f(t, x) を解く。
def f(t, x):
    return (1-x)*x

x0 = 0.1
t0 = 0
t1 = 10

for N in [100, 200, 400, 800, 1000]:
    tmp = rk1(f, x0, t0, t1, N)
    print("N =%5d, x(t1) =%15.12f" % (N, tmp[1][-1]))

# dx/dt = f(t, x) を解く。
def f(t, x):
    return (1-x)*x

x0 = 0.1
t0 = 0
t1 = 10

ndiv = 40
Ndiv = 20
N = ndiv * Ndiv

ans1 = rk1(f, x0, t0, t1, N)

# データを間引いてグラフに
T1, X1 = ans1
t1 = [T1[i] for i in range(0, len(T1), ndiv)]
x1 = [X1[i] for i in range(0, len(X1), ndiv)]

plot_list(t1, x1, is_point = True, 
          xlabel = "$t$", ylabel = "$x(t)$",
          title = "ロジスティック方程式の数値解", 
          # "ms" つまり "markersize" を小さめに
          rendering_kw={"ms":3});

# 解くべき微分方程式を eq として定義
var('t')
x = Function('x')(t)
eq = Eq(Derivative(x,t), (1-x)*x)
eq

# 初期条件を特に設定せずに微分方程式を解く
dsolve(eq, x)

# 初期条件を t=0 で x(0) = 1/10 として解く
sol1 = dsolve(eq, x, ics={x.subs(t,0): Rational(1,10)})
sol1

# ans の右辺を表示 
# rhs は rirhgt-hand-side
simplify(sol1.rhs)

p1 = plot(sol1.rhs, (t, 0, 10), "解析解", 
          xlabel = "$t$", ylabel = "$x(t)$", 
          title = "ロジスティック方程式の解", 
          rendering_kw={"color":"red"},
          show=False)

p2 = plot_list(t1, x1, "数値解", 
               is_point = True, 
               # "ms" つまり "markersize" を小さめに
               rendering_kw={"ms":3, "color":"tab:blue"}, 
               show=False)

(p1 + p2).show();

def xsol1(s):
    return sol1.rhs.subs(t,s)

xsol1(t)

def rk2(Func, x0, v0, t0, t1, N):
    """ d^2 x/dt^2 = Func(t, x, v) を
        初期条件 t = t0 で x = x0, v = v0 とし，
        [t0, t1] を N 分割して
        4次のルンゲ・クッタ法で解く
    """
    h = (t1 - t0)/N
    T = [t0]
    X = [x0]
    V = [v0]
    t = t0
    x = x0
    v = v0
    for i in range(N):
        k1 = h*v
        m1 = h*Func(t, x, v)
        k2 = h*(v+m1/2)
        m2 = h*Func(t+h/2, x+k1/2, v+m1/2)
        k3 = h*(v+m2/2)
        m3 = h*Func(t+h/2, x+k2/2, v+m2/2)
        k4 = h*(v+m3)
        m4 = h*Func(t+h, x+k3, v+m3)
        x = x + (k1 + 2*k2 + 2*k3 + k4)/6
        v = v + (m1 + 2*m2 + 2*m3 + m4)/6
        t = t + h 
        T.append(t)
        X.append(x)
        V.append(v)
    return [T, X, V]

def F(t, x, v):
    global a, b
    return -x - a*v + b*cos(t)

a = 0.5
b = 0.2
x0 = 3
v0 = 0
t0 = 0
t1 = 20

for N in [200, 400, 800, 1000]:
    T, X, V = rk2(F, x0, v0, t0, t1, N)
    print("N =%5d, x(t1) =%15.12f" % (N, X[-1]))

ndiv = 20
Ndiv = 40
N = ndiv * Ndiv

T2, X2, V2 = rk2(F, x0, v0, t0, t1, N)

t2 = [T2[i] for i in range(0, len(T2), ndiv)]
x2 = [X2[i] for i in range(0, len(X2), ndiv)]

plot_list(t2, x2, 
          is_point = True, 
          xlabel = "$t$", ylabel = "$x(t)$",
          title = "減衰＋強制振動", 
          # "ms" つまり "markersize" を小さめに
          rendering_kw={"ms":3});

t = Symbol('t')
x = Function('x')(t)
a, b = symbols('a b')

def eq(a, b):
    return Eq(Derivative(x, t, 2), 
              -x -a*Derivative(x, t) + b*cos(t))

eq(a, b)

sol0 = dsolve(eq(0, 0), x, 
         ics = {x.subs(t, 0):3, diff(x, t).subs(t,0):0})
sol0

a = Rational(5,10)
b = Rational(2,10)

sol2 = dsolve(eq(a, b), x, 
         ics = {x.subs(t, 0):3, diff(x, t).subs(t,0):0})
sol2

def xsol2(s):
    return sol2.rhs.subs(t, s)

key = "$a = %3.1f, b = %3.1f$ の解析解" % (a, b)
p4 = plot(xsol2(t), (t, 0, 20), key, 
          xlim = (0, 20), 
          xlabel = "$t$", ylabel = "$x(t)$",
          title = "減衰＋強制振動", show=False);

key = "$a = %3.1f, b = %3.1f$ の数値解" % (a, b)
p5 = plot_list(t2, x2, key,
          is_point = True, 
          title = "減衰＋強制振動", 
          # "ms" つまり "markersize" を小さめに
          rendering_kw={"ms":3}, show=False)

(p4+p5).show();