昆虫を数学する(入門)

目次

1.はじめに
2.準備
3.昆虫の数学
4.歴史・研究者
5.最後に
6.謝辞
この記事は Mathlog Advent Calendar 2023 大学数学部門 20日目の記事です。

1.はじめに

どうも、色々やる数学徒です。
今回の記事では「昆虫」にまつわる数学のお話を紹介していこうと思います。(昆虫にまつわる数学っていうよりは数理モデルのガイドみたいな感じになると思います)
僕は昔から昆虫と数学が大好きなのですがその両方が融合した興味深い分野の歴史や詳細を書いていきます。

※この記事は何か新しい試みを紹介するというよりは現在までの偉大なる先人たちの研究結果の紹介記事のようなものとなっています。
※昆虫の画像がしばしば登場いたします。苦手な方はご注意ください。

2.準備

数理モデルの世界に入るにあたって必要不可欠な知識を書いていきます。(初学者でもわかりやすく読めるように努力します)
分かる方は飛ばしてもらって構いません。

数理モデル

簡単に言ってしまえば数理モデルとは、対象となる事象(現実世界で起こるさまざまな現象)を数学の言語を用いて記述することです。(数理模型と呼ぶこともあるそうです)
この準備ではモデルを数学の言語で書き起こすための作業で重要になってくる基礎的な技術を記します。

記事中の記法・用語・注意点について

必ずこれを目に通してからお読みください

登場する物理・生物・化学の公式

以下の公式は既知として扱っちゃいます。(ここでは詳しく説明とかはしません)

すべての基本となっている超重要な式ですね。

その他高校レベルの物理公式

微分方程式

微分方程式とは中学・高校でやる方程式($x^2+x+1=0$のような)とは大きく異なります。
中学・高校でやる方程式は等式を満たすような$x$の値を見つけることを方程式を解くと言いましたが、微分方程式ではそのような条件を満たす関数を見つけなければなりません。
微分方程式は大きく分けて$2$種類です。それは「常微分方程式」と「偏微分方程式」です。($1$変数と多変数の違いです)

常微分方程式

具体的に常微分方程式では以下のような方程式を取り扱います。

<解き方>
$\displaystyle y’=Ay$
$\displaystyle A=\frac{1}{y}y’=\frac{d}{dx}\ln y(x)$
$\displaystyle \ln y(x)=\int \left[\frac{d}{dx}\ln y(x)\right]dx=\int Adx=Ax+C$($C$は積分定数)
$\displaystyle y(x)=e^{Ax+C}=e^Ce^{Ax}=De^{Ax}$
このような微分方程式の解を一般解と言いますが、初期条件(ある一点の関数の値)を与えると任意定数$D$が得られそれを特殊解と言います。
特殊解の具体例
$\displaystyle y’(x)=Ay(x),y(0)=1⇒y=e^{Ax}$

常微分方程式の解き方

微分方程式を解くにあたってスケール変換という用語が登場するので少し紹介しておきます。
スケール変換とは、ある特定の数$j$と$k$があって、$c$を任意の定数とするとき、$x→c^jx,y→c^ky$のように変換することをスケール変換といいます。常微分方程式がこのスケール変換に対しても不変であれば解をスケール変換したものも解になりますね。($x→x+c$みたいな変換は並進変換というらしい)
一般に物理に登場する方程式はこのような不変性を持ちます。そして、不変性に対応する保存則が成立します。(ネーターの定理)

力学の運動方程式では慣性の法則が関わってくるため$2$階微分方程式になってしまいますが簡単な法則であれば$1$階になることが多いです。

見たらわかる通り$x$について微分しているだけの形なので両辺を積分すればいいですね。

例)$\displaystyle y’=3(x+4)$
<解き方>
両辺に$dx$をかけます.
$\displaystyle y=\int 3(x+4)dx$
$\displaystyle =3\left(\frac{x^2}{2}+4x\right)+C$

