「い」さんからいただいた質問の回答

「い」さんから頂いた質問

 

質問なのですが ルンゲクッタのdo k=2,3のループ内ではdk（1,0）やdk(1,nx)などを使っているように見えるのですが

dk(1,i)のループを見る限るそれがない気がしています．私の勘違いかもしれないですが教えていただけないでしょうか？
https://nekodamashi-math.blog.so-net.ne.jp/2017-11-24-3

 

【回答】

ルンゲ・クッタ法を用いる計算式（漸化式）は、

　　

です。

なお、ここで、

　　

です

 

プログラム中では、はdk(i,j)で、に対応するdk(1,j)です。

したがって、

　　

となるので、

プログラムでは

　　dk(1,i) = dt*f(t(i-1,j-1),t(i,j-1),t(i+1,j-1),c)

です。

 

dk(2,j)、dk(3,j)、dk(4,j)の計算ではそれぞれ、一つ前のステップの計算結果を必要としますが、

dk(1,j)の計算に関しては、

初期条件の温度と、

ルンゲ・クッタ法を用いてその都度更新される温度

　　

を使い、

　　

と計算します。

ですから、直接的ではなく、（３）式を通じて間接的に入っているというわけです。

（回答終）

 

ルンゲ・クッタ法を用いても解けますが、記事中に書いてあるように、ルンゲ・クッタを用いた解法はオイラー法などと同じく陽解法なので、数値解を振動させないためには、安定条件を満足させる必要があり、２次元、３次元の計算の場合、この条件を満足させるのが厳しくて、偏微分方程式の初期値問題の数値解法でルンゲ・クッタ法が使われることはありません。
しかも、
１次元の場合でさえ、このようにプログラムが長く、複雑になりますし、２次元、３次元になれば、さらにプログラムは複雑化し、その上、メモリーの容量も大きくなりますので、余程の物好きでなければ、ルンゲ・クッタ法を用いて解いたりしません。
高精度で計算したいならば、無条件安定のクランク・ニコルソン法を使った方がいいんですよ。

 

 

第２種境界条件の入れ方

 

断熱条件

　　

のプログラムへの実装の仕方。

という仮想の値を設け、条件（１）を中心差分を用いて表現すると、

　　

したがって、

　　

とすればよい。

 

より一般的の第２種境界条件

　　

の場合は、

　　

なので、

　　

 

そして、i=0を含めて、ルンゲ・クッタ法で解けばよい。

 

 

コメント 3

質問なのですが ルンゲクッタのdo k=2,3のループ内ではdk（1,0）やdk(1,nx)などを使っているように見えるのですが
dk(1,i)のループを見る限るそれがない気がしています.
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
ここにdk（1,0）やdk(1,nx)は定義されていない様子です。 例えば、dk(2,1)=dt*f(t(0,0)+dk(1,0)/2,t(1,0)+dk(1,1)/2,t(2,0)+dk(1,2)/2
dk（1,0）の与えはわかりません。
私の勘違いかもしれないですが教えていただけないでしょうか？

by PAN (2019-12-05 17:21) 

コメント、ありがとうございます。

☆ここにdk（1,0）やdk(1,nx)は定義されていない様子です。 例えば、dk(2,1)=dt*f(t(0,0)+dk(1,0)/2,t(1,0)+dk(1,1)/2,t(2,0)+dk(1,2)/2
dk（1,0）の与えはわかりません。
私の勘違いかもしれないですが教えていただけないでしょうか？
◇これは、
　初期条件：T(x,0)=x(4−x)
　境界条件:T(0,t)=0, T(4,t)=0
から、
　T(0,0)=0, T(nx,0)=0
　T(0,j)=0, T(nx,j)=0
になるためです。

この部分は、
プログラム中のT(i,j)の初期化、初期条件、ならびに、境界条件の箇所で
! 初期化
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
として与えています。

ルンゲ・クッタ法で実際に計算する部分は、i=1〜nx−1、j=1〜ntで、
t(0,j)=0、t(nx,j)=0
と、この部分は不変なので、計算をする必要がないんですよ。

境界条件が
　∂T/∂x(0,t)=0, ∂T/∂x(4,t)=0
といった（断熱）条件や
　∂T/∂x(0,t)=q₀, ∂T/∂x(4,t)=q₄
といった条件が入っている場合は、時間と共にT(0,t)、T(4,t)の値が変化するので、t(0,j)、t(nx,j)の部分も含めて（ルンゲ・クッタ法で）計算しないといけませんけどね。


by nemurineko (2019-12-05 20:24) 

なるほど！分かりました、ありがとうございます。
by PAN (2019-12-05 23:12)