用Runge Kutta法近似贝塞尔函数的二阶微分方程

Question

我正在尝试用龙格-库塔方法来近似贝塞尔函数。我这里用的是RK4，因为我不知道对于一组常微分方程，RK2的公式是什么……基本上，这个方法没有奏效，或者说效果不是很好，我也不知道为什么。

总之，这段代码是用Python写的。代码如下：

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()  

Answer 1

(为了给这个老问题一个至少部分的答案。) 直接的错误在于RK4方法的实现不正确。例如，应该用

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

而不是

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) 

也就是说，斜率需要乘以步长，才能正确更新变量。