数理モデリング 2021　高須

2021/12/14, 21,

Mathematicaの練習

式もしくはコマンドを入力して、シフトキー＋Enterキー、で実行する。

1+3

3/4+1/2

myList={1,2,3,4,5}

Length[myList]

myList[[1]]

myList[[6]]

Last[myList]

AppendTo[myList,9]

Last[myList]

?Last

?Do

Do[
Print[i],{i,1,10}
]

tmp={}
Do[
AppendTo[tmp,i],{i,5,10}
]
tmp

関数の定義

fn[x_]:=x x-2

fn[3]

?Plot

Plot[fn[x],{x,-2,2},PlotRange->All]

?D

D[fn[x],x]

?Integrate

Integrate[fn[x],{x,0,5}]

mat={{1,2},{3,4}}

MatrixForm[mat]

Eigenvalues[mat]

Transpose[mat]

MatrixForm[Transpose[mat]]

非負の実数としての集団サイズに注目したモデル（差分式、微分方程式）

差分式の数値解析

離散時間の個体群動態モデル: N_{t+1} = f( N_t ) の形の差分式

指数増加・減少モデル

fn=.

fn[x_]:=r x/.r->1.1

fn[10]

fn[2]

Plot[fn[x],{x,0,2}]

seq=NestList[fn,1,50]

seq

ListPlot[seq]

ListPlot[seq,Joined->True,PlotMarkers->Automatic]

