定積分の矩形公式と台形公式の誤差 [数値解析]
定積分の矩形公式と台形公式の誤差
[a,b]をN等分した次の点列があるとする。
分割幅をhとすると、
となるが、
このとき、
と積分
を近似したときの誤差について考える。
その前に準備として、f(x)がC¹級であるとき、
となるξがa<ξ<bに存在することを以下に示す。
とし、左辺を部分積分すると、
になる。
b−aは[a,b]で非負なので、積分の第一平均値の定理より
となるξがa<ξ<bに存在する。
したがって、
となる。
よって、
ここで、
とすれば、
そして、これが、
と近似したときの誤差の評価式になる。
本当にそうなのか、
とし、
の場合について確かめてみる。
上のグラフは横軸に分割数N、縦軸に誤差の絶対値をとり、対数グラフで表したもので、このグラフを見ると、直線の勾配が約１であり、分割数Nを１０倍、つまり、分割幅hを1/10にする誤差が約1/10になることがわかる。
このことより、この誤差の評価式の妥当性を確かめることができる。
また、
となるので、同様の議論から、
と、定積分を近似したときの誤差の評価式が
で与えられることがわかる。
問題１ f(x)を[a,b]でC²級の関数とする。このとき、
となるξが存在することを示せ。
【解】
部分積分を２回用いることにより、
(b−x)(x−a)は[a,b]で非負なので、積分の平均値の定理より
となるξがa<ξ<bに存在する。
また、
だから、
（解答終）
そして、
と台形則で定積分を近似したときの誤差が
となることがわかる。
問題２ 次の等式が成立することを確かめよ。
問題３（積分の第１平均値の定理）
fが閉区間[a,b]で連続、gが[a,b]で非負連続ならば、
を満たすξが存在することを示せ。
テーラー展開を使って微分本方程式の初期値問題を解く２ [数値解析]
テーラー展開を使って微分本方程式の初期値問題を解く２
f(x)を
級の関数とする。このとき、f(x+h)はxで次のようにテーラー展開が可能である。
特にn=1のとき、
で、これは平均値の定理と呼ばれる。
これから、
となり、f(x+h)を
と近似したときの誤差が
のオーダー、すなわち、
であることがわかる。
また、f(x+h)をxのまわりで１次の項までテーラー展開すると、
となる。
したがって、h≡0のとき、
となり、
と近似したときの誤差がO(h)であることがわかる。
次に、
という常微分方程式の初期値問題について考えよう。
最もシンプルな方法は、
を用い、これを左辺におき、
と近似するものであろう。上述の議論から、
と近似したときの（打ち切り）誤差がh²程度であることがわかる。
と、点列
を等間隔hに配置すれば、次の漸化式を得ることができる。
これがオイラー法と呼ばれるもので、オイラー法の局所的な打ち切り誤差はO(h²)、すなわち、h²程度である。
次の常微分方程式の初期値問題
を、h=0.1としてオイラー法で解いてみることにする。
このとき、オイラー法は
となり、
だから、
となる。
初期条件は
となるので、
を用いて、
の値を定めることができる。
表計算ソフトを用いてx=2まで計算した結果は次のとおり。
オイラー法を用いて計算を進めるにつれ、厳密解
との誤差が増大してゆくことがわかる。
そこで、テーラー展開を利用して、より精度よく計算する手法を考えることにする。
であるから、合成関数の微分法よりy''は
となる。
一方、y(x+h)のxにおけるテーラー展開を２次の項までとると、
したがって、
これを微分方程式の左辺に代入すると、
となる。
そして、これから、次の漸化式を得ることができる。
この漸化式の局所打ち切り誤差はO(h³)なので、オイラー法よりも精度よく計算できるはずである。
この漸化式を用いて、常微分方程式の初期値問題
の近似解を求めることにする。
このとき、
となるので、
したがって、
h=-0.1として、表計算ソフトを用いて解いた結果は次の通りである。
この表とグラフを見ると――テーラー展開による方法の数値解と厳密解はほとんど一致しているので、このグラフでは同一の曲線に見える――、劇的に計算精度が向上していることがわかる。
この問題の場合、高次微分が
と求められるので、より高精度の
漸化式を得ることができる。
ちなみに、この漸化式の局所打ち切り誤差は
。
h⁵の項までとれば、４次のルンゲ・クッタ方よりも高精度に計算することが可能である。
ネムネコ法、２次のルンゲクッタ法に勝つ！！（差分の基礎） [数値解析]
ネムネコ法、２次のルンゲクッタ法に勝つ！！（差分の基礎）
ネムネコが新たに提案した、差分を用いた１階常微分方程式の初期値問題に対する数値解法は、実際、どれくらい使い物になるのか、次の非線形常微分方程式の初期値問題に使ってみた。
ネムネコが提案した方法は、上の微分方程式を
と差分方程式に置き換えて解く方法。
この方法をネムネコ法と、ひとまず、仮称することにする。
（２）式の右辺第１項の
は、前進差分法（陽解法）
に由来するものであり、第２項の
は後退差分（陰解法）
に由来するもので、（２）式の打ち切り誤差はO(h³）で２次のルンゲ・クッタ法や修正オイラー法の誤差と同程度と考えられる。
ネムネコ法は、陽解法である修正オイラー法や２次ルンゲ・クッタ法とは異なり、
の値を求めるためには超越方程式（２）を解く必要がある陰解法である。
（２）を（数値的に）解く必要があるので、修正オイラー法や２次のルンゲ・クッタ法のように表計算ソフトなどを使って簡単に微分方程式（１）の数値的な近似解を求められないという欠点がある。
なのですが、自分でプログラムを作ってこの問題を解くという場合には、これは必ずしも欠点にならない。というのは、（２）はニュートン法などを用いることによって簡単に数値的に解くことができるからだ。
そこで、プログラムを作り、解いてみた。
そして、Δx=0.1としてx=0〜1まで計算してみた結果は次のとおり。
厳密解とよく合っているじゃないのよ。我ながら驚きの結果！！
グラフにしてしまえば、厳密解との差がわからないほどよく合っているからね〜。
しかも、この計算結果は、Δx=0.1として２次のルンゲ・クッタ法を用いて解いた計算結果よりも精度が高いんだケロよ。
CASIOさんの高精度計算サイトの計算結果はコチラ。
スゴイにゃ、（仮称）ネムネコ法！！
そして、どうも、（仮称）ネムネコ法は、線形常微分方程式よりも非線形常微分方程式で強みを発揮するみたいだね。
この計算に使用したプログラムは次のとおり。
! ネムネコ法(^^ゞ
parameter (n=10)
real x(0:n) , y(0:n)
eps = 1.e-6
x = 0. ; y=0. ! 初期化
x(0)=0.; y(0)=0. ! 初期条件
dx =0.1 ! 分割の幅
do i = 1, n
x(i)= i*dx
yo = y(i-1);
! ニュートン法で非線形方程式解く
do k=1, 10
! yn = yo - (yo - dx*sin(yo)-y(i-1)-dx)/(1-dx*cos(yo))
yn= yo- (yo-0.5*dx*(sin(yo)+sin(y(i-1)))-y(i-1)-dx)/(1-0.5*dx*cos(yo))
err = abs(yn-yo)
if (abs(yn-yo).lt.eps*abs(yo)) exit ! 収束判定
yo=yn
end do
y(i)= yn ! 計算結果をセット
end do
! 計算結果の出力
write(*,*) ' x y(calc) y(exact)'
do i=0, n
write(*,100) x(i), y(i) , 2.*atan(x(i)/(2-x(i)))
end do
100 format(f10.7,1x,f10.7,1x,f10.7)
end
ちなみに、（１）の厳密解は、
などと求めることができる。
なお、
画像元：YouTubeの下の動画
お前らに質問（前進差分と後退差分）の答もどき [数値解析]
お前らに質問（前進差分と後退差分）の答もどき
常微分方程式の初期値問題
の解をy=y(x)で表すことにする。
この解は
ですが・・・。
右の図を見て欲しいのですが、前進差分による微分方程式（１）の解法というのは、（１）の解曲線y=y(x)の点(x₀,y(x₀))における接線を使ってx=x₁におけるy(x₁)の値を推測するもの。
右の図では赤い色のy₁がその推測値。
対して、後退差分による微分方程式（１）の解法は、点(x₁,y(x₁))における解曲線y=y(x)の接線の傾きをもった点(x₀,y(x₀))を通る直線をもとに、x=x₁の値y(x₁)を推測するもの。右の図では青いy₁で示されている。
というわけで、解曲線y=y(x)が下に凸であれば――下に凸ならば必ず接線は曲線の下側にある――、右の図に示してあるように、前進差分を使って推測したy₁は真の値y(x₁)より小さくなる。そして、後退差分を用いて推測したy₁は真の値y(x₁）よりも大きくなるというわけ。
逆に、解曲線y=y(x)が上に凸であれば――下に凸ならば必ず接線は曲線の上側にある――、前進差分を使って推測したy₁は真の値y(x₁)より大きくなり。そして、後退差分を用いて推測したy₁は真の値y(x₁）よりも小さくなるんだね〜。
というわけで、開曲線の（局所的な）凸凹が大小関係を決定するというわけなんですね〜。
二階微分y''が存在すれば、曲線の凸凹はy''(x)の正負の符号で判定できるから、y''の符号を調べれば大小関係は判定できる！！
なお、この話は、計算領域で解曲線y=y(x)の凹凸が変わらないことを前提としたお話。計算領域で凹凸が変わったら通用しない議論なのでその点は注意するにゃ。
このことは、
前進差分
の漸化式からも確かめられると思うにゃ。
これは点(x₀,y₀)における接線の方程式に他ならないんだから。
一方、後退差分の漸化式
は、(x₁,y(x₁))におけるy=y(x)の接線の傾きと等しい傾きを持つ点(x₀,y₀)を通る直線の方程式を表しているにゃ。
２階微分y''は曲線y=y(x)の曲がり具合、曲率、接線（直線）とのズレ具合、離れ具合を重要なあらわす因子だから、とっても重要なものなんだケロよ。
なお、x₂、y₂以降は、次の図のように順次計算してゆく。
なお、上の図はあくまで計算の仕組みを説明するための図で、実際の計算点の間隔はこの図のように大きくないので、この図のように誤差は大きくないケロ。
この点だけはくれぐれも誤解なきように！！
非線形の微分方程式だって解けるケロ（差分の基礎編） [数値解析]
非線形の微分方程式だって解けるケロ
差分法を用いた、次の非線形微分方程式の初期値問題について考える。
前進差分を用いて（１）を差分方程式に近似すると次のようになる。
これから、
という漸化式を得ることができ、x₀=0におけるyの値y₀を計算の出発点にし、前進的にの値を次々に求めることができる。
これに対して、後退差分を用いると、（１）の微分方程式に関する差分方程式は
となり、次の２次方程式を得ることができる。
この２次方程式を解くと、
という漸化式が得られ、（３）式を用いて前進的にの値を次々に求めることができる。
（１）の場合、後退差分は２次方程式を解く必要があるんだケロよ。２次方程式を解くというひと手間が加わるだけではなく、（２）と（３）を比較すればわかるように、漸化式が複雑になる。
だったら、３次方程式を解かないといけない。３次方程式ならば解の公式があるので代数的に解けないことはないけれど、n≧5のとき
この式から得られる
というのn（≧５）次方程式を代数的に解くことは一般にできないので、ニュートン法などを用いてを解く必要がある。
ということで、後退差分（陰解法？）は前進差分（陽解法？）とは違って難しいんだケロよ。だから、常微分方程式の初期値問題で後退差分を用いた解法（陰解法？）は一般に使われない。
陰解法が威力を発揮するのは、
といった偏微分方程式の数値解法なんだケロよ。
ではありますが、Δx=0.1とし、（２）式と（３）式を使って、（１）の初期値問題を解いた結果は次のとおり。
微分方程式（１）の（厳密）解は簡単に求まりそうに見えるだろうが、これはリッカチ形の微分方程式で、この解は三角関数、指数関数、対数関数を用いて表すことができない、つまり、解けないんだケロ。
なので、比較参考のために、４次のルンゲ・クッタ法を用いてΔx=0.1で計算した値を厳密解のかわりにあわせて図示している。
この計算結果を見ると、陽解法による計算結果は厳密解よりも小さく、陰解法による数値計算結果は厳密解よりも大きいこと、そして、厳密解がこの間に位置していることがわかる。
ということで、陽解法による計算結果を陰解法による数値計算結果を足して２で割る、つまり、両者の平均をとると、厳密解とかなりよく一致することがわかる。
このことは、数値計算結果の絶対誤差を図示したグラフを見るとよくわかると思う。
前回に引き続き、陽解法と陰解法の計算結果を足して２で割ると、誤差のプラスマイナスが互いに消し合って、誤差が小さくなって、厳密解に近い値が出るようになるのであった。
Δx=0.1のときに、誤差がおよそ1/10になっており、また、陽解法、陰解法の誤差の程度がO((Δx)²)であるので、この解法の誤差はO((Δx)³)と推測される。だから、修正オイラー法、２次のルンゲ・クッタ法と誤差のオーダーと同じであるに違いない。
より進んだ解法は、
①と②の両辺を足しあわせて、両辺を２で割ると、
となるにゃ。
についての２次方程式（４）を解き、こうして新たに得られた漸化式を用いて計算してみるといい。
（４）を使った計算結果は、（２）と（３）を使って求めた値の平均値と同じ値になると思うかもしれないけれど、これが一致しないんだね〜。
（４）を使ったもののほうが精度はさらに向上するはずなんだよね〜。
嘘だと思うならば、の２次方程式（４）を解き、の漸化式を求め、その漸化式で計算してみな。
というか、お前ら、これくらいはやれよな。
何故、（４）式は、（２）や（３）より精度よく計算できるのかについては、尤もらしい説明をすることはできなくはない。
というのは、
と
の中点を
と表すことにし、それに対応するyの値を
で表し、
③と④は、それぞれ、
となる。
⑤と⑥を足すと、
が消えて、（４）式になる。
③は、Δxより小さいΔx/2という分割幅で仮想のの値を求めるいることを意味し、④式はこうして求められたを元にΔxよりも小さい分割幅でを求めることを意味している。
分割幅を小さくすると一般に計算精度が向上する。
また、仮想の点
とそれに対応するyの値を考慮に入れると、（４）の左辺の
はにおける中心差分とみなすこともできる。
したがって、（４）は、（２）や（３）よりも高精度に計算できる・・・。
感覚的にこう理解すればいいのではないか。
これにテーラー展開の話を交え、「（２）や（３）の誤差の程度がO(h²)であるのに対し（４）式はO(h³)で修正オイラー法、２次のルンゲ・クッタと誤差の程度が同オーダーである」なんて尤もらしい説明もできるけれど、かなり胡散臭い話だから、ここではしないにゃ。
（２）、（３）、（４）式は、０≦φ≦１である重みφを使うと、
と、１つの式で表すことができる。
φ=0のときは陽解法、φ=1/2のときは半陰解法、φ=1のときは純陰解法と分類分けすることも可能。
そして、プログラムを作ってこの問題を解く場合、0≦φ≦1のすべてのφについて解けるものを作った方がいいにゃ。それが賢い方法だにゃ。
お前らに質問（前進差分と後退差分） [数値解析]
お前らに質問（前進差分と後退差分）！！
次の微分方程式の初期値問題がある。
この微分方程式の解が
であることは言うまでもない。
ところで、この微分方程式を差分方程式で近似すると、
となる。
一方、後退差分を用いて（１）を差分方程式に書き換えると、
Δx=0.1とし、（３）と（５）を用いて計算した結果は次のとおり。
計算結果の大小関係は
後退差分≧厳密解≧前進差分
になる。
次に、
について考える。
この解が
であることは明らかだろう。
前進差分の場合、
後退差分の場合、
h=0.1とし、（７）、（８）を用いて計算すると、
となり、大小関係が逆転し、
前進差分≧厳密解≧後退差分
の順になる。
これは偶然ですかい？
それとも、この大小関係の逆転、そして、厳密解が前進差分と後退差分による数値解の間に来ることには、何か、深い数学的な理由があるのですか？
お前らに、この深遠な（？）問題について答えてもらおうじゃないか！！
さらに新しい解法を提示 [数値解析]
さらに新しい解法を提示
§１ 平均してみた
次の１階常微分方程式の初期値問題がある。
前進差分（陽解法？）を用いてこの微分方程式を差分方程式に書き換えると、
これから、次の漸化式が得られる。
初期条件からx₀=0のときのy=y₀=0を計算の起点に、（２）を用いてy₁、y₂、・・・と前進的にx=x₁、x₂、・・・におけるyの値を求めることができる。
対して、後退差分を用いる場合（陰解法？）は、微分方程式を差分方程式
に書き換え、次の漸化式が得られる。
前進差分の場合と同様に微分方程式の数値解を求めることができる。
Δx=0.1について解いたものは次のとおり。
前進差分を用いた数値解は厳密解より大きく、後退差分の場合は厳密解の値よりも小さくなることがわかる。
絶対誤差も同程度なので、前進差分と後退差分を用いた結果を足して２で割れば、厳密解により近くなのではないかと予想できる。
予想通り、よく一致しているにゃ。
計算結果をグラフで表すと、この微分方程式の厳密解
と完全に一致しているように見える(^^)
数値計算って不思議で奥が深い、面白いと思わないかい？
§２ 少しだけ数学的な話
微分方程式（４）を前進差分を用いて差分方程式で近似すると、
これから、
という漸化式が得られる。
y₀=1なので、この漸化式から与えられる一般項は
になる。
[0,1]をn等分したものを
とおき、i=nとすれば、
n→∞のときの極限値を求めると、
となり、厳密解のy(1)=1/eと一致する。
このことは、コンピュータを用いて数値的に計算する場合、コンピュータによる計算に特有な丸め誤差などがあるのでそのとおりにならないけれど、分割数を大きくすればするほど、すなわち、分割の幅を小さくすればするほど、数値解と厳密解の誤差は小さくなり、終局的に、この両者は完全に一致する、ということを表している。数値解析の分野では、このことを適合性（consitency）があると表現する。
適合性または整合性 (consistency)
空間および時間を離散化した時の格子幅を限りなく0に近づけたときに、離散化方程式と元の微分方程式の差が0に収束することである。この差は一般には格子点についてのテイラー展開によって評価される。
https://goo.gl/2FFiDU
大きく脱線してしまった。ネムネコはこんなことを話したいわけではない。
漸化式（５）からの一般項は
になる。
したがって、
になる。
だ・か・ら、数値解の値が正しいかどうかは別にして、前進差分の数値解が意味を持つのは0<Δx<2の場合なんだケロ。
0<Δx<2という条件を満たさないと、この計算法は安定じゃない。
対して、後退差分を用いた解法は、
で、
が成立するので、Δxの値にかからず、
となり、数値的に安定している。
後退差分を用いたこの解法は無条件安定なんだケロ。
常微分方程式の初期値問題によく使われるルンゲ・クッタ法も前進差分を用いているので、Δxをあまり大きく取り過ぎると数値的に不安定になり、数値解が激しく振動したりするんだにゃ。知っていたケロか。
非線形でもないかぎり、常微分方程式の初期値問題の数値解法で振動解に見舞われることはないと思うけれど・・・。
計算に使用したスプレッドシートは以下のところに公開
エクセル版
ウェブ版
課題の解答例 [数値解析]
課題の解答例
課題 中心差分を用いて次の微分方程式の初期値問題を解く方法を考えなさい。
その方法を用い、Δx=0.1としてx=1まで計算し、厳密解と数値解とを比較しなさい。
【解答例】
の両辺を微分すると、
①式を用いて②式から
を消去すると、
①式がx=0においても成立するとすると、
したがって、この初期値問題は
となり、
差分法を用いた微分方程式の初期値問題の解法の計算結果を流用できる＼(^o^)／
（解答例終）
一階の微分方程式のままだと中心差分はうまく使えないので、中心差分を使えるように、この微分方程式を微分し２階の微分方程式にすればいいと言うわけなんですよ(^^)
修正オイラー法の計算結果は右図のとおり。
修正オイラー方を用いたほうが精度よく計算できる。
修正オイラー法の誤差の挙動が面白いね。
x=1.8の付近で最小値をとったあと、誤差が上昇しているからね。
何故、こうなるかについてはわからない(^^ゞ
これで終わってもいいのですが、クランク・ニコルソン法をヒントに、差分法を用いた一階常微分方程式の初期値問題の新しい（？）解法を思いついたので紹介するにゃ。
微分方程式
を前進差分を用いて近似すると、
後退差分を用いて近似すると、
（１）と（２）を足すと、
これをについて解くと、
という漸化式が得られる。
x₀=0におけるyの値y₀=0を計算の起点にし、漸化式（４）を用いて前進的に微分方程式（１）の数値解を求めることができる。
なお、
ね。
この解法の誤差の程度はたぶんO(h³）で、修正オイラー法や２次のルンゲ・法と同程度の誤差をもっているものと推定される。
（４）式を用いて計算した結果は次のとおり。
誤差の点で修正オイラー法に劣るものの、グラフにしてしまえば厳密解と新しい数値解は重なっており、良好な結果が得られていることがわかる。
また、２階常微分方程式直して解いたものの最大誤差が0.0047746874であり、新しい解法のそれは0.0013790814と精度的にも向上していることがわかる。
数値解、厳密解と、その誤差の挙動が一致しているあたり、なかなかオシャレだと思うにゃ。素性がいいね、この解法は。
この記事の計算で使用したスプレッドシートはコチラ。
オイラー法、修正オイラー法との比較
エクセル版
Web版
新しい解法（？）
エクセル版
Web版
差分法を用いた微分方程式の初期値問題の解法（基礎編） [数値解析]
差分法を用いた微分方程式の初期値問題の解法（基礎編）
次の２階線形微分方程式の初期値問題の差分法を用いた解法について考える。
p=dy/dxとおけば、上の微分方程式は
となり、オイラー法やルンゲ・クッタ法を用いて解くことができる。
オイラー法、ルンゲ・クッタも差分（正確にはテーラー展開）を元にして導出されるので、これらも差分法の一種と考えることができるが、ここでは異なる解法について考えることにする。
計算の起点となる点aをx₀、そして、x₋₁=a−Δx、x₁=a+Δxとし、それに対応するyの値をそれぞれy₀、y₋₁、y₁とすると、x=0における微分は中心差分を用いて、
と近似することができる。
したがって、初期条件より、
となる。
これを２回微分の式に代入すると、
となり、これを微分方程式（１）に代入すると、
となり、この式を用いて未知数であるy₁を求めることができる。
におけるyの推測値については、微分方程式（１）の差分方程式
から、
となり、この漸化式を用いて前進的に解くことができる。
新たに得た漸化式を用いて、次の初期値問題を解くことにする。
このとき、漸化式は
となる。
これならば、表計算ソフトを使って数値解を求めることができる。
Δx=0.1とし、x=3まで計算した結果は次のとおり。
この解法の局所的な打ち切り誤差はO(h³)であり、誤差の程度は修正オイラー法や２次のルンゲ・クッタ法とほぼ同程度と考えられる。
差分法を知っていれば、（修正）オイラー法や２次のルンゲ・クッタ法を知らなくても、線形の微分方程式の初期値問題も解くことができるのであった。
なお、上の初期値問題の厳密解は
である。
この計算に使用したスプレッドシートはコチラ↓
計算結果だけを知りたいヒトはWeb版のコチラ↓を
差分法と差分法を用いた微分方程式の解法（基礎編） [数値解析]
差分法と差分法を用いた微分方程式の解法（基礎編）
差分法でナビエ・ストークス方程式（略してNS方程式という）を解くことになりましたというコメントを、今日、いただいたので、差分法を用いた微分方程式の解き方について簡単に説明することにするにゃ。
何回でも微分可能な関数f(x)があるとする。
すると、h>0として、f(a+h)をテーラー展開すると、
といった式が得られる。
（１）より、
となるので、hが１に対して十分小さいとき、
と近似することができる。（３）式による近似は、hのオーダーの誤差を含んでいるので、このことをO(h)と表して、
と書いたりする。そして、これを前進差分という。
同様に、（２）式からは
というf'(a)の近似式が得られる。これを後退差分という。
さらに、（１）と（２）の差をとると、うまい具合にhの偶数次数の項が消えて、
で、hの３次以上の項をすべて寄せ集めたものをO(h³)と表わせば、
これを中心差分という。
前進差分、中心差分のどちらが精度よく計算できるか、h=0.1とし、f(x)=sinxの場合について見てみることにする。
f(x)=sinxの導関数f'(x)=cosxの中心差分による近似は極めて良好であるのに対し、前進差分の場合は誤差が大きいことがわかるだろう。
前進差分は(a,f(a))と(a+h,f(a+h))の２点を使ってその傾きを、中心差分は(a−h,f(a−h))と(a+h,f(a+h))の２点を使ってその傾きを求め、それをf'(a)に近似しているのだけれど、中心差分と前進差分による近似には月とスッポンくらいの開きがあるんだケロよ。
分割の幅hを変化させて、x=0.5における前進差分と中心差分の誤差を求めたものは次のとおり。
前進差分はhを1/10にすると誤差は約1/10になるのに対し、中心差分は誤差が約1/10²になるんだケロ。
これが記号O(h)とO(h²)の意味するところ。
つまり、
という謎の記号は、hを1/10にしたら、誤差は
になるということを表しているんだケロよ。
２階の微分係数f''(a)は、（１）式と（２）を足すと、うまい具合にhの奇数次の項が消えて、
になる。
h⁴以上の高次の項を寄せ集めたものをO(h⁴)で表すと、
となる。
そして、これは中心差分を用いた２階の微分係数f''(a)の近似式ということになる。
誤差のオーダーはO(h²)だから、hの２乗オーダーと見積もることができる。
でだ、次のような２階の微分方程式の境界値問題があるとする。
計算領域[0,1]を
と等分割h=1/nに分割し配置したとすると、この微分方程式に含まれる微分は、中心差分を用いて、
と近似できる。
で表すことにし、これを微分方程式に代入すると、
というn−1元の連立方程式を得ることができる。
で既知で、未知数は
のn−1個に対して、連立方程式の数はn−1だから、
連立方程式②を解くことによって、微分方程式①の近似解を求めることができる。
これが差分法を用いた微分方程式の境界値問題の最も基本的な解法の例である。
今日、コメント、質問をくださった翼さんは、NS方程式という非線形の連立偏微分方程式を差分法を用いて解かなければならなくなったので、これよりはずっとずっと複雑ですが、実は、この２階線形常微分方程式の解法が基礎になっている。
NS方程式のような２階の偏微分方程式の境界値問題であろうが、２階の常微分方程式の境界値問題であろうが基本は同じ。
翼さんが何をなさろうとしているのか、私にはわからないので何とも言えませんが、たぶん、研究室にはNS方程式を解く商用の汎用プログラムなどがあると思うんだよね〜。
そうではなく、１からプログラムを作らなければならないとしたら、最低でも１年はプログラム作りに専念しなければならないと思うにゃ。
NS方程式というのは、
たとえば、
といった形の非線形の連立微分方程式。
実は、これだけではNS方程式は解けなくて、連続の式（質量保存則）
を加え、偏微分方程式を閉じさせて解くんだよね〜。
NS方程式は非線形性が強くてただでも解きにくいのに加え、圧力pを解くのが非常に難しい。このすぐれた解法を開発するだけで数値解析のアカデミックな論文になり、長く歴史に名を残すことになる。１９６０年代から１９７０年代に開発された化石級の計算スキームが現在でも使われているくらいだから（笑）。
