So-net無料ブログ作成
  • ブログをはじめる
  • ログイン
ネコ騙し数学 ブログトップ
前の10件 | -

ニュートン法による非線形連立方程式の近似解法2 [ネコ騙し数学]

ニュートン法による非線形連立方程式の近似解法2

 

ニュートン法を用いると、非線形方程式

の近似解(の一つ)を

と求めることができる。

 

前回、紹介したプログラムでは、偏導関数を与えて、における偏微分係数の計算をした。

この方法では、非線形方程式(1)の関数が変わるたびに、プログラムの利用者は、偏導関数を設定しなおさなければならず、何かと面倒。

 ――偏導関数の計算を間違う場合だってある(^^ゞ――

 

ニュートン法は、近似解法であること、そして、我々が欲しいのは、偏微分係数の正確な値ではなく、を修正するために必要な

であるのだから、偏微分係数

の値は近似値で十分なはずである。

したがって、十分小さなhを選び、中心差分をもちいて

と計算することにする。

h=1/10⁶程度にとれば、(6)式で計算した打切誤差は、せいぜい、1/10¹⁰程度で、その差はほとんど無視できる。

丸め誤差の範囲の程度だにゃ。

 

この方針に従ってあらたにプログラムを作ることにする。

 

 

! ニュートン・ラプソン法で非線形連立方程式f(x,y,z)=0,g(x,y,z)=0,h(x,y,z)=0を解く
! 偏微分係数を中心差分で代用
parameter (n=3)
real a(n,n+1)

limit = 20  ! 繰り返し計算の上限
del=1.e-6
eps=1.e-6 ! 収束判定に使用する何か
a=0. ! 行列の初期化

write(*,*) 'input x0,y0,z0'
read(*,*) x,y,z ! 計算開始時の初期値設定
write(*,*)

do iter=1,limit
    ! 偏微分係数の計算
    fx=(f(x+del,y,z)-f(x-del,y,z))/(2.*del)
    fy=(f(x,y+del,z)-f(x,y-del,z))/(2.*del)
    fz=(f(x,y,z+del)-f(x,y,z-del))/(2.*del)
    gx=(g(x+del,y,z)-g(x-del,y,z))/(2.*del)
    gy=(g(x,y+del,z)-g(x,y-del,z))/(2.*del)
    gz=(g(x,y,z+del)-g(x,y,z-del))/(2.*del)
    hx=(h(x+del,y,z)-h(x-del,y,z))/(2.*del)
    hy=(h(x,y+del,z)-h(x,y-del,z))/(2.*del)
    hz=(h(x,y,z+del)-h(x,y,z-del))/(2.*del)
    
    ! 連立方程式の係数の設定
    a(1,1)=fx; a(1,2)=fy; a(1,3)=fz; a(1,4)=-f(x,y,z)
    a(2,1)=gx; a(2,2)=gy; a(2,3)=gz; a(2,4)=-g(x,y,z)
    a(3,1)=hx; a(3,2)=hy; a(3,3)=hz; a(3,4)=-h(x,y,z)
    
    call gauss(a,n,iflag) ! ガウス消去法で連立方程式を解く
    
    if (iflag.eq.0) then ! エラー判定
        write(*,*) 'エラー!! 連立方程式を解けない'
        stop
    end if
    
    dx=a(1,4); dy=a(2,4) ; dz=a(3,4) ! 修正量 Δx,Δyをセット
    x=x+dx; y=y+dy; z=z+dz ! 修正
    err=amax1(abs(dx),abs(dy),abs(dz))
    
    write(*,100) 'iter=',iter, 'x=',x,'y=',y,'z=',z
    
    if (err.lt.eps) exit ! 収束したかの判定
end do

write(*,*)

if (iter.le.limit) then
    write(*,*) '***solution***'
    write(*,*) 'x=',x,'y=',y, 'z=',z
else
    write(*,*) '収束しない'
end if

write(*,*)
write(*,*) 'check'
write(*,*) 'f(x,y,z)=',f(x,y,z),'g(x,y,z)=',g(x,y,z),'h(x,y,z)=',h(x,y,z)

100 format(a,i3,2x,2(a,f10.6,2x),a,f10.6)
end

! 利用者は連立方程式に合せてここの部分だけを変更
function f(x,y,z)  ! f(x,y,z)=0の関数定義
f=x**3-2*y-2  ! f(x,y,z)=x³-2y-2=0
end

function g(x,y,z)  ! g(x,y,z)=0の関数定義
g=x**3-5*z**2-7  ! g(x,y,z)=x³-5z²-7=0
end 

function h(x,y,z)  ! h(x,y,z)=0 の関数定義
h=y*z**2-1   ! h(x,y,z)=yx²-1=0
end

! 変更してよいのはここまで


subroutine Gauss(a,n,iflag) ! ガウス消去法
real a(n,n+1)

eps = 1.e-10
iflag = 1

! 前進消去
do k=1, n-1
    ! ピボット選択
    i_max = k
    do i=k+1, n
        if (abs(a(i,k)).gt.abs(a(i_max,k))) i_max=i
    end do
    if (i_max.ne.k) then
        do j=k,n
            w=a(k, j)
            a(k, j) = a(i_max,j)
            a(i_max,j) = w
        end do
    end if
    if (abs(a(k,k)).lt.eps) then
        iflag = 0
        return
    end if
    ! ピボット選択終わり

    do i=k+1, n
        t=a(i,k)/a(k,k)
        do j=k+1, n+1
            a(i,j)=a(i,j)-t*a(k,j)
        end do
    end do
end do
if (abs(a(n,n)).lt.eps) iflag = 0

! 後退代入
do i=n,1, -1
    d = a(i,n+1)
    do j= n, i+1, -1
        d=d- a(i,j)*a(j,n+1)
    end do
    a(i,n+1)=d/a(i,i)
end do
end

 

上のプログラムでは

とし、この連立方程式の解(の1つ)を求めている。

 

計算開始の初期値(x⁰,y⁰,z⁰)=(1,1,1)として計算した結果は次の通り。

 

 


nice!(0)  コメント(0) 

ニュートン法による非線形連立方程式の解法 [ネコ騙し数学]

ニュートン法による非線形連立方程式の解法

 

f(x,y)g(x,y)を2回微分可能な関数とする。

  

の近似解法について考えることにする。

 

(x,y)を連立方程式の解、(x⁰,y⁰)を(1)の解の推測値とする。

f(x,y)g(x,y)は仮定より微分可能なので、次のようにテーラー展開することが可能である。

  

したがって、(x⁰,y⁰)(x,y)の近傍の点とすると、

  

と近似することができる。

(x,y)は連立方程式(1)の解なので、f(x,y)=0g(x,y)=0であるから、(4)、(5)式より、(1)式は

  

と書き換えることが可能で、この連立方程式を解くことによって、(1)の解の推測値(x⁰,y⁰)をもとに解(x,y)の推測値を求めることができる。

 

なおここで、

  

である。

 

(6)式を行列を用いて書き換えると、

  

ここで、

  

とおくと、(8)の行列式|J|≠0のとき、(7)式より

  

と解くことができる。

 

(9)式で得られた連立方程式の解(x,y)をあらたに(x⁰,y⁰)にし、(8)、(9)を再計算し、(x,y)を求める。

この操作を繰り返すほど、(1)の近づくことが予想される。

 

この方法は、そのまま、より一般の

  

に拡張することが可能で、この場合は次のようになる。

  Newton2-001.png 

(12)から(非線形)連立方程式の近似解を求めることができるが、(12)を利用する場合、ヤコビ行列Jの逆行列J⁻¹を求める必要がある。しかし、連立方程式(11)を解くよりヤコビ行列Jの逆行列を求めることの方が大変なので――連立方程式を直接解く場合と比較すると、Jの逆行列を求める計算は計算量が約n倍!!――なので、連立方程式

  

をガウス消去法などで解き、

  

と計算するとよい。

 

このように(非線形)連立方程式の近似解を求める方法をニュートン(・ラプソン)法という。

 

参考までに、非線形連立方程式

  

をニュートン法で解くプログラムを示す。

 

 

! ニュートン・ラプソン法で非線形連立方程式f(x,y)=0,g(x,y)=0を解く
parameter (n=2)
real a(n,n+1)

limit = 10  ! 繰り返し計算の上限
eps=1.e-6 ! 収束判定に使用する何か
a=0. ! 行列の初期化

write(*,*) 'input x0,y0'
read(*,*) x,y ! 計算開始時の初期値設定
write(*,*)

do iter=1,limit
    call partial(f,fx,fy,g,gx,gy,x,y)
    ! 連立方程式の係数の設定
    a(1,1)=fx; a(1,2)=fy; a(1,3)=-f
    a(2,1)=gx; a(2,2)=gy; a(2,3)=-g
   
    call gauss(a,n,iflag) ! ガウス消去法で連立方程式を解く
   
    if (iflag.eq.0) then ! エラー判定
        write(*,*) 'エラー!! 連立方程式を解けない'
        stop
    end if
   
    dx=a(1,3); dy=a(2,3) ! 修正量 Δx,Δyをセット
    x=x+dx; y=y+dy ! 修正
    err=amax1(abs(dx),abs(dy))
   
    write(*,100) 'iter=',iter, 'x=',x,'y=',y
   
    if (err.lt.eps) exit ! 収束したかの判定
end do

write(*,*)

if (iter.le.limit) then
    write(*,*) '***solution***'
    write(*,*) 'x=',x,'y=',y
else
    write(*,*) '収束しない'
end if

100 format(a,i3,2x,a,f10.6,2x,a,f10.6)
end

subroutine partial(f,fx,fy,g,gx,gy,x,y) ! 偏微分などの計算
f=x*x+y*y-4 ! f(x,y)=x²+y²−4=0
g=x+2*y-3   ! g(x,y)=x+2y-3=0
fx=2.*x     !∂f/∂x
fy=2.*y     !∂f/∂y
gx=1.       !∂g/∂x
gy=2.       !∂g/∂y
end

subroutine Gauss(a,n,iflag) ! ガウス消去法
real a(n,n+1)

eps = 1.e-10
iflag = 1

! 前進消去
do k=1, n-1
    ! ピボット選択
    i_max = k
    do i=k+1, n
        if (abs(a(i,k)).gt.abs(a(i_max,k))) i_max=i
    end do
    if (i_max.ne.k) then
        do j=k,n
            w=a(k, j)
            a(k, j) = a(i_max,j)
            a(i_max,j) = w
        end do
    end if
    if (abs(a(k,k)).lt.eps) then
        iflag = 0
        return
    end if
    ! ピボット選択終わり

    do i=k+1, n
        t=a(i,k)/a(k,k)
        do j=k+1, n+1
            a(i,j)=a(i,j)-t*a(k,j)
        end do
    end do
end do
if (abs(a(n,n)).lt.eps) iflag = 0

! 後退代入
do i=n,1, -1
    d = a(i,n+1)
    do j= n, i+1, -1
        d=d- a(i,j)*a(j,n+1)
    end do
    a(i,n+1)=d/a(i,i)
end do
end


nice!(0)  コメント(0) 

対流項(移流項)を含む偏微分方程式の数値解法 [ネコ騙し数学]

対流項(移流項)を含む偏微分方程式の数値解法

 

問題 u(x,t)=f(x−ct)であるとき、

  tai-000.png

であることを示せ。ただし、cは定数とする。

【解】

ξ=x−ctとおき、合成関数の微分法を用いると、

  tai-001.png

したがって、

  tai-002.png

(解答終)

 

したがって、fを微分可能な任意の関数とすると、u(x,t)=f(x−ct)

  tai-000.png

の(一般)解になる。

f(x)=sinxとおくと、

  

となり、これが(1)の(特殊)解になることから分かるように、(1)は波の方程式の一つである。

 

(1)式の左辺第1項を

  

と、時間微分に関しては前進差分を用いて近似し、(2)式の左辺の第2項を

  

と、空間微分に関しては中心差分を用いて近似すると、(1)式は

  

と近似することができ、これから

  

という漸化式を得ることができる。

  

とおくと、

  

となる。

(6)式の無次元数(物理の単位をもたない数のこと)をクーラン数という。

 

nami_no_kai.pngしたがって、計算の起点となる、t=t₀におけるuの初期条件x=x₀における初期条件が与えられていれば、(7)式を用い逐次的に(1)式の近似解を求めることができる。

 

ということで、初期条件を

  

さらに、c=0.2Δx=1Δt=1として、(1)を解いた結果は次のようになる。

FTCS.png 

 

赤い曲線はx=1、黄色の曲線はx=2、緑はx=3、茶色はx=4におけるuの数値解の時間tにおける変化を表している。

c=0.2だからx=1に高さ1の波が到達するのはt=5なのにも関わらず、t=Δt=1に既に波の一部(?)が到達しているだけでなく、t=5になっても波の高さは1にならない。しかも波の高さは、約t=20以降、正弦曲線のような曲線を描き、振動する。

x=2x=3x=4と下流にいけばいくほど、波の最大の高さは高くなり、より激しく振動するという、物理的にありえない、数値的に振動する解を(7)式は求めてしまう。

つまり、(1)式の時間微分に前進差分、空間部分に中心差分を用いるFTCSForward in Time Center in Space)は、数値的に不安定な解法で、この解法を(1)式に用いると物理的にありえない振動解を求めてしまうということ。

 

