微分方程式で数学モデルを作ろう　2章

２章の「薬の吸収」と「水の過熱と冷却」についてmatlabで書いてみました。
まずは「薬の吸収」からです。

function dYdt=medicine_model(t,Y)
k=0.5
dYdt=[-k*Y];

%一定の濃度幅を維持し続ける場合
k=0.5;
y0=100;
T=0.5;
Ys=y0/(1-exp(-k*T));
Ym=Ys*exp(-k*T);
for i=[0:1:29];
    [t,Y]=ode45(@medicine_model,[i*T (i+1)*T],Ys);
hold on;
plot(t,Y,'b');
title('薬剤投与モデル');
xlabel('時間（day）');
ylabel('投与量(mg)');
axis([0,15,0,500]);
set(gca,'XTick',[1:1:15]);
plot([i*T i*T], [Ym Ys], 'r');
%一定の投与量を維持し続ける場合
Yi=y0*sum(exp([0:1:i]*(-k*T)));
[t,Y]=ode45(@medicine_model,[i*T (i+1)*T],Yi);
plot(t,Y,'g');
plot([i*T i*T], [(Yi-y0) Yi], 'r');
hold off;
end;
saveas(figure(1),'medicine_model','png');

次に「水の過熱と冷却」についてです。
液体の入った容器をヒーターを使って温めるケースを想定しています。
液体の比熱、質量をcw,mw,容器のをcc,mc、さらにC(t):時刻tにおける周囲との温度差とすると、
容器が時刻t=0から得た総熱量h(t)について
h(t)=(cw*mw+cc*mc)*C(t)　　　　　・・・・・・(1)
と表せます。
さらに、j:伝熱係数、a:容器の表面積とすると、
\frac{dh(t)}{dt}=-jaC(t)　　　　　・・・・・・(2)
と表せます。
(1)の微分形と(2)から
\frac{dC}{dt}=[\frac{-ja}{cw*mw+cc*mc}]C　　　　　・・・・・・(3)
と導けます。
これから、
C(t)=C(0)*e^{\(-k_{1}*t\)} \(k_{1}=[\frac{-ja}{cw*mw+cc*mc}]\)
となります。
本では初期値45から何もせず冷まして、最後に一気にヒーターで加熱した場合と、42まで下がったら45まで加熱する場合を考え、
どちらのほうが熱量が少なくてすむかを計算しています。
結果は一気に加熱した場合のほうが安くすむそうです。

function dCdt=water_model(t,C);
a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
dCdt=-k1*C;


a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
[t,C]=ode45(@water_model,[0,9],45);
plot(t,C,'b');
title('水の加熱と冷却');
xlabel('時刻');
ylabel('周囲より高い温度');
axis([0,15,0,50]);
%最後に温める
hold on;
Cm1=45*exp(-9*k1);
x=9:14;
y=12*(x-9)+Cm1;
plot(x,y,'r');
%一定の温度帯を保つ
T1=5*log(15/14);
T2=3/12;
for i=[0:1:20];
    [t,C]=ode45(@water_model,[i*(T1+T2) i*(T1+T2)+T1],45);
    plot(t,C,'g');
    plot([i*(T1+T2)+T1 (i+1)*(T1+T2)],[42 45],'r');
end;
hold off;
saveas(figure(1),'water_model3','png');


ただ、本のモデルだと加熱中の熱の放出を無視しているので、それは具体的にコスト計算してるのにどうなのかと。
加熱中はq:ヒーターの単位時間当たりの加熱量が加わるので、
(2)式が
\frac{dh(t)}{dt}=-jaC(t)+q
と変わるので、(3)式は
\frac{dC}{dt}=[\frac{-ja}{cw*mw+cc*mc}]C+[\frac{q}{cw*mw+cc*mc}]
となるので、
C(t)=C(0)*e^{\(-k_{1}*t\)}+\frac{k_{2}}{k_{1}} \(k_{2}=[\frac{q}{cw*mw+cc*mc}]\)
ということで、それを組み込んでもう一度matlabで書いてみました。

function dCdt=water_model(t,C);
a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
dCdt=-k1*C;

function dCdt=water_model2(t,C);
a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
q=120000;
k2=q/(cw*mw+cc*mc);
dCdt=-k1*C+k2;


a=100;j=20;
cw=10;mw=900;
cc=5;mc=200;
k1=j*a/(cw*mw+cc*mc);
[t,C]=ode45(@water_model,[0,9],45);
plot(t,C,'b');
title('水の加熱と冷却');
xlabel('時刻');
ylabel('周囲より高い温度');
axis([0,15,0,50]);
%最後に温める
hold on;
q=120000;
k2=q/(cw*mw+cc*mc);
Cm1=45*exp(-k1*9);
[t,C]=ode45(@water_model2,[9,15],Cm1);
plot(t,C,'r');
%一定の温度帯を保つ
T1=5*log(15/14);
T2=5*log(6/5);
for i=[0:1:10];
    [t,C]=ode45(@water_model,[i*(T1+T2) i*(T1+T2)+T1],45);
    plot(t,C,'g');
   [t,C]=ode45(@water_model2,[i*(T1+T2)+T1 (i+1)*(T1+T2)],42);
    plot(t,C,'r');
end;
hold off;
saveas(figure(1),'water_model2','png');

２つの図を比べると、結構違いがあったようです。
２つの具体的な違いの計算と、腎臓のモデルは後日あげます。