plot of "slightly complicated" functions

[imakn634]

I have the following troubles when trying to plot some "slightly complicated" functions. Is there any restriction for the functions for plot()?

from sympy import *
from spb import *

def f(x):
    if x < 0:
        return -x
    else:
        return x
x=Symbol('x')
# plot(abs(x)) works.  But, 
plot(f(x))

returns

TypeError: cannot determine truth value of Relational

def g(x):
    return nsolve(sin(t)-x, t, 0)
x=Symbol('x')
# plot(asin(x), (x, -1, 1)) works.  But, 
plot(g(x), (x, -1, 1))

returns ValueError: expected a one-dimensional and numerical function

[Davide-sd]

First example

TypeError: cannot determine truth value of Relational

That happens when executing f(x), with x being a complex symbol (a symbol with no assumption).

However, the module is able to evaluate python's functions, or lambda functions. For example:

def f(x):
    if x < 0:
        return -x
    else:
        return x

x = Symbol('x')
plot(f, (x, -2, 2), force_real_eval=True)

Note that I passed the function f to the plot command, not the application of to something, f(x). In this case, we also need to use force_real_eval=True. By default, spb evaluates functions using complex numbers, which means each x inside f will be a complex number. Since you are using comparisons inside f, for example x < 0, x needs to be a real number.

Let's dig a little deeper on how the module works:

• if f is a symbolic expression, it will generale a lambda function with lambdify. If f is already a python function, it will be evaluated as it is.
• it generates an array of coordinates. In the above case, something like x = np.linspace(-2, 2).

• it evaluates f(x) (with Numpy, if not specified otherwise). Note that the array x is passed into f. The above function is coded in such a way that x should be a single number. Let's see what really happens in our case:

  import numpy as np
  xx = np.linspace(-2, 2)
  f(xx)
  # ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

  The module is coded to handle this case. When this error is detected, it automatically switch to another evaluation strategy that loops over xx, passing only a single number to f. Hence, plot(f, (x, -2, 2), force_real_eval=True) produces the expected result. Alternatively, you could code the function slightly differently, using vectorization:

  # using numpy vectorization
  def f(x):
     y = x.copy()
     idx = x < 0
     y[idx] = -y[idx]
     return y
  
  x = Symbol('x')
  plot(f, (x, -2, 2), force_real_eval=True, n=10)

  This doesn't raise the aforementioned error, and the evaluation is faster. However, it only works if we are using basic Python operations and Numpy functions.

Second example

Here, things are a little bit different. You have a function, g, that uses sympy's nsolve. For this particular case, the initial guess to nsolve should be one number. But we have already seen that the plotting module passes an array to g:

xx = np.linspace(-1, 1)
g(xx)
# TypeError: object of type 'Symbol' has no len()

Looking at the traceback, you'll see that the error is raised inside of nsolve: since xx is an array, nsolve thinks it is an array of initial guesses and expects t to be a list of symbols of equal length to xx. In this case, the function g should be evaluated on one single number. We can force that behavior with the decorator np.vectorize:

import numpy as np
@np.vectorize
def g(x):
    return nsolve(sin(t)-x, t, 0)
t = Symbol('t') 
plot(g, (x, -1, 1))

[Davide-sd]

I should probably code to module to detect that TypeError in order to force it to switch the evaluation strategy as in the first example.

[imakn634]

Thank you for the detailed explanation. For the moment, the quick fix might be the following:

from sympy import *
from sympy.abc import *
from spb import *

def f(x):
    var('x', real = True)
    if x < 0:
        return -x
    else:
        return x

# instead of 
# plot(f(x), (x, -2, 2))
plot(f, (x, -2, 2), force_real_eval=True);

def g(x):
    return nsolve(sin(t)-x, t, 0)

# instead of
# plot(g(x), (x, -1, 1));
plot(g, (x, -1, 1));

I confirmed the code works in my environment. This fix will only work for the one variable functions, and will not for, for example, F(x, a, b).

The original plot() of sympy, without importing spb, fails to plot both f(x) and g(x), it is really the advantage of your spb. It would be nicer if we could use plot(f(x), (x, -2, 2))...

[Davide-sd]

This fix will only work for the one variable functions, and will not for, for example, F(x, a, b).

Can you come up with an example?

The original plot() of sympy, without importing spb, fails to plot both f(x) and g(x), it is really the advantage of your spb. It would be nicer if we could use plot(f(x), (x, -2, 2))...

Sorry, plot(f, (x, -2, 2)) is the best we can achieve. It's how Python works. If f is your custom python function, the moment you write f(x) you are evaluating the function and plotting the results of its evaluation. You need to understand that Python is not a mathematical language, but a programmatic one.

[imakn634]

You need to understand that Python is not a mathematical language, but a programmatic one.

I understand. Thank you for pointing it out.

Can you come up with an example?

I just wanted to say:

from sympy import *
from sympy.abc import *
from spb import *
# Both work. 
plot(sin(x), (x, 0, 1));
plot(sin, (x, 0, 1));

def F(x, a):
    return sin(x)+a

