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

[放物型偏微分方程式にシンプレクティック法は使えない!(普通にやろう(^^))] [数値解析]

[放物型偏微分方程式にシンプレクティック法は使えない!(普通にやろう(^^))]

 

 放物型にシンプレクティック法を使ってみたら、とんでもない目にあったので、今度は普通に重み付き残差法の弱形式によるガラーキン法をやってみます。支配方程式は、

                

 

 

 

です。重み関数をu(x)として、

                                  

 (2)中辺の2項目を弱形式に持ち込みます。

                       

 

 

 

従って(2)は、

                            

 

 

になります。ここでL1要素の長さ。(4)右辺1項目(弱形式)は、x1階微分の積分になってるので、線形近似くらいで十分です。要素端における温度値をT1T2とすれば、

                                

 

とあっさりと形状関数まで出てきます。

              

 

 

 

とすると、(4)右辺1項目の離散化近似は、

                   

 同様に(4)左辺の離散化近似は、

                                  

になります。

 (4)λ0とし、境界項は物理量の連続性から0になるとみなすと、(7)(8)から有限要素の一要素について、

                      

 

 

 

が得られます。左辺のConsistent Massは、恐らく実用的にはLamped Massで代用してやっても十分です。

                     

 

 

 

 いま解析領域全体をn等分割し、節点未知数をT1T2,・・・,Tn+1として、n個の要素全てについて(10)を重ね合わせた姿を想像すると、

                        

です。(11)1行目とn+1行目は、境界条件でおきかえ可です。そして2行目~n行目は、2階差分の定式化そのものじゃないですか!(^^)。こうして重み付き残差法の弱形式によるガラーキン法と、差分法がほぼ同等である事を導けました。ただ2次元以上の場合は、差分用の構造格子でなくても、有限要素用の非構造格子で同じ定式化が可能なので、重み付き残差法の弱形式ガラーキン法の方が若干便利かも知れません。

 

 もう一つわかる事は、熱源密度λがデルタ関数のa倍の場合以前は、к・∂T/∂xがデルタ関数の特異点xjで-aの不連続ジャンプを起こすとして、

                

 

 

 

と特異点両側の温度値から特異点での温度値を推定する差分法を用いましたが、定式化(4)によれば、

                  

 

 

になるので、xxjに対応する節点で(11)右辺にa/Lを加えるだけで良い事がわかります。ただしu(x)は要素端点で値1をとる重み関数(形状関数)だとします。実際にそうやって1次のオイラー法で計算してみると、(12)を用いた結果とほとんど変わらないものが出ます。その理由は、(5)による近似はC0級なので、要素境界で1階微分が不連続変化しても対処できるからです。

 というかもともと(5)の近似では要素境界で1階微分は不連続に求まります。しかし境界項は物理量の連続性から0になるとみなすという条件を暗に付けたので、不連続性を強制するような付加条件がない限り、不連続性を出来るだけ小さくするような解が決まる仕掛けになります。不連続性を強制する付加条件があれば、自然にそれに応じた不連続性が現れます。だから差分法でも(13)を、そのまま使えば良いんですよ(^^)

 やっぱり数値解法って、こうでなけりゃいけませんね(^^)

 

(執筆 ddt³さん)

 


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

[放物型偏微分方程式にシンプレクティック法は使えないか?(使えない!)] [数値解析]

[放物型偏微分方程式にシンプレクティック法は使えないか?(使えない!)]

 ddt^3です。具体的にシンプレクティック法プログラムで計算してみました。例題はネコ先生と同じです。そして結論は、使えない!となりました。以下はドツボにハマり悪戦苦闘したあげく、大敗した経緯です(^^;)

 まずシンプレクティック法を強引に使うために導いた時間に関する2階微分方程式、

                   

 

 

 

は、熱伝導方程式、

             

 

 

 

