オイラー法、修正オイラー法、そして、２次のルンゲ・クッタ法 [微分積分]
オイラー法、修正オイラー法、そして、２次のルンゲ・クッタ法
§１ オイラー法とその誤差
微分方程式の初期値問題
があり、
とする。
オイラー法は、微分方程式を積分した
の右辺の積分を
と近似し、これから得られる
を
とし、i=0を起点として、この漸化式から得られる
を微分方程式の解
の近似解としたものである。
また、これは次のように考えることができる。
関数y(x)を点
でテーラー展開すると
が得られる。
右辺の最後の項は２次なので
とおくと（記号Oはランダウのビッグオーと呼ばれるものだが、誤差のオーダーと考えてもよい）、
となる。
ここで、
という演算規則を使っている。
O(h²)は２次だから、これをhで割れば１次減って２−１＝１次なのでO(h)になると考えてよい。
これと同様に次の関係が成立する。
より一般に、
（４）を微分方程式（１）の右辺に代入すると、
さらに、
と置くと
となり、
この式からこの近似式の誤差の程度、誤差のオーダーがh²であることがわかる。
なお、この計算では、ビッグオーの計算規則
を使っている。
微分方程式の初期値問題
を、h=0.1とし、オイラー法を用いて数値的に解いた結果は次の通り。
§２ 修正オイラー法
§１で書いたように、オイラー法は、微分方程式
の左辺の微分を
で近似し、これから
あるいは
を得たもので簡単ではあるが、反面、精度が悪い。
そこで、より計算精度のよい方法について考察する。
微分方程式（１）の両辺を積分すると
となるが、
を満たすα、Βが存在すると仮定する。
y(x)を点で２次までテーラー展開すると
であるから、両辺をxで微分すると
これらを（７）式に代入すると、
（６）式の中に現れる
を点
でテーラー展開すると
これを代入すると（２）式に代入すると
（８）と（９）より、
したがって、
そして、これから次の漸化式を得る。
微分方程式の初期値問題
を、h=0.1とし、修正オイラー法を用いて数値的に解いた結果は次の通り。
微分方程式の解
とよく一致してーーグラフ上では厳密解と修正オイラー法による数値解の曲線が重なっているーー、計算精度が向上していることがわかる。
§３ まとめ、さらに、２次のルンゲ・クッタ法
何やら難しいことを書いたが、実は、簡単に（１０）、（１１）式を導くことができる。
y=y(x)の関数なので、
もxの関数である。
そこで、
とおくと、微分方程式（１）は次のようになる。
この両辺を積分すれば、
オイラー法はこの右辺の積分を
と近似したものだった。
これに対して、修正オイラー法は台形公式を用いて積分を近似したものだ。すなわち、
は未知数なので（オイラー法で）、
と近似すれば、
になるってわけ。
また、積分を台形公式でなく中点公式で
で近似し、さらに
と近似すれば、
そして、これから、次の漸化式を得る。
これを２次のルンゲ・クッタ法という。
（局所的な）誤差の程度は、修正オイラー法と同じくO(h³)である。
ネムネコ法（笑・実はただのテーラー展開）に明るい未来はあるか [お前らに質問]
ネムネコ法（笑・実はただのテーラー展開）に未来はあるか
ネムネコ法（笑）を用いて、次の初期値問題をh=0.1で解くことを考える。
とおく。
２次のネムネコ法（笑）は、（１）式と、 テーラー展開
から得られる次の漸化式
を用いて前進的に解いてゆく方法である。
そのため、
における１階微分
と２階微分
の値が必要になる。
における１階微分の近似値
については（１）から
さらに、（１）を微分すると、
となるのでにおける２階微分の近似値
は
これを（２）に代入すると、
となり、i=0を計算の起点にし、
と前進的に次々と計算することが出来る。
h=0.1で、初期条件が
だから、x₀=0なので、
といった具合に微分方程式（１）の初期値問題を数知的に解くことが出来る。
オイラー法の場合、
とおくと、
となるので、
となる。
微分方程式（１）の解は
だから、
ネムネコ法（笑）とオイラー法との計算精度の差は歴然だケロ。
しかも、ネムネコ法（笑）は簡単に精度の向上を図ることが出来る。
のテーラー展開を３次の項まで取ると
である。
３次導関数は
だから、漸化式は
したがって、
そして、
とすれば、（多分）あの有名な４次のルンゲ・クッタを越える高精度の計算が可能になるのであった。
このとき、５次のネムネコ法（笑）の計算結果は
CASIOさんの高精度計算サイトの結果によると、４次のルンゲクッタは
で、５次のネムネコ法（笑）の計算結果に軍配が上がるのであった。
このように、ネムネコ法（笑）はいくらでも簡単に高精度化が可能で、（しかも、打ち切り）誤差の解析ーーただのテーラー展開なのだから当たり前！！ーーも容易。
しかし、その反面、問題ごとに使用する漸化式を自分で作らなければならず、プログラムの自動化には向かない。この点は素直に認めなければならい。
しか〜し、微分方程式の初期値問題の表計算ソフトを用いた数値解法には向いているのかもしれない。
では、例によって、お前らに問題。
問題１ 微分方程式
に関して、次の問に答えよ。
（１） 一般解を求めよ。
（２） x=0のときy=0とする特殊解を求めよ。
問題２ 微分方程式
に関して、次の問に答えよ。
（１） 微分方程式の解を求めよ。
（２） h=0.1とし２次のネムネコ法（笑）を用いて、x=0.1〜0.5まで数値的に解け。そして、計算結果を比較せよ。
できたら、オイラー法による数値解も求め、計算結果を比較せよ。
ちなみに、CASIOさんのところで４次のルンゲ・クッタ法を用いて解くと、次の答えになる。
男なら、必ず、やれよな！！
このブログの訪問者に女性はいないと思うけれどーーリケジョはいるかもしれないが、ネムネコの頭の中で、リケジョは女性と別のカテゴリーに分類される(^^ゞ。だって、ネムネコは、リケジョは自分と同じ種族だと思っているから・・・ーー、
女もそうさ、見てるだけじゃ、始まらない♪
ネムネコ法（笑）による斜方投射の近似計算 [ねこ騙し物理]
ネムネコ法（笑。実はただのテーラー展開）による斜方投射の近似計算
前回、連立微分方程式
を、γ=0.3、θ=π/3の場合について、オイラー法を用いて近似的に解いた。
その計算結果は、次のとおり。
この計算結果を見ると、最高到達高さHの近辺でオイラー法による計算結果と厳密解との差異が少し大きい。
そこで、オイラー法
を改良し、近似計算の精度向上を目指すことにする。
まず、u(t+Δt)を
とテーラー展開する。
微分方程式より
したがって、
と近似でき、これから次の漸化式を得る。
同様に、
これから、次の漸化式を得る。
改めて、まとめて書くと、
k=0,1,2,・・・について
が連立微分方程式（１）、（２）の数値解になる。
この方針にしたがって、先と同じ条件で計算した結果は次の通り。
近似計算の精度が劇的に向上し、良好な結果が得られていることがわかる。
オイラー法を用いた斜方投射の近似計算 [ねこ騙し物理]
オイラー法を用いて斜方投射のシミュレーションをしてみた。
角度θで球（砲弾）を初速度
で射出したとすると、運動方程式は
初期条件は
とする。
ここでu、vは、それぞれ、x軸とy軸の正の方向の速度であり、また、mは質量、gは重力加速度、そしてk>0は空気抵抗の比例定数である。
m、g、kといった文字があると計算が厄介。
そこで、適当な無次元化を行うと、上の微分方程式は次の形の微分方程式になる。
このとき初期条件は
この微分方程式の解は、簡単な計算から
と求められ、オイラー法を用いて、この連立微分方程式をわざわざ数値的に解く必要はないのだが、それを言ったら身も蓋もない。
だから、このような批判をしてはいけないにゃ。
さて、
を積分すると
になる。
この右辺の積分を
と近似すると
などとやるのが正式なのだろうが、この方法はすこし面倒なので、ここでは採用しないことにする。
テーラー展開から
になる。
この微分方程式の場合
となるので、この結果を代入すると
だから、下の式の右辺で
とu(t+Δt)を近似したときの誤差は
であることがわかる。
でもあるので、（４）式は
さてさて、
とし、
とし、
から得られる、
の値で
における微分方程式の解
の値で近似する。
において
であるとき、（３）と（４）式から
これを、離散化した時の（局所的な）打ち切り誤差などと呼ぶことがあるようだにゃ。
同様に
という漸化式が得られ、
を起点にし、前進的に
などの前進的にその値を次々と求めることができる。
これがオイラー法を用いた連立微分方程式（１）、（２）の数値的な解法である。
γ=0.3、θ=π/3=60°、Δt=0.1として解いた結果は次のとおり。
uとxはよく合っているのですが、vとyが少し合わない。
Δt=0.1という粗い計算でここまで合ってれば十分といえば十分。
着弾点はほとんど一致しているので、敵を大砲で吹き飛ばすには十分な結果だろう。
でも、なんか気にらないケロ。
そこで、近似計算法を改良すると――仮称：ネムネコ法（笑）――、次のように良好な結果が得られる。
精度が劇的に向上し、ネムネコ法（笑）による計算結果と厳密解との違いを視覚的に見つけることはできないケロ！！
なお、投射角θ=60°なのに、そう見えないのは、x軸とy軸の目盛の間隔が１：１でないため。もしこれが１：１ならば、６０°で投げ出されていることがわかる。
この計算モデルの場合、
無次元重力加速（初速、質量、空気抵抗の比例定数の違いは全てここに反映される）
の場合、θ=６０°で砲弾を撃つと、砲弾はほとんど真上から落ちてくるんだね。知らなかった。
ちなみに、この無次元重力加速度γと投射角度θが同一であれば、同じ軌道を描く！！
それはそれとして、この計算結果を見ると、uは0に、vは−γ=−0.3に収束することがわかる。
この値は微分方程式
に、du/dt=0、dv/dt=0を代入すれば出てくるんだケロよ。
つまり、終端速度を求めるだけならば、微分方程式を解く必要がない。
知っていたケロか？
このままでは癪なんで、修正オイラー法、ルンゲ・クッタなんて高級な計算方法――計算法の導出が１変数の微分積分の範囲を越える！！――は使わずに、計算精度を向上させる方法について考え、スプレッドシートを作るにゃ。
斜方投射の動画をBloggerにアップ！！ ついでに数値的に解いてみた [ねこ騙し物理]
試しにオイラー法で解いてみた
「この微分方程式くらいならば、オイラー法で十分じゃねぇ」と思ったので、
今日の記事で取り上げた
の近似解（数値解）を、最も簡単なオイラー法で求めてみた。
その結果は、コチラ↓。
とりあえず、これだけ合っていれば十分だろう。
上の連立微分方程式の場合、オイラー法を用いると、、時刻tにおけるu(t)、x(t)の値を元に、時刻t+Δtのu(t+Δt)とx(t+Δt)が
という簡単な式で表される。
Δtを一定にし、
とすると
ここからは余談だけれど、
①は次のように変形すると、
となり、そして、Δtは一定と仮定しているので、数列
は初項
、公比1−Δtの等比数列であることわかる。
したがって、
で、u₀=1とし、区間[0,1]をn等分すると、
となるので、
そした、n→∞の極限を取ると
で、
の解は
だから、t=1とすると、
となり、この両者は一致するのであった。
ということで、
もうひとつの微分方程式
も、オイラー法で、結構、いい精度で解くことが出来そうだね。
こちらの場合は、
話を簡単にするために、
としたけれど、
角度θで投げ出す斜方投射の場合、初期値は
だケロよ。
そして、γは無次元化された重力加速度であることに注意。
この無次元重力加速度に質点の質量、空気抵抗の比例定数などがすべて集約されているんだケロよ。
あと、Bloggerの方に、斜方投射のアニメーション動画をアップしておいたので、興味のある奴は見るといいにゃ。
リンク先はこちら↓
http://nemneko.blogspot.com/2019/08/blog-post.html
そしていつもどおり、自画自賛をして、記事を結ぶのであった。
斜方投射と到達距離 [ねこ騙し物理]
斜方投射と到達距離
北朝鮮が、また、ロケット花火を打ち上げたので、極初歩的な弾道軌道の計算をしてみることにするにゃ。
地球は球形だけれど、問題を簡単化するために、地球は平ら、重力加速度gは一定と仮定する。
§１ 放物運動 （空気抵抗がない場合）
時刻t=0で速さ
で地平線（x軸）となす角度θで、質量mの質点が投げ出されたとする。
このとき、（ニュートンの）運動方程式は次のようになる。
初期条件は
とおき、①の両辺をmで割ると、
t=0のとき
だから
これを積分すると
初期条件x(0)=0から
②の両辺をmで割ると
とおき、これを書き換えると
重力加速度gは定数だから、両辺を積分すると
初期条件
より
さらに積分すると
初期条件y(0)=0より
したがって、
が解になる（高校の物理の公式）。
次に、質点の運動の軌道を求めることにする。
（３）より、
これを（４）式に代入すると
したがって、最高の高さHは
到達距離Lは
0<θ<９０°とすると、sin2θ=1になるのはθ=45°の時だから、最大到達距離
である。
もっと簡単に解けるけれど、高校の物理などとの兼ね合いで、こう解いてみたにゃ。
問 高校の物理らしく、最高の高さHと到達距離を求めよ。
【解】
v=0のときに、yは最大になる。
したがって、（２）より
これを（４）式に代入すると、
（４）式でy=0になるtの値を求めると、
したがって、
（解答終）
§２ 空気抵抗がある場合
以上の議論は空気抵抗がない場合の話。空気抵抗が速度に比例する場合、運動方程式は次のようになる。
kは比例定数だケロ。
初期条件は空気抵抗がない場合と同一とする。
③の一般解は見た瞬間
とわかる。
初期条件はだから
これを積分すると
t=0のときx=0だから
④はじっと見ると
初期条件はだから
これを積分すると
t=0のときy=0だから
よって、
したがって、微分方程式の解は
（１１）をtについて解き、それを（１２）に代入すれば、軌道をy=f(x)という形で表すことができるけれど、形がオドロオドロなりすぎるにゃ。
問 変数分離で④を解け。
地平線（x軸）にぶつかれば跳ね返るけれど、跳ね返ることなくもし永遠に落下し続けるとすれば、終端速度は
そして、
ここで、お前らに問題。
問題１ t≧0とするとき、（１２）式で与え与えられるyの最大値Hを求めよ。
【ヒント】
yが最大になるとき、
到達距離Lを求めようと無謀なことは考えないほうがいいにゃ。
だって、
の解を求める必要があるからだにゃ。
まぁ、t=0という自明な解はすぐに求まるが・・・。
問題２ km/t>0が１に対して非常に小さいとき、空気抵抗がない時の解（１）〜（４）は、空気抵抗を考慮した解（９）〜（１２）の近似になるはずである。
このとことを示せ。
【ヒント】
ところで、ネムネコは微分方程式
を次元（物理単位）をもつまんま解いたけれど、この微分方程式は何らかの方法で無次元化して解くべきだね。
たとえば、無次元化された速度
は無次元化された時間
を導入する。
すると、（無次元化された）微分方程式は次のようになるのかな。
ここで、γは無次元化された重力加速度で
こうすれば、「メートル」、「フィート」や、「キログラム」、「ポンド」といった単位によらないより普遍的な方程式になり、微分方程式を解くのも楽になる。
こうすると、初期条件は
ネムネコが考えるに、
ddt³さんが、きっと、こういった話をしてくれるに違いない。
さらに、初期値問題⑨を解くスプレッドシートを作ってくれるに違いない(^^ゞ
定積分の近似と誤差の限界 [微分積分]
定積分の近似と誤差の限界
f(x)を[a,b]でC¹級、すなわち、f'(x)が[a,b]で連続であるとする。
を
と近似したときの誤差について考える。
部分積分すると
仮定より、[a,b]において(b−x)は非負連続なので、積分の（第一）平均値の定理より
を満たすξが存在する。
したがって、
これが、積分の値を（１）式の右辺で近似したときの誤差を与える式になる。
次に、を
と近似した時の誤差について考える。
部分積分すると
積分の平均値の定理より、
したがって、
f(x)=x、a=0、b=1の場合について計算してみると、
また、f'(x)=1なので、
となり、（２）、（４）を満たしていることがわかる。
定理
f(x)が[a,b]でC¹級であるとき、
を満たすξが存在する。
また、
とおき、（２）、（４）の剰余項が消えるようにLとRの平均を取ると、
と台形公式が得られる。
区間[a,b]をn等分し、
となるように点
を配する。
すなわち、
そして、[a,b]の小区間
を作り、
とすると、定理より
が成立する。
また、
という関係があるので、
の最大値をMとすると、
したがって、
と、を近似した誤差の限界は
になる。
だから、分割数を１０倍にすれば誤差は約1/10に、分割数を１００倍にすれば約1/100になる、
台形公式の場合、
とおくと、
なので、
同様の議論をすると、
とを近似したときの誤差の限界は
ここで、Mは
の最大値である。
つまり、台形公式（６）の場合、分割数を１０倍にすると誤差は約1/100、分割数を１００倍にすると誤差は約1/10000になるってワケ。
台形公式と誤差 [微分積分]
台形公式とその誤差
有界閉区間[a,b]で定義された関数f(x)の積分
を
で近似する計算法を台形公式という。
f(x)が１次以下の多項式（関数）である場合、
が成立する。
f(x)=x²として、
の近似値を（１）で求めてみると、
になる。
したがって、この時の誤差は
問題１ a、bが任意の実数のとき、
はf(x)が１次以下の多項式であれば成り立つが、f(x)が２次以上の多項式であると、一般に、成り立たない。
このとき、次の問に答えよ。
（１） f(x)=x²のとき、（A)の右辺の値を左辺の積分の近似値とみれば、その誤差は(b−a)³に比例することを示せ。
（２） 積分
を求めるのに、積分の区間0≦x≦1をn個の区間に等分し、各区間ごとに近似式（A)を用い、それらの近似値を加えて出す。このようにして得られる積分の誤差を1/10⁴より小さくするには、nをいくつにすればよいか。
【解】
（１）
だから、
（２） [0,1]をn等分し
とおき、
とすると、（１）より
したがって、誤差は
になる。
この誤差を1/１０⁴より小さくすれば良いので、
したがって、n≧41にとればよい。
（解答終）
本当にこの答えが正しいかどうか、自分が過去に作ったプログラム（スクリプト）でn=40、41について計算してみると、n=40のときは誤差が1/10⁴を越え、n=41のとき誤差が1/10⁴より小さくなっていて、n=41が求める答えであることがわかる。
定理１ （台形公式の誤差）
f(x)は[a,b]でC²級とすると、
を満たすξが存在する。
【証明】
とおくと、部分積分を２度用いると
[a,b]で(b−x)(x−a)は非負連続なので、積分の平均値の定理より
を満たすξが存在する。
一方、
だから、これを（B）に代入すると、
（証明終）
定理 積分の（第一）平均値の定理
f(x)が[a,b]で連続、g(x)が[a,b]で非負連続ならば、
を満足するξが存在する。
問題２ f(x)を[a,b]でC²級の関数とするとき、
となることを示せ。
ただし、
とする。
【解】
部分積分より
一方、
となるので、これを代入すると、
（解答終）
さてさて、f(x)=x²、a=0、b=1とすると、
となるので、定理から
となり、先に求めた結果と一致している。
また、この定理から、問題１の（１）の答が
であることがわかる。
ところで、g(x)=1とすると、積分の平均値の定理から
となるξが存在する。
このξの代わりに
と、aとbの中点をとって、を次のように近似する。
これを中点公式という。
そして、中点公式も台形公式と同じく、f(x)が１次以下の多項式（関数）であれば、
と等号が成立し、（２）の右辺でを近似した誤差は(b−a)³に比例する。
ここで問題になるのは、台形公式と中点公式のどちらが精度よく計算できるかだ。
試しに、を中点公式で近似すると、
したがって、誤差は
となり、台形公式の誤差の1/2倍になっていることがわかる。
定理２ （中点公式の誤差）
f(x)は[a,b]でC²級とすると、
というわけで、一般に、中点公式の方が台形公式よりも精度よくを近似することができる。
また、
とすると、シンプソンの公式、台形公式、中点公式の間には
という関係が成り立ち、シンプソンの公式は、台形公式と中点公式の剰余校が消えるような重みつき平均になっていることがわかる（※）。
なお、図形的に言うと、台形公式のTの値は下図の台形ABGF、中点公式Mの値はEDGF、そして、シンプソンの公式のSの値は、y=f(x)を２次関数で近似した放物線（青色で示してある）とx軸、直線AF、BCで囲まれた面積を表す。
（※） 定理１、定理２の式中に出てくるa<ξ<bの値は、一般的に異なることに注意。
お前らに質問 （８月２２日 定積分の近似 台形公式） [お前らに質問]
今回は、簡単な問題！！
問題 a、bが任意の実数のとき、
はf(x)が１次以下の多項式であれば成り立つが、f(x)が２次以上の多項式であると、一般に、成り立たない。
このとき、次の問に答えよ。
（１） f(x)=x²のとき、（A)の右辺の値を左辺の積分の近似値とみれば、その誤差は(b−a)³に比例することを示せ。
（２） 積分を求めるのに、積分の区間0≦x≦1をn個の区間に等分し、各区間ごとに近似式（A)を用い、それらの近似値を加えて出す。このようにして得られる積分の誤差を1/10⁴より小さくするには、nをいくつにすればよいか。
簡単な問題のはずだから、さくっと解いてもらおうじゃないか。
このブログの数学の記事を大学受験生が見ているとは思わないが、ひょっとしたら、どこかの大学で似たような問題が出るかもしれないにゃ。
だから、（２）の答は４１じゃ〜なかろうか。（このスクリプトによると、n=40だと、1/10⁴より大きい・・・）
いやまぁ、これを計算するプログラムを自分で作り、そして、計算してこの答（？）を見つけ出したというのであれば、解答として認めてもいいが・・・。
問題の解答とさらなる質問 （定積分の近似計算） [お前らに質問]
問題２ a、bを任意の実数とするとき、等式
はf(x)が３次以下の多項式であれば成り立つが、４次以上の多項式であると成り立たない（a≠b）。
f(x)=x⁴のとき、(A)の右辺の値を左辺の積分の近似値とみれば、その誤差は(b−a)⁵に比例することを示せ。
【略解】
よって、誤差は(b−a)⁵に比例する。
（略解終）
（A）で定積分を近似した場合、
その誤差は(b−a)⁵に比例する
という結論は（数値計算の）理論上重要な事実。
この簡単な計算から、シンプソンの公式の（局所的な）誤差が(b−a)⁵に比例するということがわかるんだよね。
昨日の数学の記事で、
「a、bを任意の実数とするとき、等式
はf(x)が３次以下の多項式であれば成り立つ」
と書いた。
そして、公式（A）の右辺はf(x)を
の３点を用いて２次関数で近似することによって得られるという話をした。
だとすれば、（A）式で正確に計算できるのはf(x)が２次以下の多項式（関数）の場合であるはずだ。
しかし、
f(x)=x³、さらに、a=0、b=1とし、（A）式の右辺を用いて、
を近似的に（？）計算すると、
となり、
の正しい値を与える。
考えれみれば、妙な話だ。
というわけで、お前らに次の命題を証明してもらおうじゃないか。
命題
a、bを任意の実数とするとき、等式
はf(x)が３次以下の多項式であれば成り立つ
もちろん、
とし、
と、
を計算し、この両者の値が等しくなることで、証明してもらって結構だにゃ。
さらに、
【Β（べーた）関数を用いた証明】
と置き、置換積分をすると、
だから、
である。
したがって、
（証明終）
もちろん
だから、
