python - 龙格库塔逼近贝塞尔函数，二阶微分方程

Question

我正在尝试用 Runge-Kutta 逼近 Bessel 函数。我在这里使用 RK4，因为我不知道 ODE 系统的 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) 

也就是说，斜率必须按步长缩放才能获得变量的正确更新。