# Math 2250 # Project 2 # October, 1998 # # Introduction: This project has two parts. In Part A you will use # Euler's method for an example problem. This will illustrate why # ``nice'' first order systems of ordinary differential equations should # have unique solutions for initial value problems. We won't be precise # about what ``nice'' means - you can look in one of the appendices to # our text for more details. Our reasoning in this section is like the # disucussion of existence and uniqueness for solutions to initial value # problems which we had previously, in Chapter 1 and the first Maple # project: For first order DE's the solution graph is tangent to the # slope field determined by f. For first order systems of differential # equations, the solution curve is tangent to the vector field # determined by f. # In section 4.1 of the text it is indicated how a single nth order # differential equation can be reformulated as a system of n first order # systems differential equations: getting to specify the initial values # of the function and its first (n-1) derivatives at a given t-value # corresponds to getting to specify an initial value for the vector made # out of the function and its first (n-1) derivatives, in the # reformulated first-order system. Therefore existence and uniqueness # for initial value problems in the context of systems of first order # differential equations, explains existence and uniqueness initial # value problems in the context of nth order differential equations, # i.e. the Theorem we discussed in section 3.1. # After doing Euler to help you understand the theory of systems of # first order differential equations, you will use the Runge-Kutta # algorithm on the same model problem, to illustrate the fact that this # algorithm is much more effective computationally than Euler. # In Part A of the project you will be led along by the commands # which are already written for you. Therefore, it is strongly # recommended that you get the web version of this document by saving it # from the address http://www.math.utah.edu/~korevaar/2250proj2.txt # After you save it to your directory as a *.txt document you can open # it from Maple, if you tell Maple to look for a ``Maple Text'' # document. You make this choice from the bottom of the ``open'' dialog # box. Then when you want to execute several lines of code at once, as # you will want to do with the do-loops for Euler and Runge Kutta, you # blacken them with your mouse, go up to the ``Edit'' menu choice, and # pick the ``join execution groups'' option. This will put a big # bracket around your black-lighted section; then you can execute # everything in it at once. # In Part B of the project you will use Runge-Kutta to study the # Duffing Equation: This is a forced oscillation problem with a # non-linear spring, with force function -(kx + bx^3) instead of -kx. # It is discussed in section 6.5 of the text, starting at page 390. The # only difference between this equation and what we already understand # is the bx^3 term. This turns out to be a big change, as you will see. # In Part B you will be expected to provide more thinking and more code # on your own. # # Part A : Numerical techniques # Geometrically a system of first order DE's is telling you that # your tangent vector is determined by where you are at a certain # instant. In Euler for systems you compute the tangent vector from # where you are at stage i, and then for small time step ``h'' you use # this constant tangent. So if you were at the point x[i] at time t[i], # and if you DE is dx/dt=f(t,x), then at time t[i]+h=t[i+1] you will be # at x[i]+h*f(t[i],x[i])=x[i+1]. This makes for a short do-loop # iteration. We will try Euler on the elementary initial value problem [ dx ] [----] [ dt ] [y ] [ ] = [ ] [ dy ] [-x] [----] [ dt ] [x(0)] [1] [ ] = [ ] [y(0)] [0] # Notice that if you convert the undamped, unforced, spring # equation # 2 d x ---- + x = 0 2 dt # with initial conditions x(0)=1, dx/dt(0)=0, to a first order system in # which the first function is x(t) and the second function y(t)=dx/dt, # you get exactly the initial value problem and system which we are # studying. Therefore, since we can solve the spring initial value # problem and deduce that x(t)=cos(t), we know that the solution to our # system of first order DE's is x(t)=cos(t), y(t)=-sin(t). If we draw # the resulting path (x(t),y(t)) in the plane we will get the # unit-radius circle, going around clockwise. # Notice, (and check), that no matter which initial value problem # we took for this system, our solution curve would be (Acos(t-a), # -Asin(t-a)), i.e. a circle going clockwise around the origin. Its # radius is exactly the amplitude of the corresponding solution to the # spring problem. # A picture of several of the solution curves, together with a # sample of the tangent vectors at various points in the plane would be # called a phase diagram for the spring problem. MOre generally for a # system of two first order autonomous (no t-dependence) differential # equations for x(t) and y(t), the x-y plane is called the phase plane, # and a phase diagram would consist of the tangent vector field and # pictures of several solution curves (x(t),y(t)). Notice that the use # of the word ``phase'' here does not seem to have a whole lot to do # with its previous use in chapter 3. That is unfortunate and leads to # confusion. Chapter 6 deals in more depth with phase planes. # Let's return to our model problem. We know what our solution # curves should look like, let's see how Euler does in comparison: > restart:with(linalg):with(plots):with(DEtools): > #load up the tools! > f:=(x,y,t)->y; #f is the first comp. in our > #autonomous DE, in our case > #dx/dt=y, so f(x,y,t)=y f := (x, y, t) -> y > g:=(x,y,t)->-x; #g is the 2nd comp. in our DE, > #i.e. dy/dt=-x g := (x, y, t) -> -x > n:=10;T:=2.0*Pi;h:=T/n; #we will do Euler > #for t=0..2*Pi, 10 time steps n := 10 T := 2.0 Pi h := .2000000000 Pi > x:=vector(n+1);y:=vector(n+1); > #to hold x and y coords at step i x := array(1 .. 11, []) y := array(1 .. 11, []) > x[1]:=1.0;y[1]:=0.0; #initial values x[1] := 1.0 y[1] := 0 # Now comes the simple Euler do-loop: > for i from 1 to n do > ti:=h*(i-1): #the time is ti > x[i+1]:=evalf(x[i]+h*f(x[i],y[i],ti)): > y[i+1]:=evalf(y[i]+h*g(x[i],y[i],ti)): > od: # We have inserted the evalf command above to prevent MAPLE from trying # to do exact arithmetic; otherwise MAPLE would try to define everything # recursively and if n was too large, say n=100, x[101] would be # defined in terms of x[100] which would defined in terms of x[99], etc, # etc, etc, and your computer would run out of memory. With the evalf # command MAPLE makes a decimal approximation at each stage. # # We will plot the tangent field, the approximate solution, and the # exact solution: > approxsol:=pointplot({seq([x[i],y[i]],i=1..n+1)},connect=true): > #connect=true makes Maple connect the points > #you get from Euler with line segments. It only > #works well with small numbers of points. > exactsol:=plot([cos(t),-sin(t),t=0..2*Pi]): > velocfield:=fieldplot([x2,-x1],x1=-2..2,x2=-2..2): > #We also think of the tangent field as a > #velocity field, hence the particular choice > #of name. > display({approxsol,exactsol,velocfield}); # #A1) Explain the picture above. Can you explain geometrically why # your Euler approximation is spiraling away from the exact circle # solution? By the way, every now and then Maple seems to connect the # approximate solution points incorrectly, so that the last point is # connected to the first one. Hopefully that doesn't happen to you. # # #A2) Redo Euler with n=100 time steps for the same time interval. # Create a display like you did with n=10. When you make the pointplot, # do NOT set connect=true: MAPLE seems to have real trouble with this # option when there are a lot of points . You will see the spiraling # still, but it won't be as bad. Compare the Euler approximation for # [x(2*Pi),y(2*Pi)] with the correct values (1,0) by evaluating # x[101],y[101]. # # #A3) Here is the code for Runge-Kutta. As with the Euler code we # wrote it for a general first order system of two DE's, even though our # current f and g functions are autonomous (don't use the variable t). # In Part B of the project you will want to allow t-dependence, however. # Use it and the rest of the commands to create a picture with n=20, # and compare your answer to that of Euler with n=100. In particular, # compute your approximation for [x(2*Pi),y(2*Pi)], and compare it to # [1,0]. Here are the commands to do most of what is asked: > restart;with(linalg):with(DEtools):with(plots): > #restarting clears out memory. It's not > #a bad idea to do this when you start a new > #portion of a project, especially if you > #want to use the same variable names for > #different things. > f:=(x,y,t)->y; #same f > > g:=(x,y,t)->-x; #same g > n:=20;T:=2*Pi;h:=evalf(T/n); > #20 steps from 0 to 2*Pi > x:=vector(n+1);y:=vector(n+1); > x[1]:=1;y[1]:=0; #initial conditions # # Here's Runge-Kutta. Compare to chapter 1: it's just been vectorized. # The evalf's have been added to avoid recursion in some places and just # to speed up the code in others. > for i from 1 to n do > ti:=(i-1)*h: #current time > xi:=evalf(x[i]):yi:=evalf(y[i]): #current x and y > k1:=evalf(f(xi,yi,ti)): #left-hand slope > l1:=evalf(g(xi,yi,ti)): #left-hand slope > k2:=evalf(f(xi+h*k1/2,yi+h*l1/2,ti+h/2)): #1st midpoint slope > l2:=evalf(g(xi+h*k1/2,yi+h*l1/2,ti+h/2)): #other component > k3:=evalf(f(xi+h*k2/2,yi+h*l2/2,ti+h/2)): #2nd midpoint slope > l3:=evalf(g(xi+h*k2/2,yi+h*l2/2,ti+h/2)): #other component > k4:=evalf(f(xi+h*k3,yi+h*l3,ti+h)): #right-hand slope > l4:=evalf(g(xi+h*k3,yi+h*l3,ti+h)): #other component > k:= evalf((k1+2*k2+2*k3+k4)/6): #x-comp of displacement > l:= evalf((l1+2*l2+2*l3+l4)/6): #y-comp of displacement > x[i+1]:=evalf(x[i]+h*k): #updated x-coord > y[i+1]:=evalf(y[i]+h*l): #updated y-coord > od: > approxsol:=pointplot({seq([x[i],y[i]],i=1..n+1)}): > exactsol:=plot([cos(t),-sin(t),t=0..2*Pi]): > display({approxsol,exactsol}); # Notice how closely the points match the circle. Now go back and # answer the rest of A3. # # Part B: Duffing's Equation and Period Doubling in dynamical systems # Because you are going to use a lot of computing in this section, # and because our lab machines don't have great memories, it is # recommended that you create a new file to do Part B. One way to do # that, if you are working from the web version, is to use the File menu # item to now ``Save As'' another, new filename. Then your work up # until now will be saved in your old-named file, and you will have a # current copy to work with in you newly- named file. You can then # delete the Part A stuff from your new file that you don't want for # part B. Make sure to keep the stuff related to Runge Kutta. # Hopefully you are thinking about the particular mass-wire model # which the book describes on page 390, figure 6.5.9, as a possible # model for the Duffing ``spring''. We pick b=-1, m=c=k=1 as simple # values of the parameters, so we are studying the equation (14) on page # 391: 2 d x dx 3 ---- + ---- - x + x = F0 cos(t) 2 dt dt # Notice that if we consider the unforced system (F0=0) there are now # three equilibrium solutions (equilibrium solutions are constant # solutions which solve the DE, just like in Chapter 2), namely the # roots of -x + x^3 =0. These are x=0, x=1, x=1. This is qualitatively # what we would expect from the configuration in figure 6.5.9; in fact # the three equilibrium solutions are represented in that picture. # Since you would expect the middle one to be unstable (configurations # starting near it will not stay near it!), x=0 is unstable and x=1,-1 # are stable equilibria, we expect. In fact, if we ``linearize'' near # x=0, considering small displacements x, we have a negative spring # constant in Hooke's law! Well, this is consistent with the # configuration: small displacements lead to forces which pull in the # direction of the displacement. # We will actually focus near the stable solution x=1, for this # problem. First let's try to linearize near this equilibrium. We # linearize by changing from the function x(t) to the function # z(t)=x(t)-1. You do this change because if you are interested in # x(t)'s near 1 you are interested in z(t)'s near 0. You should verify # that if you plug x(t)=z(t)+1 into our equation above you get 2 d z dz 2 3 ---- + ---- + 2 z + 3 z + z = F0 cos(t) 2 dt dt # Thus, for solutions near the equilbrium x=1 (z=0), we expect behaviour # similar to that of solutions to the linearized equation 2 d z dz ---- + ---- + 2 z = F0 cos(t) 2 dt dt # #B1) Find the page and formula number (and record them here) which # leads us to deduce that the steady-state solution to the linear DE we # just derived is F0 cos(t - 1/4 Pi) z := t -> ------------------ sqrt(2) # (Hint: This was covered in section 3.6) You should verify that our # formula above is correct, by indicating the values of the various # terms section 3.6 formula, in our case. # # #B2) Convert the second order equation at the start of Part B ((14) on # page 391 of text) into a first order system, by setting y=dx/dt. # Verify that this leads to the first-order system [ dx ] [----] [ dt ] [ y ] [ ] = [ ] [ dy ] [ 3 ] [----] [x - x - y + F0 cos(t)] [ dt ] # # #B3) Since there is damping in the problem, we expect to get # steady-state solutions as t->infinity. At least that's what happened # in section 3.6 for the linear version of this problem. Actually, that # will still turn out to be true when the forcing amplitude F0 is small # enough, but it will turn out to be gloriously false at other values of # F0. # Let's see how Runge-Kutta does for a small F0, and compare our # numerical solution to the linearized solutions from B1. The fact that # we have damping means that our equilibrium solution (if it exists) is # ``stable'', and this has the effect that Runge-Kutta will work well; # the kind of pathology illustrated at the end of the first Maple # project (text, page 123) will not occur for us. We will pick the # equilibrium solution x=1, y=dx/dt=0 as our initial value. The forcing # term has period 2*Pi; we will let our code run for 100 time units to # reach steady state. Our plot will tell us if this was long enough. # You are to obtain the pictures which are displayed below, using # the Runge-Kutta code from Part A. Of course you must redefine the # functions f and g to reflect the current system. These pictures were # obtained with F0=0.1, T=200, n=1000, h=T/n=0.2. Only the last 500 # points were plotted, to correspond with t=100..200. In the plot which # shows the oscillation x(t) at time t, notice that t=i*h. If the # forcing amplitude is small we expect to get solution curves which look # like circles, with radii of about 0.7 times the forcing amplitude. # (From your work in #B1, since 0.7 is a crude approximation to # 1/sqrt(2)). Thus in case F0=0.1 we expect approximate circles of # radius about 0.07. # In the interests of honesty, we point out that one can use Maple # commands like DEsolve to get numerical solutions for the Duffing # Equation; we are making you use Runge-Kutta in order to maximize your # learning in this project. DEsolve would be using a similar algorithm # anyway. # > approxsol:=pointplot({seq([x[i],y[i]],i=n/2..n+1)}): > #We assume n is even, we are taking data from > #the last half of our time interval. > linapprox:=plot([1+cos(t)/(10*sqrt(2)),sin(t)/(10*sqrt(2)),t=0..2*Pi]) > : > display({approxsol,linapprox}); > justx:=pointplot({seq([i*h,x[i]],i=n/2..n+1)}): > display({justx}); > # So, as we hoped, we get a slightly imperfect circle, with radius right # around 0.07, which compares well with the linearized problem. The # second graph, which is showing how x(t) alone is changing, shows that # the response steady state vibration has period 2*Pi, just like the # driving force. This is what we expected based upon our understanding # of the linear problem. # We emphasize that there is no closed-form way of solving the # nonlinear problem, so it must be solved numerically. However, if we # are near an equilibrium solution we expect the linearized problem to # be a good predictor of actual behavior. # # # #B4) Complete PROJECT 3 on page 392 of the text, using the Runge-Kutta # code as above.