绘制使用 RK4 方法绘制的微分方程的结果

Question

我有一个求解 Lane-Emden 方程的代码，它在求解后给出了正确的值，但是当我使用 matplotlib 绘制这些列表时，我得到的是空白图像而不是正确的图。

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()

正确的剧情应该是this。谢谢。

编辑执行下面建议的编辑后，我收到此错误输出。

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.

EDIT2

Answer 1

总结一下导致解决方案的原始 post 下的评论：

考虑在 for 循环后仅调用 plt.plot() 和 plt.show() 一次，并使用 plt.plot(Xi[0],Phi[0]) 而不是单个绘制 Xi 和 Phi 的各个子列表就像你在这里做的那样 plt.plot(xi, phi).

绘制使用 RK4 方法绘制的微分方程的结果

Plotting the results of the differential equation plotted suing RK4 method

matplotlib

runge-kutta