そこで、なぜ、FTCS法は不安定なのか、フォン・ノイマンの安定判定法を用いて調べることにする。

  tai-006.png

とし、(7)式に代入すると、

  

となり、

  

したがって、フォン・ノイマンの安定判定法によると、クーラン数に関わらず、FTCSは無条件不安定!!

 

ということで、FTCSは微分方程式(1)の数値解を求めることに使えないことが明らかになった。

そこで、(1)式の第2項の空間微分を

  tai-007.png

風上化上流化)すると、c>0のとき、(1)は

  

したがって、

  

となる。

(8)式を(13)式に代入すると、

  

よって、

  

したがって、

の時に安定になる。

風上(上流)差分を用いて、c=0.2Δx=1Δt=1、すなわち、クーラン数C=0.2として計算した結果は、次の通り。

 

Up-Wind.png

 

高さ1の波の到達時刻は実際とは異なるけれど、今度は振動することなく計算できていることが分かる。

また、下図を見ると、この波の解は

  tai-005.png

でなければならないのに、風上差分を用いた(13)の数値解は、波の形が崩れていることが分かるだろう。

UpWind-001.png

参考までに、同一の条件でFTCSを用いて計算した結果を下図に示す。

 

FTCS-002.png

 


nice!(0)  コメント(0) 

