変数を分離しよう
★方向ベクトルたちをばらまく
０．微分方程式は関数さがし
このワークシートはMath by Codeの一部です。
＜微分方程式の解とは？＞
普通、「方程式」と言えば、等式を満たす数を求めるための情報のことだ。
平たく言えば「逆算」だね。
（マスコミがよく使う「勝利の方程式」のような「法則、公式」という意味ではない。）
たとえば、「かけ算の方程式」を解く、つまり、かけ算の逆算をするには、かけ算の逆、割り算をする。
だから、「微分の方程式」を解く、つまり、微分の逆算をするには、積分をすることになりそうだね。
代数方程式ならば、「代数の等式」という情報を満足させる「数」を求めることがゴールになる。
微分方程式ならば、「微分の等式」という情報を満足させる「解」を求めることがゴールになるでしょう。「微分方程式の解」は「数」ではなく、「関数」です。
代数方程式を解くためには解の公式、因数定理、因数分解、置換による次元下げなど、
数と数、数式と数式の「関係性を利用した変形」が重要だったね。
微分方程式を解くためには、関数と関数の「関係性を利用した変形」が重要になるでしょう。
たとえば、積分、部分積分、置換積分、級数展開、フーリエ・ラプラスの変換（と逆変換）などが微分積分の関係性を利用した変形だね。
たとえば、f(x)=１という定数関数を不定積分すると、F(x)=x+Cという関数になる。
だから、「dy/dx=1」は微分方程式だ。f(x)=x+Cが解だというのは形式的にわかる。
ということで、積分の基礎知識が必要になりそうだから、確認しておこう。
＜微積の基本公式の確認＞
（微分は->、不定積分は=>として定数Cは省略する）
べき関数 xn=>1/(n+1) * xn+1
三角関数 sinx => -cos x , cos x =>sinx
sin(ax) => -1/a cos(ax) , cos(ax)=> 1/a sin (ax)
指数対数系 1/x => loge|x|
e ^ax => 1/a * e ^ax
・合成の微分 f(g(x)) -> df/dg *dg/dx
からの置換積分
g=g(x) とすると、f f(g(x)) dg/dx dx = ∫f(g) dg => F(g)
(例)g=ax+b f(g)=e^gとみると、 ∫ae^(ax+b) dx = e^g + C
・積の微分 fg -> f'g + fg' から、 fg=∫(f' *g + fg')
からの積の積分（部分積分） ∫f' *g dx = - ∫f g' dx + f*gまたは、∫f *g' dx = f*g- ∫f' g dx
g',gを格上げしてg,Gとかくと、∫f *g dx =f*G - ∫f' G dx
積f*gの片方だけ、積分形にして、もう片方は積分の中で微分する。
テーラー展開
e^x=Σ1/n! * x^n
sin x = Σ (-1)^n/ (2n+1)! *x ^(2n+1)
cos x = Σ (-1)^n/ (2n)! *x ^(2n)
オイラー等式
e^ix= sin x + i cos x
＜微分方程式を解くことの意味は？＞
では、これってどういう意味なのだろうか。
それが、上のアプレットからわかる。
図形的にみると、「dy/dx」は関数y=f(x)の曲線の点(x,y)における接線を表す関数y=f'(x)だった。
つまり、これって、ｘｙ平面でのf(x)の変換の仕方についての一般的な情報ともいえるね。
この意味からすると、「dy/dx=1」は接線の傾きが一定の１になるのが関数y=f(x)ということになる。
つまり、f(x)はどこでも傾き１になる曲線、つまり、傾き１の直線ということだから、
一般的な解はf(x)=x+Cになるね。この関数群が一般解だ。
ｘ＝０のときのf(0)=0ならば、f(x)=xと限定できる。これが特殊解になる。
出発点がｘ＝０で、変化の仕方がdy/dx=1という特徴をもつ関数がf(x)=xだと限定できると、
ｘ＝０以外でのf(x)のふるまいが予測できるというのが微分方程式を解くことのよさになるでしょう。
代数方程式の場合は、1次、2次、3次と次数をメインに種類分けし、
種別に対する解の公式や解法を探るということをしてきた。
微分方程式はどうやってタイプ分けし、どう対処するのでしょうか。
それを探っていきましょう。
１．微分方程式のタイプ
微分方程式についていろいろなネーミングやタイプ分けがあるので、
アレルギーが起きないように、さらっと外観しておこう。
＜種別＞
（階数）
代数方程式では、最高の次数がNならば、N次方程式と呼んだように、
微分方程式では、最高の階数（微分の回数）がNならば、N階方程式と呼ぶ。
・常微分方程式(ODE)は、yをｘで何回か微分したものを要素としている。
yが1変数ｘだけの関数のときだね。
・偏微分方程式は、それに対して多変数関数を微分した方程式。
yを3回微分したものをy(3)と書いたり、y'''と書いたりする。時間ｔの関数を微分3回微分したときは。
ｔの関数にドットを３つ乗せて表したりすることもあるね。
（例）
・yのｎ階微分をy(n)と書くことにすると、y(3)+y(2)+y(1)=f(x)+g(y)は3階
・y(2)=-ayは2階。yは単振動を表すことがある。
・上記のy(1)=-1は1階。
（線形／非線形、同次／非同次、定数）
・yのｎ階微分をy(n)と書き、Akがxの関数、ΣAk*y(k)=0のとき、線形同次。
Akが定数なら、特に、定数線形同次というね。
・yのｎ階微分をy(n)と書き、Akがyの関数だったり、ΣAk*(y(k))j=0のように、
yの微分の項で1次（線形）でないものがあるとか、yの微分同士の積のが、非線形同次。
・さらに、＝０ではなく、=B(x)のとき、それぞれ非同次という。だから、
線形非同次、非線形非同次などもあるね。
（例）
（線形同次) y(2)+5y(1)=0, xy(2)+x2y(1)=y
（線形非同次）xy(2)+x2y(1)=x
（非線形同次) y(2)+yy(1)=0, (y(1))2+x2y=0
（非線形非同次）y(1)y(3)+xy(2)=3x2
＜基本型＞
直接型、変数分離型、同次型、その他、などタイプ分けがある。
・直接型dy/dx=f(x)
右辺の不定積分で、y=F(x)＋Cと一般解を求めることができるね。
・変数分離型dy/dx=f(x)g(y)
1/g(y) dy = f(x) dxのように、変数を等式の左右に分離できる型。
両辺に積分記号をつけた∫ 1/g(y)dy= ∫f(x)dxから、G(y)=F(x)+C。
・同次型dy/dx=f(y/x)
u=y/xを代入することで、変数分離形に還元できる。
微分方程式。変数分離y'=-2xy
微分方程式。変数分離y'=ｘ(2y-1)
微分方程式。変数分離y'=-4x/y
微分方程式。同次型y'=(x^2+y^2)/2xy
微分方程式。同次型y'=(2x-y)/(x-2y)
２．変数分離による解法
＜変数分離型dy/dx=f(x)/g(y)＞
dy/dx=f(x)/g(y)を代数的に変数分離したf(x)dx=g(y)dyの両辺積分∫f(x)dx＝∫g(y) dyから解が求められる。
dy/dx=f(x)/g(y)でｙ＝h(x)とすると、dy/dx=h'(x)=f(x)/g(h(x))だからf(x)=h'(x)g(h(x))
したがって積分は∫f(x)dx=∫h'(x)g(h(x))dx。
置換積分dy= dy/dx dx= h'(x)dx より、右辺＝∫ g(y)h'(x)dx= ∫g(y) dy。だから、∫f(x)dx＝∫g(y) dy。
（例）
「dy/dx=-2xyの一般解と(x,y)=(0,1)を通る特殊解」は？
1/ydy=-2xdxの両辺積分∫1/ydy=∫-2xdxから、log|y|=-x2＋C |y|=e(-x^2+C)
だから、y=±e(-x^2+C) 。一般解はy=De(-x^2)、(x,y)=(0,1)からD=1で特殊解はy=e(-x^2)
（例）
「dy/dx=x(2y-1)の一般解と(x,y)=(0,1)を通る特殊解」は？
1/(2y-1)dy=xdxの両辺積分1/2∫2/(2y-1)dy=∫xdxから、1/2log|2y-1|=1/2x2＋C |2y-1|=e(x^2+2C)
だから、y=±1/2e(x^2+2C) +1/2。一般解はy=De(x^2)+1/2、(x,y)=(0,1)からD=1/2で特殊解はy=(e(x^2)+1)/2
（例）
「dy/dx=-4x/yの一般解と(x,y)=(0,2)を通る特殊解」は？
ydy=-4xdxの両辺積分∫ydy=∫-4xdxから、 1/2・y2=-2x2＋C 一般解は(x/1)2+(y/2)2=C/2。
(x,y)=(0,2)からC=2で、特殊解は(x/1)2+(y/2)2=12の楕円。ｙ＝０除く
＜同次型dy/dx=f(y/x)＞
「dy/dx=f(y/x)ならu=y/xとして変数分離で両辺積分∫(1/x) dx =∫1/(f(u)-u) duで、u=y/xに戻せる。」
u=y/xとおくとdy/dx=f(u)。y=xuの両辺微分dy/dx=(xu)'=(x)'u+x(u)'=u+x・du/dx。
・2式からf(u)=u+du・x/dx 。 x/dx=(f(u)-u)/du。(1/x) dx =1/(f(u)-u) du
両辺積分∫(1/x) dx =∫1/(f(u)-u) duして、uをy/xに戻せば一般解になるね。
・ y=mx型の解があるときはm=y/x=u（定数）。dy/dx=f(m)=mを満たすmをy=mxに代入すればよい。
（例）
「dy/dx=(x2+y2)/2xyの一般解」は？
u=y/xとおくとdy/dx=(1+(y/x)2)/(2y/x)=(1+u2)/2u。y=xuの両辺微分dy/dx=(xu)'=u+x・du/dx。
２式からu+x・du/dx=(1+u2)/2u。x/dx={(1-u2)/2u} /duこれから2u/(u2-1) du=-1/x dxと変数分離。
から、 となり、
だから、y2-x2=Cx。C=0の場合、y＝±xも解（特異解）となる。C=0も含めた一般解はy2-x2=Cx
（例）
「dy/dx=(2x − y)/(x − 2y)の一般解」は？
u=y/xとおくとdy/dx=(2-y/x)/(1-2(y/x))=(2-u)/(1-2u)。y=xuの両辺微分dy/dx=(xu)'=u+x・du/dxから
u+x・du/dx=(2-u)/(1-2u)。x/dx={(2-u-u+2u2)/(1-2u} /duから(1-2u)/2(1-u+u2) du=1/x dxと変数分離。 から、 となり、 だから、x2-xy+y2=C (C=0の場合も含む）
３．実装例と計算例
質問：係数分離で解ける微分方程式をコードで解くにはどうしたらよいでしょうか。
科学技術の計算と視覚化の用途では、pythonかjuliaが適しています。
pythonユーザーが多いので、pythonを使います。
javascrptはrustやC++で作った部品をwebassenble化してweb表示するときには役立ちますが、高速計算や
大量データ処理が必要なければ、pythonにsympyをimportするのがよいでしょう。
sympyは微積分の機能があり、sympy.abcは変数を関数として設定する機能があります。
y = Function('y')
これで、yの文字や関数として読み取ります。
しかし、dy/dxはこのままでは読み取れないのでDerivative(y(x), x)とかきましょう。
また、微分方程式がA=Bという等式だとしたら、
Eq(A,B)で渡します。それに名前をつけましょう。たとえば、deq。
微分方程式の解はdsolve(方程式、y(x))とかくと、「y=」の形式の解を返してくれます。
ただし、注意点があります。
A,Bによませる式では、yはy(x)とかいて、xの関数であることをいちいち明示することです。
やってみると、変数分離型でも同次型でも解けます。
y^2=Pとなる関数が解になるときは、y=の形の回答となり、わざわざ+sqrt(P),-sqrt(P)の２つを返してくるので見た目が期待とちがったりしますね。
[IN]=======================================
from sympy.abc import *
from sympy import *
y = Function('y') # y=f(x)
def dydxEq(fx):
Eq(左辺,右辺)
deq = Eq(Derivative(y(x), x), fx)
return dsolve(deq, y(x))
[IN]=======================================
dydxEq(-2*x*y(x)) #一般解はy=De(-x^2)
[OUT]=====================================
[IN]=======================================
dydxEq(x*(2*y(x)-1)) #一般解はy=De(x^2)+1/2
[OUT]=====================================
[IN]=======================================
dydxEq(-4*x/y(x)) #一般解は(x/1)^2+(y/2)^2=C/2
[OUT]=====================================
[Eq(y(x), -sqrt(C1 - 4*x**2)), Eq(y(x), sqrt(C1 - 4*x**2))]
[IN]=======================================
dydxEq((x**2+(y(x))**2)/(2*x*y(x))) #一般解はy^2-x^2=Cx
[OUT]=====================================
[Eq(y(x), -sqrt(x*(C1 + x))), Eq(y(x), sqrt(x*(C1 + x)))]
変数分離による解法とコードをみてきましたが、
現実問題としてどんな場面で
変数分離できる例が登場するかをみておきましょう。
＜日本の人口＞
時刻 t におけるある国の総人口を N(t) とする。
短い時間間隔 dt における出生数と死亡数は総人口と時間間隔に比例するから，
出生数=αNdt, 死亡数=βNdtとおける。
時間間隔 dt での人口の増加 dN=αNdt–βNdt=γNdtとかけるね。
これは微分方程式dN/dt=γN
∫1/N dN= ∫γ dt から、logN=λt+C
日本の令和５年(t=2023)出生率6.0/1000, 死亡率13.0/1000 ,人口N=121 193 394
γ=(6.0-13.0)/1000=-0.007。
log121193394=-0.007*2023+ C より、C=log121193394+0.007*2023=
import math
math.log(121193394)+0.007*2023
32.77389812516354
logN=λt+32.77とすると、N=e^(-0.007t) *e^32.77
N(2050)=math.exp(-0.007* 2050)*math.exp(32.77)
=99931948.77012336
≒99931948
日本の人口は、2023年の1億2千万人から、2050年は9993万人になる？
そんなことはないかもしれません。
国土審議会によると9515万人と推定されています。
すでに高齢化社会なので、死亡率が高い人の割合が多く含まれるからでしょう。
＜放射性物質の減少率＞
放射性物質の量をN、時間ごとの減少率をγとすると、
人口モデルと同型の微分方程式dN/dt=-γN(γ>0)ができる。
∫1/N dN= -∫γ dt から、logN=-λt+c、N＝e^c*e^(-λt)=Ce^(-γt)となる。
半減期は、Ce^(-γx)=1/2Ce^(-γt) から、e^(-γx)=2^-1 *e^(-γt)。対数をとると、-yx=-log2- γt
log2=γ(x-t) から、半減する年数x-t=log2/γ
2011東日本大震災で福島第一原子力発電所から放出された主な放射性物質ヨウ素131、セシウム134、セシウム137、ストロンチウム90の4種類の半減期は順に約8日間、2年、30年、29年とします。
このうち、セシウム137の半減期が一番長いので、この物質が2011年に放出した量をN(2011)として、
2050にどれだけの量が残っているかを計算してみよう。
半減する年数x-t=log2/γ=30とする。γ=log2/30となる。N(2011)=Ce^(-log2/30 * 2011).
これから、C=N(2011)/e^(-log2/30 * 2011)
N(2050)=Ce^(-log2/30 * 2050)=N(2011)/e^(-log2/30 * 2011) * e^(-log2/30 * 2050)
=N(2011)*e^{(-log2/30 * 2050)-(-log2/30 * 2011)}
=N(2011)*e^{(-log2/30 * 39)}
import math
math.exp(-math.log(2)/30*39)
=0.4061261981781178
まだ、約40%残るという計算です。
ただ、その場所にとどまっているわけではないので、拡散して濃度が薄まれば数字よりは付近にある量は
もっと減っているでしょう。
＜地球脱出速度＞
地球の表面から発射された飛行体が地球の重力から逃れるための最小の初速ｖを求めたい。
地球、飛行体の質量をM,mとして、地球の中心からの飛行体への距離をrとしよう。
地球の半径をR,重力をg、万有引力定数をGとすると、
この飛行体が受ける力F=-ma=GMm/r2から、加速度a=-GM/r2となる。
発射時は、加速度はaを-gにおきかえ、rをRにおきかえて、g=GM/R2から、GM=gR2となる。
２式をまとめると、a=-(gR2)/r2=-g(R/r)2となるね。
ここで、v=dr/dt, a=dv/dtだから、微分方程式はa=dv/dt=dv/dr*dr/dt=
v *dv/dr =-gR2/r2
変数分離して、v dv=-gR2/r2 dr だから両辺積分∫v dv=-gR2∫1/r2 dr
1/2v2=gR2/r+c これから、c=1/2v2-gR2/r
発射時点では、v=v0, r=Rだから、c=1/2 v02-gR
無限遠点r=∞で、v=0になるから、c=1/2v2-gR2/r→0-0=0。c=0
まとめると、c=1/2 v02-gR=0から、v0=√(2gR)
R：地球の平均半径 6372(km)
g：地球の重力 0.00980 (km/s2)
R=6372
g=0.00980
math.sqrt(2*g*R)
11.175473144346059
v0=11.2(km/s)の初速で、地球脱出ができるね。
質問：ここまでの現実的な問題の解に対して、具体的なイメージがわくような視覚化はどうやればよいでしょう。
微分方程式の解となる関数をgeogebraで表示しましょう。
独立変数となるｔやｒをｘにした式にすることが必要だね。
独立変数をスライダーなどで変化させることで、求めた計算結果などの確認ができるでしょう。