<解き方>
右辺の邪魔な$g(y)$を両辺で割ってキレイにします。
$\displaystyle \frac{1}{g(y)}\frac{dy}{dx}=f(x)$
$x$について積分します。
$\displaystyle \int\frac{dy}{g(y)}=\int f(x)dx$
もちろん不定積分なので積分定数はでてきてしまいますが、両辺の積分定数をまとめちゃって大丈夫です。
また、両辺を$g(y)$で割っちゃったわけですが$g(y)=0$についても考えねばなりませんね。
$g(y)=0$になる点を$y=a$とおくと上ではその場合を除外していることになりますね。$y=a$であった場合、両辺が$0$ということで微分方程式を満たします。この解は一般解に含まれていませんが、一般解の任意定数の極限値として再現できるので特異解とは言いません。

例)$\displaystyle y’=6x^2y$
<解き方>
$\displaystyle y’=6x^2y⇒\frac{1}{y}y’=6x^2$($y≠0$)
両辺に$dx$をかける
$\displaystyle \int\frac{1}{y}dy=\int6x^2dx$
$\displaystyle \log|y|=2x^3+C$
$\displaystyle |y|=e^{2x^3+C}⇒y=±e^C×e^{2x^3}$
($±e^C$をAとおく)
$\displaystyle y=Ae^{2x^3}$
忘れてはいけないのが$y=0$のときです。
$\displaystyle y’=0$となり微分方程式を満たします。

まず、同次スケール変換とは$j=k=1$のスケール変換のことです。$\displaystyle \frac{dy}{dx}$はこの変換に対して不変なので上のように書けるわけです。$y=xv$という変換を行うと$xv’+v=f(v)$となります。
これは$\displaystyle v’=\frac{1}{x}(f(v)-v)$に帰着します。

$f(x)=0$ならば変数分離形であるため$\displaystyle y=\exp\left(-\int p(x)dx\right)$
ここで$\displaystyle y=v\exp\left(-\int p(x)dx\right)$という変数変換を行う。
$\displaystyle P’(x)=-p(x)P(x)$

$\displaystyle P(x)v’=f(x)$

したがって、$\displaystyle v=\int\frac{f(x)}{P(x)}dx$

$\displaystyle y=P(x)\int\frac{f(x)}{P(x)}dx$
よって解は$\displaystyle =\exp(-\int p(x))・\int[f(x)\exp(\int p(x)dx)]dx$となります。

$z=F(x,y)$の全微分は$\displaystyle dx=\frac{\partial F}{\partial x}dx+\frac{\partial F}{\partial y}dy$
微分方程式を$\displaystyle \frac{\partial F}{\partial x}dx+\frac{\partial F}{\partial y}dy=0$のように表示できれば解は$F(x,y)=C$となります。
このような微分方程式を完全微分方程式といいます。
与えられた微分方程式$P(x,y)dx+Q(x,y)dy=0$が完全微分方程式であるための必要十分条件は$\displaystyle \frac{\partial P(x,y)}{\partial y}=\frac{\partial Q(x,y)}{\partial x}$です。
例) $(4xy+9x^2)dx+(x^2+2y)dy=0$
解は$2x^2y+4x^3+y^2=A$なります。

$y’$に対して明示的に解かれていない微分方程式のことを非正規型微分方程式といいます。
最も簡単な例としてクレーローの微分方程式が挙げられます。$y=y’x+f(y’)$
この一般解は$y=Cx+f(C)$です。

例) $f(y,y’,y’’)=0$を考えてみましょう。
$\displaystyle y’=\frac{dy}{dx}=\left(\frac{dx}{dy}\right)^{-1}=\frac{1}{x’}$
$\displaystyle y’’=\frac{dy’}{dx}=\frac{dy}{dx}\frac{d}{dy}\left(\frac{1}{x’}\right)=\frac{1}{x’}(-1)\frac{x’’}{x’^2}=-\frac{x’’}{x’^3}$
より$\displaystyle f\left(y,\frac{1}{x’},-\frac{x’’}{x’^3}\right)=0$に帰着できます。

