# rk1: one step in the standard: Runge-Kutta method for the ODE # # dy/dx = f(x,y) # # with initial conditions y(x0) = y0. (see also file "desolvers") # # Example: # # f := (x,y) -> y # rk1( 0.0, 1.0, 0.1, f); # # Usage: # # rk( x0, y0, n, f ); # rk1:=proc(x0,y0,h,f) local k1,k2,k3,k4,s: k1:=h*f(x0,y0): k2:=h*f((x0+h/2),(y0+k1/2)): k3:=h*f((x0+h/2),(y0+k2/2)): k4:=h*f((x0+h),(y0+k3)): s:=y0+(k1+2*k2+2*k3+k4)/6: RETURN(s): end: # rk: repetitively executes rk1 # # Example: # # f := (x,y) -> y # rk( 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 the Runge-Kutta # method with 100 subdvisions. Note that to # compute \$e\$ using Maple, do this: # # evalf(E); # # Usage: # # rk( 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. # rk := 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 := rk1( x, y, h, f ): # compute new y-value x := x + h: # compute new x-value od: RETURN(y); # value of rk(...) end: rkplot := 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 := rk1( 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: