ベタ書きのコード
先ほどの説明の関数の文字を置き換えます。
- $T$→$u$
- $u$→$c$
１次元の移流方程式
\begin{align*}\frac{\partial u(x,t)}{\partial t}+c\frac{\partial u(x,t)}{\partial x}=0\tag{7}\end{align*}
※関数$u(x)$が速度$c$で伝播する方程式
C++で移流方程式をそのまま後退差分で実装すると以下のようになります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 | #include <iostream> #include <vector> using namespace std; using std::vector; int main() { int i, nt; int imax = 100, ntmax = 500; double dt = 0.001, xmax=2.0, c=1.0; double dx,cfl; vector<double> x(imax, 0.0); vector<double> U(imax, 0.0); dx = xmax / imax; cout << "dx = " << dx << endl; cfl = c * dt / dx; //initial state for (i = 0; i < x.size(); i++) { x[i] = dx * i; cout << "x[" << i << "] = " << x[i] << endl; if ((x[i] >= 0.5) && (x[i] <= 1.0)) { U[i] = 1.0; } else { U[i] = 0.0; } } // initial state for (i = 0; i < x.size(); i++) { cout << "U[" << i<< "] = " << U[i] << endl; } // Time Step for (nt = 0; nt < ntmax; nt++) { for (i = 1; i < x.size() - 1; i++) { U[i] = U[i] - cfl * (U[i] - U[i - 1]); } U[0] = U[1]; U[x.size() - 1] = U[x.size() - 2]; } // final state for (i = 0; i < x.size(); i++) { cout << "U[" << i << "] = " << U[i] << endl; } return 0; } |
解きたい方程式が簡単なのでコード量も少ないですね。
しかし、このようなべた書きコードを書いたら1次元から3次元に変更しようと思うと大がかりな改造する必要で1回限りの使い捨てコードになってしまいます。
メインの差分法は、
1 | U[i] = U[i] - cfl * (U[i] - U[i - 1]); |
ですが、直感的には、
1 | U = U - cfl * grad(U); |
のように書けると嬉しいですよね。
本記事を通してそのように書くことが可能になります。
割と直感的なコード。
移流方程式って難しい。 pic.twitter.com/QGcPUHTXMG— カマキリ🐲Python頑張る昆虫 (@t_kun_kamakiri) July 14, 2022
オブジェクト指向型プログラミング
では、ここからオブジェクト指向でコードを修正していきます。
段階を踏んでコードを書いていくことにします。
Follow @t_kun_kamakiri