例) $f(y,xy’,x^2y’’)=0$
<解き方>
$x=e^u$とおく
$\displaystyle f\left(y,\frac{dy}{du},\frac{d^2y}{du^2},-\frac{dy}{du}\right)=0$となる
また、$x→cx$と$y\to cy$を同時に同次スケール変換のもとに不変な微分方程式$\displaystyle f(x^{-1}y,y’,xy’’,\cdots,x^{n-1}y^{(n)})$は$y=xv$とおけば$x$の1変数スケール変換不変な場合に帰着する.

ここで$A_n(x),…,A_0(x)$は$x$のみの関数であり、この解をいくつか互いに線形独立なものがあるならばそれらの線形結合も解となります。

$n$次元の縦ベクトル$\mathbf y:=(y,y’,\cdots,y^{(n-1)})^t$と$n$次元の横ベクトル$\mathbf p:=(p_n,p_{n-1},\cdots,p_1)$とおけば上の微分方程式は$y^{(n)}+(\mathbf p・\mathbf y)=0$と書けますね。
この微分方程式に対して重ね合わせの原理が成り立つため一般解は$y=C_1y_1+C_2y_2+\cdots+C_ny_n$となります。

$n$個の$n$次元ベクトル$\mathbf y_1,\cdots ,\mathbf y_n$を成分表示してできるロンスキー行列式という次のような行列式を考えます。
$\displaystyle W(y_1,\cdots,y_n):=\textup{det}(\mathbf y _1,\mathbf y_2,\cdots,\mathbf y_n)$
$:=\begin{vmatrix} y_1 & y_2 & \cdots &y_n \\ y_1’& y_2’&\cdots &y_n’ \\ \vdots&\vdots&\ddots&\vdots\\ y_1^{(n-1)}&y_2^{(n-1)}&\cdots&y_n^{(n-1)} \end{vmatrix}$
この$n$個の解が一次独立であるための必要十分条件はロンスキー行列式が$0$でないことです。
$y_k^{(n)}=-p_1y_k^{(n-1)}-p_2y_k^{(n-2)}-\cdots-p_{n-1}y_k’-p_ny_k$
$\displaystyle \frac{d}{dx}W=-p_1W$
$\displaystyle W=\exp\left(-\int p_1(x)dx\right)$
こんなふうに表示できるってわけですね。

実際に$2$階常微分方程式を解いてみましょう。
$2$階以上になるとちょっと工夫しなければなりません。(これはもうテクニックを覚えるしかないです)

まず、$\displaystyle y’’+f(x)y’+g(x)y=0$を考えます。

解の形を$y=e^{\alpha x}$と仮定して、上の式の各項に代入します。
$\displaystyle y’’=\alpha^2e^{\alpha x},f(x)y’=f(x)\alpha e^{\alpha x},g(x)y=g(x)e^{\alpha x}$

$\displaystyle (\alpha^2+f(x)\alpha+g(x))e^{\alpha x}=0⇒\alpha^2+\alpha f(x)+g(x)=0$
($e^{\alpha x}>0$より上のようになります)
また、上の式は$\alpha$の$2$次方程式であるため次のように解がえられます。
$\displaystyle \alpha=\frac{-f(x)±\sqrt{f(x)^2-4g(x)}}{2}$
$\displaystyle \alpha_1=\frac{-f(x)+\sqrt{f(x)^2-4g(x)}}{2},\alpha_2=\frac{-f(x)-\sqrt{f(x)^2-4g(x)}}{2}$

$\displaystyle y(x)=c_1e^{\alpha_1 x}+c_2e^{\alpha_2 x}$

要は$n$階の定数係数線形微分方程式は$n$次方程式に帰着するってわけですね。
上の解説文の中に語弊を招く記述が見られたので読み飛ばしてください(今後訂正いたします)
例)$\displaystyle \frac{d^2y}{dx^2}=-y$
<解き方>
$\displaystyle y=c_1\sin x+c_2\cos x$
$\displaystyle y=C\sin(x+\alpha)$
(三角関数の合成を用いれば特殊な場合の微分方程式もワンパターンで書けることがしばしばあります)