放物型偏微分方程式の解法 ルンゲ・クッタ法 [ネコ騙し数学]

放物型偏微分方程式の解法 ルンゲ・クッタ法

 

次の熱伝導方程式

  

の右辺を差分法を用いて次のよう

  

と近似し、

とおくと、ルンゲ・クッタ法を用いて、(1)式の近似解を求めることができる。

すなわち、

  

 

例によって、(2)を用いて、次の問題を解くと、右のグラフのようになる。

 

問題

  

を、初期条件

  

境界条件

  

のもとで、Δt=1Δx=1として解け。

 

 

Runge-Kutta-graph-001.png

 

Δt=1Δx=1という粗い計算でも、良好な計算結果が得られていることが分かる。

 

Δt=0.1Δx=0.4とすると、さらに良好な結果が得られ、厳密解との差はほとんど認められなくなる。

 

Runge-Kutt-Graph-003.png

 

なのですが、ルンゲ・クッタ法も陽解法の一種なので、

  

のとき数値解が振動することがあるので、ルンゲ・クッタ法を用いて熱伝導方程式を解くとき、

  

になるようにΔtΔxを取ったほうがよい。

このことを示すために、Δt=0.36Δx=0.5のとき、すなわち、

として計算すると、40ステップのt=14.4の前後で振動を始め、t=18ではさらに激しく振動する。

Runge-Kutta-graph-004.png

Runge-Kutta-graph-005.png

