restart;Runge-Kutta stencils to solve LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUkmbWZyYWNHRiQ2KC1JI21vR0YkNi5RMCZEaWZmZXJlbnRpYWxEO0YnLyUlc2l6ZUdRIzE4RicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLyUmZmVuY2VHUSZ1bnNldEYnLyUqc2VwYXJhdG9yR0Y6LyUpc3RyZXRjaHlHRjovJSpzeW1tZXRyaWNHRjovJShsYXJnZW9wR0Y6LyUubW92YWJsZWxpbWl0c0dGOi8lJ2FjY2VudEdGOi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRkktRiM2KEYuLUYvNi5RIn5GJ0YyRjUvRjlRJmZhbHNlRicvRjxGUi9GPkZSL0ZARlIvRkJGUi9GREZSL0ZGRlJGR0ZKLUkjbWlHRiQ2JlEieEYnRjIvJSdpdGFsaWNHUSV0cnVlRicvRjZRJ2l0YWxpY0YnRjIvJTBmb250X3N0eWxlX25hbWVHUSdOb3JtYWxGJ0Y1LyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zkby8lKWJldmVsbGVkR0ZSLUZaNiZRInlGJ0YyRmduRmpuLUkobWZlbmNlZEdGJDYlLUYjNiZGWUYyRlxvRjVGMkY1Rk4tRi82LlEiPUYnRjJGNUZRRlNGVEZVRlZGV0ZYL0ZIUSwwLjI3Nzc3NzhlbUYnL0ZLRmVwRk4tRlo2JlEiZkYnRjJGZ25Gam4tRl1wNiUtRiM2KEZZLUYvNi5RIixGJ0YyRjVGUS9GPEZpbkZURlVGVkZXRlhGRy9GS1EsMC4zMzMzMzMzZW1GJ0Zpby1GXXA2JS1GIzYlRllGMkY1RjJGNUYyRjVGMkY1RjJGNQ==The purpose of this worksheet is to introduce and derive the Runge-Kutta stencils to solve the first order equation LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkmbWZyYWNHRiQ2KC1JI21vR0YkNi1RMCZEaWZmZXJlbnRpYWxEO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R0Y3LyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGNy8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZGLUYjNiZGLi1GLzYtUSJ+RidGMi9GNlEmZmFsc2VGJy9GOUZPL0Y7Rk8vRj1GTy9GP0ZPL0ZBRk8vRkNGT0ZERkctSSNtaUdGJDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvRjNRJ2l0YWxpY0YnRjIvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl5vLyUpYmV2ZWxsZWRHRk8tRlc2JVEieUYnRlpGZ24tSShtZmVuY2VkR0YkNiQtRiM2JEZWRjJGMi1GLzYtUSI9RidGMkZORlBGUUZSRlNGVEZVL0ZFUSwwLjI3Nzc3NzhlbUYnL0ZIRl9wLUZXNiVRImZGJ0ZaRmduLUZnbzYkLUYjNidGVi1GLzYtUSIsRidGMkZOL0Y5RmZuRlFGUkZTRlRGVUZEL0ZIUSwwLjMzMzMzMzNlbUYnRmNvRmZvRjJGMi1GVzYjUSFGJ0Yy, and study some of their properties.General form of the stencilsThe generic LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=-stage Runge-Kutta stencil to solveLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkmbWZyYWNHRiQ2KC1JI21vR0YkNi1RMCZEaWZmZXJlbnRpYWxEO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R0Y3LyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGNy8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZGLUYjNidGLi1GLzYtUSJ+RidGMi9GNlEmZmFsc2VGJy9GOUZPL0Y7Rk8vRj1GTy9GP0ZPL0ZBRk8vRkNGT0ZERkctSSNtaUdGJDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvRjNRJ2l0YWxpY0YnLyUwZm9udF9zdHlsZV9uYW1lR1EnTm9ybWFsRidGMi8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGYW8vJSliZXZlbGxlZEdGTy1GVzYlUSJ5RidGWkZnbi1JKG1mZW5jZWRHRiQ2JC1GIzYlRlZGaW5GMkYyLUYvNi1RIj1GJ0YyRk5GUEZRRlJGU0ZURlUvRkVRLDAuMjc3Nzc3OGVtRicvRkhGYnAtRlc2JVEiZkYnRlpGZ24tRmpvNiQtRiM2KEZWLUYvNi1RIixGJ0YyRk4vRjlGZm5GUUZSRlNGVEZVRkQvRkhRLDAuMzMzMzMzM2VtRidGZm9GaW9GaW5GMkYyRmluRjI= involves defining LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= quantities LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn as follows:LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUSI9RicvRjZRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZILyUpc3RyZXRjaHlHRkgvJSpzeW1tZXRyaWNHRkgvJShsYXJnZW9wR0ZILyUubW92YWJsZWxpbWl0c0dGSC8lJ2FjY2VudEdGSC8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRlctRi82JVEiZkYnRjJGNS1JKG1mZW5jZWRHRiQ2JC1GIzYuLUYsNiUtRi82JVEieEYnRjJGNS1GIzYlLUYvNiVRImlGJ0YyRjVGMkY1Rj0tRkE2LVEiK0YnRkRGRkZJRktGTUZPRlFGUy9GVlEsMC4yMjIyMjIyZW1GJy9GWUZqby1GLDYlLUYvNiVRImNGJ0YyRjVGOEY9LUYvNiVRImhGJ0YyRjUtRkE2LVEiLEYnRkRGRi9GSkY0RktGTUZPRlFGUy9GVlEmMC4wZW1GJy9GWVEsMC4zMzMzMzMzZW1GJy1GLDYlLUYvNiVRInlGJ0YyRjVGYW9GPUZmb0ZhcC1JK211bmRlcm92ZXJHRiQ2Jy1GQTYtUSYmU3VtO0YnRkQvRkdRJnVuc2V0RicvRkpGaHEvRkxGNC9GTkZocS9GUEY0L0ZSRjQvRlRGaHFGaHAvRllRLDAuMTY2NjY2N2VtRictRiM2Ji1GLzYlUSJtRidGMkY1RkAtSSNtbkdGJDYkUSIxRidGREZELUYjNiUtRi82JVEiTUYnRjJGNUYyRjVGUy8lLGFjY2VudHVuZGVyR0ZILUYsNiUtRi82JVEnJiM5NDU7RicvRjNGSEZELUYjNidGOkZkcEZjckYyRjVGPS1GLDYlRi4tRiM2JUZjckYyRjVGPUZERkRGRA==In this expression, LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRImlGJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUS0mVGlsZGVUaWxkZTtGJy9GNlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkgvJSlzdHJldGNoeUdGSC8lKnN5bW1ldHJpY0dGSC8lKGxhcmdlb3BHRkgvJS5tb3ZhYmxlbGltaXRzR0ZILyUnYWNjZW50R0ZILyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGV0YuLUkobWZlbmNlZEdGJDYkLUYjNiQtRiw2JS1GLzYlUSJ4RidGMkY1RjhGPUZERkRGRA== represents our numeric approximation to solution of the differential equation at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JJW1zdWJHRiQ2JUYrLUYjNiUtRiw2JVEiaUYnRi9GMkYvRjIvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0Y5 as usual. The coefficients LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiY0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEnJiM5NDU7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYnLUYvNiVRIm5GJy9GM1EldHJ1ZUYnL0Y2USdpdGFsaWNGJy1JI21vR0YkNi1RIixGJ0Y1LyUmZmVuY2VHRjQvJSpzZXBhcmF0b3JHRj4vJSlzdHJldGNoeUdGNC8lKnN5bW1ldHJpY0dGNC8lKGxhcmdlb3BHRjQvJS5tb3ZhYmxlbGltaXRzR0Y0LyUnYWNjZW50R0Y0LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRi82JVEibUYnRj1GP0Y9Rj8vJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0Y1 are constants that we are going to choose shortly, and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEiaEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JJW1zdWJHRiQ2JS1GLDYlUSJ4RidGL0YyLUYjNictRiw2JVEiaUYnRi9GMi1GNjYtUSIrRidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjIyMjIyMjJlbUYnL0ZORmhuLUkjbW5HRiQ2JFEiMUYnRjlGL0YyLyUvc3Vic2NyaXB0c2hpZnRHUSIwRictRjY2LVEoJm1pbnVzO0YnRjlGO0Y+RkBGQkZERkZGSEZnbkZpbi1GUDYlRlItRiM2JUZXRi9GMkZeb0Y5 is the stepsize. Once one calculates LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn given LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkobWZlbmNlZEdGJDYkLUYjNiYtSSVtc3ViR0YkNiUtSSNtaUdGJDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiUtRjQ2JVEiaUYnRjdGOkY3RjovJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIixGJy9GO1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjkvJSlzdHJldGNoeUdGTS8lKnN5bW1ldHJpY0dGTS8lKGxhcmdlb3BHRk0vJS5tb3ZhYmxlbGltaXRzR0ZNLyUnYWNjZW50R0ZNLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjE2JS1GNDYlUSJ5RidGN0Y6Rj1GQkZJRklGSQ==, the value of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUYvNiVRImlGJ0YyRjUtSSNtb0dGJDYtUSIrRicvRjZRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZFLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRlQtSSNtbkdGJDYkUSIxRidGQUYyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1GPjYtUS0mVGlsZGVUaWxkZTtGJ0ZBRkNGRkZIRkpGTEZORlAvRlNRLDAuMjc3Nzc3OGVtRicvRlZGXG9GLi1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiUtRi82JVEieEYnRjJGNUY4RmVuRkFGQS1GLzYjUSFGJ0ZB is given byLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUYvNiVRImlGJ0YyRjUtSSNtb0dGJDYtUSIrRicvRjZRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZFLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRlQtSSNtbkdGJDYkUSIxRidGQUYyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1GPjYtUSI9RidGQUZDRkZGSEZKRkxGTkZQL0ZTUSwwLjI3Nzc3NzhlbUYnL0ZWRlxvLUYsNiVGLi1GIzYlRjpGMkY1RmVuRj0tRi82JVEiaEYnRjJGNS1JK211bmRlcm92ZXJHRiQ2Jy1GPjYtUSYmU3VtO0YnRkEvRkRRJnVuc2V0RicvRkdGXHAvRklGNC9GS0ZccC9GTUY0L0ZPRjQvRlFGXHAvRlNRJjAuMGVtRicvRlZRLDAuMTY2NjY2N2VtRictRiM2Ji1GLzYlUSJuRidGMkY1RmhuRldGQS1GIzYlLUYvNiVRIk1GJ0YyRjVGMkY1RlAvJSxhY2NlbnR1bmRlckdGRS1GLDYlLUYvNiVRImJGJ0YyRjUtRiM2JUZpcEYyRjVGZW4tRiw2JS1GLzYlUSJrRidGMkY1RmhxRmVuLUY+Ni1RIi5GJ0ZBRkNGRkZIRkpGTEZORlBGY3AvRlZGZHBGQQ==In this expression, the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiYkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRictRi82I1EhRicvRjZRJ25vcm1hbEYn are constants. Here is what these equations look like for a particular choice of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= (which you can change):unassign(alpha,b,c):
M := 4;
for n from 1 to M do:
k_def[n] := k[n] = f(x[i]+c[n]*h,y[i]+h*add(alpha[n,m]*k[m],m=1..M)):
od;
evolve := y[i+1] = y[i] + h*add(b[n]*k[n],n=1..M);A given Runge-Kutta scheme is defined when we assign definate values to the constants LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkobXN1YnN1cEdGJDYoLUkobWZlbmNlZEdGJDYmLUYjNigtSSVtc3ViR0YkNiUtSSNtaUdGJDYlUScmIzk0NTtGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUYjNictRjc2JVEibkYnL0Y7USV0cnVlRicvRj5RJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiLEYnRj0vJSZmZW5jZUdGPC8lKnNlcGFyYXRvckdGRi8lKXN0cmV0Y2h5R0Y8LyUqc3ltbWV0cmljR0Y8LyUobGFyZ2VvcEdGPC8lLm1vdmFibGVsaW1pdHNHRjwvJSdhY2NlbnRHRjwvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GNzYlUSJtRidGRUZHRkVGRy8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRkktRjQ2JS1GNzYlUSJiRidGRUZHLUYjNiVGQkZFRkdGXm9GSS1GNDYlLUY3NiVRImNGJ0ZFRkdGZm9GXm9GPUY9LyUlb3BlbkdRInxmckYnLyUmY2xvc2VHUSJ8aHJGJy1GIzYpRkJGSUZbby1GSjYtUSI9RidGPUZNL0ZQRjxGUUZTRlVGV0ZZL0ZmblEsMC4yNzc3Nzc4ZW1GJy9GaW5GanAtSSNtbkdGJDYkUSIxRidGPUZFRkctRiM2JS1GNzYlUSJNRidGRUZHRkVGRy8lMXN1cGVyc2NyaXB0c2hpZnRHRmBvRl5vL0krbXNlbWFudGljc0dGJFEsW25vbmUsbm9uZV1GJy1GNzYjUSFGJ0Y9. In the literature, these are often presented as a Butcher tableau, which is a matrix arranged in the following way:A := Matrix([seq([seq(alpha[n,m],m=1..M)],n=1..M)]):
B := Vector[row]([seq(b[n],n=1..M)]):
C := Vector[column]([seq(c[n],n=1..M)]):
Butcher := Matrix([[C,A],[``,B]]);You might notice an immediate barrier to implementing this scheme: the equations (1.1) actually involve the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn's on both the lefthand and righthand sides. Hence, it will not be in general possible to analytically solve for the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn's except for very specific forms of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGQUZBRkE=. In this sense, the most general Runge-Kutta methods are implicit schemes: the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn's are defined implicitly. However, we can obtain an explicit scheme by setting a number of the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEnJiM5NDU7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYnLUYvNiVRIm5GJy9GM1EldHJ1ZUYnL0Y2USdpdGFsaWNGJy1JI21vR0YkNi1RIixGJ0Y1LyUmZmVuY2VHRjQvJSpzZXBhcmF0b3JHRj4vJSlzdHJldGNoeUdGNC8lKnN5bW1ldHJpY0dGNC8lKGxhcmdlb3BHRjQvJS5tb3ZhYmxlbGltaXRzR0Y0LyUnYWNjZW50R0Y0LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRi82JVEibUYnRj1GP0Y9Rj8vJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0Y1 coefficients equal to zero:for n from 1 to M do:
for m from n to M do:
alpha[n,m] := 0:
od;
od;
Butcher; With these choices, we see that the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn's are now defined explicitly:for n from 1 to M do;
k_def[n];
od;That is, the first equation allows you to calculate LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0Y+, which can then be used in the second to give you LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMkYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1GLzYjUSFGJ0Y+, etc. This is the general form of an explicit Runge-Kutta LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=-stage stencil. It is traditional (but not necessary) to also set LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiY0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIj1GJ0Y+LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZJLyUpc3RyZXRjaHlHRkkvJSpzeW1tZXRyaWNHRkkvJShsYXJnZW9wR0ZJLyUubW92YWJsZWxpbWl0c0dGSS8lJ2FjY2VudEdGSS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRlgtRjs2JEZCRj5GPg==, resulting in the following Butcher tableau:c[1] := 0;
Butcher, map(x->if (x=0) then `` else x end if,Butcher);Notice that we have written the tableau in two ways, the second just has blanks where all the zeros reside and is how the tableau is usually written in books. Now, how do we choose the remaining coefficients in the Butcher tableau? The idea is to select them in such a way that the local error in the stencil is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSShtZmVuY2VkR0YkNiQtRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJoRicvRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2Jy1GLDYlUSJMRidGQEZCLUkjbW9HRiQ2LVEiK0YnRjIvJSZmZW5jZUdGMS8lKnNlcGFyYXRvckdGMS8lKXN0cmV0Y2h5R0YxLyUqc3ltbWV0cmljR0YxLyUobGFyZ2VvcEdGMS8lLm1vdmFibGVsaW1pdHNHRjEvJSdhY2NlbnRHRjEvJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0Znbi1JI21uR0YkNiRRIjFGJ0YyRkBGQi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGMkYyLUYsNiNRIUYnRjI= where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= is the desired global error in the solution; that isLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkobWZlbmNlZEdGJDYkLUYjNi8tSSNtaUdGJDYlUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYsNiQtRiM2Ji1GMTYlUSJ4RidGNEY3LUkjbW9HRiQ2LVEiK0YnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGSS8lKXN0cmV0Y2h5R0ZJLyUqc3ltbWV0cmljR0ZJLyUobGFyZ2VvcEdGSS8lLm1vdmFibGVsaW1pdHNHRkkvJSdhY2NlbnRHRkkvJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0ZYLUYxNiVRImhGJ0Y0RjdGRUZFLUZCNi1RKCZtaW51cztGJ0ZFRkdGSkZMRk5GUEZSRlRGVkZZRjAtRiw2JC1GIzYkRj5GRUZFLUZCNi1RKiZ1bWludXMwO0YnRkVGR0ZKRkxGTkZQRlJGVEZWRllGZW4tSSttdW5kZXJvdmVyR0YkNictRkI2LVEmJlN1bTtGJ0ZFL0ZIUSZ1bnNldEYnL0ZLRmlvL0ZNRjYvRk9GaW8vRlFGNi9GU0Y2L0ZVRmlvL0ZXUSYwLjBlbUYnL0ZaUSwwLjE2NjY2NjdlbUYnLUYjNiYtRjE2JVEibkYnRjRGNy1GQjYtUSI9RidGRUZHRkpGTEZORlBGUkZUL0ZXUSwwLjI3Nzc3NzhlbUYnL0ZaRl1xLUkjbW5HRiQ2JFEiMUYnRkVGRS1GIzYkLUYxNiVRIk1GJ0Y0RjdGRUZULyUsYWNjZW50dW5kZXJHRkktSSVtc3ViR0YkNiUtRjE2JVEiYkYnRjRGNy1GIzYkRmZwRkUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JJm1mcmFjR0YkNigtRiM2JC1GW3I2JS1GMTYlUSJrRidGNEY3RmByRmJyRkUtSSltcGhhbnRvbUdGJDYkLUYjNiYtRjE2I1EhRictRiM2JkZkcy1GIzYnRmRzRj4tRkI2LUZbcUZFRmhvRmpvL0ZNRmlvRlxwL0ZRRmlvL0ZTRmlvRl9wRlxxRl5xLUYxNiVRImFGJ0Y0RjdGRUZkc0ZFRmRzRkUvJSxjb25zdHJhaW50c0dRLGhlaWdodC1vbmx5RicvJS5saW5ldGhpY2tuZXNzR0Zkci8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0ZqdC8lKWJldmVsbGVkR0ZJLUZCNi1RKSZ2ZXJiYXI7RidGRUZob0Zqb0ZbcEZccEZedEZfdEZfcC9GV1EsMC4xMTExMTExZW1GJy9GWkZjdS1GZnI2KC1GYHM2JC1GIzYmRmRzLUYjNiVGZHMtRjE2JVElZih4KUYnRjRGN0ZFRmRzRkVGY3QtRiM2Ji1GLDYkLUYjNiYtRltyNiVGPi1GIzYkLUYxNiVRImlGJ0Y0RjdGRUZici1GQjYtUSIsRidGRUZHL0ZLRjZGTEZORlBGUkZURmBwL0ZaUSwwLjMzMzMzMzNlbUYnLUZbcjYlRjBGaHZGYnJGRUZFRmlwLUYsNiQtRiM2J0Y+Rl13RjBGW29GRUZFRkVGZnRGaHRGW3VGXXVGRUZFRmlwLUYxNiVRIk9GJy9GNUZJRkUtRiw2JC1GIzYkLUklbXN1cEdGJDYlRmVuLUYjNiYtRjE2JVEiTEYnRjRGN0ZBRl9xRkUvJTFzdXBlcnNjcmlwdHNoaWZ0R0ZkckZFRkUtRkI2LVEiLkYnRkVGR0ZKRkxGTkZQRlJGVEZgcC9GWkZhcEZFMore concretely, to determine the coefficients we need to:Construct an expression for the local error by first subtracting the RHS from the LHS of evolve, subbing in the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn definitions, and then making the substitutions LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkobWZlbmNlZEdGJDYmLUYjNiktSSVtc3ViR0YkNiUtSSNtaUdGJDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiUtRjQ2JVEiaUYnRjdGOkY3RjovJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIixGJy9GO1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjkvJSlzdHJldGNoeUdGTS8lKnN5bW1ldHJpY0dGTS8lKGxhcmdlb3BHRk0vJS5tb3ZhYmxlbGltaXRzR0ZNLyUnYWNjZW50R0ZNLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjE2JS1GNDYlUSJ5RidGN0Y6Rj1GQkZFRkUtRjE2JUZcby1GIzYnRj8tRkY2LVEiK0YnRklGSy9GT0ZNRlBGUkZURlZGWC9GZW5RLDAuMjIyMjIyMmVtRicvRmhuRmhvLUkjbW5HRiQ2JFEiMUYnRklGN0Y6RkJGSUZJLyUlb3BlbkdRInxmckYnLyUmY2xvc2VHUSJ8aHJGJy1GRjYtUSJ+RidGSUZLRmZvRlBGUkZURlZGWEZaL0ZobkZmbi1GRjYtUTAmUmlnaHRUZWVBcnJvdztGJ0ZJRktGZm8vRlFGOUZSRlRGVkZYL0ZlblEsMC4yNzc3Nzc4ZW1GJy9GaG5GXXEtRiw2Ji1GIzYqRjNGRUZcby1GLDYkLUYjNiRGM0ZJRklGRUZcby1GLDYkLUYjNiZGM0Zjby1GNDYlUSJoRidGN0Y6RklGSUZJRklGXnBGYXBGSQ==Expand the result in a Taylor series in LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiaEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= and simplify derivatives using the differential equation LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIidGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4xMTExMTExZW1GJy8lJ3JzcGFjZUdRJjAuMGVtRictSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RidGL0YyRjlGOS1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORmZuLUYsNiVRImZGJ0YvRjItRlE2JC1GIzYnRlUtRjY2LVEiLEYnRjlGOy9GP0YxRkBGQkZERkZGSC9GS0ZPL0ZOUSwwLjMzMzMzMzNlbUYnRitGUEY5RjlGOQ==.Select LUkobWZlbmNlZEc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXJvd0dGJDYoLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEoJmFscGhhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictRiw2Ji1GMjYlUSJuRicvRjZRJXRydWVGJy9GOVEnaXRhbGljRictSSNtb0dGJDYtUSIsRidGOC8lJmZlbmNlR0Y3LyUqc2VwYXJhdG9yR0ZBLyUpc3RyZXRjaHlHRjcvJSpzeW1tZXRyaWNHRjcvJShsYXJnZW9wR0Y3LyUubW92YWJsZWxpbWl0c0dGNy8lJ2FjY2VudEdGNy8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYyNiVRIm1GJ0ZARkJGOC8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRkQtRi82JS1GMjYlUSJiRidGQEZCLUYsNiRGPUY4RmluRkQtRi82JS1GMjYlUSJjRidGQEZCRmFvRmluRjhGOC8lJW9wZW5HUSJ8ZnJGJy8lJmNsb3NlR1EifGhyRic= in such a way that the coefficients of all termsLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbW9HRiQ2LVEvJlByb3BvcnRpb25hbDtGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRjQvJSlzdHJldGNoeUdGNC8lKnN5bW1ldHJpY0dGNC8lKGxhcmdlb3BHRjQvJS5tb3ZhYmxlbGltaXRzR0Y0LyUnYWNjZW50R0Y0LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGQy1JJW1zdXBHRiQ2JS1JI21pR0YkNiVRImhGJy8lJ2l0YWxpY0dRJXRydWVGJy9GMFEnaXRhbGljRictRiM2JS1GSjYlUSJwRidGTUZQRk1GUC8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGLw== are zero (LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEicEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJSZsZTtGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJMRidGL0YyLUY2Ni1RIilGJ0Y5L0Y8RjFGPi9GQUYxRkJGREZGRkgvRktRLDAuMTY2NjY2N2VtRicvRk5GWEY5. Note that our choices have to work for all possible forms of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGQUZBRkE=.Obviously, we want to choose LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= as big as possible in order minimize the one-step error. It turns out for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJSZsZTtGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjRGJ0Y5Rjk= there are enough coefficents in the Butcher tableau to achieve LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJMRidGL0YyRjk=, but for higher LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= the global error cannot be made that small. In the following sections, we will construct explicit Runge-Kutta stencils for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjJGJ0Y5Rjk= and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjRGJ0Y5Rjk= and hence choose LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJMRidGL0YyRjk=.The LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JlEiTUYnLyUlYm9sZEdRJmZhbHNlRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LlEiPUYnRi8vRjZRJ25vcm1hbEYnLyUmZmVuY2VHRjEvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGMS8lKnN5bW1ldHJpY0dGMS8lKGxhcmdlb3BHRjEvJS5tb3ZhYmxlbGltaXRzR0YxLyUnYWNjZW50R0YxLyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTi1JI21uR0YkNiVRIjJGJ0YvRjwvJSdmYW1pbHlHUTBUaW1lc35OZXd+Um9tYW5GJy8lJXNpemVHUSMxMkYnRi9GPA== Huen (modified Euler) methodrestart;In this section, we will construct an example of an LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjJGJ0Y5Rjk= explicit Runge-Kutta stencil. The relevant equations are (see previous section):M := 2;
for n from 1 to M do:
for m from n to M do:
alpha[n,m] := 0:
od;
od;
c[1] := 0;
for n from 1 to M do:
k_def[n] := k[n] = f(x[i]+c[n]*h,y[i]+h*add(alpha[n,m]*k[m],m=1..M)):
od;
evolve := y[i+1] = y[i] + h*add(b[n]*k[n],n=1..M);
A := Matrix([seq([seq(alpha[n,m],m=1..M)],n=1..M)]):
B := Vector[row]([seq(b[n],n=1..M)]):
C := Vector[column]([seq(c[n],n=1..M)]):
Butcher := map(x->if (x=0) then `` else x end if,Matrix([[C,A],[``,B]]));It will be useful to have a set indicating what we want to solve for:unknowns := {seq(seq(alpha[n,m],m=1..n-1),n=2..M),seq(b[n],n=1..M),seq(c[n],n=2..M)};We could have also obtained this by using the indets function to find all the indeterminant quantities (all symbols that have not been assigned a value) in the Butcher tableau:unknowns := indets(Butcher);Actually, this returns the set we want with an extra element `` corresponding to all the blank elements in the Butcher tableau. This can be removed using the minus operator:unknowns := unknowns minus {``};We now need to construct an expression for the local error. First, we form the LHS-RHS of evolve:Error := (lhs-rhs)(evolve);Then, let us substitute in the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn definitions (we convert k_def to a list first so we can put them all in at once:Error := subs(convert(k_def,list),Error);Actually, we see in this expression that there is still a LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0Y+ (due to the fact that the definition of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMkYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0Y+ has LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1GLzYjUSFGJ0Y+ in it), so we need to do this again:Error := subs(convert(k_def,list),Error);This expression has no reference to the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRIm5GJ0YyRjVGMkY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvRjZRJ25vcm1hbEYn's. Note we could have generated the error in a single line by recursively subbing in k_def using the $ operator. It is not hard to convince yourself that the number of times k_def needs to be subsituted in for an LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=-stage method is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=, so we take that into account in our code:Error := subs(convert(k_def,list)$M,(lhs-rhs)(evolve));Now, we need to replace the numeric data with the exact function values they are meant to represent:Subs := {y[i] = y(x),y[i+1] = y(x+h), x[i] = x};
Error := subs(Subs,Error);We now expand the Error in a Taylor series of order LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=:Error := series(Error,h,M+1);We have not yet used the fact that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRInhGJ0YvRjIvRjNRJ25vcm1hbEYnRj1GPQ== is the solution to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIidGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4xMTExMTExZW1GJy8lJ3JzcGFjZUdRJjAuMGVtRictSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RidGL0YyRjlGOS1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORmZuLUYsNiVRImZGJ0YvRjItRlE2JC1GIzYnRlUtRjY2LVEiLEYnRjlGOy9GP0YxRkBGQkZERkZGSC9GS0ZPL0ZOUSwwLjMzMzMzMzNlbUYnRitGUEY5RjlGOQ==. We can use this to remove the derivatives of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRInhGJ0YvRjIvRjNRJ25vcm1hbEYnRj1GPQ== appearing above:derivative[1] := D(y)(x) = f(x,y(x));
for j from 2 to M do:
derivative[j] := subs(derivative[1],convert(diff(derivative[j-1],x),D)):
od;
Error := subs(convert(derivative,list),Error);Now, to determine what conditions the coefficients need to satisfy for the error to be LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSShtZmVuY2VkR0YkNiQtRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJoRicvRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2JS1JI21uR0YkNiRRIjNGJ0YyRkBGQi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGMkYyRjI=, we can rearrange this expression to be a polynomial in LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkobWZlbmNlZEdGJDYmLUYjNjAtSSNtaUdGJDYlUSJoRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiLEYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GMTYlUSJmRidGNEY3RjotSSZtZnJhY0dGJDYoLUY7Ni1RKyZQYXJ0aWFsRDtGJ0Y+L0ZBUSZ1bnNldEYnL0ZERmluL0ZGRmluL0ZIRmluL0ZKRmluL0ZMRmluL0ZORmluRk8vRlNGUS1GIzYoLUYxNiNRIUYnRmVuLUknbXNwYWNlR0YkNiYvJSdoZWlnaHRHUSYwLjBleEYnLyUmd2lkdGhHUSYwLjJlbUYnLyUmZGVwdGhHRltwLyUqbGluZWJyZWFrR1ElYXV0b0YnLUYxNidRInhGJ0Y0LyUrZm9yZWdyb3VuZEdRLFsyMDAsMCwyMDBdRicvJSxwbGFjZWhvbGRlckdGNkY3RmNvRj4vJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmFxLyUpYmV2ZWxsZWRHRkItRjs2LVEifkYnRj5GQC9GREZCRkVGR0ZJRktGTUZPRmBvRlVGOi1GWTYoRmVuLUYjNiYtRjs2LUZnbkY+RkBGaXFGRUZHRklGS0ZNRk9GYG9GZnEtRjE2JVEieUYnRjRGN0Y+RlxxRl9xRmJxRmRxRlVGOi1GMTYlUSJPRicvRjVGQkY+LUYsNiQtRiM2JC1JJW1zdXBHRiQ2JUYwLUYjNiUtSSNtbkdGJDYkUSIzRidGPkY0RjcvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRj5GPkY+Rj4vJSVvcGVuR1EifGZyRicvJSZjbG9zZUdRInxockYnLUY7Ni1RIi5GJ0Y+RkBGaXFGRUZHRklGS0ZNRk9GYG9GPg== This can be accomplishs using the collect command, but we first need to write down a set corresponding to the terms we want to collect together. To do this, first lets get a set with all the indeterminants in Error:vars := indets(Error);Now, we subtract out the unknown coefficients:vars := vars minus unknowns;Finally, there is no need to factor over LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= or LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRInhGJ0YvRjIvRjNRJ25vcm1hbEYnRj0tRiw2I1EhRidGPQ==, so we take these out, and the indets function does not count the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSShtZmVuY2VkR0YkNiQtRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJoRicvRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2JS1JI21uR0YkNiRRIjNGJ0YyRkBGQi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGMkYyLUYsNiNRIUYnRjI= term as an indeterminant, so we add this back in:vars := vars minus {x,y(x)};
vars := vars union {O(h^(M+1))};Now we can rewrite Error as a polynomial in all the unique combinations of the terms in vars using collect. However, collect works best on simple sums (i.e., objects of Maple type `+`) as opposed to series (i.e., objects of Maple type series), so we first convert Error before collecting:whattype(Error);
Error := convert(Error,`+`);
whattype(Error);
Error := collect(Error,vars,'distributed');The 'distributed' option ensures that the coefficients of the resulting polynomial will only involve elements of the unknowns set. We can put each of the coefficients in the polynomial into a set using coeffs:sys := {coeffs(Error,vars)};Actually, this set includes the coefficient of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSShtZmVuY2VkR0YkNiQtRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJoRicvRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2JS1JI21uR0YkNiRRIjNGJ0YyRkBGQi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGMkYyLUYsNiNRIUYnRjI=, which is just 1. We're not really interested in this, so we remove it from the list. After this, we set each coefficient equal to zero to yield a system of equations for LUkobWZlbmNlZEc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXJvd0dGJDYoLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEoJmFscGhhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictRiw2Ji1GMjYlUSJuRicvRjZRJXRydWVGJy9GOVEnaXRhbGljRictSSNtb0dGJDYtUSIsRidGOC8lJmZlbmNlR0Y3LyUqc2VwYXJhdG9yR0ZBLyUpc3RyZXRjaHlHRjcvJSpzeW1tZXRyaWNHRjcvJShsYXJnZW9wR0Y3LyUubW92YWJsZWxpbWl0c0dGNy8lJ2FjY2VudEdGNy8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYyNiVRIm1GJ0ZARkJGOC8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRkQtRi82JS1GMjYlUSJiRidGQEZCLUYsNiRGPUY4RmluRkQtRi82JS1GMjYlUSJjRidGQEZCRmFvRmluRjhGOC8lJW9wZW5HUSJ8ZnJGJy8lJmNsb3NlR1EifGhyRic=.sys := sys minus {1};
sys := map(x->x=0,sys);We can try to attempt to solve this system:sol := solve(sys);We actually fail to get a definate solution, we rather get a one-parameter family of solutions. This is because the system of equations is actually underdetermined; there are many possible solutions. All of these solutions will generate a valid, explicit 2-stage Runge-Kutta stencil, which can be seen if we put sol into Error:subs(sol,Error);Stated another way: there is no unique 2-stage explicit Runge-Kutta method. However, to use the method we do need actual values for LUkobWZlbmNlZEc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXJvd0dGJDYoLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEoJmFscGhhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictRiw2Ji1GMjYlUSJuRicvRjZRJXRydWVGJy9GOVEnaXRhbGljRictSSNtb0dGJDYtUSIsRidGOC8lJmZlbmNlR0Y3LyUqc2VwYXJhdG9yR0ZBLyUpc3RyZXRjaHlHRjcvJSpzeW1tZXRyaWNHRjcvJShsYXJnZW9wR0Y3LyUubW92YWJsZWxpbWl0c0dGNy8lJ2FjY2VudEdGNy8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYyNiVRIm1GJ0ZARkJGOC8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRkQtRi82JS1GMjYlUSJiRidGQEZCLUYsNiRGPUY4RmluRkQtRi82JS1GMjYlUSJjRidGQEZCRmFvRmluRjhGOC8lJW9wZW5HUSJ8ZnJGJy8lJmNsb3NlR1EifGhyRic=, so we need to make a choice. The customary additional assumption is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiYkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMkYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIj1GJ0Y+LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZJLyUpc3RyZXRjaHlHRkkvJSpzeW1tZXRyaWNHRkkvJShsYXJnZW9wR0ZJLyUubW92YWJsZWxpbWl0c0dGSS8lJ2FjY2VudEdGSS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRlgtSSZtZnJhY0dGJDYoLUY7NiRRIjFGJ0Y+RjgvJS5saW5ldGhpY2tuZXNzR0Zqbi8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0Zfby8lKWJldmVsbGVkR0ZJRj4= , which yields a definitive answer for the other constants:assumptions := {b[2]=1/2};
solution := solve(subs(assumptions,sys));
answer := assumptions union solution;This is known as the Heun or modified Euler stencil. Its Butcher tableau issubs(answer,Butcher);Finally, the actual equations one uses in a numerical method are:for n from 1 to M do:
subs(answer,k_def[n]):
od;
subs(answer,evolve);And the local error issubs(answer,Error);Notice that the Huen method gives the same one-step accuracy as the trapezoidal method, but it is fully explicit.The classic fourth order Runge-Kutta schemerestart;
Digits := 14;We now turn our attention to an explicit 4-stage Runge-Kutta method. Note that the code from the previous section works for any value of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=, so we could just re-run it to obtain the equations satisfied by LUkobWZlbmNlZEc2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXJvd0dGJDYoLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEoJmFscGhhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictRiw2Ji1GMjYlUSJuRicvRjZRJXRydWVGJy9GOVEnaXRhbGljRictSSNtb0dGJDYtUSIsRidGOC8lJmZlbmNlR0Y3LyUqc2VwYXJhdG9yR0ZBLyUpc3RyZXRjaHlHRjcvJSpzeW1tZXRyaWNHRjcvJShsYXJnZW9wR0Y3LyUubW92YWJsZWxpbWl0c0dGNy8lJ2FjY2VudEdGNy8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYyNiVRIm1GJ0ZARkJGOC8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnRkQtRi82JS1GMjYlUSJiRidGQEZCLUYsNiRGPUY4RmluRkQtRi82JS1GMjYlUSJjRidGQEZCRmFvRmluRjhGOC8lJW9wZW5HUSJ8ZnJGJy8lJmNsb3NlR1EifGhyRic= when LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjRGJ0Y5Rjk=. But we can take this opportunity to re-write the code a little more compactly:M := 4;
for n from 1 to M do:
for m from n to M do:
alpha[n,m] := 0:
od;
od:
c[1] := 0:
for n from 1 to M do:
k_def[n] := k[n] = f(x[i]+c[n]*h,y[i]+h*add(alpha[n,m]*k[m],m=1..M)):
od:
evolve := y[i+1] = y[i] + h*add(b[n]*k[n],n=1..M):
unknowns := {seq(seq(alpha[n,m],m=1..n-1),n=2..M),seq(b[n],n=1..M),seq(c[n],n=2..M)}:
derivative[1] := D(y)(x) = f(x,y(x)):
for j from 2 to M do:
derivative[j] := subs(derivative[1],convert(diff(derivative[j-1],x),D)):
od:
Subs := {y[i] = y(x),y[i+1] = y(x+h), x[i] = x}:
Error := subs(convert(derivative,list),convert(series(subs(convert(k_def,list)$M,Subs,(lhs-rhs)(evolve)),h,M+1),`+`)):
vars := indets(Error) minus unknowns minus {x,y(x)}:
sys := {coeffs(collect(Error,vars,'distributed'),vars)} minus {1}:
sys := map(x->simplify(x=0),sys);
This is the set of equations that must be satisfied by LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkobWZlbmNlZEdGJDYmLUYjNigtSSVtc3ViR0YkNiUtSSNtaUdGJDYlUSgmYWxwaGE7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYmLUY0NiVRIm5GJy9GOFEldHJ1ZUYnL0Y7USdpdGFsaWNGJy1JI21vR0YkNi1RIixGJ0Y6LyUmZmVuY2VHRjkvJSpzZXBhcmF0b3JHRkMvJSlzdHJldGNoeUdGOS8lKnN5bW1ldHJpY0dGOS8lKGxhcmdlb3BHRjkvJS5tb3ZhYmxlbGltaXRzR0Y5LyUnYWNjZW50R0Y5LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjQ2JVEibUYnRkJGREY6LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGRi1GMTYlLUY0NiVRImJGJ0ZCRkQtRiM2JEY/RjpGW29GRi1GMTYlLUY0NiVRImNGJ0ZCRkRGY29GW29GOkY6LyUlb3BlbkdRInxmckYnLyUmY2xvc2VHUSJ8aHJGJy1GNDYjUSFGJ0Y6 to have a valid 4-stage stencil. We can now ask Maple to solve this system of equations. I have supressed the output as it is very long. (Warning: the following step takes about 30 seconds on my machine. If you execute it on your computer, the number printed out is how much CPU time you used. See ?time() for more details.)st := time():
rk4_solutions := {solve(sys)}:
time()-st;The number of classes of solutions returned by solve can be obtained using the nops command:NN := nops(rk4_solutions);And here is an example of class of solutions:rk4_solutions[1];This is a parameteric solution; i.e., just as in the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjJGJ0Y5Rjk= case, there are an infinite number of solutions for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkobWZlbmNlZEdGJDYmLUYjNigtSSVtc3ViR0YkNiUtSSNtaUdGJDYlUSgmYWxwaGE7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYmLUY0NiVRIm5GJy9GOFEldHJ1ZUYnL0Y7USdpdGFsaWNGJy1JI21vR0YkNi1RIixGJ0Y6LyUmZmVuY2VHRjkvJSpzZXBhcmF0b3JHRkMvJSlzdHJldGNoeUdGOS8lKnN5bW1ldHJpY0dGOS8lKGxhcmdlb3BHRjkvJS5tb3ZhYmxlbGltaXRzR0Y5LyUnYWNjZW50R0Y5LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjQ2JVEibUYnRkJGREY6LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGRi1GMTYlLUY0NiVRImJGJ0ZCRkQtRiM2JEY/RjpGW29GRi1GMTYlLUY0NiVRImNGJ0ZCRkRGY29GW29GOkY6LyUlb3BlbkdRInxmckYnLyUmY2xvc2VHUSJ8aHJGJy1GNDYjUSFGJ0Y6. The "classic fourth order Runge-Kutta method" involves the following assumptions:assumptions := {alpha[4, 1] = 0, alpha[3, 1] = 0, alpha[4, 2] = 0, b[4] = b[1], b[3] = b[2]};This leads to the following solution for the coefficients:sol := solve(subs(assumptions,sys));
answer := sol union subs(sol,assumptions);The Butcher tableau for this stencil isA := Matrix([seq([seq(alpha[n,m],m=1..M)],n=1..M)]):
B := Vector[row]([seq(b[n],n=1..M)]):
C := Vector[column]([seq(c[n],n=1..M)]):
Butcher = subs(answer,map(x->if (x=0) then `` else x end if,Matrix([[C,A],[``,B]])));Here are the basic equationsrk4 := subs(answer,[seq(k_def[n],n=1..M),evolve]);We can confirm that the one-step error is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSShtZmVuY2VkR0YkNiQtRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJoRicvRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2Jy1GLDYlUSJNRidGQEZCLUkjbW9HRiQ2LVEiK0YnRjIvJSZmZW5jZUdGMS8lKnNlcGFyYXRvckdGMS8lKXN0cmV0Y2h5R0YxLyUqc3ltbWV0cmljR0YxLyUobGFyZ2VvcEdGMS8lLm1vdmFibGVsaW1pdHNHRjEvJSdhY2NlbnRHRjEvJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0Znbi1JI21uR0YkNiRRIjFGJ0YyRkBGQi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGMkYyLUYsNiNRIUYnRjI= as expected:simplify(subs(answer,Error));We conclude by writing a procedure that makes use of the classic Runge-Kutta method to solve LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIidGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4xMTExMTExZW1GJy8lJ3JzcGFjZUdRJjAuMGVtRictSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RidGL0YyRjlGOS1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORmZuLUYsNiVRImZGJ0YvRjItRlE2JC1GIzYnRlUtRjY2LVEiLEYnRjlGOy9GP0YxRkBGQkZERkZGSC9GS0ZPL0ZOUSwwLjMzMzMzMzNlbUYnRitGUEY5RjlGOQ==. We first transform the rk4 stencil into a set of mappings:for n from 1 to M+1 do:
K[n] := unapply(rhs(rk4[n]),h,x[i],y[i],seq(k[m],m=1..n-1));
od;Here is a procedure that uses these to calculate a numeric solution in the interval LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJSZpbjtGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JKG1mZW5jZWRHRiQ2Ji1GIzYmLUkjbW5HRiQ2JFEiMEYnRjktRjY2LVEiLEYnRjlGOy9GP0YxRkBGQkZERkZGSC9GS1EmMC4wZW1GJy9GTlEsMC4zMzMzMzMzZW1GJy1GVTYkUSIxRidGOUY5RjkvJSVvcGVuR1EiW0YnLyUmY2xvc2VHUSJdRidGOQ== using LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= steps and with initial data LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUkjbW5HRiQ2JFEiMEYnL0YzUSdub3JtYWxGJ0Y+Rj4tSSNtb0dGJDYtUSI9RidGPi8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGRi8lKXN0cmV0Y2h5R0ZGLyUqc3ltbWV0cmljR0ZGLyUobGFyZ2VvcEdGRi8lLm1vdmFibGVsaW1pdHNHRkYvJSdhY2NlbnRHRkYvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZVLUYsNiVRI3kwRidGL0YyLUZBNi1RIjpGJ0Y+RkRGR0ZJRktGTUZPRlFGU0ZWRj4=RK4 := proc(y0,N)
local h, x, y, i, k:
h := evalf(1/N):
x[0] := 0:
y[0] := y0:
for i from 1 to N do:
x[i] := x[i-1] + h:
k[1] := K[1](h,x[i-1],y[i-1]):
k[2] := K[2](h,x[i-1],y[i-1],k[1]):
k[3] := K[3](h,x[i-1],y[i-1],k[1],k[2]):
k[4] := K[4](h,x[i-1],y[i-1],k[1],k[2],k[3]):
y[i] := K[5](h,x[i-1],y[i-1],k[1],k[2],k[3],k[4]):
od:
[seq([x[i],y[i]],i=0..N)]:
end proc:We'll compare the fourth order Runge-Kutta results with those obtained by the forward Euler stencil as calculated by EulerSol:Euler := proc(y0,N)
local h, x, y, i:
h := evalf(1/N):
x[0] := 0:
y[0] := y0:
for i from 1 to N do:
x[i] := x[i-1] + h:
y[i] := y[i-1] + h*f(x[i-1],y[i-1]):
od:
[seq([x[i],y[i]],i=0..N)]:
end proc:Here is a plot of the output generated by the two stencils compare to the analytic solution of the equation for a particular choice of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGQUZBRkE=.f := (x,y) -> -y^2+2*x;
y0 := 2;
N := 5;
ode := diff(y(x),x)=f(x,y(x));
analytic_sol := rhs(dsolve({ode,y(0)=y0}));
plot([Euler(y0,N),RK4(y0,N),analytic_sol],x=0..1,style=[point$2,line],legend=[`Forward Euler`,`Runge Kutta 4th order`,`analytic solution`],axes=boxed);
We see for even a moderate number of steps, the agreement between the Runge-Kutta method and the analytic solution is remarkable. We can quantify just how much better the Runge-Kutta stencil does by defining a measure of the global error LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEnJiM5NDk7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJ0Yy as the magnitude of the discrepancy between the numerical and actual values of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUkjbW5HRiQ2JFEiMUYnL0YzUSdub3JtYWxGJ0Y+Rj5GPg== for each stencil:epsilon := (y0,N,Method) -> evalf(abs(Method(y0,N)[N+1][2]-eval(analytic_sol,x=1)));Here is a log-log plot of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEnJiM5NDk7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJ0Yy as a function of the number of steps LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=:data[RK4] := [seq([N,epsilon(y0,N,RK4)],N=4..50,2)]:
data[Euler] := [seq([N,epsilon(y0,N,Euler)],N=4..50,2)]:
plots[loglogplot]([data[Euler],data[RK4]],legend=[`Forward Euler`,`Runge Kutta 4th order`],axes=boxed,labels=[`number of steps`,`global error`]);Clearly, the error for the Runge-Kutta method is several orders of magnitude lower. Futhermore, the global error curves both look linear on the log-log plot, which suggests that there is a power law dependence of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEnJiM5NDk7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJ0Yy on LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiTkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEnJiM5NDk7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RKCZUaWxkZTtGJ0YyLyUmZmVuY2VHRjEvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGMS8lKnN5bW1ldHJpY0dGMS8lKGxhcmdlb3BHRjEvJS5tb3ZhYmxlbGltaXRzR0YxLyUnYWNjZW50R0YxLyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGSS1JJW1zdWJHRiQ2JS1GLDYlRi4vRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2JS1JI21uR0YkNiRRIjBGJ0YyRlFGUy8lL3N1YnNjcmlwdHNoaWZ0R0ZaLUklbXN1cEdGJDYlLUYsNiVRIk5GJ0ZRRlMtRiM2Ji1GNjYtUSomdW1pbnVzMDtGJ0YyRjlGO0Y9Rj9GQUZDRkUvRkhRLDAuMjIyMjIyMmVtRicvRktGY28tRiw2JVEiYUYnRlFGU0ZRRlMvJTFzdXBlcnNjcmlwdHNoaWZ0R0ZaLUY2Ni1RLyZQcm9wb3J0aW9uYWw7RidGMkY5RjtGPUY/RkFGQ0ZFRkdGSi1GaG42JS1GLDYlUSJoRidGUUZTLUYjNiVGZW9GUUZTRmhvRjI=. We can determine the power LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= by fitting a power law to the data using Statistics[PowerFit]:`numerically determined global error (Forward Euler) = ` || O(h^(-Statistics[PowerFit](convert(data[Euler],Matrix))[2]));
`numerically determined global error (Runge-Kutta 4th order) = ` || O(h^(-Statistics[PowerFit](convert(data[RK4],Matrix))[2]));This matches our expectation that if the one-step error in a stencil is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSShtZmVuY2VkR0YkNiQtRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJoRicvRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2Jy1GLDYlUSJwRidGQEZCLUkjbW9HRiQ2LVEiK0YnRjIvJSZmZW5jZUdGMS8lKnNlcGFyYXRvckdGMS8lKXN0cmV0Y2h5R0YxLyUqc3ltbWV0cmljR0YxLyUobGFyZ2VvcEdGMS8lLm1vdmFibGVsaW1pdHNHRjEvJSdhY2NlbnRHRjEvJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0Znbi1JI21uR0YkNiRRIjFGJ0YyRkBGQi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGMkYyRjI=, the global error ought to be LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEiT0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSShtZmVuY2VkR0YkNiQtRiM2JC1JJW1zdXBHRiQ2JS1GLDYlUSJoRicvRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2JS1GLDYlUSJwRidGQEZCRkBGQi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGMkYyRjI= for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiTkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJSZHdDtGJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjFGJ0Y5Rjk= [see (3.10) above].NOTE: the above syntax for Statistics[PowerFit] only works in Maple 15. For previous versions, you need to do this:with(LinearAlgebra):
M1,M2 := convert(data[Euler],Matrix),convert(data[RK4],Matrix):
X1,X2 := Column(M1,1),Column(M2,1):
Y1,Y2 := Column(M1,2),Column(M2,2):
`numerically determined global error (Forward Euler) = ` || O(h^(-Statistics[PowerFit](X1,Y1)[2]));
`numerically determined global error (Runge-Kutta 4th order) = ` || O(h^(-Statistics[PowerFit](X2,Y2)[2]));TTdSMApJN1JUQUJMRV9TQVZFLzQzMDcxNDI3NjBYLCUpYW55dGhpbmdHNiI2IltnbCEiJSEhISM6IiYiJiUhRyMiIiIiIiNGKEYpRidGJ0YoIiIhRisjRikiIidGJ0YnRihGKyNGKSIiJEYnRidGJ0YpRi5GJ0YnRidGJ0YsRiU=