Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It only takes a minute to sign up.
Sign up to join this community
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?
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$.
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.
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.