Return to SymPy で理工系の数学C

SymPy で常微分方程式

必要なパッケージを import します。

In [1]:
from sympy.abc import * 
from sympy import *

# SymPy Plotting Backends (SPB): グラフを描く際に利用
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

変数分離法

例題

$$ \frac{dy}{dx} = – 2 x\, y $$

In [2]:
# y を x の関数として宣言
y = Function('y')

# 解くべき微分方程式
eq = Eq(Derivative(y(x), x), -2*x*y(x))
eq
Out[2]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} = – 2 x y{\left(x \right)}$
In [3]:
dsolve(eq, y(x))
Out[3]:
$\displaystyle y{\left(x \right)} = C_{1} e^{- x^{2}}$

上の出力で $C_1$ は積分定数です。

マルサスの人口モデル

$$\frac{dN}{dt} = \gamma \, N$$

を初期条件 $t = t_0$ で $N(t_0) = N_0$ として解く。

In [4]:
N = Function('N')
var('gamma')

eqm = Eq(Derivative(N(t), t), gamma*N(t))
eqm
Out[4]:
$\displaystyle \frac{d}{d t} N{\left(t \right)} = \gamma N{\left(t \right)}$
In [5]:
var('t0, N0')

dsolve(eqm, N(t), ics={N(t).subs(t,t0):N0})
Out[5]:
$\displaystyle N{\left(t \right)} = N_{0} e^{\gamma t} e^{- \gamma t_{0}}$

参考：米国の人口とマルサスモデル

上記によれば，1790年から1930年のアメリカの人口は以下のようになっている（人口の単位は百万人）。

In [6]:
usa = [
# 西暦, 人口（百万人）
  [1790,  3.9], 
  [1800,  5.3],
  [1810,  7.2],
  [1820,  9.6],
  [1830, 12.9],
  [1840, 17.1],
  [1850, 23.2],
  [1860, 31.4],
  [1870, 38.6],
  [1880, 50.2],
  [1890, 62.9],
  [1900, 76.0],
  [1910, 92.0],
  [1920,106.5],
  [1930,123.2]
]

$t_0 = 1790$ （年）とすると，$N_0 = 3.9$ 。Python のリスト要素を取り出すインデックスは 0 ゼロはじまりであるから…

In [7]:
t0 = usa[0][0]
N0 = usa[0][1]

残りのパラメーター $\gamma_1$ は，別の時刻 $t_1$ における $N_1 = N_m(t_1)$ から求める。

\begin{eqnarray}
N_1 &=& N_0 e^{\gamma_1 (t_1 – t_0)} \\
\log \frac{N_1}{N_0} &=& \gamma_1 (t_1 – t_0) \\
\therefore\ \ \gamma_1 &=& \frac{1}{t_1 – t_0} \log \frac{N_1}{N_0}
\end{eqnarray}

たとえば，$t_1 = 1830$ とすると（Python のインデックスは 0 ゼロはじまりだから）…

In [8]:
t1 = usa[4][0]  # 5行目の1列目 
N1 = usa[4][1]  # 5行目の2列目 

gamma1 = 1/(t1-t0) * log(N1/N0)
gamma1
Out[8]:
$\displaystyle 0.0299062689558006$

したがって，マルサスの人口モデルを $N_m(t)$ とすると

In [9]:
def Nm(t):
    return N0 * exp(gamma1*(t-t0))

このようにして求められた $N(t)$ を，人口データ usa と共グラフにしてみます。

In [10]:
# SymPy Plotting Backends (SPB) で陽関数を描く

p1 = plot(Nm(t), (t, 1790, 1940), "マルサスモデル", 
          legend = True, 
          xlabel = "年", ylabel = "人口（単位：百万人）");

SPB で点，$x$ 座標 $y$ 座標の数値データをグラフにする際の書式は

plot_list([x1, ..., xn], [y1, ..., yn])

一方，リスト usa には

