Numeric Integration using Python and Pylab

Pylab is Python module that contain NumPy and MatPlotlib.

Let f(x) = x^2 and we want to integrate it from 0 to 1.

Numerically we could use square method or trapezoidal (which is almost better).

Delphi.

Still using Delphi?

Yup, it teach me about programming discipline, :)

Flappy Bird Like in Delphi

Remember the infamous Flappy Bird? Yup, I will create the program based on that algorithm.

Newton Polynomial in Python

I wrote code for this in Delphi. This time I want to rewrite it in Python based on this wiki.

I use this set of data point
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)

and I use xc=3 for the test data.

It's obvious that these sets of data points have quadratic form and f(xc) must have value of 9.

The heart of code lay on this

def n(j,xc,x):
    n = 1
    for i in arange(j):
        n *= (xc-x[i])

    return n
def a(j,l,x,y):
    if j==0:
        return y[0]
    elif j-l==1 :
        return (y[j]-y[l])/(x[j]-x[l])
    else:
        return (a(j,l+1,x,y)-a(j-1,l,x,y))/(x[j]-x[l])
        

def N(xc,x,y):
    N = 0
    for j in range(len(x)):
        N += a(j,0,x,y)*n(j,xc,x)
    
    return N

Look at the function a(j,l,x,y), that's recursive function to obtain Newton divided difference value.

the whole code is this

from pylab import *

def n(j,xc,x):
    n = 1
    for i in arange(j):
        n *= (xc-x[i])

    return n
def a(j,l,x,y):
    if j==0:
        return y[0]
    elif j-l==1 :
        return (y[j]-y[l])/(x[j]-x[l])
    else:
        return (a(j,l+1,x,y)-a(j-1,l,x,y))/(x[j]-x[l])
        

def N(xc,x,y):
    N = 0
    for j in range(len(x)):
        N += a(j,0,x,y)*n(j,xc,x)
    
    return N

x = []
y = []
#initial value
x.append(0)
x.append(1)
x.append(2)
x.append(4)
x.append(5)

y.append(0)
y.append(1)
y.append(4)
y.append(16)
y.append(25)

#for testing
xc = 3

yc = N(xc,x,y)

print ''
print xc, yc
#plot
t = linspace(-7,7,100)
u = N(t,x,y)
plot(t,u)
grid(True)
show()

.


here's the graphics

with another sets of data points, I have another result
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)

Newton Polinomial.

Here's code for Newton's divided differences interpolation polynomial (quite mouthful huh, :) ).

The purpose of this method is to create a function (polynomial) that passes through given set of data points.




I read data point from several edit box.

procedure TForm1.Button3Click(Sender: TObject);
var i:integer;
begin
  for i:=0 to n do begin
    x[i]:=strToFloat(kx[i].Text);
    y[i]:=strToFloat(ky[i].Text);
  end;
  xc:=strToFloat(kxc.Text);
  yc:=fn(xc);
  kyc.Text:=floatToStr(yc);
  gambarNewton;
end;

.
kx and ky is tEdit created when button1 is clicked

procedure TForm1.Button1Click(Sender: TObject);
var i:integer;
begin
  button2.Enabled:=true;
  button3.Enabled:=true;
  button4.Enabled:=true;
  button5.Enabled:=true;

  n:=strToInt(edit1.Text);
  kxc:=tEdit.Create(form1);       kyc:=tEdit.Create(form1);
  kxc.Parent:=form1;              kyc.Parent:=form1;
  kxc.Left:=36;                   kyc.Left:=72;
  kxc.Width:=36;                  kyc.Width:=36;
  kxc.Text:='0,5';

  for i:=0 to n do begin
    kx[i]:=tEdit.Create(form1);     ky[i]:=tEdit.Create(form1);
    kx[i].Parent:=form1;            ky[i].Parent:=form1;
    kx[i].Top:=36+36*i;             ky[i].Top:=36+36*i;
    kx[i].Left:=36;                 ky[i].Left:=72;
    kx[i].Width:=36;                ky[i].Width:=36;
    kx[i].Text:=intToStr(i);        ky[i].Text:=intToStr(i);
  end;
end;

xc is x coordinate where the corresponding y (yc) is obtained using Newton method by calling it

yc=fn(xc)

function tform1.fn(xs:real):real;
var i:integer;fs:real;
begin
  fs:=0;
  for i:=0 to n do begin
    fs:=fs+b(i,0)*c(xs,i);
  end;
  fn:=fs;
