python - sympy 方程中的未知数学错误(python)
问题描述
我试图通过使用lambdify将它们转换为numpy格式来绘制r(0,20)范围内的函数j0,j1和j10。我使用了以下代码:
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
from ipywidgets.widgets import interact
sym.init_printing(use_latex="mathjax")
x, y, z, t = sym.symbols('x y z t')
r = sym.symbols("r", positive=True)
j0 = (sym.diff(((sym.cos(sym.sqrt(r**2-2*r*t)))/r),t)).subs({t:0})
j1 = (sym.diff(((sym.cos(sym.sqrt(r**2-2*r*t)))/r),t,2)).subs({t:0})
j10 = (sym.diff(((sym.cos(sym.sqrt(r**2-2*r*t)))/r),t,11)).subs({t:0})
k = sym.lambdify(r,j0)
l = sym.lambdify(r,j1)
m = sym.lambdify(r,j10)
myr = np.linspace(0,20,1000)
plt.plot(myr,k(myr),label="$j_{0}(r)$")
plt.plot(myr,l(myr),label="$j_{1}(r)$")
plt.plot(myr,m(myr),label="$j_{10}(r)$")
plt.ylim(-1,1)
plt.legend()
plt.xlabel("r")
plt.ylabel("$j_{n}(r)$")
我得到了这个输出:
这似乎至少部分正确,但是我也收到了我以前从未见过的错误消息:
/anaconda3/lib/python3.6/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in true_divide
"""
/anaconda3/lib/python3.6/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in true_divide
"""
/anaconda3/lib/python3.6/site-packages/numpy/__init__.py:1: RuntimeWarning: divide by zero encountered in true_divide
"""
/anaconda3/lib/python3.6/site-packages/numpy/__init__.py:1: RuntimeWarning: invalid value encountered in true_divide
"""
/anaconda3/lib/python3.6/site-packages/numpy/__init__.py:1: RuntimeWarning: divide by zero encountered in true_divide
"""
我怀疑这与使用 .subs({t:0}) 有关,但是经过大量修改和修改代码后,我发现如果不使用 .subs,我无法获得我想要的 j0、j1 和 j10 的公式。我认为这个错误会产生连锁反应,因为当我尝试将 j10 的公式替换为以下等式(应该为 0)时,我得到一个错误引用“不正确的语法”:
(r**2)*sym.diff(m,r,2) + (2*r)*sym.diff(m,r) + (r**2 - 10*(10+1))*m
其中 m 是 j10 的 numpy 版本。
任何帮助将非常感激。
解决方案
您的问题是由除以零引起的,即使限制r->0
可能是有限的,这在数字上也很难处理。对于这个问题,我会有两个(略有不同)的解决方案。
1)用数学精确结果替换问题点。在您的示例中,这将意味着一些事情(是您首先在纸上得出limit
的函数的确切解决方案):r->0
myr = np.linspace(0,20,1000)
k_noerror = np.concatenate([[limit], k(myr[1:])])
plt.plot(myr,k_noerror,label="$j_{0}(r)$")
2)如果您自己无法计算限制,您可以通过将零替换为非常小的值来解决问题,即:
myr = np.linspace(0,20,1000)
myr[0] = 1e-3
plt.plot(myr,k(myr),label="$j_{0}(r)$")
推荐阅读
- javascript - 在javascript中将时间戳转换为特定的日期格式
- javascript - 如何有效地检查连续数字列表是否缺少任何元素
- c# - MySQL c# 连接字符串故障转移
- ios - 更改 ATS 设置会导致我的应用程序挂起 - Metro 构建但应用程序无法启动
- vba - 查找所有最早的日期并在消息框中显示它们
- keras - 如何在 lstm keras 中获得非负输出
- python - 熊猫比较两列写入另一列
- visual-studio - TFS 工作区错误
- reactjs - 无法使用 jsx 导入手动放置的节点模块
- cluster-computing - 使用 bing ClusterLayer (clientWidth) 时出错