usa = [
 [x1, y1],
 [x2, y2],
 ...,
 [xn, yn]
]

のように数値データが格納されているため，$x$ 座標のみ，$y$ 座標のみのリストを作成するために少しだけ工夫が必要です。やり方を2つほど…

  1. for 文を使ってリストを作成する例。
  2. numpy.array() で配列に変換して「列」成分を指定する例。
In [11]:
# SPB で2次元リスト usa を plot_list() する例 1
# for 文で x 成分リスト usat と y 成分リスト usaN を作成
usat = []
usaN = []
for dat in usa:
    usat.append(dat[0])
    usaN.append(dat[1])

p2 = plot_list(usat, usaN, "アメリカの人口", 
               # 線で繋がす点で, 凡例を表示
               is_point = True, legend = True, 
               xlabel = "年", ylabel = "人口（単位：百万人）");

In [12]:
# SPB で2次元リスト usa を plot_list() する例 2
# NumPy 配列に変換して1列目 [:,0] 2列目 [:,1] を指定
import numpy as np
npusa = np.array(usa)

p3 = plot_list(npusa[:,0], npusa[:,1], "アメリカの人口", 
               # 線で繋がす点で, 凡例を表示
               is_point = True, legend = True, 
               xlabel = "年", ylabel = "人口（単位：百万人）");

SPB で上記2つのグラフを重ねて表示させる例。

In [13]:
p4 = p1 + p2
p4.show();

ヴェアフルストによる修正人口モデル

$$
\frac{dN}{dt} = \gamma N \left(1 – \frac{N}{N_{\rm max}}\right)
$$

In [14]:
N = Function('N')
var('gamma Nmax')

eqv = Eq(Derivative(N(t), t), gamma*N(t)*(1-N(t)/Nmax))
eqv
Out[14]:
$\displaystyle \frac{d}{d t} N{\left(t \right)} = \gamma \left(1 – \frac{N{\left(t \right)}}{Nmax}\right) N{\left(t \right)}$
In [15]:
var('t0 N0')

ansv = dsolve(eqv, N(t), ics = {N(t).subs(t,t0):N0})
ansv.simplify()
Out[15]:
$\displaystyle N{\left(t \right)} = \frac{N_{0} Nmax e^{\gamma t}}{N_{0} e^{\gamma t} – \left(N_{0} – Nmax\right) e^{\gamma t_{0}}}$

参考：米国の人口とヴェアフルストモデル

初期条件を $t_0 = 1790$（年）のとき $N(t_0) = N_0 = 3.9$（百万人）とします。

$t_1$（年）と $t_2$（年）の値を使って $\gamma$ と $N_{\rm max}$ を求めます。指数関数を含む連立方程式はなかなか解いてくれないかもしれないので，簡単な代数方程式の形にして解きます。

$N(t)$ の分子分母を $e^{\gamma t}$ で割り，さらに

\begin{eqnarray}
n_0 &\equiv& \frac{N_0}{N_{\rm max}}
\end{eqnarray}

とします。

In [16]:
t0 = usa[0][0]
N0 = usa[0][1]
var('n0 gamma2')

def Nv1(t):
    return N0/(n0 + (1-n0)*exp(gamma2*(t0-t)))
Nv1(t)
Out[16]:
$\displaystyle \frac{3.9}{n_{0} + \left(1 – n_{0}\right) e^{\gamma_{2} \cdot \left(1790 – t\right)}}$

$t_1 = 1850$（年）と $t_2 = 1910$（年）の値を入れて連立方程式の形にします。簡単な連立方程式にするため，さらに $T \equiv e^{-60\, \gamma_2}$ とし，$n_0$ と $T$ に関するシンプルな連立方程式の形にして，solve() で解きます。

In [17]:
t1 = usa[6][0];
eq1 = Eq(usa[6][1], Nv1(t1)) 