# This is normal and works.
plot(F(x, 1), (x, 0, 1));
# This not.
plot(F, (x, 0, 1));

[Davide-sd]

Of course it doesn't work: what's a supposed to be? What kind of plot would you like to get? Here is a couple of examples of what you can do with that function.

We can let a be a number. Then, creates multiple functions by changing that number:

a, x = var("a x")

def F(p1, p2):
    return sin(p1) + p2

functions = [F(x, p) for p in range(3)]
plot(*functions, (x, 0, 1))

In the second example, we can let a be a parameter and create a widget-like plot (an interactive application). Note the slider in the following plot. This only works if you are using Jupyter and you have installed the complete plotting module:

%matplotlib widget
plot(F(x, a), (x, -1, 1), params={a: (1, -3, 3)}, ylim=(-5, 5))

If you are not using Jupyter, you can try the following which creates an interactive application and load it on a new browser page:

plot(F(x, a), (x, -1, 1), params={a: (1, -3, 3)}, ylim=(-5, 5), imodule="panel", servable=True)

[imakn634]

Thank you for many example solutions. I just wanted to say the following:

if I want to plot the numerical solution $t$ of, say, $y(t, x) \equiv \sin t - x = 0$, the code is

def g(x):
    return nsolve(y(t, x), t, 0)

# instead of
# plot(g(x), (x, -1, 1));
plot(g, (x, -1, 1));

Let us call it as "function name only method (without parentheses)" or simply "The Method".

The Method works. However, if we want to plot the numerical solution $t$ of $y(t, x, a) = 0$, the function will be

def G(x, a):
    return nsolve(y(t, x, a), t, 0)

In this case,

plot(G(x, 1), (x, -1, 1));

will fail because of the reason you explained.

The Method suggests (function name only!)

plot(G, (x, -1, 1));

Of course it doesn't work. That's all.

The Method works only for the single variable functions. So I think the change of the function might work.

def GG(x):
    global a
    return nsolve(y(t, x, a), t, 0)

a=1
plot(GG, (x, -1, 1));

I'm just trying to plot Trajectory of a projectile with air resistance.

• Projectile motion - Wikipedia

The vertical position is (something like),

$$ y(t, \theta, \beta)= \frac{1-e^{-\beta t}}{\beta} \sin\theta + \frac{1 -\beta t – e^{-\beta t}}{\beta^2} $$

Time to reach ground $t_d$ is the solution of

$$y(t_d, \theta, \beta) = 0,$$

which can be obtained by td = nsolve(y(t, theta, beta), t, ..)...

[Davide-sd]

How about letting python figure out what the parameters are? I don't see the benefits of over complicating things...

Here A*sin(w*t) is a mathematical function of 3 parameters. I can set two of the outside the python function:

x, t = symbols('x t') 
A = 2
w = 3
def g(x):
    return nsolve(A*sin(w*t)-x, t, 0)
plot(g, (x, -2, 2))

[imakn634]

Thanks. I agree with your suggestion. No need to declare "global a" in the definition of the funcion.

[imakn634]

If we could use plot(g(x, a), (x, -2, 2)), then plotting multiple expressions would be simply,

plot(g(x,1), g(x,2), ..., (x, -2, 2))

For the moment, because we can not plot(g(x, a), (x, -2, 2)), we will need to do

a = 1
p1 = plot(g, (x, -2, 2), show=False)
a = 2
p2 = plot(g, (x, -2, 2), show=False)
...
(p1+p2+...).show()

[Davide-sd]

You can already do that:

x, t = symbols('x t') 
w = 4
def g(A, w):
    return lambda x: nsolve(A*sin(w*t)-x, t, 0)

plot(g(2, w), g(3, w), (x, -2, 2))

First, let's understand what is going on. g(A, w) is a Python function, with two parameters, A, w. Once it is invoked, it returns a new lambda function which uses those parameters. This lambda function requires only one argument, x, so the plotting module is able to evaluate it.

Let's look at the picture: as you can see, only one function was plotted. That's a bug: for reasons to long to explain, the plotting module is able to deal with only one lambda-python function per call. The other functions are discarded. I don't intend to fix this bug as it requires way too much time.

You can plot them like this:

A = [2, 3, 4]
plots = [plot(g(a, w), (x, -2, 2), show=False) for a in A]
sum(plots).show()

Alternatively, you can use the graphics module:

graphics(
    line(g(2, w), (x, -2, 2)),
    line(g(3, w), (x, -2, 2)),
)

[imakn634]

Thank you. The information is much helpful!

[Davide-sd]

New release is out, fixing the problems shown on this issue.

from sympy import *
from spb import *

def f(x):
    if x < 0:
        return -x
    else:
        return x

x = Symbol('x')
plot(f, (x, -2, 2))

def g(x):
    return nsolve(sin(t)-x, t, 0)
x, t = symbols('x, t') 
plot(g, (x, -1, 1))

w = 4
def h(A, w):
    return lambda x: nsolve(A*sin(w*t)-x, t, 0)

A = [2, 3, 4]
plots = [plot(h(a, w), (x, -2, 2), show=False) for a in A]
sum(plots).show()

I'm going to close this issue.