はじめに
下記の記事 1（以下、元記事）にインスパイアされ、本記事を書きました。
大学入試問題をpythonで解く #Python - Qiita
ただ、元記事では解答にまで至っていなかったので、本記事では答えを求めてみます。
また、勉強を兼ねて Python から GNU Octave 2（以下、Octave）にコードを書き換えています。Octave は MATLAB 3 とほぼ互換性があるらしいです。
使用している Octave のバージョンは 5.2.0 です。
symbolic パッケージ 4 と geometry パッケージ 5 を使用しています。
SVG 画像を PNG 画像に変換するため、Inkscape 6 を利用しています。
Octave の紹介
Octave の「よいところ」と「よくないところ」を簡単に紹介します。※個人の感想です
詳細を知りたい方は、下記の公式ドキュメントを参照してください。
https://docs.octave.org/latest/index.html
よいところ
- 数値計算に適している（MATLAB と互換性がほぼある）
- ベクトル、行列を簡単に記述できる
- 線型代数、微分方程式、数値積分などの分野でよく利用する関数が用意されている
- プロットの GUI 表示、画像ファイル出力も簡単
- ライセンスが GPL （プロプライエタリでない）
上記「よいところ」の一部をコードで示します。
% 線型代数
b = [4; 9; 2]; % 列ベクトル. `'` が共軛転置演算子なので `[4 9 2]'` でも OK
A = [ 3 4 5; % 行列. 直感的に記述できる
1 3 1;
3 5 9 ];
x = A \ b % A を係数行列とする線型方程式系 Ax = b を解く. `A \ b` は `inverse(A) * b` と同値
%%% 出力
%% x =
%%
%% -1.5000
%% 4.0000
%% -1.5000
% 数値積分
f = @(x) 2 * sin(x); % f(x) = 2sin(x)
quad(f, 0, pi) % F(π) - F(0)
%%% 出力
%% ans = 4
% プロット
x = -10 : 0.1 : 10; % [-10 -9.9 -9.8 ... 9.9 10]. -10 から 10 までステップが 0.1 の行ベクトル
y = sin (x); % [sin(-10) sin(-9.9) sin(-9.8) ... sin(9.9) sin(10)]
plot (x, y);
title ("Simple 2-D Plot");
xlabel ("x");
ylabel ("sin (x)");
図 2.1-1 Octave の REPL にて、グラフをプロットしている様子（OS は Ubuntu）
よくないところ
- ドキュメントがわかりにくい
- 公式ドキュメントは英語
- 関数の使用例が少ない
- 日本語で使い方を解説してくれているサイトもあるが、情報が古い
- 実装されていない MATLAB の関数がある
-
polyshape()など
-
上記「よくないところ」の対応策です。
- ドキュメントがわかりにくい
- → MATLAB の日本語ドキュメントを併用する
- → Wiki とフォーラムを漁る
- → 公式ドキュメントとソースコードをがんばって読む
- 実装されていない MATLAB の関数がある
- → 代替可能な Octave のパッケージを探す
- → 自力で実装する。有料サポートを受ける
入試問題
\begin{align}
|x| + |y| - 2 & \leq 0 \tag{1} \\
x^2 + 2x - 2y - 1 & \geq 0 \tag{2} \\
\end{align}
\begin{array}{c}
(1), (2) を同時に満たす点(x, y)の全体からなる領域を D とする。 \\
D を図示し、その面積を求めよ。 \\
\end{array}
元記事では $(1)$ が $|x| + |y| - 2 \geq 0$ となっていますが、有限の紙面では $D$ の作図が出来ず、その面積も無限大になってしまうので変更を加えています。
解答
準備
式変形
問題を解きやすくするため、式変形しておきます。
まずは、$(1)$ から
\begin{align}
|x| + |y| - 2 & \leq 0 \\
|y| & \leq -|x| + 2\\
\end{align}
ここで、$y$ について場合分けします。
\begin{cases}
y \leq -|x| + 2 & (y \geq 0) &&& (1a) \\
y \gt |x| - 2 & (y \lt 0) &&& (1b) \\
\end{cases}
さらに、$(1a)$ を $x$ について場合分けして、
\begin{cases}
y \leq -x + 2 & (x \geq 0 \land y \geq 0) &&& (1.1) \\
y \leq x + 2 & (x \lt 0 \land y \geq 0) &&& (1.2) \\
\end{cases}
$(1b)$ も $x$ について場合分けします。
\begin{cases}
y \gt x - 2 & (x \geq 0 \land y \lt 0) &&& (1.3) \\
y \gt -x - 2 & (x \lt 0 \land y \lt 0) &&& (1.4) \\
\end{cases}
つづいて、$(2)$ を変形します。
\begin{align}
x^2 + 2x - 2y - 1 & \geq 0 \\
x^2 + 2x - 1 & \geq 2y \\
y & \leq \frac{1}{2}x^2 + x - \frac{1}{2} \tag{2.1} \\
\end{align}
スクリプト
入試問題中で定義された数式や先程変形した数式、他のスクリプトから呼ばれるヘルパ関数などを定義しておきます。
各スクリプトの詳細は、下記ページを参照してください。
https://github.com/reodon/akbrobot/tree/main/advent-calendar/2023_theory/04th/octave
Octave を実行するワーキングディレクトリのファイル構造を下記に示します。
.
├── out.d/
├── private/
│ ├── arg2num.m
│ ├── conf_axes.m
│ ├── f_1.m
│ ├── f_11.m
│ ├── f_12.m
│ ├── f_13.m
│ ├── f_14.m
│ ├── f_2.m
│ ├── f_21.m
│ ├── fill_poly.m
│ ├── fprint_hr.m
│ ├── fprint_lf.m
│ ├── plot_f1_f2.m
│ └── save_plot.m
├── README.md
├── area_integral_numeric.m
├── area_montecarlo.m
├── area_polyarea.m
├── plot_d.m
├── plot_integral.m
├── plot_manual.m
├── plot_montecarlo.m
├── plot_ready.m
└── put_integral_symbolic.m
プロット
$(2.1)$ のグラフの形状を把握するため、その右辺 $\frac{1}{2}x^2 + x - \frac{1}{2}$ に対して、二次方程式の解の公式を適用します。
\frac{-b \pm \sqrt{b^2 - 4ac}}{2a} = -1 \pm \sqrt{2}
さらに、最小値を求めるために微分します。
\begin{align}
\biggl( \frac{1}{2}x^2 + x - \frac{1}{2} \biggr)' & = x + 1 \\
x & = -1 \\
\end{align}
$(2.1)$ について、$x = -1 \pm \sqrt{2}$ のとき $y = 0$ となり、$x = -1$ のとき $y$ が最小となることがわかりました。
念の為、$(2.1)$ の根と最小値を Octave で計算しておきます。
% (2.1) の根と最小値の検算
pkg load symbolic;
syms f_21(x) x;
f_21(x) = (sym(1)/2)*x.^2 + x - (sym(1)/2);
roots(sym2poly(f_21)) % 根
%%% 出力
%% ans =
%%
%% -2.41421
%% 0.41421
[-1 - sqrt(2), -1 + sqrt(2)]' % 手計算の根
%%% 出力
%% ans =
%%
%% -2.41421
%% 0.41421
fminbnd(function_handle(f_21(x)), -3, 3) % 最小値. 範囲は -3 <= x <= 3
%%% 出力
%% ans = -1
手計算は間違っていないようです。
大体の形状が把握できたところで、下記スクリプトを実行してグラフ画像を生成します。
function plot_ready()
clf; hold on;
plot_f1_f2(true); conf_axes;
pin_axis_x(-1 + sqrt(2), "-1 + √2");
pin_axis_x(-1 - sqrt(2), "-1 - √2");
save_plot(dbstack.name);
endfunction
function pin_axis_x(x, t)
line(repelem(x, 2), [-.1 .1], "color", [128 128 128]./256, "linewidth", 1.5);
text(x, -.1, t, "fontsize", 14, "verticalalignment", "top", "horizontalalignment", "center");
endfunction
$ octave plot_ready.m
上記のコマンドを実行して、生成された画像がこちらです。
図 4.1.3-1 plot_ready.png
Octave によるグラフのプロットは、問題なさそうです。
作図
$D$ を図示します。
fill() という関数を利用します。
塗りつぶされた多角形領域をプロットします。
-
fill()
function fill_poly(poly, fillopt = {}, centroid_text = [])
pkg load geometry;
ccw_poly = orientPolygon(poly);
fill(ccw_poly(:, 1), ccw_poly(:, 2), fillopt{:}); % <--- fill()
if (ischar(centroid_text))
c = shapeCentroid(polygon2shape(ccw_poly));
text(c(1), c(2), centroid_text, "fontsize", 14, "horizontalalignment", "center");
endif
endfunction
function plot_d()
clf; hold on;
plot_f1_f2; conf_axes;
x{1} = -1 : .01 : 1;
x{2} = 1 : .01 : 2;
x{3} = 2 : -.01 : 0;
x{4} = 0 : -.01 : -1;
y{1} = f_21(x{1});
y{2} = f_11(x{2});
y{3} = f_13(x{3});
y{4} = f_14(x{4});
xs = cell2mat(x);
ys = cell2mat(y);
fill_poly([xs' ys'], {"y"}, "D");
save_plot(dbstack.name);
endfunction
$ octave plot_d.m
図 4.2-1 plot_d.png
作図は大丈夫そうです。
面積
手計算
Octave で計算をする前に、手計算で答えを出しておきます。
下記に定義する $D_{man}$ を定積分して、残りの正方形部分の面積（$D \setminus D_{man}$）を足す方針にします。
\begin{align}
x & \leq 1 \tag{3} \\
y & \geq -1 \tag{4} \\
\end{align}
\begin{array}{c}
(2), (3), (4) を同時に満たす点(x, y)の全体からなる領域を D_{man} とする。
\end{array}
$ octave plot_manual.m
図 4.3.1-1 plot_manual.png
\begin{align}
D_{man} & = \int^1_{-1} \Big( \frac{1}{2}x^2 + x - \frac{1}{2} \Big) - (-1) dx \\
& = \frac{1}{2}\left[ \frac{1}{3}x^3 + x^2 + x\right]^1_{-1} \\
& = \frac{1}{2}\left[\Big(\frac{1}{3} + 1 + 1\Big) - \Big(-\frac{1}{3} + 1 - 1\Big) \right] \\
& = \frac{4}{3} = 1.3333\ldots \\
\end{align}
$D$ の面積は
\begin{align}
D & = D_{man} + (D \setminus D_{man}) \\
& = \Big(\frac{4}{3}\Big) + (\sqrt{2})^2 \\
& = \frac{10}{3} = 3.3333\ldots
\end{align}
計算間違いが怖いので、Octave で $D_{man}$ の面積を計算しておきます。
% D_{man} の検算
pkg load symbolic;
syms f(x) x;
f(x) = ((sym(1)/2)*x.^2 + x - sym(1)/2) - (-1);
quad(function_handle(f), -1, 1)
%%% 出力
%% ans = 1.3333
4/3
%%% 出力
%% ans = 1.3333
計算間違いは無かったようです。
モンテカルロ法
モンテカルロ法で面積を求めます。
function area_montecarlo(n = 10^4)
n = arg2num(n);
limit = 3; % x, y の絶対値の最大値
x = (rand(n, 1)*2 - 1) * limit;
y = (rand(n, 1)*2 - 1) * limit;
condition = @(x, y) f_1(x, y) <= 0 && f_2(x, y) >= 0;
cnt = length(find(arrayfun(condition, x, y)));
p = cnt / n;
D = p * (limit*2)^2;
fprintf(stderr, "%f\n", D); % 標準出力に不要な文字列が表示されるため、面積は標準エラー出力に表示
endfunction
試行回数（n）が 1,000, 10,000, 100,000 のとき、それぞれどのような値になるか見ていきます。
n = 1,000 (1k)
$ n=$((10*3)); time (for i in {1..5}; do octave area_montecarlo.m $n > /dev/null; done)
2.772000
2.916000
2.808000
3.024000
3.996000
real 0m2.393s
user 0m2.718s
sys 0m1.192s
$ n=$((10**3)); octave plot_montecarlo.m $n
図 4.3.2-1 n = 1k のときのモンテカルロ法のプロット
n = 10,000 (10k)
$ n=$((10**4)); time (for i in {1..5}; do octave area_montecarlo.m $n > /dev/null; done)
3.124800
3.236400
3.502800
3.236400
3.416400
real 0m4.592s
user 0m4.809s
sys 0m1.297s
$ n=$((10**4)); octave plot_montecarlo.m $n
図 4.3.2-2 n = 10k のときのモンテカルロ法のプロット
n = 100,000 (100k)
$ n=$((10**5)); time (for i in {1..5}; do octave area_montecarlo.m $n > /dev/null; done)
3.290760
3.286080
3.355920
3.343320
3.363840
real 0m31.827s
user 0m30.705s
sys 0m1.736s
$ n=$((10**5)); octave plot_montecarlo.m $n
図 4.3.2-3 n = 100k のときのモンテカルロ法のプロット
考察
試行回数が多くなるほど、誤差は小さく、実行時間は大きくなっているようです。
積分
$D$ を下記の3つの領域（$D_{1}$, $D_{2}$, $D_{3}$）に分割して、それぞれを定積分して足し合わせる方針にします。
$ octave plot_integral.m
図 4.3.3-1 定積分で面積を求めるにあたり、分割した3つの領域
数値積分
quad() を利用します。
数値積分の結果を得ます。
-
quad()
function area_integral_numeric()
D(1) = quad(@(x) f_21(x) - f_14(x), -1, 0);
D(2) = quad(@(x) f_21(x) - f_13(x), 0, 1);
D(3) = quad(@(x) f_11(x) - f_13(x), 1, 2);
fprintf(stderr, "%f\n", sum(D));
endfunction
$ octave area_integral_numeric.m > /dev/null
3.333333
記号積分
symbolic パッケージを利用します。
数式の記号的操作を行うためのパッケージです。
https://ja.wikipedia.org/wiki/数式処理システム
-
symbolicパッケージ
function put_integral_symbolic()
pkg load symbolic;
syms x;
s.g(1) = f_21(x) - f_14(x);
s.g(2) = f_21(x) - f_13(x);
s.g(3) = f_11(x) - f_13(x);
s.D(1) = int(s.g(1), -1, 0);
s.D(2) = int(s.g(2), 0, 1);
s.D(3) = int(s.g(3), 1, 2);
fd = fopen(strcat("./out.d/", dbstack.name, ".out"), "w");
fprint_struct(fd, s); % 出力をフォーマットするため、ファイルに書き出します
fclose(fd);
endfunction
function fprint_struct(fd, s)
fprint_lf(fd);
fprintf(fd, "D ="); fdisp(fd, pretty(sum(s.D)));
for i = 1:length(s.g)
fprint_hr(fd);
fprintf(fd, strcat("D", num2str(i), " =")); fdisp(fd, pretty(s.D(i)));
fprintf(fd, strcat("g", num2str(i), " =")); fprint_lf(fd); fprint_lf(fd);
fdisp(fd, pretty(s.g(i)));
endfor
endfunction
$ octave put_integral_symbolic.m >& /dev/null
$ cat out.d/put_integral_symbolic.out
D = 10/3
-----------------------
D1 = 2/3
g1 =
2
x 3
── + 2⋅x + ─
2 2
-----------------------
D2 = 5/3
g2 =
2
x 3
── + ─
2 2
-----------------------
D3 = 1
g3 =
4 - 2⋅x
その他
polyarea()
polyarea() を利用します。
多角形の面積を求める関数です。
-
polyarea()
function area_polyarea(step = .01)
step = arg2num(step);
x{1} = -1 : step : 1;
x{2} = 1 : step : 2;
x{3} = 2 : -step : 0;
x{4} = 0 : -step : -1;
y{1} = f_21(x{1});
y{2} = f_11(x{2});
y{3} = f_13(x{3});
y{4} = f_14(x{4});
D = polyarea(cell2mat(x), cell2mat(y));
fprintf(stderr, "%f\n", D);
endfunction
$ octave area_polyarea.m .1 > /dev/null
3.335000
$ octave area_polyarea.m .001 > /dev/null
3.333333
step を小さくするほど、誤差も小さくなるようです。
おわりに
Qiita で思い通りに数式（tex）が書けず、Octave と MATLAB の違いに苦しみましたが、色々と勉強になりました。
手計算を何度も間違えたあとで、symbolic パッケージの存在を知りましたが、もっと早く見つけていれば効率的に記事を書けたかもしれません。
時間があれば、Julia 7 で書き直そうと思います。