# Before you begin, read this first. # Please click # here # to connect to the answer form for this worksheet. After doing so, # put the cursor on the next restart line to begin the worksheet. # When you get to a homework question, enter the answer both on the # worksheet and on the answer form. You may use cut and paste # operations to simplify this procedure. Whenever you have completed # entering the answer, make sure you place the cursor on the next # command line to continue with the worksheet. Save the worksheet # with your answers when you are done to help you study > restart; > with( plots ): > with( DEtools ): # Math 221-0 # Numerical Methods # The worksheet applies the three numerical methods, Euler, Improved # Euler, and Runge-Kutta, to a differential equation. # The function, > f# , that is the right-hand side of the differential equation > f := (t,y) -> y*(1-y^2); / 2\ f := (t, y) -> y \1 - y / # The stepsize > h := 0.1; h := 0.1 # The number of steps to be taken > N := 10; N := 10 # The initial condition > T[0] := 0.0; > Y[0] := 2.0; T[0] := 0. Y[0] := 2.0 # The slope at the initial point (x=1, y=4.0) is > f( T[0], Y[0] ); -6.000 # the approximate value of the solution at the first point (t=0.1) is # approximated to be > T[1] := T[0] + h: > Y[1] := Y[0] + h*f( T[0], Y[0] ): > printf(`t=%1.1f \t y=%1.4f\t k=%1.4f`,T[1], Y[1], f(T[1],Y[1])); # t=.1 y=1.4000 k=-1.3440 # Repeating the same steps four more times yields to give the first # five steps: > T[2] := T[1] + h: > Y[2] := Y[1] + h*f( T[1], Y[1] ): > printf(`t=%1.1f \t y=%1.4f \t k=%1.4f`,T[2], Y[2],f(T[2],Y[2])); # t=.2 y=1.2656 k=-.7616 > T[3] := T[2] + h: > Y[3] := Y[2] + h*f( T[2], Y[2]): > printf(`t=%1.1f \t y=%1.4f \t k=%1.4f`,T[3], Y[3], f(T[3],Y[3])); > T[4] := T[3] + h: > Y[4] := Y[3] + h*f(T[3],Y[3]): > printf(`t=%1.1f \t y=%1.4f \t k=%1.4f`,T[4], Y[4], f(T[4],Y[4])); > T[5] := T[4] + h: > Y[5] := Y[4] + h*f(T[4],Y[4]): > printf(`t=%1.1f \t y=%1.4f \t k=%1.4f`, T[5], Y[5], f(T[5],Y[5])); # Combining the solutions into a table and calculating 10 terms up to # t=1.0 > h := 0.1: > N := 10: > T[0] := 0.0: > Y[0] := 2.0: > printf(` t \t y \t\t\t\t k\n`); > K[0] := f(T[0],Y[0]): > printf(`%1.2f \t%1.5f \t%1.5f \n`,T[0], Y[0],K[0]); > for n from 1 to N do > T[n] := T[n-1] +h: > Y[n] := Y[n-1] + h * K[n-1]: > K[n] := f(T[n],Y[n]): > printf(` %01.2f \t%1.5f \t%01.5f \n`,T[n], Y[n],K[n]); > od: # To understand what has just been computed, let's look at some plots. # First, create the direction field for this problem. > ODE := diff( y(t), t ) = y(t)*(1-y(t)^2); > dir_plot := DEplot( ODE, y, t=0..1,[[y(0)=2]],linecolour=CYAN, y=1..2 > ): # The points computed in Euler's Method can also be plotted > euler_pts := [[T[0],Y[0]], [T[1],Y[1]], [T[2], Y[2]],[T[3],Y[3]], > [T[4],Y[4]], [T[5],Y[5]], [T[6],Y[6]], [T[7],Y[7]], > [T[8],Y[8]], > [T[9],Y[9]], [T[10],Y[10]]]; > euler_plot := plot( euler_pts, color=MAGENTA, thickness=2 ): > display( { euler_plot, dir_plot }, > title=`Euler's Method, h=0.1` ); # Observe how the computed solution closely follows the arrows in the # direction field. # The main deviation comes in the first step from t=0 to 0.2 where # the Euler solution # is not near the direction of the direction field for t=0.09. # We now implement the improved Euler method for the same differential # equation > h := 0.1: > N := 10: > T[0] := 0.0: > IY[0] := 2.0: > printf(` t \t y \t\t\t\t k1 \t k2 \n`); > K1[0] := f(T[0],IY[0]): > Z[0] := IY[0] + h*K1[0]: > T[1] :+ T[0] + h: > K2[0] := f(T[1],Z[0]): > printf(`%1.2f \t%1.5f \t%1.5f \t%1.5f \n`,T[0], IY[0],K1[0],K2[0]); > for n from 1 to N do > IY[n] := IY[n-1] + 0.5*h*(K1[n-1]+K2[n-1]): > K1[n] := f(T[n],IY[n]): > Z[n] := IY[n] + h*K1[n]: > T[n+1] := T[n] +h: > K2[n] := f(T[n+1],Z[n]): > printf(` %01.2f \t%01.5f \t%01.5f \t%01.5f \n`, > T[n],IY[n],K1[n],K2[n]); > od: # Check the above steps to convice yourself that this is the improved # Euler method. # Plot of the solution using improved Euler > improved_euler_pts := [[T[0],IY[0]], [T[1],IY[1]], [T[2],IY[2]], > [T[3],IY[3]], [T[4],IY[4]], [T[5],IY[5]], > [T[6],IY[6]], > [T[7],IY[7]], [T[8],IY[8]], [T[9],IY[9]], [T[10],IY[10]]]: > improved_euler_plot := plot(improved_euler_pts, color=GREEN, > thickness=2 ): > display( { dir_plot,euler_plot, improved_euler_plot }, > title=`Improved Euler's Method, h=0.1` ); # Runge-Kutta for the same differential equation: > h := 0.1: h2 := 0.5*h: > N := 10: > T[0] := 0.0: > RKY[0] := 2.0: > printf(` t \t y \t\t\t\t k1 \t k2 \t k3 \t k4 \n`); > K1[0] := f(T[0],RKY[0]): > Z1[0] := RKY[0] + h2*K1[0]: > K2[0] := f(T[0]+h2,Z1[0]): > Z2[0] := RKY[0] + h2*K2[0]: > K3[0] := f(T[0]+h2,Z2[0]): > Z3[0] := RKY[0] + h*K3[0]: > T[1] :+ T[0] + h: > K4[0] := f(T[1],Z3[0]): > printf(`%1.2f \t%1.5f \t%1.5f \t%1.5f \t%1.5f \t%+1.4f \n`, > T[0], RKY[0],K1[0],K2[0],K3[0],K4[0]); > for n from 1 to N do > RKY[n] := RKY[n-1] + h*(K1[n-1]+2*K2[n-1]+2*K3[n-1]+K4[n-1])/6.0: > K1[n] := f(T[n],RKY[n]): > Z1[n] := RKY[n] + h2*K1[n]: > K2[n] := f(T[n]+h2,Z1[n]): > Z2[n] := RKY[n] + h2*K2[n]: > K3[n] := f(T[n]+h2,Z2[n]): > Z3[n] := RKY[n] + h*K3[n]: > T[n+1] :+ T[n] + h: > K4[n] := f(T[n+1],Z3[n]): > printf(`%01.2f \t%01.5f \t%01.5f \t%01.5f \t%01.5f \t%01.5f \n`, > T[n], RKY[n],K1[n],K2[n],K3[n],K4[n]); > od: # Check the above steps to convice yourself that this is the # Runge-Kutta method # The outputs of the three methods can be compared with the exact # solution. # For the above example the exact solution is > sqrt(C*exp(2*T[n])/(C*exp(2*T[n]-1.0) # where > C = Y[0]^2/(Y[0]^2 - 1.0) > C := Y[0]^2/(Y[0]^2-1): > printf(` t \t Euler \t Improved \t Runge-Kutta \t Exact > Solution\n`); > for n from 0 to N do > Ex[n] := sqrt(C*exp(2*T[n])/(C*exp(2*T[n])-1.0)): > printf(`%1.2f \t %1.5f \t %1.5f \t %1.5f \t %1.5f\n`, > T[n], Y[n], IY[n], RKY[n], Ex[n]); > od: # Plot the solution using Runge-Kutta and the solution from DEplot: > RK_pts := [ [T[0],RKY[0]], [T[1],RKY[1]], [T[2],RKY[2]], > [T[3],RKY[3]], [T[4],RKY[4]], [T[5],RKY[5]], > [T[6],RKY[6]], [T[7],RKY[7]], [T[8],RKY[8]], > [T[9],RKY[9]], [T[10],RKY[10]] ]: > RK_plot := plot(RK_pts, 0..1, 1..2, colour=BLUE, thickness=2 ): > display( {RK_plot, dir_plot},title=`Runge-Kutta Method, h=0.1` ); # Notice that the two solution are so close together that you can not # tell them apart. # Now we display the outputs from the three methods together. > display( { dir_plot, euler_plot, improved_euler_plot, RK_plot }, > title=`Three Numerical Methods, h=0.1` ); # Homework. Take the differential equation y''=(1 + y^2)/4 and apply # the three numerical methods for y(0)=1 to approximate the solution up # to t=0.5. Use step size h=0.1. For the Euler method, also use # step size h=0.05. Compare the numerical and output of the three # methods. # Also compare the solution with the exact equation. # (To solve use separation of variables. The solution will involve an # arctan and eventually # y(t) = tan((t+evalf(Pi))/4). ) Hint: you will have to redefine the # function f, the ODE, the initial condition every place it is used, # and the formula for the exact solution. # Note above in the construction of the exact solution for the first # equation C is defined in terms of the initial condition. Then # Ex[n] is the exact solution at time T[n] in terms of C. # This line which defines C needs to be removed because for Y[0]=1.0 # the denominator in the definition will be zero. Merely define E[n] # using the formula using tan given above. # Now compare the three methods graphically for this new equation. # Note that if the solutions are too close together only one graph will # be displayed. # Problem 1. Briefly describe what you discovered about the three # different numerical methods for this # example and how they compared with the exact solution.