So-net無料ブログ作成
数値解析 ブログトップ
前の10件 | -

放物型偏微分方程式の解法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:nx,0:nt)
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になるのはΔt=0のときだから、これを除くと(7)が得られる。

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

  

に従うとする。

すると、

  

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

  

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

そして、r>0だから、

  

 

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

 

 

問題

  

を、初期条件

  

境界条件

  

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

【解】

Δt

  

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

  

となる。

  

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

 

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

 ht-graph-001.png

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

 

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

  

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

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

 

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

 

 


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

[ガラーキン法がいつも最良って訳じゃないんですよ] [数値解析]

[ガラーキン法がいつも最良って訳じゃないんですよ]

 

※どうでも良いトリビア:

 ロシア(ソ連)では、カントロヴィッチの方法っていう時もあるみたいです。

※次に愚痴:

 ネコ先生だって語れるはずなのに。でもお声がかかると嬉しくなって書いちゃう(^^)

 

 重み付き残差法は、変分原理のない系にも有限要素法を簡単に定式化するために開発された方法なので、いちおう変分原理とは別の根拠に基づきます。関数空間の理論です。関数空間って何なのさ?って話になる訳ですが、具体的にはフーリエ級数なんかでいつもやってる事です。

 任意の連続関数は、収束半径内で任意の独立な関数系で展開できる事がわかってます。ここで関数系とは、単に関数の集合{f0(x)f1(x),・・・}の事で、独立とは各fk(x)kを除く{fj(x)}で展開できない事を言います。また関数f(x)の展開とは適当な定数a0a1,・・・を用いて、

  

と出来ることです。独立な関数系は必ず可算無限個(有限個または自然数程度の無限個のこと)の関数を持ちます。

 定義として述べるとなんか鬱陶しいですが、良く目にするのがテーラー級数です。

  

 式(2)の関数系は{1xx2,・・・}で可算無限個の関数を持ち、xkkを除く{xj}で展開できないのは明らかです。また、

  

という事になります。式(2)の収束半径はf(x)の性質によって決まります。ここはちょっと不便です。

 

 不便なので収束半径無限大の関数系はないものか?とさがします。そんな便利なものはない訳ですが、実質的に収束半径無限大とみなせるものならあります。三角関数系、

  gala-001.png

です。(4)j0を除き、fj(x)(2π)-1/2 sin(2πjx/L)と表せます。(4)による展開は、

  gala-002.png

になります。式(5)の右辺は0≦x≦Lの範囲で、必ず収束するのがわかってます。ここでLは任意です。

 実用上-∞≦x≦+∞というf(x)の定義域はあり得ません。よって任意の0≦x≦L(5)が保証されるなら、実質的に収束半径無限大と同じです(^^)

 しかし(5)では、2つの事が不明です。三角関数系の独立性と、係数a0a1,・・・の与え方です。これらは三角関数系の正規直交性という性質から解決されます。三角関数系のメンバーを最初と同じように{f0(x)f1(x),・・・}で表すと、

  

なのは、直接計算してすぐわかります(三角関数の積和、倍角公式を使用)。δkjはクロネッカーのデルタ。

  gala-009.png

 そうすると式(5)の左辺がfk(x)である時、両辺にj≠kとなるfj(x)をかけて区間[0L]で積分すれば、

  

なので、(6)の関係から左辺は0,右辺はajの項のみ残ってaj0です。jjkでなければ、任意でした。これは、fk(x)kを除く{fj(x)}(1)のように展開しようとしても、係数ajが全て0になってしまい、fk(x){fj(x)}で展開できなかった事を意味します。従って三角関数系は独立。

 次に(5)に戻って、同様に(5)の両辺にfk(x)をかけて区間[0L]で積分すれば、

  

なので右辺にはakのみ残り、

  

が得られます。kは任意なので、係数a0a1,・・・の与え方もわかりました。

 ちなみに式(6)で∫fk(x)fk(x)dx1になる事を正規性,j≠kで∫fk(x)fj(x)dx0になる事を直交性と呼びます。

 

 ところで式(8)からの(9)の導出を、もうちょっと一般化できるのはすぐわかります。正規性と直交性を外して関数系に独立性だけを要求します。そうすると(8)は、

  

と書けます。ただし、

  

の事だと定義します。関数空間では(11)の形を内積と呼びます(←理由はきかないで(^^;))。(fkfj)(fkf)は、fk(x)fj(x)f(x)も既知関数なので(11)で必ず具体的数値になります。そしてk012,・・・なのでした。kはどこまでも続きますが、実用上そうなる事はあり得ません。必ずあるnで打ち切ります。従って(10)は、

  gala-005.png

