2023年12月21日
今日の活動
気分が浮かれているのでCoCo壱番屋で骨付きチキンがはいったスープカレー食べたんだけど、チキンがスープにつかって、これどうやって食べんねんと大変だった。
スープカレーは美味しかった。
今日のわんこ
寝たふりしたら何するんだろうと思って試してみたら、おもちゃもって走り込んできた。
Advent Calendar 21日目
数学IIIを学んでるとき、数IIIの時点でもうラスボスなのに、微分方程式となるとラスボスが一段階変化した、その上で高校から一歩先に出たところに存在するので裏ボスな存在に感じなかった？
おさらい
\(y'(x) = y\)となる微分方程式を考える。初期値は\(y(0)=1\)とするとき、\(y(1)\)を求めたい。
答え: \(y(x) = e^x\)より\(y(1) = e \approx 2.718 \)
MyClassで微分方程式を解くと言っても、昨日のがほぼ全てである。というのも昨日のは微分方程式を加減乗除だけで計算していた。これはMyClassでも計算できる。
\(x=0\)のときから変化した\(x=0+\Delta x\)を考える。MyClassでいうと、MyClass([0,1,0,0,0])だ。
もう一度ルンゲクッタ法を利用する。
x_o = MyClass([0,1,0,0,0])
k_1 = x_o
k_2 = x_o + 0.5*k_1
k_3 = x_o + 0.5*k_2
k_4 = x_o + k_3
print((k_1 + 2*k_2 + 2*k_3 + k_4)*(1/6))
# [0, 1.708333, 0, 0, 0]
これは、\(x=0\)から\(\Delta x\)だけ変化した場合、yはy(0)から\(1.708333\Delta x\)だけ変化しますよと言っている。
\(\Delta x\)がどれぐらいだと正確かはわからないが、y(0)から1変化させたy(1)の値が知りたいので、\(y(1) = y(0) + 1.708333\)である。
つまり、\(y(1) = 2.708333\)で、まあ昨日やったルンゲクッタ法と同じ答えだ。
これだと昨日と何もかわらないから、もう一個別の方法でやろう。せっかく2日前に積分を定義したのだ、微分方程式をピカールの逐次近似法により積分方程式に書き換える。なお何故逐次近似法を使っていいのか(リプシッツ条件がどうとか連続関数かどうとか)はうまいこと考えてほしい。僕には成立するような数式しか扱わないので興味がないのだが、知りたい人には本の情報を教える。
x_o = MyClass([1,0,0,0,0])
def f(x_o, x_n):
return x_o + (x_n).integrate()
f(x_o, x_o)
# [1, 1, 0, 0, 0]
これは、\(x=0\)から\(\Delta x\)だけ変化したとき、yはy(0)から\(1\Delta x\)だけ変化しますよと言っている。つまりy(1)=2とこのコードは主張している。が、当然そんな精度を今更受け入れることはできない。
ところで逐次近似法は漸化式として、もう一度積分方程式をつくるとよりよい近似になるという話であった。漸化式ということは、再帰的に書けると聞いた記憶があるので、こうなる。
x_o = MyClass([1,0,0,0,0])
def f(x_o, x_n, n):
if n == 1:
return x_o + (x_n).integrate()
else:
return x_o + f(x_o, x_n, n-1).integrate()
f(x_o, x_o, 1), f(x_o, x_o, 2), f(x_o, x_o, 3), f(x_o, x_o, 4), f(x_o, x_o, 5)
# [1, 1, 0, 0 ,0]
# [1, 1, 0.5, 0, 0]
# [1, 1, 0.5, 0.166667, 0]
# [1, 1, 0.5, 0.166667, 0.0416667]
# [1, 1, 0.5, 0.166667, 0.0416667]
n=4とn=5は同じ答えになったがとりあえず今は気にせず、[1, 1, 0.5, 0.16666667, 0.0416667]について考える。
これは、\(x=0\)から\(\Delta x\)だけ変化したとき、yはy(0)から\(1 \Delta x + 0.5 \Delta x^2 + 0.16666667 \Delta x^3 + 0.04166667\Delta x^4\)だけ変化しますよと言っている。
今回\(\Delta x = 1\)であり、初項の1はy(0)を意味しているので、結局、\(y(1) = np.sum([1, 1, 0.5, 0.16666667, 0.0416667]) = 2.708333\)となる。
さっきのルンゲクッタ法と同じ答えになったが、あのルンゲクッタ法は4次4段と呼ばれているもので、要は4次の精度で計算していた。今回のMyClassも4次まで計算したので、同じ答えになった。という話である。
よって、
x_o = MyClass([1,0,0,0,0,0])
def f(x_o, x_n, n):
if n == 1:
return x_o + (x_n).integrate()
else:
return x_o + f(x_o, x_n, n-1).integrate()
f(x_o, x_o, 5)
# [1, 1, 0.5, 0.166667, 0.0416667, 0.00833333]
であり、
np.sum(f(x_o, x_o, 5).arr)
# 2.716666
と、精度が上がるわけである。また、3次や4次の係数が目に見えるのも結構嬉しい。よね？
ちょっとさっきの再帰関数を弄る。
x_o = MyClass([1,0,0,0,0])
def f(x_o, x_n, n):
if n == 1:
return x_o + sp.symbols("x")*(x_n).integrate()
else:
return x_o + sp.symbols("x")*f(x_o, x_n, n-1).integrate()
f(x_o, x_o, 4)
# array([1, x, x**2/2, x**3/6, x**4/24], dtype=object)
これ、この微分方程式の答えが\(1+x+x^2/2+x^3/6+x^4/24\)であると主張している。ちょうど\(y(x) = e^x\)をテイラー展開した形になる。直接\(e^x\)を求めることはできないが、数値解として十分なことはできている。
類題。\(y'(x) = y^2, y(0) = 1\)のとき。聡明な読者は見ただけで\(y(x) = \frac{1}{1-x}\)と気づいていらっしゃるかもしれないが、
x_o = MyClass([1,0,0,0,0])
def f(x_o, x_n, n):
if n == 1:
return x_o + sp.symbols("x")*((x_n)**2).integrate()
else:
return x_o + sp.symbols("x")*(f(x_o, x_n, n-1)**2).integrate()
f(x_o, x_o, 4)
# array([1, x, x**2, x**3, 1.0*x**4], dtype=object)
より、やっぱり\(y(x) = \frac{1}{1-x}\)をテイラー展開した形を求めることができる。嬉しいね。
全体を通しての参考文献
平山 弘, Vol.2019-HPC-169 No.2(pdf)
ここまで3週間、読者の皆様お疲れ様でした。僕もお疲れ様でした。明日からはクリスマスが終わった前提のブログを書きます。
業務に直接関係ない、書ける範囲でつらつら書いたわけだけど、毎日帰宅して1時間近く考えるのは自分の理解にだいぶ役だった。読者には、「微積なんて卒業したら一生使わないと思ってたけど、未だに使ってる人って実在するんだ」と思ってもらえればそれでいい。
コメントを投稿
注意事項- 承認を受けるまでコメントは公開されません。
- 非公開を希望した場合、コメントは公開されません。