end;

the fn function call the two other function. The b function, a recursive contain divided difference like this

function tform1.b(i,j:integer):real;
begin
  if i=0 then b:=y[0]
    else if (i-j)=1 then
      b:=(y[i]-y[j])/(x[i]-x[j])
        else
          b:=(b(i,j+1)-b(i-1,j))/(x[i]-x[j]);
end;

and c function, a recursive function (or you could rewrite it using simple for command)

function tform1.c(xs:real;i:integer):real;
begin
  if i=0 then c:=1
    else c:=(xs-x[i-1])*c(xs,i-1);
end;

and finally, draw the data and the function on image1

fprocedure tform1.gambarNewton;
var i,x0,y0:integer;px,py:real;
begin
  x0:=image1.Width div 2;               y0:=image1.Height div 2;
  image1.Canvas.Brush.Color:=clLime;
  image1.Canvas.Rectangle(0,0,image1.Width,image1.Height);
  image1.Canvas.Brush.Color:=clWhite;
  image1.Canvas.Pen.Color:=clBlack;
  image1.Canvas.MoveTo(0,y0);           image1.Canvas.LineTo(image1.Width,y0);
  image1.Canvas.MoveTo(x0,0);           image1.Canvas.LineTo(x0,image1.Height);
  for i:=-300 to 300 do begin
    px:=i/skala;                          py:=skala*fn(px);
    image1.Canvas.Pixels[x0+i,y0-round(py)]:=clGreen;
  end;
  for i:=0 to n do begin
    px:=x0+skala*x[i];                        py:=y0-skala*y[i];
    image1.Canvas.Ellipse(round(px)-7,round(py)-7,round(px)+7,round(py)+7);
  end;
  px:=x0+skala*xc;                        py:=y0-skala*yc;
  image1.Canvas.Brush.Color:=clred;
  image1.Canvas.Ellipse(round(px)-7,round(py)-7,round(px)+7,round(py)+7);
  image1.Canvas.Brush.Color:=clwhite;
end;

Short Function.

Here's my implementation of function according to The Power of 10;

"Restrict functions to a single printed page."

As bonus, I didn't use global variable if possible. So if a function or procedure need a variable from others, it have to be passed using parameter on that function.

If we look at the code below, we know that it can be rewritten using a long single procedure or function. But according The Power of Ten, a function should be as short as possible so it could be printed in a single page.

So, instead one long multiple page function, I write/break it as several short-single-printed-page functions. :)

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;

type
  TForm1 = class(TForm)
    Edit1: TEdit;
    Edit2: TEdit;
    Button1: TButton;
    Edit3: TEdit;
    Edit4: TEdit;
    procedure proses;
    function konversi(a:real;c,d:char):string;
    function konversiC(a:real;d:char):string;
    function konversiF(a:real;d:char):string;
    function konversiR(a:real;d:char):string;
    function konversiK(a:real;d:char):string;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.Button1Click(Sender: TObject);
begin
  proses;
end;

procedure tform1.proses;
var
  a:real;
  s,b:string;
  c,d:char;
begin
  a:=strToFloat(edit1.Text);
  s:=edit2.Text;
  c:=s[1];
  s:=edit4.Text;
  d:=s[1];
  b:=konversi(a,c,d);
  edit3.Text:=b;
end;

function tform1.konversi(a:real;c,d:char):string;
begin
  case c of
    'C':konversi:=konversiC(a,d);
    'F':konversi:=konversiF(a,d)
    'R':konversi:=konversiR(a,d)
    'K':konversi:=konversiK(a,d)
    else konversi:='error';
  end;
end;

function tform1.konversiC(a:real;d:char):string;
begin
  case d of
    'C':konversiC:=floatToStr(a);
    'F':konversiC:=floatToStr(a*9/5+32);
    'R':konversiC:=floatToStr(a*4/5);
    'K':konversiC:=floatToStr(a+273);
    else konversiC:='Error';
  end;
end;

function tform1.konversiF(a:real;d:char):string;
begin
  case d of
    'C':konversiF:=floatToStr((a-32)*5/9);
    'F':konversiF:=floatToStr(a);
    'R':konversiF:=floatToStr((a-32)*4/9);
    'K':konversiF:=floatToStr((a-32)*5/9+273);
    else konversiF:='Error';
  end;