計算に使用したプログラムは次の通り。

 

 

! ルンゲ・クッタ法による解法
parameter (nt=20, nx=10)
real t(0:nx, 0:nt)
real dk(4,0:nx)
real l,kappa

l=4.; kappa=0.5
dt=0.2
dx = l/nx
dx2=dx*dx
c=kappa/dx2

dk= 0
   
! 初期化
do i=0, nx
    do j=0,nt
        t(i,j)=0.
    end do
end do

! 初期条件
do i=0,nx
    x=i*dx
    t(i,0)=x*(l-x)
end do

! 境界条件
do j=0,nt
    t(0,j)=0.
    t(nx,j)=0.
end do

! Runge-Kutta Method
do j=1,nt
    do i=1,nx-1
        dk(1,i) = dt*f(t(i-1,j-1),t(i,j-1),t(i+1,j-1),c)
    end do
   
    do k=2,3
        do i=1,nx-1
            dk(k,i)=dt*f(t(i-1,j-1)+dk(k-1,i-1)/2,t(i,j-1)+dk(k-1,i)/2,t(i+1,j-1)+dk(k-1,i+1)/2,c)
        end do
    end do

    do i=1,nx-1
        dk(4,i)=dt*f(t(i-1,j-1)+dk(3,i-1),t(i,j-1)+dk(3,i),t(i+1,j-1)+dk(3,i+1),c)
    end do
   
    do i=1, nx-1
        t(i,j) = t(i,j-1)+(dk(1,i)+2.*dk(2,i)+2.*dk(3,i)+dk(4,i))/6
    end do
end do

write(6,*) '*** Resutl ***'
do j=0,nt
    write(6,100) j*dt,(t(i,j),i=0,nx)
end do

write(*,*)
write(6,*) '*** Exact ***'
do j=0,nt
    write(6,100) j*dt,(exact(i*dx,j*dt),i=0,nx)
end do

100 format(12(f8.5,1x))

end

function f(x,y,z,c)
f=c*(x-2.*y+z)
end

! 厳密解
function exact(x,t)
pi = acos(-1.0)
s=0.
do i=1,30
    s=s+128/((2*i-1)**3*pi**3)*exp(-(2*i-1)**2*pi**2/32*t)*sin((2*i-1)*pi/4*x)
end do
exact = s
end


nice!(0)  コメント(0) 

[陽解法でも収束するかも] [ネコ騙し数学]

[陽解法でも収束するかも]

 

 ddt^3です。ちょっと面白いなぁ~と思って、差分法に関する安定条件をググってみると、フォン・ノイマンの安定条件ってのが出てきます(註1。

  http://dr-asa.hatenablog.com/entry/2017/08/31/154806

 

 上記HPによると、放物型偏微分方程式の陽解法が安定であるためには、任意の整数0≦mに対して、

 

 

 

 

だそうです。ここでGは時間方向への誤差の拡大率,βmは、

 

 

 

であり、Lは解析領域の長さです。この条件は解の大きさに誤差も比例するという前提なので、mは解を空間方向にフーリエ展開した場合の波数(モード次数)と考えられます。

 к1/2Δt0.1Δx0.4L4として、m012,・・・に対してG(m)をプロットしてみると、図-1になります。

 

kaiketsu-graph-001.png

 

 波数20以上は周期的な繰り返しです。解析結果のグラフをみてみると、どう考えても1次モードのみ励起されるはずです。

 

 

です。

 

 Gは誤差の拡大率です。初期誤差をεとし、誤差は解の大きさに比例するという前提を信じれば、

 

 

 

 

で、時間ステップnにおける誤差を評価できる事になります(たぶん最悪の予想)。T(xt)は解,tnnステップ目の時刻です。

 

kaiketsu-tab-001.png 解析結果は放物線が一様に圧縮されてくような感じなので、T(xtn)T(x0)x2の値の変化で代表させます。紙面から読み取った値は、

 ・・・表-1くらいかな?。

これをグラフに起こして強引に指数関数近似したのが図-2です。相関係数がR21とウルトラ良いので(^^)、図-2の補間結果から、

 

 

で近似する事にします(註2)。

 

kaiketsu-graph-002.png

 

 次に初期誤差εですが、T(x0)x(4x)=-x24xなので、t0でのx方向への2階差分の絶対値はa(0)1のはず。

 時間差分は1階で、空間差分は2階だから、

 

 

くらいかな?(^^;)

 

 以上を使って、誤差の積算Ε(n)を計算してやると、図-3になります。

 

kaiketsu-graph-003.png

 

「ネコの旦那!。陽解法で永遠に数値積分しても、最大誤差0.0102くらいで収束しそうでっせ!」

 

 陰解法は、長時間積分に対して陽解法より平均的に精度が良いだけの話で、局所的な精度は陽解法が上回る場合もあり得ると、自分は思っています。例えばあまり長くない時間積分であれば、無条件安定なシンプレクティック法よりルンゲクッタ法の方が、ずう~っとずう~っと精度が良い。

 シンプレクティック法は陽解法のくせに(まさに陰に)陰解法の特性を持つのに対して、ルンゲクッタ法による超長時間積分では、解は0に減衰します。という訳で・・・。

 

予想1)

 陽解法でも発散しないケースだったから、陽解法の精度がバカみたいに良かった?。

 

 でもですね、n100で既に0.001~最大0.01程度の誤差は見込まれます。単精度ですよね?。問題のスケール規模は1001のオーダーで単精度の有効数字は6桁。倍精度でやりたいなぁ~(^^;)

 

