Solving for motion in one dimension using the second order the Runge-Kutta method.
D. E. Soper
13 October 2003

Off[General::"spell"];Off[General::"spell1"]

2nd order Runge-Kutta operator

The function rk2 has inputs
   t = total length of  time for which the equation m*d^2x/dt^2 = f[y] is to be solved.
   n = number of intervals into which this total proper time is to be divided. Thus the
         time step is Δt = t/n.
   x0 = starting value of x
   v0 = starting value of v.
   
The output is a list of x values: xlist = {x[[1]], x[[2]], ...,x[[n]]}, where each x[[j]] is x  at time step j.

Required functions:
   force[x]
   friction[v]

rk2[t_,n_,x0_,v0_]:= Module[
(* The internal variables *)
{Δt,v,x,vmid,xmid,txlist},  
Δt = t/n; (* the time step *)
v = v0;   (* initial value of v *)
x = x0;   (* initial value of x *)
txlist = {{0,x}};  (* starting output list *)
Do[  (* repleat steps below for j from 1 to n *)
(* Approximate x and v at midpoint of step, *)
(* using the differential equation *)
xmid = x + v * Δt/2 //N; 
vmid = v + friction[v]/m * Δt/2
         + force[x]/m * Δt/2      //N;
(* Use these x and v values to go from begining *)
(* of the step to the end according to the      *)
(* differential equation.                       *)
x = x + vmid * Δt //N;
v = v + friction[vmid]/m * Δt
      + force[xmid]/m * Δt       //N;
(* Append the t and x values to result list *)
txlist = Append[txlist,{j*Δt,x}],
{j,1,n}];
txlist  (* report the result list *)
]

Define mass and force

First, we define a global variable for the mass but simply set it to 1. This setting
could be changed latter if desired.

m = 1;

Here is an example of a force

force[x_] := - Sin[x];

Here is an example of a frictional force

friction[v_] := 0;

Try it

Plot[{force[x],-x},{x,-1,1}]

txlist = rk2[10,30,0,1];

plot1 = ListPlot[txlist,
PlotStyle -> {Hue[0.6],PointSize[0.02]}]

xlist = rk2[10,200,0,1];

plot2 = ListPlot[xlist, PlotStyle -> Hue[0.9]]

Show[plot1,plot2]

Try it with linear force

force[x_] := - x;

xlist = rk2[10,200,0,1];

plot3 = ListPlot[xlist, PlotStyle -> Hue[0.5]]

Show[plot2,plot3]

See if the initial conditions matter

force[x_] := - x;

xlist = rk2[10,200,0,0.1];

plot4 = ListPlot[xlist, PlotStyle -> Hue[0.5]]

force[x_] := - Sin[x];

xlist = rk2[10,200,0,0.1];

plot5 = ListPlot[xlist, PlotStyle -> Hue[0.9]]

Show[plot4,plot5]

Linear force with friction

force[x_] := - x;

friction[v_] := - 1.0 v;

xlist = rk2[10,200,0,1];

plot3 = ListPlot[xlist, 
PlotRange -> All,
PlotStyle -> Hue[0.8]]