end;

function tform1.konversiR(a:real;d:char):string;
begin
  case d of
    'C':konversiR:=floatToStr(a*5/4);
    'F':konversiR:=floatToStr(a*9/4+32);
    'R':konversiR:=floatToStr(a);
    'K':konversiR:=floatToStr(a*5/4+273);
    else konversiR:='Error';
  end;
end;

function tform1.konversiK(a:real;d:char):string;
begin
  case d of
    'C':konversiK:=floatToStr(a-273);
    'F':konversiK:=floatToStr((a-273)*9/5+32);
    'R':konversiK:=floatToStr((a-273)*4/5);
    'K':konversiK:=floatToStr(a);
    else konversiK:='Error';
  end;
end;

end.

.

Another Turtle in Circle

There's always another way to solve something.

So, I have another code for "Turtle in Circle" code, :)

In the script below, I use turtle position to determine if it's still inside circle or not. If it outside circle, instead of send it to zero position, I send it to random position inside circle.

import turtle
from random import uniform
import numpy as np

turtle.shape("turtle")
#turtle.speed(1)
x = 0
y = 0
rmax=40

for i in range (1,1000):
    a = uniform (-90,90) #angle
    turtle.left(a)
    d = uniform (-75,75) #distance
    x = turtle.xcor()+d*np.cos(a*np.pi/180)
    y = turtle.ycor()+d*np.sin(a*np.pi/180)
    r = np.sqrt(x*x+y*y)
    if r>rmax:
        turtle.setx(uniform(-rmax,rmax))
        turtle.sety(uniform(-rmax,rmax))
        x = 0
        y = 0
    else:
        turtle.forward(d)

turtle.exitonclick()



 

.

Turtle in Circle

I use previous code and improve it so the turtle could only move at certain circle area.

import turtle
from random import uniform
import numpy as np

turtle.shape("turtle")
#turtle.speed(1)
x = 0
y = 0

for i in range (1,1000):
    a = uniform (-90,90) #angle
    turtle.left(a)
    d = uniform (-75,75) #distance
    x += d*np.cos(np.pi*a/180)
    y += d*np.sin(np.pi*a/180)
    r = np.sqrt(x*x+y*y)
    if r>40:
        turtle.setx(0)
        turtle.sety(0)
        x = 0
        y = 0
    turtle.forward(d)

turtle.exitonclick()



 

.

Random Turtle Movement.

I use turtle module, the standard module, in Python.

The turtle movement has random direction (angle), and random distance (forward/backward).

import turtle
from random import uniform

turtle.shape("turtle")
turtle.speed(1)
for i in range (1,100):
    #random angle
    a = uniform (-90,90)
    turtle.left(a)
    #random move
    d = uniform (-100,100)
    turtle.forward(d)

turtle.exitonclick()



 

.

Prime Number on Python

We've done it using Delphi, how about Python? Easy, :)

for i in range (2,200):
    prime = True
    for j in range (2,i):
        if i%j==0: 
            prime=False
    if prime==True:
        print i

.

Prime Number on Delphi

Using Delphi to generate prime number? Okay, there's some code out there with

if (i mod 2 <>0) and (i mod 3 <>0) and (i mod 5 <>0) and (i mod 7 <>0) then i is prime.

Don't use that.

That's just for prime number below 100.

Use this instead.

var i,j:integer;
   prime:boolean;
begin
  for i:=2 to 200 do begin
    prime:=true;
    for j:=2 to i-1 do begin
      if i mod j=0 then prime:=false;
    end;
    if prime=true then memo1.Lines.Append(intToStr(i));
  end;
end;

.

Collision.

 Here's the Code

#code
from visual import *
from random import uniform

