みけCATのにっき(仮)
つれづれなるまゝに、日くらし、PCにむかひて、心に移りゆくよしなし事を、そこはかとなく書きつくれば、あやしうこそものぐるほしけれ。
(本当か!?)
出典

波動方程式の差分近似

アバター
みけCAT
記事: 6734
登録日時: 13年前
住所: 千葉県
連絡を取る:

波動方程式の差分近似

投稿記事 by みけCAT » 9年前

講義で波動方程式の差分近似による解法を紹介された時、ウェーブマシンのようなモノを描画したい。を思い出しました。
俺「波動ってことはウェーブじゃね?ということはもしかしたら使えるんじゃね?」

まずは、波動方程式を差分近似で解いて漸化式を出します。
hadouhouteisiki.png
波動方程式
hadouhouteisiki.png (24.19 KiB) 閲覧数: 811 回
そして、実装してみました。

CODE:

#include 
#include 

class Key {
	private:
		static const int KEY_NUM = 256;
		int keyStatus[KEY_NUM];
	public:
		Key() {
			for (int i = 0; i < KEY_NUM; i++) keyStatus[i] = 0;
		}
		bool updateKey() {
			char keyData[KEY_NUM];
			if (GetHitKeyStateAll(keyData) != 0) return false;
			for (int i = 0; i < KEY_NUM; i++) {
				keyStatus[i]++;
				if (!keyData[i]) keyStatus[i] = 0;
			}
			return true;
		}
		int getKey(int id) const {
			if (id < 0 || KEY_NUM <= id) return 0;
			return keyStatus[id];
		}
		int operator[](int id) const {
			return getKey(id);
		}
};

static bool getFloatValue(float *out, int x, int y) {
	char input[128] = "abc";
	int ret = KeyInputSingleCharString(x, y, sizeof(input) - 1, input, TRUE);
	if (ret == 1) {
		float buf;
		char *end;
		buf = strtof(input, &end);
		if (*end != '\0') return false;
		if (out != NULL) *out = buf;
		return true;
	} else {
		return false;
	}
}

int WINAPI WinMain(HINSTANCE,HINSTANCE,LPSTR,int) {
	if (ChangeWindowMode(TRUE) != DX_CHANGESCREEN_OK) return 1;
	SetMainWindowText("波動方程式の差分近似");
	if (DxLib_Init() != 0) return 1;
	SetDrawScreen(DX_SCREEN_BACK);
	Key key;

	const int WIDTH = 640;
	const int HEIGHT = 480;

	float h_x = 1.0f;
	float h_t = 1.0f;
	float s = 1.0f;
	bool isFixedEnd = true;

	float wave[3][WIDTH] = {{0.0f}};

	while (ProcessMessage() == 0 && ClearDrawScreen() == 0 && key.updateKey()) {
		// 波をローテートする(logrotate的な意味で)
		for (int i = 0; i < WIDTH; i++) {
			wave[2][i] = wave[1][i];
			wave[1][i] = wave[0][i];
		}
		// 次の波の位置を計算する
		double p = h_t * h_t / (s * s * h_x * h_x);
		for (int i = 1; i < WIDTH - 1; i++) {
			wave[0][i] = p * wave[1][i - 1] + 2.0f * (1.0f - p) * wave[1][i]
				+ p * wave[1][i + 1] - wave[2][i];
		}
		// 波の端を操作する
		int mx, my;
		GetMousePoint(&mx, &my);
		wave[0][0] = my - HEIGHT / 2;
		if (isFixedEnd) {
			// 固定端なので位置は0
			wave[0][WIDTH - 1] = 0.0f;
		} else {
			// 自由端 : よくわからないので、とりあえず右の値として左の値を用いる
			int idx = WIDTH - 1;
			wave[0][idx] = 2.0 * p * wave[1][idx - 1] + 2.0f * (1.0f - p) * wave[1][idx] - wave[2][idx];
		}

		// パラメータの設定を行う
		if (key[KEY_INPUT_F1] == 1) {
			float val = 0.0f;
			if (getFloatValue(&val, 10, 400)) h_x = val;
		} else if (key[KEY_INPUT_F2] == 1) {
			float val = 0.0f;
			if (getFloatValue(&val, 180, 400)) h_t = val;
		} else if (key[KEY_INPUT_F3] == 1) {
			float val = 0.0f;
			if (getFloatValue(&val, 350, 400)) s = val;
		} else if (key[KEY_INPUT_F4] == 1) {
			isFixedEnd = !isFixedEnd;
		} else if (key[KEY_INPUT_ESCAPE] == 1) {
			for (int i = 0; i < WIDTH; i++) {
				wave[0][i] = wave[1][i] = wave[2][i] = 0.0f;
			}
			SetMousePoint(mx, HEIGHT / 2);
		}

		// パラメータを描画する
		int clWhite = GetColor(255, 255, 255);
		DrawFormatString(10, 10, clWhite, "h_x = %f", h_x);
		DrawFormatString(180, 10, clWhite, "h_t = %f", h_t);
		DrawFormatString(350, 10, clWhite, "s = %f", s);
		DrawString(530, 10, isFixedEnd ? "固定端" : "自由端", clWhite);
		// 操作方法を描画する
		DrawString(10, 440, "F1 : h_x設定", clWhite);
		DrawString(180, 440, "F2 : h_t設定", clWhite);
		DrawString(350, 440, "F3 : s設定", clWhite);
		DrawString(520, 440, "F4 : 端設定", clWhite);
		// 波を描画する
		for (int i = 1; i < WIDTH; i++) {
			DrawLine(i - 1, (int)(wave[0][i - 1] + HEIGHT / 2.0),
				i, (int)(wave[0][i] + HEIGHT / 2.0), clWhite);
		}

		ScreenFlip();
	}

	DxLib_End();
	return 0;
}
結果は、h_x = h_t = s = 1.0の時はいい感じに波が伝わり、定常波のような現象も観測できましたが、
パラメータを変えると、波形が上下に飛んでいってしまうなどの微妙な感じになりました。
dousarei.png
動作例
dousarei.png (22.25 KiB) 閲覧数: 824 回
添付ファイル

[拡張子 zip は無効化されているため、表示できません]


アバター
GRAM
記事: 164
登録日時: 13年前

Re: 波動方程式の差分近似

投稿記事 by GRAM » 9年前

多分伝播速度が粒子間距離を上回るからじゃないかな(というようなことが質問ページの僕の主張)
CFDとかにおけるクーラン数条件に引っかかってる。多分。まじめにかんがえたわけじゃないけど。
要は陽に微分方程式を解くならば、現象の伝播速度は粒子間距離より小さくないといけない。
この例だと粒子間距離は1ドットなので、1フレームの伝播速度は1ドット未満である必要がある。
陰に説くならこの条件はなくなるけど普通にめんどくさい。
伝播速度を上げたいならば1秒間の計算回数を増やせばいいはず。
あとは計算誤差の問題だと思う。