を時間tで微分して作ったものなので、(1)(2)に対する必要条件です。不用意な条件でやるとあらぬ解が出そうですが、(1)の時間積分手順を考えてみると、原理的には初期条件さえ決めればシンプレクティック法で以後の全ての挙動は決まるので、(1)の解が(2)の解と一致するためには、初期条件の設定に敏感そうだと想像できます。

 初期条件は0≦x≦4T(x0)=-x24x,境界条件はT(0t)T(04)0です。またλ0。式(1)は、重み付き残差法の弱形式を通じて有限要素法に持ち込み、節点未知数として図-1の節点温度Tと温度勾配wを取ったので((3)は形状関数,ξx/4)初期条件として、初期温度と温度勾配分布および初期温度速度と勾配速度分布を与える必要があります。それには熱伝導方程式(2)を用いるのが妥当です。一要素の節点は両側の隣接要素と共有されるので、式(3)C1級の3次近似になります。

 

 

 

 

 

 

 

 

 

ddt³-smp-fig-01.png

 

 初期温度は明らかにT(x0)=-x24xです。初期温度勾配は次のように決めました。解析領域の両端x04では温度一定なので(2)より、T(xt)xでの2階微分は0になる必要があります。この条件が成立するようにx04の傾きw1(0)w2(4)を決めます。中間節点位置xxjでは要素間でxでの2階微分が連続になるようにw(xj)を定めます。

 解析領域をn分割するとすると(n要素)、節点はn+1個あり、一節点2自由度なので未知数は2(n+1)個。そのうち初期温度分布からn+1個の温度値が決まり、両端での2階微分の値と中間での2階微分の連続性から、勾配値n+1個に対する条件も得られので、初期温度分布さえ与えれば温度分布と温度勾配分布を計算できるようになります。ただしこれは、2(n+1)個の節点未知数に関する連立一次方程式を解く事になります。結果はC2級の3次近似です(← すでにドツボの予感(^^;))。

 温度速度分布については(2)から、xでの2階微分がわかれば良い訳ですが、さっき2階微分が連続になるように温度勾配を与えたので、(3)ξで(xで)2階微分すれば温度速度分布が得られます。

 勾配速度分布については、(2)xで微分したものから得られるので、xでの3階微分です。これも(3)ξ3階微分すれば得られますが、(3)3次関数なので一般には不連続です。でも(2)には、xでの3階微分は直接的には出てこないので、こんなところでどうだろう?と始めました。





 図-1の横軸は時間,縦軸は解析領域中央(x2)における温度変化です。d-オイラー1次のラインは、通常の差分を行い1次のオイラー法で時間積分した結果で、ここでは厳密解とみなします。

 初期条件をC2級近似しシンプレクティック法に持ち込んだ結果を、s-離散-Non.Corで表します。御覧のようにいともあっさり発散します。t10では2.86×10204に達しました(^^;)。離散はLamped Massの意味。Consistent Massでもやってみましたが結果は変わりません。

 ここで気づいた事があります。式(1)λ0とすると、

            

 

 

 

ですが、これは正のフィードバック系です。

 自分が良く扱う梁の曲げ振動は、

          

 

 

 

という形をしてますが、これは負のフィードバックです。右辺にマイナス符号があるので、Tが大きくなり過ぎると逆向きの運動をさせようという加速度(復元力)が働らくので、安定した振動解になります。しかし(2)は正のフィードバック系なので、Tがちょっとでも大きくなり過ぎるとどこまでも加速して行きます(どんだけぇ~(^^))。つまり数値誤差にも非常に弱い、という事になります。

 

 という訳で修正項(安定化機構)が必要でないか?と考えました(←ドツボへの入り口)。そこで重み付き残差法のもともとの形に戻ると、

                               

 

