clear;

// *** 共通部分 ***
// 入力パラメータ
k = 1.0;
w = %pi / 2;
t = 80;
d = 0.2;
f = 1.0;
t = 80;			// 時間の最大値
n = t / 0.1 + 1;	// 時間の刻み数
T = linspace(0,t,n);
// 初期値
x0 = 1.0;
dx0 = 0.5;
X0 = [x0; dx0];

// *** 常微分方程式ソルバによる数値解 ***
// 微分方程式の定義
function dx = dump(t,x)
  dx(1) = x(2);
  dx(2) = - d * x(2) - k * x(1) + f * cos(w * t);
endfunction
// 常微分方程式ソルバ
X = ode(X0,0,T,dump);

// *** 解析解 ***
kdw = (k - w ^ 2) ^ 2 + (d * w) ^ 2;
ac = x0 - f * (k - w ^ 2) / kdw;
as = - (dx0 + d / 2 * ac - f * d * w ^ 2 / kdw) / sqrt(k - d ^ 2 / 4);
th = atan(as / ac);
a = ac / cos(th);
ph = - atan(d * w, k - w ^ 2);
AX = a * exp(- d / 2 * T) .* cos(sqrt(k - d ^ 2 / 4) * T + th) ..
     + f * cos(w * T + ph) / sqrt(kdw);

// *** グラフのプロット ***
// 時間-位置プロット
xsetech([0,0,1,0.45]);
plot(T,X(1,:),'or');	// 数値解
plot(T,AX,'--g')	// 解析解
// グラフの体裁
xlabel("t");
ylabel("x(t)");
legend("Numerical","Analytical");
// 位置-速度(位相図)プロット
xsetech([0,0.5,1,0.45]);
plot(X(1,:),X(2,:),'-r');
xlabel("x");
ylabel("v");