予想2)

 倍精度計算したら、劇的に状況が変化するかも。あくまで「するかも」。

 

 以上、他人のプログラムなんて頼まれても読みたくないddt^3の、他人事のような発言でした(^^;)

 

(執筆:ddt³さん)

 

 

(註1)フォン・ノイマン法による陽解法の安定性判定については、ねこ騙し数学でも過去に「陽解法による拡散方程式の安定性」という記事で取り上げている。

http://nekodamashi-math.blog.so-net.ne.jp/2016-11-23

 

(註2)微分方程式の厳密解は、

そして、n=1のとき

となるので、約5%くらいの誤差で

と近似できる。

 


nice!(0)  コメント(0) 

放物型偏微分方程式の解法4 まとめ [ネコ騙し数学]

放物型偏微分方程式の解法4 まとめ

 

熱伝導方程式

  

を、陽解法を用いて離散化すると、

  

(純)陰解法を用いると、

  

となる。

そして、クランク・ニコルソン法の場合は、(2)、(3)式の辺々を足しあわせ2で割った

  

 

である。

ここで、

  

 

また、(2)、(3)、(4)式は

と一つの式表すことができ、

  

になる。

(6)式を

  

と変形し、

  

とおくと、(7)式は

  

となり、陽解法、陰解法、クランク・ニコルソン法と解法の区別をすることなく、三重対角行列を係数に持つ連立方程式として(1)の近似解を求めることができる。

 

この方針にしたがって作ったプログラムは以下の通り。

 

parameter (nx_max=20, nt_max=100)
real t(0:nx_max,0:nt_max)
real kappa,l

nt = 10; nx = 8
l=4.0
dt = 0.25
dx = l/nx
dx2 = dx*dx
kappa=0.5


write(6,*) '解法を入力:陽解法 1, 陰解法 2, クランク・ニコルソン 3'
read(5,*) iflag

if(iflag.eq.1) then
    alpha=0.
else if (iflag.eq.2) then
    alpha = 1.
else
    alpha=0.5
end if


do i=0, nx
    do j=0, nt
        t(i,j)=0.
    end do
end do

! initial conditon
do i=0, nx
    x=i*dx
    t(i,0)=x*(l-x)
end do

! boundary condition
do j=1, nt
    t(0,j)=0.
    t(nx,j)=0.
end do

r= kappa*dt/dx2

call solver(t,nx,nt,r,alpha)

write(*,*) '*** Result ***'
do j=0,nt
    write(*,100) j*dt, (t(i,j), i=0,nx)
end do

write(*,*)
write(*,*) '*** exact solution ***'
do j=0,nt
    write(*,100) j*dt, (exact(i*dx,j*dt), i=0,nx)
end do

!close(1)
100 format(12(f8.5, 1x))
end

! 厳密解
function exact(x,t)
pi = acos(-1.0)
s=0.
do i=1,30
    s=s+128/((2*i-1)**3*pi**3)*exp(-(2*i-1)**2*pi**2/32*t)*sin((2*i-1)*pi/4*x)
end do
exact = s
end


! 【注意】 ここから下はいじってはいけない
! 差分法を用いて解く本体
subroutine solver(t,nx,nt,r,alpha)
real t(0:20,0:100)
real a(nx), b(nx), c(nx), d(nx)

do j=1, nt
    do i=1,nx-1
        a(i)=-alpha*r
        b(i)=1+2.*alpha*r
        c(i)=-alpha*r
        d(i)=t(i,j-1)+(1-alpha)*r*(t(i-1,j-1)-2*t(i,j-1)+t(i+1,j-1))
    end do
    d(1)=d(1)-a(1)*t(0,j)
    d(nx-1)=d(nx-1)-c(nx-1)*t(nx,j)
    call tdma(a,b,c,d,nx-1)

    do i=1,nx-1
        t(i,j)=d(i)
    end do
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

 

陽解法の時にも連立方程式を解いているので、陽解法のとき、少し無駄な計算をすることになるけれど、陽解法、陰解法、クランク・ニコルソン法といった解法の区別をすることなく計算できる。

また、αを一種の重みと考え、α01/21だけではなく、0≦α≦1の任意のαを選び計算することもできる。

ただ、

にとると、

  

の値によって、伝播誤差のために、振動解が得られる。

 

例によって、(1)のκ=1/2とし、条件

  

とし、Δt=0.1Δx=0.4に選んで計算した結果を以下に示す。

陽解法
Explicit.png

陰解法

Implicit.png

クランク・ニコルソン法

Crank-Nicolson method.png

 

ΔtΔxをこの程度にとれば、陽解法、陰解法、クランク・ニコルソン法ともに、(1)の厳密解とよく一致することが分かる。

 

また、Δt=0.2Δx=0.4とし、陽解法で解いた結果は次の通り。

 

Explicit2.png

 

このとき、

  

となるので、この条件で陽解法を用いて解こうとすると、伝播誤差のために、計算を進めるたびに、数値解の振動が激しくなり、厳密解からの逸脱が大きくなることが分かる。

 


nice!(0)  コメント(0) 

放物型偏微分方程式の解法3 クランク・ニコルソン法 [ネコ騙し数学]

放物型偏微分方程式の解法3 クランク・ニコルソン法

 

微分方程式

  

を、陽解法を用いて離散化すると、

  