でした。(4)右辺の23項目は物理量の連続性より0としました。でも今は初期条件を除きC1級近似なので、数値上は23項目も0ではありません。(3)などでは物理量の連続性から、これらを0とする強制条件を与えた方が上手く行くのですが、(2)が非常に数値儀差に対して敏感な系だとわかったので、これらの値も考慮すれば安定するのでは?と考えました。(4)右辺の23項目は、要素端でのいわば不釣り合い力を表すからです。それらを考慮するという事は、数学的には偏微分方程式をかなり人為的に操作する事にもなるのですが、系が超発散しやすいので、仕方ないだろうという訳です。

 

 最初に言いますがみなさん、こんな事は絶対に考えてはいけませんよ。数値解法が誤差に非常に弱いとわかった瞬間にあきらめるべきです。でもやってみたかったんだよなぁ~(^^;)。という訳で時間ステップの各段で(4)の境界項を考慮した結果が、図-1の赤ライン:s-離散-Cor.です。

 結果はある意味で安定しました。温度はt0から時間に比例して下がり続け、t10では-6.0382.86×10204と比較すれば大進歩です。・・・やっぱり不釣り合い力が影響してるのだ。という事は不釣り合い力が出ないようにしてやれば・・・と、ますますドツボへ。

 初期条件では温度分布からC2級近似してるから、(4)右辺の2階微分項は0になります。だったら時間ステップの各段で得られた温度分布から、同じように温度勾配,温度速度,勾配速度をC2級近似になるように決め直したら?。結果は水色点線のs-離散-Cor.-C2です。オ~~、s-離散-Cor.と全く同じ結果だ。やはり不釣り合い力だ!。

 (4)右辺の3階微分項も0になるようにC3級近似しよう。それには最低4次近似必要だが、(3)の形状関数に4階微分値の定数項を未知数として含めるだけだから、なんとかなるだろう。温度分布は事前に与えられるとして、温度勾配は最初と同じように2階微分の両端での値と要素間での連続性から決め、3階微分も要素間での連続性から決める。・・・とやってみると、なんと条件が1個足りなくなりました。3階微分の値の1個が不定になるのです。普通ならこの辺りでやめるのですが・・・(^^;)

 

 ところがドツボ状態だと、都合の良い方へ都合の良い方へと確証バイアスが働きます。3階微分の値の1個が不定にならならそれは自由に選べるから、最適近似を選べるという事じゃないか。これは良い事だ!。最初に従来の方法でC2級近似を行い、得られた温度勾配分布との差が最も小さくなるように、最小2乗法で3階微分の値を選んでやればどうだろう?。やってみると、なかなか綺麗なC3級近似が得られます。その計算結果が、s-離散-C3です。t105.23×1096。最初に戻っちゃった・・・(^^)

 

 C3級近似すると、(4)の境界項はいつも0です。予想ですが、少なくとも空間方向の近似が厳密解と同等程度でなければ、この方法は上手く行かないのでしょう。1次元ならば、それは空間方向の境界要素法を併用する事でほぼ可能ですが、仮に上手く行ったとしても、非常に不安定な系を解いてる事に変わりはありません。ちょっと激しい変動をする熱源密度なんかを付けただけで発散しない保証は0です。また2次元以上では境界要素法でも、空間近似の影響が入ります。ここで万策尽きました。

 これだけ手をかけてもまともな応えが得られないという事は、放物型へのシンプレクティック法の使用は、明らかに失敗です。理論と数値計算は違うのだ、と言い訳しておきます。

 

 撤収ゅ~~!(^^;)

 


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

差分法による2次元非定常熱伝導方程式の解法3 クランク・ニコルソン法 [数値解析]

差分法による2次元非定常熱伝導方程式の解法3 クランク・ニコルソン法

 

2次元非定常熱伝導方程式

に対する陽解法の漸化式は

であり、純陰解法による漸化式は

である。

(2)、(3)式は重み係数αを用いると、

と表すことができる。

陽解法はα=0であり、純陰解法はα=1の場合である。そして、(4)式の重み係数α=1/2とした

ものがクランク・ニコルソン法(半陰解法)である。

陽解法(2)、陰解法(3)の打切誤差がであるのに対し、クランク・ニコルソン法の打切誤差はなので、Δx,ΔyΔtが同一の場合、陽解法、陰解法よりも精度のよい計算結果が得られ、また、陰解法と同様に無条件安定である。

 

以下に、

をクランク・ニコルソン法で解くプログラムを示す。

 

! 2次元否定情熱伝導方程式の解法 (クランク・ニコルソン法)
parameter (nx=10,ny=10,nt=10)
real t(0:nx,0:ny,0:nt)
real kappa

alpha=0.5  !α=0のとき陽解法、α=1/2のときクランク・ニコルソン法,α=1のとき純陰解法