ListLogPlot[seq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","N_t"}]

Plot[{fn[x],x},{x,0,4},PlotRange->{{0,4},{0,4}},AspectRatio->1,AxesLabel->{"N_t","N_t+1"}]

ロジスティックモデル（離散時間）

fn[x_]:=r (1-x)x/.r->3.2

fn[0.2]

Plot[{fn[x],x},{x,0,1},AspectRatio->1]

seq=NestList[fn,0.01,50]

ListPlot[seq,Joined->True,PlotMarkers->Automatic,PlotRange->{0,1},AxesLabel->{"t","N_t"}]

ListLogPlot[seq,Joined->True,PlotMarkers->Automatic,PlotRange->All]

gCobweb=Plot[{fn[x],x},{x,0,1},AspectRatio->1,AxesLabel->{"N_t","N_t+1"}]

convertData[x_List]:=Module[
{len=Length[x],seqCob={}},
AppendTo[seqCob,{x[[1]],0}];
Do[
AppendTo[seqCob,{x[[i]],x[[i+1]]}];
AppendTo[seqCob,{x[[i+1]],x[[i+1]]}],{i,1,len-1}
];
Return[seqCob];
]

seq

seq2=convertData[seq]

gCobwebSim=ListPlot[seq2,Joined->True,PlotStyle->Red]

Show[gCobweb,gCobwebSim]

Rickerロジスティックモデル（離散時間）

fn[x_]:=Exp[r (1-x)]x/.r->3

fn[0.2]

seq=NestList[fn,0.01,200]

ListPlot[seq,Joined->True,PlotMarkers->Automatic]

ListLogPlot[seq,Joined->True,PlotMarkers->Automatic]

gPhase=Plot[{fn[x],x},{x,0,2},AspectRatio->1]

seq

seq2=convertData[seq];

gCob=ListPlot[seq2,Joined->True,PlotStyle->Red]

Show[gPhase,gCob]

gPhase2=Plot[{fn[fn[x]],x},{x,0,2},AspectRatio->1]

Show[gPhase2,gCob]

Nicholson - Baileyモデル （2 種系寄生モデル）

nextState[{h_,p_}]:=Module[{hh,pp,next},
hh=r Exp[-a p]h;
pp=(1-Exp[-a p])h;
next={hh,pp}/.{r->1.2,a->1};
Return[next];
]

nextState[{1,0.2}]

seq=NestList[nextState,{1,0.2},100]

ListPlot[seq,Joined->True,PlotMarkers->Automatic]

ListPlot[seq,Joined->True,PlotMarkers->Automatic,PlotRange->All,AxesLabel->{"H","P"}]

MatrixForm[seq]

seqT=Transpose[seq]

MatrixForm[seqT]

ListPlot[seqT,Joined->True,PlotMarkers->Automatic]

MatrixForm[seq]

MatrixForm[Transpose[seq]]

{seqh,seqp}=Transpose[seq];

ListPlot[{seqh,seqp},Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","H, P"}]

Nicholson - Baileyモデル + Host Ricker model （2 種系寄生モデル）

nextState[{h_,p_}]:=Module[{hh,pp,next},
hh=Exp[r(1-h/k)] Exp[-a p]h;
pp=(1-Exp[-a p])h;
next={hh,pp}/.{r->1.2,a->1,k->5};
Return[next];
]

nextState[{1,0.2}]

seq=NestList[nextState,{1,0.2},50];

ListPlot[seq,Joined->True]

ListPlot[seq,Joined->True,PlotMarkers->Automatic,PlotRange->All,AxesLabel->{"H","P"}]

seqT=Transpose[seq];

ListPlot[seqT,Joined->True,PlotMarkers->Automatic]

{seqh,seqp}=Transpose[seq];

ListPlot[{seqh,seqp},Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","H, P"}]

常微分方程式の数値解

指数モデル

deq={x'[t]==r x[t],x[0]==x0}

DSolve[deq,x[t],t]

?DSolve

?Plot

deqN=deq/.{r->1,x0->0.01}

tEnd=10

sol=NDSolve[deqN,x[t],{t,0,tEnd}]

Plot[Evaluate[x[t]/.sol],{t,0,tEnd},PlotRange->All]

LogPlot[Evaluate[x[t]/.sol],{t,0,tEnd}]

ロジスティックモデル （連続時間）

deq={x'[t]==r(1-x[t])x[t],x[0]==x0}

DSolve[deq,x[t],t]

deqN=deq/.{r->1,x0->0.01}

tEnd=10

sol=NDSolve[deqN,x[t],{t,0,tEnd}]

Plot[Evaluate[x[t]/.sol],{t,0,tEnd}]

LogPlot[Evaluate[x[t]/.sol],{t,0,tEnd}]

Lotka Volterra競争モデル

deq = {n1'[t]==r1(1-(n1[t] + a12 n2[t])/k1)n1[t],
n2'[t]==r2(1-(a21 n1[t] + n2[t])/k2)n2[t],
n1[0]==n10,n2[0]==n20}

DSolve[deq,{n1[t],n2[t]},t]

para = {r1->1, k1->1, a12->0.2, r2->1, k2->1, a21->1.3}

deqN1 =deq/.para/.{n10->0.1, n20->0.1}

NDSolveで微分方程式を数値的に解く

tEnd=50;
sol1 = NDSolve[deqN1,{n1[t],n2[t]},{t,0,tEnd}]

Plot[ Evaluate[{n1[t],n2[t]}/.sol1],{t,0,tEnd},AxesLabel->{"t","n1, n2"},PlotRange->All]

g1 = ParametricPlot[ Evaluate[{n1[t],n2[t]}/.sol1],{t,0,50},PlotRange->All,AxesLabel->{"n1","n2"},AxesOrigin->{0,0}]

deqN2 =deq/.para/.{n10->1, n20->0.1}

sol2 = NDSolve[deqN2,{n1[t],n2[t]},{t,0,tEnd}]

g2 = ParametricPlot[ Evaluate[{n1[t],n2[t]}/.sol2],{t,0,50},PlotRange->All,AxesLabel->{"n1","n2"},AxesOrigin->{0,0}]

deqN3 =deq/.para/.{n10->0.1, n20->1}

sol3 = NDSolve[deqN3,{n1[t],n2[t]},{t,0,tEnd}]

g3 = ParametricPlot[ Evaluate[{n1[t],n2[t]}/.sol3],{t,0,50},PlotRange->All,AxesLabel->{"n1","n2"},AxesOrigin->{0,0}]

deqN4 =deq/.para/.{n10->1, n20->1}

sol4 = NDSolve[deqN4,{n1[t],n2[t]},{t,0,tEnd}]

g4 = ParametricPlot[ Evaluate[{n1[t],n2[t]}/.sol4],{t,0,50},PlotRange->All,AxesLabel->{"n1","n2"},AxesOrigin->{0,0}]

Show[g1,g2,g3,g4,PlotRange->{{0,1.5},{0,1.5}}]

?VectorPlot

vec={r1(1-n1-a12 n2)/k1,r2(1-a21 n1-n2)/k2}/.para

gvec=VectorPlot[vec,{n1,0,1.5},{n2,0,1.5}]

Show[g1,g2,g3,g4,gvec,PlotRange->{{0,1.5},{0,1.5}}]

gcont=ContourPlot[vec=={0,0},{n1,0,1.5},{n2,0,1.5}]

Show[g1,g2,g3,g4,gvec,gcont,PlotRange->{{0,1.5},{0,1.5}}]

Lotka Volterra捕食モデル

deq = {h'[t]==rH h[t] - a h[t]p[t],
p'[t]==-rP p[t]+b h[t]p[t],
h[0]==h0,p[0]==p0}

para = {rH->1,a->1, rP->0.8, b->1}

deqN1 =deq/.para/.{h0->1, p0->0.5}

NDSolveで微分方程式を数値的に解く

tEnd=50;
sol1 = NDSolve[deqN1,{h[t],p[t]},{t,0,tEnd}]

Plot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},PlotRange->{0,2},AxesLabel->{"t","h[t], p[t]"}]

g1 = ParametricPlot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},AspectRatio->1,PlotRange->{{0,2},{0,2}},AxesLabel->{"h[t]", "p[t]"}]

Lotka Volterra捕食モデル + Logistic growth of the prey

deq = {h'[t]==rH( 1-h[t]/k)h[t] - a h[t]p[t],
p'[t]==-rP p[t]+b h[t]p[t],
h[0]==h0,p[0]==p0}

para = {rH->1,a->1, rP->0.8, b->1,k->8}

deqN1 =deq/.para/.{h0->1, p0->0.5}

tEnd=100;
sol1 = NDSolve[deqN1,{h[t],p[t]},{t,0,tEnd}]

Plot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},PlotRange->{0,2},AxesLabel->{"t","h[t], p[t]"}]