となる。

陰解法を用いて離散化すると、

  

(2)と(3)の平均を取ると、

  

になる。

  

とおくと、

  

になる。

ここで、

とおくと、

   

と3重対角行列を係数に持つ連立方程式が得られ、三重対角行列アルゴリズムTDMAを用いてこの連立方程式を解くことが可能になる。

 

大胆に、(6)式の以外の項を落とすと、

  

となり、r>0のとき、

  

となり、クランク・ニコルソン法は無条件安定であることが予想できる。

 

例によって、次の問題をクランク・ニコルソン法で解いてみる。

 

問題

  

を、初期条件

  

境界条件

  

のもとで、Δt=1Δx=1として解け。

 

 

ht-graph-004.png

 

Δt=1Δx=1という粗い計算でも、良好な計算結果が得られていることが分かる。

 

Δt=0.5Δx=0.4にすると、厳密解とクランク・ニコルソン法による近似解はほとんど等しくなる。

 

ht-graph-010.png

ht-graph-011.png

ht-graph-012.png

 

この計算に使用したプログラムは次の通り。

 

! crank-nickolson method
parameter (nt = 20, nx = 10)
real t(0:nx,0:nt)
real a(nx),b(nx),c(nx),d(nx)
real kappa,l

l=4.0
dt = 0.5
dx = l/nx
dx2 = dx*dx
kappa=0.5

do i=0, nx
    do j=0, nt
        t(i,j)=0.
    end do
end do

! initial conditon
do i=0, nx
    x=i*dx
    t(i,0)=x*(l-x)
end do

! boundary condition
do j=1, nt
    t(0,j)=0.
    t(nx,j)=0.
end do

r= kappa*dt/dx2
do j=1, nt
    do i=1,nx-1
        a(i)=-r/2
        b(i)=1+r
        c(i)=-r/2
        d(i)=t(i,j-1)+r/2*(t(i-1,j-1)-2*t(i,j-1)+t(i+1,j-1))
    end do
    d(1)=d(1)-a(1)*t(0,j)
    d(nx-1)=d(nx-1)-c(nx-1)*t(nx,j)
    call tdma(a,b,c,d,nx-1)

    do i=1,nx-1
        t(i,j)=d(i)
    end do
end do

write(*,*) '*** Crank-Nicolson Method ***'
do j=0,nt
    write(*,100) j*dt, (t(i,j), i=0,nx)
end do

write(*,*)
write(*,*) '*** exact solution ***' 
do j=0,nt
    write(*,100) j*dt, (exact(i*dx,j*dt), i=0,nx)
end do

100 format(12(f8.5, 1x))

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

! 厳密解
function exact(x,t)
pi = acos(-1.0)
s=0.
do i=1,10
    s=s+128/((2*i-1)**3*pi**3)*exp(-(2*i-1)**2*pi**2/32*t)*sin((2*i-1)*pi/4*x)
end do
exact = s
end

 


nice!(0)  コメント(0) 

放物型偏微分方程式の数値解法2 [ネコ騙し数学]

放物型偏微分方程式の数値解法2

 

前回に引き続き、差分法を用いた熱伝導方程式の数値解法について述べることにする。

  

この方程式の時間微分を

  

右辺の空間微分を

  

と近似すると、(1)式は

  

となる。

  

とおくと、

  

という3重対角行列を係数を持つ連立方程式が得られる。

 

(1)の解をT=T(x,t)とし、初期条件を

  

境界条件を

  

とする。

そして、

  

とすると、初期条件より

  

境界条件より

  

となる。

したがって、

  

という連立方程式が得られ、逐次的に、3重対角行列アルゴリズム(TDMA)を用いて解くことによって、微分方程式(1)の近似解を求めることができる。

 

このような解法を陰解法という。

陰解法は、前回紹介した陽解法とは異なり、無条件安定である。

 

陽解法との比較参照のために、前回と同じ問題

 

問題

   

を、初期条件

  

境界条件

  

のもとで、Δx=1Δt=1として解け。

 

を解かせたものは次の通り。

ht-graph-002.png 

この問題の場合、厳密解と比較すると、陰解法による近似解は厳密解よりも大きいことがわかる。

さらに、前回の陽解法による結果を加え比較したものが次のグラフである。

 

ht-graph-003.png

込み入っていて、少しわかりづらいと思うのだけれど、厳密解は陰解法と陽解法の間に位置している。

ということで、陰解法と陽解法によって得られた数値解を足して2で割ると厳密解に近い結果が得られることが予想できる。

この予想は正しく、陽解法と陰解法の中間的な解法、クランク・ニコルソン法という解法が存在する。

 

ht-graph-004.png

 

Δt=1Δx=1という非常に粗い計算でありながら、クランク・ニコルソン法は、厳密解とよく一致していることが分かる。

 

ただし、この事実をもって、陽解法や陰解法は精度が悪いと思わないで欲しい。

Δt=0.1程度にとれば、陽解法、陰解法ともに良好な結果を得ることができる。

陽解法はともかく、陰解法は熱や流れなどの数値計算ではよく使われている、優れた解法である。

 

問題の初期条件を

  

境界条件を

  

と変え、陰解法を用いてΔt=0.5Δx=0.4として計算した結果を以下に示す。

このようにあらい格子でも、厳密解をよく再現していることが分かる。

ht-graph-007.png

ht-graph-008.png

ht-graph-009.png

 

 

そして、t→∞の定常解

  

