Math 321 Differential Equations Class04 Euler's Method
The differential equation NiQtSSVtcm93RzYjL0krbW9kdWxlbmFtZUc2IkksVHlwZXNldHRpbmdHSShfc3lzbGliR0YoNiUtSSNtaUdGJTY5USFGKC8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYoLyUlc2l6ZUdRIzEwRigvJSVib2xkR1EmZmFsc2VGKC8lJ2l0YWxpY0dRJXRydWVGKC8lKnVuZGVybGluZUdGOC8lKnN1YnNjcmlwdEdGOC8lLHN1cGVyc2NyaXB0R0Y4LyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GKC8lK2JhY2tncm91bmRHRkQvJSdvcGFxdWVHRjgvJStleGVjdXRhYmxlR0Y4LyUpcmVhZG9ubHlHRjgvJSljb21wb3NlZEdGOC8lKmNvbnZlcnRlZEdGOC8lK2ltc2VsZWN0ZWRHRjgvJSxwbGFjZWhvbGRlckdGOC8lMGZvbnRfc3R5bGVfbmFtZUdRKjJEfk1hdGhfMUYoLyUqbWF0aGNvbG9yR0ZELyUvbWF0aGJhY2tncm91bmRHRkQvJStmb250ZmFtaWx5R0YyLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGKC8lKW1hdGhzaXplR0Y1LUkmbWZyYWNHRiU2Ki1GJDYjLUYtNjlRI2R5RihGMEYzRjZGOUY8Rj5GQEZCRkVGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbkZobkZbby1GJDYjLUYtNjlRI2R0RihGMEYzRjZGOUY8Rj5GQEZCRkVGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbkZobkZbby8lLmxpbmV0aGlja25lc3NHUSIxRigvJStkZW5vbWFsaWduR1EnY2VudGVyRigvJSludW1hbGlnbkdGX3AvJSliZXZlbGxlZEdGOC8lK2ZvcmVncm91bmRHRkQvJStiYWNrZ3JvdW5kR0ZERiw3IzYjKiZJI2R5R0YoIiIiSSNkdEdGKCEiIg== = NiQtSSVtcm93RzYjL0krbW9kdWxlbmFtZUc2IkksVHlwZXNldHRpbmdHSShfc3lzbGliR0YoNiUtSSNtaUdGJTY5USFGKC8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYoLyUlc2l6ZUdRIzEwRigvJSVib2xkR1EmZmFsc2VGKC8lJ2l0YWxpY0dRJXRydWVGKC8lKnVuZGVybGluZUdGOC8lKnN1YnNjcmlwdEdGOC8lLHN1cGVyc2NyaXB0R0Y4LyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GKC8lK2JhY2tncm91bmRHRkQvJSdvcGFxdWVHRjgvJStleGVjdXRhYmxlR0Y4LyUpcmVhZG9ubHlHRjgvJSljb21wb3NlZEdGOC8lKmNvbnZlcnRlZEdGOC8lK2ltc2VsZWN0ZWRHRjgvJSxwbGFjZWhvbGRlckdGOC8lMGZvbnRfc3R5bGVfbmFtZUdRKjJEfk1hdGhfMkYoLyUqbWF0aGNvbG9yR0ZELyUvbWF0aGJhY2tncm91bmRHRkQvJStmb250ZmFtaWx5R0YyLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGKC8lKW1hdGhzaXplR0Y1LUYkNiYtRi02OVEiZkYoRjBGM0Y2RjlGPEY+RkBGQkZFRkdGSUZLRk1GT0ZRRlNGVUZYRlpGZm5GaG5GW28tSSNtb0dGJTYzUTAmQXBwbHlGdW5jdGlvbjtGKC8lJWZvcm1HUSZpbmZpeEYoLyUmZmVuY2VHRjgvJSpzZXBhcmF0b3JHRjgvJSdsc3BhY2VHUSQwZW1GKC8lJ3JzcGFjZUdGX3AvJSlzdHJldGNoeUdGOC8lKnN5bW1ldHJpY0dGOC8lKG1heHNpemVHUSlpbmZpbml0eUYoLyUobWluc2l6ZUdRIjFGKC8lKGxhcmdlb3BHRjgvJS5tb3ZhYmxlbGltaXRzR0Y4LyUnYWNjZW50R0Y4LyUwZm9udF9zdHlsZV9uYW1lR0ZXLyUlc2l6ZUdGNS8lK2ZvcmVncm91bmRHRkQvJStiYWNrZ3JvdW5kR0ZELUYkNiUtRmNvNjNRIihGKC9GZ29RJ3ByZWZpeEYoL0Zqb0Y7RltwL0ZecFEudGhpbm1hdGhzcGFjZUYoL0ZhcEZjci9GY3BGO0ZkcEZmcEZpcEZccUZecUZgcUZicUZkcUZmcUZocS1GJDYlLUYtNjlRInRGKEYwRjNGNkY5RjxGPkZARkJGRUZHRklGS0ZNRk9GUUZTRlVGWEZaRmZuRmhuRltvLUZjbzYzUSIsRihGZm9GaW8vRlxwRjtGXXAvRmFwUTN2ZXJ5dGhpY2ttYXRoc3BhY2VGKEZicEZkcEZmcEZpcEZccUZecUZgcUZicUZkcUZmcUZocS1GLTY5USJ5RihGMEYzRjZGOUY8Rj5GQEZCRkVGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbkZobkZbby1GY282M1EiKUYoL0Znb1EocG9zdGZpeEYoRmFyRltwRmJyL0ZhcFEydmVyeXRoaW5tYXRoc3BhY2VGKEZlckZkcEZmcEZpcEZccUZecUZgcUZicUZkcUZmcUZocUYsRiw3IzYjLUkiZkdGKDYkSSJ0R0YoSSJ5R0Yo = NiQtSSVtcm93RzYjL0krbW9kdWxlbmFtZUc2IkksVHlwZXNldHRpbmdHSShfc3lzbGliR0YoNiUtSSNtaUdGJTY5USFGKC8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYoLyUlc2l6ZUdRIzEwRigvJSVib2xkR1EmZmFsc2VGKC8lJ2l0YWxpY0dRJXRydWVGKC8lKnVuZGVybGluZUdGOC8lKnN1YnNjcmlwdEdGOC8lLHN1cGVyc2NyaXB0R0Y4LyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GKC8lK2JhY2tncm91bmRHRkQvJSdvcGFxdWVHRjgvJStleGVjdXRhYmxlR0Y4LyUpcmVhZG9ubHlHRjgvJSljb21wb3NlZEdGOC8lKmNvbnZlcnRlZEdGOC8lK2ltc2VsZWN0ZWRHRjgvJSxwbGFjZWhvbGRlckdGOC8lMGZvbnRfc3R5bGVfbmFtZUdRKjJEfk1hdGhfM0YoLyUqbWF0aGNvbG9yR0ZELyUvbWF0aGJhY2tncm91bmRHRkQvJStmb250ZmFtaWx5R0YyLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGKC8lKW1hdGhzaXplR0Y1LUYkNiUtSSNtb0dGJTYzUSomdW1pbnVzMDtGKC8lJWZvcm1HUSdwcmVmaXhGKC8lJmZlbmNlR0Y4LyUqc2VwYXJhdG9yR0Y4LyUnbHNwYWNlR1EkMGVtRigvJSdyc3BhY2VHRlxwLyUpc3RyZXRjaHlHRjgvJSpzeW1tZXRyaWNHRjgvJShtYXhzaXplR1EpaW5maW5pdHlGKC8lKG1pbnNpemVHUSIxRigvJShsYXJnZW9wR0Y4LyUubW92YWJsZWxpbWl0c0dGOC8lJ2FjY2VudEdGOC8lMGZvbnRfc3R5bGVfbmFtZUdGVy8lJXNpemVHRjUvJStmb3JlZ3JvdW5kR0ZELyUrYmFja2dyb3VuZEdGRC1GJDYoLUkjbW5HRiU2OVEiMkYoRjBGM0Y2L0Y6RjhGPEY+RkBGQkZFRkdGSUZLRk1GT0ZRRlNGVUZYRlpGZm4vRmluUSdub3JtYWxGKEZbby1GYG82M1ExJkludmlzaWJsZVRpbWVzO0YoL0Zkb1EmaW5maXhGKEZmb0Zob0Zqb0ZdcEZfcEZhcEZjcEZmcEZpcEZbcUZdcUZfcUZhcUZjcUZlcS1GLTY5USJ0RihGMEYzRjZGOUY8Rj5GQEZCRkVGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbkZobkZbb0Zgci1JJW1zdXBHRiU2JS1GLTY5USJ5RihGMEYzRjZGOUY8Rj5GQEZCRkVGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbkZobkZbb0ZpcS8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRihGLEYsRiw3IzYjLCQqKCIiIyIiIkkidEdGKEZmcylJInlHRihGZXNGZnMhIiI=.
f := (t, y) -> -2*t*y^2;
The initial condition NiQtSSVtcm93RzYjL0krbW9kdWxlbmFtZUc2IkksVHlwZXNldHRpbmdHSShfc3lzbGliR0YoNiUtSSNtaUdGJTY5USFGKC8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYoLyUlc2l6ZUdRIzEwRigvJSVib2xkR1EmZmFsc2VGKC8lJ2l0YWxpY0dRJXRydWVGKC8lKnVuZGVybGluZUdGOC8lKnN1YnNjcmlwdEdGOC8lLHN1cGVyc2NyaXB0R0Y4LyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GKC8lK2JhY2tncm91bmRHRkQvJSdvcGFxdWVHRjgvJStleGVjdXRhYmxlR0Y4LyUpcmVhZG9ubHlHRjgvJSljb21wb3NlZEdGOC8lKmNvbnZlcnRlZEdGOC8lK2ltc2VsZWN0ZWRHRjgvJSxwbGFjZWhvbGRlckdGOC8lMGZvbnRfc3R5bGVfbmFtZUdRKjJEfk1hdGhfNEYoLyUqbWF0aGNvbG9yR0ZELyUvbWF0aGJhY2tncm91bmRHRkQvJStmb250ZmFtaWx5R0YyLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGKC8lKW1hdGhzaXplR0Y1LUYkNiYtRi02OVEieUYoRjBGM0Y2RjlGPEY+RkBGQkZFRkdGSUZLRk1GT0ZRRlNGVUZYRlpGZm5GaG5GW28tSSNtb0dGJTYzUTAmQXBwbHlGdW5jdGlvbjtGKC8lJWZvcm1HUSZpbmZpeEYoLyUmZmVuY2VHRjgvJSpzZXBhcmF0b3JHRjgvJSdsc3BhY2VHUSQwZW1GKC8lJ3JzcGFjZUdGX3AvJSlzdHJldGNoeUdGOC8lKnN5bW1ldHJpY0dGOC8lKG1heHNpemVHUSlpbmZpbml0eUYoLyUobWluc2l6ZUdRIjFGKC8lKGxhcmdlb3BHRjgvJS5tb3ZhYmxlbGltaXRzR0Y4LyUnYWNjZW50R0Y4LyUwZm9udF9zdHlsZV9uYW1lR0ZXLyUlc2l6ZUdGNS8lK2ZvcmVncm91bmRHRkQvJStiYWNrZ3JvdW5kR0ZELUYkNiUtRmNvNjNRIihGKC9GZ29RJ3ByZWZpeEYoL0Zqb0Y7RltwL0ZecFEudGhpbm1hdGhzcGFjZUYoL0ZhcEZjci9GY3BGO0ZkcEZmcEZpcEZccUZecUZgcUZicUZkcUZmcUZocS1GJDYjLUkjbW5HRiU2OVEiMEYoRjBGM0Y2L0Y6RjhGPEY+RkBGQkZFRkdGSUZLRk1GT0ZRRlNGVUZYRlpGZm4vRmluUSdub3JtYWxGKEZbby1GY282M1EiKUYoL0Znb1EocG9zdGZpeEYoRmFyRltwRmJyL0ZhcFEydmVyeXRoaW5tYXRoc3BhY2VGKEZlckZkcEZmcEZpcEZccUZecUZgcUZicUZkcUZmcUZocUYsRiw3IzYjLUkieUdGKDYjIiIh = NiQtSSVtcm93RzYjL0krbW9kdWxlbmFtZUc2IkksVHlwZXNldHRpbmdHSShfc3lzbGliR0YoNiUtSSNtaUdGJTY5USFGKC8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYoLyUlc2l6ZUdRIzEwRigvJSVib2xkR1EmZmFsc2VGKC8lJ2l0YWxpY0dRJXRydWVGKC8lKnVuZGVybGluZUdGOC8lKnN1YnNjcmlwdEdGOC8lLHN1cGVyc2NyaXB0R0Y4LyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GKC8lK2JhY2tncm91bmRHRkQvJSdvcGFxdWVHRjgvJStleGVjdXRhYmxlR0Y4LyUpcmVhZG9ubHlHRjgvJSljb21wb3NlZEdGOC8lKmNvbnZlcnRlZEdGOC8lK2ltc2VsZWN0ZWRHRjgvJSxwbGFjZWhvbGRlckdGOC8lMGZvbnRfc3R5bGVfbmFtZUdRKjJEfk1hdGhfNUYoLyUqbWF0aGNvbG9yR0ZELyUvbWF0aGJhY2tncm91bmRHRkQvJStmb250ZmFtaWx5R0YyLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGKC8lKW1hdGhzaXplR0Y1LUklbXN1YkdGJTYmLUYtNjlRInRGKEYwRjNGNkY5RjxGPkZARkJGRUZHRklGS0ZNRk9GUUZTRlVGWEZaRmZuRmhuRltvLUYkNiMtSSNtbkdGJTY5USIwRihGMEYzRjYvRjpGOEY8Rj5GQEZCRkVGR0ZJRktGTUZPRlFGU0ZVRlhGWkZmbi9GaW5RJ25vcm1hbEYoRltvLyUvc3Vic2NyaXB0c2hpZnRHUSIwRigvJSxwbGFjZWhvbGRlckdGOEYsNyM2IyZJInRHRig2IyIiIQ==, the time increment, and the number of increments.
t[0] := 0; y[0] := 1; Delta_t := 0.3; max_k := 10;
First step of Euler's method.
t[1] := t[0] + Delta_t;
y[1] := y[0] + f(t[0],y[0])*Delta_t;
LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2OVEhRicvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSp1bmRlcmxpbmVHRjcvJSpzdWJzY3JpcHRHRjcvJSxzdXBlcnNjcmlwdEdGNy8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJStiYWNrZ3JvdW5kR1EuWzI1NSwyNTUsMjU1XUYnLyUnb3BhcXVlR0Y3LyUrZXhlY3V0YWJsZUdGOi8lKXJlYWRvbmx5R0Y3LyUpY29tcG9zZWRHRjcvJSpjb252ZXJ0ZWRHRjcvJStpbXNlbGVjdGVkR0Y3LyUscGxhY2Vob2xkZXJHRjcvJTBmb250X3N0eWxlX25hbWVHUSkyRH5JbnB1dEYnLyUqbWF0aGNvbG9yR0ZDLyUvbWF0aGJhY2tncm91bmRHRkYvJStmb250ZmFtaWx5R0YxLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy8lKW1hdGhzaXplR0Y0t[2] := t[1] + Delta_t;
y[2] := y[1] + f(t[1],y[1])*Delta_t;
Many steps of Euler's method.
for k from 0 to (max_k - 1) do
t[k+1] := t[k] + delta_t;
y[k+1] := y[k] + f(t[k],y[k])*delta_t;
od;
points := [seq([t[k],y[k]],k=0..max_k)];
plot(points, style=point);
Library needed for the display function that will superimpose two or more types of plots.
with(plots):
Superimpose Euler approximate solution on known analytic solution.
points := [seq([t[k],y[k]],k=0..max_k)]:
p1 := plot(points,style=point):
soln := t -> 1/(1+t^2):
p2 := plot(soln(t), t=t[0]..t[max_k]):
display([p1,p2]);
View and calculate the errors.
errors := seq(abs(y[k]-soln(t[k])),k=0..max_k):
points := [seq([t[k],[errors][k+1]],k=0..max_k)]:
plot(points,style=point);
max(errors);
errors[max_k];
Putting it all together.
f := (t, y) -> -2*t*y^2; t[0] := 0; y[0] := 1;
delta_t := 0.001; max_k := 3000;
for k from 0 to max_k-1 do
t[k+1] := t[k] + delta_t;
y[k+1] := y[k] + f(t[k],y[k])*delta_t;
od:
points := [seq([t[k],y[k]],k=0..max_k)]:
p1 := plot(points,style=point):
soln := t -> 1/(1+t^2):
p2 := plot(soln(t), t=t[0]..t[max_k]):
display([p1,p2]);
errors := seq(abs(y[k]-soln(t[k])),k=0..max_k):
points := [seq([t[k],[errors][k+1]],k=0..max_k)]:
plot(points,style=point);
max(errors);
errors[max_k];