var('T')
eq1 = eq1.subs(exp(-60*gamma2), T)
eq1
Out[17]:
$\displaystyle 23.2 = \frac{3.9}{T \left(1 – n_{0}\right) + n_{0}}$
In [18]:
t2 = usa[12][0];
eq2 = Eq(usa[12][1], Nv1(t2))
eq2 = eq2.subs(exp(-60*gamma2), T)
eq2
Out[18]:
$\displaystyle 92.0 = \frac{3.9}{T^{2} \cdot \left(1 – n_{0}\right) + n_{0}}$
In [19]:
ans = solve([eq1, eq2], [n0, T])
ans
Out[19]:
[(0.0200125277046207, 0.151115116017121)]

$n_0, \ T$ からもとの $N_{\rm max}, \ \gamma_2$ の値になおすと…

In [20]:
n0, T = ans[0]
Nmax = N0/n0
gamma2 = -1/60*log(T)

def NV(t):
    return N0*Nmax*exp(gamma2*t)/(
        N0*exp(gamma2*t)-(N0-Nmax)*exp(gamma2*t0))
In [21]:
p5 = plot(NV(t), (t, 1790, 1940), "ヴェアフルストモデル", 
          legend = True, 
          xlabel = "年", ylabel = "人口（単位：百万人）");

In [22]:
p6 = p1+p2+p5
p6.show();

1階線形微分方程式と積分因子法

SymPy では，1階線形微分方程式は，特に積分因子法を coding しなくても，dsolve() 関数で解くことができます。

例題

$$ \frac{dy}{dx} + \frac{y}{x} = \frac{\sin x}{x} $$

SymPy での三角関数 $\sin x, \cos x, \tan x$ は sin(x), cos(x), tan(x) と書きます。

In [23]:
y = Function('y')
var('x')
eq = Eq(Derivative(y(x), x)+y(x)/x, sin(x)/x)
eq
Out[23]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} + \frac{y{\left(x \right)}}{x} = \frac{\sin{\left(x \right)}}{x}$
In [24]:
dsolve(eq, y(x))
Out[24]:
$\displaystyle y{\left(x \right)} = \frac{C_{1} – \cos{\left(x \right)}}{x}$

参考：あえて積分因子法の公式で…

この例題をあえて積分因子法の公式に従ってやってみる。

一般的な1階線形微分方程式（非同次方程式の格好しているもの）を

$$\frac{dy}{dx} + P(x) y = Q(x)$$

と書くと，積分因子 $g(x)$ は

$$g(x) = \exp \left\{\int^x P(x’)\,dx’ \right\}$$

解は

$$y = \frac{1}{g(x)} \left\{\int^x g(x’) Q(x’) \, dx’ + C \right\}$$

この例題では

$$P(x) = \frac{1}{x}, \ Q(x) = \frac{\sin x}{x}$$

In [25]:
def P(x):
    return 1/x

def Q(x):
    return sin(x)/x

g = exp(integrate(P(x), x))
g
Out[25]:
$\displaystyle x$
In [26]:
y = 1/g * (integrate(g * Q(x), x) + C)
y
Out[26]:
$\displaystyle \frac{C – \cos{\left(x \right)}}{x}$

最も簡単な定数係数2階微分方程式

最も簡単な定数係数2階線形微分方程式 $y” + y = 0$ と $y” – y = 0$ を解く。

In [27]:
y = Function('y')
var('x')
eq1 = Derivative(y(x), x, 2) + y(x)
eq1
Out[27]:
$\displaystyle y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)}$
In [28]:
dsolve(eq1, y(x))
Out[28]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(x \right)} + C_{2} \cos{\left(x \right)}$
In [29]:
eq2 = Derivative(y(x), x, 2) - y(x)
eq2
Out[29]:
$\displaystyle – y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)}$
In [30]:
dsolve(eq2, y(x))
Out[30]:
$\displaystyle y{\left(x \right)} = C_{1} e^{- x} + C_{2} e^{x}$

参考：脊髄反射によらずに解く