limit = 100
omega=1.0

dx=1./nx; dx2=dx*dx
dy=1./ny; dy2=dy*dy
dt=0.1

kappa=0.025

pi=acos(-1.0)

t=0.

! initial condition
do i=0,nx
    do j=0,ny
        t(i,j,0)=sin(pi*i*dx)*sin(pi*j*dy)
    end do
end do

! caliculate coeficients
aw=alpha*kappa/dx2
ae=alpha*kappa/dx2
as=alpha*kappa/dy2
an=alpha*kappa/dy2
ap0=1./dt
ap=aw+ae+as+an+ap0
aw0=(1.0-alpha)*kappa/dx2
ae0=(1.0-alpha)*kappa/dx2
as0=(1.0-alpha)*kappa/dy2
an0=(1.0-alpha)*kappa/dy2


do k=1,nt
    do iter=1,limit
        r_max=0.0
        do i=1,nx-1
            do j=1,ny-1
                s=aw0*t(i-1,j,k-1)+ae0*t(i+1,j,k-1)+as0*t(i,j-1,k-1)+an0*t(i,j+1,k-1) - (aw0+ae0+as0+an0)*t(i,j,k-1)
                res = aw*t(i-1,j,k)+ae*t(i+1,j,k)+as*t(i,j-1,k)+an*t(i,j+1,k) + ap0*t(i,j,k-1) -ap*t(i,j,k) + s
                cor=omega*res/ap
                r_max=amax1(r_max,abs(cor))
                t(i,j,k)=t(i,j,k)+cor
            end do
        end do
        if (r_max.lt.1.0e-6) exit
    end do
    write(*,*) k, iter, r_max
end do

do k=0,nt
    write(*,*) 'time=',k*dt
    do j=0,ny,1
        write(*,100) (t(i,j,k),i=0,nx,1)
    end do
    write(*,*)
end do

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

end

 

なお、このプログラムは、(5)式に基づくものではなく、(4)式をもとにコーディングしているので、重み係数α01/21にすることによって、陽解法、クランク・ニコルソン法、純陰解法のそれになるようにしてある。

 

計算結果

 


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

差分法による2次元否定情熱伝導方程式の解法2 純陰解法 [数値解析]

差分法による2次元否定情熱伝導方程式の解法2 純陰解法

 

2次元費定常熱伝導方程式

  

に、次の差分法を用いると、

  

次の差分方程式を得る。

  

したがって、この連立(差分)方程式を解くことによって、(1)の数値的な近似解を求めることができる。

このように(1)の近似解を求める方法を純陰解法という。

前回紹介した陽解法による2次元非定常熱伝導方程式の解法の安定条件は

  

であったが、純陰解法は無条件安定であり、空間分割ΔxΔy、時間分割Δtを自由に選ぶことができる。

また、(3)式の係数を見ると、κ>0Δt>0なので、

  

と優対角条件を満たしており、連立(差分)方程式は、ガウス・ザイデル法やSOR法といった反復法を解くことができる。

 

純陰解法を用いて

  

を用いて解くプログラムを以下に示す。

 

parameter (nx=10,ny=10,nt=10)
real t(0:nx,0:ny,0:nt)
real kappa

limit = 100
omega=1.3

dx=1./nx; dx2=dx*dx
dy=1./ny; dy2=dy*dy
dt=0.1

kappa=0.025

pi=acos(-1.0)

t=0.

! initial condition
do i=0,nx
    do j=0,ny
        t(i,j,0)=sin(pi*i*dx)*sin(pi*j*dy)
    end do
end do

! caliculate coeficients
aw=kappa/dx2
ae=kappa/dx2
as=kappa/dy2
an=kappa/dy2
ap0=1./dt
ap=aw+ae+as+an+ap0

do k=1,nt
    do iter=1,limit
        r_max=0.0
        do i=1,nx-1
            do j=1,ny-1
                res = aw*t(i-1,j,k)+ae*t(i+1,j,k)+as*t(i,j-1,k)+an*t(i,j+1,k)+ap0*t(i,j,k-1)-ap*t(i,j,k)
                cor=omega*res/ap
                r_max=amax1(r_max,abs(cor))
                t(i,j,k)=t(i,j,k)+cor
            end do
        end do
        if (r_max.lt.1.0e-6) exit
    end do
    write(*,*) k, iter, r_max
