Return to Maxima で理工系の数学C

Maxima で常微分方程式

変数分離法

例題

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

解くべき方程式を，変数 eq に代入します。diff() はMaxima で微分を行う関数ですが，微分方程式を定義するときは ' で始まるのが大事です。

'diff(y, x) ではなく，diff(y, x) と ' を最初につけるのを忘れるとどうなるかというと…

In [1]:
diff(y, x);
Out[1]:
\[\tag{${\it \%o}_{1}$}0\]
In [2]:
/* 解くべき方程式を，変数 eq に代入。*/

eq: 'diff(y, x) = -2*x*y;
Out[2]:
\[\tag{${\it \%o}_{2}$}\frac{d}{d\,x}\,y=-2\,x\,y\]
In [3]:
/* 微分方程式を解く関数は ode2() 。*/
ode2(eq, y, x);
Out[3]:
\[\tag{${\it \%o}_{3}$}y={\it \%c}\,e^ {- x^2 }\]

上記の答えの %c は積分定数 $C$ を表します。

マルサスの人口モデル

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

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

In [4]:
eq2: 'diff(N, t) = gamma * N;
Out[4]:
\[\tag{${\it \%o}_{4}$}\frac{d}{d\,t}\,N=N\,\gamma\]
In [5]:
ode2(eq2, N, t);
Out[5]:
\[\tag{${\it \%o}_{5}$}N={\it \%c}\,e^{t\,\gamma}\]
In [6]:
/* 1階微分方程式の初期条件 */

ic1(%, t = t0, N = N0);
Out[6]:
\[\tag{${\it \%o}_{6}$}N=N_{0}\,e^{t\,\gamma-t_{0}\,\gamma}\]

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

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

In [7]:
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$ 。したがって，マルサスの人口モデルを $N_m(t)$ とすると

In [8]:
t0: usa[1][1]; /* 1行目の1列目 */
N0: usa[1][2]; /* 1行目の2列目 */

kill(t, gamma)$

Nm(t) := N0 * exp(gamma1*(t-t0));
Out[8]:
\[\tag{${\it \%o}_{8}$}1790\]
Out[8]:
\[\tag{${\it \%o}_{9}$}3.9\]
Out[8]:
\[\tag{${\it \%o}_{11}$}{\it Nm}\left(t\right):=N_{0}\,\exp \left(\gamma_{1}\,\left(t-t_{0}\right)\right)\]

残りのパラメーター $\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$ とすると…

In [9]:
t1: usa[5][1]; /* 5行目の1列目 */
N1: usa[5][2]; /* 5行目の2列目 */

gamma1: 1/(t1-t0) * log(N1/N0);
Out[9]:
\[\tag{${\it \%o}_{12}$}1830\]
Out[9]:
\[\tag{${\it \%o}_{13}$}12.9\]
Out[9]:
\[\tag{${\it \%o}_{14}$}0.02990626895580065\]

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

Maxima でグラフを作成する関数として plot2d() などのいわゆる plot 系と，draw2d() などのいわゆる draw 系があります。パラメータ設定等の書式が異なりますので，それぞれの例をあげておきます。

In [10]:
/* plot2d() によるグラフ */

plot2d([[discrete, usa], Nm(t)], [t, 1790, 1940], 
       [style, [points, 1.5], lines], 
       [legend, "アメリカの人口", "マルサスモデル"], 
       [xlabel, "年"], [ylabel, "人口（単位：百万人）"])$

In [11]:
/* draw2d() によるグラフ */
draw2d(
  color = "red",
  key = "マルサスモデル", 
  explicit(Nm(t), t, 1790, 1940), 
  
  color = "blue", 
  point_type = 7,
  point_size = 0.7,
  key = "アメリカの人口", 
  points(usa), 
  
  xlabel = "年", ylabel = "人口（単位：百万人）",
  yrange = [-20, 410], xaxis = true
)$

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

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

In [12]:
/* 念のため，変数を初期化 */
kill(N, t, Nmax, N0, t0)$

eq3: 'diff(N, t) = gamma2 * N * (1 - N/Nmax);
Out[12]:
\[\tag{${\it \%o}_{20}$}\frac{d}{d\,t}\,N=N\,\left(1-\frac{N}{{\it Nmax}}\right)\,\gamma_{2}\]
In [13]:
ans: ode2(eq3, N, t);
Out[13]:
\[\tag{${\it \%o}_{21}$}-\frac{\log \left(N-{\it Nmax}\right)-\log N}{\gamma_{2}}=t+{\it \%c}\]

