Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|cite|improve this question
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
up vote 2 down vote accepted

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

Here's the output:

Plots of $P_{n+1}(x)$ using SAGE

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

share|cite|improve this 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
assume(x > 0)
assume(m > 0)
y = e^(m*x^2/2)
G += plot(p[-1].subs(m=1),-xw,xw)
G += plot(p[-1].subs(m=1),-xw,xw,color='green')
for k in range(3):
    G += plot(p[-1].subs(m=1),-xw,xw,color=(1,k/3,k/3)),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()
    assume(x > 0)
    y = e^(m*x^2/2)
    G += plot(p[-1],-x_max,x_max)
    G += plot(p[-1],-x_max,x_max,color='green')
    for k in range(n):
        G += plot(p[-1],-x_max,x_max,color=(1,k/n,k/n)),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()
    assume(x > 0)
    y = e^(m*x^2/2)
    G += plot(p[-1],-x_max,x_max)
    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
        G += plot(z,-x_max,x_max,color=(1,k/n,k/n))

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:


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

share|cite|improve this answer
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;
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);

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:

Plotted using "plot"

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

enter image description here

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

Plotted 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:

Plotted negative part using "semilogy"

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.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.