偏微分方程式

基礎的な物理や数理モデルでは常微分方程式で話は終わりますが、複雑なモデルになってくると話は変わりってきます。
例)$\displaystyle \frac{\partial u}{\partial x}-\frac{\partial u}{\partial y}=x^2-y^2$
都合の良い形に移項する。
$\displaystyle \frac{\partial u}{\partial x}-x^2=\frac{\partial u}{\partial y}-y^2=a$
積分する
$\displaystyle u=\int(x^2+a)dx+\int(y^2+a)dy$
$\displaystyle \;\;\;=\frac{1}{3}(x^3+y^3)+a(x+y)+b$

↑
こんな感じのものを偏微分方程式って言います。(解ける時の方が少ないですね)
実用例を紹介しておきます。
熱伝導方程式の導出です。

フーリエ級数展開

数学を学ぶ最初の段階に遭遇すれば大興奮不可避ですね。
区間$[-\pi,\pi]$で定義された有界で積分可能な関数$f(x)$を考えます。
このとき、$\displaystyle a_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\cos nxdx$,$\displaystyle b_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\sin nxdx$とし$\displaystyle S(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}(a_n\cos nx+b_n\sin nx)$で定義される級数を$f(x)$のフーリエ級数といい$f(x)$と$f’(x)$が区分的に連続なとき$f(x)$は区分的になめらかであるといいます。

フーリエ変換

上のフーリエ級数展開を周期的でない関数にも適用できるように導入されたのがフーリエ変換です。
完全正規直交関数列$\displaystyle u_n(x)=\frac{1}{\sqrt{2l}}\exp\left({\frac{in\pi x}{l}}\right)$,$n\in\mathbb Z$
を考えます。
$\displaystyle f(x)=\frac{1}{\sqrt{2l}}\sum_{n=-\infty}^{\infty}C_n\exp\left(\frac{in\pi x}{l}\right)$
ここで$\displaystyle C_n=(f,u_n)=\frac{1}{\sqrt{2l}}\int_{-l}^lf(x)\exp\left(\frac{-in\pi x}{l}\right)dx$
(以降$\displaystyle F_n=\frac{1}{\pi}C_n$というものを考える)
すると、
$\displaystyle f(x)=\frac{1}{\sqrt{2\pi}}\sum_{n=-\infty}^{\infty}\frac{\pi}{l}F_n\exp\left(\frac{in\pi x}{l}\right)$
$\displaystyle F_n=\frac{1}{\sqrt{2\pi}}\int_{-l}^lf(x)\exp\left(\frac{-in\pi x}{l}\right)dx$
$l\to\infty$を考えてみましょう。

$\displaystyle k_n:=\frac{n\pi}{l}$($n\in \mathbb Z$)
$\displaystyle f(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}F(k)\exp(ikx)dk$

$\displaystyle F(k)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(x)\exp(-ikx)dx$
したがって$l\to\infty$としたとき
$\displaystyle f(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}F(k)\exp(ikx)dk$となります。
また、$\displaystyle \frac{F(k)}{\sqrt{2\pi}}=\mathcal F(k)$とおけば

$\displaystyle f(x)=\int_{-\infty}^{\infty}\mathcal F(k)\exp(ikx)dx$
$\displaystyle \mathcal F(k)=\frac{1}{2\pi}\int_{-\infty}^{\infty}f(x)\exp(-ikx)dx$
と簡潔に表せます。
特にこの記事内でいうフーリエ変換は以下で定義したものです。(上はただの導入です)

ラプラス変換

具体的に色々な関数をラプラス変換していきましょう。

実際にラプラス変換を用いて微分方程式を解いてみましょう。

例1)$3y’+4y=0,y(0)=1$
<解き方>
両辺ラプラス変換します。
$3(sY-y(0))+4Y=0$
$(3s+4)Y=3$
$\displaystyle Y=\frac{3}{3s+4}=\frac{1}{s+\frac{4}{3}}$
$\displaystyle y=e^{-\frac{4}{3}t}$