g1 = ParametricPlot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},AspectRatio->1,PlotRange->{{0,2},{0,2}},AxesLabel->{"h[t]", "p[t]"}]

Lotka Volterra捕食モデル+ Logistic growth of the 
prey + Satiation of the predator

deq = {h'[t]==rH( 1-h[t]/k)h[t] - a h[t]p[t]/(1+c h[t]),
p'[t]==-rP p[t]+b h[t]p[t]/(1+c h[t]),
h[0]==h0,p[0]==p0}

para = {rH->1,a->1, rP->0.8, b->1,k->8,c->0.9}

deqN1 =deq/.para/.{h0->4, p0->4}

tEnd=200;
sol1 = NDSolve[deqN1,{h[t],p[t]},{t,0,tEnd}]

Plot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},PlotRange->{0,10},AxesLabel->{"t","h[t], p[t]"}]

gDyna = ParametricPlot[ Evaluate[{h[t],p[t]}/.sol1],{t,0,tEnd},AspectRatio->1,PlotRange->{{0,10},{0,10}},AxesLabel->{"h[t]", "p[t]"}]

ローレンツ方程式 Lorenz

deq={x'[t]==-p x[t]+p y[t],
y'[t]==-x[t]z[t]+r x[t]-y[t],
z'[t]==x[t]y[t]-b z[t],
x[0]==x0,y[0]==y0,z[0]==z0}

para={p->10,r->28,b->8/3}

deqN=deq/.para/.{x0->0.1,y0->0.1,z0->-1}

tEnd=100;
sol=NDSolve[deqN,{x[t],y[t],z[t]},{t,0,tEnd}]

