/*
* dynamics.c
*/
/**********************************************************************
**********************************************************************
*
* 2 次元の定数係数線形常微分方程式
* x'(t) = a x + b y
* y'(t) = c x + d y
* に初期条件
* x(0)=x0
* y(0)=y0
* を課した常微分方程式の初期値問題を解いて、相図を描く。
*
* このプログラムは次の4つの部分から出来ている。
* main()
* 主プログラム
* 行列の係数の入力、ウィンドウのオープン等の初期化をした後に、
* ユーザーにメニュー形式でコマンドを入力してもらう。
* 実際の計算のほとんどは他の副プログラムに任せている。
* draworbit(x0,y0,h,tlimit)
* (x0,y0) を初期値とする解の軌道を描く。
* 刻み幅を h、追跡時間を tlimit とする。
* 近似解の計算には Runge-Kutta 法を用いる。
* fx(x,y)
* 微分方程式の右辺の x 成分
* fy(x,y)
* 微分方程式の右辺の y 成分
**********************************************************************/
#include <stdio.h>
#include <math.h>
#include "fplot.h"
/* 係数行列 A の成分（common 文により function fx,fy と共有する）*/
double a, b, c, d;
/* ウィンドウに表示する範囲（common 文により draworbit と共有する）*/
double xleft, xright, ybottom, ytop;
void draworbit(double, double, double, double);
int main()
{
/* 初期値 */
double x0, y0;
/* 時間の刻み幅、追跡時間 */
double h, tlimit, newh, newT;
/* メニューに対して入力されるコマンドの番号 */
int cmd;
/* マウスのボタンの状態 */
int lbut, mbut, rbut;
/* ウィンドウに表示する範囲の設定 */
xleft = -1.0;
xright = 1.0;
ybottom = -1.0;
ytop = 1.0;
/* 時間刻み幅、追跡時間（とりあえず設定） */
h = 0.01;
tlimit = 10.0;
/* 行列の成分を入力 */
printf(" a,b,c,d=");
scanf("%lf %lf %lf %lf", &a, &b, &c, &d);
/* ウィンドウを開く*/
openpl();
fspace2(xleft, ybottom, xright, ytop);
/* x 軸、y 軸を描く */
fline(xleft, 0.0, xright, 0.0);
fline(0.0, ybottom, 0.0, ytop);
xsync();
/* メイン・ループの入口 */
while (1) {
/* メニューを表示して、何をするか、番号で選択してもらう */
printf("したいことを番号で選んで下さい。\n");
printf(" -1:メニュー終了, 0:初期値のキー入力, 1:初期値のマウス入力,");
printf(" 2:刻み幅 h,追跡時間 T 変更(現在 h=%7.4f, T=%7.4f)\n",
h, tlimit);
scanf("%d", &cmd);
/* 番号 cmd に応じて、指示された仕事をする */
if (cmd == 0) {
/* 初期値の入力*/
printf(" 初期値 x0,y0=");
scanf("%lf %lf", &x0, &y0);
draworbit(x0, y0, h, tlimit);
} else if (cmd == 1) {
while (1) {
printf("マウスの左ボタンで初期値を指定して下さい (右ボタンで中止)\n");
fmouse(&lbut, &mbut, &rbut, &x0, &y0);
if (lbut == 1) {
printf("(x0,y0)=(%g,%g)\n", x0, y0);
draworbit(x0, y0, h, tlimit);
} else if (mbut == 1) {
printf("(x0,y0)=(%g,%g)\n", x0, y0);
draworbit(x0, y0, -h, tlimit);
} else {
printf("マウスによる初期値の入力を打ち切ります。\n");
break;
}
}
} else if (cmd == 2) {
/* 時間刻み幅、追跡時間の変更 */
printf(" 時間刻み幅 h, 追跡時間 T= ");
scanf("%lf %lf", &newh, &newT);
if (newh != 0 && newT != 0) {
h = newh;
tlimit = newT;
printf("新しい時間刻み幅 h = %g, 新しい追跡時間 T = %g\n",
h, tlimit);
} else {
printf(" h=%g, T=%g は不適当です。\n", newh, newT);
}
} else if (cmd == -1) {
/* 終了 -- メイン・ループを抜ける */
break;
}
}
mkplot("dynamics.plot");
printf("fplotウィンドウを左ボタンでクリックして下さい\n");
closepl();
return 0;
}
/***********************************************************************/
/***********************************************************************/
/* 指示された初期値に対する解軌道を描く*/
void draworbit(double x0, double y0, double h, double tlimit)
{
double x, y, fx(), fy();
double k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y, t;
/* 時刻を 0 にセットする */
t = 0.0;
/* 初期値のセット */
x = x0;
y = y0;
/* 初期点を描く */
fpoint(x, y);
/* ループの入口 */
do {
/* Runge-Kutta 法による計算 */
/* k1 の計算 */
k1x = h * fx(x, y);
k1y = h * fy(x, y);
/* k2 の計算 */
k2x = h * fx(x + k1x / 2.0, y + k1y / 2.0);
k2y = h * fy(x + k1x / 2.0, y + k1y / 2.0);
/* k3 の計算 */
k3x = h * fx(x + k2x / 2.0, y + k2y / 2.0);
k3y = h * fy(x + k2x / 2.0, y + k2y / 2.0);
/* k4 の計算 */
k4x = h * fx(x + k3x, y + k3y);
k4y = h * fy(x + k3x, y + k3y);
/* (Xn+1,Yn+1) の計算 */
x += (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
y += (k1y + 2.0 * k2y + 2.0 * k3y + k4y) / 6.0;
/* 解軌道を延ばす */
fcont(x, y);
/* 時刻を 1 ステップ分進める */
t += h;
/* まだ範囲内かどうかチェック */
}
while (abs(t) <= fabs(tlimit));
/* サブルーチンを抜ける */
/* 確実に描画が終るのを待つ */
xsync();
}
/**********************************************************************/
/**********************************************************************/
/* 微分方程式の右辺のベクトル値関数 f の x 成分 */
double fx(double x, double y)
{
return a * x + b * y;
}
/* 微分方程式の右辺のベクトル値関数 f の y 成分 */
double fy(double x, double y)
{
return c * x + d * y;
}