4.5魚の個体群の資源開発

魚の漁獲量を、魚の個体の大きさwと、魚の群れの規模nから考えたい。
個体の大きさwについて
ベルタランフィモデルによると、ある定数\alpha、\betaを用いて
\frac{dw}{dt}=\alpha w^{\frac{2}{3}}-\beta w
となる微分方程式ができる。これは、ベルヌーイの微分方程式
\frac{dy}{dx}+f(x)y=g(x)y^n
の特別な場合らしい。変数変換y^{1-n}=v、つまりv=w^{\frac{1}{3}}とすることで
\frac{dw}{dt}+\frac{\beta}{3}v =\frac{\alpha}{3}
という線型1階微分方程式に帰着できる。これは、積分因子e^{\int \frac{\beta}{3}dt}=e^{ \frac{\beta}{3}t}を用いて
w=(\frac{\alpha}{\beta})^3 (1+\frac{A\beta}{\alpha}e^{-\frac{\beta}{3}t})^3
となる。
ここで、t \to \inftyのときw\to w_{\infty}とすれば、w_{\infty}=(\frac{\alpha}{\beta})^3である。
また、w(0)=0とするとA=-\frac{\alpha}{\beta}であるから、結局
w(t)=w_{\infty}(1-e^{-\frac{\beta}{3}t})^3
となる。
　
群れの規模について
稚魚である初期状態n_0の群れが、時間経過でどうなるかは、
n(t)=n_0 e^{-Rt}
と考えていいだろう。
ここで、漁獲努力Fは行うけれども、ある時刻M才以下の魚は小さいので獲らない、とすれば、
t=Mのとき、n(M)=n_0 e^{-RM}（t=Mでの新しい初期値）
t>Mのとき、n(t)=n_0 e^{-RM}e^{-(R+F)(t-M)}（魚の減少はR\to R+Fとより早くなる）
となる。漁獲高Yを求めるには、t\leq Mをすべて足し合わせればいいから、t才で獲られてしまった魚の数をN_tとすれば
\begin{align}Y&=\int _{\hspace{5}M}^{\hspace{25}\infty} w(t)N_t dt\\&=\int _{\hspace{5}M}^{\hspace{25}\infty} w_{\infty}(1-e^{-\frac{\beta}{3}t})^3 Fn_0 e^{-RM}e^{-(R+F)(t-M)} dt\end{align}
という、FとMを変数とする関数になる。
ただし、魚の供給は考えていない。

R <- 0.003 #魚の群れの減衰定数
n0 <- 10 #魚の群れの初期値
beta <- 0.25 #魚の大きさの定数
winf <- 10000 #魚の最大の大きさ
Mmin <- 20 #漁獲していい最小の年齢

Gyokaku <- function(M, Fg){
	time <- seq(M, 100000, length=1)
	result <- winf*(1 - exp((-beta/3)*time))*Fg*n0*exp(-R*M)*exp(-(R + Fg)*(time - M))
	
	return(sum(result))
}

len <- 300
M <- seq(1, 300, length=len)
Fg <- seq(0, 10, length=len)

GyokakuRes <- matrix(0, nr=length(M), nc=length(Fg))
for(i in 1:nrow(GyokakuRes)){
	for(j in 1:ncol(GyokakuRes)){
		GyokakuRes[i, j] <- Gyokaku(M[i], Fg[j])
	}
}
GyokakuRes[1:Mmin, ] <- 0 
image(M, Fg, GyokakuRes, xlab="Fish age", ylab="Fish catch effort")

library(rgl)
col.slice <- 15 #濃度段階の設定。
colpalet <- rainbow(col.slice)[GyokakuRes/(max(GyokakuRes)) * col.slice + 1]
persp3d(GyokakuRes, col=colpalet, xlab="Fish age", ylab="Fish catch effort"
, zlab="Fish catch")