# Math 2250 # Project 2 # February 1999 # # 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 Math 2250 Maple home page # http://www.math.utah.edu/~korevaar/2250spring99.html 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 should make sure they # are part of a single execution group. You can see if they are by # checking that there is a bracket on the left of the worksheet which # groups them all together. If this is not the case, you blacken the # commands in question 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 be asked to do several problems # from section 4.3 of the text. # # Part A : Numerical techniques on a model probem # 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}); # Problem 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. ) # # Problem 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]. # # 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). # The next sequence of commands will create an approximate solution with # n=20 time steps, for the same intial value problem we just studied. # As you hit the return key through the various commands try to slow # yourself down and understand what each command is doing. > 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 > f := (x, y, t) -> y > g:=(x,y,t)->-x; #same g g := (x, y, t) -> -x > n:=20;T:=evalf(2*Pi);h:=evalf(T/n); > #20 steps from 0 to 2*Pi n := 20 T := 6.283185308 h := .3141592654 > x:=vector(n+1);y:=vector(n+1); x := array(1 .. 21, []) y := array(1 .. 21, []) > x[1]:=1;y[1]:=0; #initial conditions x[1] := 1 y[1] := 0 # # Here's Runge-Kutta. Compare to chapter 1: it's just been vectorized. # The evalf's have been added to avoid recursion and to speed up the # code. > for i from 1 to n do > ti:=(i-1)*h: #current time > xi:=x[i]: yi:=y[i]: #current x and y > k1:=f(xi,yi,ti): #left-hand slope > l1:=g(xi,yi,ti): #left-hand slope > k2:=f(xi+h*k1/2,yi+h*l1/2,ti+h/2): #1st midpoint slope > l2:=g(xi+h*k1/2,yi+h*l1/2,ti+h/2): #other component > k3:=f(xi+h*k2/2,yi+h*l2/2,ti+h/2): #2nd midpoint slope > l3:=g(xi+h*k2/2,yi+h*l2/2,ti+h/2): #other component > k4:=f(xi+h*k3,yi+h*l3,ti+h): #right-hand slope > l4:=g(xi+h*k3,yi+h*l3,ti+h): #other component > k:= (k1+2*k2+2*k3+k4)/6: #x-comp of displacement > l:= (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}); # Problem A3) Compare the accuracy of Euler with n=10, Euler with # n=100, and Runge-Kutta with n=20, in terms of how well they # approximate the exact values of x and y when T=2*Pi. You should # exhibit the numerical values of the three approximations, as well as # the exact solution [1,0] corresponding to that T. # # Part B: Some Problems from 4.3: Your instructor may pick (some of) # these problems, or others. # # Problem B1) Do #9 page 250. In addition to what the text asks for, # create a plot which contains the tangent vector field, the exact # solution curve for t between 0 and 1, and the Runge-Kutta # approximation points when h=0.05. (See the commands above for the # various pieces you will need.) # # Problem B2) Do #11 page 250. Notice this system is not autonomous, # so the tangent vector field is changing in time as well as in space. # Make a table comparing the exact solution to the h=0.1 and h=0.05 # Runge-Kutta approximations, at t-values 0, 0.2, 0.4, 0.6, 0.8, 1.0. # # Problem B3) Do #13 page 250. You will have to convert this second # order equation into a first order system, like we did in part A. # # Problem B4) Use the methods section 2.3 page 88 to find the exact # solution to B3), and compare with your numerical work.