波動方程式の数値解法

数値解法

\frac { 1 }{ { c }^{ 2 } } \frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } =\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ x }^{ 2 } } +\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }y^{ 2 } } 

今回は、二次元波動方程式を差分法にて数値積分を試みます。
二次元程度でしたら、一般的なコンピュータならリアルタイムにその時間発展を観察することができるので、とても楽しいです。

導出についてはこちら
波動方程式の導出（弦の運動） - Qiita

式変形

コンピュータに数値計算させるためには、ある程度時間(t)と空間(x)を細かく区切って進めるような式を準備してあげる必要があります。
まずは空間のほうですが、以下の考え方が使えます。これを差分法といいます。

\frac { { \partial  }u(x,y,t) }{ { \partial  }{ x } } \cong \frac { u(x,y,t)-u(x-⊿x,y,t) }{ ⊿x } \\ \frac { { \partial  }u(x,y,t) }{ { \partial  }{ x } } \cong \frac { u(x+⊿x,y,t)-u(x,y,t) }{ ⊿x } \\ \frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ x }^{ 2 } } \cong \frac { \frac { u(x+⊿x,y,t)-u(x,y,t) }{ ⊿x } -\frac { u(x,y,t)-u(x-⊿x,y,t) }{ ⊿x }  }{ ⊿x } =\frac { u(x-⊿x,y,t)-2u(x,y,t)+u(x+⊿x,y,t) }{ { ⊿x }^{ 2 } } 

一番下の式も、上の式と結局のところ同じことをやっているだけです。
これにより、式から位置に対する二階微分を消去します。

\frac { 1 }{ { c }^{ 2 } } \frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } =\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ x }^{ 2 } } +\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }y^{ 2 } } \\ \frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } ={ c }^{ 2 }(\frac { u(x-⊿x,y,t)-2u(x,y,t)+u(x+⊿x,y,t) }{ { ⊿x }^{ 2 } } +\frac { u(x,y-⊿y,t)-2u(x,y,t)+u(x,y+⊿y,t) }{ { ⊿y }^{ 2 } } )

ここで、⊿x=⊿yだということにして、右辺を整理します。

