LoginSignup
0
0

大学入試問題をoctaveで解く

Last updated at Posted at 2024-05-15

はじめに

下記の記事 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)");

screenshot_octave_repl_ubuntu.png
図 2.1-1 Octave の REPL にて、グラフをプロットしている様子(OS は Ubuntu)

よくないところ

  • ドキュメントがわかりにくい
    • 公式ドキュメントは英語
    • 関数の使用例が少ない
    • 日本語で使い方を解説してくれているサイトもあるが、情報が古い
  • 実装されていない MATLAB の関数がある

上記「よくないところ」の対応策です。

入試問題

\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

手計算は間違っていないようです。

大体の形状が把握できたところで、下記スクリプトを実行してグラフ画像を生成します。

plot_ready.m
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

上記のコマンドを実行して、生成された画像がこちらです。

plot_ready.png
図 4.1.3-1 plot_ready.png

Octave によるグラフのプロットは、問題なさそうです。

作図

$D$ を図示します。

fill() という関数を利用します。
塗りつぶされた多角形領域をプロットします。

private/fill_poly.m
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
plot_d.m
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

plot_d.png
図 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

plot_manual.png
図 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

計算間違いは無かったようです。

モンテカルロ法

モンテカルロ法で面積を求めます。

area_montecarlo.m
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

plot_montecarlo.1000.png
図 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

plot_montecarlo.10000.png
図 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

plot_montecarlo.100000.png
図 4.3.2-3 n = 100k のときのモンテカルロ法のプロット

考察

試行回数が多くなるほど、誤差は小さく、実行時間は大きくなっているようです。

積分

$D$ を下記の3つの領域($D_{1}$, $D_{2}$, $D_{3}$)に分割して、それぞれを定積分して足し合わせる方針にします。

$ octave plot_integral.m

plot_integral.png
図 4.3.3-1 定積分で面積を求めるにあたり、分割した3つの領域

数値積分

quad() を利用します。
数値積分の結果を得ます。

area_integral_numeric.m
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/数式処理システム

put_integral_symbolic.m
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() を利用します。
多角形の面積を求める関数です。

area_polyarea.m
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 で書き直そうと思います。

  1. https://qiita.com/tonkoarimoto/items/5da77c98fe33b53b3c43

  2. https://octave.org/

  3. https://jp.mathworks.com/products/matlab.html

  4. https://gnu-octave.github.io/packages/symbolic/

  5. https://gnu-octave.github.io/packages/geometry/

  6. https://inkscape.org/

  7. https://julialang.org/

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0