Plot[Evaluate[{x[t],y[t],z[t]}/.sol],{t,0,tEnd},AxesLabel->{"t","x, y, z"}]

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol],{t,0,tEnd},BoxRatios->{1,1,1},AxesLabel->{"x","y","z"},PlotPoints->1000]

Van der Pol oscillator

Rössler system

空間分布を考慮したモデル

生物個体を粒子と見なして、空間上 （2 次元空間） を移動する様子をモデル化する。

ランダム拡散モデル（ランダムウォーク）

各粒子が独立に一定のルールで移動を繰り返すモデルを考える

具体例0：単位時間毎に同期して上下左右に一定の距離 Lだけ移動する。離散時間のモデル。
具体例1：単位時間毎に等方的に同期して一定の距離 Lだけ移動する。離散時間のモデル。
具体例2：単位時間毎に等方的に同期して距離Lだけ移動するが、Lはある確率変数。離散時間のモデル。
具体例3：各粒子が等方的に非同期で距離Lだけ移動するが、Lはある確率変数。連続時間のモデル。

ランダム拡散をシミュレーションし、空間上の点の分布を可視化する。
簡単なプログラミングで可視化を可能とするフリーのライブラリ openFrameworks を用いる。

URL: http://openframeworks.cc

Mathematicaによるランダムウォークシミュレーション　具体例0

points=Table[{0,0},{100}]

ListPlot[points,AspectRatio->1]

moveAPoint[{x_,y_}]:=Module[{distance=0.01,rv,xt,yt,dx,dy},
rv=Random[];
If[rv<0.25,dx=distance;dy=0,
If[rv<0.5,dx=-distance;dy=0,
If[rv<0.75,dx=0;dy=-distance,
dx=0;dy=distance
]
]
];
xt=x+dx;
yt=y+dy;
Return[{xt,yt}];
]

points[[1]]

moveAPoint[points[[1]]]

moveAPoint[moveAPoint[points[[1]]]]

NestList[moveAPoint,points[[1]],10]

ListPlot[NestList[moveAPoint,points[[1]],10]]

movePoints[pts_]:=Module[{len=Length[pts],pt,ptsAfter={}},
Do[
pt=moveAPoint[pts[[i]]];
AppendTo[ptsAfter,pt],{i,1,len}
];
Return[ptsAfter];
]

movePoints[points]

result=NestList[movePoints,points,100];

ListPlot[result[[5]],AspectRatio->1,PlotRange->{{-1,1},{-1,1}}]

glist={};
Do[
g=ListPlot[result[[i]],AspectRatio->1,PlotRange->{{-1,1},{-1,1}}];
AppendTo[glist,g],{i,1,Length[result]}
]

glist

ListAnimate[glist]

Mathematicaによるランダムウォークシミュレーション　具体例1

points=Table[{0,0},{100}];

ListPlot[points,AspectRatio->1]

moveAPoint[{x_,y_}]:=Module[{distance=0.01,angle,xt,yt,dx,dy},
angle=Random[]*2*Pi;
dx=distance*Cos[angle];
dy=distance*Sin[angle];
xt=x+dx;
yt=y+dy;
Return[{xt,yt}];
]

points[[1]]

moveAPoint[points[[1]]]

movePoints[pts_]:=Module[{len=Length[pts],pt,ptsAfter={}},
Do[
pt=moveAPoint[pts[[i]]];AppendTo[ptsAfter,pt],{i,1,len}
];
Return[ptsAfter];
]

movePoints[points];

result=NestList[movePoints,points,50];

ListPlot[result[[2]],AspectRatio->1,PlotRange->{{-0.1,0.1},{-0.1,0.1}}]

ListPlot[result[[3]],AspectRatio->1,PlotRange->{{-0.1,0.1},{-0.1,0.1}}]

glist={};
Do[
g=ListPlot[result[[i]],AspectRatio->1,PlotRange->{{-1,1},{-1,1}}];
AppendTo[glist,g],{i,1,Length[result]}
]

ListAnimate[glist]