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

【问题标题】：Runge Kutta approximation Bessel function, second order differential equation龙格库塔逼近贝塞尔函数，二阶微分方程
【发布时间】：2014-03-20 04:49:58
【问题描述】：

我正在尝试使用 Runge-Kutta 逼近 Bessel 函数。我在这里使用 RK4，因为我不知道 RK2 对于 ODE 系统的方程是什么.....基本上它不起作用，或者至少它不是很好，我不知道为什么。

反正它是用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()

【问题讨论】：

您能否详细说明“工作不太好”？它会产生异常等吗？干杯。
它应该重现此页面上的情节（或至少非常相似的内容）mathworld.wolfram.com/BesselFunctionoftheFirstKind.html，但事实并非如此。它做了一些奇怪的事情。你能运行它吗？只要你有我要导入的 python 包，它就是完全独立的。
它没有？它产生什么？还是根本不起作用？
它适用于 m = 0 和 m=1，并产生以下（分别）imgur.com/U0J4y8U,aVZ7CMi#0 请注意，在这两种情况下它都不能“正确启动”，看看它有多大对于 m = 1。它对 m>1 不起作用，因为边界条件状态 J_m(0) = 0 和 Z_m(0) = 0 所以一切总是零，我不知道你应该如何解决这个问题.
@GeneralPancake 你能否提供一些它应该通过但没有通过的单元测试。这可能是传达失败的地点和方式的“最佳”方式。

标签： python runge-kutta

【解决方案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)

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

【讨论】：