# Problem:  Euler's Method
# 
# Euler's algorithm approximates the solution to the differential equation dy/dx = f(x,y) by 
# moving from point (x_n, y_n) to point (x_n+1, y_n+1) along the line with slope f(x_n, y_n).  
# This implies that
# 
# y_n+1 - y_n 
# -----------------  = f(x_n,y_n)     
#  x_n+1 - x_n                         and    y_n+1 = y+n + h*f(x_n,y_n)                                                   
# where h = x_n+1 - x_n
# 
# 
# h is called the step size.  In the standard Euler's method, step size is kept constant.
# 
# The following examples display the approximations generated by the Euler algorithm for 
# three first order equations.
# 
# Example 1
#                          dy/dx = y + x       y(0) = 0
# 
#                     Because the differential equation is of a special type, first order linear, it can 
# be solved in closed form.  The solution is y=e^x -x-1.  This is the equation that                     
# was used to generate the blue curve shown bleow. 
> with(plots):


> ExactSolution := (x) -> exp(x) - x - 1;


> picture1 := plot(ExactSolution(x), x=0..4, color=COLOR(RGB,0,0,1)):


> h:=0.5;  # the step size
> oldx:=0;  
> oldy:=0;  # the initial condition
> L:=[];
> slope:=(x,y) -> x+ y;
> counter:=0;
> while counter <= 4/h do L:=[op(L),oldx]; L:=[op(L),oldy];oldy:=oldy+h*slope(oldx
> ,oldy);counter:=counter+h;oldx:=counter*h; od;


> picture2:=plot(L,x=0..4,style=point,color=COLOR(RGB,1,0,0)):


> display({picture1,picture2});


# The graph displays the blue solution curve and a red approximation curve (the points are 
# where the aproximations were made).  Try changing the step size, h, to be 0.5, 0.25, and 
# 0.125 (and then reevaluate the region).  How close can you get to the blue curve?
# 
# Example 2
# 
# dy                                                  Although this differential equation is not much more
# ---  = y^2 + x^2       y(0)= 0         complicated than the one in Example 1, it can be shown 
# dx                                                  that it does not have any closed form solutions.
# 
# 
# 
> 


> h:=1;
> oldx:=0;
> oldy:=0; # the initial condition
> L:=[];
> slope:=(x,y) -> y^2 - x^2;
> counter:=0;
> while counter <= 4/h do L:=[op(L),oldx]; L:=[op(L),oldy];oldy:=oldy+h*slope(oldx
> ,oldy);counter:=counter+h;oldx:=counter*h; od;
> plot(L,style=line,color=COLOR(RGB,1,0,0));


# Reduce the step size as above and reevaluate the region.  Notice how drastically the 
# solution changes with a smaller step size.  It's very easy to be misled as to the form of a 
# solution with an inapropriate step size.
# 
# 
# Example 3
# 
# dy
# ---  = sin(y) - x*y         y(0) = 1    This equation does not have any closed form solutions.
# dx
# 
> 


> h:=1;
> oldx:=0;
> oldy:=1; #the initial conditon
> L:=[];
> slope:=(x,y) -> sin(y) - x*y;
> counter:=0;
> while counter <= 4/h do L:=[op(L),oldx]; L:=[op(L),oldy];oldy:=oldy+h*slope(oldx
> ,oldy);counter:=counter+h;oldx:=counter*h; od;
> plot(L, style=line,color=COLOR(RGB,1,0,0));


# Reduce step size as above, and reevaluate region.  No matter what happens to the rest of 
# the solution, the initial condition is always met by this method.
> 


>