例2)$y’-6y-2=0,y(0)=1$
<解き方>
両辺ラプラス変換する
$\displaystyle sY-y(0)-6Y-\frac{2}{s}=0$
移行して初期値を代入する
$\displaystyle sY-1-6Y=\frac{2}{s}$
$\displaystyle (s-6)Y=\frac{2}{s}+1$
$\displaystyle Y=\frac{2}{s(s-6)}+\frac{1}{s-6}$
$BBB$(部分分数分解)します
$\displaystyle Y=\frac{2+s}{s(s-6)}$
これを逆ラプラス変換します。
$\displaystyle y=\frac{1}{3}(4e^{6x}-1)$

$2$階常微分方程式も同じようにして解くことができます。

積分方程式

正直積分方程式は微分方程式ほど種類も多くなく解法も簡単なのでいくつか例を示すだけにとどまります。(追記12/20:ここで触れる積分方程式は簡易なものだけ、という意図でしたが上手く表現することができていませんでした。)

例)$\displaystyle f(x)=\int_0^5tf(t)dt+x^2$を満たす関数$f(x)$
ここで気づくべきことは$\displaystyle \int_0^5tf(t)dt$が$x$の関数ではなく定数だという点です。$A$とでも置いておきましょう。
そうすると$\displaystyle f(x)=x^2+A$とかけ、
積分内の$f(t)$にこれを代入すると解が導けます。
$\displaystyle A=\int_0^5t^3+At$
$\displaystyle A=\left[\frac{t^4}{4}+\frac{t^2}{2}A\right]_0^5$
$\displaystyle A=-\frac{625}{46}$
$\displaystyle f(x)=x^2-\frac{625}{46}$
ちょっと数字が汚くなっちゃいました。

例)$\displaystyle \int_0^xf(t)dt=x^2+a$を満たす関数$f(x)$と定数$a$
両辺微分します。
$f(x)=2x$が得られ元の式に$x=0$を代入します。
$a=0$が得られました。

反応拡散系

ここに形は載せていますが、詳しく解説しちゃうとこの記事がとてつもなく長くなってしまうので個人的にわかりやすかった解説サイトを載せておきます。(追記12/18:結局書きました)
解説サイト

具体的な数理モデル

上で述べた考え方を用いて具体的な数理モデルを考えましょう。(有名なやつばかりで面白くないのでとばしてもらって構いません)

単振動

フックの法則や運動方程式に従うと次のように立式できます。
$-ky=my’’$
初期値$y(0)=0$より
$\displaystyle y=A\cos \sqrt{\frac{k}{m}}t$

マルサスの人口モデル

微分というものの本質を考えれば容易に以下が導けます。
$\displaystyle \frac{dP}{dt}=mP$

ロジスティック方程式

このロジスティック方程式は見てわかる通りマルサスモデルの強化版みたいなものです。
ここで人口の増加だけでなく病気などに感染した数も考えましょう(全体の人口$N$は一定)
このとき、切片$k$,$N$のときは$0$になるような比例係数を考えればいいです。
したがって比例係数を$K$とでも置けば
$\displaystyle K=k\left(1-\frac{y}{N}\right)$
とかけます。
$y$の変化率が$y$に比例するとすれば
$\displaystyle \frac{dy}{dt}=k\left(1-\frac{y}{N}\right)$
この微分方程式をロジスティック方程式というそうです。(数理モデルで一番有名ですね)
$\displaystyle $

ロトカ・ヴォルテラの方程式

長くなってしまうのでここでは詳しく途中式は載せませんが以下のような微分方程式となることが知られています。
$\displaystyle \frac{dx}{dt}=ax-bxy$
$\displaystyle \frac{dy}{dt}=cxy-dy$
$x$は被食者の個体数、$y$は捕食者の個体数、$t$は経過した時間を表します。($a,b,c,d$は実数のパラメータ)