end do

do k=0,nt
    write(*,*) 'time=',k*dt
    do j=0,ny,1
        write(*,100) (t(i,j,k),i=0,nx,1)
    end do
    write(*,*)
end do

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

end

 

計算結果は次の通り。

 





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

差分法による2次元非定常熱伝導方程式の解法1 陽解法 [数値解析]

差分法による2次元非定常熱伝導方程式の解法1 陽解法

 

差分法による2次元非定常熱伝導方程式

  

の数値解法について考えることにする。

  

だから、(1)式は

  

となり、

  

のとき、(3)式を用いて前進的に安定して解くことができる。

Δx=Δy=hのとき、(3)式は

  

となる。

そこで、

  

となるようにΔthを選ぶと、(3)式は次のようになる。

  

 

ところで、ラプラス方程式

  

の数値解法で

  

という漸化式が出たが、(7)式はこれと同一の形をしていることが分かる。

そして、(7)式を用いてタイムステップを1つ進めることと、連立方程式(8)を反復法の一つであるヤコビ法

  

を用いて解くときの反復回数が1回増えることとは同じことになる。

 ――(9)のと、uの右肩に付いている記号(n)は反復計算の反復回数を表す――

 

 

問題 (1)の初期条件、境界条件が

  

であるとき、Δx=0.1Δy=0.1κ=0.025とし、(1)の数値解を求めよ。

【解】

h=Δx=Δy=0.1とおくと、

  

したがって、(7)式を用いて前進的に解くことができる。

 

計算結果

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

 

parameter (nx=10,ny=10,nt=10)
real t(0:nx,0:ny,0:nt)

pi=acos(-1.0)
dx=0.1; dy=0.1; dt=0.1

t=0.0
! initial condition
do i=0,nx
    do j=0,ny
        t(i,j,0)=sin(pi*dx*i)*sin(pi*dy*j)
    end do
end do

do k=1, nt
    do i=1,nx-1
        do j=1,ny-1
            t(i,j,k)=0.25*(t(i-1,j,k-1)+t(i+1,j,k-1)+t(i,j-1,k-1)+t(i,j+1,k-1))
        end do
    end do
end do

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

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

end

 

 

(3)式を用いるより汎用的なプログラムを作ることも可能であるが、安定的に解くことのできる条件が

  

と厳しいので、これはひとまず保留だにゃ。

 

 


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

ラプラス方程式の数値解法2 [数値解析]

ラプラス方程式の数値解法2

 

ラプラス方程式

  

に差分法を適用すると、

  

という連立差分方程式が得られる。

前回、この連立方程式をガウス・ザイデル法とSOR法といった反復法を用いて解いた。しかし、ガウス・ザイデル法、SOR法ともに、1回の反復で計算の影響が隣接する格子にしか及ばないので、収束が遅く、収束解を得るまでに多くの反復計算を要する。そこで、この反復計算を減らす方法について、今回、考えることにする。

Lapgrid.pngその前に、Patankharの流儀にならって(3)式を

  

と書き換えることにする。

ここで、

  

であり、

  

である。

WWestEEastSSouthNNorthを省略したものくらいの意味。

 

ガウス・ザイデル反復法は、計算領域の各点に対して

  

と計算する方法である。このため、反復計算ごとの、の更新値には、これに隣接するの4点の値しか影響を与えない。このため、ガウス・ザイデル法の収束は遅くなる。

そこで、

  

と変形する。すると、左辺は3項係数行列の方程式形となり、

  

とすることによって、TDMA(トマスのアルゴリズム)を用いてj列の(の推測値)を一気に計算することができる。

j列のすべてについてこのように計算することによって、1回の反復計算の収束速度を速め、反復回数を減らすことができる。

 

一方、SOR法は計算領域の各点に対して

  

と計算する。

