# -*- coding: utf-8 -*-
"""
Created on Sun Dec 09 18:28:40 2012

@author: ваха
"""
    
from numpy import *
import pylab
from matplotlib import mlab
from pylab import *

a=2
b=3
m=1
F=13
t0=0
tmax=10
epsilon=0.01
y0=[0,0,math.pi/3,0,0,0]
h=0.1

def rkm45(f,t0,tmax,y0,epsilon,h,i):
    t=t0
    n=int((tmax-t0)/h)
    s1=(6,5)
    k=zeros(s1)
    s2=(6,n+1)
    y=zeros(s2)
    yb=zeros(s2)
    yp=zeros(6)
    y[:,0]=y0[:]
    k[:,0]=f(t,y[:,i])
    k[:,1]=f(t+h/3,y[:,i]+h/3*k[:,0])
    k[:,2]=f(t+h/3,y[:,i]+h/6*k[:,0]+h/6*k[:,1])
    k[:,3]=f(t+h/2,y[:,i]+h/2*k[:,0]+h/8*3*k[:,1])
    k[:,4]=f(t+h/2,y[:,i]+h/2*k[:,0]-3/2*h*k[:,2]-2*h*k[:,3]) 
    yp=y[:,i]+h/2*(k[:,0]+3*k[:,2]+4*k[:,3])
    y[:,i+1]=y[:,i]+h/6*(k[:,0]+4*k[:,3]+k[:,4])
    t=t+h
    R=0.2*math.fabs(y[2,i+1]-yp[2])       
    if R>=epsilon:
        i=i+1
        h=h/2.0
        rkm45(func,t0,tmax,y0,epsilon,h,i)
    else:
        if R<=epsilon/64.0:
            i=i+1            
            h=2*h
            rkm45(func,t0,tmax,y0,epsilon,h,i)
        else: 
            if i<3:
                while(i<3):
                    t=t+h                
                    k[:,0]=f(t,y[:,i])
                    k[:,1]=f(t+h/3,y[:,i]+h/3*k[:,0])
                    k[:,2]=f(t+h/3,y[:,i]+h/6*k[:,0]+h/6*k[:,1])
                    k[:,3]=f(t+h/2,y[:,i]+h/2*k[:,0]+h/8*3*k[:,1])
                    k[:,4]=f(t+h/2,y[:,i]+h/2*k[:,0]-3/2*h*k[:,2]-2*h*k[:,3]) 
                    y[:,i+1]=y[:,i]+h/6*(k[:,0]+4*k[:,3]+k[:,4])
                    i=i+1
            while i<n-1:
                yb[:,i+1]=y[:,i]+h/24*[55*f(t,y[:,i])-59*f(t-h,y[:,i-1])+37*f(t-2*h,y[:,i-2])-9*f(t-3*h,y[:,i-3])]
                y[:,i+1]=y[:,i]+h/24*[9*f(t+h,yb[:,i+1])+19*f(t,y[:,i])-5*f(t-h,y[:,i-1])+f(t-2*h,y[:,i-2])] 
                i=i+1
    print y
    return y

def func(t,y):
    nu=zeros(6)
    nu=y
    dk1=nu[1]
    dk3=nu[3]
    dk5=nu[5]
    dk2=F*math.cos(nu[4])/m
    dk4=F*math.sin(nu[4])/m
    dk6=(3*a*F)/((a^2+b^2)*m)
    res=zeros(5)
    res=[dk1,dk2,dk3,dk4,dk5,dk6]
    return res
    


nice=rkm45(func,t0,tmax,y0,epsilon,h,0)