どうやら，Maxima はちゃんと $N= \dots$ のように解いてはくれないようです。あとは人間の手で整理していくことになります。

上式 ans の両辺に $\gamma_2$ をかけた式を ans2 とし，

In [14]:
ans2: gamma2 * ans;
Out[14]:
\[\tag{${\it \%o}_{22}$}\log N-\log \left(N-{\it Nmax}\right)=\gamma_{2}\,\left(t+{\it \%c}\right)\]

両辺の exponential をとった式を ans3 とする。lhs() は左辺（left-hand-side），rhs() は右辺（right-hand-side）をとる関数。

In [15]:
ans3: exp(lhs(ans2)) = exp(rhs(ans2));
Out[15]:
\[\tag{${\it \%o}_{23}$}\frac{N}{N-{\it Nmax}}=e^{\gamma_{2}\,\left(t+{\it \%c}\right)}\]

この ans3 式を N について解く。方程式を解くのは solve() 関数。

In [16]:
sol: solve(ans3, N);
Out[16]:
\[\tag{${\it \%o}_{24}$}\left[ N=\frac{{\it Nmax}\,e^{\gamma_{2}\,t+{\it \%c}\,\gamma_{2}}}{e^{\gamma_{2}\,t+{\it \%c}\,\gamma_{2}}-1} \right] \]

最後に初期条件，$t = t_0$ で $N = N_0$ を入れる。（一般に solve() の解は複数かも知れないのでリスト形式 [... ] で出力される。1個目をとるのが sol[1]。）

In [17]:
sol1: ic1(sol[1], t = t0, N = N0);
Out[17]:
\[\tag{${\it \%o}_{25}$}N=\frac{N_{0}\,{\it Nmax}\,e^{\gamma_{2}\,t}}{\left({\it Nmax}-N_{0}\right)\,e^{\gamma_{2}\,t_{0}}+N_{0}\,e^{\gamma_{2}\,t}}\]
In [18]:
define(Nv(t), rhs(sol1));
Out[18]:
\[\tag{${\it \%o}_{26}$}{\it Nv}\left(t\right):=\frac{N_{0}\,{\it Nmax}\,e^{\gamma_{2}\,t}}{\left({\it Nmax}-N_{0}\right)\,e^{\gamma_{2}\,t_{0}}+N_{0}\,e^{\gamma_{2}\,t}}\]

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

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

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

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

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

とします。

In [19]:
t0: usa[1][1]; /* 1行目の1列目 */
N0: usa[1][2]; /* 1行目の2列目 */

/* n0 = N0/Nmax とおいた */

Nv1(t):= N0/(n0 + (1-n0)*exp(gamma2*(t0-t)));
Out[19]:
\[\tag{${\it \%o}_{27}$}1790\]
Out[19]:
\[\tag{${\it \%o}_{28}$}3.9\]
Out[19]:
\[\tag{${\it \%o}_{29}$}{\it Nv}_{1}\left(t\right):=\frac{N_{0}}{n_{0}+\left(1-n_{0}\right)\,\exp \left(\gamma_{2}\,\left(t_{0}-t\right)\right)}\]

$t_1 = 1850$（年）と $t_2 = 1910$（年）の値を入れて連立方程式の形にします。

In [20]:
t1: usa[7][1];
t2: usa[13][1];

eq1: usa[7][2] = Nv1(t1);
eq2: usa[13][2] = Nv1(t2);
Out[20]:
\[\tag{${\it \%o}_{30}$}1850\]
Out[20]:
\[\tag{${\it \%o}_{31}$}1910\]
Out[20]:
\[\tag{${\it \%o}_{32}$}23.2=\frac{3.9}{n_{0}+e^ {- 60\,\gamma_{2} }\,\left(1-n_{0}\right)}\]
Out[20]:
\[\tag{${\it \%o}_{33}$}92.0=\frac{3.9}{n_{0}+e^ {- 120\,\gamma_{2} }\,\left(1-n_{0}\right)}\]

簡単な連立方程式にするため，さらに $T \equiv e^{-60\, \gamma_2}$ とし，$n_0$ と $T$ に関するシンプルな連立方程式の形にして，solve() で解きます。