ここで、ωは収束を早めるための加速係数で1<ωであり、は反復前のの値である。

(8)式は

  

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

したがって、(7)式の係数を

  

とすることによって、さらに収束を早めることができる。

 

 

parameter (nx=100,ny=100)
real u(0:nx,0:ny)
real aw(nx,ny),ae(nx,ny),as(nx,ny),an(nx,ny),ap(nx,ny)
real a(nx),b(nx),c(nx),d(nx)

u=0.
a=0.; b=0.; c=0.; d=0;

aw=0. ; ae=0. ; as=0. ; an=0. ; ap=0.

pi=acos(-1.0)
dx=1./nx; dx2=dx*dx
dy=1./ny; dy2=dy*dy

limit=400
alp=1.3

! 境界条件 (y=1)
do i=0,nx
    u(i,ny)=sin(pi*i*dx)
end do

! 係数の計算
do i=1,nx-1
    do j=1,ny-1
        aw(i,j)=1./dx2
        ae(i,j)=1./dx2
        as(i,j)=1./dy2
        an(i,j)=1./dy2
        ap(i,j)=aw(i,j)+ae(i,j)+as(i,j)+an(i,j)
    end do
end do

do iter=1, limit
    e_max=0.0
    do j=1,ny-1
        do i=1,nx-1
            a(i)=-aw(i,j)
            b(i)=ap(i,j)/alp
            c(i)=-ae(i,j)
            d(i)=as(i,j)*u(i,j-1)+an(i,j)*u(i,j+1)+(1.-alp)*ap(i,j)/alp*u(i,j)
        end do
        d(1)=d(1)+aw(1,j)*u(0,j)
        d(nx-1)=d(nx-1)+ae(nx-1,j)*u(nx,j)
       
        call tdma(a,b,c,d,nx-1)
       
        do i=1,nx-1
            cor=u(i,j)-d(i)
            e_max=amax1(e_max,abs(cor))
            u(i,j)=d(i)
        end do
    end do
   
    if (e_max.lt.1.e-6) exit
end do

write(*,*) iter,e_max

do j=0,ny,10
    write(*,100) (u(i,j),i=0,nx,10)
end do

open(1,file='Laplace.dat', status = 'replace')
do i=0,nx
    do j=0,ny
        write(1,110) i*dx, j*dy, u(i,j)
    end do
    write(1,*)
end do
close(1)

100 format (10(f8.5,1x),f8.5)
110 format(2(f8.5,1x),f8.5)

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

 


 

  

を、xy方向ともに100分割し、Δx=Δy=0.01とし、ω=1.3、さらに収束条件をとした場合、SOR法の反復回数は3376回であるのに対し、直接法TDMASOR法を組み合わせた計算法の反復回数は359回と反復回数を約1/10に抑えることができる。

 

PCの計算速度がべらぼうに速くなった現在でも、これくらいの計算量になると、計算速度の差を十分に体感できる!!

ただし、計算格子の数が50×50以下ならば、反復回数はSOR法の方が多いけれど、サブルーティンを使っていないのでむしろSOR法の方が計算時間は短いのではないだろうか(^^

 

今回紹介した計算法よりもさらに収束を早めるADI法が存在するけれど、基本的な考え方は今回紹介した方法と同様であり、プログラムの複雑さの割にたいして速くなるわけでもないので、割愛するにゃ。

 

 


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

楕円形方程式(ラプラス方程式)の数値解法 [数値解析]

楕円形方程式(ラプラス方程式)の数値解法

 

ラプラス方程式

  

の差分法を用いた数値解法を考える。

  

だから、(1)は

  

と近似できる。

(3)式は優対角条件

  

を満たしており、連立方程式(3)を反復法(ガウス・ザイデル法)を用いて解くことができる。

 

特にΔx=Δyのとき(2)式は

  

 

  

を、Δx=Δy=0.1として、(4)を用いて数値的に解け。

 

計算結果

 


 

計算に使用したプログラム

 

parameter (nx=10,ny=10)
real u(0:nx,0:ny)

u=0.0
limit = 1000
pi=acos(-1.0)

dx=1./nx
dy=1./ny

! 境界条件 (y=1)
do i=0,nx
    u(i,ny)=sin(pi*i*dx)
end do

do iter=1, limit
    res_max=0.
    do i = 1, nx-1
        do j=1, ny-1
            res=u(i,j)-0.25*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)) ! 残差
            res_max=amax1(abs(res),res_max)
            u(i,j)= 0.25*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))
        end do
    end do
    if (res_max.lt.1.e-6) exit ! 収束判定
