差分法と差分法を用いた微分方程式の解法（基礎編）

差分法と差分法を用いた微分方程式の解法（基礎編）

 

 

差分法でナビエ・ストークス方程式（略して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を解くのが非常に難しい。このすぐれた解法を開発するだけで数値解析のアカデミックな論文になり、長く歴史に名を残すことになる。１９６０年代から１９７０年代に開発された化石級の計算スキームが現在でも使われているくらいだから（笑）。

 

 

２階線形微分方程式の境界値問題

　　

を解くプログラムのソースコードは次のとおり。ただし、このプログラムでは計算に使用する点の間隔をdxとしてある。

 

 

! 差分法を用いて y''+f(x)y'+g(x)y=h(x) の境界値問題（ディリクレ条件）を解くプログラム
! ただし、固有値問題は解けない!!
parameter (ns=100)
real x(0:ns), y(0:ns)

x = 0.; y =0;

n=10

a=0.; b=1. ! 境界のx座標
x(0)=a; x(n)=b
y(0)=1.; y(n)=exact(1.0) ! 境界条件

call Solver(x,y,n)

write(*,*) ' i      x         y       exact'
do i=0, n
    write(*,100) i, x(i), y(i), exact(x(i))
end do

100 format(i3,1x,f9.5,1x,f9.5,1x,f9.5)
end

function exact(x)
e=exp(1.)
exact=1/(1-e)*exp(-x)- e/(1-e)*exp(-2.*x)
end

function f(x)
f=3.
end

function g(x)
g=2.
end

function h(x)
h=0.
end

! ここより下はいじると危険
! ブラックボックスとして使うべし


subroutine Solver(x,y,n)
real a(n),b(n),c(n),d(n)
real x(0:n), y(0:n)

n1=n-1
dx=(x(n)-x(0))/n

do i=1,n1
    x(i)=x(0)+i*dx
end do

! 差分法によって得られる連立方程式の係数の計算
do i=1,n1
    a(i)=1-dx/2.*f(x(i))
    b(i)=-(2-dx*dx*g(x(i)))
    c(i)=1+dx/2.*f(x(i))
    d(i)=dx*dx*h(x(i))
end do
    d(1)=d(1)-a(1)*y(0) ! x=aにおける境界条件の処理
    d(n1)=d(n1)-c(n1)*y(n) ! x=bにおける境界条件の処理

call tdma(a,b,c,d,n1) ! TDMAで連立方程式を解く

do i=1,n1
    y(i)=d(i) ! 計算結果をセット
end do

end

!  TDMA 連立方程式を解くサブルーティン
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

 

 

　　

について解いた計算結果。