# Need (free) software to plot recursively defined functions.

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?

-
sage, perhaps even online. –  bgins Apr 7 '12 at 23:04
You can use octave... –  Quimey Apr 7 '12 at 23:21
why do you need these? what's your real problem? what are you trying to achieve? –  bgins Apr 7 '12 at 23:50
wait, does $\{\cdot\}$ mean the fractional part? –  Rahul Apr 8 '12 at 5:07
It's customary to put the subscript by the letter and before the parentheses: $P_{n+1}(x)$, and to cite your previous work or otherwise give the mathematical context. It helps to know that you are trying to characterize the behavior of, and that the $P_n$ are the particular solutions to, the differential equation $y'-mxy=x^n$ with initial condition $(0,0)$. –  bgins Apr 8 '12 at 10:01

I was thinking about this post again and decided a SAGE script may help, and it did! The graphs look much better. It's just Python using SAGE defined functions and methods:

x = var('x');
t = var('t');
G = Graphics();
m = 5;
forget() # Only run this command in an isolated instance!
assume(x > 0)  # Doesn't really matter if x >= 0 or x < 0, I checked.
y = e^(.5*m*(x^2));
P = [0 for l in range(11)];
P[0] = y*integrate(e^(-.5*m*(t^2)), (t, 0, x));
P[1] = (y - 1)/m;

for l in xrange(1, 10):
P[l + 1] = (P[l - 1] - ((x^l)/l))*l/m;

k = 5;
for n in range(k):
G += plot(P[n], xmin = -pi/2, xmax= pi/2, ymin = -10, ymax = 10, color = (1, 5/4 + .12*n, 5/4 + .12*n), legend_label = 'Plot of P[%s]'%(str(n)))
G.show()


Here's the output:

You can always turn the script into a function of $m$ and $k$ as well if you'd like.

-

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)


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)


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$.

-
How do you know which function is which? –  Pedro Tamaroff Apr 7 '12 at 23:31
p starts as an empty array. then i start appending. when the array has some elements, indexing starts at 0 for the first element. the index -1 is a shorthand for the last element in the array, just as if it were computing a remainder modulo the array size. this is based on python but perl behaves the same way. (you can't do this in c/c++ and i think not in java but you can in java's scripting dialect, groovy.) –  bgins Apr 7 '12 at 23:45
I have no idea what you're saying. I know nothing about programming. Can´t this be ploted without so many coding? –  Pedro Tamaroff Apr 7 '12 at 23:46
How many values of $m$ do you need and over what range of $x$? What do you really need here? –  bgins Apr 7 '12 at 23:49
I´d simply like to see how the functions behave for various $m$, and I'm thinking I'd prefer $m<0$ –  Pedro Tamaroff Apr 7 '12 at 23:51

Try Octave first. It is an excellent free open-source alternative to MATLab. As far as programming goes, you're going to need to do some programming. Though, a discussion of this would be better suited on StackExchange Programmers, scripting of the sort encountered in Octave (or MATLab) may not exactly be termed "programming;" but, that's beside the point. Here's a small Octave/MATLab script for accomplishing what you want, though if you're not familiar with MATLab or Octave, search for some tutorials for MATLab (since the two are almost identical, 97% $\pm$ 2%). I'm going to assume $m$ is constant, if you need a variable instead, I recommend setting up two arrays for $x \text{ and } m$ (perhaps equal, but certainly of equal size), and using meshgrid and surf, or simply plot3 if you're looking for a parametric representation:

x = linspace(-10, 10, 500);
k = 1:500; # This is the indexing array for P. Notice, x and k are equal
# in length.
N = 20;
m = 5;
figure
hold on

P(1, k) = exp(.5*m.*(x.^2)).*trapz(x.*exp(-.5*m.*(x.^2)));
P(2, k) = exp(.5*m.*(x.^2))./m;

for n = 2:N
P(n + 1, k) = (n/m).*(P(n - 1, k) - (x.^n)./n);
end

for n = 1:N
plot(x, P(n, k)) # You may opt for "semilogy" here instead of
end                  # "plot" though that would only show you the
# positive y-axis.

axis square
axis([-5 5 -100 100])


Here's the visual output for the script:

Here's the same output, with much greater y-axis limits:

Here's the log-plot for the positive part using "semilogy":

And, here's the log-plot for the negative part using "semilogy" It's been rotated 180° in the plane, however:

Don't take my word here, though, verify it!

Finally, I have to warn you, there is a good reason MATLab costs thousands of dollars. The commercially developed software can run circles around free software like Octave whenever resource-intensive computations are involved. I've limited the arrays x and n here to be rather small, which in turn may greatly affect the quality of the plots produced. If you really need that horsepower, however, you can either invest \$100, if you're a student, for a MATLab student license or the time to learn SAGE, as it's the next best thing for this sort of work. As far as numerical analysis goes, I'm not sure about Mathematica; though, it's stellar for most other applications. Finally, I'm not endorsing any product; I have simply provided advice based on experience and what I've picked up after spending a good chunk of time in Engineering and Mathematics departments.

-