restart;with(plots):with(DEtools):Numerical Methods, Scalar EquationThe worksheet applies the three numerical methods, Euler, Improved Euler, and Runge-Kutta, to a scalar differential equation.The function, LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2OVEiZkYnLyUnZmFtaWx5R1EwVGltZXN+TmV3flJvbWFuRicvJSVzaXplR1EjMTJGJy8lJWJvbGRHUSZmYWxzZUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUqdW5kZXJsaW5lR0Y3LyUqc3Vic2NyaXB0R0Y3LyUsc3VwZXJzY3JpcHRHRjcvJStmb3JlZ3JvdW5kR1EoWzAsMCwwXUYnLyUrYmFja2dyb3VuZEdRLlsyNTUsMjU1LDI1NV1GJy8lJ29wYXF1ZUdGNy8lK2V4ZWN1dGFibGVHRjcvJSlyZWFkb25seUdGNy8lKWNvbXBvc2VkR0Y3LyUqY29udmVydGVkR0Y3LyUraW1zZWxlY3RlZEdGNy8lLHBsYWNlaG9sZGVyR0Y3LyUwZm9udF9zdHlsZV9uYW1lR1ErMkR+Q29tbWVudEYnLyUqbWF0aGNvbG9yR0ZDLyUvbWF0aGJhY2tncm91bmRHRkYvJStmb250ZmFtaWx5R0YxLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy8lKW1hdGhzaXplR0Y0, that is the right-hand side of the differential equationf := (t,y) -> y*(1-y^2);The stepsizeh := 0.1;The initial conditionT[0] := 0.0; Y[0] := 2.0;The slope at the initial point (x=1, y=4.0) isf( T[0], Y[0] );The approximate value of the solution at the first point (t=0.1) is approximated to beT[1] := T[0] + h:
Y[1] := Y[0] + h*f( T[0], Y[0] ):
printf(`\t t=%1.1f \t y=%1.4f\t k=%1.4f`,T[1], Y[1], f(T[1],Y[1]));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 t=%1.1f \t y=%1.4f \t k=%1.4f`,T[2], Y[2],f(T[2],Y[2]));T[3] := T[2] + h:
Y[3] := Y[2] + h*f( T[2], Y[2]):
printf(`\t 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 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 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.0h := 0.1:
N := 10:
T[0] := 0.0:
Y[0] := 2.0:
printf(` \t Euler method\n`);
printf(` t_n \t y_n \t k_n\n`);
for n from 0 to N do
K[n] := f(T[n],Y[n]):
printf(` % 1.2f \t%1.5f \t%01.5f \n`,T[n], Y[n],K[n]);
T[n+1] := T[n] +h:
Y[n+1] := Y[n] + h * 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 plottedeuler_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]]];Now displace both the phase portrait created by maple and the approximation using Euler's methodeuler_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.We now implement the improved Euler method for the same differential equationh := 0.1:
N := 10:
T[0] := 0.0:
IY[0] := 2.0:
printf(` \t Improved Euler method \n`);
printf(` t_n \t y_n \t k1_n \t k2_n \n`);
for n from 0 to N do
T[n+1] := T[n] +h:
K1[n] := f(T[n],IY[n]):
Z[n] := IY[n] + h*K1[n]:
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]);
IY[n+1] := IY[n] + 0.5*h*(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 in addition to the two previous plots. 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:
print(`Runge-Kutta method`);
printf(` t_n \t y_n \t k1_n \t z1_n \t k2_n \t\t z2_n \t k3_n \t z3_n \t k4_n \n`);
for n from 0 to N do
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 \t%01.5f \t%01.5f \t%01.5f \n`,
T[n], RKY[n],K1[n],Z1[n],K2[n],Z2[n],K3[n],Z3[n],K4[n]);
RKY[n+1] := RKY[n] + h*(K1[n]+2*K2[n]+2*K3[n]+K4[n])/6.0:
od:Check the above steps to convice yourself that this is the Runge-Kutta methodThe outputs of the three methods can be compared with the exact solution.
For the above example the exact solution is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2OVEhRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSp1bmRlcmxpbmVHRjcvJSpzdWJzY3JpcHRHRjcvJSxzdXBlcnNjcmlwdEdGNy8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJStiYWNrZ3JvdW5kR1EuWzI1NSwyNTUsMjU1XUYnLyUnb3BhcXVlR0Y3LyUrZXhlY3V0YWJsZUdGNy8lKXJlYWRvbmx5R0Y3LyUpY29tcG9zZWRHRjcvJSpjb252ZXJ0ZWRHRjcvJStpbXNlbGVjdGVkR0Y3LyUscGxhY2Vob2xkZXJHRjcvJTBmb250X3N0eWxlX25hbWVHUSsyRH5Db21tZW50RicvJSptYXRoY29sb3JHRkMvJS9tYXRoYmFja2dyb3VuZEdGRi8lK2ZvbnRmYW1pbHlHRjEvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLyUpbWF0aHNpemVHRjQtSSZtc3FydEdGJDYlLUkmbWZyYWNHRiQ2Ki1GIzYmLUYsNjlRIkNGJ0YvRjJGNUY4RjtGPUY/RkFGREZHRklGS0ZNRk9GUUZTRlVGWEZaRmZuRmhuRltvLUkjbW9HRiQ2M1ExJkludmlzaWJsZVRpbWVzO0YnLyUlZm9ybUdRJmluZml4RicvJSZmZW5jZUdGNy8lKnNlcGFyYXRvckdGNy8lJ2xzcGFjZUdRJDBlbUYnLyUncnNwYWNlR0ZlcC8lKXN0cmV0Y2h5R0Y3LyUqc3ltbWV0cmljR0Y3LyUobWF4c2l6ZUdRKWluZmluaXR5RicvJShtaW5zaXplR1EiMUYnLyUobGFyZ2VvcEdGNy8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJTBmb250X3N0eWxlX25hbWVHRlcvJSVzaXplR0Y0LyUrZm9yZWdyb3VuZEdGQy8lK2JhY2tncm91bmRHRkYtSSVtc3VwR0YkNiUtRmlvNjNRImVGJy9GXXBGLkZfcEZhcEZjcEZmcEZocEZqcEZccUZfcUZicUZkcUZmcUZocUZqcUZcckZeci1GIzYlLUZpbzYzUSIoRicvRl1wUSdwcmVmaXhGJy9GYHBGOkZhcC9GZHBRLnRoaW5tYXRoc3BhY2VGJy9GZ3BGYHMvRmlwRjpGanBGXHFGX3FGYnFGZHFGZnFGaHFGanFGXHJGXnItRiM2Ji1JI21uR0YkNjlRIjJGJ0YvRjJGNS9GOUY3RjtGPUY/RkFGREZHRklGS0ZNRk9GUUZTRlVGWEZaRmZuL0ZpblEnbm9ybWFsRidGW29GaG8tSSVtc3ViR0YkNiYtRiw2OVEiVEYnRi9GMkY1RjhGO0Y9Rj9GQUZERkdGSUZLRk1GT0ZRRlNGVUZYRlpGZm5GaG5GW28tRiM2Iy1GLDY5USJuRidGL0YyRjVGOEY7Rj1GP0ZBRkRGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbkZobkZbby8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnLyUscGxhY2Vob2xkZXJHRjdGKy1GaW82M1EiKUYnL0ZdcFEocG9zdGZpeEYnRl5zRmFwRl9zL0ZncFEydmVyeXRoaW5tYXRoc3BhY2VGJ0Zic0ZqcEZccUZfcUZicUZkcUZmcUZocUZqcUZcckZeci8lMXN1cGVyc2NyaXB0c2hpZnRHRml0RistRiM2JUYrLUYjNiZGZW9GaG8tRmFyNiVGY3ItRiM2JUZpci1GIzYmRitGY3MtRmlvNjNRKCZtaW51cztGJ0ZccEZfcEZhcC9GZHBRMG1lZGl1bW1hdGhzcGFjZUYnL0ZncEZjdkZocEZqcEZccUZfcUZicUZkcUZmcUZocUZqcUZcckZeci1GZnM2OVEkMS4wRidGL0YyRjVGaXNGO0Y9Rj9GQUZERkdGSUZLRk1GT0ZRRlNGVUZYRlpGZm5GanNGW29GXHVGY3VGK0YrLyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zddy8lKWJldmVsbGVkR0Y3RlxyRl5yRlxyRl5yRis=where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2OVEhRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSp1bmRlcmxpbmVHRjcvJSpzdWJzY3JpcHRHRjcvJSxzdXBlcnNjcmlwdEdGNy8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJStiYWNrZ3JvdW5kR1EuWzI1NSwyNTUsMjU1XUYnLyUnb3BhcXVlR0Y3LyUrZXhlY3V0YWJsZUdGNy8lKXJlYWRvbmx5R0Y3LyUpY29tcG9zZWRHRjcvJSpjb252ZXJ0ZWRHRjcvJStpbXNlbGVjdGVkR0Y3LyUscGxhY2Vob2xkZXJHRjcvJTBmb250X3N0eWxlX25hbWVHUSsyRH5Db21tZW50RicvJSptYXRoY29sb3JHRkMvJS9tYXRoYmFja2dyb3VuZEdGRi8lK2ZvbnRmYW1pbHlHRjEvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLyUpbWF0aHNpemVHRjQtRiM2Ji1GLDY5USJDRidGL0YyRjVGOEY7Rj1GP0ZBRkRGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbkZobkZbby1JI21vR0YkNjNRIj1GJy8lJWZvcm1HUSZpbmZpeEYnLyUmZmVuY2VHRjcvJSpzZXBhcmF0b3JHRjcvJSdsc3BhY2VHUS90aGlja21hdGhzcGFjZUYnLyUncnNwYWNlR0ZfcC8lKXN0cmV0Y2h5R0Y3LyUqc3ltbWV0cmljR0Y3LyUobWF4c2l6ZUdRKWluZmluaXR5RicvJShtaW5zaXplR1EiMUYnLyUobGFyZ2VvcEdGNy8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJTBmb250X3N0eWxlX25hbWVHRlcvJSVzaXplR0Y0LyUrZm9yZWdyb3VuZEdGQy8lK2JhY2tncm91bmRHRkYtSSZtZnJhY0dGJDYqLUYjNiVGKy1JKG1zdWJzdXBHRiQ2Jy1GLDY5USJZRidGL0YyRjVGOEY7Rj1GP0ZBRkRGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbkZobkZbby1GIzYjLUkjbW5HRiQ2OVEiMEYnRi9GMkY1L0Y5RjdGO0Y9Rj9GQUZERkdGSUZLRk1GT0ZRRlNGVUZYRlpGZm4vRmluUSdub3JtYWxGJ0Zbby1GaHI2OVEiMkYnRi9GMkY1RltzRjtGPUY/RkFGREZHRklGS0ZNRk9GUUZTRlVGWEZaRmZuRlxzRltvLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy8lL3N1YnNjcmlwdHNoaWZ0R0Zjc0YrLUYjNiVGKy1GIzYmRitGXXItRmNvNjNRKCZtaW51cztGJ0Zmb0Zpb0ZbcC9GXnBRMG1lZGl1bW1hdGhzcGFjZUYnL0ZhcEZedEZicEZkcEZmcEZpcEZccUZecUZgcUZicUZkcUZmcUZocS1GaHI2OVEkMS4wRidGL0YyRjVGW3NGO0Y9Rj9GQUZERkdGSUZLRk1GT0ZRRlNGVUZYRlpGZm5GXHNGW29GKy8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGaHQvJSliZXZlbGxlZEdGN0ZmcUZocUYrRis=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` );LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2OVEhRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSp1bmRlcmxpbmVHRjcvJSpzdWJzY3JpcHRHRjcvJSxzdXBlcnNjcmlwdEdGNy8lK2ZvcmVncm91bmRHUSpbMjU1LDAsMF1GJy8lK2JhY2tncm91bmRHUS5bMjU1LDI1NSwyNTVdRicvJSdvcGFxdWVHRjcvJStleGVjdXRhYmxlR0Y6LyUpcmVhZG9ubHlHRjcvJSljb21wb3NlZEdGNy8lKmNvbnZlcnRlZEdGNy8lK2ltc2VsZWN0ZWRHRjcvJSxwbGFjZWhvbGRlckdGNy8lMGZvbnRfc3R5bGVfbmFtZUdRKTJEfklucHV0RicvJSptYXRoY29sb3JHRkMvJS9tYXRoYmFja2dyb3VuZEdGRi8lK2ZvbnRmYW1pbHlHRjEvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLyUpbWF0aHNpemVHRjQ=