带RK4的二阶ODE

Question

我需要用RK4绘制一个显示2个变量的图形，到目前为止我已经完成了

from numpy import arange
from pylab import plot,xlabel,ylabel,show
Qger = 400
K = 20
T1 = 150
T2 = 60
N = 1000
h = (T2-T1)/N
rpoints = arange(6.0,8.0,h)
xpoints = []
x = 423
def df(s,t):
    dTdt = -Qger*t/(2*K) + 172.8/t
    return dTdt

for r in rpoints:
    xpoints.append(x)
    k1 = h*df(x,r)
    k2 = h*df(x+0.5*k1,r+0.5*h)
    k3 = h*df(x+0.5*k2,r+0.5*h)
    k4 = h*df(x+k3,r+h)
    x += (k1+2*k2+2*k3+k4)/6
pylab.plot(rpoints,xpoints)
pylab.xlabel("Raio")
pylab.ylabel("Temperatura")
pylab.show

但这是一个一阶颂歌的RK4，因为我不知道，也不集成 hand，但我不能这样做，也不能使用scipy，所以有人能给我解释一下如何集成这个函数或者如何使用RK4和二阶ODE。功能如下。在

这是函数，只有T和r是变量，其余的是0

Answer 1

你应该能够把上面的结果写成“半离散”形式，也就是说，只对r的偏导数表示dT/dT。如果你能找到一个数值或其他方法来近似等于dT/dT的项，即dT/dT=df（r，…）的RHS，那么显式的RK4就可以适用了。在

在这种方法中，时间步进法（RK4）只适用于温度对时间的一阶导数。在