mems.html

mems.mws

> restart;

> ode := diff(w(eta),eta,eta)+(n-1)/eta*diff(w(eta),eta)=1/w(eta)^2;

> n := 2;

> L0:=1e-3:

>

> pp := dsolve({ode, w(L0)=1, D(w)(L0)=0}, numeric):

> pp(1);

> with(plots):

> odeplot(pp, [eta^2/w(eta)^3, 1/w(eta)], eta=L0..100, numpoints=1300);pic1:=%:

> odeplot(pp, [eta^2/w(eta)^3, 1/w(eta)], eta=100..4000, numpoints=300, color=blue);pic2:=%:

>

> display(pic1,pic2, title="Bif diagram in 2D", labels=["lambda", "u(0)"], plot([4/9,y,y=0..1], color=black, linestyle=dot));

> n := 1: pp1 := dsolve({ode, w(L0)=1, D(w)(L0)=0}, numeric):

> odeplot(pp1, [eta^2/w(eta)^3, 1/w(eta)], eta=L0..100, numpoints=1000, title="Bif diagram in 1D", labels=["lambda", "u(0)"]);

> pic:=NULL:
b:=0.1:
to 20 do
odeplot(pp, [eta/b, eval(1/w(eta), pp(b))*w(eta)], eta=L0..b, numpoints=300, title=sprintf("lambda=%a", b^2*eval(1/w(eta)^3, pp(b))));
pic:=pic,%:
b:=b*2:
end:

> display(pic, insequence=true);

> plot(eta^(2/3), eta=0..1, color=black, thickness=1);

> display(%,pic): display(%, title="");

>