対流項（移流項）を含む偏微分方程式の数値解法 [ネコ騙し数学]
対流項（移流項）を含む偏微分方程式の数値解法
問題 u(x,t)=f(x−ct)であるとき、
であることを示せ。ただし、cは定数とする。
【解】
ξ=x−ctとおき、合成関数の微分法を用いると、
したがって、
（解答終）
したがって、fを微分可能な任意の関数とすると、u(x,t)=f(x−ct)は
の（一般）解になる。
f(x)=sinxとおくと、
となり、これが（１）の（特殊）解になることから分かるように、（１）は波の方程式の一つである。
（１）式の左辺第１項を
と、時間微分に関しては前進差分を用いて近似し、（２）式の左辺の第２項を
と、空間微分に関しては中心差分を用いて近似すると、（１）式は
と近似することができ、これから
という漸化式を得ることができる。
とおくと、
となる。
（６）式の無次元数（物理の単位をもたない数のこと）をクーラン数という。
したがって、計算の起点となる、t=t₀におけるuの初期条件
とx=x₀における初期条件
が与えられていれば、（７）式を用い逐次的に（１）式の近似解
を求めることができる。
ということで、初期条件を
さらに、c=0.2、Δx=1、Δt=1として、（１）を解いた結果は次のようになる。
赤い曲線は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=2、x=3、x=4と下流にいけばいくほど、波の最大の高さは高くなり、より激しく振動するという、物理的にありえない、数値的に振動する解を（７）式は求めてしまう。
つまり、（１）式の時間微分に前進差分、空間部分に中心差分を用いるFTCS（Forward in Time Center in Space）は、数値的に不安定な解法で、この解法を（１）式に用いると物理的にありえない振動解を求めてしまうということ。
そこで、なぜ、FTCS法は不安定なのか、フォン・ノイマンの安定判定法を用いて調べることにする。
とし、（７）式に代入すると、
となり、
したがって、フォン・ノイマンの安定判定法によると、クーラン数に関わらず、FTCSは無条件不安定！！
ということで、FTCSは微分方程式（１）の数値解を求めることに使えないことが明らかになった。
そこで、（１）式の第２項の空間微分を
と風上化（上流化）すると、c>0のとき、（１）は
したがって、
となる。
（８）式を（１３）式に代入すると、
よって、
したがって、
の時に安定になる。
風上（上流）差分を用いて、c=0.2、Δx=1、Δt=1、すなわち、クーラン数C=0.2として計算した結果は、次の通り。
高さ１の波の到達時刻は実際とは異なるけれど、今度は振動することなく計算できていることが分かる。
また、下図を見ると、この波の解は
でなければならないのに、風上差分を用いた（１３）の数値解は、波の形が崩れていることが分かるだろう。
参考までに、同一の条件でFTCSを用いて計算した結果を下図に示す。
放物型偏微分方程式の解法 ルンゲ・クッタ法 [ネコ騙し数学]
放物型偏微分方程式の解法 ルンゲ・クッタ法
次の熱伝導方程式
の右辺を差分法を用いて次のよう
と近似し、
とおくと、ルンゲ・クッタ法を用いて、（１）式の近似解を求めることができる。
すなわち、
例によって、（２）を用いて、次の問題を解くと、右のグラフのようになる。
問題
を、初期条件
境界条件
のもとで、Δt=1、Δx=1として解け。
Δt=1、Δx=1という粗い計算でも、良好な計算結果が得られていることが分かる。
Δt=0.1、Δx=0.4とすると、さらに良好な結果が得られ、厳密解との差はほとんど認められなくなる。
なのですが、ルンゲ・クッタ法も陽解法の一種なので、
のとき数値解が振動することがあるので、ルンゲ・クッタ法を用いて熱伝導方程式を解くとき、
になるようにΔt、Δxを取ったほうがよい。
このことを示すために、Δt=0.36、Δx=0.5のとき、すなわち、
として計算すると、４０ステップのt=14.4の前後で振動を始め、t=18ではさらに激しく振動する。
計算に使用したプログラムは次の通り。
! ルンゲ・クッタ法による解法
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
［陽解法でも収束するかも］ [ネコ騙し数学]
［陽解法でも収束するかも］
ddt^3です。ちょっと面白いなぁ～と思って、差分法に関する安定条件をググってみると、フォン・ノイマンの安定条件ってのが出てきます（註１。
http://dr-asa.hatenablog.com/entry/2017/08/31/154806
上記HPによると、放物型偏微分方程式の陽解法が安定であるためには、任意の整数0≦mに対して、
だそうです。ここでGは時間方向への誤差の拡大率，βmは、
であり、Lは解析領域の長さです。この条件は解の大きさに誤差も比例するという前提なので、mは解を空間方向にフーリエ展開した場合の波数（モード次数）と考えられます。
к＝1/2，Δt＝0.1，Δx＝0.4，L＝4として、m＝0，1，2，・・・に対してG(m)をプロットしてみると、図-1になります。
波数20以上は周期的な繰り返しです。解析結果のグラフをみてみると、どう考えても1次モードのみ励起されるはずです。
です。
Gは誤差の拡大率です。初期誤差をεとし、誤差は解の大きさに比例するという前提を信じれば、
で、時間ステップnにおける誤差を評価できる事になります（たぶん最悪の予想）。T(x，t)は解，tnはnステップ目の時刻です。
解析結果は放物線が一様に圧縮されてくような感じなので、T(x，tn)／T(x，0)をx＝2の値の変化で代表させます。紙面から読み取った値は、
・・・表-1くらいかな？。
これをグラフに起こして強引に指数関数近似したのが図-2です。相関係数がR2＝1とウルトラ良いので(^^)、図-2の補間結果から、
で近似する事にします（註２）。
次に初期誤差εですが、T(x，0)＝x(4－x)＝－x2＋4xなので、t＝0でのx方向への2階差分の絶対値はa(0)＝1のはず。
時間差分は1階で、空間差分は2階だから、
くらいかな？(^^;)。
以上を使って、誤差の積算Ε(n)を計算してやると、図-3になります。
「ネコの旦那！。陽解法で永遠に数値積分しても、最大誤差0.0102くらいで収束しそうでっせ！」
陰解法は、長時間積分に対して陽解法より平均的に精度が良いだけの話で、局所的な精度は陽解法が上回る場合もあり得ると、自分は思っています。例えばあまり長くない時間積分であれば、無条件安定なシンプレクティック法よりルンゲクッタ法の方が、ずう～っとずう～っと精度が良い。
シンプレクティック法は陽解法のくせに（まさに陰に）陰解法の特性を持つのに対して、ルンゲクッタ法による超長時間積分では、解は0に減衰します。という訳で・・・。
予想1)
陽解法でも発散しないケースだったから、陽解法の精度がバカみたいに良かった？。
でもですね、n＝100で既に0.001～最大0.01程度の誤差は見込まれます。単精度ですよね？。問題のスケール規模は100＝1のオーダーで単精度の有効数字は6桁。倍精度でやりたいなぁ～(^^;)。
予想2)
倍精度計算したら、劇的に状況が変化するかも。あくまで「するかも」。
以上、他人のプログラムなんて頼まれても読みたくないddt^3の、他人事のような発言でした(^^;)。
（執筆：ddt³さん）
（註１）フォン・ノイマン法による陽解法の安定性判定については、ねこ騙し数学でも過去に「陽解法による拡散方程式の安定性」という記事で取り上げている。
http://nekodamashi-math.blog.so-net.ne.jp/2016-11-23
（註２）微分方程式の厳密解は、
そして、n=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
陽解法の時にも連立方程式を解いているので、陽解法のとき、少し無駄な計算をすることになるけれど、陽解法、陰解法、クランク・ニコルソン法といった解法の区別をすることなく計算できる。
また、αを一種の重みと考え、αを0、1/2、1だけではなく、0≦α≦1の任意のαを選び計算することもできる。
ただ、
にとると、
の値によって、伝播誤差のために、振動解が得られる。
例によって、（１）のκ=1/2とし、条件
とし、Δt=0.1、Δx=0.4に選んで計算した結果を以下に示す。
陽解法
陰解法
クランク・ニコルソン法
Δt、Δxをこの程度にとれば、陽解法、陰解法、クランク・ニコルソン法ともに、（１）の厳密解とよく一致することが分かる。
また、Δt=0.2、Δx=0.4とし、陽解法で解いた結果は次の通り。
このとき、
となるので、この条件で陽解法を用いて解こうとすると、伝播誤差のために、計算を進めるたびに、数値解の振動が激しくなり、厳密解からの逸脱が大きくなることが分かる。
放物型偏微分方程式の解法３ クランク・ニコルソン法 [ネコ騙し数学]
放物型偏微分方程式の解法３ クランク・ニコルソン法
微分方程式
を、陽解法を用いて離散化すると、
となる。
陰解法を用いて離散化すると、
（２）と（３）の平均を取ると、
になる。
とおくと、
になる。
ここで、
とおくと、
と３重対角行列を係数に持つ連立方程式が得られ、三重対角行列アルゴリズムTDMAを用いてこの連立方程式を解くことが可能になる。
大胆に、（６）式の
と
以外の項を落とすと、
となり、r>0のとき、
となり、クランク・ニコルソン法は無条件安定であることが予想できる。
例によって、次の問題をクランク・ニコルソン法で解いてみる。
問題
を、初期条件
境界条件
のもとで、Δt=1、Δx=1として解け。
Δt=1、Δx=1という粗い計算でも、良好な計算結果が得られていることが分かる。
Δt=0.5、Δx=0.4にすると、厳密解とクランク・ニコルソン法による近似解はほとんど等しくなる。
この計算に使用したプログラムは次の通り。
! 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
放物型偏微分方程式の数値解法２ [ネコ騙し数学]
放物型偏微分方程式の数値解法２
前回に引き続き、差分法を用いた熱伝導方程式の数値解法について述べることにする。
この方程式の時間微分を
右辺の空間微分を
と近似すると、（１）式は
となる。
とおくと、
という３重対角行列を係数を持つ連立方程式が得られる。
（１）の解をT=T(x,t)とし、初期条件を
境界条件を
とする。
そして、
とすると、初期条件より
境界条件より
となる。
したがって、
という連立方程式が得られ、逐次的に、３重対角行列アルゴリズム（TDMA）を用いて解くことによって、微分方程式（１）の近似解を求めることができる。
このような解法を陰解法という。
陰解法は、前回紹介した陽解法とは異なり、無条件安定である。
陽解法との比較参照のために、前回と同じ問題
問題
を、初期条件
境界条件
のもとで、Δx=1、Δt=1として解け。
を解かせたものは次の通り。
この問題の場合、厳密解と比較すると、陰解法による近似解は厳密解よりも大きいことがわかる。
さらに、前回の陽解法による結果を加え比較したものが次のグラフである。
込み入っていて、少しわかりづらいと思うのだけれど、厳密解は陰解法と陽解法の間に位置している。
ということで、陰解法と陽解法によって得られた数値解を足して２で割ると厳密解に近い結果が得られることが予想できる。
この予想は正しく、陽解法と陰解法の中間的な解法、クランク・ニコルソン法という解法が存在する。
Δt=1、Δx=1という非常に粗い計算でありながら、クランク・ニコルソン法は、厳密解とよく一致していることが分かる。
ただし、この事実をもって、陽解法や陰解法は精度が悪いと思わないで欲しい。
Δt=0.1程度にとれば、陽解法、陰解法ともに良好な結果を得ることができる。
陽解法はともかく、陰解法は熱や流れなどの数値計算ではよく使われている、優れた解法である。
問題の初期条件を
境界条件を
と変え、陰解法を用いてΔt=0.5、Δx=0.4として計算した結果を以下に示す。
このようにあらい格子でも、厳密解をよく再現していることが分かる。
そして、t→∞の定常解
に収束することが確認できる。
誤差は最大で１０％くらい。
Δt=0.1、Δx=0.4にすると、誤差は最大で約３％。
陰解法よりより高精度な計算ができる、クランク・ニコルソン法は、次回ということで。
計算に使用した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
放物型偏微分方程式の数値的解法１ [ネコ騙し数学]
放物型偏微分方程式の数値的解法１
放物型偏微分方程式の熱伝導方程式
を例に取り、差分法を用いた数値的な解法について考えることにする。
（１）の左辺の時間（偏）微分
右辺の２階の空間（偏）微分を
と、差分法を用いて近似すると、
ここで、
とおくと、
となる。
を計算する際、
はすべて既知なので、（５）式を用いて逐次的にを計算することができる。
とくに、r=1/2のとき、（５）式は
になる。
このような解法を陽解法という。
陽解法は代数方程式を解くことなく簡単に計算できるが、
のとき、（５）を用いて求めた近似解は安定ではなく、振動解が得られる。
したがって、
となるようにΔt、Δxを選ばないといけない。
（５）式の右辺を
と書き換え、右辺第２項を無視すると――なんと大胆な(^^ゞ――
誤差がこれにしたがって伝播するとする、安定であるためには、
でなければならい。
したがって、安定であるためには
・・・。
もう少し正確な議論をすると、誤差
が
に従うとする。
すると、
だから、少なくとも、安定であるためには
でなければならないに違いない。この条件を満たさないと、近似解は振動したり、発散するだろう。
そして、r>0だから、
正確な議論をするためには、von Neumannの判定法などを用いる必要があるが、それは厄介なので、簡易的にこの関係を求めてみたにゃ。
問題
を、初期条件
境界条件
のもとで、Δx=1として解け。
【解】
Δtを
となるように、Δt=1にとると、
となる。
以下、同様に計算すると、次の表が得られる。
この結果をグラフにすると、次のようになる。
定量的にはともかく、このような粗い計算であっても、指数関数的な現象を捉えており、定性的には正しい結果が得られていることがわかる。
より精度よく計算するために、Δx=0.1とすると、条件（７）より
よって、Δtを0.01以下にとる必要があり、計算量が大きく増えてしまう。
したがって、実際は、熱伝導方程式を陽解法を用いて解くことはない。
こうした制約のない陰解法を次回紹介することにする。
平行平板ポアズイユ流れの熱伝達のプログラム（最終版） [ネコ騙し数学]
! 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
陰的オイラー法 [ネコ騙し数学]
陰的オイラー法
微分方程式
の解は
で与えられる。
この式の右辺の定積分を
と近似し、
と逐次的に微分方程式の近似解
を求める方法を陰的オイラー法という。
等間隔の場合、すなわち、
のとき、
である。
このとき、陰的オイラー法は
となり、この漸化式を解くと
になる。
また、誤差
も漸化式に従うとすれば、
のとき、すなわち、
のとき、陰的オイラー法は安定した解法になる。
λ=3、h=1のとき、λh=3>2となり陰的オイラー法は安定であり、
に収束する。
しかし、この条件のときの微分方程式の解は
であり、陰的オイラー法による近似解は厳密解の挙動を表していない。陰的オイラー法による近似解が安定でかつ意味を持つのは、λh<０のときである。
数値計算でいう安定とは収束解が得られるという意味であり、その収束解が微分方程式の解であることを意味しないことに注意。
陰的オイラー法は、陽的オイラー法と同様に打切誤差がO(h²)程度で、しかも、方程式を解かなければならないので、実際に使われることはない。
たとえば、
の場合、
を解かなければならない。
この方程式は二分法やニュートン法を用いて数値的に解くことができるが、各段階でこの計算を行わなければならないので、計算量が多くなってしまう。理論的は話ではともかく、陰的オイラー法の打切誤差はO(h²)程度なので、とてもじゃないけれど、こんな計算は馬鹿らしくてやってられない。
参考までに、λ=1、h=0.1のときに、陽的オイラー法と陰的オイラー法を用いて（６）を数値的に解いた計算結果を示す。
この場合、陽的オイラー法の方が陰的オイラー法よりも精度よく計算できていることが分かる。
λ=−1、h=0.1のときは、次のようになる。この場合は、陰的オイラー法の方が精度よく計算できていることが分かるだろう。
また、一般的な傾向として、厳密解と比較した時、陰的オイラー法は大きめの値を出し、陽的オイラー法は小さめの値を出し、厳密解はこの２つの間に位置し、したがって、陽的、陰的オイラ法の近似解の平均をとれば厳密解に近い値を求められそうである。
そこで、λ=1、h=0.1の場合について計算してみると、次のようになる。
予想通り、劇的に精度が向上したにゃ(^^)
オイラー法の誤差の内訳 [ネコ騙し数学]
オイラー法の誤差の内訳
この記事を書いた私自身、このあたりの話はよく分かっていないのだけれど・・・。
微分方程式を
とする。
（１）式の点
における近似解、厳密解を
とするとき、で発生した誤差
は、オイラー法の場合、
次のステップの
の近似解
に伝播していく。
の場合、
だから、
である。
そして、この伝播誤差
の他に、k+1段の積分計算にともなう打切誤差
が加わり、
となり、k+1段に
として伝わってゆくのでしょう。
この考えに基づき、誤差の蓄積の内訳を求めてみた。
λ=1、h=0.1の場合は次の通り。
オイラー法の場合、伝播誤差をまったく含まないx=0.1を除き、誤差の大部分を伝播誤差が占めており、この伝播誤差のために、計算を進めるたびに厳密解との乖離が大きくなってゆくことが分かる。
安定であるλ=−1、h=0.1の場合は、次のようになる。
この場合、伝播誤差（の絶対値）はx=1.2の付近まで増加するのは、伝播誤差の減衰よりも計算のたびに加わる打切誤差の（絶対値）の増加が大きいため。
この場合も、誤差の多くを占めていることが分かる。
なお、この計算においては、打切誤差は、相対誤差と伝播誤差の差から求めてある。
ネムネコは数値計算の素人だから、こうした計算法の是非についてはわからない。そして、数値計算屋さんと伝播誤差を違う意味で捉えているのかもしれない。
だから、この方面に詳しいヒトや数値計算屋さんに、こうした計算の是非などについて教えて欲しい。
