restart; with( plots ): with( DEtools ): Numerical Methods, System of Equations The worksheet applies the three numerical methods, Euler, Improved Euler, and Runge-Kutta, to a differential equation. The function F(x,y) = (F1(x,y),F2(x,y)) is the right-hand side of the differential equation F1 := (x,y) -> y; F2 := (x,y) -> -4*x; The stepsize h := 0.1; The initial condition T := 0.0: X := 0.0: Y := 1.0: printf(`\t t=%1.1f \t x=%1.4f\t y=%1.4f`,T, X, Y); The vector at the initial point (x0, y0) is (F1( X, Y),F2(X,Y)); the approximate value of the solution at the first point (t=0.1) is approximated to be T := T + h: X := X + h*F1(X,Y): Y := Y + h*F2(X,Y ): printf(`\t t=%1.1f \t x=%1.4f\t y=%1.4f`,T, X, Y); Repeating the same steps four more times yields to give the first five steps: T := T + h: X := X + h*F1(X,Y): Y := Y + h*F2(X,Y): printf(`\t t=%1.1f \t x=%1.4f \t y=%1.4f`,T, X,Y); T := T + h: X := X + h*F1(X,Y): Y := Y + h*F2(X,Y): printf(`\t t=%1.1f \t x=%1.4f \t y=%1.4f`,T, X,Y); T := T + h: X := X + h*F1(X,Y): Y := Y + h*F2(X,Y): printf(`\t t=%1.1f \t x=%1.4f \t y=%1.4f`,T, X,Y); T := T + h: X := X + h*F1(X,Y): Y := Y + h*F2(X,Y): printf(`\t t=%1.1f \t x=%1.4f \t y=%1.4f`,T, X,Y); Combining the solutions into a table and calculating N terms up to t=N*h Step size and number of iterates h := 0.1; N := 10; Initial conditions of time, x, and y: (EX and EY are the x and y coordinated for the Euler method.) T := 0.0; EX := 0.0; EY := 1.0; printf(` t_n \t x_n \t y_n \t k_(1,n)\t k_(2,n)\n`); for n from 0 to N do K_1[n] := F1(EX[n],EY[n]): K_2[n] := F2(EX[n],EY[n]): printf(` %01.2f \t%1.5f \t%01.5f \t%1.5f\t%1.5f\n`,T[n],EX[n],EY[n],K_1[n],K_2[n]); EX[n+1] := EX[n] + h * K_1[n]: EY[n+1] := EY[n] + h * K_2[n]: T[n+1] := T[n] +h: od: To understand what has just been computed, let's look at some plots. First, create the phase portrait and direction field for this problem. SysDE := [diff(x(t),t) = y(t), diff(y(t),t) = -4*x(t)]; sys_plot := DEplot(SysDE, [x(t),y(t)], t=0..1, [[x(0)=EX,y(0)=EY]], linecolour=CYAN, x=-1..1, y=-1..1, stepsize=0.1, method=classical[rk4],scene=[x,y]): The points computed in Euler's Method can also be plotted euler_pts := [[EX,EY], [EX,EY], [EX,EY],[EX,EY],[EX,EY], [EX,EY], [EX,EY], [EX,EY], [EX,EY], [EX,EY], [EX,EY]]; Now display both the plot of the true solution and the approximation using Euler's method euler_plot := plot( euler_pts, color=MAGENTA, thickness=2 ): display( { euler_plot, sys_plot }, title=`Euler's Method, h=0.1` ); Observe how the computed solution closely follows the arrows in the direction field. Also note that the numerical solution increases in size more than the true solution. We now implement the improved Euler method for the same differential equation. IX and IY are the x and y coordinates for the Improved Euler method. h := 0.1: N := 10: T := 0.0: IX := 0.0: IY := 1.0: printf(` t_n \t x_n \t y_n \t k1_(1,n)\t k1_(2,n) \t k2_(1,n) \t k2_(2,n)\n`); for n from 0 to N do K1_1[n] := F1(IX[n],IY[n]): K1_2[n] := F2(IX[n],IY[n]): Z_1[n] := IX[n] + h*K1_1[n]: Z_2[n] := IY[n] + h*K1_2[n]: K2_1[n] := F1(Z_1[n],Z_2[n]): K2_2[n] := F2(Z_1[n],Z_2[n]): printf(` %01.2f \t%1.5f \t%01.5f \t%1.5f\t%1.5f \t%01.5f \t%01.5f\n`,T[n],IX[n],IY[n],K1_1[n],K1_2[n],K2_1[n],K2_2[n]); IX[n+1] := IX[n] + 0.5*h * (K1_1[n]+K2_1[n]): IY[n+1] := IY[n] + 0.5*h * (K1_2[n]+K2_2[n]): T[n+1] := T[n] +h: od: Check the above steps to convice yourself that this is the improved Euler method. Plot of the solution using improved Euler in addition to the phase portrait created by maple. improved_euler_pts := [[IE,IY], [IX,IY], [IX,IY], [IX,IY], [IX,IY], [IX,IY], [IX,IY], [IX,IY], [IX,IY], [IX,IY], [IX,IY]]: improved_euler_plot := plot(improved_euler_pts, color=GREEN, thickness=2 ): display( { sys_plot, improved_euler_plot }, title=`Improved Euler's Method, h=0.1` ); Runge-Kutta for the same differential equation. RKX and RKY are the x and y coordinates for the Runge-Kutta method. h := 0.1: h2 := 0.5*h: N := 10: T := 0.0: RKX := 0.0: RKY := 1.0: printf(` t_n \t x_n \t y_n \t k1_(1,n)\t k1_(2,n) \t k2_(1,n) \t k2_(2,n)\t k3_(1,n) \t k3_(2,n) \t k4_(1,n) \t k4_(2,n)\n`); for n from 0 to N do K1_1[n] := F1(RKX[n],RKY[n]): K1_2[n] := F2(RKX[n],RKY[n]): Z1_1[n] := RKX[n] + h2*K1_1[n]: Z1_2[n] := RKY[n] + h2*K1_2[n]: K2_1[n] := F1(Z1_1[n],Z1_2[n]): K2_2[n] := F2(Z1_1[n],Z1_2[n]): Z2_1[n] := RKX[n] + h2*K2_1[n]: Z2_2[n] := RKY[n] + h2*K2_2[n]: K3_1[n] := F1(Z2_1[n],Z2_2[n]): K3_2[n] := F2(Z2_1[n],Z2_2[n]): Z3_1[n] := RKX[n] + h*K3_1[n]: Z3_2[n] := RKY[n] + h*K3_2[n]: K4_1[n] := F1(Z3_1[n],Z3_2[n]): K4_2[n] := F2(Z3_1[n],Z3_2[n]): printf(` %01.2f \t%1.5f \t%01.5f \t%1.5f\t%1.5f \t %01.5f \t%1.5f \t %1.5f \t%1.5f \t %1.5f \t%01.5f\n`,T[n],RKX[n],RKY[n],K1_1[n],K1_2[n],K2_1[n],K2_2[n],K3_1[n],K3_2[n], K4_1[n],K4_2[n]); RKX[n+1] := RKX[n] + h*(K1_1[n]+2*K2_1[n]+2*K3_1[n]+K4_1[n])/6: RKY[n+1] := RKY[n] + h*(K1_2[n]+2*K2_2[n]+2*K3_2[n]+K4_2[n])/6: T[n+1] := T[n] +h: od: Check the above steps to convice yourself that this is the Runge-Kutta method Compare the solutions. The exact solution has y(t) = cos(2*t). printf(` t \t Euler \t Improved \t Runge-Kutta \t Exact Solution\n`); for n from 0 to N do printf(`%1.2f \t %1.5f \t %1.5f \t %1.5f \t %1.5f\n`, T[n], EY[n], IY[n], RKY[n], cos(2*T[n])); od: Plot the solution using Runge-Kutta and the solution from DEplot: RK_pts := [ [RKX,RKY], [RKX,RKY], [RKX,RKY], [RKX,RKY], [RKX,RKY], [RKX,RKY], [RKX,RKY], [RKX,RKY], [RKX,RKY], [RKX,RKY], [RKX,RKY] ]: RK_plot := plot(RK_pts, colour=BLUE, thickness=2 ): display( {sys_plot, RK_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( { sys_plot, euler_plot, improved_euler_plot, RK_plot }, title=`Three Numerical Methods, h=0.1` ); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2OVEhRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSp1bmRlcmxpbmVHRjcvJSpzdWJzY3JpcHRHRjcvJSxzdXBlcnNjcmlwdEdGNy8lK2ZvcmVncm91bmRHUSpbMjU1LDAsMF1GJy8lK2JhY2tncm91bmRHUS5bMjU1LDI1NSwyNTVdRicvJSdvcGFxdWVHRjcvJStleGVjdXRhYmxlR0Y6LyUpcmVhZG9ubHlHRjcvJSljb21wb3NlZEdGNy8lKmNvbnZlcnRlZEdGNy8lK2ltc2VsZWN0ZWRHRjcvJSxwbGFjZWhvbGRlckdGNy8lMGZvbnRfc3R5bGVfbmFtZUdRKTJEfklucHV0RicvJSptYXRoY29sb3JHRkMvJS9tYXRoYmFja2dyb3VuZEdGRi8lK2ZvbnRmYW1pbHlHRjEvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLyUpbWF0aHNpemVHRjQ=