三体問題(おまけ)

物理だとこんな問題もあります。(特殊な場合の解は有名ですね)
他にもSIRモデルなど実用的な式は大体数理モデルです。

3.昆虫の数学

では数理モデルについてある程度理解が深まったと思うので今回の本題、「昆虫の数学」についていくつかのテーマを紹介しようと思います。

素数ゼミ

みなさんは素数ゼミという一風変わったセミを知っていますか？このセミはなんと驚くべきことに素数年周期で地上に出てきます。($13$,$17$年周期の$2$種類らしい)
近年、なぜこのセミは素数年周期なのかという研究が進んでいます。(数理モデルを用いた)

昆虫の目

こちらはかなり面白い話ですよ。
知っている方も多いかもしれませんが、ここではチューリング・パターンについて記します。
提供:竹さん
上の図$1$はアブの$1$種なのですが、複眼を少し拡大したものを見てみましょう。 拡大された複眼
図$2$の複眼を見たらわかるように動物の角膜にも何か数学的性質があるように思えませんか？
このように生物に現れる模様の規則性は拡散方程式を用いて表現できます。
まず、拡散方程式を導入する際に重要となってくるのは生物の模様も化学反応の結果だと考えることです。
このとき、ある化学反応を考える上で登場する$n$種類の物質とその濃度を考えます。($u_i$とおきます。$i=1,2,…n$)この$n$種の物質の反応速度を$\displaystyle \frac{\partial u_i}{\partial t}=f_i(u_1,…,u_n)$と書くことにします。
また、$n$種の物質の拡散方程式は$\displaystyle \frac{\partial u_i}{\partial t}=D_i \Delta^2u_i$と書けます。
ここで$D$は各物質の拡散係数を対角成分に持つ対角行列です。
そして反応も計算に含み考える場合は$\displaystyle \frac{\partial u_i}{\partial t}=f_i(u_1,…,u_n)+D_i\Delta^2u_i$
これを反応拡散方程式といいます。この方程式を基としてパターンが生じる条件や実際に生じるパターンを考えるそうです。

次に$Turing$不安定性というものを導入します。
$Turing$不安定性とは「拡散がないときに一様定常状態が安定であるような状況」のことを指すそうです。このとき濃度の時間発展方程式は$u_t=\gamma f(u,v)$などと書けます。(参考文献でパラメータを$\gamma$とおいてあったのでそれし従います)
この$Turing$不安定性によって生成されるパターンを$Turing$パターンと呼びます。
$\displaystyle \frac{\partial u}{\partial t}=\gamma f(u,v)+\Delta^2u$

