# Math 2280 # Maple Project 2 # Summer 2001 # # # Instruction: You are expected to read relevant sections from the book # and work through the project. The only work which you have to turn in # consists of those questions which are clearly labeled as PROBLEMS. # There are many other times when you are asked to verify certain # computations or facts, but these verifications are for your own benefit # in understanding the material. Problem 1-4 are warm-up questions. # Save a copy of the file using your web browser, and then open it from # Maple as a ``Maple text'' document. # # In this project we will study numerical methods for # approximating solutions to first order differential equations. We # will see soon how higher order differential equations can be converted # into first order systems of differential equations # We will be working through material from sections 2.4-2.6 of the # text. You may want to read the text and consult the Computing # Projects manual, projects for chapter 2. # if you find explanations here to be deficient. # The most basic method of approximating solutions to differential # equations is called Euler's method, after the 1700's mathematician who # first formulated it. If you want to approximate the solution to the # initial value problem dy/dx = f(x,y), y(x0)=y0, first pick a step size # ``h''. Then for x between x0 and x0+h, use the constant slope # f(x0,y0). At x-value x1:=x0+h your y-value will therefore be y1:=y0 + # f(x0,y0)h. Then for x between x1 and x1+h you use the constant slope # f(x1,y1), so that at x2:=x1+h your y-value is y2:=y1+f(x1,y1)h. You # continue in this manner. It is easy to visualize if you understand # the slope field concept we've been talking about; you just use the # slope field where you are at the end of each step to get a slope for # the next step. It is straightforward to have the computer do this # sort of tedious computation for you. In Euler's time such # computations would have been done by hand! # # A good first example to illustrate Euler's method is our favorite # DE from the time of Calculus, namely dy/dx =y, say with initial value # y(0)=1, so that y=exp(x) is the solution. Let's take h=0.2 and try to # approximate the solution of the x-interval [0,1]. Since the # approximate solution will be piecewise affine, we only need to know # the approximations at the discrete x values x=0,0.2,0.4,0.6,0.8,1. # Here's a simple ``do loop'' to make these computations. # > restart; #clear any memory from earlier work > x0:=0.0; xn:=1.0; y0:=1.0; n:=5; h:=(xn-x0)/n; # specify initial > # values, number of steps, step size. > x0 := 0 xn := 1.0 y0 := 1.0 n := 5 h := .2000000000 > f:=(x,y)->y; #this is the function f(x,y) in dy/dx = f(x,y) f := (x, y) -> y > x:=x0; y:=y0; #initialize x,y for the do loop x := 0 y := 1.0 > for i from 1 to n do > k:= f(x,y): #current slope, use : to suppress output > y:= y + h*k: #new y value via Euler > x:= x + h: #updated x-value: > print(x,y,exp(x)); #display current values, > #and compare to exact solution. > od: #``od'' ends a do loop .2000000000, 1.200000000, 1.221402758 .4000000000, 1.440000000, 1.491824698 .6000000000, 1.728000000, 1.822118800 .8000000000, 2.073600000, 2.225540928 1.000000000, 2.488320000, 2.718281828 # Notice your approximations are all a little too small, in particular # your final approximation 2.488... is short of the exact value of # exp(1)=e=2.71828.. The reason for this is that because of the form of # our f(x,y) our approximate slope is always less than the actual slope. # We can see this graphically using plots: > with(plots):with(linalg): > xval:=vector(n+1);yval:=vector(n+1); #to collect all our points xval := array(1 .. 6, []) yval := array(1 .. 6, []) > xval[1]:=x0; yval[1]:=y0; #initial values xval[1] := 0 yval[1] := 1.0 > #paste in the previous work, and modify for plotting: > for i from 1 to n do > x:=xval[i]: #current x > y:=yval[i]: #current y > k:= f(x,y): #current slope > yval[i+1]:= y + h*k: #new y value via Euler > xval[i+1]:= x + h: #updated x-value: > od: #``od'' ends a do loop > approxsol:=pointplot({seq([xval[i],yval[i]], i=1..n+1)},connect=true): > exactsol:=plot(exp(t),t=0..1,`color`=`black`): > #used t because x was already used above > display({approxsol,exactsol}); # The picture you just created is like the one in figure 2.4.3, on page # 100 of the text. In your picture the upper graph is of the # exponential function, the lower graph is of your Euler approximation. # # It should be that as your step size ``h'' gets smaller, your # approximations get better to the actual solution. This is true if # your computer can do exact math (which it can't), but in practice you # don't want to make the computer do too many computations because of # problems with round-off error and computation time, so for example, # choosing h=0.0000001 would not be practical. But, trying h=0.01 should # be instructive. # # PROBLEM 1. Change the n-value to 100. Keep the other data the same. # Only print out your approximations when h is a multiple of 0.1. The # commands you need are below. > x0:=0.0; xn:=1.0; y0:=1.0; n:=100; h:=(xn-x0)/n; > x:=x0; y:=y0; > for i from 1 to n do > k:= f(x,y): #current slope > y:= y + h*k: #new y value via Euler > x:= x + h: #updated x-value: > if frac(i/10)=0 > then print(x,y,exp(x)); > fi; #use the ``if'' test to decide when to print; > #the command ``frac'' computes the remainder > #of a quotient, it will be zero for us if i > #is a multiple of 10. > od: .1000000000, 1.104622127, 1.105170918 .2000000000, 1.220190041, 1.221402758 .3000000000, 1.347848916, 1.349858808 .4000000000, 1.488863735, 1.491824698 .5000000000, 1.644631823, 1.648721271 .6000000000, 1.816696699, 1.822118800 .7000000000, 2.006763370, 2.013752707 .8000000000, 2.216715220, 2.225540928 .9000000000, 2.448632678, 2.459603111 1.000000000, 2.704813834, 2.718281828 # PROBLEM 2. How close is your estimate to e now, based on n=100? # # Actually, considering how many computations you did in problem # B.1, you are not so close to the exact solution. In more complicated # problems it is a very serious issue to find relatively efficient ways # of approximating solutions. An entire field of mathematics, # ``numerical analysis'' deals with such issues for a variety of # mathematical problems. The book talks about some improvements to # Euler in sections 2.5 and 2.6. If you are interested in this # important field of mathematics you should read 2.5 and 2.6 carefully. # Let's summarize some highlights below. # For any time step h the fundamental theorem of calculus asserts # that, since dy/dx = f(x,y(x)), > y(x+h)=y(x) + Int(f(t,y(t)),t=`x`..`x+h`); x+h / | y(x + h) = y(x) + | f(t, y(t)) dt | / x # The problem with Euler is that we always approximated this integral by # h*f(x,y(x)), i.e. we used the left-hand endpoint as our approximation # of the ``average height''. The improvements to Euler depend on better # approximations to that integral. ``Improved Euler'' uses an # approximation to the Trapezoid Rule to approximate the integral. The # trapezoid rule for the integral approximation would be # (1/2)*h*(f(x,y(x))+f((x+h),y(x+h)). Since we don't know y(x+h) we # approximate it using unimproved Euler, and then feed that into the # trapezoid rule. This leads to the improved Euler do loop below. Of # course before you use it you must make sure you initialize everything # correctly. # > x:=x0; y:=y0; n:=5; h:=(xn-x0)/n; x := 0 y := 1.0 n := 5 h := .2000000000 > for i from 1 to n do > k1:=f(x,y): #left-hand slope > k2:=f(x+h,y+h*k1): #approximation to right-hand slope > k:= (k1+k2)/2: #approximation to average slope > y:= y+h*k: #improved Euler update > x:= x+h: #update x > print(x,y,exp(x)); > od: > .2000000000, 1.220000000, 1.221402758 .4000000000, 1.488400000, 1.491824698 .6000000000, 1.815848000, 1.822118800 .8000000000, 2.215334560, 2.225540928 1.000000000, 2.702708163, 2.718281828 # Notice you almost did as well with n=5 as you did with n=100 in # unimproved Euler. # # PROBLEM 3. Use improved Euler with n= 100 to see how well you do # for the same initial value problem. Use the template above, modified # to print out data only at x-multiples of 0.1 # # One can also use Taylor approximation methods to improve upon # Euler; by differentiating the equation dy/dx = f(x,y) one can solve # for higher order derivatives of y in terms of the lower order ones, # and then use the Taylor approximation for y(x+h) in terms of y(x). # See the book for more details of this method, we won't do it here. # In the same vein as ``improved Euler'' we can use the Simpson # approximation for the integral instead of the Trapezoid one, and this # leads to the Runge-Kutta method which is widely used in real software. # (You may or may not have talked about Simpson's Rule in Calculus, it # is based on a quadratic approximation to the function f, whereas the # Trapezoid rule is based on a first order approximation.) Here is the # code for the Runge-Kutta method. The text explains it in section 2.6. # Simpson's rule approximates an integral in terms of the integrand # values at each endpoint and at the interval midpoint. Runge-Kutta # uses two different approximations for the midpoint value. # Before you use the loop below you must initialize your values, # like you did above. # # > > x:=x0; y:=y0; n:=5; h:=(xn-x0)/n; x := 0 y := 1.0 n := 5 h := .2000000000 > for i from 1 to n do > k1:=f(x,y): #left-hand slope > k2:=f(x+h/2,y+h*k1/2): #1st guess at midpoint slope > k3:=f(x+h/2,y+h*k2/2): #second guess at midpoint slope > k4:=f(x+h,y+h*k3): #guess at right-hand slope > k:=(k1+2*k2+2*k3+k4)/6: #Simpson's approximation for the integral > x:=x+h: #x update > y:=y+h*k: #y update > print(x,y,exp(x)); #display current values > od: > .2000000000, 1.221400000, 1.221402758 .4000000000, 1.491817960, 1.491824698 .6000000000, 1.822106456, 1.822118800 .8000000000, 2.225520825, 2.225540928 1.000000000, 2.718251136, 2.718281828 # Notice how close Runge-Kutta gets you to the correct value of e, with # n=5. # # PROBLEM 4. Do Runge-Kutta for the same initial value problem, with # n=20. Compare to your results using Euler and improved Euler, using # n=100. There should be a huge improvement. # # As we know, solutions to non-linear DE's can blow up, and there are # other interesting pathologies as well, so if one is doing numerical # solutions there is a real need for care. The text has a number of # examples of this. You are to do the following two problems: Don't # forget to insert textual explanations to go with your MAPLE # computations. # # PROBLEM 5. Section 2.4 page 103, #25. (illustrates the danger of # missing blow-up) # # PROBLEM 6. Section 2.6 page 134: Work through and understand # example #4. Recreate the data in the last two columns of the table on # the top of page 134, using your Runge-Kutta routine (i.e. the # Runge-Kutta approximation with h=0.05 and the actual y value). # # Understand (for yourself) how mathematical instability caused by # the solution to the homogeneous problem can lead to numerical # instability like that illustrated in this important example of Problem # B.6, even when you're using a ``good'' numerical technique like Runge # Kutta. This example illustrates why people who create numerical # solutions to problems must also be well-grounded in the mathematical # theory behind the problems. # #