display(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

Ball    = sphere(radius=2, length=4, opacity=.3)

Bola    = []
n       = 5
for i in arange (n):
    bola        = sphere(color=color.green,radius=uniform(.2,.73))
    bola.pos    = vector(uniform(-1.5,1.5),uniform(-1.5,1.5),uniform(-1.5,1.5))
    bola.v      = vector(uniform(-1,1),uniform(-1,1),uniform(-1,1))
    Bola.append(bola)
    
dt          = 1./16

def pantul():
    global Bola
    for bola in Bola:
        r   = bola.pos
        v   = bola.v
        if mag(r)>=Ball.radius:
            r   = 1.9*norm(r)
            vp  = (dot(v,norm(r)))*norm(r)
            vr  = v-vp
            v   = vr - vp
            bola.r  = r
            bola.v  = v
    for i in arange (n-1):
        for j in range(i+1, n):
            ri  = Bola[i].pos
            rj  = Bola[j].pos
            vi  = Bola[i].v
            vj  = Bola[j].v
            rc  = rj-ri
            if Bola[i].radius+Bola[j].radius>mag(rc):
                vpi = dot(vi,norm(rc))*norm(rc)
                vri = vi-vpi
                vpj = -dot(vj,norm(rc))*norm(rc)
                vrj = vj-vpj

                vi  = vpj+vri
                vj  = vpi+vrj

                Bola[i].v   = vi
                Bola[i].v   = vj

def proses():
    for bola in Bola:
        r   = bola.pos
        v   = bola.v
        a   = vector(0,0,0)
        v   += a*dt
        r   += v*dt

        bola.pos  = r

    pantul()

while 1:
    rate(37)
    proses()

.

Bouncing Ball inside a Cone

  I use vector projection and rejection to calculate velocity after bouncing the side of cone, :)

#code
from visual import *
from random import uniform