In [21]:
kill(T, n0)$

eq1: subst(exp(-60*gamma2)=T, eq1);
eq2: subst(exp(-120*gamma2)=T**2, eq2);

ans: solve([eq1, eq2], [n0, T]);

rat: replaced 23.2 by 116/5 = 23.2

rat: replaced -3.9 by -39/10 = -3.9

rat: replaced 92.0 by 92/1 = 92.0

rat: replaced -3.9 by -39/10 = -3.9
Out[21]:
\[\tag{${\it \%o}_{35}$}23.2=\frac{3.9}{n_{0}+T\,\left(1-n_{0}\right)}\]
Out[21]:
\[\tag{${\it \%o}_{36}$}92.0=\frac{3.9}{n_{0}+T^2\,\left(1-n_{0}\right)}\]
Out[21]:
\[\tag{${\it \%o}_{37}$}\left[ \left[ n_{0}=\frac{87477}{4371112} , T=\frac{3354}{22195} \right] \right] \]

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

In [22]:
Nmax: float(N0/ev(n0, ans[1][1]));
gamma2: float(-1/60*log(ev(T, ans[1][2])));
Out[22]:
\[\tag{${\it \%o}_{38}$}194.8779313419527\]
Out[22]:
\[\tag{${\it \%o}_{39}$}0.03149522291474969\]
In [23]:
/* plot2d() によるグラフ */

plot2d([Nm(t), [discrete, usa], Nv(t)], [t, 1790, 1940], 
       [point_type,bullet],
       [style, lines, [points, 1.5], lines], 
       [color, red, blue, green],
       [legend, "マルサスモデル", "アメリカの人口", "ヴェアフルストモデル"], 
       [xlabel, "年"], [ylabel, "人口（単位：百万人）"])$

In [24]:
/* draw2d() によるグラフ */
draw2d(
  color = "red",
  key = "マルサスモデル", 
  explicit(Nm(t), t, 1790, 1940), 
  
  color = "blue", 
  point_type = 7,
  point_size = 0.7,
  key = "アメリカの人口", 
  points(usa), 
  
  color = "green",
  key = "ヴェアフルストモデル",
  explicit(Nv(t), t, 1790, 1940), 
  
  xlabel = "年", ylabel = "人口（単位：百万人）",
  yrange = [-20, 410], xaxis = true
)$

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

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

例題

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

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

