RK4方法求解Lane-Emden方程后matplotlib绘图空白及报错求助

我编写了使用RK4方法求解Lane-Emden方程的Python代码，求解得到的数值结果正确，但使用matplotlib对输出的列表数据绘图时，仅生成空白图像，无法得到标准的Lane-Emden方程解曲线。按照他人建议修改代码后又出现维度不匹配报错。

原始代码 #

import numpy as np
import matplotlib.pyplot as plt

n = 14
theta_0 = 1
phi_0 = 0
h = 0.01
xi_0 = 0
xi_max = 100
   
theta = theta_0
phi = phi_0
xi = xi_0 + h
   
Theta = [[] for i in range(n)]
Phi = [[] for i in range(n)]
Xi = [[] for i in range(n)]
    
for i in range(n):
    Theta[i].append(theta)
    Phi[i].append(phi)
    Xi[i].append(xi)
     
def dThetadXi(phi,xi): #r1
    return -phi/xi**2
        
def r2(phi,xi):
    return dThetadXi(phi+h,xi+h*dThetadXi(phi,xi))
    
def r3(phi,xi):
    return dThetadXi(phi+h,xi+h*r2(phi,xi))
      
def r4(phi,xi):
    return dThetadXi(phi+h,xi+h*r3(phi,xi))
    
def dPhidXi(theta,xi,n): #k1
    return theta**(n)*(xi**2)
    
        
def k2(theta,xi,n):
    return dPhidXi(theta+h,xi+h*dPhidXi(theta,xi,n),n)
    
def k3(theta,xi,n):
    return dPhidXi(theta+h,xi+h*k2(theta,xi,n),n)
        
def k4(theta,xi,n):
    return dPhidXi(theta+h,xi+h*k3(theta,xi,n),n)
        
for i in range(n):  
    while xi < xi_max:
        if theta < 0:
            break
        dTheta = (step/6)*(dThetadXi(phi,xi)+2*r2(phi,xi)+2*r3(phi,xi)+r4(phi,xi))
        dPhi = (step/6)*(dPhidXi(theta,xi,i/2.)+2*k2(theta,xi,n)+2*k3(theta,xi,n)+k4(theta,xi,n))
        theta = theta+ dTheta
        phi = phi +dPhi
        xi = xi + h
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)
    print (i/2., round(xi,2), round(dThetadXi(phi,xi),2), round(xi/3./dThetadXi(phi,xi),2), round(1./(4*np.pi*(i/2.+1))/dThetadXi(phi,xi)**2,2))
    theta = theta_0
    phi = phi_0
    xi = xi_0 + h
    plt.plot(phi, xi)
    plt.show()

修改后报错信息 #

usr/local/lib/python3.7/dist-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order, subok=True)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: float() argument must be a string or a number, not 'list'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
<ipython-input-20-df5bae47b875> in <module>()
     62     phi = phi0
     63     xi = xi0 + step
---> 64     plt.plot(Phi, Xi)
     65 

7 frames
/usr/local/lib/python3.7/dist-packages/numpy/core/_asarray.py in asarray(a, dtype, order)
     81 
     82     """
---> 83     return array(a, dtype, copy=False, order=order)
     84 
     85 

ValueError: setting an array element with a sequence.

• 变量名不统一：原始代码中定义步长为h，但RK4更新逻辑里使用了未定义的step变量，会直接触发变量不存在报错
• 绘图数据维度错误：Theta、Phi、Xi都是嵌套列表，外层对应不同的多态指数分组，内层是每个分组对应的求解序列。修改代码时直接把整个嵌套列表传给plt.plot，matplotlib无法解析嵌套结构，就会触发类型转换报错
• 原始绘图逻辑错误：原始代码循环中只绘制了单次迭代的单个数值点，当然只会生成空白图像，且每次循环都调用plt.show()，无法把多个多态指数的曲线绘制在同一张图上
• RK4公式错误：原始代码中的r2_r4、k2k4的增量计算不符合RK4标准格式，会导致求解结果存在偏差

import numpy as np
import matplotlib.pyplot as plt

n_group = 14  # 计算0到7共7组多态指数
theta_0 = 1
phi_0 = 0
h = 0.01  # 统一使用h作为步长变量
xi_0 = 0
xi_max = 100
   
theta = theta_0
phi = phi_0
xi = xi_0 + h
   
Theta = [[] for i in range(n_group)]
Phi = [[] for i in range(n_group)]
Xi = [[] for i in range(n_group)]
    
for i in range(n_group):
    Theta[i].append(theta)
    Phi[i].append(phi)
    Xi[i].append(xi)
     
def dThetadXi(phi,xi): 
    if xi == 0:
        return 0  # 单独处理xi=0的情况避免除零错误
    return -phi/xi**2

# 修正为标准RK4计算格式
def rk4_theta(phi, xi, h):
    k1 = dThetadXi(phi, xi)
    k2 = dThetadXi(phi + k1 * h/2, xi + h/2)
    k3 = dThetadXi(phi + k2 * h/2, xi + h/2)
    k4 = dThetadXi(phi + k3 * h, xi + h)
    return (k1 + 2*k2 + 2*k3 + k4) * h /6

def dPhidXi(theta,xi,n): 
    return (theta**n) * (xi**2)

def rk4_phi(theta, xi, n, h):
    k1 = dPhidXi(theta, xi, n)
    k2 = dPhidXi(theta + k1 * h/2, xi + h/2, n)
    k3 = dPhidXi(theta + k2 * h/2, xi + h/2, n)
    k4 = dPhidXi(theta + k3 * h, xi + h, n)
    return (k1 + 2*k2 + 2*k3 + k4) * h /6
        
for i in range(n_group):
    n_poly = i/2  # 当前组对应的多态指数n
    while xi < xi_max:
        if theta < 0:
            break
        dTheta = rk4_theta(phi, xi, h)
        dPhi = rk4_phi(theta, xi, n_poly, h)
        theta += dTheta
        phi += dPhi
        xi += h
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)
    print(n_poly, round(xi,2), round(dThetadXi(phi,xi),2), round(xi/3./dThetadXi(phi,xi),2), round(1./(4*np.pi*(n_poly+1))/dThetadXi(phi,xi)**2,2))
    # 取当前组的序列绘图，添加图例
    plt.plot(Xi[i], Theta[i], label=f'n={n_poly:.1f}')
    # 重置变量计算下一组
    theta = theta_0
    phi = phi_0
    xi = xi_0 + h

# 所有曲线绘制完成后再设置显示参数
plt.xlabel(r'$\xi$')
plt.ylabel(r'$\theta_n(\xi)$')
plt.xlim(0, 10)
plt.ylim(0, 1.1)
plt.legend()
plt.grid(alpha=0.3)
plt.show()

运行上述代码即可得到符合预期的Lane-Emden方程解曲线，不会再出现空白图或者维度报错问题。

内容的提问来源于stack exchange，提问作者Isha Upadhyay