blog

木星と火星

木星を観測しやすい季節になりました.また,月末には火星が接近します.口径80 mmの屈折式望遠鏡にカメラを取り付けてリレーレンズ方式で撮影しました.

jupiter16

mars16

バイオメカニズム学会誌特集

バイオメカニズム学会誌の特集「柔軟物の硬さ」を担当しました.Vol. 40, No. 2として5月1日に発行されました.特集では,

  • 内山,永岡:押し込み型の硬度計,バイオメカニズム学会誌,Vol. 40, No. 2, pp. 97-102, (2016)

を寄稿しました.押し込み型の硬度計の計測原理と,ウレタンフォームなど様々な柔軟物を市販の硬度計で計測したときの指示値間の関係について解説しました.用いた硬度計は,アスカーゴム硬度計E型(高分子計器株式会社),PEK-1(株式会社井元製作所),TDM-Z1(有限会社トライオール),FGRT-5(日本電算シンポ株式会社),TK-03C(株式会社特殊計測)の5機種です.

_MG_3803
アスカーゴム硬度計E型
_MG_3809
PEK1
_MG_3827
TDM-Z1
_MG_3830
FGRT-5
_MG_3815
TK-03C

餌と捕食者の関係

物理情報数学Bで,行列の累乗の応用例として,餌と捕食者の関係を表す行列差分方程式を説明しています.

「第k年のイワシの数をx_1(k),クジラの数をx_2(k)とする.クジラがいないときの1年後のイワシの数を\alpha x_1(k)とする.つまり,第k年のイワシの数に比例して変化する.クジラに1年間に食べられるイワシの数を\beta x_2(k)とする.つまり,クジラの数に比例するイワシが食べられる.イワシがいないときの1年後のクジラの数を\gamma x_2(k)とする.また,イワシを食べて1年間に増加するクジラの数を\delta x_1(k)とする.」

上記の関係を行列差分方程式で表すと

   \left[\begin{array}{c}x_1(k+1) \\x_2(k+2)\end{array}\right]   =   \left[\begin{array}{rr}\alpha & -\beta \\ \delta & \gamma\end{array}\right]   \left[\begin{array}{c}x_1(k) \\x_2(k)\end{array}\right]

になります.一般解は

\boldsymbol{x}(k) = \left[\begin{array}{rr}\alpha & -\beta \\ \delta & \gamma\end{array}\right]^k\boldsymbol{x}(0)

になります.

この行列差分方程式を講義では手計算で解けるように,\alpha=1.5\beta=0.25\gamma=1.1\delta=0.07として解きます.しかし,この行列の固有値は1以上の実数が2つですので,餌と捕食者の数の変化に期待される振動的な振る舞いにはなりません.初期値のクジラの数が多いと,下図のようにイワシの数は減少に転じます.

Graph0
餌と捕食者の数(x1(0)=100,x2(0)=150)

餌と捕食者の関係のモデルとしては,ロトカ・ボルテラのモデルが知られています.このモデルでは,

   \displaystyle\frac{dx_1(t)}{dt} = \alpha x_1(t) - \beta x_1(t)x_2(t) \\   \displaystyle\frac{dx_2(t)}{dt} = \delta x_1(t)x_2(t) - \gamma x_2(t)

で表されます.非線形微分方程式ですから,x_1(t)x_2(t)を解析的に求めることはできません.適当な値を代入して数値解を求めると,x_1(t)x_2(t)が振動的な振る舞いをすることがわかります.

また,平衡点はdx_1(t)/dt=0dx_2(t)/dt=0から,原点とx_1=\gamma/\deltax_2=\alpha/\betaにあります.平衡点近傍では,x_1(t)x_2(t)の項が小さいので無視できると仮定すると,(0, 0)の平衡点近傍では,

   \displaystyle\frac{d}{dt}\left[\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right]   =   \left[\begin{array}{rr}\alpha & 0 \\   0 & -\gamma\end{array}\right]   \left[\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right]

となり,行列の固有値は実数で1つは正でもう1つは負です.したがって,(0, 0)は鞍点です.

また,もう1つの平衡点(\gamma/\delta, \alpha/\beta)の近傍では,x_1(t)x_1+\gamma/\deltax_2(t)x_2+\alpha/\betaとおいて,x_1(t)x_2(t)の項が小さいので無視できると仮定すると

   \displaystyle\frac{d}{dt}\left[\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right]   =   \left[\begin{array}{rr}0 & -\frac{\beta\gamma}{\delta} \\   \frac{\alpha\delta}{\beta} & 0\end{array}\right]   \left[\begin{array}{c}x_1(t) \\x_2(t)\end{array}\right]

と表されます.したがって,固有値が純虚数となり,振動する解になります.

\alpha=10\beta=0.1\gamma=10\delta=0.1とし,x_1(0)=10x_2(0)=100 で数値解を求めると下図のようになります.上のパネルはx_1(t)x_2(t)の時間発展です. 下のパネルはx_1(t)x_2(t)の関係です.

x1(t)とx2(t)の時間発展
lv2
x1(t)とx2(t)の関係(t=0で(10, 100)とx1=0のヌルクラインに近く,x2が急速に減少するが,原点が鞍点であるのである程度減少したところでx1が急速に増加する.)

初期値を平衡点((100, 100))に近い点(99, 99)とすると下図のように円に近い軌道を描きます.

lv3
初期値を(99, 99)と平衡点(100, 100)に近い値にしたので円に近い軌道を描く

MATLABのコードです.

  • lv.m
tspan = [0:0.001:5];
[t y] = ode45(@LotkaVolterra,tspan, [99 99]);
figure(1);
plot(t, y(:,1), 'b', t, y(:,2), 'r');
figure(2);
plot(y(:,1), y(:,2));
  • LotkaVolterra.m
function dx = LotkaVolterra(t, x)
alpha = 10;
beta = 0.1;
gamma = 10;
delta = 0.1;
dx = zeros(2,1);
dx(1) = alpha * x(1) - beta * x(1) * x(2);
dx(2) = delta * x(1) * x(2) - gamma * x(2);
end

第55回日本生体医工学会大会

富山で開催されました第55回日本生体医工学会大会にて

  • K. Saito, T. Uchiyama: ‘Estimation of muscle stiffness during one cycle of pedaling exercises using system identification’
  • K. Suzuki, T. Uchiyama: ‘Influence of footwear on stiffness of tibialis anterior muscle during walking’
  • T. Kato, T. Uchiyama: ‘System identification of mechanomyogram in voluntary contraction of abductor digiti minimi muscle’

以上,3件を発表しました.

入学式

本日は慶應義塾大学の入学式です.矢上キャンパスにも塾旗が掲揚されています.

_DSC8226

卒業式

本日は卒業式でした.日吉の記念館で卒業式の後,矢上キャンパスで学科に分かれて学位記を授与されます.

_MG_5694

研究室の卒業生5名の内,3名は四月から社会人です.

_MG_5576