display(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

Cone        = cone(pos = (0,0,0), axis=(0,5,0), radius = 3, opacity = .2) 


bola        = sphere(color=color.green,radius=.2)
bola.y      = 1
bola.x      = -1
bola.z      = 1

v           = vector(1,-1,0)
dt          = 1./16
r           = bola.pos
rc          = Cone.radius
h           = vector(Cone.axis)

def pantul():
    global r,v
    #tumbukan dengan lantai
    if r.y<0:
        r.y = 0
        v.y *= -1

    rp      = vector(r.x,0,r.z)
    hb      = h.y - r.y
    rmaks   = hb/h.y*rc
    c       = h-rmaks*norm(rp)    #vektor garis singgung 
                                  #selimut kerucut dengan bidang singgung
    #tumbukan dengan selimut kerucut
    if mag(rp)>rmaks:
        rp  = norm(rp)*rmaks
        r   = vector(rp.x,r.y,rp.z)
        vp  = dot(v,norm(c))*norm(c)
        v   = 2*vp-v
    print v



def proses():
    global r,v
    a   = vector(0,0,0)
    v   += a*dt
    r   += v*dt

    bola.pos  = r

    pantul()

while 1:
    rate(37)
    proses()

.

Bouncing Ball inside Sphere

from visual import *
from random import uniform

display(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

Ball    = sphere(radius=2, length=4, opacity=.3)


bola        = sphere(color=color.green,radius=.2)
bola.y      = 1
bola.x      = -1
bola.z      = 1

v           = vector(2,1,0)
dt          = 1./16
r           = bola.pos
rc          = Ball.radius

def pantul():
    global r,v
    if mag(r)>=rc:
        r   = 1.99*norm(r)
        vp  = (dot(v,norm(r)))*norm(r)
        vr  = v-vp
        v   = vr - vp

def proses():
    global r,v
    a   = vector(0,0,0)
    v   += a*dt
    r   += v*dt

    bola.pos  = r

    pantul()

while 1:
    rate(37)
    proses()

.

Bouncing inside Cylinder

from visual import *
from random import uniform

display(center=(0,2,0),background=(1,1,1), autoscale=False, range=4.5,        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

silinder    = cylinder(radius=2, length=4, opacity=.3)
silinder.rotate(angle=pi/2, axis=(0,0,1),origin=(0,0,0))


bola        = sphere(color=color.green,radius=.2)
bola.y      = 0
bola.x      = 0
bola.z      = 1

v           = vector(2,0,0)
dt          = 1./16
r           = bola.pos
rc          = silinder.radius

def pantul():
    global r,v
    if mag(r)>=rc:
        r   = 1.9*norm(r)
        vp  = (dot(v,norm(r)))*norm(r)
        vr  = v-vp
        v   = vr - vp

def proses():
    global r,v
    a   = vector(0,0,0)
    v   += a*dt
    r   += v*dt

    bola.pos  = r

    pantul()

while 1:
    rate(7)
    proses()

.

Bounce Over Spherical Surface

#code

from visual import *
from random import uniform


display(center=(0,2,0),background=(1,1,1), autoscale=False, range=7.5,        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)

ball        = sphere(radius=2, color=color.red, opacity = .5)
r2          = ball.radius
silinder    = cylinder(radius=2, length=4, opacity=.3)
silinder.rotate(angle=pi/2, axis=(0,0,1),origin=(0,0,0))


bola        = sphere(color=color.green,radius=.2)
bola.y      = 3
bola.x      = uniform(-1,1)
bola.z      = uniform(-1,1)


v           = vector(0,2,0)
dt          = 1./8.
r           = bola.pos

def pantul():
    global r,v
    print v
    if mag(r)<r2:
        print mag(r)
        arah    = norm(r)
        dv      = dot(v,arah)
        v       -= dv*arah
        r       = (r2+.2)*arah
    
    

def proses():
    global r,v
    a   = vector(0,-1,0)
    v   += a*dt
    r   += v*dt

    bola.pos  = r

    pantul()

    

while 1:
    rate(11)
    proses()

.

N-Spring System

 Using Visual Python

 I like the result, :)

#code
from visual import *

n       = 13
display(center=(n/2,0,0),background=(1,1,1), autoscale=False, range=(7),        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.orange)
dt      = 1./8.
dx      = []
k       = []
m       = []
l0      = []
l       = []
v       = []
x       = []
for i in arange(n):
    dx.append(.1)
    k.append(1.)
    m.append(1.)
    l0.append(1.)
    l.append(1.)
    v.append(0.)
    x.append(0.)
pegas   = []
kotak   = []

for i in arange(n):
    spring  = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
    pegas.append(spring)
    ko      =  box(pos=(0,0,0), width=.5, height=.5, length= .5, color= color.green)
    kotak.append(ko)

box(pos=(-1,.64,0), width=n, height=2, length= 2, color=color.black)
box(pos=(n/2.,-.36,0), width=n, height=.2, length= n, color=color.white,opacity=.9)

#usikan
l[0]        = 1.2 
l[n-1]      = .9
#posisi x
position    = 0
for i in arange(n):
    position    +=l[i]
    x[i]        = position
    kotak[i].x  = x[i]


def updatePegas():
    global l
    for i in arange(n):
        if i!=0:
            pegas[i].x  = x[i-1]
            l[i]        = x[i]-x[i-1]
        else:
            l[i]        = x[i]

        kotak[i].x      = x[i]
        pegas[i].length = l[i]
        
def proses():
    for i in arange(n):
        dx[i]   = l[i]-l0[i]
        f0      = -k[i]*dx[i]
        if i<n-1:
            dx[i+1] = l[i+1]-l0[i+1]
            f1      = -k[i+1]*dx[i+1]
            a       = (f0-f1)/m[i]
        else:
            a       = f0/m[i]
        v[i]    += a*dt
        x[i]    += v[i]*dt
    updatePegas()
while 1:
    rate (19)
    proses()

.

Here's the Culprit

 In Visual Python, helix object will generate error if helix.length = some array like the code below. I use dl, an array, for the length value

from visual import *

display(center=(1,0,0),background=(1,1,1), autoscale=False, range=(2,2,2),        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)
n       = 1
dt      = 1./8.
dl      = ones(1)

pegas   = []

print dl


for i in arange(n):
    spring  = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
    pegas.append(spring)

def updatePegas(l):
    pegas[0].length = l
    
    

def proses():
    global dl
    l   = pegas[0].length
    if l>2:
        l   = 2
        dl[0]  *= -1
    elif l<.5:
        l   = .5
        dl[0]  *= -1
    l   += dl[0]
    updatePegas(l)
    
while 1:
    rate (19)
    proses()

    
    
    
    

    

.



And the result is

Traceback (most recent call last):
  File "springList.py", line 40, in 
    proses()
  File "springList.py", line 36, in proses
    updatePegas(l)
  File "springList.py", line 22, in updatePegas
    pegas[0].length = l
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/VPython-6.11-py2.7-macosx-10.6-intel.egg/visual_common/primitives.py", line 850, in set_length
    self.__frame.axis = self.__axis.norm()
AttributeError: 'numpy.ndarray' object has no attribute 'norm'

If we change dl from array to list, like the code below, everything is suddenly OK, :)

from visual import *

display(center=(1,0,0),background=(1,1,1), autoscale=False, range=(2,2,2),        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)
n       = 1
dt      = 1./8.
#dl      = ones(1)
dl      = []
dl.append(.1)
pegas   = []

print dl


for i in arange(n):
    spring  = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
    pegas.append(spring)

def updatePegas(l):
    pegas[0].length = l
    
    

def proses():
    global dl
    l   = pegas[0].length
    if l>2:
        l   = 2
        dl[0]  *= -1
    elif l<.5:
        l   = .5
        dl[0]  *= -1
    l   += dl[0]
    updatePegas(l)
    
while 1:
    rate (19)
    proses()

    
    
    
    

    

.

Gonna rewrite the code.

Double Spring System

Akhirnya.

With Visual Python module

I couldn't use list for spring length since it'll trigger some error for helix object. It's very unfortunate because it would come handy as we expand the number of spring and mass.

Anyway, here's the code

#code
from visual import *
from random import uniform,random
from visual.controls import *

display(center=(1,0,0),background=(1,1,1), autoscale=False, range=(2,2,2),        
        width=600, height=600,  forward=(-.4,-.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)
n       = 2
dt      = 1./8.
pegas   = []
kotak   = []
box(pos=(-1,0,0), width=2, height=2, length= 2, color=color.black)
box(pos=(0,-.36,0), width=2, height=.2, length= 5, color=color.black,opacity=.3)
for i in arange(n):
    spring  = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
    pegas.append(spring)
    ko      = box(pos=(0,0,0), width=.5, height=.5, length= .5, color= color.green)
    kotak.append(ko)

k0      = 1.
k1      = 1.

m0      = 1.
m1      = 1.

l00     = 1.
l01     = 1.

l0      = 1.
l1      = 1.1

x0      = l0
x1      = l0+l1

v0      = 0.
v1      = 0.


y = 1.
print x1
def updatePegas():
    global x0,x1
    kotak[0].x      = x0
    pegas[0].length = l0
    
    kotak[1].x      = x1
    pegas[1].x      = l0
    pegas[1].length = l1

    

def proses():
    global l0,v0,x0,l1,v1,x1
    #untuk m0
    dx0     = l0-l00
    f0     = -k0*dx0
    dx1     = l1-l01
    f1     = -k1*dx1
    a0      = (f0-f1)/m0
    v0      += a0*dt
    x0      += v0*dt
    l0      = x0
    #untuk m1
    a1      = f1/m1
    v1      += a1*dt
    x1      += v1*dt
    l1      = x1-x0
    
    
    updatePegas()
    
while 1:
    rate (39)
    y   += .1
    proses()

    
    
    
    

    

.

Normal Mode

 Planned to write coupled oscillator, using double spring.

 For some reason, Vpython refused to set length of helix with some value from array or list or any value derived from them. It only accept plain number or  number from simple variable (like a = 3. ). It hugely messed the whole script as I used l = [] for helix length.

 So, rewriting the code, create l0 and l1 manually.

 Didn't have energy to code the rest. So, at the moment, just call it "normal mode" coupled oscillator, heheh...

#code

from visual import *
from random import uniform,random
from visual.controls import *

display(center=(0,0,0),background=(1,1,1), autoscale=False, range=(5,5,3),        
        width=600, height=600,  forward=(-1.4,-1.3,-1)) #arah kamera

distant_light(direction=(1,1,1), color=color.red)
n       = 2
dt      = 1./8.
pegas   = []
kotak   = []
w       = ones(n)
w       /= 2.
for i in arange(n):
    spring  = helix(pos=(0,0,0), axis=(5,0,0), radius=0.2, color=color.red, length=1.)
    pegas.append(spring)
    ko      = box(pos=(0,0,0), width=.5, height=.5, length= .5, color= color.green)
    kotak.append(ko)

k0      = 1.
m0      = 1.

l00     = 1.
l01     = 1.

l0      = 1.1
l1      = 1.

x0      = l0+w[0]/2

v0      = 0.
v1      = 0.


y = 1.
print pegas[0].length
def updatePegas():
    global x0
    kotak[0].x      = x0
    pegas[0].length = l0
    
    x1              = l0+w[0]+l1+w[1]/2
    kotak[1].x      = x1
    pegas[1].x      = l0+w[0]

    

def proses():
    global l0,v0,x0
    dx      = l0-l00
    f0      = -k0*dx
    a0      = f0/m0
    v0      += a0*dt
    l0      += v0*dt
    x0      = l0+w[0]/2

    
    updatePegas()
    
while 1:
    rate (19)
    y   += .1
    proses()

    
    
    
    

    

.