に収束することが確認できる。

誤差は最大で10%くらい。

Δt=0.1Δx=0.4にすると、誤差は最大で約3%。

 

陰解法よりより高精度な計算ができる、クランク・ニコルソン法は、次回ということで。

 

計算に使用したFortranプログラムは次の通り。

 

parameter (nt = 20, nx = 10)
real t(0:nx,0:nt)
real a(nx),b(nx),c(nx),d(nx)
real kappa,l

l=4.0
dt = 0.5
dx = l/nx
dx2 = dx*dx
kappa=0.5

do i=0, nx
    do j=0, nt
        t(i,j)=0.
    end do
end do

! initial conditon
do i=0, nx
    x=i*dx
    t(i,0)=x*(l-x)+x/l
end do

! boundary condition
do j=1, nt
    t(0,j)=0.
    t(nx,j)=1.
end do

do j=1, nt
    do i=1,nx-1
        a(i)=-kappa/dx2
        b(i)=2.*kappa/dx2 + 1/dt
        c(i)=-kappa/dx2
        d(i)=t(i,j-1)/dt
    end do
    d(1)=d(1)-a(1)*t(0,j)
    d(nx-1)=d(nx-1)-c(nx-1)*t(nx,j)
    call tdma(a,b,c,d,nx-1)

    do i=1,nx-1
        t(i,j)=d(i)
    end do
end do

write(*,*) ' *** result *** '
do j=0,nt
    write(*,100) j*dt, (t(i,j), i=0,nx)
end do

write(*,*)
write(*,*) 'exact solution' 
do j=0,nt
    write(*,100) j*dt, (exact(i*dx,j*dt), i=0,nx)
end do

write(*,*)
write(*,*) 'error(percent)' 
do j=1,nt
    write(*,100) j*dt, ((t(i,j) - exact(i*dx,j*dt))/exact(i*dx,j*dt)*100, i=1,nx-1)
end do
100 format(12(f8.5, 1x))
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

function exact(x,t)

pi = acos(-1.0)
s=0.
do i=1,10
    s=s+128/((2*i-1)**3*pi**3)*exp(-(2*i-1)**2*pi**2/32*t)*sin((2*i-1)*pi/4*x)
end do
exact = s+x/4.
end


nice!(0)  コメント(0) 

放物型偏微分方程式の数値的解法1 [ネコ騙し数学]

放物型偏微分方程式の数値的解法1

 

grid--001.png放物型偏微分方程式の熱伝導方程式

  

を例に取り、差分法を用いた数値的な解法について考えることにする。

 

(1)の左辺の時間(偏)微分

  

右辺の2階の空間(偏)微分を

  

と、差分法を用いて近似すると、

  

ここで、

  

とおくと、

  

となる。

を計算する際、はすべて既知なので、(5)式を用いて逐次的にを計算することができる。

とくに、r=1/2のとき、(5)式は

  

になる。

 

このような解法を陽解法という。

陽解法は代数方程式を解くことなく簡単に計算できるが、

  

のとき、(5)を用いて求めた近似解は安定ではなく、振動解が得られる。

したがって、

  

となるようにΔtΔxを選ばないといけない。

 

(5)式の右辺を

  