$y” + y = 0$ を脊髄反射によらずに解く。両辺に $2 y’$ をかけて整理すると

$$\left(\left(y’\right)^2 + y^2\right)’ = 0$$

微分してゼロということは，かっこの中身は定数であるということなので，

$$ \left(y’\right)^2 + y^2 = \mbox{const.} \equiv a^2$$

とおける。つまり

$$\frac{dy}{dx} = \pm \sqrt{a^2 – y^2}$$

を解けばよい。まずは $+$ の式から…

In [31]:
# 一般性を失うことなく第1の積分定数の平方根 a を正にできる
var('a', positive = True)

eq1 = Eq(Derivative(y(x), x), sqrt(a**2 - y(x)**2))
eq1
Out[31]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} = \sqrt{a^{2} – y^{2}{\left(x \right)}}$
In [32]:
dsolve(eq1, y(x))
Out[32]:
$\displaystyle y{\left(x \right)} = a \sin{\left(C_{1} + x \right)}$
In [33]:
eq2 = Eq(Derivative(y(x), x), -sqrt(a**2 - y(x)**2))
eq2
Out[33]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} = – \sqrt{a^{2} – y^{2}{\left(x \right)}}$
In [34]:
dsolve(eq2, y(x))
Out[34]:
$\displaystyle y{\left(x \right)} = a \sin{\left(C_{1} – x \right)}$

$y(x) = a \sin\left(C_1 – x\right) = -a \sin(x – C_1)$
ということで 2 つの解をまとめると，積分定数あらためて $\pm a \Rightarrow C_2, \pm C_1 \Rightarrow C_1$ とおいて

\begin{eqnarray}
y &=& C_2 \sin(x + C_1) \\
&=& C_2 \left(\sin x \cos C_1 + \cos x \sin C_1 \right) \\
&=& (C_2 \cos C_1) \sin x + (C_2 \sin C_1) \cos x \\
&\equiv& A \cos x + B \sin x
\end{eqnarray}

のように書けるということがわかる。

最も簡単な定数係数2階微分方程式：続き

$y” + K y = 0$ あるいは移項して $y” = – K y$ を解く。

$K > 0$ の場合

In [35]:
y = Function('y')
var('x')
var('K', positive = True)

eq1 = Eq(Derivative(y(x), x, 2), -K*y(x))
eq1
Out[35]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = – K y{\left(x \right)}$
In [36]:
dsolve(eq1, y(x))
Out[36]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(\sqrt{K} x \right)} + C_{2} \cos{\left(\sqrt{K} x \right)}$

$K < 0$ の場合

$K \equiv – K_p, \ K_p > 0$ として…

In [37]:
var('K', negative = True)

eq2 = Eq(Derivative(y(x), x, 2), -K*y(x))
eq2
Out[37]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = – K y{\left(x \right)}$
In [38]:
dsolve(eq2, y(x))
Out[38]:
$\displaystyle y{\left(x \right)} = C_{1} e^{- x \sqrt{- K}} + C_{2} e^{x \sqrt{- K}}$

$K = 0$ の場合

In [39]:
eq3 = Eq(Derivative(y(x), x, 2), 0)
eq3
Out[39]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = 0$
In [40]:
dsolve(eq3, y(x))
Out[40]:
$\displaystyle y{\left(x \right)} = C_{1} + C_{2} x$

人類の至宝：オイラーの公式

SymPy がオイラーの公式を知っているか確認する。

SymPy では虚数単位 $i$ は I です。以下のように import しています。

In [41]:
# π，ネイピア数，虚数単位
from sympy import pi, E, I

指数関数 $e^{i x}$ を .rewrite(cos) をつけて三角関数で書き直してみます。

In [42]:
Eq(exp(I*x), 
   exp(I*x).rewrite(cos))
Out[42]:
$\displaystyle e^{i x} = i \sin{\left(x \right)} + \cos{\left(x \right)}$

オイラーの等式

