4.1はじめに

\frac{dy}{dx}+P(x)y=Q(x)
となる微分方程式を線型1階微分方程式というらしい。
積分因子と呼ばれるe^{\int P(x)dx}というものを求め、これを使うと解ける。
両辺に各々かけて
e^{\int P(x)dx}\frac{dy}{dx}+e^{\int P(x)dx}P(x)y=e^{\int P(x)dx}Q(x)
となる。ところで、e^{\int P(x)dx}yをxで微分すると上式の左辺そのものだから、
\frac{d}{dx}(e^{\int P(x)dx}y)=e^{\int P(x)dx}Q(x)
と書き直せる。あとは積分することで、解は
y=e^{-\int P(x)dx}\int e^{\int P(x)dx} Q(x)dx + Ce^{-\int P(x)dx}
となる。
　
例題
\frac{dy}{dx}+xy=2x
を解く。P(x)=x、Q(x)=2xであり、積分因子はe^{\int xdx}。
e^{\frac{x^2}{2}}\frac{dy}{dx}+xe^{\frac{x^2}{2}}y=2xe^{\frac{x^2}{2}}
\frac{d}{dx}(e^{\frac{x^2}{2}}y)=2xe^{\frac{x^2}{2}}
y=2+Ce^{-\frac{x^2}{2}}

4.2広告に対する売上げ反応

商品は広告しないと売上が落ちる。
今、売上の定常的減少の仕方が直線的だったので、販売速度S、時刻t、定数\lambda、muを用いて
\log S =-\lambda t +\mu
とできる。これは、人口問題でも扱った
\frac{dS}{dt}=-\lambda S…(1)
と同じである。
　
さて、ここで広告普及率A、市場飽和度Mとして、広告の効果を考える。
飽和度は人口問題の修正や技術革新でも扱った。
売上の増減は、広告普及率Aと市場飽和率\frac{M-S}{M}に影響を受けるから、(1)式は
\frac{dS}{dt}=r A\frac{M-S}{M}-\lambda S（rは定数）
となる。整理すると
\frac{dS}{dt}+(\frac{rA}{M}+\lambda)S=rA
となる。b=\frac{rA}{M}+\lambdaとすれば、積分因子はe^{\int bdt}=e^{bt}となるので、解くと
S=\frac{rA}{b}+Ce^{-bt}
広告開始した時刻[tex=0]で売上S=S_0とすれば、C=S_0-\frac{rA}{b}なので、最終的に
S(t)=\frac{rA}{b}+(S_0-\frac{rA}{b})e^{-bt}
となる。
広告終了後は、(1)式の挙動を示す。

Stime <- 1 #観察期間
Stimes <- seq(0, Stime, length=1000)
S0 <- 100 #広告開始時の市場に出ている具合
r <- 0.1 #広告係数
lambda <- 0.001 #減少係数
A <- 100 #広告速度
M <- 1000 #市場最大数

S <- numeric(length(Stimes) * 2)
S[1] <- S0

for(i in 1:(length(Stimes) - 1)){
	dS <- r * A - (r * A / M + lambda) * S[i]
	S[i + 1] <- S[i] + dS
}

for(i in length(Stimes):(2 * length(Stimes) - 1)){
	dS <- - lambda * S[i]
	S[i + 1] <- S[i] + dS
}
S <- S / M #市場独占度に変換
plot(S, ylim=c(0, 1), type="l", xlab="time", ylab="market ratio")
abline(v=length(Stimes), lty=2, col=2)

4.3美術品の贋作

やりたいことは、放射性元素でやったことの応用である。
y(t)：時刻tにおける通常の鉛1gごとのPb^{210}の量
y(t_0)=y_0：t_0は製造開始時刻
r(t)：通常の鉛の中における毎分1gごとのRa^{226}の崩壊数
として、
\frac{dy}{dx}=-\lambda y +r(t)
とできる。積分因子はe^{\int \lambda dt}=e^{\lambda t}である。これより、
y=e^{-\lambda t}\int e^{\lambda t}r(t)dt + Ae^{-\lambda t}
となる。
ここで、取り扱っている問題の都合上、r(t)は定数とみなしてよいらしく、積分定数も考えると
y(t)=\frac{r}{\lambda}+(y_0 - \frac{r}{\lambda})e^{-\lambda (t-t_0)}
となる。

4.4電気回路

読者に任せようと言っている
R\frac{dQ}{dt}+\frac{Q}{C}=E
をやろう。
\frac{dQ}{dt}+\frac{1}{RC}Q=\frac{E}{R}
となる。積分因子はe^{\int \frac{1}{RC}dt}=e^{\frac{1}{RC}t}であるので、
Q(t)=e^{-\frac{1}{RC}t}\int e^{\frac{1}{RC}t}\frac{E(t)}{R}dt +Ae^{-\frac{1}{RC}t}
E(t)が定数なら、
Q(t)=\frac{E}{R}+(Q_0-\frac{E}{R})e^{-\frac{1}{RC}t}
である。
交流電流で、例えばE=E_0 \cos(\omega t)と与えられるとすると、
\int e^{\frac{1}{RC}}\cos(\omega t)dt=\frac{(RC)^2}{1+(\omega RC)^2}e^{\frac{1}{RC}t}\{ \frac{1}{RC}\cos(\omega t)+\omega \sin(\omega t)\} + Ae^{-\frac{1}{RC}t}
となる。時刻t=0、でQ=Q_0とすると
A=Q_0 + \frac{RC}{1+(\omega RC)^2}E_0
となるので、
E_0\frac{(RC)^2}{1+(\omega RC)^2}\{ \frac{1}{RC}\cos(\omega t)+\omega \sin(\omega t)\}+(Q_0 + \frac{RC}{1+(\omega RC)^2}E_0)e^{-\frac{1}{RC}}
となる。

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")


4.7五大湖の汚染

湖のような、ある巨大な容器の中の水が汚染されてしまい、それをどうにか浄化する作業を考える。
条件として、
・降雨量と蒸発量は互いに釣り合っている
・水が湖に流入するときは、完全な混合が起こり、汚染物は一様に分布する
・汚染物は流出によってしか、湖から除去されない
として考える。
V：湖の容積
P_l：湖の汚染度
P_i：湖への流入物の汚染度
r：流率
とすると、微小時間dtの間で全汚染物の正味の変化は
d(VP_l)=(P_i-P_l)(rdt)
となる。これは
\frac{dP_l}{dt}+\frac{r}{V}P_l=\frac{rP_i}{V}
となる線型1階微分方程式である。積分因子はe^{\int\frac{r}{V}dt}=e^{\frac{r}{V}t}であり、これを用いると
P_l(t)=e^{-\frac{r}{V}t}P_l(0)+e^{-\frac{r}{V}t}(\frac{r}{V})\int _{\hspace{5}0}^{\hspace{25}t}e^{\frac{r}{V}t}P_i dt
となる。
汚染物の流入が止まったとき、つまりP_i=0のとき
t=\frac{V}{r}\log (\frac{P_l(0)}{P_l(t)})
となる。これは、現在規定されている水準のパーセンテージまで汚染を減らすのに要する時間らしい。
t=\frac{V}{r}はRainy値と言われるものらしく、五大湖それぞれ既にある。
初期状態P_l(0)から何パーセント汚染が減るか、それにかかる時間が計算できる。