という展開係数(a0a1,・・・,an)に関する連立方程式を表している事になります。この連立方程式の左辺の係数マトリックスが特異でないのを決める条件が、じつは関数系の独立性です。

 ところで(12)で決定された係数(a0a1,・・・,an)で展開されたf(x)の近似、

  

は収束するんでしょうか?。しますよ。「有限個」しか足してないんだから(^^)。これの含むところは、展開範囲[0L]を本来の収束半径を越えて設定すると、近似次数nを上げるほど精度が悪くなる事もあり得る、という事です。言いたいのは、近似次数はほどほどに、展開範囲もほどほどに・・・(^^;)

 そういう訳で実用上、収束半径を余り気にする事無く、正規直交でない関数系{1xx2,・・・}なんかにも同様な計算が可能になります。

 

 ・・・これってガラーキン法ですよね?。そして明らかにもっと一般化できますよね。関数系{fj(x)}に対して別の独立な関数系{gj(x)}を持ってきて、

  gala006.png

を作れば、とにかく展開係数(a0a1,・・・,an)は計算できます。・・・これって重み付き残差法ですよね?(^^)。重み付き残差法は関数解析における、最も一般的形と思って間違いではないと思います。

 gをデルタ関数にとれば選点法,{xj}にとればモーメント法,{∂L/∂aj}にとれば最小二乗法です(Lfに対する微分作用素)。実用的には積極的反対理由がない限り、場合に応じて一番便利なものを使います(^^;)

 

 ちなみに関数系の独立性は、一般にはさっきみたいに簡単にはいきません。またここまで可算無限個のメンバーを持つ独立な関数系は、任意の連続関数を展開できるとしてやってきましたが、この性質を完全性と言います。完全性の証明はさらに困難です。なので普通は気にせずに重み付き残差法を使います。というか、昔の人たちが苦労して証明してくれた正規完全直交関数系の数々が、「岩波公式集」なんかに特殊関数として山のように載ってます(^^)

 

 最後に微分方程式に対して変分原理が存在する場合です。変分原理を使うためにはラグラジアンF(汎関数ともいう)を見つける必要があります。

  

ならラグラジアンFは、

 

  

です。※なれれば少なくとも線形系に関しては、すぐFを書けるようになります。変分原理は、

  

と書けます。以下の変形は本質的に、境界要素法入門の最初と同じです。まず変分を取ります(u'du/dx)。

  

 変分法の定石は、変分δuの微分を部分積分で消去する、でした。

  gala-siki-19.png

 ここでu(x)を式(13)の形で近似するとします。たいがい境界条件を満たす{fj(x)}を使うので、普通は(19)の一項目は0になります。実際にやるとわかりますが、Fとして(16)をとった場合、(19)()内は(15)に戻ります。それをL(u)x0と書くと、

  

です。次に変分δuの具体的な姿を考えます。u(x)は式(13)の形で近似すると「決めた」ので、じつは任意の変分と言いながら実際には、任意に動かせるのは展開係数(a0a1,・・・,an)しかありません。

  

 各δajには「微小」という制限は付きますが、その範囲では完全に自由です。自由なので、jkを除いてδaj0とする事も可能です。もしそうすれば、

  

ですが、δakは定数なので両辺をδakで割れます。従って、

  

  ※ k=1,2,・・・,n

となり、ガラーキン法の条件が自然に導かれます。

 

 これが「積極的反対理由」です。変分原理が使える場合ガラーキン法で近似すると、変分原理を概ね満たしながら収束する解になります。そして変分原理とそれから導かれる微分方程式は、うるさい事を言わなければ同値です。一方ガラーキン法以外では、必ずしも変分原理を満たさない解の収束になります。実用上完全に解を収束させる事はできません(どこかで打ち切るから)。従って変分原理が使える場合はガラーキン法の精度が一番良いだろうと、予想できます。もちろん厳密な関数解析による証明はありますよ。でも余りに面倒で自分は読んでません(^^;)

 

 という訳で、変分原理を使えないケースではガラーキン法がいつも最良って訳じゃありません。理想的には問題に適した近似法を選択するべきですが、どれを使っても経験上そんなにかけ離れた結果にはならず、面倒なので定式化と計算が一番簡単なガラーキン法を使っちゃいますね(^^)

 


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) 

陰的オイラー法 [数値解析]

陰的オイラー法

 

微分方程式

  

の解は

  

で与えられる。

この式の右辺の定積分を

  

と近似し、

  

と逐次的に微分方程式の近似解を求める方法を陰的オイラー法という。

 

等間隔の場合、すなわち、

  

のとき、

  

である。

 

  