$\theta = \pi$ のとき，$e^{i \pi} + 1 = 0$。$\pi$ は SymP では pi です。

In [43]:
exp(I*pi) + 1
Out[43]:
$\displaystyle 0$

オイラーの公式からみた三角関数と双曲線関数の関係

三角関数 $\cos x, \sin x$ を exponentialize() で（複素）指数関数表示します。

In [44]:
cos(x).rewrite(exp)
Out[44]:
$\displaystyle \frac{e^{i x}}{2} + \frac{e^{- i x}}{2}$
In [45]:
sin(x).rewrite(exp)
Out[45]:
$\displaystyle – \frac{i \left(e^{i x} – e^{- i x}\right)}{2}$

双曲線関数 $\cosh x, \sinh x$ も指数関数表示します。

In [46]:
cosh(x).rewrite(exp)
Out[46]:
$\displaystyle \frac{e^{x}}{2} + \frac{e^{- x}}{2}$
In [47]:
sinh(x).rewrite(exp)
Out[47]:
$\displaystyle \frac{e^{x}}{2} – \frac{e^{- x}}{2}$

さて，三角関数や双曲線関数の変数（引数）が虚数でもいいのだと拡張すると…

In [48]:
cosh(I*theta)
Out[48]:
$\displaystyle \cos{\left(\theta \right)}$
In [49]:
sinh(I*theta)
Out[49]:
$\displaystyle i \sin{\left(\theta \right)}$
In [50]:
cos(I*x)
Out[50]:
$\displaystyle \cosh{\left(x \right)}$
In [51]:
sin(I*x)
Out[51]:
$\displaystyle i \sinh{\left(x \right)}$

例題

微分方程式 $y” + K y = 0$ の $K > 0$ の場合の解は

\begin{eqnarray}
y &=& A \cos\left( \sqrt{K} x\right) + \frac{B}{\sqrt{K}} \sin\left( \sqrt{K} x\right)
\end{eqnarray}

のように書ける。このことを使って，$K<0$ の場合と $K = 0$ の場合の解を上記の解から直接求める。

In [52]:
var('A B x')
var('K', positve =True)

y1 = A*cos(sqrt(K)*x) + B/sqrt(K)*sin(sqrt(K)*x)
y1
Out[52]:
$\displaystyle A \cos{\left(\sqrt{K} x \right)} + \frac{B \sin{\left(\sqrt{K} x \right)}}{\sqrt{K}}$
$K < 0$ の場合

$K = – |K|$ とすると…

In [53]:
y2 = y1.subs(K, -abs(K))
y2
Out[53]:
$\displaystyle A \cosh{\left(x \sqrt{\left|{K}\right|} \right)} + \frac{B \sinh{\left(x \sqrt{\left|{K}\right|} \right)}}{\sqrt{\left|{K}\right|}}$
$K = 0$ の場合
In [54]:
Eq(Limit(y1, K, 0), 
   Limit(y1, K, 0).doit())
Out[54]:
$\displaystyle \lim_{K \to 0^+}\left(A \cos{\left(\sqrt{K} x \right)} + \frac{B \sin{\left(\sqrt{K} x \right)}}{\sqrt{K}}\right) = A + B x$

定数係数2階線形同次方程式

定数係数2階線形微分方程式（同次方程式）は以下のように書ける。

$$ \frac{d^2 y}{dx^2} + 2 b \frac{dy}{dx} + cy = 0$$

$b, c$ は定数。一般解 $y$ は以下のように場合分けして…

In [55]:
y = Function('y')
var('x')
var('b c')

eq = Derivative(y(x), x, 2) + 2*b*Derivative(y(x), x) + c*y(x)
eq
Out[55]:
$\displaystyle 2 b \frac{d}{d x} y{\left(x \right)} + c y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)}$

特に仮定なしで eq = 0 を dsolve() したときの解。ルートの中身 $b^2 -c >0$ を仮定したときの解になっている。

