clear; // ************************************ // // 共通部分 // // ************************************ // *** 入力パラメータ *** a = 1; k1 = input("k1 = "); k2 = input("k2 = "); // 時間ベクトル T = linspace(0,10); // ************************************ // // 数値解 // // ************************************ // 解くべき微分方程式 function dx = unreversible(t,x) dx = k1 * (a - x) ^ 2 - k2 * x; endfunction // 初期値 x0 = 0; // 微分方程式の数値解 X = ode(x0,0,T,unreversible); // ************************************ // // 解析解 // // ************************************ b1 = (2 * a + k2 / k1 + sqrt(4 * a * k2 / k1 + (k2 / k1) ^ 2)) / 2; b2 = (2 * a + k2 / k1 - sqrt(4 * a * k2 / k1 + (k2 / k1) ^ 2)) / 2; d = sqrt(4 * a * k2 / k1 + (k2 / k1) ^ 2); Xa = b2 * (exp(k1 * d * T) - 1) ./ (exp(k1 * d * T) - b2 / b1); // ************************************ // // プロット // // ************************************ plot(T,X,'or'); plot(T,Xa,'-g'); xtitle(strcat(["Chemical Reaction: ","a = ",string(a),", k1 = ",string(k1),", k2 = ",string(k2)])); xlabel("$t [s]$"); ylabel("$N_x$"); legend(["Numerical","Analytical"],4);