# Heun1: one step in Heun'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 # Heun1( 0.0, 1.0, 0.1, f); # # Usage: # # Heun( x0, y0, h, f ); # Heun1:=proc(x0,y0,h,f) local s,s0: s0 := f(x0,y0): s := y0+h*(s0+f((x0+h),(y0+h*s0)))/2: RETURN(s): end: # Heun: repetitively executes Heun1 # # Example: # # f := (x,y) -> y # Heun( 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 Heun's method # with 100 subintervals. Note that to # compute \$e\$ using Maple, do this: # # evalf(E); # # Usage: # # Heun( a, b, Y, h, 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. # Heun := 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 := Heun1( x, y, h, f ): # compute new y-value x := x + h: # compute new x-value od: RETURN(y); # value of Heun(...) end: Heunplot := 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 := Heun1( 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: