Search code examples
pythondifferential-equationsrunge-kutta

numpy.float64 error in my Runge-Kutta 4th-order method for a series of differential equations

This is a 4th order Runge-Kutta method I've made to eventually graph some differential equations.

The goal is to create a 4 by 100,000x.1 array that gives me the value of x, y, dx, dy at every point in the timestep, so that I can graph any equation with those 4 parameters. 

#Assumptions
x0, y0 = -.250, .433
x1, y1 = -.250,-.433
x2, y2 =  .500, .000
R      =  .2
C      =  .5
d      =  .25

#Imports
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as intgr 
import math

#ag = [[   x0,  y0], [   x1,   y1], [  x2,  y2]]
mag = [[-.250,.433], [-.250,-.433], [.500,.000]]

def der( xin, t ):
    mag = [[-.250,.433],[-.250,-.433],[.500,.000]]
    x   = xin[0]
    y   = xin[1]
    vx  = xin[2]
    vy  = xin[3]
    dx  = vx
    dy  = vy
    vx2 = 0
    vy2 = 0
    vx1 = -R * vx - C * x
    vy1 = -R * vy - C * y

    for i in range( mag.__len__() - 1 ):
        vx2 = vx2 + (   ( mag[i][0] - x )
                    / ( ( mag[i][0] - x )**2
                      + ( mag[i][1] - y )**2
                      +                 d**2
                        )**1.5
                      )
        vy2 = vy2 + (   ( mag[i][1] - y )
                    / ( ( mag[i][0] - x )**2
                      + ( mag[i][1] - y )**2
                      +                 d**2
                         )**1.5
                      )
        vx3 = vx1 + vx2
        vy3 = vy1 + vy2

        array = [dx,dy,vx3,vy3]
        return array

dt    = .1
t     = np.arange( 0, 100000, dt )
xzero = [.2, .2, 0, 0]

def RK4( func, xzero, t ):
    rows     = xzero.__len__()
    columns  =     t.__len__()
    x        = np.zeros( ( rows, columns ) )
    x_t      = 0
    ind      = 0
    x[:,ind] = xzero
    dt       = t[1] - t[0] 

    for time in t[0:len( t ) - 1]:
        ind      = ind + 1
        K1       = dt * func( x[:,ind-1],           time )
        K2       = dt * func( x[:,ind-1] + .5 * K1, time + .5 * dt )
        K3       = dt * func( x[:,ind-1] + .5 * K2, time + .5 * dt )
        K4       = dt * func( x[:,ind-1] +      K3, time +      dt )
        x[:,ind] = x[:,ind-1] + ( 1.0 / 6.0 ) * (     K1
                                                + 2 * K2
                                                + 2 * K3
                                                +     K4
                                                  )
    return x

print( RK4( func = der, xzero = xzero, t = t ) )

Produces a numpy float 64 error

I'm not exactly sure why but some variable in my code isn't being interpreted as a number?

Thanks for the help in advance and let me know if I should provide more code or a larger context.

Solution

  • The error message:

    You are trying to multiply a floating point number with an instance of a list.

    This kind of operation is actually well defined for integers, where you get the concatenation of multiple copies of the input list ( Given a = [1, 2, 3]; print( 2*a ) returns [1, 2, 3, 1, 2, 3] ). Thus the error message.

    Solution:

    You will want to use numpy consistently and especially the vector arithmetic that its array object provides.

    As a first point, the return of the ODE function from RK4() should be rather articulated as: 

    return np.array( [dx, dy, vx3, vy3] )