Modeling Colchea
A System with localized Resonnances

E.A. Bartnik
Institute of Theoretical Physics
University of Warsaw
eabartnik@gmail.com

Mathematica Day
     Warsaw, 2017

Cochlea models

ClearAll["Global`*"]
$PrePrint=MatrixForm

Ear anatomy

image_ear=

image3=

Below (Subscript[x, n]) is the displacement of the n-th mass, (Subscript[v, n]) - it's velocity.
(Subscript[M, n]) represents the fluid displaced by the membrane composed od masses (Subscript[x, n]).
(Subscript[k, n]) are spring parameters (stiffness)

Parameter setup
Mass and Stiffness matrices definition

nmax=64
h=1.0/nmax
xx=Table[x[n],{n,nmax+1}];
m[n_]=0.00035;
mucont=m[1] nmax;
Print["mu cont = ",mucont]

mm1=1.0;
mm1cont=mm1/nmax;
Print["mm1 cont = ",mm1cont]
mm2=0.1;
mm2cont=mm2/nmax;
Print["mm2 cont = ",mm2cont]
bm=(mm2-mm1)/(Sqrt[nmax]-1.);
am=mm1-bm;
mm[n_]=am+Sqrt[n] bm;

k[n_]=1.0;
kcont=k[1]/nmax;
Print["k cont = ",kcont]
x[nmax+1]=0;
ekin1=Sum[ m[n] x[n]^2/2,{n,nmax}];
mass1=Table[D[ekin1,x[n],x[m] ],{n,nmax},{m,nmax}];
ekin2=Sum[ mm[n](x[n+1]-x[n])^2/2,{n,nmax}];
mass2=Table[D[ekin2,x[n],x[m] ],{n,nmax},{m,nmax}];
pot=Sum[ k[n](x[n+1]-x[n])^2/2,{n,nmax}];
stiff=Table[D[pot,x[n],x[m] ],{n,nmax},{m,nmax}];

Finding eigenfrequencies and eigenmodes

{vals,vecs}=Eigensystem[{stiff,mass1+mass2}];
vals;
ListPlot[vals,PlotLabel->"Eigenfrenquencies"]
ListPlot[Log[vals],PlotLabel->"Eigenfrequencies
(log scale)"]
ListPlot[vecs^2,Joined->True ,PlotRange->All,
PlotLabel->"Eigenfunctions"]

ListPlot[vecs[[nmax/2]],Joined->True ,PlotRange->All,
PlotLabel->"Typical eigenfunction"]

Continous case:
New class of wave equations

Lagrange density

ℒ = (m/2) ((( ((∂ϕ)/(∂t))))^2) + ((M(x))/2)((((((∂^2)ϕ)/(∂x∂t)) ))^2) - ((K(x))/2)((( ((∂ϕ)/(∂x))))^2)

Euler-Lagrange Equation

m ((((∂^2)ϕ)/(∂(t^2))) ) = (∂/(∂x))(M(x) (((∂^3)ϕ)/(∂x∂(t^2))))- (∂/(∂x))(K(x) ((∂ϕ)/(∂x)))

Stationary solution

ϕ(t,x)=(e^(i ω t)) u(x)

Equation for u(x)

((ω^2)M(x)-K(x))  (((d^2)u(x))/(dx^2)) = (-(ω^2) ((dM(x))/dx)+((dK(x))/dx) ) ((du(x))/dx)+m (ω^2)u(x)

ClearAll["Global`*"]
$PrePrint=MatrixForm

Adjusting parameters to get logarytmic frequency response

xmin=0.0;
xmax=1.0;

wmax=10.0; (* for xinr *)
wmin=1.0;   (* for xmax *)

mu=1.0;
mzer=1.0;

alf=Log[wmin/wmax]/xmax
wx=wmax Exp[alf x]
Plot[{wx,0},{x,xmin,xmax},
PlotLabel->"w(x)"]
Plot[Log[wx],{x,xmin,xmax},
PlotLabel->"Log(w(x))"]

k[x_]:=mzer wmax^2 Exp[2 alf x];





Plot[k[x],{x,xmin,xmax},
PlotRange->All,
PlotLabel->"k(x)"];

Plot[{k[x],0},{x,0.9 xmax,xmax},
PlotRange->All,
PlotLabel->"Log(k(x)}"];

Solving equations for a mode with frequency w

A small imaginary part is added to frequency w to awoid the singularity

eps=0.01;
w=2.5+ I eps;
a=k[x]/w^2-mzer-0.0 I eps;
Plot[{Re[a],Im[a]},{x,xmin,xmax}];
b=D[k[x],x]/w^2;
c=mu;

eqs={a u''[x]+b u'[x]+c u[x]==0,
u[xmax]==0,
u'[xmax]==1}

sol=NDSolve[eqs,u,{x,xmax,xmin}]
Plot[Re[u[x] /. sol],{x,xmin,xmax},
PlotRange->All];
Plot[Im[u[x] /. sol],{x,xmin,xmax},
PlotRange->All];
Plot[Abs[u[x] /. sol],{x,xmin,xmax},
PlotRange->All,
PlotLabel->"TYPICAL CONTINOUS EIGENFUNCTION"]

Facit

It is necessary to use a system with non diagonal mass matrix.
The new wave equation is sharply peaked (has a singularity) when the coefficient
of the second space derivative vanishes.