FitzHugh-Nagumo model (Ch 14, page 235)

Clear["Global`*"]
Off[General::spell1]

<<VectorFieldPlots`;

System

Parameters

{(Subscript[x, min]),(Subscript[x, max])}={-0.8,1.2};
{(Subscript[y, min]),(Subscript[y, max])}={-0.2,0.4};

f[x_,y_,a_,b_,ξ_]:=(1/ξ)(x(1-x)(x+a)-y+b)

g[x_,y_,ϵ_]:=x-ϵ y

Vector field

V={f[x,y,a,b,ξ],g[x,y,ϵ]}

Fixed Points & Nullclines

(Subscript[x, 1])=.
(Subscript[y, 1])=.

fp=Solve[{f[(Subscript[x, 1]),(Subscript[y, 1]),0.1,0.01,0.01]==0,g[(Subscript[x, 1]),(Subscript[y, 1]),0.5]==0},{(Subscript[x, 1]),(Subscript[y, 1])}]

fp1=Sort[fp[[1]]];

{(Subscript[x, 1]),(Subscript[y, 1])}={Extract[fp1,{1}][[2]],Extract[fp1,{2}][[2]]};
fixedpoint1=Show[Graphics[{PointSize[0.025],Hue[0.33],Point[{(Subscript[x, 1]),(Subscript[y, 1])}]},PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}}]];

nullclines=ContourPlot[{f[x,y,0.1,0.01,0.01]==0},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},ContourStyle->{{Green,Thickness[0.005]}}];

Show[nullclines,fixedpoint1]

Jacobian

a11=Simplify[D[f[x,y,a,b,ξ],x]]

a12=Simplify[D[f[x,y,a,b,ξ],y]]

a21=Simplify[D[g[x,y,ϵ],x]]

a22=Simplify[D[g[x,y,ϵ],y]]

J={{a11,a12},{a21,a22}}

MatrixForm[J]

Acceleration

γ=Factor[J . V]

Fixed Points Stability - Eigenvalues

Eigenvalues[J]/.{x-> (Subscript[x, 1]),y->(Subscript[y, 1])}

Phase Portrait

solution=NDSolve[{x'[t]==f[x[t],y[t],0.1,0.01,0.01],x[0]==1,

y'[t]==g[x[t],y[t],0.5],y[0]==1},{x,y},

{t,0,100},Method->"ExplicitRungeKutta", MaxSteps->Infinity,StartingStepSize->0.001];phaseplot=ParametricPlot[Evaluate[{x[t],y[t]}/.solution],{t,50,100},PlotPoints->1000,PlotRange->{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},AxesLabel->{"X","Y"},PlotStyle->{{Red,Thickness[0.005]}},AspectRatio->1,Frame->True]

Show[phaseplot,fixedpoint1,nullclines]

Flow Curvature Manifold

ϕ=FullSimplify[Det[{V,γ}]]

ϕL=Factor[Det[{V,γ}]/.{ξ-> 0.01,ϵ-> 0.5,a-> 0.1,b-> 0.01}]

flowmanifoldL=ContourPlot[{ϕL==0},{x,(Subscript[x, min]),0},{y,(Subscript[y, min]),(Subscript[y, max])},ContourStyle->{{Blue,Thickness[0.005]}},AspectRatio->1,Frame->True]

ϕR=Factor[Det[{V,γ}]/.{ξ-> 0.01,ϵ-> 0.5,a-> 0.1,b-> 0.01}]

flowmanifoldR=ContourPlot[{ϕR==0},{x,0.6,(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},ContourStyle->{{Blue,Thickness[0.005]}},AspectRatio->1,Frame->True]

Flow Curvature Invariance

Normal Vector

η=Factor[{D[ϕ,x],D[ϕ,y]}]

Invariance

Ψ=Factor[V . η]

K=Factor[Tr[J]ϕ -((2 (-1+a+3 x) (-x+y ϵ))/(ξ^3))(((b+(1-x) x (a+x)-y))^2)]

FullSimplify[Ψ - K]

Ψ=Factor[V . η]/.{ξ-> 0.01,ϵ-> 0.5,a-> 0.1,b-> 0.01}

invariantpart=ContourPlot[{Ψ==0},{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])},PlotPoints->50,ContourStyle->{{Magenta,Thickness[0.005],Dashing[{.025}]}},Frame->True,AspectRatio->1]

Show[phaseplot,fixedpoint1,nullclines,flowmanifoldL,flowmanifoldR,invariantpart,ImageSize->{600,600}]