clear; // *** 共通部分 *** // 入力パラメータ d = 0.05; f = 7.5; b = 1.0; t = 2 * %pi; // 時間の最大値 n = t / 0.01 + 1; // 時間の刻み数 T = linspace(0,t,n); // 初期値 x0 = 0.0; dx0 = 0.0; X0 = [x0; dx0]; // *** 常微分方程式ソルバによる数値解 *** // 微分方程式の定義 function dx = duffin(t,x) dx(1) = x(2); dx(2) = - d * x(2) - b * x(1) ^ 3 + f * cos(t); endfunction for k = 1:5000; if k - fix(k / 100) * 100 == 0 then k // 100ごとに数値をカウント表示 end X = ode(X0,0,T,duffin); // 常微分方程式ソルバ X0 = [X(1,n); X(2,n)]; XV(1,k) = X(1,n); XV(2,k) = X(2,n); end plot(XV(1,:),XV(2,:),'.r','markersize',1); xlabel("x"); ylabel("v");