I am a newcomer to Python, my knowledge of the programming language is still in its infancy, so I copied the Runge-Kutta Python script shown here and modified it for my purposes. Here is my current script:
import numpy as np
from matplotlib import pyplot as plt
a=0
b=np.pi
g=9.8
l=1
N=1000
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 ) )
from math import sqrt
dy = RK4(lambda t, y: -y)
t, y, dt = 0., 1., np.divide((b-a),float(N))
i=0
T=np.zeros((N+1,1))
DY=T
Y=T
while t < (b-dt):
T[i]=t
DY[i]=dy(t,y,dt)
t, y = t + dt, y + dy( t, y, dt )
Y[i]=y
i=i+1
plt.figure(1)
plt.plot(T,Y)
plt.show()
You can ignore the g and l variables in the first few lines, I was going to solve the problem of the simple pendulum, but then I remembered this is a first order ODE solver, so now my ODE is dy/dx=-y. I have been running this in IPython. I was expecting T to be an array, basically the equivalent to linspace(0,pi,N+1) in MATLAB. So T is to be a set of N+1 evenly-spaced values between (and including) 0 and pi, but instead this is a sample of its contents (which came as an output from running T):
In [101]: T
Out[101]:
array([[ 0.99686334],
[ 0.99373651],
[ 0.9906195 ],
...,
[ 0.04334989],
[ 0.04321392],
[ 0. ]])
(including the Input and Output lines to give some context as to what I am referring to here, in case it is unclear). Oh and by-the-way, if you're wondering why I didn't use T=np.linspace(a,b,num=N+1) instead of defining it in this loop well it is because this gives similarly unusual T arrays.
With
DY=T
Y=T
you are only copying the references, they all point to the same array object. Essentially, the contents of all will be Y as it is the last assigned.
Use T.copy() to get a separate array object.