波動方程式の数値解

質問

今、波動方程式を差分法の陰解法を用いて数値的に解いています。
しかし解いた結果はなぜか時間と共に減衰しています。
周期は厳密解とほぼ一致しています。
プログラム的には問題ないと思うのですが。
一応、使用した陰の差分方程式は以下のとおりです。

(u[n+1][j]-2*u[n][j]+u[n-1][j])/(dt^2)=(u[n+1][j+1]-2*u[n+1][j]+u[n+1][j-1])/(dx^2)
とし、ここでnは時間、jは距離の変化を表していてdt/dx=rと置きます。上式を変形して
r^2*u[n+1][j-1]-(1+2*r^2)u[n+1][j]+r^2*u[n+1][j+1]=u[n-1][j]-2*u[n][j]
この式を行列で表して連立一次方程式を解くことにより解を得ています。

なんらか原因が分かる方お願いいたします。

回答No.2

nta

>繰り返し計算することにより丸め誤差が大きくなって減衰してしまったのではないかと思うのですが、その可能性はあるでしょうか？

もちろんその可能性はあります。解析対象のスケールLと固有振動周波数fの積fL が１に近ければrの値は妥当なものが得られると思いますが、極端に小さいとか大きい場合には２乗の効果がによって丸め誤差はかなり大きなものになります。初期値を見る限りｆ＝１でL=1となっていますか。u'(x,0)=0に対応できるような初期化が正しく与えられているでしょうか。

お礼

ご回答ありがとうございます。
これといって丸め誤差の解決作は見つかりませんが、頑張ってみます。

回答No.1

nta

１次元モデルらしいことは分かりますが、初期条件・境界条件も励振源の性質も示されていませんから、減衰すればいいのか、振幅が増大するべきなのか判断ができません。境界条件のみを与えて固有値を求めるようにプログラムすれば過渡的なものを考える必要がないので煩わしくはないと思いますが、初期条件がモデルに会わないときには振動が分岐していくため挙動はかなり複雑だと思います。減衰というのはそういうことを意味しているのかもしれません。ただ、自作のプログラムに問題がないことを証明するのはかなりいろいろな問題を解いて実績を積んでからですから、なんどもなんどもプログラムコードを見直す必要があります。

補足

ご回答ありがとうございます。
質問の仕方がお粗末で申し訳ございません。
指摘のように上の差分方程式は1次元の弦の振動を表しています。
境界条件は両端でu(0,t)=u(l,t)=0,初期条件はu(x,0)=sin(πx),u'(x,0)=0で与えています。
全く同様な問題を陽解法で解いたところ、減衰せず発散もせずに一定の波形が現れました。私的には、繰り返し計算することにより丸め誤差が大きくなって減衰してしまったのではないかと思うのですが、その可能性はあるでしょうか？
よろしくお願いいたします。