4
$\begingroup$

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?

  • 0
    [tag:reference-request] should not be used as a standalone tag; see [meta](http://meta.math.stackexchange.com/questions/2498/the-meta-tags). (I am not that sure that reference-request is a suitable tag here.)2012-08-09

3 Answers 3

3

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:

Plots of <span class=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.

6

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

  • 0
    I think $x$ in $-5,5$ would suffice.2012-04-08
6

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:

Plotted using

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

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

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.