clear; // *** 解析解とソルバ解の共通部分 *** // 初期振幅の入力 g = 9.8; // 重力加速度 l = g / (2 * %pi) ^ 2; // 糸の長さ // 初期条件 th0 = input("th0 = "); // 初角度 q0 = 0; // 初角速度 X0 = [th0; q0]; // *** 常微分方程式ソルバによる解 *** // 解くべき微分方程式の定義 function dx = pend(t,x) dx(1) = x(2); dx(2) = - g / l * sin(x(1)); endfunction T = linspace(0,2,200); // 時間ベクトル TH = ode(X0, 0, T, pend); // 常微分方程式ソルバ plot(T,TH(1,:),'-r'); // プロット // *** 解析解 *** THanaly = linspace(- th0, th0, 20); Tanaly = - sqrt(l / (2 * g)) * integrate('1 ./ sqrt(cos(th) - cos(th0))','th',th0,THanaly); plot(Tanaly,THanaly,'og'); // *** グラフの体裁 *** legend("o.d.e","integration",1); xlabel("t[s]"); ylabel("$\theta \mathrm{[rad]}$");