> #OK to post #Timothy Nasralla, HW9, 10/03/2021 > read "/Users/tan88/OneDrive - Rutgers University/M9.txt" > Help9() Orb(f,x,x0,K1,K2), Orb2D(f,x,x0,K) , FP(f,x) , SFP(f,x) , Comp(f,x) > #Question 1: Find the stable points using ORB and SPF (i) 2x(1-x) , (ii) > 2.5*x*(1-x) , (iii) (3, 1)x(1-x) , (iv) (4+x)/(3+x), (v) (3+x)/(4+x), (vi) > (3+x+x^2)/(4+x+2x2) > SFP((2*x*(1-x)),x) ; Orb((2*x*(1-x)),x,.5,1000,1020) [0.5000000000] [0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000, 0.5000000000] > > #The stable fixed points for i are 0.50 > SFP((2.5*x*(1-x)),x) ; Orb((2.5*x*(1-x)),x,.5,1000,1020) [0.6000000000] [0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000, 0.6000000000] > > #The stable fixed points for ii are 0.60 > SFP((x*(1-x)),x) ; SFP((3*x*(1-x)),x) [] [] > Orb((x*(1-x)),x,.5,1000,1020) ; Orb((3*x*(1-x)),x,.5,1000,1020) [0.0009913908610, 0.0009904080051, 0.0009894270971, 0.0009884481311, 0.0009874711014, 0.0009864960022, 0.0009855228278, 0.0009845515726, 0.0009835822308, 0.0009826147968, 0.0009816492650, 0.0009806856297, 0.0009797238854, 0.0009787640265, 0.0009778060475, 0.0009768499429, 0.0009758957071, 0.0009749433347, 0.0009739928202, 0.0009730441582, 0.0009720973432, 0.0009711523700] [0.6591746334, 0.6739903083, 0.6591821178, 0.6739831602, 0.6591895800, 0.6739760328, 0.6591970200, 0.6739689264, 0.6592044378, 0.6739618410, 0.6592118337, 0.6739547760, 0.6592192077, 0.6739477317, 0.6592265598, 0.6739407081, 0.6592338903, 0.6739337046, 0.6592411992, 0.6739267215, 0.6592484865, 0.6739197585] > #There are seemingly no stable fixed points for iii. > SFP(((4+x)/(3+x)),x) ; Orb(((4+x)/(3+x)),x,0.5,1000,1020) [1.236067977] [1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978, 1.236067977, 1.236067978] > #For iv, the stable fixed point is 1.236067977. > SFP(((3+x)/(4+x)),x) ; Orb(((3+x)/(4+x)),x,0.5,1000,1020) [0.791287848] [0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475, 0.7912878475] > #The stable fixed point for v is 0.791287848. > SFP(((3+x+x*x)/(4+x+2*x*x)),x) ; Orb(((3+x+x*x)/(4+x+2*x*x)),x,0.5,1000,1020) [0.7351392587] [0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591, 0.7351392591] > #The stable fixed points for vi are 0.7351392591. > #Question 2: For the system x --> (x+a)/(x+b), experiment with a = 1, b = 2, a > = 2, b = 3, and a = 12, b = 17 to find any SFP and FP's then compare to ORB, > SFP, and FP function. > f := ((x+a)/(x+b)) ; solve((f=x),x) (1/2) (1/2) 1 1 1 / 2 \ 1 1 1 / 2 \ - - b + - + - \b + 4 a - 2 b + 1/ , - - b + - - - \b + 4 a - 2 b + 1/ 2 2 2 2 2 2 > diff(((x+a)/(x+b)),x) 1 x + a ----- - -------- x + b 2 (x + b) > #Generally speaking, the solution from the solve function shows all x values > that are fixed values and the solution from the diff function shows the > function that will be used to determine whether or not any of the fixed points > are stable. > #Case 1, a = 1, b= 2 > f := (x+1)/(x+2) ; solve((f=x),x) 1 1 (1/2) 1 (1/2) 1 - - - - 5 , - 5 - - 2 2 2 2 > test1 := evalf(-1/2-sqrt(5)/2) ; test2 := evalf(-1/2 + sqrt(5)/2) > evalf(1/(test1+2)-(test1+1)/((test1+2)^2)) ; > evalf(1/(test2+2)-(test2+1)/((test2+2)^2)) 6.854101940 0.1458980339 > Orb(((x+1)/(x+2)),x,.5,1000,1020) [0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888, 0.6180339888] > FP((x+1)/(x+2),x) [-1.618033988, 0.6180339880] > SFP((x+1)/(x+2),x) [0.6180339880] > #For the situation where a = 1 and b = 2, the only SFP is 0.618. > #Case 2, a = 2 and b = 3 > f := (x+2)/(x+3) ; solve((f=x),x) (1/2) (1/2) -1 - 3 , 3 - 1 > test1 := evalf(-1-sqrt(3)) ; test2 := evalf(-1+ sqrt(3)) > evalf(1/(test1+3)-(test1+2)/((test1+3)^2)) ; > evalf(1/(test2+3)-(test2+2)/((test2+3)^2)) 13.92820327 0.0717967697 > Orb(((x+2)/(x+3)),x,.5,1000,1020) [0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076, 0.7320508076] > FP((x+2)/(x+3),x) [-2.732050808, 0.732050808] > SFP((x+2)/(x+3),x) [0.732050808] > #For set 2 where a = 2 and b=3, the only SFP is 0.732. > #Case 3, a = 12 and b= 17 > f := (x+12)/(x+17) ; solve((f=x),x) (1/2) (1/2) -8 - 2 19 , -8 + 2 19 > test1 := evalf(-8-2*sqrt(19)) ; test2 := evalf(-8+ 2*sqrt(19)) > evalf(1/(test1+17)-(test1+12)/((test1+17)^2)) ; > evalf(1/(test2+17)-(test2+12)/((test2+17)^2)) 62.78407369 0.01592760650 > Orb(((x+12)/(x+17)),x,.5,1000,1020) [0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871, 0.7177978871] > FP((x+12)/(x+17),x) [-16.71779789, 0.717797888] > SFP((x+12)/(x+17),x) [0.717797888] > #As shown above, the only fixed points are x = 0.717797888. > #Question 3: For an arbitrary k between 1&4, x--> k*x*(1-x) prove that x = 0 > is never stable. For what values of k is the oher fixed point stable? And what > is the first bifurcation point? > #For all values where k = a number from 1-4, setting x = 0 is always going to > be a fixed point such that x = k*x*(1-x), in which 0 = 0. When finding if > they're stable, the derivative of the function will be k - 2*k*x. Solving for > x in x = k*x*(1-x) turns to be k*x-x-k*x^2, which when simplified is 0= > x(k-1-kx) and the fixed points will therefore be 0 = x or (k-1)/k #Whenever x > = 0, the derivative f'(0) will be equal to k which for all k = [1,4], the > value will always be greater than or equal to 1, which means the point is no > stable. diff((k*x-k*x^2),x); Error, (in solve) cannot solve for an unknown function with other operations in its arguments -2 k x + k > #For all values k > 1 and 3 > k, the other fixed point is stable. When k = 1, > the only fixed points are x = 0 which has f'(0) = 1, making it a non stable > fixed point. When k = 3, x = 0 or 2/3 for fixed points. Since we know x = 0 > won't work, I tried 2/3 which results in f'(2/3) = -1, which is also a non > stable fixed point and this applies to all values for 3 < k < 4. Therefore, > the only stable fixed points are inbetween 1 < k < 3. solve((x = x*(1-x)),x); > solve((x = 2*x*(1-x)),x); solve((x = 3*x*(1-x)),x); 0, 0 1 0, - 2 2 0, - 3 > 1 - 2*0*1;-2*2*1/2 + 2; -2*3*2/3 + 3; 1 0 -1 > #Based on the information above the bifurcation point is going to be 3. > #Question 4: For an arbitrary k between 1,4 find the equilibrium points of > x--> f(f(x)) > #f(x) = k*x*(1-x), f(f(x)) = k*k*x*(1-x)*(1-k*x*(1-x)) > simplify(k*k*x*(1-x)*(1-k*x*(1-x))) 2 / 2 \ -k x (-1 + x) \k x - k x + 1/ > solve((-k^2*x*(x-1)*(k*x^2-k*x+1)) = x, x) (1/2) (1/2) 1 1 1 / 2 \ 1 1 1 / 2 \ - k + - + - \k - 2 k - 3/ - k + - - - \k - 2 k - 3/ k - 1 2 2 2 2 2 2 0, -----, -------------------------------, ------------------------------- k k k > #When k<3, the last 2 equilibrium points are imaginary so I tested 3 in which > the last 3 equilibrium points are all equivalent, giving only two fixed points > and period 2 orbit. solve((-3^2*x*(x-1)*(3*x^2-3*x+1)) = x, x) 2 2 2 0, -, -, - 3 3 3 > solve((-3.1^2*x*(x-1)*(3.1*x^2-3.1*x+1)) = x,x); > diff((-3.1^2*x*(x-1)*(3.1*x^2-3.1*x+1)),x) 0., 0.6774193548, 0.5580141252, 0.7645665200 / 2 \ / 2 \ -9.61 (-1 + x) \3.1 x - 3.1 x + 1/ - 9.61 x \3.1 x - 3.1 x + 1/ - 9.61 x (-1 + x) (6.2 x - 3.1) > #Based on the fact that 3 gives period 2 orbit, it seems as though k>3 gives > period 4 orbit. > Orb:=proc(f,x,x0,K1,K2) local x1,i,L: x1:=x0: for i from 1 to K1 do > x1:=subs(x=x1,f): #we don't record the first values of K1, since we are > interested in the long-time behavior of the orbit od: L:=[x1]: for i from K1 > to K2 do x1:=subs(x=x1,f): #we compute the next member of the orbit > L:=[op(L),x1]: #we append it to the list od: L: #that's the output end: > Orb(3.1*3.1*x*(1-x)*(1-3.1*x*(1-x)),x,0.5,1000,1020); #My previous prediction > was seemingly wrong since there is no oscillation between two points in this > situation. [0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252, 0.5580141252] > #However, for some reason k = 3.5 seems to be the point of bifurcation based > on the orb function. > Orb(3.4*3.4*x*(1-x)*(1-3.4*x*(1-x)),x,0.5,1000,1020) [0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476, 0.4519632476] > Orb(3.5*3.5*x*(1-x)*(1-3.5*x*(1-x)),x,0.5,1000,1020) [0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828, 0.5008842103, 0.3828196828]