CODEWITHSUNDEEP
pythonnumpymatplotlibsympy
yerman
yerman

Reputation: 287

Unknown mathematical error in sympy equation (python)

I am attempting to plot the functions j0,j1 & j10 in the range r(0,20) by converting them to numpy format using lambdify. I used the following code:

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)$")

I got this output:

enter image description here

Which seems at least partially correct, however I also got this error message which I've never seen before:

/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
  """

I suspect this has something to do with using .subs({t:0}) however after much revising and reworking the code I find I'm unable to get the formulas I want for j0, j1 and j10 without using .subs. I think this error has a knock-on effect as I get an error quoting "incorrect syntax" when I try to substitute the formula for j10 into the following equation (which is supposed to go to 0):

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

Where m is the numpy version of j10.

Any help would be much appreciated.

Upvotes: 1

Views: 138

Answers (1)

Thomas Kühn
Thomas Kühn

Reputation: 9840

Your problem is caused by a division by zero, which is numerically hard to deal with, even though the limit of r->0 may be finite. I would have two (slightly different) solutions to the problem.

1) Replace the problematic point with the mathematical exact result. In your example this would mean something along the lines (limit is the exact solution of the function for r->0 that you derive first on paper) :

myr = np.linspace(0,20,1000)
k_noerror = np.concatenate([[limit], k(myr[1:])])
plt.plot(myr,k_noerror,label="$j_{0}(r)$")

2) If you cannot calculate the limit yourself, you may be able to fix the problem by replacing your zero with a very small value, i.e.:

myr = np.linspace(0,20,1000)
myr[0] = 1e-3
plt.plot(myr,k(myr),label="$j_{0}(r)$")

Upvotes: 1

Related Questions