In [56]:
ans1 = dsolve(eq, y(x))
ans1.simplify()
Out[56]:
$\displaystyle y{\left(x \right)} = \left(C_{1} e^{x \sqrt{b^{2} – c}} + C_{2} e^{- x \sqrt{b^{2} – c}}\right) e^{- b x}$

ひとつだけの解ですますには…

以下の表現ひとつで3つの場合の全てに対応した解になることを示しておきます。

In [57]:
def y(mu, x):
    return exp(-b*x)*(A*cos(mu*x) + B/mu * sin(mu*x))
y(mu, x)
Out[57]:
$\displaystyle \left(A \cos{\left(\mu x \right)} + \frac{B \sin{\left(\mu x \right)}}{\mu}\right) e^{- b x}$

$c-b^2 > 0$ の場合は $\mu \equiv \sqrt{c-b^2}$ として…

In [58]:
y(sqrt(c - b**2), x)
Out[58]:
$\displaystyle \left(A \cos{\left(x \sqrt{- b^{2} + c} \right)} + \frac{B \sin{\left(x \sqrt{- b^{2} + c} \right)}}{\sqrt{- b^{2} + c}}\right) e^{- b x}$

$c-b^2 < 0$ の場合は $\mu \equiv \sqrt{c-b^2} = i \sqrt{b^2-c}$ として…

In [59]:
y(I*sqrt(b**2-c), x)
Out[59]:
$\displaystyle \left(A \cosh{\left(x \sqrt{b^{2} – c} \right)} + \frac{B \sinh{\left(x \sqrt{b^{2} – c} \right)}}{\sqrt{b^{2} – c}}\right) e^{- b x}$

最後に $c-b^2 = 0$ の場合は $\mu \rightarrow 0$ の極限をとって…

In [60]:
limit(y(mu, x), mu, 0)
Out[60]:
$\displaystyle \left(A + B x\right) e^{- b x}$

定数係数2階線形非同次方程式

人力で解く際には，同次方程式と非同次方程式とでは，解く手間がずいぶん違ったが，SymPy ではどちらも同じ。dsolve() を使う。

例題

非同次方程式 $y” + a^2 y = \sin b x$ を解く。

In [61]:
y = Function('y')
var('x')
var('a b', positive = True)

eq = Eq(Derivative(y(x), x, 2) + a**2 * y(x), sin(b*x))
eq
Out[61]:
$\displaystyle a^{2} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = \sin{\left(b x \right)}$
In [62]:
dsolve(eq, y(x))
Out[62]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(a x \right)} + C_{2} \cos{\left(a x \right)} + \frac{\sin{\left(b x \right)}}{a^{2} – b^{2}}$

$a = b$ の場合は上記の解の分母がゼロになってしまうので，別途計算する必要がある。

この状況は，力学では固有振動数 $a$
に等しい振動数の外力が加えられた時に起こる「共鳴（共振）」と呼ばれる現象である。

In [63]:
b = a
eq2 = Eq(Derivative(y(x), x, 2) + a**2 * y(x), sin(b*x))
eq2
Out[63]:
$\displaystyle a^{2} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = \sin{\left(a x \right)}$
In [64]:
dsolve(eq2, y(x))
Out[64]:
$\displaystyle y{\left(x \right)} = C_{2} \sin{\left(a x \right)} + \left(C_{1} – \frac{x}{2 a}\right) \cos{\left(a x \right)}$

積分定数 $C_1, C_2 $がつくのは同次方程式の一般解。非同次方程式の特解
$\displaystyle -\frac{x}{2a} \cos (a x)$ は，振幅が $x$ に比例して単調増加していく。

あえてロンスキアンを使って特殊解を求める

非同次方程式も dsolve() で解がもとまったわけだが，そこをあえてロンスキアンを使って特殊解を求めてみる。

まず，同次方程式は…