$\displaystyle \frac{\partial v}{\partial t}=\gamma g(u,v)+d\Delta^2v$
を満たす$2$種物質反応拡散方程式を考えます。
さらにここで考えている領域$B$において次のような条件を課します。(ゼロフラックス境界条件というらしい)
$\mathbf n\in\partial B$において$\displaystyle (\mathbf n・\Delta)\mathbf u=0$,$\mathbf u=\begin{bmatrix}u\\v\end{bmatrix}$
ここで$\partial B$は領域$B$の境界、$\mathbf n$は法線ベクトルです。参考文献によるとこの条件は「領域の端で物質の出入りがないよ」ということを意味しているそうです。
($Turing$パターンも良いですが、拡散がある場合を考えているのでした)
以上の前提を踏まえて体表パターンについて入りましょう。
シュミレーションをするためには上で述べたように具体的な化学反応が必要なのですが、ここでは参考文献通りに$Thomas$の系という化学反応系を用います
$\displaystyle f(u,v)=a-u-h(u,v),g(u,v)=\alpha(b-v)-h(u,v),h(u,v):=\frac{ρuv}{1+u+Ku^2}$
この化学反応の反応拡散方程式は、
$\displaystyle \frac{\partial u}{\partial t}=\gamma f(u,v)+\Delta^2 u,$
$\displaystyle \frac{\partial v}{\partial t}=\gamma g(u,v)+d\Delta^2v.$
と書けます。
この系を用いるパターン形成のシミュレーションは次のようになるらしいです。
$\displaystyle \mathbf u_{i,j}^n:=\mathbf u(i\Delta x,j\Delta y,n\Delta t)(i=0,1,2,…,N_x-1,j=0,1,…,N_y-1)$
$\displaystyle u_{i,j}^{n+1}=u_{i,j}^n+\Delta t\gamma f(u_{i,j}^n,u_{i,j}^n)+\frac{\Delta t}{(\Delta x)^2}(u_{i+1,j}^n+u_{i,j+1}^n+u_{i-1,j}^n+u_{i,j-1}^n-4u_{i,j}^n),$
$\displaystyle v_{i,j}^{n+1}=v_{i,j}^n+\Delta t\gamma g(u_{i,j}^n,u_{i,j}^n)+\frac{\Delta t}{(\Delta x)^2}d(v_{i+1,j}^n+v_{i,j+1}^n+v_{i-1,j}^n+v_{i,j-1}^n-4v_{i,j}^n).$
ここで$x$方向に対するゼロフラックス境界条件は、$\mathbf u_{-1,j}^n:=\mathbf u_{0,j}^n,\mathbf u_{N_x,j}^n:=\mathbf u_{N_x-1,j}^n.$
周期境界条件は、
$\displaystyle \mathbf u_{-1,j}^n:=\mathbf u_{N_x-1,j}^n,\mathbf u_{N_x,j}^n:=\mathbf u_{0,j}^n$となる。
どのように導出するのかは参考文献とした こちら をご覧ください。
しかし、これを実際にシミュレートするのは大変な労力であったため基本的な反応拡散方程式のパラメータをイジってアブの複眼っぽいものを生成してもらいました。(僕はPythonすら扱えないのでうぃりあむさんに丸投げしてしまいました。笑)
シミュレートしたものが次の図です。
うぃりあむさんに作ってもらった画像
実際にアブの複眼と似ていますね。(このシミュレートを適応させるには生物学的な制限があるらしいですが僕はあまり理解が及んでいないので気になる方は「チューリンング・パターン」で調べてみてください)おそらく複眼の水晶体の境界条件や生物学的にモルフォゲンが〜という細かい制約がたくさん必要みたいです。改めてこの分野を開拓している方々の凄さがわかりました。

翅

この話はあまり有名じゃないかもしれません。
皆さんはトンボなどの不均翅亜目の仲間の翅の構造は幾何学的だと思ったことはありませんか？

提供:竹さん
実はトンボの翅脈はボロノイ構造に酷似しているのです
ボロノイ構造の定義をしておきます。(wikiより)

距離空間$(X,d)$ 内の有限な部分集合 $P\in X$が与えられたとき、各点 $p\in P$ を母点またはサイトと呼び、これに対して、$X$の中で「$P$の点の中で$p$が最も近い」点の集合
を$p$の(ボロノイ)領域と呼び、$P$ の全ての点の領域を集めた集合（の誘導するセル複体）をボロノイ図と呼ぶ。
ボロノイ領域の境界をボロノイ境界と呼び、各々のボロノイ境界の交点をボロノイ点と呼ぶ。

このボロノイ構造がトンボの翅脈などの自然界で現れるのはその性質の良さなどの工学的な理由があるみたいです。
さらに面白い話としてハネカクシやテントウムシの翅の構造も工学的観点から研究されているそうです。→ これに関する論文

はばたき

ここでは昆虫が飛ぶときにどのように数学的に記述できるかについて述べるのですが、実は近年までマルハナバチのように胴体が翅に比べて大きい種類は流体力学的に解明できていませんでした。しかし、近年流体力学などの活用により研究が加速度的に進歩しているようです。
はばたきに関する論文
僕はまだ全然流体力学のことを理解していないのでここでは詳しく紹介できなそうです。(基礎でいきなり物体の大きさを考えない、みたいな話が出てきて全く理解できませんでした笑)
提供:竹さん

