# Euler1: one step in Euler's method for the ODE # # dy/dx = f(x,y) # # with initial conditions y(x0) = y0. (See also file "desolvers") # # Example: # # f := (x,y) -> y # Euler1( 0.0, 1.0, 0.1, f); # # Usage: # # Euler( x0, y0, h, f ); # Euler1 := ( x0, y0, h, f) -> y0 + h*f( x0,y0 ): # Euler: repetitively executes Euler1 # # Example: # # f := (x,y) -> y # Euler( 0.0, 1.0, 1.0, 100, f ) # # This computes an approximation to $e$ by solving # dy/dx = y with y(0.0) = 1.0 using Euler's method # with 100 subdivisons of the interval [0,1], returning # y(1.0). Note that to compute $e$ using Maple, do this: # # evalf(E); # # Usage: # # Euler( a, b, Y, n, f) # # The value computed is an approximation to y(b) where # y(x) is a solution to dy/ds = f(x,y) with y(a) = Y. # Euler := proc( x0, x1, y0, n, f) local x, y, h, i: # set up internal variables h := evalf((x1 - x0)/n): # compute step size x := x0; y := y0: # initialize for i from 1 to n do y := Euler1( x, y, h, f ): # compute new y-value x := x + h: # compute new x-value od: RETURN(y); # value of Euler(...) end: Eulerplot := proc( x0, x1, y0, n, f) local i, x, y, h, data: # set up internal variables h := evalf((x1 - x0)/n): # compute step size x := x0; y := y0: # initial point data := array(0..n); # array to hold the data data[0] := [x0,y0]; for i from 1 to n do y := Euler1( x, y, h, f ): # compute new y-value x := x + h: # compute new x-value data[i] := [x,y] od: plot( convert( data, list ) ); # convert array to list & plot end: