'diff(x,t)=k*x*(1-x);

ode2(%,x,t);

ic1(%,t=0,x=x[0]),logcontract;

solve(%,x);

define(x(t),rhs(first(%o4)))$

wxplot2d(makelist(subst([k=d*0.25,x[0]=0.1],x(t)),d,1,7),
[t,0,15],cons(legend,makelist(k=d*0.25,d,1,7)),
[gnuplot_preamble,"set key right bottom"]);

solve(r*c*(1-c)=c,c);

F(r,x):=r*x*(1-x)$

load(dynamics)$

set_random_state(make_random_state (true))$

x[0]:random(1.0);

evolution(F(0.25,x),x[0],15,[y,0,1]);

evolution(F(1.3,x),x[0],15,[y,0,1]);

evolution(F(3.5,x),0.3,25,[y,0,1]);

staircase(F(3.5,x),x[0],500,[x,0,1],[y,0,1]);

orbits(F(r,x), x[0], 150, 200, [r, 2.5, 4], [style, dots]);

orbits(F(r,x),x[0],150,200,[r,3.5,3.6],[x,0.3,0.4],
[style,dots]);

evolution(F(4,x),3/4,50,[style,[lines,2]]);

eps:random(0.000000000001);

evolution(F(4,x),3/4+eps,50,[style,[lines,2]]);

evolution(F(4,x),0.2,50000,[style,dots]);

numer:false$

latractor: [10*y-10*x, -x*z+28*x-y, x*y-8*z/3]$

linitial: [-8, 8, 27]$

lsolution:rk(latractor,[x,y,z],linitial,[t,0,50,0.01])$

lpoints: map(lambda([x], rest(x)), lsolution)$

load(draw)$

draw3d(point_type=none,points_joined=true,color=orange,
xlabel="x(t)",ylabel="y(t)",zlabel="z(t)",
xtics=10,ytics=10,ztics=10,points(lpoints));

remfunction(x)$

gradef(x(t),%sigma*(y(t)-x(t)))$

gradef(y(t),x(t)*(%rho -z(t))-y(t))$

gradef(z(t),x(t)*y(t)-%beta *z(t))$

x(t)^2+y(t)^2+(z(t)-(%sigma +%rho))^2=R^2*(%sigma +%rho)^2;

facsum(diff(%,t),x(t),y(t),z(t));

factor(solve((sqrt(2*%beta)*z(t)-u)^2-u^2=
2*%beta*z(t)^2-2*%beta*(%sigma+%rho)*z(t),u));

solve(1/a = sqrt(2*%sigma/(%beta*(%sigma +%rho)^2/2)),a);

solve(1/b = sqrt(2/(%beta*(%sigma +%rho)^2/2)),b);

solve(1/c = sqrt(2*%beta/(%beta*(%sigma +%rho)^2/2)),c);

assume(%epsilon>0,N>0)$

remvalue(x[0])$

solve(abs(f[N](x[0]+%epsilon)-f[N](x[0]))=
%epsilon*exp(N*%lambda),%lambda);

gradef(f(x),f\'(x))$

atvalue(f(y),y=x[0],x[1])$

atvalue('diff(f(y),y),y=x[0],f\'(x[0]))$

at(diff(f(f(x)),x),x=x[0]);

define(dF(r,x),diff(F(r,x),x))$

r:3$

maxiter:50000$

x[1]:random(1.0);

for j:2 thru maxiter do (x[j]:F(r,x[j-1]))$

L:0$

for j:1 thru maxiter do L:L+log(abs(dF(r,x[j])))$

'%lambda=L/maxiter;

'V(x)=integrate(x^3/4-x,x);

wxplot2d(integrate(x^3/4-x,x),[x,-3,3],[ylabel,"V(x)"]);

duff:[v,-v/10+x-x^3/4]$

icduff:[3,10]$

solduff:rk(duff,[x,v],icduff,[t,0,100,0.1])$

curveduff:map(lambda([x],rest(x,-1)),solduff)$

pointsduff:map(lambda([x],rest(x)),solduff)$

wxdraw2d(point_type=none,points_joined=true,
xlabel="t",ylabel="x(t)",points(curveduff));

wxdraw2d(point_type=none,points_joined=true,
xlabel="x(t)",ylabel="v(t)",points(pointsduff));

duffing:[v,-v/10+x-x^3/4+sin(t)]$

icduffing:[0,0]$

sduffing:rk(duffing,[x,v],icduffing,[t,0,100,0.1])$

cduffing:map(lambda([x],rest(x,-1)),sduffing)$

pduffing:map(lambda([x],rest(x)),sduffing)$

wxdraw2d(point_type=none,points_joined=true,color=magenta,
xlabel="t",ylabel="x(t)",points(cduffing));

wxdraw2d(point_type=none,points_joined=true,color=magenta,
xlabel="x(t)",ylabel="v(t)",points(pduffing));

duffing1:[v,-v/10+x-x^3/4+2.5*sin(2*t)]$

icduffing1:[0,0]$

sduffing1:rk(duffing1,[x,v],icduffing1,[t,0,100,0.1])$

pduffing1:map(lambda([x],rest(x)),sduffing1)$

wxdraw2d(point_type=none,points_joined=true,color=coral,
xlabel="x(t)",ylabel="v(t)",points(pduffing1));

miter:25$

%tau:bfloat(%pi)$

sduffing2:rk(duffing,[x,v],icduffing,
[t,0,miter*2*%tau,%tau/30])$

pduffing2:create_list(sduffing2[i],i,makelist(i*60,i,1,miter))$

poinduffing2:map(lambda([x],rest(x)),
makelist(pduffing2[i],i,1,miter))$

wxdraw2d(point_type=filled_circle,color=magenta,xtics=1,
ytics=1,xrange=[-4,1.5],yrange=[0,4],points(poinduffing2));

transient:map(lambda([x],rest(x)),
makelist(pduffing2[i],i,5,miter))$

wxdraw2d(point_type=filled_circle,color=magenta,xtics=1,
ytics=1,xrange=[-4,1.5],yrange=[0,4],points(transient));

duffing3:[v,-v/10+x-x^3/4+2.5*sin(2*t)]$

icduffing3:[0,0]$

maxiter:1000$

sduffing3:rk(duffing3,[x,v],icduffing3,
[t,0,maxiter*%tau,%tau/30])$

pduffing3:
create_list(sduffing3[i],i,makelist(i*30,i,1,maxiter))$

poinduffing3:map(lambda([x],rest(x)),
makelist(pduffing3[i],i,1,maxiter))$

wxdraw2d(point_size=0.3,point_type=circle,
xrange=[-5,5],yrange=[-7,3],xtics=1,ytics=1,
color=coral,points(poinduffing3),grid=true);

resolution:20$

for n:1 thru resolution do
(Z[n]:substpart("[",zeromatrix(10*n,10*n),0),
boxcount[n]:0,
for k:1 thru maxiter do (
ix:floor(n*(poinduffing3[k][1]+5)),
iy:floor(n*(poinduffing3[k][2]+7)),
if is(Z[n][ix][iy]=0) then
(Z[n][ix][iy]:1,boxcount[n]:boxcount[n]+1)))$

fitdim:makelist(
[log(n)/log(10),log(boxcount[n])],n,1,resolution),numer$

wxdraw2d(point_size=1,point_type=filled_circle,
color=dark_violet,points(fitdim));

fitdimension:rest(rest(fitdim,7),-2)$

load(stats)$

model:linear_regression(fitdimension);

fiteqn:take_inference('b_estimation,model);

fractaldimension:second(%);