個体群

次に社会性昆虫のコロニー内でのカーストなどの数理モデルを紹介します。まず、社会性昆虫とは

社会性昆虫（しゃかいせいこんちゅう、英語: social insect）とは、ハチやシロアリのように、集団を作り、その中に女王や働き蟻（蜂）のような階層があるような生活をしているなど、人間のそれに似た社会的構造を備える昆虫を指す(Wikipediaより)

人間のように集団で生活する昆虫を指します。
提供:高橋さん
提供:高橋さん
以下のようなモデルが提案されているっぽいです。

ここで${}^jp_i$は個体$j$の 表現型$i$を，${}^jq_i$は個体$j$の幼若ホルモン $i$の濃度を表しています。
これらも反応拡散方程式として考えられているようです。
詳しくは 参考にした論文 を参照してください。(僕もまだ全体理解が及んでいません)

寄生に関する数理モデル

では、ついに最後のテーマです。
冒頭で僕は昆虫・数学が好きだと言いましたが、実はもう一つ好きな分野があります。そうです！僕のアイコンにもなっている寄生蜂です。(寄生虫ほど高等な生物はいないでしょう)
調べてみると、そんな寄生蜂にも数理モデルがあるそうなのです！
おそらくイチジクコバチの仲間
一風変わったアシブトコバチの仲間
上の画像の個体たちは僕が採集した中で特にお気に入りのやつだったので少し自慢。キレイですよね！
昆虫の寄生関係では次のようなモデルが成り立つそうです。(昆虫は繁殖が同期しているので離散時間モデルを用いるそう)

($H_t$は時刻$t$におけるホスト個体密度、$P_t$も同様にパラサイト個体密度、$f$はホストが寄生を免れる確率、$R,c$はそれぞれホストの増殖率、寄生し成長したパラサイト数)
Twitterのフォロワーさんによると寄生蜂などの産卵管の長さもモデル化できるそうです。

4.歴史・研究者

歴史

最初は数理モデルの歴史は案外浅いのではないかと思っていたが、考えてみたらそもそも物理学などの分野は事象をモデル化して数学の言葉で述べているのだから浅いわけがありませんでした。しかし、昆虫などのモデル化というのは最近ようやく研究され始めたのでこれから発展していくことでしょう。

研究者

内田　俊郎
合原　一幸
大崎 浩一
久保田　耕平
宮竹　貴久
アラン・チューリング
高須　夫悟
その他にも僕の知識不足で知らない巨人たちがいるのでぜひ皆さん一度は調べてみてください。

5.最後に

お疲れ様でした〜:;(∩´﹏`∩);:
今までの僕の投稿した記事の中で最長になってしまいました。
高校の物理なども微分方程式の知識があるとある程度説明がついたりするのでこれを機にちょっと考えてみるのもありですね。(逆に力学的エネルギー保存則とかはバリバリの大学物理レベルの証明なので感謝だけして勝手に使っちゃいましょう、先生に質問したらめちゃくちゃ渋い顔されました笑)
この記事をきっかけに数学徒も数学初学者も数理の世界を楽しんで頂けたのなら光栄です。
今回の記事は本当に色々な人に支えられてようやく作成することができました。(特に昆虫チャットの皆様、塵芥さん、うぃりあむさん)ありがとうございました。
提供:高橋さん
追記12/18 保存則(経路が違っても変わらないこと)を使って証明するそうです。ちなみに運動方程式は経験則らしいです。(＞人＜;)(詳しく調べてないからそうではないかも)

6.謝辞

本記事の作成にあたり、多くの方々にご指導ご鞭撻を賜りました。物理のT先生、並びに代数のT先生、本記事の作成にあたり、適切なご助言を賜りました。感謝申し上げます。
昆虫チャットのみなさまには、多くの情報と写真を提供頂きました。厚く御礼申し上げます。
最後に、友人・先輩方には本記事の遂行にあたり多大なご助言、ご協力頂きました。ここに誠意を表します。