このとき、陰的オイラー法は

  

となり、この漸化式を解くと

  

になる。

また、誤差も漸化式に従うとすれば、

  

のとき、すなわち、

  

のとき、陰的オイラー法は安定した解法になる。

 

λ=3h=1のとき、λh=3>2となり陰的オイラー法は安定であり、

  

に収束する。

しかし、この条件のときの微分方程式の解は

  

であり、陰的オイラー法による近似解は厳密解の挙動を表していない。陰的オイラー法による近似解が安定でかつ意味を持つのは、λh<0のときである。

数値計算でいう安定とは収束解が得られるという意味であり、その収束解が微分方程式の解であることを意味しないことに注意。

陰的オイラー法は、陽的オイラー法と同様に打切誤差がO(h²)程度で、しかも、方程式を解かなければならないので、実際に使われることはない。

たとえば、

  

の場合、

  

を解かなければならない。

この方程式は二分法やニュートン法を用いて数値的に解くことができるが、各段階でこの計算を行わなければならないので、計算量が多くなってしまう。理論的は話ではともかく、陰的オイラー法の打切誤差はO(h²)程度なので、とてもじゃないけれど、こんな計算は馬鹿らしくてやってられない。

 

参考までに、λ=1h=0.1のときに、陽的オイラー法と陰的オイラー法を用いて(6)を数値的に解いた計算結果を示す。

 

 

 

この場合、陽的オイラー法の方が陰的オイラー法よりも精度よく計算できていることが分かる。

 

λ=−1h=0.1のときは、次のようになる。この場合は、陰的オイラー法の方が精度よく計算できていることが分かるだろう。

 

 

 

また、一般的な傾向として、厳密解と比較した時、陰的オイラー法は大きめの値を出し、陽的オイラー法は小さめの値を出し、厳密解はこの2つの間に位置し、したがって、陽的、陰的オイラ法の近似解の平均をとれば厳密解に近い値を求められそうである。

そこで、λ=1h=0.1の場合について計算してみると、次のようになる。

 

 

 

予想通り、劇的に精度が向上したにゃ(^^)

 

 


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

オイラー法の誤差の内訳 [数値解析]

オイラー法の誤差の内訳

 

この記事を書いた私自身、このあたりの話はよく分かっていないのだけれど・・・。

 

微分方程式を

  

とする。

(1)式の点における近似解、厳密解をとするとき、で発生した誤差

  

は、オイラー法の場合、

  

次のステップのの近似解に伝播していく。

 

の場合、

  

だから、

  

である。

そして、この伝播誤差の他に、k+1段の積分計算にともなう打切誤差が加わり、

  

となり、k+1段に

  

として伝わってゆくのでしょう。

 

この考えに基づき、誤差の蓄積の内訳を求めてみた。

λ=1h=0.1の場合は次の通り。

 

 

Euler-gosa-005.png 

Euler-gosa-006.png

 

オイラー法の場合、伝播誤差をまったく含まないx=0.1を除き、誤差の大部分を伝播誤差が占めており、この伝播誤差のために、計算を進めるたびに厳密解との乖離が大きくなってゆくことが分かる。

 

 

安定であるλ=−1h=0.1の場合は、次のようになる。

 

 

Euler-gosa-007.png

 

Euler-gosa-008.png

 

この場合、伝播誤差(の絶対値)はx=1.2の付近まで増加するのは、伝播誤差の減衰よりも計算のたびに加わる打切誤差の(絶対値)の増加が大きいため。

この場合も、誤差の多くを占めていることが分かる。

 

λ=1h=0.1のとき、こうして見積もった打切誤差の積算値を求めると、x=2のときに、およそ、−0.314。対して、厳密解とオイラー法を用いて求めた近似解との相対誤差は−0.662で2倍以上になっている。このことからも、積分計算をするたびに発生する打切誤差が増幅されて伝播していることが分かる。

 

 

 

なお、この計算においては、打切誤差は、相対誤差と伝播誤差の差から求めてある。(計算法の詳細は補足を参照!!)

ネムネコは数値計算の素人だから、こうした計算法の是非についてはわからない。そして、数値計算屋さんと伝播誤差を違う意味で捉えているのかもしれない。

だから、この方面に詳しいヒトや数値計算屋さんに、こうした計算の是非などについて教えて欲しいにゃ。

 

 

補足


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

誤差の蓄積 [数値解析]

誤差の蓄積

 

  

の解は

  gosa-001.png

であり、Euler法は(2)式の積分を

  gosa-002.png

と近似し、

  

を計算し、以下、逐次的に

  

と計算し、微分方程式(1)の近似解を求める方法である。

特に、

  