と書き換え、右辺第2項を無視すると――なんと大胆な(^^ゞ――

  

誤差がこれにしたがって伝播するとする、安定であるためには、

  

でなければならい。

したがって、安定であるためには

  

・・・。

もう少し正確な議論をすると、誤差

  

に従うとする。

すると、

  

だから、少なくとも、安定であるためには

  

でなければならないに違いない。この条件を満たさないと、近似解は振動したり、発散するだろう。

そして、r>0だから、

  

 

正確な議論をするためには、von Neumannの判定法などを用いる必要があるが、それは厄介なので、簡易的にこの関係を求めてみたにゃ。

 

 

問題

  

を、初期条件

  

境界条件

  

のもとで、Δx=1として解け。

【解】

Δt

  

となるように、Δt=1にとると、

  

となる。

  

以下、同様に計算すると、次の表が得られる。

 

この結果をグラフにすると、次のようになる。

 ht-graph-001.png

定量的にはともかく、このような粗い計算であっても、指数関数的な現象を捉えており、定性的には正しい結果が得られていることがわかる。

 

より精度よく計算するために、Δx=0.1とすると、条件(7)より

  

よって、Δt0.01以下にとる必要があり、計算量が大きく増えてしまう。

したがって、実際は、熱伝導方程式を陽解法を用いて解くことはない。

 

こうした制約のない陰解法を次回紹介することにする。



nice!(0)  コメント(0) 

平行平板ポアズイユ流れの熱伝達のプログラム(最終版) [ネコ騙し数学]

! 2平板間の流れ(平板ポアズイユ流れ、壁温一定)を解くプログラム
! NS eq. U(∂U/∂X)+V(∂U/∂Y)=-dP/dX+1/Re(∂²U/∂²Y)
! 連続の式 ∂U/∂X+∂V/∂Y=0
! エネルギー方程式 U(∂T/∂X)+V(∂T/∂Y)=1/(RePr)(∂²U/∂²Y)
parameter (m=10,n=5,n1=6) ! n1=n+1
real u(0:m,0:n),v(0:m,0:n), dpdx(m)
real a(n),b(n),c(n),d(n+1),e(n),f(n+1)
real t(0:m,0:n)
real tb(m), qnu(m)
real dummy_v(n);
iter_limit = 200 ! 流れ場の反復計算の反復回数の上限
eps = 1.0e-6 ! 流れ場の収束判定

xl=10
re=100.; pr =0.7
dy=0.5/n
dy2=dy*dy
dx=xl/m

! 速度u,vの初期化
do i=0,m
    do j=0,n
        u(i,j) =0.0 ; v(i,j)=0
    end do
end do

! 流入速度(X=0)の速度
do j=1,n
    u(0,j)=1.0
end do

do i=0, m ! 温度の初期化
    do j=1,n
        t(i,j)=0
    end do
end do

do i=1,m ! 壁の温度
    t(i,0)=1
end do


do i=1,m
    x= x+dx
    do iter=1, iter_limit
    ! NS方程式の係数行列の計算と値のセット
        do j=1,n
            a(j)=-v(i,j)/(2.*dy)-1./(re*dy2)
            b(j)=u(i,j)/dx+2./(re*dy2)
            if (j.ne.n) then
                c(j)=v(i,j)/(2*dy) -1./(re*dy2)
            else
                c(j)=0.
            end if
            d(j)=u(i-1,j)/dx*u(i,j)
            e(j)=1.
            f(j)=dy
        end do
        a(n)=a(n)-1./(re*dy2)
        d(n1)=0.5
        f(n)=dy/2
        f(n1)=0.

        ! 速度と圧力を解くプログラムを呼び出す
        call calcup(a,b,c,d,e,f,n,n1)
        do j=1,n
            u(i,j)=d(j)
        end do
        err=abs(1-dpdx(i)/d(n1))
        dpdx(i)=d(n1)
       
        ! 連続の式からvを計算 (ここは改良すべき・・・)
        do j=1,n-1
!            v(i,j) = v(i,j-1)+dy/dx*(u(i-1,j)-u(i,j))
            um=0.5*(u(i-1,j) + u(i-1,j-1))
            up=0.5*(u(i,j) + u(i,j-1))
            v(i,j)=v(i,j-1)+dy/dx*(um-up) ! 少し改良・・・
        end do
        if (err.lt.eps) exit
    end do
end do


! エネルギー方程式の連立方程式の係数の計算
do i=1, m
    do j=1,n
        a(j)=-v(i,j)/(2.*dy)-1./(re*pr*dy2)
        b(j)=u(i,j)/dx+2./(re*pr*dy2)
        if (j.ne.n) then
            c(j)=v(i,j)/(2*dy) -1./(re*pr*dy2)
        else
            c(j)=0.
        end if
        d(j)= u(i,j)/dx*t(i-1,j)
    end do
    d(1)=d(1)+1./(re*pr*dy2)*t(i,0)
    a(1)=0.
    a(n)=2*a(n)
    c(n)=0.

    call tdma(a,b,c,d,n) ! 連立方程式を解く
   
    do j=1,n
        t(i,j)=d(j)
    end do
end do

! 台形公式でバルク温度(混合平均温度)の計算
    do i=1,m
        s=0.5*u(i,n)*t(i,n)
        do j =1,n-1
            s = s + u(i,j)*t(i,j)
        end do
        tb(i)=2*s*dy;
    end do

! 局所の熱伝導量、並びに局所ヌセルト数を計算
    do i=1, m
        dtdy= (4*t(i,1)-t(i,2)-3*t(i,0))/(2*dy)
        qnu(i)=-dtdy/(1-tb(i))
    end do

write(6,*) ' ***** 計算結果 ***** '
write(6,*) 'U'
do i=1, m
    write(6,100) i*dx,(u(i,j),j=1,n), dpdx(i)
end do
write(6,*)
write(6,*) 'V'
do i=1, m
    write(6,100) i*dx,(v(i,j),j=1,n)
end do

write(6,*)

write(6,*)
write(6,*) 'T'
do i=1, m
    write(6,100) i*dx,(t(i,j),j=1,n), qnu(i)
end do

100 format(f8.5,1x, 12(f8.5,1x))
110 format(a, 1x, 5(f8.5,1x),a)

end

! UとdP/dXを解くサブルーティン
subroutine calcup(a,b,c,d,e,f,n,n1)
real a(n), b(n), c(n),d(n1),e(n) ,f(n1)

! 前進消去
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)
    d(n1)=d(n1)-f(i)/b(i)*d(i)
    e(i+1)=e(i+1)-ratio*e(i)
    f(i+1)=f(i+1)-f(i)/b(i)*c(i)
    f(n1)=f(n1)-f(i)/b(i)*e(i)
end do
f(n1)=f(n1)-f(n)/b(n)*e(i)
d(n1)=d(n1)-f(n)/b(n)*d(n)

! 後退代入
d(n1)=d(n1)/f(n1)
d(n)=(d(n)-e(n)*d(n1))/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1)-e(i)*d(n1))/b(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

出力部は自分で作るといいにゃ。そこまでは面倒を見きれないケロよ。

mはX方向の分割数、nはY方向の分割数、n1=n+1だケロ。Xの最大値は10だから、m=100、n=10くらいが丁度いいんじゃないかな。mを1000、nを100にしても、m×nオーダーの計算法だから、エンターキーを押した瞬間に終わると思うけれど、U、V、T合せて30万のデータ数になるから、mとnをあまり大きくすのはやめたほうがいいと思うにゃ。


タグ:数値解析
nice!(0)  コメント(0) 
前の10件 | - ネコ騙し数学 ブログトップ

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。