end do

do j=0,ny
    write(*,100) (u(i,j),i=0,nx)
end do

100 format (10(f8.5,1x),f8.5)
end


 

Δx≠Δyの場合でも計算できる、より汎用的なプログラムは次の通り。なお、このプログラムでは、反復計算の収束を早めるためにSOR法を用いている。

 

parameter (nx=10,ny=10)
real u(0:nx,0:ny)
real aw(nx,ny),ae(nx,ny),as(nx,ny),an(nx,ny),ap(nx,ny)

u=0.
aw=0. ; ae=0. ; as=0. ; an=0. ; ap=0.

pi=acos(-1.0)
dx=1./nx; dx2=dx*dx
dy=1./ny; dy2=dy*dy
limit=1000

omega=1.3

! 境界条件 (y=1)
do i=0,nx
    u(i,ny)=sin(pi*i*dx)
end do

! 係数の計算
do i=1,nx-1
    do j=1,ny-1
        aw(i,j)=1./dx2
        ae(i,j)=1./dx2
        as(i,j)=1./dy2
        an(i,j)=1./dy2
        ap(i,j)=aw(i,j)+ae(i,j)+as(i,j)+an(i,j)
    end do
end do

do iter=1, limit
    res_max=0.0
    do i=1,nx-1
        do j=1,ny-1
            sum = aw(i,j)*u(i-1,j)+ae(i,j)*u(i+1,j)+as(i,j)*u(i,j-1)+an(i,j)*u(i,j+1)
            res = sum-ap(i,j)*u(i,j)
            res_max=amax1(res_max,abs(res)/ap(i,j))
            u(i,j)=u(i,j)+omega*res/ap(i,j)
        end do
    end do
    if (res_max.lt.1.e-6) exit
end do

do j=0,ny,1
    write(*,100) (u(i,j),i=0,nx,1)
end do

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

end

 


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

差分法を用いた、波動方程式の解法 陽解法 [数値解析]

差分法を用いた、波動方程式の解法 陽解法

 

差分法を用いた、波動方程式

  

の数値的解法について考える。

  

だから、(1)式は

  

ここで、

  

とおき、について解くと

  

あるいは、

  

となる。

 

  

となるようにΔxΔtをとると、(6)式は、さらに

  

と、簡単な形になる。

 

また、初期条件より、

  

となり、が与えられる。

t=0のときも(7)式が成り立つとすると仮定し、(8)式と初期条件を(7)に代入すると、

  

となり、(9)を計算の出発点とし、以後、(7)を用いて、前進的に解くことができる。

 

 

問題 境界条件と初期条件が与えられている。c=2の場合の(1)の解を、Δx=1、Δt=1/2として計算せよ。

【解】

このとき、C=cΔt/Δx=1となるので、(7)式と(9)式を使うことができる。

したがって、

 

(解答終)

 

この計算結果から、問題の解が周期が4であることがわかる。(下図参照)

 

 

計算プログラムを以下に示す。

 

parameter(nx=40,nt=200)
real u(0:nx,0:nt)
real l

u = 0.

l=4.
cu=2.
dx=l/nx
dt=0.05

CC = cu*dt/dx

! initial condition
do i =1,nx
    u(i,0)=f(i*dx,l)
end do
   
! boundary condtion
do j=1,nt
    u(0,j)=0.
    u(nx,j)=0.
end do

c1 = cc*cc
c2 = 2*(1.-cc*cc)
c3 = cc*cc

r= cc*cc

do j=1,nt
    do i=1,nx-1
        if (j.ne.1) then
            u(i,j)=r*(u(i-1,j-1)-2.*u(i,j-1)+u(i+1,j-1))+2*u(i,j-1)-u(i,j-2)
        else
            u(i,j)=0.5*(u(i-1,j-1)+u(i+1,j-1))
        end if
    end do