\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } ={ c }^{ 2 }(\frac { u(x-⊿x,y,t)-2u(x,y,t)+u(x+⊿x,y,t) }{ { ⊿x }^{ 2 } } +\frac { u(x,y-⊿x,t)-2u(x,y,t)+u(x,y+⊿x,t) }{ { ⊿x }^{ 2 } } )\\ \frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } ={ c }^{ 2 }\frac { u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) }{ { ⊿x }^{ 2 } } \\ \frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } =\frac { { c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} 

とだいぶ式が簡単になりました。
ここでさらに左辺を同じ方法で整理します。

\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } =\frac { { c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} \\ \frac { u(x,y,t-⊿t)-2u(x,y,t)+u(x,y,t+⊿t) }{ { ⊿t }^{ 2 } } =\frac { { c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} \\ u(x,y,t-⊿t)-2u(x,y,t)+u(x,y,t+⊿t)=\frac { { { ⊿t }^{ 2 }c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} \\ u(x,y,t+⊿t)=2u(x,y,t)-u(x,y,t-⊿t)+\frac { { { ⊿t }^{ 2 }c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} 

これで準備は整いました。

プログラミング

さて、上記の式展開によって、u(x,y,t)の二次元配列と、u(x,y,t-⊿t)の二次元配列があれば、
適当な初期値を定め、どんどん時間発展させることが可能です。
端っこの処理が気になりますが、今回は簡単に、端っこ一マスは処理しないで常に0であることにします。
この端っこの処理を考えるときの規則を境界条件と呼び、常に0の初期値に固定してしまうようなものをディリクレ境界条件と呼ぶそうです。一方初期状態を決める規則を初期条件と呼びます。（この辺あまり自信がありません）
数学的には積分定数の確定が目的です。

これを踏まえてC++コードに素直に落とし込むと、以下のようになります。
数式をしっかり整理することで、効率の良いコードに落とし込みやすくなります。

#include <glm/glm.hpp>
#include <glm/ext.hpp>

class WaveSolver {
public:
	WaveSolver() {
		_u_cur.resize(kWaveGrid * kWaveGrid);
		_u_new.resize(kWaveGrid * kWaveGrid);
		_u_pre.resize(kWaveGrid * kWaveGrid);
		this->reset();
	}
	void reset() {
		auto gauss = [](float x, float sigma) {
			return 1.0f / glm::sqrt(glm::two_pi<float>()) * sigma * glm::exp(-x * x / (2.0f * sigma * sigma));
		};

		int cx1 = kWaveGrid / 3;
		int cy1 = kWaveGrid / 3;
		int cx2 = kWaveGrid * 2 / 3;
		int cy2 = kWaveGrid * 2 / 3;
		for (int x = 1; x < kWaveGrid - 1; ++x) {
			for (int y = 1; y < kWaveGrid - 1; ++y) {
				int index = this->valueIndex(x, y);
				float value = 0.0f;

				{
					float norm = glm::distance(glm::vec2(x, y), glm::vec2(cx1, cy1));
					value += gauss(norm, 3.0f) * 20.0f;
				}
				{
					float norm = glm::distance(glm::vec2(x, y), glm::vec2(cx2, cy2));
					value += gauss(norm, 3.0f) * 20.0f;
				}
				
				_u_cur[index] = value;
			}
		}

		_u_pre = _u_new = _u_cur;
	}
	void step() {
		float deltaX = float(kWaveWidth) / float(kWaveGrid);
		float deltaT = 1.0f / 60.0f;
		float c = 2.0f;
		float mul = deltaT * deltaT * c * c / (deltaX * deltaX);
		for (int x = 1; x < kWaveGrid - 1; ++x) {
			for (int y = 1; y < kWaveGrid - 1; ++y) {
				int index = this->valueIndex(x, y);

				float uL = _u_cur[this->valueIndex(x - 1, y)];
				float uR = _u_cur[this->valueIndex(x + 1, y)];
				float uT = _u_cur[this->valueIndex(x, y - 1)];
				float uB = _u_cur[this->valueIndex(x, y + 1)];

				float u_pre = _u_pre[index];
				float u = _u_cur[index];
				_u_new[index] = u + u - u_pre + mul * (-4.0f * u + uL + uR + uT + uB);
			}
		}

		std::swap(_u_pre, _u_cur);
		std::swap(_u_cur, _u_new);
	}

	enum {
		kWaveGrid = 300,
		kWaveWidth = 100
	};

	int valueIndex(int x, int y) const {
		return y * kWaveGrid + x;
	}
	float value(int x, int y) const {
		return _u_cur[valueIndex(x, y)];
	}
	int widthN() const {
		return kWaveGrid;
	}
	int heightN() const {
		return kWaveGrid;
	}

	float width() const {
		return kWaveWidth;
	}
	float height() const {
		return kWaveWidth;
	}
private:
	std::vector<float> _u_pre;
	std::vector<float> _u_cur;
	std::vector<float> _u_new;
};

本質的なコードはstep関数の中の一部の行だけで、あとはさして重要ではありません。

これをビジュアライズすると、このようになります。
きれいですね！

減衰

よくありふれた現象として、減衰があります。減衰を考慮にいれると、外力などが加わらない限り、いずれ波の高さが減っていき最後には0に収束していきます。
今回は物理的な根拠はありませんが、速度に比例する減衰項（定数kとして）を入れてみます。粘性抵抗のような考え方です。
これをほとんど上記と変わりませんが、最後の項だけ少し変わってきます。

\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } ={ c }^{ 2 }(\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ x }^{ 2 } } +\frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }y^{ 2 } } )-k\frac { { \partial  }u(x,y,t) }{ { \partial  }t } \\ \frac { { \partial  }^{ 2 }u(x,y,t) }{ { \partial  }{ t }^{ 2 } } =\frac { { c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} -k\frac { { \partial  }u(x,y,t) }{ { \partial  }t } \\ \frac { u(x,y,t-⊿t)-2u(x,y,t)+u(x,y,t+⊿t) }{ { ⊿t }^{ 2 } } =\frac { { c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} -k\frac { u(x,y,t)-u(x,y,t-⊿t) }{ ⊿t } \\ u(x,y,t-⊿t)-2u(x,y,t)+u(x,y,t+⊿t)=\frac { { { ⊿t }^{ 2 }c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} -k⊿t\{ u(x,y,t)-u(x,y,t-⊿t)\} \\ u(x,y,t+⊿t)=2u(x,y,t)-u(x,y,t-⊿t)+\frac { { { ⊿t }^{ 2 }c }^{ 2 } }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} -k⊿t\{ u(x,y,t)-u(x,y,t-⊿t)\} 

プログラム側も該当箇所だけアップデートします。

int index = this->valueIndex(x, y);

float uL = _u_cur[this->valueIndex(x - 1, y)];
float uR = _u_cur[this->valueIndex(x + 1, y)];
float uT = _u_cur[this->valueIndex(x, y - 1)];
float uB = _u_cur[this->valueIndex(x, y + 1)];

float u_pre = _u_pre[index];
float u = _u_cur[index];

float k = 0.1f;
float damp = - k * deltaT * (u - u_pre);
_u_new[index] = u + u - u_pre + mul * (-4.0f * u + uL + uR + uT + uB) + damp;

こちらもビジュアライズすると確かにいいかんじに減衰していることが確認できました。

是非みなさんも色々数式やらコードやらをいじって楽しんでみてはいかがでしょうか