Question

Ordinary Differential Equations (a) Write a Python function implementing the 4'th order Runge-Kutta method. (b) Solve the following amusing variation on a pendulum problem using your routine. A pendulum is suspended from a sliding collar as shown in the diagram below. The system is at rest when an oscillating motion y(t) = Y sin (omega t) is imposed on the collar, starting at t = 0. The differential equation that describes the pendulum motion is given by: d^2 theta/dt^2 = -g/L sin theta + omega^2/L Y cos theta sin omega t Plot theta versus t for t = 0 to 40 s and determine the largest value of 0 during this time. Use g = 9.80665 m/s^2, L = 1.0 m, y = 0.25 m, and omega = 2.5 rad/s.

Answer 1

a) answer

def RK4(f):
   return lambda t, y, dt: (
           lambda dy1: (
           lambda dy2: (
           lambda dy3: (
           lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
           )( dt * f( t + dt , y + dy3 ) )
                   )( dt * f( t + dt/2, y + dy2/2 ) )
                   )( dt * f( t + dt/2, y + dy1/2 ) )
                   )( dt * f( t      , y        ) )

def theory(t): return (t**2 + 4)**2 /16

from math import sqrt
dy = RK4(lambda t, y: t*sqrt(y))

t, y, dt = 0., 1., .1
while t <= 10:
   if abs(round(t) - t) < 1e-5:
                print("y(%2.1f)\t= %4.6f \t error: %4.6g" % ( t, y, abs(y - theory(t))))
    t, y = t + dt, y + dy( t, y, dt )

output:

y(0.0)  = 1.000000      error:    0

y(1.0)  = 1.562500      error: 1.45722e-07

y(2.0)  = 3.999999      error: 9.19479e-07

y(3.0)  = 10.562497     error: 2.90956e-06

y(4.0)  = 24.999994     error: 6.23491e-06

y(5.0)  = 52.562489     error: 1.08197e-05

y(6.0)  = 99.999983     error: 1.65946e-05

y(7.0)  = 175.562476    error: 2.35177e-05

y(8.0)  = 288.999968    error: 3.15652e-05

y(9.0)    = 451.562459      error: 4.07232e-05