end do

write(*,*) '*** result***'
do j=0,nt,20
    write(6,100) j*dt, (u(i,j),i=0,nx,10)
end do

write(*,*)
write(*,*) '*** exact ***'

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

write(*,*)
       
100 format(11(f8.5,1x),f8.5)

end

function f(x,l)
real l

f=x*(l-x)
end

function exact(x,t)
real l

pi = acos(-1.0)
c=2.
l=4.

sum = 0.
do n=1,30
    sum = sum+1./(2*n-1)**3*cos((2*n-1)*pi*c*t/l)*sin((2*n-1)*pi*x/l)
end do

exact = 8*l*l/pi**3*sum
end

 

 

計算結果

 

 

 

 


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

バーガース方程式の数値解法 計算結果の比較 [数値解析]

バーガース方程式の数値解法 計算結果の比較

 

 

ν=0.01Δx=0.025Δt=0.025で、風上差分(UpWind Scheme)、中心差分(Central Difference Scheme)、ハイブリッド法(Hybrid Method)、そして、有限体積法(Finite Volume Method)を用いて解いた計算結果(t=0.5)を以下に示す。

 

 

 

 

 

定性的にはともかく、計算スキームによって計算結果がかなり異なることが分かると思う。

 

 


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

バーガース方程式の数値解法8 有限体積法(Finite Volume Method) [数値解析]

バーガース方程式の数値解法8 有限体積法(Finite Volume Method

 

バーガース方程式

  

  

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

 

Grid-002.png

 

そこで、(2)式を検査体積(Control Volume)で

  

と積分すると、

  

となる。

(3)式の各項を

  

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

  

と近似できる。

は、それぞれ、の境界面の速度で未知数であるが、

  

と近似することにより、(4)式は2次の連立方程式となり、近似的に(1)式を解くことができる。

なお、ここで、max{a,b}abの小さくない数を表す。

 

以下に、この方法を用いて

   

を解くプログラムを示す。

 

parameter (nx=20, nt=40)
real u(0:nx,0:nt)
real a(nx),b(nx),c(nx),d(nx)
real nyu

pi=acos(-1.0)
a=0. ; b=0. ; c=0.; d=0
u=0.

nyu=0.01
dx=1./nx; dt =1./nt
dx2=dx*dx

! initial condition
do i=0,nx
    u(i,0)=sin(2*pi*i*dx)
end do
u(nx,0)=0.

! boudary condition

! caliculate coefficients
do j=1,nt
    do iter=1,100
        do i=1, nx-1
            u(i,j)=u(i,j-1)
        end do
        do i=1,nx-1
            fe=0.5*(u(i+1,j) + u(i,j))
            fw=0.5*(u(i,j)+u(i-1,j))
            aw=nyu/dx2+amax1(fw,0.0)/(2.*dx)
            ae=nyu/dx2+amax1(-fe,0.0)/(2.*dx)
            a(i) = -aw
             b(i)=2.*nyu/dx2+amax1(-fw,0.0)/(2.*dx)+amax1(fe,0.0)/(2*dx)+1/dt
            c(i)=-ae
            d(i)=u(i,j-1)/dt
        end do
        d(1)=d(1)-a(1)*u(0,j)
        d(nx-1)=d(nx-1)-c(nx-1)*u(nx,j)
        call tdma(a,b,c,d,nx-1)
        do i=1,nx-1
            u(i,j)=d(i)
        end do
    end do
end do

do j=0,nt,4
    write(6,100) j*dt, (u(i,j),i=0,nx,2)
end do

write(*,*)

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

 

計算結果を見ると、この方法では、対流項を1次精度の上流差分で近似しているので、この方法の打切誤差は、陰解法と同程度のと考えられるが、ΔxΔtを粗くとっても、次のように精度よく計算ができていることが分かる。

――実線は厳密解(もどき)――

 

FVM-graph-002.png

FVM-graph-001.png

 

 

 

ネムネコ、一押しの計算法である。

 

 


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