1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
| muvect:=[2.5, 3.0, 3.1, 3.56, 3.57, 4.0];
mu:=muvect[1];
f:=x*exp(mu*(1-x));
dfdx:=diff(f,x);
itmax:=10000: # number of iterations to calculate
data:=array(0..itmax): # save the iterates in data
x0:=0.1: # initial condition
for i from 1 by 1 to 1000 do # iterate a few times to get past transients
x0:=eval(f,x=x0):
od:
x0; # see where we are after transients
data[0]:=x0: # compute the orbit
for i from 1 by 1 to itmax do
x0:=eval(f,x=x0):
data[i]:=x0:
od:
Lyapunov:=0: # compute the Lyapunov exponent
for i from 1 by 1 to itmax do
Lyapunov:=Lyapunov + ln(abs(eval(dfdx,x=data[i]))):
od:
Lyapunov:=Lyapunov/itmax; |
Partager