[Math] Need (free) software to plot recursively defined functions.

math-softwarereference-request

I need software to plot

$$P_{n+1}(x) = \frac{n}{m}\left\{ P_{n-1}(x)-\frac{x^{n}}{n} \right\} $$

given that

$$P_0(x) = \exp\left(\displaystyle \frac{mx^2}{2}\right)\int_0^x \exp\left(-\displaystyle \frac{mt^2}{2}\right)dt$$

$$P_1(x) =\frac{1}{m}\left\{ \exp\left(\displaystyle \frac{mx^2}{2}\right)-1 \right\}$$

Any good reference?

Best Answer

It's not terribly pretty, but here's a start in sage:

p = []
G = Graphics()
xw = 4  # maximum/window "radius" for x
yw = 10 # maximum/window "radius" for y
var('x,m')
assume(x > 0)
assume(m > 0)
y = e^(m*x^2/2)
p.append(y*integrate(1/y,(x,0,x)))
G += plot(p[-1].subs(m=1),-xw,xw)
p.append((y-1)/m)
G += plot(p[-1].subs(m=1),-xw,xw,color='green')
for k in range(3):
    p.append(p[k]*(k+1)/m-x^(k+1)/m)
    G += plot(p[-1].subs(m=1),-xw,xw,color=(1,k/3,k/3))
G.show(ymin=-yw,ymax=yw,figsize=4)

enter image description here

Here's two versions of the above code, made into a function to play around with. This first version requires the $x$ and $y$ limits. It also shows you the functions it calculated.

def myplot(m,n,x_max,y_max):
    p = []
    G = Graphics()
    var('x')
    assume(x > 0)
    y = e^(m*x^2/2)
    p.append(y*integrate(1/y,(x,0,x)))
    G += plot(p[-1],-x_max,x_max)
    p.append((y-1)/m)
    G += plot(p[-1],-x_max,x_max,color='green')
    for k in range(n):
        p.append(p[k]*(k+1)/m-x^(k+1)/m)
        G += plot(p[-1],-x_max,x_max,color=(1,k/n,k/n))
    G.show(ymin=-y_max,ymax=y_max)
    return p

This second version does not require $y$ limits, but scales automatically in the $y$ dimension.

def my_plot(m,n,x_max):
    p = []
    G = Graphics()
    var('x')
    assume(x > 0)
    y = e^(m*x^2/2)
    p.append(y*integrate(1/y,(x,0,x)))
    G += plot(p[-1],-x_max,x_max)
    p.append((y-1)/m)
    G += plot(p[-1],-x_max,x_max,color='green')
    for k in range(n):
        z = p[k]*(k+1)/m-x^(k+1)/m
        p.append(z)
        G += plot(z,-x_max,x_max,color=(1,k/n,k/n))
    G.show()

If you can, try playing around with these and let me know what you find. There are a few potential problems. (1) is it calculating a formula that works for $x<0$? (2) I draw $P_0$ in blue, $P_1$ in green, and the iterates in successively lighter shades of red. However, I didn't seem to see as many shades of red as I expected. It should draw $n$ iterates, $n$ being the second function argument. (3) There are some runtime errors but it does produce a plot. For example:

myplot(-1,5,6,6)

enter image description here

Hope this helps for now! I'll try to check back within 24 hours.

Upon reflection/recollection, I noticed and was going to ask whether you are aware that $P_{n}(x)$ is the particular solution to the differential equation $y'-mxy=x^n$ with initial condition $(0,0)$, but I see now that you are!

You can also characterize the qualitative behavior for a fixed $m$ from the vector field plot, which easily justifies your intuition that $y$ is bounded when $m<0$ and grows rapidly when $m>0$.

Related Question