使用 RK4（硬编码）的 Lane-emden 解决方案

Question

我正在努力寻找以1/2 为间隔的值n=[0,6] 的Lane-Emden 方程的绘图解。我是 Python 新手，似乎不知道如何使用 RK4 来完成这项工作。请帮忙！

当前进度。

TypeError: unsupported operand type(s) for Pow: 'int' and 'list' on line 37 in main.py

我在定义为r2、r3、r4 和k2、k3、k4 的方程中添加后才出现错误。

import numpy as np
import matplotlib.pyplot as plt

n = [0,1,2,3,4,5,6,7,8,9,10,11,12,13]
theta0 = 1
phi0 = 0
step = 0.01
xi0 = 0
xi_max = 100

theta = theta0
phi = phi0
xi = xi0 + step

Theta = [[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
Phi = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
Xi = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]

for i in n:
    Theta[i].append(theta)
    Phi[i].append(phi)
    Xi[i].append(xi)
 
def dTheta_dXi(phi,xi): #r1
    return -phi/xi**2
    
def r2(phi,xi):
    return dTheta_dXi(phi+step,xi+step*dTheta_dXi(phi,xi))

def r3(phi,xi):
    return dTheta_dXi(phi+step,xi+step*r2(phi,xi))
    
def r4(phi,xi):
    return dTheta_dXi(phi+step,xi+step*r3(phi,xi))

def dPhi_dXi(theta,xi,n): #k1
    return theta**(n)*xi**2
    
def k2(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*dPhi_dXi(theta,xi,n),n)

def k3(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*k2(theta,xi,n),n)
    
def k4(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*k3(theta,xi,n),n)
    
for i in n: 
    while xi < xi_max:
        if theta < 0:
            break
        dTheta = (step/6)*(dTheta_dXi(phi,xi)+2*r2(phi,xi)+2*r3(phi,xi)+r4(phi,xi))
        dPhi = (step/6)*(dPhi_dXi(theta,xi,i/2.)+2*k2(theta,xi,n)+2*k3(theta,xi,n)+k4(theta,xi,n))
        theta += dTheta
        phi += dPhi
        xi += step
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)
    print i/2., round(xi,4), round(dTheta_dXi(phi,xi),4), round(xi/3./dTheta_dXi(phi,xi),4), round(1./(4*np.pi*(i/2.+1))/dTheta_dXi(phi,xi)**2,4)
    theta = theta0
    phi = phi0
    xi = xi0 + step

Answer 1

将 RK4 用于耦合系统

如果将一阶系统理解为处理向量状态的向量值系统，那么 RK4 的标量版本

k1 = f(x,y)
k2 = f(x+0.5*h, y+0.5*h*k1)
k3 = f(x+0.5*h, y+0.5*h*k2)
k4 = f(x+h, y+h*k3)
x,y = x+h, y+h/6*(k1+2*k2+2*k3+k4)

也可以直接用于向量的情况。有时，以组件方式实现此功能似乎很有教育意义。虽然在数学课本中最好使用单字母变量名，可能带有下标或上标，但程序代码中的变量通常是多字母的。因此，使用 k2_Theta 和 k2_Phi 代替 r2 和 k2 会更具描述性。

然后，用于评估k3 组件的状态具有theta+0.5*step*k2_Theta 和phi+0.5*step*k2_Phi 参数变得相当直观。

k2_Xi等对于自变量总是1，所以第三阶段的值就是xi+0.5*step。

实施细节 RK4

k1 等的值在步骤内部是固定的，并且是导数函数评估的结果。将它们本身声明为函数绝对没有意义。也就是说，专门针对这种情况的 RK4 步骤变成了只是

def RK4_update(theta, phi, xi, step, n):
    k1_Theta = dTheta_dXi(phi, xi)
    k1_Phi = dPhi_dXi(theta, xi, n)
    k2_Theta = dTheta_dXi(phi+0.5*step*k1_Phi, xi+0.5*step)
    k2_Phi = dPhi_dXi(theta+0.5*step*k1_Theta, xi+0.5*step, n)
    k3_Theta = dTheta_dXi(phi+0.5*step*k2_Phi, xi+0.5*step)
    k3_Phi = dPhi_dXi(theta+0.5*step*k2_Theta, xi+0.5*step, n)
    k4_Theta = dTheta_dXi(phi+step*k3_Phi, xi+step)
    k4_Phi = dPhi_dXi(theta+step*k3_Theta, xi+step, n)
    dTheta = (step/6)*(k1_Theta+2*k2_Theta+2*k3_Theta+k4_Theta)
    dPhi = (step/6)*(k1_Phi+2*k2_Phi+2*k3_Phi+k4_Phi)
    return dTheta, dPhi

关于 Lane-Emden 方程的奇异性

对于存在于xi=0 的解决方案，至少需要phi ~ xi^k 和k>=2。这使得theta 几乎是恒定的，这反过来又导致整合phi = theta0^n*xi^3/3，然后在另一个等式中给出theta = theta0 - theta0^n*xi^2/6。这允许在不使用数值方法的情况下迈出远离奇点的第一步。

xi = step 
theta, phi = theta0 - theta0**n*xi**2/6, theta0**n*xi**3/3
Xi[i] = [0, xi]
Theta[i] = [theta0, theta]
Phi[i] = [0, phi]

那么主循环可以写成

for i in range(N): 
    n = i/2
    xi = step 
    theta, phi = theta0 - theta0**n*xi**2/6, theta0**n*xi**3/3
    Xi[i] = [0, xi]
    Theta[i] = [theta0, theta]
    Phi[i] = [0, phi]
    while xi < xi_max:
        if theta < 0:
            break
        dTheta, dPhi = RK4_update(theta,phi,xi,step,n)
        theta += dTheta
        phi += dPhi
        xi += step
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)

然后绘制

for i in range(N):
    plt.plot(Xi[i],Theta[i], label=f"n={i/2}")
plt.grid(); plt.legend(); plt.show()

结果

使用的技巧：为避免负值的有理幂，请将 theta**n 替换为 theta*abs(theta)**(n-1) 或类似的延续。

---

旧内容

您应该再次探索更新在哪里进行。 xi 是自变量，因此只得到更新 0.5*step 和 step，theta 的更新使用导数 dTheta_dXi 和类似地 phi 使用斜率更新 dPhi_dXi

def r2(phi,xi):
    return dTheta_dXi(phi+0.5*step*dPhi_dXi(theta,xi,n),xi+0.5*step)

def k2(theta,xi,n):
    return dPhi_dXi(theta+0.5*step*dTheta_dXi(phi,xi),xi+0.5*step,n)

def r3(phi,xi):
    return dTheta_dXi(phi+0.5*step*k2(theta,xi,n),xi+0.5*step)

等等

现在可以看到，由于等式的耦合性质，您需要theta 和phi 作为任何地方的参数。此外，即使最终可行，您最终也会多次计算许多值，而在一个循环中组装所有内容只需要一次计算。