最も簡単な波動方程式ー伝達方程式を解いてみましょう.

\[\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0\]

これは, 曲線\[\frac{\text{d}x(t)}{\text{d}t}=c\]つまり\[x(t)=ct+x_0\]に沿っての方向微分, つまり関数値$u(t,x(t))$の変化を考えると,

\[\frac{\text{d}u(t,x(t))}{\text{d}t}=\frac{\partial u}{\partial t}(t,x(t))+\frac{\partial u}{\partial x}(t,x(t))\frac{\text{d}x(t)}{\text{d}t}=0\]となることから, 解は\[u(t,x(t))=\text{Const.}\]となることがわかります. この曲線(今は直線か・・・)を特性曲線と呼びます. 例えば$t=0$で初期値$f(x)$が与えられているのであれば, \[u(0,x(0))=u(0,x_0)=f(x_0)\]であり, これが時間によらないので\[u(t,x(t))=f(x_0)=f(x(t)-ct)\]が成り立ちます．$xt$平面を特性曲線が埋め尽くしているので, $t,x$を独立変数と考えると, \[u(t,x)=f(x-ct)\]が解です.

陽的差分法

差分法でこれを解いてみましょう．空間$x$及び時間$t$を, 整数$n,i$を用いて以下のように離散化します:\[t_n=n\Delta t,\quad x_i=i\Delta x\]. 差分法では, この格子上の値\[ u_i^{(n)}=u(t_n,x_i)\]の変化について差分式を考え, それを解いて解を求めます．

差分法で用いる格子点の選択は重要です．

図の右上の3点を用いると, 時間微分と空間微分が差分近似できる:\[\frac{\partial u}{\partial t}\sim\frac{u_N-u_+}{\Delta t}\]\[\frac{\partial u}{\partial x}\sim\frac{u_+-u_-}{\Delta x}\]これを用いると，差分式\[\frac{u_N-u_+}{\Delta t}+c\frac{u_+-u_-}{\Delta x}=0\]が得られるので, これを解くと\[u_N=(1-K)u_+ +K u_-,\quad K=\frac{c\Delta t}{\Delta x}\]

これを用いると, $n=1$のデータを図の1,2,...の順に求めていくことができます. $n=1$のデータがわかれば, $n=2$がわかる・・・

こうして，時間発展を追うことができます. $u_N$が全て既知の値で表されている点に注目して，陽的差分法と呼びます.

時間発展を追う問題では, 初心者の想像とは異なり, 時間$t$の方向に配列を用意することはありません．1層づつ解いていくので, 例えば$n=10$を解くときには, 既に$n=1,2,\ldots, 8$くらいまでは全く関係がなくなっています．もちろん時間方向の配列を用意できるなら用意すれば良いのですが，それでは，時間配列を用いない同業他者様は, 同じメモリーサイズを用いて，ずっと細かい空間格子で計算してしまう(つまり競争に負ける)からです．

だいたいは，こんな感じのプログラムになるでしょう:

#include <string>

main():
...前略...
CONV mySolver(nx,dt,"pde_10.data");//空間nx分割, 時間ステップや出力ファイル名を指定
mySolver.Write();                  //初期値も保存しよう
for(int nt=1;nt<=ntmax;nt++)       //ntmax回時間発展
{
  mySolver.Step();                      //dt進む
  if (nt%ntsave==0) mySolver.Write();   //たまに保存する
}

conv.hpp:
...前略...
std::vector<double> U_VAL;
std::vector<double> U_OLD;
double dt;
...略...
CONV(int N,double dt,std::string filename):dt(dt)
{
  ...略...
  dx=xmax/N;
  U_VAL.resize(N+1);           //格子点は0,1,2,...NまでのN+1個！
  U_OLD.resize(U_VAL.size());  //Nを使うと間違えるから
  CFL=c*dt/dx;
  ...略...
  ofs.open(filename);
}
void STEP()
{
  for(int i=0;i<U_VAL.size();i++) U_OLD[i]=U_VAL[i];
  for(int i=1;i<U_VAL.size();i++)
    U_VAL[i]=(1.-CFL)*U_OLD[i]+CFL*U_OLD[i-1];
  time+=dt;
}

std::string filenameというものが使われていますが, これは, int だと整数, doubleだと浮動小数点数であるのと同じように, std::string だと「文字列」を保存する変数を意味しています．例えば

std::string AA="なすびはeggplant";

これで, AAの値が「なすびはeggplant」になります．利用するには, 先頭の#include <string> が必要ですが, 非常によく使われるものなので，省略してもどこかでincludeされ, 利用できるケースも多いです. さらに文法的には，この部分が新顔です．これはメンバー変数初期化リストといって，メンバー変数(引数の変数)と書きます. これで，メンバー変数の初期値が引数で与えられたものに設定されます. この例では, 青字のdtがメンバー変数ですので, STEP()の中ではメンバー変数dtが利用されます. ところがコンストラクタには引数dtが与えられているので, コンストラクタの中では引数dtが利用されます．こういうふうに, 変数には有効範囲(スコープ)があります．有効範囲は{}で区別されます:

例えば

{
    ここで double DX=1.0; と定義すると, DXは{}の中でだけ有効で浮動小数点数である
   {
   ここで int DX; と定義すると, DXは整数になってしまう
   }
   でもここではDXは浮動小数点数に戻るし値は1.0である. 上の整数DXは，無かったことになる
}

コンストラクタでは広い範囲で有効なメンバー変数dtが, コンストラクタの{}内では引数dtで隠されている，と考えることができます（例）．

数値計算的には, STEPのこの辺では, 前時間ステップの分布をU_OLDにコピーしてから, 新時刻のU_VALを計算していることに注意してくださいね．

プログラムを完成させると，こんな感じになります. 厳密解を点線で表示しています:

$c>0$であるため, 正方向に波が進んでいく様子がわかりますね！・・・でも，なんや，波形がだれまくっている気も, 多少しますね．

皆さんの仕事は，これをもう少しマシに解くことです．

プログラムが動くようになったら，以下の検討を行ってみてください：

$K$の値をCourant-Friedrichs-Lewy (CFL)条件という. この解法は, CFL < 1 の場合にしか利用できない.
- CFL > 1 では何が起こるか，確かめてみよ.
- CFLを色々に変えて, 厳密解と比較してみよ.
- $\Delta t, \Delta x$を小さくすると厳密解に近づくのか?
差分式$u_N=(1-K)u_+ +K u_-$は内挿・外挿の公式のような気がする. どのような内挿・外挿であるか, 考えてみよ.
初期値が滑らかな関数である場合, 初期値が不連続な関数である場合について, 計算を行ってみよ.
$\Delta x, \Delta t$をさまざまに変え, 誤差がどの程度であるか調べてみよ.
$c>0$の場合, 下図の差分では, 全く解けないことを確かめよ. $c<0$の場合には, どうか?
研究においては，上の挿絵のようなものを描く必要がある．描けるようになっておけ.
- PowerPointを使う人も多いが, レイヤーもないし, 精密に描けないようである．好みのものを見つけておく必要がある．LibreOffice/Draw, Visio, Dia, Inkscape, CorelDraw Blenderなどなど．
- LibreOffice/DrawとInkscapeは, graph で作成したPDFファイルを編集することができるので便利

陰的差分法

陽的差分法では, CFL条件によって, $\Delta t$を大きく選ぶことができません. 多くの問題では, 長時間経過して定常解すなわち\[\frac{\partial}{\partial t}=0\]となる状態が速く知りたいので, 困ります. そこで,

空間微分を, 解くべき時刻で評価する

ことを考えてみましょう. この場合, 解くべき時刻の関数値$u_N$に対する連立方程式が得られます(下図). 差分方程式の解が陽的に与えられず, 連立方程式の解として与えられるので, 陰的差分法と呼びます.

しかしこの場合，連立方程式といっても, 壁面から順番に解けば, 解が定まっていく特例になっているので, 簡単に試すことができるでしょう.

実際に$c>0$の場合についてプログラムを作成し, $\Delta t$が大きくても時間発展を求められることを確かめてください.
$c<0$の場合は, 差分スキームはどのようになるでしょうか?
差分スキームをよくみると, 内挿になっています. どの点とどの点を内挿したことになっていますか? 特性方向を用いて述べてください.

高精度差分

ここまでは, 2点のデータ$u_{i}$, $u_{i-1}$ を用いて最低精度の差分式を用いていました. もっと点数をたくさん選べば, 精度の高い差分方程式が得られます.

例えば$u_{i}$, $u_{i-1}$, $u_{i-2}$の3点を用いる場合, \[\frac{\partial u}{\partial x}=\alpha u_{i} + \beta u_{i-1} +\gamma u_{i-2}\]となるように$\alpha, \beta, \gamma$を選べば良い. 具体的には, $x_i$周りにTaylor展開して\[u_{i}=u_{i}\], \[u_{i-1}=u_{i}-\frac{\partial u}{\partial x}(x_i)\Delta x+\frac{1}{2}\frac{\partial^2 u}{\partial x^2}(x_i)\Delta x^2+\cdots\] \[u_{i-2}=u_{i}-\frac{\partial u}{\partial x}(x_i)2\Delta x+\frac{1}{2}\frac{\partial^2 u}{\partial x^2}(x_i)4\Delta x^2+\cdots\]となるので, $\alpha, \beta, \gamma$を選べば \[\frac{\partial u}{\partial x}=\alpha u_{i} + \beta u_{i-1} +\gamma u_{i-2}+O(\Delta x^3)\]となるように選ぶことができます．

