Amath 301 Homework 4: Due Friday 8/12 at 4pm Summer 2016 Warning: This homework is long, do not wait until the last minute to start or you will not finish. Also, given its length, it is probably a good idea to clear all between exercises to free up memory Exercise 1: Chemical Differential Equations Consider the chemical reaction k1 [C] X Y. The molecule X is converted to Y . There is a catalyst [C] that controls the rate of the reaction, such that the rate of conversion of X to Y is k1 [C] molecules/second. The catalyst itself, however, is affected by the amount of X present. Namely, [C] = [X] . k2 + [X] Finally, suppose we have an external machine that is adding or removing X periodically, where the time-dependent rate is given by k3 sin(t). Then the differential equation governing the concentration of [X] is: d[X] = k1 [C][X] + k3 sin(t). dt Suppose k1 = 2, k2 = 0.25, and k3 = 1.5. Suppose also that we begin with [X] = 1 at t = 0. We wish to compute [X] as a function of t from t = 0 to t = 1. (a) Integrate the differential equation to t = 1 using the forward Euler's method with a time step of h = 0.1. Save the resulting values ([X](t = 0), [X](t = 0.1), . . . , [X](t = 1))T as a column vector in A1.dat. (b) Integrate the differential equation to t = 1 using the fourth order Runge-Kutta method with a time step of h = 0.1. Save the resulting values ([X](t = 0), [X](t = 0.1), . . . , [X](t = 1))T as a column vector in A2.dat. (c) Repeat exercises (a) and (b) with a time step of h = 0.01. Save the results in a two-column matrix where the first column contains the solution from the forward Euler's method and the second column contains the solution from the fourth order Runge-Kutta method. Save this matrix in A3.dat (d) Integrate the differential equation to t = 1 using the Matlab function ode45 with tspan=[0;1] and default settings otherwise. Save the output t and [X], in that order, in a two-column matrix in A4.dat. 1 Amath 301 Homework 4: Due Friday 8/12 at 4pm Summer 2016 Try plotting your results on the same plot to see how they compare. No need to submit comments, but explore a bit! Exercise 2: A Second Order Differential Equation Consider the damped pendulum equation, which describes the motion of a pendulum that has some impediment hindering its movement (the hinge may need oil, or there might be friction in the air, etc.). The motion of this pendulum is described by: g x00 (t) = sin(x(t)) x0 (t), L where x(t) is the angle of displacement, g is the acceleration due to gravity, L is the length of the pendulum, and is a damping coefficient. Note that this equation is nonlinear. However, if we assume that x \u001c 1, then we can approximate sin(x(t)) with x(t) (think about the Taylor series for sin(x)). This gives us: g x00 (t) = x(t) x0 (t). L Now that the equation is linear, it can be written in terms of a linear system by introducing v(t) = x0 (t), i.e. \u0015 \u0015 \u0014 \u0015\u0014 \u0014 d x(t) A11 A12 x(t) . = A21 A22 v(t) dt v(t) (a) Find the matrix A if g = 10, L = 8, and = 3. Save it in A5.dat. (b) Recall that the forward Euler method for a linear system of differential equations can be written as an iteration: \u0014 \u0015 \u0014 \u0015 xn+1 x = (I + tA) n , vn+1 vn and that the stability of this method depended upon the spectral radius (maximum size of the eigenvalues) of M = I + tA. Use these concepts, and the matrix A found above, to determine the maximum time step tmax we can take and still have the forward Euler method be stable. You should do this by hand. Save the value you find in A6.dat. (c) Solve the linear system of differential equations, with the matrix A above, using the forward Euler method from t = 0 to t = 50. For initial conditions, use [x(0), v(0)]T = [1, 0]T , and use a timestep of 0.8tmax . Plot the displacement x(t) to verify that your solution is stable (i.e. it does not blow up to very large values). Store the resulting x(t) in a row vector. Save this as A7.dat. 2 Amath 301 Homework 4: Due Friday 8/12 at 4pm Summer 2016 (d) Now solve the same linear system with the same initial conditions using forward Euler again, but use t = 1.05tmax instead (only five percent above the maximum). Plot the results to verify that the scheme is no longer stable. Store the resulting x(t) in a row vector. Save this as A8.dat. (e) Solve the same problem with the same initial conditions using the built-in solver ode45. As before, solve this from t = 0 to t = 50. Use default settings otherwise. Save the computed x(t) values as a column vector in A9.dat. Save the computed v(t) values as a column vector in A10.dat. Try comparing these result with your results from (c) and (d). Does your forward Euler scheme perform better if you reduce the time step even further? Explore further, but do not worry about submitting commentary. 3 Amath 301 Homework 4: Due Friday 8/12 at 4pm Summer 2016 Exercise 3: A Chaotic ODE The following three-dimensional autonomous chaotic ODE was proposed and analyzed in a recent paper, Analysis of a novel three-dimensional chaotic system (Li et al., 2013): x = ax + f yz y = cy dxz z = bz + ey 2 Where x, y, z are the state variables, and a, b, c, d, e, f are positive constant parameters. You will now model this system, reproducing a few of their results (the paper is avaliable at http://www.sciencedirect.com/science/article/pii/S0030402612002823?np=y, you should check it out). Much of this work will parallel that done in video lecture 23. (a) You will begin by creating a cube of initial values. The cube should be oriented in the natural manner; that is, with edges parallel the x, y, and z axes. Each edge should be .4 long. The spacing between grid points should be .05. The center of the cube is at (x, y, z) = (.3, 2, 12). Using meshgrid() create three three-dimensional arrays containing the x-coordinates, y-coordinates, and z-coordinates, respectively. To avoid dealing with high dimensional arrays (as is done in the online video lectures), rescale each of these matrices as a vector (a matrix, A , can be rescaled as a vector by V = A(:)). Store the resultant vectors as the rows in a three row matrix with the x, y, and z coordinate vectors in the 1st, 2nd, and 3rd rows, respectively. Save the matrix as A11.dat. To check you work, try the command plot3(x,y,z,'.') and verify that you have a cube of points spaced as described above. (b) Let a = 16, b = 5, c = 10, d = 6, e = 18, f = .05. Model the above system with Runge-Kutta 4 (you can reuse some of your code from Homework 4, Problem 3), from time 0 to time 7, using a time step t = .01, and with the initial condition cube created in (a). Structure your solution in the same way as your initial condition, however you should now have a 3-dimensional array (size = 701, 3, 729). If your output is called A, then A(1,:,:) will be your initial condition cube, A(2,:,:) will correspond to the positions of those initial points after t seconds, and so on. Be sure to properly vectorize your code, or it will run extremely slowly (and will probably terminate prior to completion on scorelator). Once completed, use the squeeze() function to make the 3-D array of solutions into 3 2-D matrices containing the x, y, or z coordinates at the different time steps as the rows of these matrices. Type help squeeze to learn how to use this command. Save the matrix of x-coordinates as A12.dat, the matrix of y-coordinates as A13.dat, and the matrix of z coordinates as A14.dat. The size of each of these should be 701 729 (the x, y, and z coordinates of the 729 points from time 0 to time 7 in steps of 0.01). 4 Amath 301 Homework 4: Due Friday 8/12 at 4pm Summer 2016 (c) You will not submit anything more for this problem, but you should construct some cool plots! Plot your solutions at the last time using the code: plot3(A12(end,:),A13(end,:),A14(end,:),'.'); where A12, A13, and A14 are what you created in part b). Verify that your solution looks similar to the plots of the system found in the paper referenced above (this is also a good way to help you troubleshoot if things are not working out). Another interesting plot you should make is of the tracks of a few of your initial positions: plot3(A12(:,1:5),A13(:,1:5),A14(:,1:5)); You can increase the number of tracks you plot (you have one for each of the 93 starting positions, but you only), but note that this will quickly become quite taxing on your computer. Finally, try watching the full evolution of the system by typing: for i=1:length(t) plot3(A2(i,:),A3(i,:),A4(i,:),'.'); axis([-.5,.5,-10,10,0,50]) drawnow pause(h) end Make sure to comment out your plots before submitting the assignment. 5 Amath 301 Homework 4: Due Friday 8/12 at 4pm Summer 2016 Exercise 4: A Chaotic Map A map is a system which progresses in discrete steps according to a function f , such that we have xn+1 = f (xn ) (this is opposed to the ODEs we have studied, which advance continuously according to a differential equation). In this class, we have seen maps in the context of solving ODE's, where we have approximated the solution to the continuous systems with similarly constructed maps (such as the forward Euler or Runge Kutta methods). But maps are worth studying in their own right. In this problem we will be looking at the famous logistic map: xn+1 = rxn (1 xn ) Despite its simplicity, it exhibits complex and chaotic behavior. a) With x0 = .7 calculate 30 iterations of the logistic map where r takes on the values 2.5, 3.2, 3.52, and 4. Store these as the columns of a 31 4 matrix (x0 is included, so the first row of this matrix should be equal to 0.7 in every entry). Each column of the matrix corresponds to one of the r values in the order given. Save this matrix as A15.dat Try plotting each of these solutions (make sure to comment this out before submitting). For r = 2.5, 3.2, and 3.52, you should notice that the solutions eventually begin to repeat themselves (in the language of dynamical systems, they \"converge to a periodic orbit\"). We will call the \"period\" of an orbit the number of iterations it takes to return to where it was before. For example, the sequence \"1,2,1,2,1,. . . \" would have period two, while the sequence \"1,2,6,4,5,1,2,6,4,5,...\" would have period five. Define PR as the period of the orbit the sequence of the logistic map approaches when r = R. Looking at the plots of your results, identify the periods when r = 2.5 (that is, P2.5 ), when r = 3.2, and when r = 3.52. Save these in a column vector [P2.5 ; P3.2 ; P3.52 ] as A16.dat. Notice that with r = 4, the solution does not appear to have any regular behavior. b) Now find the first 500 iterations of the map with x0 = .7 under a whole range of r values: r = 2.5:.001:4. These should be stored in a 501 1501 matrix, in which each column represents the sequence from x0 to x5 00 for a different r. Make sure your code is properly vectorized (there should only be a single for loop). Save the last 100 iterations (a 100 1501 matrix) in A17.dat. To see the nature of the long-term solutions with respect to r, plot these last hundred iterations using this code: plot(R,A7,'k.','Markersize',1); You have made what is known as the bifurcation diagram of the logistic map, a celebrated image from the history of chaos theory. Notice how the map first becomes fractal-like, repeatedly splitting, before eventually giving rise to chaos as r increases. Make sure to comment out the plot before submitting. 6