In [25]:
eq: 'diff(y, x) + y/x = sin(x)/x;
Out[25]:
\[\tag{${\it \%o}_{44}$}\frac{d}{d\,x}\,y+\frac{y}{x}=\frac{\sin x}{x}\]
In [26]:
ode2(eq, y, x);
Out[26]:
\[\tag{${\it \%o}_{45}$}y=\frac{{\it \%c}-\cos x}{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 [27]:
P(x):= 1/x;
Q(x):= sin(x)/x;

define(g(x), exp(integrate(P(x), x)));

y = 1/g(x) * (integrate(g(x)*Q(x), x) + C);
Out[27]:
\[\tag{${\it \%o}_{46}$}P\left(x\right):=\frac{1}{x}\]
Out[27]:
\[\tag{${\it \%o}_{47}$}Q\left(x\right):=\frac{\sin x}{x}\]
Out[27]:
\[\tag{${\it \%o}_{48}$}g\left(x\right):=x\]
Out[27]:
\[\tag{${\it \%o}_{49}$}y=\frac{C-\cos x}{x}\]

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

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

In [28]:
eq1: 'diff(y, x, 2) + y = 0;
Out[28]:
\[\tag{${\it \%o}_{50}$}\frac{d^2}{d\,x^2}\,y+y=0\]
In [29]:
ode2(eq1, y, x);
Out[29]:
\[\tag{${\it \%o}_{51}$}y={\it \%k}_{1}\,\sin x+{\it \%k}_{2}\,\cos x\]
In [30]:
eq2: 'diff(y, x, 2) - y = 0;
Out[30]:
\[\tag{${\it \%o}_{52}$}\frac{d^2}{d\,x^2}\,y-y=0\]
In [31]:
ode2(eq2, y, x);
Out[31]:
\[\tag{${\it \%o}_{53}$}y={\it \%k}_{1}\,e^{x}+{\it \%k}_{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 [32]:
assume(a > 0)$
eq1: 'diff(y, x) = sqrt(a**2 - y**2);
Out[32]:
\[\tag{${\it \%o}_{55}$}\frac{d}{d\,x}\,y=\sqrt{a^2-y^2}\]
In [33]:
ode2(eq1, y, x);
Out[33]:
\[\tag{${\it \%o}_{56}$}\arcsin \left(\frac{y}{a}\right)=x+{\it \%c}\]

$y = $ の形に解きます。

In [34]:
solve(%, y);
Out[34]:
\[\tag{${\it \%o}_{57}$}\left[ y=a\,\sin \left(x+{\it \%c}\right) \right] \]

$-$ の式についても同様に…

In [35]:
eq2: 'diff(y, x) = -sqrt(a**2 - y**2);
ode2(eq2, y, x);
solve(%, y);
Out[35]:
\[\tag{${\it \%o}_{58}$}\frac{d}{d\,x}\,y=-\sqrt{a^2-y^2}\]
Out[35]:
\[\tag{${\it \%o}_{59}$}-\arcsin \left(\frac{y}{a}\right)=x+{\it \%c}\]
Out[35]:
\[\tag{${\it \%o}_{60}$}\left[ y=-a\,\sin \left(x+{\it \%c}\right) \right] \]

ということで 2 つの解をまとめると，積分定数 $C_1, C_2$ を使って

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

のように書けるということ。

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

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

In [36]:
eq: 'diff(y, x, 2) = - K * y;
Out[36]:
\[\tag{${\it \%o}_{61}$}\frac{d^2}{d\,x^2}\,y=-K\,y\]

$K > 0$ の場合

In [37]:
ode2(eq, y, x);

/* positive; と回答する */

XIs K positive, negative or zero?
\Is K positive, negative or zero?
positive;
Out[37]:
\[\tag{${\it \%o}_{62}$}y={\it \%k}_{1}\,\sin \left(\sqrt{K}\,x\right)+{\it \%k}_{2}\,\cos \left(\sqrt{K}\,x\right)\]

$K < 0$ の場合

In [38]:
ans3: ode2(eq, y, x);

/* negative; と回答する */

XIs K positive, negative or zero?
\Is K positive, negative or zero?
negative;
Out[38]:
\[\tag{${\it \%o}_{63}$}y={\it \%k}_{1}\,e^{i\,\sqrt{K}\,x}+{\it \%k}_{2}\,e^ {- i\,\sqrt{K}\,x }\]

上記の解は間違いではないかも知れないが，ややこしい。

$K \equiv -|K|$ としてみる。

In [39]:
ev(%, K=-abs(K));
Out[39]:
\[\tag{${\it \%o}_{64}$}y={\it \%k}_{2}\,e^{\sqrt{\left| K\right| }\,x}+{\it \%k}_{1}\,e^ {- \sqrt{\left| K\right| }\,x }\]

$K = 0$ の場合

In [40]:
ode2(eq, y, x);

/* zero; と回答する */

XIs K positive, negative or zero?
\Is K positive, negative or zero?
zero;
Out[40]:
\[\tag{${\it \%o}_{65}$}y={\it \%k}_{2}\,x+{\it \%k}_{1}\]

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

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

Maxima では虚数単位 $i$ は %i です。

In [41]:
exp(%i * theta);
Out[41]:
\[\tag{${\it \%o}_{66}$}e^{i\,\vartheta}\]

demoivre() 関数は複素指数関数を三角関数で表します。

In [42]:
exp(%i * theta) = demoivre(exp(%i * theta));
Out[42]:
\[\tag{${\it \%o}_{67}$}e^{i\,\vartheta}=i\,\sin \vartheta+\cos \vartheta\]

オイラーの等式

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

In [43]:
exp(%i * %pi) + 1;
Out[43]:
\[\tag{${\it \%o}_{68}$}0\]

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

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

In [44]:
cos(x) = exponentialize(cos(x));
sin(x) = exponentialize(sin(x));
Out[44]:
\[\tag{${\it \%o}_{69}$}\cos x=\frac{e^{i\,x}+e^ {- i\,x }}{2}\]
Out[44]:
\[\tag{${\it \%o}_{70}$}\sin x=-\frac{i\,\left(e^{i\,x}-e^ {- i\,x }\right)}{2}\]

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

In [45]:
cosh(x) = exponentialize(cosh(x));
sinh(x) = exponentialize(sinh(x));
Out[45]:
\[\tag{${\it \%o}_{71}$}\cosh x=\frac{e^{x}+e^ {- x }}{2}\]
Out[45]:
\[\tag{${\it \%o}_{72}$}\sinh x=\frac{e^{x}-e^ {- x }}{2}\]

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

In [46]:
cosh(%i * theta);
sinh(%i * theta);
Out[46]:
\[\tag{${\it \%o}_{73}$}\cos \vartheta\]
Out[46]:
\[\tag{${\it \%o}_{74}$}i\,\sin \vartheta\]
In [47]:
cos(%i * x);
sin(%i * x);
Out[47]:
\[\tag{${\it \%o}_{75}$}\cosh x\]
Out[47]:
\[\tag{${\it \%o}_{76}$}i\,\sinh x\]

例題

微分方程式 $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$ の場合の解を上記の解から直接求める。

$K > 0$ の場合
In [48]:
y = y1: A * cos(sqrt(K) * x) + B/sqrt(K) * sin(sqrt(K) * x);
Out[48]:
\[\tag{${\it \%o}_{77}$}y=\frac{B\,\sin \left(\sqrt{K}\,x\right)}{\sqrt{K}}+A\,\cos \left(\sqrt{K}\,x\right)\]
$K < 0$ の場合
In [49]:
y = ev(y1, K=-abs(K));
Out[49]:
\[\tag{${\it \%o}_{78}$}y=\frac{B\,\sinh \left(\sqrt{\left| K\right| }\,x\right)}{\sqrt{\left| K\right| }}+A\,\cosh \left(\sqrt{\left| K\right| }\,x\right)\]
$K = 0$ の場合
In [50]:
y = limit(y1, K, 0);
Out[50]:
\[\tag{${\it \%o}_{79}$}y=B\,x+A\]

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

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

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

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

In [51]:
eq: 'diff(y, x, 2) + 2 * b * 'diff(y, x) + c * y = 0;
Out[51]:
\[\tag{${\it \%o}_{80}$}\frac{d^2}{d\,x^2}\,y+2\,b\,\left(\frac{d}{d\,x}\,y\right)+c\,y=0\]

$c – b^2 > 0$ の場合

In [52]:
ode2(eq, y, x), factor;

/* c - b**2 は positive; と答える */

X        2
Is c - b  positive, negative or zero?
\        2
Is c - b  positive, negative or zero?
positive;
Out[52]:
\[\tag{${\it \%o}_{81}$}y=e^ {- b\,x }\,\left({\it \%k}_{1}\,\sin \left(\sqrt{c-b^2}\,x\right)+{\it \%k}_{2}\,\cos \left(\sqrt{c-b^2}\,x\right)\right)\]

$c – b^2 < 0$ の場合

In [53]:
ode2(eq, y, x), factor;

/* c - b**2 は negative; と答える */

X        2
Is c - b  positive, negative or zero?
\        2
Is c - b  positive, negative or zero?
negative;
Out[53]:
\[\tag{${\it \%o}_{82}$}y=\left({\it \%k}_{1}\,e^{2\,\sqrt{b^2-c}\,x}+{\it \%k}_{2}\right)\,e^{-\sqrt{b^2-c}\,x-b\,x}\]

$c – b^2 = 0$ の場合

In [54]:
ode2(ev(eq, c=b**2), y, x);

XIs b zero or nonzero?
\Is b zero or nonzero?
nonzero;
Out[54]:
\[\tag{${\it \%o}_{83}$}y=\left({\it \%k}_{2}\,x+{\it \%k}_{1}\right)\,e^ {- b\,x }\]

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

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

In [55]:
y(mu):= exp(-b*x) * (A*cos(mu * x) + B/mu*sin(mu * x));
Out[55]:
\[\tag{${\it \%o}_{84}$}y\left(\mu\right):=\exp \left(\left(-b\right)\,x\right)\,\left(A\,\cos \left(\mu\,x\right)+\frac{B}{\mu}\,\sin \left(\mu\,x\right)\right)\]

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

In [56]:
y = y(sqrt(c-b**2));
Out[56]:
\[\tag{${\it \%o}_{85}$}y=e^ {- b\,x }\,\left(\frac{B\,\sin \left(\sqrt{c-b^2}\,x\right)}{\sqrt{c-b^2}}+A\,\cos \left(\sqrt{c-b^2}\,x\right)\right)\]

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

In [57]:
y = y(%i * sqrt(b**2-c));
Out[57]:
\[\tag{${\it \%o}_{86}$}y=e^ {- b\,x }\,\left(\frac{B\,\sinh \left(\sqrt{b^2-c}\,x\right)}{\sqrt{b^2-c}}+A\,\cosh \left(\sqrt{b^2-c}\,x\right)\right)\]

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

In [58]:
y = limit(y(mu), mu, 0), factor;

XIs b x an integer?
\Is b x an integer?
no;
Out[58]:
\[\tag{${\it \%o}_{87}$}y=\left(B\,x+A\right)\,e^ {- b\,x }\]

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

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

例題

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

In [59]:
/* a, b は正（ゼロでない）と仮定する。*/
assume(a>0, b>0)$

eq: 'diff(y, x, 2) + a**2 * y = sin(b*x);
Out[59]:
\[\tag{${\it \%o}_{89}$}\frac{d^2}{d\,x^2}\,y+a^2\,y=\sin \left(b\,x\right)\]
In [60]:
ode2(eq, y, x);
Out[60]:
\[\tag{${\it \%o}_{90}$}y=-\frac{\sin \left(b\,x\right)}{b^2-a^2}+{\it \%k}_{1}\,\sin \left(a\,x\right)+{\it \%k}_{2}\,\cos \left(a\,x\right)\]

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

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

In [61]:
eq1: ev(eq, b = a); /* ev(eq) は b = a の場合を評価してくれる。*/
Out[61]:
\[\tag{${\it \%o}_{91}$}\frac{d^2}{d\,x^2}\,y+a^2\,y=\sin \left(a\,x\right)\]
In [62]:
ode2(eq1, y, x);
Out[62]:
\[\tag{${\it \%o}_{92}$}y={\it \%k}_{1}\,\sin \left(a\,x\right)-\frac{x\,\cos \left(a\,x\right)}{2\,a}+{\it \%k}_{2}\,\cos \left(a\,x\right)\]

積分定数 $\% k_1, \% k_2 $がつくのは同次方程式の一般解。非同次方程式の特解は，振幅が $x$ に比例して単調増加していく。

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

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

まず，同次方程式は…

In [63]:
eq0: 'diff(y, x, 2) + b**2 * y = 0;
Out[63]:
\[\tag{${\it \%o}_{93}$}\frac{d^2}{d\,x^2}\,y+b^2\,y=0\]
In [64]:
ode2(eq0, y, x);
Out[64]:
\[\tag{${\it \%o}_{94}$}y={\it \%k}_{1}\,\sin \left(b\,x\right)+{\it \%k}_{2}\,\cos \left(b\,x\right)\]

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

In [65]:
y1(x):= sin(a*x);
y2(x):= cos(a*x);
Out[65]:
\[\tag{${\it \%o}_{95}$}y_{1}\left(x\right):=\sin \left(a\,x\right)\]
Out[65]:
\[\tag{${\it \%o}_{96}$}y_{2}\left(x\right):=\cos \left(a\,x\right)\]

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

In [66]:
W: y1(x) * diff(y2(x), x) - diff(y1(x), x) * y2(x);
W: trigsimp(W);
Out[66]:
\[\tag{${\it \%o}_{97}$}-a\,\sin ^2\left(a\,x\right)-a\,\cos ^2\left(a\,x\right)\]
Out[66]:
\[\tag{${\it \%o}_{98}$}-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’$$

$b \neq a$ を仮定して…

In [67]:
R(x):= sin(b*x);

ys: y2(x)*integrate(R(x)*y1(x)/W, x) - y1(x)*integrate(R(x)*y2(x)/W, x);
Out[67]:
\[\tag{${\it \%o}_{99}$}R\left(x\right):=\sin \left(b\,x\right)\]
Out[67]:
\[\tag{${\it \%o}_{100}$}\frac{\sin \left(a\,x\right)\,\left(-\frac{\cos \left(\left(b+a\right)\,x\right)}{2\,\left(b+a\right)}-\frac{\cos \left(\left(b-a\right)\,x\right)}{2\,\left(b-a\right)}\right)}{a}-\frac{\cos \left(a\,x\right)\,\left(\frac{\sin \left(\left(a-b\right)\,x\right)}{2\,\left(a-b\right)}-\frac{\sin \left(\left(b+a\right)\,x\right)}{2\,\left(b+a\right)}\right)}{a}\]

expand(ys) でかっこをはずして展開して…
trigexpand(%) で三角関数の加法定理を使い…
最後に trigsimp(%) で簡単化すると…

最後に trigsimp(%) で簡単化すると…

In [68]:
expand(ys)$
trigexpand(%)$
trigsimp(%);
Out[68]:
\[\tag{${\it \%o}_{103}$}-\frac{\sin \left(b\,x\right)}{b^2-a^2}\]

無事，非同次方程式の特殊解
$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$$

とすると…

In [69]:
assume(x>0)$
ys: y2(x)*integrate(R(t)*y1(t)/W, t, 0, x) - y1(x)*integrate(R(t)*y2(t)/W, t, 0, x);
Out[69]:
\[\tag{${\it \%o}_{105}$}\frac{\cos \left(a\,x\right)\,\left(\left(b-a\right)\,\sin \left(\left(b+a\right)\,x\right)+\left(b+a\right)\,\sin \left(\left(a-b\right)\,x\right)\right)}{a\,\left(2\,b^2-2\,a^2\right)}+\frac{\sin \left(a\,x\right)\,\left(\frac{b}{b^2-a^2}-\frac{\left(b-a\right)\,\cos \left(\left(b+a\right)\,x\right)+\left(b+a\right)\,\cos \left(\left(b-a\right)\,x\right)}{2\,b^2-2\,a^2}\right)}{a}\]

なんだかんだで簡単化してやると…

In [70]:
expand(ys)$
trigexpand(%)$
trigsimp(%)$
expand(%);
Out[70]:
\[\tag{${\it \%o}_{109}$}\frac{b\,\sin \left(a\,x\right)}{a\,b^2-a^3}-\frac{a\,\sin \left(b\,x\right)}{a\,b^2-a^3}\]

… となり，$\displaystyle y_s = \frac{\sin (b x)}{a^2-b^2}$ のほかに $\sin (a x)$ に比例する項も現れます。しかし，$\sin (a x)$ に比例する項は，同次方程式の一般解に組み込まれるので，$\sin (b x)$ に比例する項が結局非同次方程式の特殊解となります。

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

$b = a$ の場合は，非同次項 $R(x) = \sin (a x)$ として…

In [71]:
R(x):= sin(a*x)$

ys: y2(x)*integrate(R(x)*y1(x)/W, x) - y1(x)*integrate(R(x)*y2(x)/W, x);
Out[71]:
\[\tag{${\it \%o}_{111}$}-\frac{\cos \left(a\,x\right)\,\left(a\,x-\frac{\sin \left(2\,a\,x\right)}{2}\right)}{2\,a^2}-\frac{\cos ^2\left(a\,x\right)\,\sin \left(a\,x\right)}{2\,a^2}\]

expand(ys) でかっこをはずして展開して…
trigexpand(%) で三角関数の加法定理を使い…
最後に trigsimp(%) で簡単化すると…

In [72]:
expand(ys)$
trigexpand(%)$
trigsimp(%);
Out[72]:
\[\tag{${\it \%o}_{114}$}-\frac{x\,\cos \left(a\,x\right)}{2\,a}\]

無事，非同次方程式の特殊解
$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$$

とすると…

In [73]:
assume(x > 0)$

ys: y2(x)*integrate(R(t)*y1(t)/W, t, 0, x) - y1(x)*integrate(R(t)*y2(t)/W, t, 0, x);
Out[73]:
\[\tag{${\it \%o}_{116}$}\frac{\cos \left(a\,x\right)\,\left(\sin \left(2\,a\,x\right)-2\,a\,x\right)}{4\,a^2}+\frac{\left(\frac{1}{2\,a}-\frac{\cos ^2\left(a\,x\right)}{2\,a}\right)\,\sin \left(a\,x\right)}{a}\]

なんだかんだで簡単化してやると…

In [74]:
expand(ys)$
trigexpand(%)$
trigsimp(%)$
expand(%);
Out[74]:
\[\tag{${\it \%o}_{120}$}\frac{\sin \left(a\,x\right)}{2\,a^2}-\frac{x\,\cos \left(a\,x\right)}{2\,a}\]

… となり，やはり非同次方程式の特殊解の項のほかに，$\sin (a x)$ に比例する項，つまり同次方程式の一般解に組み込まれる項も現れます。