In [65]:
eq0 = Eq(Derivative(y(x), x, 2) + a**2 * y(x), 0)
eq0
Out[65]:
$\displaystyle a^{2} y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} = 0$
In [66]:
dsolve(eq0, y(x))
Out[66]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(a x \right)} + C_{2} \cos{\left(a x \right)}$

上記のように同次方程式の1次独立な基本解はそれぞれ…

In [67]:
def y1(x):
    return sin(a*x)
def y2(x):
    return cos(a*x)

ロンスキアンは $W(x) \equiv y_1 y_2′ – y_1′ y_2 = \cdots$

In [68]:
W = (y1(x)*diff(y2(x), x) - diff(y1(x), x)*y2(x)).simplify()
W
Out[68]:
$\displaystyle – a$

非同次方程式の特殊解 $y_s(x)$ は，非同次項を $R(x) = \sin b x$ として

$$y_s(x) = y_2(x) \int^x \frac{R(x’) y_1 (x’)}{W(x’)} dx’ – y_1(x) \int^x \frac{R(x’) y_2(x’)}{W(x’)} dx’$$

In [69]:
var('a b', positive = True)
def R(x):
    return sin(b*x)

ys = ( y2(x) * integrate(R(x)*y1(x)/W, x) 
      -y1(x) * integrate(R(x)*y2(x)/W, x))

ys.simplify().expand()
Out[69]:
$\displaystyle \begin{cases} \frac{\sin{\left(b x \right)}}{a^{2} – b^{2}} & \text{for}\: a \neq b \\- \frac{x \cos{\left(a x \right)}}{2 a} + \frac{\sin{\left(b x \right)} \cos{\left(a x – b x \right)}}{2 a b} & \text{otherwise} \end{cases}$

上記の otherwise の場合の答えが今ひとつです。直接 $a = b$ とすれば…

In [70]:
def R(x):
    return sin(a*x)

ys = ( y2(x) * integrate(R(x)*y1(x)/W, x) 
      -y1(x) * integrate(R(x)*y2(x)/W, x))

ys.simplify().expand()
Out[70]:
$\displaystyle – \frac{x \cos{\left(a x \right)}}{2 a} + \frac{\sin{\left(a x \right)}}{2 a^{2}}$

otherwise, すなわち $a = b$ の場合の非同次方程式の特殊解を求めたはずなのに，$\sin (a x)$ に比例する項も出てきました。これは，同次方程式の一般解に吸収できる項です。

非同次方程式の特殊解 $y_s$ を求める公式には，積分の下限の任意性があるので，試しに以下のように

$$y_s(x) = y_2(x) \int_0^x \frac{R(t) y_1 (t)}{W(t)} dt – y_1(x) \int_0^x \frac{R(t) y_2(t)}{W(t)} dt$$

とすると，$a \neq b$ の場合の特殊解は…

In [71]:
def R(x):
    return sin(b*x)

ys = ( y2(x) * integrate(R(t)*y1(t)/W, (t, 0, x)) 
      -y1(x) * integrate(R(t)*y2(t)/W, (t, 0, x)))

ys.simplify().expand()
Out[71]:
$\displaystyle \begin{cases} \frac{a \sin{\left(b x \right)}}{a^{3} – a b^{2}} – \frac{b \sin{\left(a x \right)}}{a^{3} – a b^{2}} & \text{for}\: a \neq b \\- \frac{x \cos{\left(a x \right)}}{2 a} + \frac{\sin{\left(b x \right)} \cos{\left(a x – b x \right)}}{2 a b} & \text{otherwise} \end{cases}$

… となり，$ \displaystyle y_s = \frac{\sin (b x)}{a^2-b^2}$ のほかに，$\sin (a x)$ に比例する項も現れます。しかし，この項は同次方程式の一般解に吸収できる項です。

このように，ロンスキアンを使った公式で特殊解を求めると，不定積分の不定性（積分の下限の不定性）により，同次方程式の一般解に組み込まれる項があらわれることもありますので，そのへんは各自処理していただくということで…