clear; // *** 解析解とソルバ解の共通部分 *** g = 9.8; // 重力加速度 l = g / (2 * %pi) ^ 2; // 糸の長さ // 初期条件 th0 = input("th0 = "); // 初角度 q0 = 0; // 初角速度 X0 = [th0; q0]; T = linspace(0,2,100); // 時間ベクトル // *** 厳密解 *** // 解くべき微分方程式の定義 function dx = pend(t,x) dx(1) = x(2); dx(2) = - g / l * sin(x(1)); endfunction TH = ode(X0, 0, T, pend); // 常微分方程式ソルバ plot(T,TH(1,:),'-r'); // プロット // *** 近似解 *** // 解くべき微分方程式の定義 function dx = pendapp(t,x) dx(1) = x(2); dx(2) = - g / l * x(1); endfunction THapp = ode(X0, 0, T, pendapp); // 常微分方程式ソルバ plot(T,THapp(1,:),'og'); // プロット // *** グラフの体裁 *** legend("exact solution","approximate solution",1); xlabel("t[s]"); ylabel("$\theta \mathrm{[rad]}$");