TMA MS325 02 Question 2, part (d) restart: f:=x->x*exp(r*(1-x)); r:=2.50; plot({f(x),x},x=0..2); # There is a fixed point at x=1. So, the points of the period-2 orbit will lie on either side of x=1 # So we locate each point by solving the composition of f with itself equal to x x1:=fsolve((f@@2)(x)=x,x=0..1); x2:=fsolve((f@@2)(x)=x,x=1..2); # Now we can classify the orbit, by evaluating the derivative of f at each point and taking the product D(f)(x1)*D(f)(x2); # The absolute value is less than 1, hence it is a period-2 sink # Now for r=2.55, proceeding in the same manner... r:=2.55: x1:=fsolve((f@@2)(x)=x,x=0..1); x2:=fsolve((f@@2)(x)=x,x=1..2); D(f)(x1)*D(f)(x2); # The absolute value is greater than 1, so it is a period-2 source TMA MS325 02 Question 2, part (e) bifdiagram := proc(r1,r2,N::posint) local bifdiag,r,dr,x,pts,j,cond1,cond2: cond1 := style=point,symbol=point,view=[r1..r2,0..1],color=black: cond2 := labels=["r","x*"],labelfont=[TIMES,ITALIC,18],axes=BOXED: dr := evalf((r2-r1)/N); bifdiag := NULL: for r from r1 to r2 by dr do x := 0.7; for j from 1 to 100 do x := evalf(x*exp(r*(1-x))); # Removes transients end do; pts := NULL: for j from 1 to 200 do x := evalf(x*exp(r*(1-x))); # Finds attracting point pts := pts,[r,x]; # Appends attracting point # to attractor end do; bifdiag := bifdiag,plot([pts],cond1,cond2); # Plots attractor at given a and appends plot # of attractor to bifurcation diagram end do: plots:-display(bifdiag); end proc: bifdiagram(0,4,200); TMA MS325 02 Question 2, part (g) lypt_Ricker := proc(r1,r2,xint,M::posint,N::posint) local x,r,dr,h,j,pts,cond1,cond2,P1,P2: cond1:= color=black,title="Lyapunov Exponent h vs r for the Ricker Model": cond2:= labels=["r","h"],labelfont=[TIMES,ITALIC,18],axes=BOXED: dr := evalf((r2-r1)/M): pts:=NULL: for r from r1 to r2 by dr do h:=0: x := xint: for j from 1 to 100 do # Removes transients x := evalf(x*exp(r*(1-x))): end do: for j from 1 to N do x := evalf(x*exp(r*(1-x))): h := evalf(h+ln(abs(exp(r*(1-x))*(1-x*r)))): end do: h:=h/N; pts := pts,[r,h]: end do: P1 := plot([pts],cond1): P2 := plot(0,xaxis=r1..r2,color=black): plots:-display({P1,P2},cond2): end proc: lypt_Ricker(0,4,0.4,500,200);