#
# First, Euler's method is observed.  Improved Euler's method is then
# observed to be worse for this case.  This is discussed.  Runge-Kutta
# is mentioned but not discussed.
#
# Written by Doug Hensley
#
> 
> 
>  f:=(x,y)->x+sin(y); \
ezeuler:=proc(f,startx,finishx,starty,hop)\
local x,y,h,point;\
global eulerseq;\
y:=starty;h:=hop;point:=[startx,y];eulerseq:=point;\
for x from startx+h by h to finishx do \
     y:=evalf(y+h*f(x-h,y));point:=[x,y];\
     eulerseq:=eulerseq,point od:\
y end;\


                            f := (x,y) -> x + sin(y)

> ezeuler(f,1,10,0,1);
                                   44.34661583
> [eulerseq];

  [[1, 0], [2, 1.], [3, 3.841470985], [4, 6.197346360], [5, 10.11161279],

      [6, 14.47751992], [7, 21.42015681], [8, 28.96062356], [9, 36.32695232],

      [10, 44.34661583]]
> plot(");
> 
> 
# This is all well and good, but it's not all that fast. The defect of the 
# method is that it assumes, # (wrongly, but some compromise is 
# necessary) that the slope of the solution remains fixed as x goes from 
# x to x+h. A better approximation to the average value of the slope 
# over that interval is the slope at x+h/2.
> 
> impvdEuler:=\
proc(f,firstx,lastx,firsty,stepsize)\
local x,y,h,ty,yy,point;\
global impvdEulerlist;\
h:=stepsize;y:=firsty;\
yy:=ezeuler(f,firstx,firstx+stepsize,firsty,stepsize^2);\
point:=[firstx,firsty],[firstx+h,yy];\
impvdEulerlist:=[point];\
for x from firstx+2*h by h to lastx do \
      ty:=yy;\
      yy:=y+2*h*evalf(f(x-h,yy));\
      y:=ty;point:=[x,yy];\
      impvdEulerlist:=[op(impvdEulerlist),point];\
      od;\
yy end;

>  impvdEuler(f,1,10,0,1/10);\
ezeuler(f,1,10,0,1/10);\
plot({[eulerseq],impvdEulerlist});

                                   49.10131736

                                   48.80116568
> 
# Use the larger button in the upper right hand corner of your plot 
# window to blow this up to full-screen. Then select, in the style menu, 
# POINT. This way you can see how the two results are slightly 
# different.
> RK:=proc(f,firstx,lastx,firsty,stepsize)\
local h,x,y,i,k1,k2,k3,k4,point;\
global rklist;\
h:=stepsize;x:=firstx;y:=firsty;point:=[x,y];\
rklist:=[point];\
while x<=lastx-h do \
      k1:=evalf(f(x,y));\
      k2:=evalf(f(x+0.5*h,y+0.5*h*k1));\
      k3:=evalf(f(x+0.5*h,y+0.5*h*k2));\
      k4:=evalf(f(x+h,y+h*k3));\
      y:=y+(h/6)*(k1+2*k2+2*k3+k4);\
      x:=x+h;\
      point:=[x,y];\
      rklist:=[op(rklist),point];\
      od;\
y;end;

>  RK(f,1,10,0,1/10);\


                                   49.14631012
> plot(rklist);
> 
--------------------------------------------------------------------------------
> 
>  ff:=(x,y)->sin(x)+sin(y);

                         ff := (x,y) -> sin(x) + sin(y)
--------------------------------------------------------------------------------
> ezeuler(ff,1,10,0,1/4);

                                   3.380453928
--------------------------------------------------------------------------------
> e1:=[eulerseq]:plot(e1);
--------------------------------------------------------------------------------
> impvdEuler(ff,1,10,0,1/4);

                                   4.782750447
--------------------------------------------------------------------------------
> e2:=impvdEulerlist:plot({e1,e2});
# Wow! Is Improved Euler really an improvement? As you may have 
# heard in calculus, the very simplest method can often be best when 
# you have noisy data. Here, we are using too large a step size. The 
# basic Euler method takes the abuse in stride and gives its usual 
# mediocre approximation. The fine-strung Improved Euler method 
# goes berserk trying to track the DE. The theorem says that f(x,y) 
# must be continuous. At step size 1/4, the continuity of sin(x)+sin(y) 
# isn't much in evidence!
--------------------------------------------------------------------------------
> RK(ff,1,10,0,1/4);

                                   3.311381822
--------------------------------------------------------------------------------
> e3:=rklist:plot({e1,e3});
--------------------------------------------------------------------------------
> impvdEuler(ff,1,10,0,1/50);
--------------------------------------------------------------------------------
> 
# You could use this, or the built-in tool DEplot1, to investigate the 
# solutions to y'=sin(x)+sin(y).