Euler's Method - The Maple Way restart:with(plots):with(DEtools):We begin again by entering the initial value problem.deq:=diff(y(x),x) = 5/sqrt(x+25);
IC:=y(0)=0;We solve the equation, change the right-hand side to a function using unapply, and create the plot structure for our domain. Here we choose not to initialize the array structures.Y:=dsolve({deq,IC},y(x));
Y:=unapply(rhs(Y),x);
p[1]:=plot(Y(x),x=0..700):We use dsolve to implement Euler's method by setting the type to numeric, the method to classical[foreuler] (classical would be sufficient since foreuler, standing for forward Euler method, is the default), and providing a starting point for the independent variable and a stepsize.Y100:=dsolve({deq, IC}, y(x), type=numeric,
method=classical[foreuler], start=0.0,stepsize=100);The output in this case is a procedure from which the data we want can be extracted. Suppose we simply wanted to find the data for when LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkjbW5HRiQ2JFEkMjAwRidGPkY+RitGPg==.Y100(200);We can do the following to create an ordered pair representing a single data point.eval([x,y(x)],Y100(200));We can also simply extract the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= value at a data point.eval(y(x),Y100(200));We can do the following to see the whole set of data.for i from 0 to 7 do
Y100(i*100);
od;We can also choose to see this as a listing of data points.for i from 0 to 7 do
eval([x,y(x)],Y100(100*i));
od;Suppose we ask for a value, such as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkjbW5HRiQ2JFEkMjUwRidGPkY+RitGPg==, which falls between the computed points. Maple uses linear interpolation, that is, it follows the line segment joining the nearest two computed points.Y100(250);We can again choose to see this simply as a data point.eval([x,y(x)],Y100(250));We will use odeplot from the plots package to create first the point plot and then the line plot to connect the points. p[2]:=odeplot(Y100,[x,y(x)],0..700, numpoints=7, style=point,color=blue):
p[3]:=odeplot(Y100,[x,y(x)],0..700, numpoints=7, style=line,color=blue):
display(seq(p[n],n=1..3));Let's repeat the whole process for a stepsize of 50, giving 14 steps.Y50:=dsolve({deq, IC}, y(x), type=numeric,
method=classical[foreuler], start=0.0,stepsize=50):
p[4]:=odeplot(Y50,[x,y(x)],0..700, numpoints=14, style=point,color=magenta):
p[5]:=odeplot(Y50,[x,y(x)],0..700, numpoints=14, style=line,color=magenta):
display(seq(p[n],n=1..5));Finally, let's repeat it one more time with a stepsize of 10 for 70 steps.Y10:=dsolve({deq, IC}, y(x), type=numeric,
method=classical[foreuler], start=0.0,stepsize=10):
p[6]:=odeplot(Y10,[x,y(x)],0..700, numpoints=70, style=point,color=green):
p[7]:=odeplot(Y10,[x,y(x)],0..700, numpoints=70, style=line,color=green):
display(seq(p[n],n=1..7));It seems clear that the smaller the stepsize, the more accuracy one gets. However, even our smallest stepsize here is still quite large. If one does not include the stepsize command, Maple uses the minimum of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJoRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkjbW5HRiQ2JFEmMC4wMDVGJ0Y+Rj5GK0Y+ and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJoRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbWZyYWNHRiQ2KC1GIzYmRistRiM2Ji1GLDYlUSJiRidGNEY3LUY7Ni1RKCZtaW51cztGJ0Y+RkBGQ0ZFRkdGSUZLRk0vRlBRLDAuMjIyMjIyMmVtRicvRlNGXG8tSSVtc3ViR0YkNiUtRiw2JVEieEYnRjRGNy1GIzYkLUkjbW5HRiQ2JFEiMEYnRj5GPi8lL3N1YnNjcmlwdHNoaWZ0R0Zpb0Y+RitGPi1GIzYkLUZnbzYkUSIzRidGPkY+LyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0ZmcC8lKWJldmVsbGVkR0ZCRj5GK0Y+ where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= = LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMEYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHRj1GPg== .. LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiYkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=. Let's look at the line display without the points.display(p[1],p[3],p[5],p[7]);Finally, let's print a table of our values. In the table, LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEleTEwMEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=, for instance, refers to the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=-values from a step size of 100, while LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= refers to the exact function values.printf(" n x[n] y100(x[n]) y50(x[n]) y10(x[n]) y(x[n])");
printf("\134n\134n");
for n from 0 to 7 do
printf("%2d %4.0f %12.4f %12.4f %12.4f %12.4f \134n",n,eval(x,Y100(100*n)),eval(y(x),Y100(100*n)),eval(y(x),Y50(100*n)),eval(y(x),Y10(100*n)),evalf(Y(100*n)));
od;In Chapter 3, we used DEplot to draw direction fields along with solutions. Those solutions were not the exact solutions through a point, but used a numeric method called the fourth order classical Runge-Kutta method as a default. We can choose other methods such as the Euler method. It is also good to know that the default stepsize is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJoRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkmbWZyYWNHRiQ2KC1GIzYmRistRiM2Ji1GLDYlUSJiRidGNEY3LUY7Ni1RKCZtaW51cztGJ0Y+RkBGQ0ZFRkdGSUZLRk0vRlBRLDAuMjIyMjIyMmVtRicvRlNGXG8tRiw2JVEiYUYnRjRGN0Y+RitGPi1GIzYkLUkjbW5HRiQ2JFEjMjBGJ0Y+Rj4vJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRlxwLyUpYmV2ZWxsZWRHRkJGPkYrRj4= where x = a .. b. Let's compare the actual solution (red) with both Euler's method (green) and the Runga-Kutta method (blue).YY10:=dsolve({deq, IC}, y(x), type=numeric,
method=classical[rk4], start=0.0,stepsize=10):
p[8]:=odeplot(YY10,[x,y(x)],0..700, numpoints=70, style=line,color=blue):
display(p[1],p[7],p[8]);