高精度の差分を用いて, 方程式を解いてみなさい.
CFL条件は, どのように変わるか? $\Delta t $や$\Delta x$を変えて試してみてください.

時間方向に高精度の差分を用いると, 前時間ステップの分布をU_OLDにコピーしてから, 新時刻のU_VALを計算するだけでは不十分です. 前前時間ステップの分布をU_SUPEROLDにコピーしてから, 前時間ステップの分布をU_OLDにコピーしてから, 新時刻のU_VALを計算する必要があります．

データが大きくなると，コピーにも時間がかかるようになります．これを避けるために工夫しても良いでしょう．

オペレータの定義

面倒くさいので, U[nt][i] で時刻nt, 空間i のデータが出てくれれば良いのです．そこで, そのような時空間配列クラスを作成しましょう. 例えば, 時間データが3層(つまり, nt, nt-1, nt-2 にアクセスしたい)という場合

#include <array>
#include <vector>
...略...
class txData
{
private:
  static const int step=3;
  std::array<std::vector<double>,step> data;
public:
  std::vector<double>& operator[](int nt)
  {
    return data[nt%step];
  }
  txData(int N)
  {
     for(int nt=0;nt<step;nt++) data[nt].resize(N);
  }
};

ここの std::array<クラス名, 整数> というのは, std::vector<クラス名> と同じなのですが, サイズが指定した整数に固定されているものです. この例では,　要素が (doubleの配列, doubleの配列, doubleの配列)という3成分の配列になっています.　配列要素のdoubleの配列サイズがバラバラだと面倒なので, コンストラクタで全部N成分である，と統一しています．

ここの static ってのは, txDataのオブジェクトで共通の値である，という制約をかけます．const は, 予想できるように, 定数である, 変更不能である，と制約をかけるものです.

ただ単に制約をかけるだけではありません．計算が爆速になります. あなたのプログラムが遅いのは, constを入れないからだ，と断言できるくらいには速くなります.

この operator[] という関数は, オブジェクト変数名[n] という書き方をされた時に, どういう値が得られるか，を定義するものです．

具体的には, プログラム中で

txData U_VAL[nx];

と定義されたとして，U_VAL[nt] で得られる値を定めます. 上の例では, nt を3で割った余りが0の時には data[0] つまり最初の doubleの配列が得られます．ntが1増えると, 3で割った余りが1になるので data[1] つまりdoubleの配列が, さらにntが1増えると余りが2なので data[2] つまりdoubleの配列が得られます. ntがさらに1増えると, 余りが0に戻るのでdoubleの配列が得られます．このように，ntが増えるたびに巡回して参照するデータになります．こういうデータをリングメモリー, リング配列，などと呼びます．U_VAL[nt] と書けるのでプログラムは読みやすくなりますが, 実際に確保されているメモリーサイズは, ntが3個分しかないのでメモリーが節約できます．ntが増えたからといって「前前時間ステップの分布をU_SUPEROLDにコピー」する必要がなくなりますので，計算が高速化します.

U_VAL[nt] で得られるのは, doubleの配列ですので, その i 番目の要素は U_VAL[nt][i] で参照できます．結局

$u_i^{nt}$は, プログラム中では U_VAL[nt][i] と書けば良い.

となります．

オペレータは，どう定義しても良いのですが，多少の制約はあります．例えば[]は元々添字$a_i$オペレータですので, 引数が1つしかダメだとか, 制限があります．

波動方程式

2階の波動方程式\[\frac{\partial u}{\partial t}+c\frac{\partial v}{\partial x}=0,\quad\frac{\partial v}{\partial t}+c\frac{\partial u}{\partial x}=0\]ようするに\[\frac{\partial^2 u}{\partial t^2}-c^2\frac{\partial^2 u}{\partial x^2}=0\]を解いてみるのも面白いでしょう.

初期条件は,\[u(t=0,x)=\phi(x),\quad \frac{\partial u}{\partial t}(t=0,x)=\psi(x)\]でしたよね．遊ぶ子は育つってわけで，遊んでみましょう．

例1 滑らか初期値

\[u(t=0,x)=\exp\left[-5x^2\right],\quad \frac{\partial u}{\partial t}(t=0,x)=0\]

数学で勉強した通り！綺麗に左右に波が伝わっていますね！楽勝です．

例2 カクカク初期値

\[u(t=0,x)=\cases{1 &($|x|<1$)\\0&($|x|>1$)},\quad \frac{\partial u}{\partial t}(t=0,x)=0\]

ぎゃあ．不連続な初期値は厳しいわ〜. これでも$\Delta t$や$\Delta x$は，例1と同じなんですけどね．

いやあうまくいかないんでプログラム見せて〜って向きは，こちらをどうぞ．

偏微分方程式 I

陽的差分法

陰的差分法

高精度差分

オペレータの定義

波動方程式