matplotlib - 绘制使用 RK4 方法绘制的微分方程的结果
问题描述
我有一个求解 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()
正确的情节应该是这样的。谢谢。
编辑 在执行下面建议的编辑后,我得到了这个错误输出。
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.
解决方案
总结导致解决方案的原始帖子下的评论:
考虑在 for 循环之后只调用一次plt.plot()
andplt.show()
并绘制Xi
和Phi
with的各个子列表,而plt.plot(Xi[0],Phi[0])
不是像你在此处所做的那样单点plt.plot(xi, phi)
。
推荐阅读
- swagger - Nelmio Docs 正在返回 401 Undocumented
- autocad - 有没有办法从 dxf 文件手动计算实体总长度?
- arrays - 在汇编中将数组打印为双字
- php - 我想使 MySQL 查询适应 postgreSQL
- java - 在握手 SSL 错误期间收到 close_notify?
- virtualbox - DBCA 权限被拒绝
- bash - 后缀从 0 增加到 1 后设置退出状态 1
- python - 根据表单中的值重定向到视图(Django)
- php - 在php中减去mysql数据库中两个不相关表的两列的总和
- reactjs - 我使用 mobx 5.15.7 看到奇怪的行为