と、等間隔の時、オイラー法は

  

となる。

そして、(5)の(局所的な)打ち切り誤差程度である。コンピュータを用いる解法では、この打ち切り誤差の他に計算機特有の丸め誤差などが発生するが、丸め誤差は打ち切り誤差と似た性質を持っているので、打ち切り誤差に含めることが可能であろう。

 

が正確にであるとき、積分によって得られるの近似値の間には

  

という関係が成立する。ここで、は局所的な打ち切り誤差を表す。

さらに、が誤差を含めば、にはこの誤差伝播誤差と局所的な打ち切り誤差が入る。この2つの誤差によって生じる誤差は、計算をするたびに蓄積し、時に近似計算の結果を無意味なものにしてしまう。

 

微分方程式(1)、すなわち、y’=f(x,y)のオイラー法による積分(5)において、が正確であれば、

  

である。

ここで、

  

の誤差を含めば、

  

右辺第2項をテーラー展開し、は十分小さく高次の項を無視できるとすると、

  

(10)式を(9)式に代入すると、

  

したがって、伝播誤差は

  

である。

よって、では、

  gosa-004.png

となり、では、

  

となる。

のとき、計算を進めるたびに誤差は増大してゆく。そして、hが小さければ、

  

であれば、誤差の影響は積分を繰り返すたびに小さくなり、やがて消失してしまう。

 

次の微分方程式があるとする。

  

この微分方程式の解は

  

である。

(13)式にEuler法を用いると、

   

となる。

誤差も(14)式に従うとすれば、

  

したがって、誤差因子1+λh

  

すなわち、

  

であれば、(伝播)誤差は常に一定の範囲に収まる、つまり、有界で、安定である。

λ>0のとき、(16)を満たすことはできないので、(積分の)計算を進めるたびに、伝播誤差が指数関数的に増大する。

 

参考までに、λ=±1h=0.1としたときの結果をグラフに示す。

λ=1のとき、誤差が指数関数的に増加していることが分かる。

 

Euler-gosa-001.png

 

対して、λ=−1のとき、誤差は0≦x≦1では増加するが、さらに計算が進むに従い、誤差が減少してゆくことが分かる。

 

Euler-gosa-002.png

 

また、λh=−2となるように、λ=−2h=1とすると、次のような振動解が得られる。

 

Euler-gosa-003.png

 

そして、λ=−3、h=1とすると、この計算はやがて発散してしまう。

Euler-gosa-004.png


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

差分法による2階常微分方程式の境界値問題の解法 [数値解析]

差分法による2階常微分方程式の境界値問題の解法

 

次の2階常微分方程式

を、差分法を用いて解く解法について考える。

 

閉区間[a,b]a≦x≦b)をn等分し

とし、

とおき、さらに分点におけるyの値を

とおく。

 

(1)式の導関数を

と差分法を用いて近似すると、k=1,2,・・・,n−1において

というn−1本の方程式が得られる。

未知数は、n−1個なので、連立方程式(4)を解くことにより、これらの未知数を求めることができる。

 

連立方程式(4)は

とおくと

と書き換えることができる。

そして、(6)式の左辺の係数行列は、3重対角行列なので、TDMATri Diagonal Matrix Algorithm)を用いて解くことができる。

 

以下に、微分方程式

を解かせたプログラムと、分割数をn=10n=100とした場合の計算結果を示す。

 

! 差分法を用いて y''+f(x)y'+g(x)y=h(x) の境界値問題(でぃれくれ型)を解くプログラム
! ただし、固有値問題は解けない!!
parameter (ns=100,ns1=101)
real x(0:ns), y(0:ns)
data x,y/ns1*0.,ns1*0/

n=10

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

call Solver(x,y,n)

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

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

function f(x)
f=-14./x
end

function g(x)
g=x**3
end

function h(x)
h=2*x**3
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) ! 境界条件
    d(n1)=d(n1)-c(n1)*y(n) ! 境界条件

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

 

計算結果はおかしく見えるかもしれないけれど、この計算結果は正しいケロよ。

 

上の微分方程式の解は、ベッセル関数という特殊関数を用いないと表わせないうえに、ベッセル関数は組み込み関数にないので、n=100のときの数値解を厳密解だと考えて欲しいにゃ。

 

少し違うけれど、十進BASIC用のプログラムは
http://nekodamashi-math.blog.so-net.ne.jp/2016-01-24

に書いてあるので、Fortranを使えないヒト、あるいは、Fortranから他のコンピュータ言語にできないヒトは、コチラのプログラムを参考にしてほしいにゃ。

 


nice!(1)  コメント(0) 
前の10件 | - 数値解析 ブログトップ