Euler法によって計算した値は，厳密な値のグラフと比べると誤差があるのがわかる。この誤差は∆xを小さ
く取ることにより，小さくすることができる。例えば，上のプログラムでDx=0.001としてみよう。（このと
き10行目のrange(1,81)を変えて，計算のステップ数も増やす必要がある。）
46.3 微分方程式の厳密解
この節ではdesolve関数を使って微分方程式を解きます。desolveの文法とオプションは次の通りです：
微分方程式を解く(desolve)
1 d e s o l v e ( de , dvar , o p t i o n s ) ここで
• de: 微分方程式(differential equation)
• dvar: 従属変数(dependent variable:求める関数）
です。オプションoptionsは次が用意されています：
• ics: 初期条件(initial conditions)もしくは境界条件を指定するときに設定します
• ivar: 独立変数(independent variable)
• show_method: 解法に使用した手法(微分方程式の型)を表示する
• contrib_ode: Clairaut, Lagrange, Riccati型の方程式も含めて解を探す場合にTrueにする。
これらのオプションは省略可能で，指定しない場合Falseになります。
さて簡単な微分方程式
y′(x) =y(x) (17)
を解いてみましょう。これ解くには次のようにします：
1 var (’ x ’) # xは 形 式 的 な 変 数
2 y = f u n c t i o n (’ y ’, x ) # yをxの 関 数 と す る!!
3 d e s o l v e ( d i f f ( y , x ) == y , y ) # 微 分 方 程 式 y ’= y を yに つ い て 解 く
実行結果の例 _C*e^x
上のプログラムの2行目で，yはxの関数であるという宣言をしています。3行目のdiff(y,x)は微分y′(x) を意味していて，diff(y,x)==yが解きたい微分方程式ということになります。実行結果の_Cは任意定数 です。
desolveでは微分方程式を指定するときに==0を省略することができます。例えば，desolve(diff(y,x)-y,
y)とdesolve(diff(y,x)-y == 0, y)は同じ意味です。
1 var (’ x ’); y = f u n c t i o n (’ y ’, x )
2 d e s o l v e ( d i f f ( y , x ) - y , y ) # 微 分 方 程 式 y ’ - y =0 を yに つ い て 解 く
実行結果の例 _C*e^x
46.3.1 微分方程式の解法1
次の微分方程式を解く事を考えましょう。
y′(x) + 2y(x) sin(x) = 0 (18)
微分方程式を勉強した事がある人は，変数分離型の解法によって直ちにこれを解く事ができると思います。忘 れた人もいるかもしれないので念のために変数分離型の解法を復習しておきます。まず上の微分方程式を
y′
y =−2 sin(x) (19)
と変形して*14，両辺をxで積分する事により logy=
∫ y′
ydx= 2 cos(x) +c (20)
となります。これをyについて解くことにより答えは
y(x) = exp(2 cos(x) +c) (21)
となります。これも一般解といえますが，これではyが恒等的に0となる自明な解を含みません。そこで，
C=ecを新たに任意定数とすると，(18)の解は
y(x) =Ce2 cos(x) (22)
となります。さて，上の微分方程式をSageで解いてみましょう：
1 x = var (’ x ’); y = f u n c t i o n (’ y ’, x )
2 DE = d i f f ( y , x ) + 2* y * sin ( x ) == 0 # 微 分 方 程 式 をD Eと 名 付 け る
3 d e s o l v e ( DE , y ) # D Eをyに つ い て 解 く
実行結果の例 _C*e^(2*cos(x))
46.3.2 微分方程式の解法2 次の微分方程式を考えます：
y′(x) =x−y(x)
x+y(x). (23)
これも変数分離型として解く事ができますが，Sageに解かせてみましょう：
*14y(x) = 0のときはどうするんだと指摘されそうですが，とりあえずy(x)̸= 0と仮定しましょう。
46.3 微分方程式の厳密解 46 微分方程式
1 x = var (’ x ’); y = f u n c t i o n (’ y ’, x ) 2 d e s o l v e ( d i f f ( y , x ) == ( x - y )/( x + y ) , y )
実行結果の例 -1/2*x^2 + x*y(x) + 1/2*y(x)^2 == _C
微分方程式(23)の解y(x)は，微分を含まない方程式
−1
2x2+xy(x) +1
2y(x)2=C (24)
を満たす関数として陰(implicitly)に解かれました。上式はy(x)についての2次関数なので解の公式を使っ て直ちに解く事ができますが，Sageのsolveにそれをやらせてみましょう：
1 x = var (’ x ’); y = f u n c t i o n (’ y ’, x )
2 SOL = d e s o l v e ( d i f f ( y , x ) == ( x - y )/( x + y ) , y ) # 微 分 方 程 式 の 答 え をS O Lと 名 づ け る
3 s o l v e ( SOL , y ) # 方 程 式S O Lをyに つ い て 解 く
実行結果の例 [y(x) == -x - sqrt(2*x^2 + 2*_C), y(x) == -x + sqrt(2*x^2 + 2*_C)]
したがって，微分方程式の答えは，y(x) =−x±√
2x2+ 2Cとなることが分かります。
±と定数Cに応じて，微分方程式の解は無数にありますが，これをグラフに描いてみましょう：
1 p = G r a p h i c s () # 空 の グ ラ フ ィ ッ ク ス オ ブ ジ ェ ク ト を 作 る
2 for C in r a n g e(4 ) : # 定 数 Cを0か ら3ま で 変 え な が ら 解 の グ ラ フ をpに 加 え る
3 p += p l o t ( - x - s q r t (2* x ^2 + 2* C ) , ( x , -3 ,3) ,
4 c o l o r = hue ( C /4 ,1 ) , l e g e n d _ l a b e l =’ C = ’+str( C )+’ , $ + $ ’) # Cご と に 色 を 変 え る 5 p += p l o t ( - x + s q r t (2* x ^2 + 2* C ) , ( x , -3 ,3) ,
6 c o l o r = hue ( C /4 ,0.5) , l e g e n d _ l a b e l =’ C = ’+str( C )+’ ,$ - $ ’)# Cご と に 色 を 変 え る 7
8 p . s h o w () # pを 描 画
-3 -2 -1 1 2 3
-5 5
C=0, + C=0, − C=1, + C=1, − C=2, + C=2, − C=3, + C=3, −
図34 出力結果：微分方程式(23)の解のグラフ
46.3.3 微分方程式の解法３
次の微分方程式を考えます：
y′= sin(x+y) + sin(x−y) (25)
これも変数分離型で解く事ができますが，このままではSageは解く事ができません，
1 var (’ x ’); y = f u n c t i o n (’ y ’, x )
2 d e s o l v e ( d i f f ( y , x ) == sin ( x + y )+ sin ( x - y ) , y )
実行結果の例 NotImplementedError: Maxima was unable to solve this ODE. Consider to
set option contrib_ode to True.
右辺を加法定理で因数分解してやれば解く事ができます：
1 var (’ x ’); y = f u n c t i o n (’ y ’, x )
2 DE = d i f f ( y , x ) == sin ( x )* cos ( y )+ sin ( y )* cos ( x )+ sin ( x )* cos ( y ) - sin ( y )* cos ( x ) 3 d e s o l v e ( DE , y )
実行結果の例 -1/4*log(sin(y(x)) - 1) + 1/4*log(sin(y(x)) + 1) == c - cos(x)