Неизвестная математическая ошибка в уравнении симпы (python)

Question

Я пытаюсь отобразить функции j0, j1 и j10 в диапазоне r (0,20), преобразовав их в пустой формат с помощью lambdify.Я использовал следующий код:

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}), однако после значительного пересмотра и переработки кода, который я нахожуЯ не могу получить нужные формулы для j0, j1 и j10 без использования .subs.Я думаю, что эта ошибка имеет эффект зацепки, так как я получаю сообщение об ошибке «цитирую неверный синтаксис», когда пытаюсь подставить формулу для j10 в следующее уравнение (которое должно идти в 0):

(r**2)*sym.diff(m,r,2) + (2*r)*sym.diff(m,r) + (r**2 - 10*(10+1))*m

Где m - это версия n10 с клочками.

Любая помощь будет принята с благодарностью.

Answer 1

Ваша проблема вызвана делением на ноль, с которым численно трудно справиться, даже если предел 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)$")