def graf1():
    x = arange(len(nice[0,:]))
    y = nice[0,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf2():
    x = arange(len(nice[0,:]))
    y = nice[1,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf3():
    x = arange(len(nice[0,:]))
    y = nice[2,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf4():
    x = arange(len(nice[0,:]))
    y = nice[3,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf5():
    x = arange(len(nice[0,:]))
    y = nice[4,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()

def graf6():
    x = arange(len(nice[0,:]))
    y = nice[5,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
    

l1 = [7,8,9]
l2 = [6,6,6]
l3 = [(e1-e2)*2 for e1,e2 in zip(l1,l2)]
l3
Out[4]: [2, 4, 6]

# -*- coding: utf-8 -*-
"""
Created on Sun Dec 09 18:28:40 2012
@author: ваха
"""
    
from numpy import *
import pylab
from matplotlib import mlab
from pylab import *
a=2
b=3
m=1
F=13
t0=0
tmax=10
epsilon=0.01
y0=[0,0,math.pi/3,0,0,0]
h=0.1
def adams(f,t0,tmax,y0,epsilon,h,i):
    t=t0
    n=int((tmax-t0)/h)
    s1=(6,5)
    k=zeros(s1)
    s2=(6,n+1)
    y=zeros(s2)
    yb=zeros(s2)
    yp=zeros(6)
    y[:,0]=y0[:]
    k[:,0]=f(t,y[:,i])
    k[:,1]=f(t+h/3,y[:,i]+h/3*k[:,0])
    k[:,2]=f(t+h/3,y[:,i]+h/6*k[:,0]+h/6*k[:,1])
    k[:,3]=f(t+h/2,y[:,i]+h/2*k[:,0]+h/8*3*k[:,1])
    k[:,4]=f(t+h/2,y[:,i]+h/2*k[:,0]-3/2*h*k[:,2]-2*h*k[:,3]) 
    yp=y[:,i]+h/2*(k[:,0]+3*k[:,2]+4*k[:,3])
    y[:,i+1]=y[:,i]+h/6*(k[:,0]+4*k[:,3]+k[:,4])
    t=t+h
    R=0.2*math.fabs(y[2,i+1]-yp[2])       
    if R>=epsilon:
        i=i+1
        h=h/2.0
        rkm45(func,t0,tmax,y0,epsilon,h,i)
    else:
        if R<=epsilon/64.0:
            i=i+1            
            h=2*h
            rkm45(func,t0,tmax,y0,epsilon,h,i)
        else: 
            if i<3:
                while(i<3):
                    t=t+h                
                    k[:,0]=f(t,y[:,i])
                    k[:,1]=f(t+h/3,y[:,i]+h/3*k[:,0])
                    k[:,2]=f(t+h/3,y[:,i]+h/6*k[:,0]+h/6*k[:,1])
                    k[:,3]=f(t+h/2,y[:,i]+h/2*k[:,0]+h/8*3*k[:,1])
                    k[:,4]=f(t+h/2,y[:,i]+h/2*k[:,0]-3/2*h*k[:,2]-2*h*k[:,3]) 
                    y[:,i+1]=y[:,i]+h/6*(k[:,0]+4*k[:,3]+k[:,4])
                    i=i+1
            while i<n-1:
                l0=y[:,i]
                l1=f(t,y[:,i])
                l2=f(t-h,y[:,i-1])
                l3=f(t-2*h,y[:,i-2])
                l4=f(t-3*h,y[:,i-3])
                d0=y[:,i]
                d1=f(t+h,yb[:,i+1])
                d2=f(t,y[:,i])
                d3=f(t-h,y[:,i-1])
                d4=f(t-2*h,y[:,i-2])
                yb[:,i+1]=[e0+h/24*(55.0*e1-59.0*e2+37.0*e3-9.0*e4) for e0,e1,e2,e3,e4 in zip(l0,l1,l2,l3,l4)]
                y[:,i+1]=[d0+h/24*(9*d1+19*d2-5*d3+d4) for d0,d1,d2,d3,d4 in zip(d0,d1,d2,d3,d4)] 
                i=i+1
    print y
    return y
def func(t,y):
    nu=zeros(6)
    nu=y
    dk1=nu[1]
    dk3=nu[3]
    dk5=nu[5]
    dk2=F*math.cos(nu[4])/m
    dk4=F*math.sin(nu[4])/m
    dk6=(3*a*F)/((a^2+b^2)*m)
    res=zeros(5)
    res=[dk1,dk2,dk3,dk4,dk5,dk6]
    return res
    
nice=adams(func,t0,tmax,y0,epsilon,h,0)
def graf1():
    x = arange(len(nice[0,:]))
    y = nice[0,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
def graf2():
    x = arange(len(nice[0,:]))
    y = nice[1,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
def graf3():
    x = arange(len(nice[0,:]))
    y = nice[2,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
def graf4():
    x = arange(len(nice[0,:]))
    y = nice[3,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
def graf5():
    x = arange(len(nice[0,:]))
    y = nice[4,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
def graf6():
    x = arange(len(nice[0,:]))
    y = nice[5,:]
    area = 1
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
    
graf1()
graf2()
graf3()
graf4()
graf5()
graf6()