大学入試問題をoctaveで解く

## はじめに
下記の記事 [^1]（以下、元記事）にインスパイアされ、本記事を書きました。
[大学入試問題をpythonで解く #Python - Qiita](https://qiita.com/tonkoarimoto/items/5da77c98fe33b53b3c43)

ただ、元記事では解答にまで至っていなかったので、本記事では答えを求めてみます。
また、勉強を兼ねて Python から [GNU Octave](https://octave.org/) [^2]（以下、Octave）にコードを書き換えています。Octave は [MATLAB](https://jp.mathworks.com/products/matlab.html) [^3] とほぼ互換性があるらしいです。

使用している Octave のバージョンは 5.2.0 です。
symbolic パッケージ [^4] と geometry パッケージ [^5] を使用しています。
SVG 画像を PNG 画像に変換するため、[Inkscape](https://inkscape.org/) [^6] を利用しています。

## Octave の紹介
Octave の「よいところ」と「よくないところ」を簡単に紹介します。※個人の感想です
詳細を知りたい方は、下記の公式ドキュメントを参照してください。
https://docs.octave.org/latest/index.html

### よいところ
- 数値計算に適している（MATLAB と互換性がほぼある）
    - ベクトル、行列を簡単に記述できる
    - 線型代数、微分方程式、数値積分などの分野でよく利用する関数が用意されている
- プロットの GUI 表示、画像ファイル出力も簡単
- ライセンスが GPL （プロプライエタリでない）
    - https://docs.octave.org/latest/Copying.html


上記「よいところ」の一部をコードで示します。

```matlab
% 線型代数

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
```

```matlab
% 数値積分

f = @(x) 2 * sin(x); % f(x) = 2sin(x)
quad(f, 0, pi)       % F(π) - F(0)

%%% 出力
%% ans =  4
```

```matlab
% プロット

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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/73094367-2f51-9f50-a618-0c38e0dd2d41.png)
_図 2.1-1　Octave の REPL にて、グラフをプロットしている様子（OS は Ubuntu）_


### よくないところ
- ドキュメントがわかりにくい
    - 公式ドキュメントは英語
    - 関数の使用例が少ない
    - 日本語で使い方を解説してくれているサイトもあるが、情報が古い
- 実装されていない MATLAB の関数がある
    - `polyshape()` など
        - https://jp.mathworks.com/help/matlab/ref/polyshape.html

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

- ドキュメントがわかりにくい
    - → MATLAB の日本語ドキュメントを併用する
        - https://jp.mathworks.com/help/matlab/
    - → Wiki とフォーラムを漁る
        - https://wiki.octave.org/
        - https://octave.discourse.group/
    - → 公式ドキュメントとソースコードをがんばって読む
        - https://docs.octave.org/latest/index.html
        - https://github.com/gnu-octave/octave
- 実装されていない MATLAB の関数がある
    - → 代替可能な Octave のパッケージを探す
        - https://gnu-octave.github.io/packages/
    - → 自力で実装する。有料サポートを受ける
        - https://www.octave.org/missing.html


## 入試問題
```math
\begin{align}
|x| + |y| - 2 & \leq 0 \tag{1} \\
x^2 + 2x - 2y - 1 & \geq 0 \tag{2} \\
\end{align}
```
```math
\begin{array}{c}
(1), (2) を同時に満たす点(x, y)の全体からなる領域を D とする。 \\
D を図示し、その面積を求めよ。 \\
\end{array}
```
<!--
```math
\begin{split}
(1) を満たす & 点(x_1, y_1) の全体からなる領域を D_1 とし、 \\
(2) を満たす & 点(x_2, y_2) の全体からなる領域を D_2 とする。 \\
(1), (2) を同時に満たす & 点(x_3, y_3)の全体からなる領域を D_3 とする。
\end{split}
```
-->

:::note
元記事では $(1)$ が $|x| + |y| - 2 \geq 0$ となっていますが、有限の紙面では $D$ の作図が出来ず、その面積も無限大になってしまうので変更を加えています。
:::

## 解答
### 準備
#### 式変形
問題を解きやすくするため、式変形しておきます。

まずは、$(1)$ から
```math
\begin{align}
|x| + |y| - 2 & \leq 0 \\
|y| & \leq -|x| + 2\\
\end{align}
```

ここで、$y$ について場合分けします。
```math
\begin{cases}
y \leq -|x| + 2 & (y \geq 0) &&& (1a) \\
y \gt |x| -  2 & (y \lt 0) &&& (1b) \\
\end{cases}
```

さらに、$(1a)$ を $x$ について場合分けして、
```math
\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$ について場合分けします。
```math
\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)$ を変形します。
```math
\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
```

<!--
```matlab:f1_v2.m
function z = f1_v2(x, y)
  z = abs(x) + abs(y) - 2
endfunction
```
```matlab:f2_v2.m
function z = f2_v2(x, y)
  z = x.^2 + 2*x - 2*y - 1
endfunction
```

```matlab:f1_v1.m.f11a.m
%% 数式 (1.1a)
function y = f1_v1(x)
  y = -x + 2
endfunction
```

```matlab:f2_v1.m.f21.m
%% 数式 (2.1)
function y = f2_v1(x)
  y = (x.^2 + 2*x - 1) / 2
endfunction
```

```matlab:conf_axis.m.init_axes.m
%% グラフの設定
function conf_axis(limits = {-3, 3})
  [min, max] = limits{:};

  grid on;
  grid minor; % minor grid が表示されない...
  axis square;

  line([min max], [0 0], "color", "k", "linewidth", 1); % x軸
  line([0 0], [min max], "color", "k", "linewidth", 1); % y軸

  text(0, 0, "●", "fontsize", 14, "horizontalalignment", "center");
  text(max - .1, .1, "x", "fontsize", 14, "verticalalignment", "bottom", "horizontalalignment", "right");
  text(.1, max - .1, "y", "fontsize", 14, "verticalalignment", "top");

  xlim([min max]);
  ylim([min max]);

  set(gca, "xtick", min:1:max, "ytick", min:1:max);
  set(gca, "gridalpha", .8, "gridlinestyle", "-");
  % set(gca, "xaxislocation", "origin", "yaxislocation", "origin"); % 軸を原点に合わせる場合
endfunction
```
-->

#### プロット
$(2.1)$ のグラフの形状を把握するため、その右辺 $\frac{1}{2}x^2 + x - \frac{1}{2}$ に対して、二次方程式の解の公式を適用します。
```math
\frac{-b \pm \sqrt{b^2 - 4ac}}{2a} = -1 \pm \sqrt{2}
```

さらに、最小値を求めるために微分します。
```math
\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 で計算しておきます。

```matlab
% (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
```

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

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

```matlab: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
```

```shell-session
$ octave plot_ready.m
```

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

![plot_ready.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/9345beaf-4885-641e-89db-dd969e08fe77.png)
_図 4.1.3-1　plot_ready.png_

<!--
```plantuml
@startuml
<style>
root { Padding 0; Margin 0; }
caption { FontSize 28; }
</style>
title <img:https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/7f4ca885-09b9-165c-1e05-1fde357e494a.png>
caption 図１　plot_1.png
@enduml
```
-->

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

### 作図
$D$ を図示します。

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

- `fill()`
    - https://docs.octave.org/latest/Two_002dDimensional-Plots.html#index-fill
    - https://jp.mathworks.com/help/matlab/ref/fill.html

```matlab: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
```

```matlab: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
```

```shell-session
$ octave plot_d.m
```

![plot_d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/026ff947-f154-0299-6ece-a72bcb84b22f.png)
_図 4.2-1　plot_d.png_

作図は大丈夫そうです。

### 面積
#### 手計算
Octave で計算をする前に、手計算で答えを出しておきます。
下記に定義する $D_{man}$ を定積分して、残りの正方形部分の面積（$D \setminus D_{man}$）を足す方針にします。

```math
\begin{align}
x & \leq 1 \tag{3} \\
y & \geq -1 \tag{4} \\
\end{align}
```
```math
\begin{array}{c}
(2), (3), (4) を同時に満たす点(x, y)の全体からなる領域を D_{man} とする。
\end{array}
```

```shell-session
$ octave plot_manual.m
```

![plot_manual.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/540d7e9d-7e7e-cb23-224f-989c1eca1357.png)
_図 4.3.1-1　plot_manual.png_

```math
\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$ の面積は
```math
\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}$ の面積を計算しておきます。
```matlab
% 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
```

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

<!--
```math
\begin{align}
& \int^1_{\sqrt{2}-1} \frac{1}{2} (x^2 + 2x - 1) dx \\
& \qquad = \frac{1}{2}\left[\frac{1}{3}x^3 + x^2 - x \right]^1_{\sqrt{2}-1} \\
& \qquad = \frac{1}{2} \bigg\{ \Big( \frac{1}{3} \Big) - \Big( \frac{1}{3}(\sqrt{2}-1)^3 + (\sqrt{2}-1)^2 - (\sqrt{2}-1) \Big) \bigg\} \\
& \qquad = \frac{1}{6} \bigg\{(1) - (\sqrt{2}-1)\Big((\sqrt{2}-1)^2 + 3(\sqrt{2}-1) - 3 \Big) \bigg\} \\
& \qquad = \frac{1}{6} \bigg\{1 - (\sqrt{2}-1)(2 -2\sqrt{2} + 1 + 3\sqrt{2}- 3 - 3) \bigg\} \\
& \qquad = \frac{1}{6} \bigg\{1 - (\sqrt{2}-1)(\sqrt{2} - 3) \bigg\} \\
& \qquad = \frac{1}{6} \bigg\{1 - (2 - 4\sqrt{2} + 3) \bigg\} \\
& \qquad = \frac{1}{6}(4\sqrt{2} - 4) \\
& \qquad = \frac{2}{3}(\sqrt{2} - 1)
\end{align}
```

$D$ の面積は
```math
\begin{align}
D & = \frac{2}{3}(\sqrt{2} - 1) + \frac{1}{2} \\
& \simeq 0.776142375
\end{align}
```

---

```math
\begin{align}
& \int^1_{-1} \frac{1}{2}(x^2 + 2x) dx \\
& \qquad = \frac{1}{2}\left[ \frac{1}{3}x^3 + x^2 \right]^1_{-1} \\
& \qquad = \frac{1}{2}\left[\Big(\frac{1}{3} + 1\Big) - \Big(-\frac{1}{3} + 1\Big) \right] \\
& \qquad = \frac{1}{6} \\
\end{align}
```

$D$ の面積は
```math
\begin{align}
D & = \frac{1}{6} + 2 \\
& = 2.166666\ldots
\end{align}
```
---
-->

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

```matlab: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)

```shell-session
$ 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
```

```shell-session
$ n=$((10**3)); octave plot_montecarlo.m $n
```

![plot_montecarlo.1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/99a556ee-5368-e0c6-7857-c3b8e1e15fb1.png)
_図 4.3.2-1　n = 1k のときのモンテカルロ法のプロット_


##### n = 10,000 (10k)

```shell-session
$ 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
```

```shell-session
$ n=$((10**4)); octave plot_montecarlo.m $n
```
![plot_montecarlo.10000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/9c536020-069d-15d5-ffe4-8c2d9e41ee3a.png)
_図 4.3.2-2　n = 10k のときのモンテカルロ法のプロット_

##### n = 100,000 (100k)

```shell-session
$ 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
```

```shell-session
$ n=$((10**5)); octave plot_montecarlo.m $n
```

![plot_montecarlo.100000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/5088e439-0ed2-ca76-be2b-7f2cf75f3576.png)
_図 4.3.2-3　n = 100k のときのモンテカルロ法のプロット_

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

#### 積分

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

```shell-session
$ octave plot_integral.m
```

![plot_integral.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/501772/e174d717-addc-ec9e-11db-c2fb30e7e277.png)
_図 4.3.3-1　定積分で面積を求めるにあたり、分割した3つの領域_

##### 数値積分

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

- `quad()`
    - https://docs.octave.org/latest/Functions-of-One-Variable.html#index-quad
    - https://jp.mathworks.com/help/matlab/ref/quad.html

```matlab: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
```

```shell-session
$ octave area_integral_numeric.m > /dev/null
3.333333
```

##### 記号積分

`symbolic` パッケージを利用します。
数式の記号的操作を行うためのパッケージです。
https://ja.wikipedia.org/wiki/数式処理システム

- `symbolic` パッケージ
    - https://gnu-octave.github.io/packages/symbolic/
    - `@sym:int()`
        - https://octave.sourceforge.io/symbolic/function/@sym/int.html

```matlab: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
```

```shell-session
$ 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()`
    - https://docs.octave.org/latest/Voronoi-Diagrams.html#index-polyarea
    - https://jp.mathworks.com/help/matlab/ref/polyarea.html

```matlab: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
```

```shell-session
$ 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](https://julialang.org/) [^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/