3

I'm trying to approximate the Bessel function with Runge-Kutta. I'm using RK4 here cause I don't know what the equations are for RK2 for a system of ODEs.....Basically it's not working, or at least it's not working very well, I don't know why.

Anyway it's written in python. Here is code:

from __future__ import division 
from pylab import* 

#This literally just makes a matrix 
def listoflists(rows,cols): 
    return [[0]*cols for i in range(rows)] 

def f(x,Jm,Zm,m): 
    def rm(x): 
        return (m**2 - x**2)/(x**2) 
    def q(x):
        return 1/x
    return rm(x)*Jm-q(x)*Zm

n = 100 #No Idea what to set this to really 
m = 1 #Bessel function order; computes up to m (i.e. 0,1,2,...m)
interval = [.01, 10] #Interval in x 
dx = (interval[1]-interval[0])/n #Step size 
x = zeros(n+1) 
Z = listoflists(m+1,n+1) #Matrix: Rows are Function order, Columns are integration step (i.e.                 function value at xn)
J = listoflists(m+1,n+1) 

x[0] = interval[0] 
x[n] = interval[1]

#This reproduces all the Runge-Kutta relations if you read 'i' as 'm' and 'j' as 'n'
for i in range(m+1): 
    #Initial Conditions, i is m 
    if i == 0: 
        J[i][0] = 1 
        Z[i][0] = 0
    if i == 1: 
        J[i][0] = 0 
        Z[i][0] = 1/2 
    #Generate each Bessel function, j is n 
    for j in range(n):  
        x[j] = x[0] + j*dx 
        
        K1 = Z[i][j]
        L1 = f(x[j],J[i][j],Z[i][j],i)
        K2 = Z[i][j] + L1/2 
        L2 = f(x[j] + dx/2, J[i][j]+K1/2,Z[i][j]+L1/2,i) 
        K3 = Z[i][j] + L2/2 
        L3 = f(x[j] +dx/2, J[i][j] + K2/2, Z[i][j] + L2/2,i)
        K4 = Z[i][j] + L3 
        L4 = f(x[j]+dx,J[i][j]+K3, Z[i][j]+L3,i) 

        J[i][j+1] = J[i][j]+(dx/6)*(K1+2*K2+2*K3+K4) 
        Z[i][j+1] = Z[i][j]+(dx/6)*(L1+2*L2+2*L3+L4) 

plot(x,J[0][:])
show()
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
GeneralPancake
  • 535
  • 2
  • 6
  • 15
  • 1
    Can you elaborate on "not working very well"? Does it produce an Exception, etc? Cheers. – Ffisegydd Mar 20 '14 at 07:49
  • It should reproduce the plot on this page (or at least something very similar) http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html and it does not. It does something weird. Can you run it? It's completely standalone as long as you have the python packages I'm importing. – GeneralPancake Mar 20 '14 at 13:40
  • 1
    And it doesn't? What does it produce? Or doesn't it work at all? – Ffisegydd Mar 20 '14 at 13:41
  • It works for m = 0 and m=1, and produces the following (respectively) http://imgur.com/U0J4y8U,aVZ7CMi#0 Notice that it doesn't "start right" in either case, and look at how huge it gets for m = 1. It flatout does not work for m>1 because the boundary conditions state J_m(0) = 0 and Z_m(0) = 0 so everything is always zero and I don't know how you are supposed to fix that. – GeneralPancake Mar 20 '14 at 13:44
  • @GeneralPancake Can you provide some unit tests that it should pass but do not. That is probably the "best" way to communicate the where and how it's failing. – wheaties Mar 20 '14 at 13:49
  • What do you mean by unit tests? – GeneralPancake Mar 20 '14 at 13:52

1 Answers1

0

(To close this old question with at least a partial answer.) The immediate error is that the RK4 method is incorrectly implemented. For instance instead of

K2 = Z[i][j] + L1/2 
L2 = f(x[j] + dx/2, J[i][j]+K1/2, Z[i][j]+L1/2, i)

it should be

K2 = Z[i][j] + L1*dx/2 
L2 = f(x[j] + dx/2, J[i][j]+K1*dx/2, Z[i][j]+L1*dx/2, i)

that is, the slopes have to be scaled by the step size to get the correct update for the variables.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51