Menu
▾
▴
visualpython-users
|
From: Andy D. <dou...@la...> - 2002-01-07 21:11:59
|
I've been considering using Vpython in a course to illustrate the
propagation of wave packets with dispersion. The problem is it's just too
slow. I'm hoping someone on this list with experience with the Numeric
module might be able to help me speed it up. I think I'm spending way
too much time explicitly looping over arrays.
Sorry for the length of what follows, but it might be helpful to be
specific:
Here's my original naive attempt:
#!/usr/bin/python
# Show propagation of a wave pulse with dispersion.
from visual.graph import *
def w(k):
"Dispersion relation w(k), gives v ~ 1/sqrt(k)."
return sqrt(k)
def yt(x,t):
"Total wave obtained by superimposing sinusoids with k from 1 to 4."
ytotal = 0.0
for k in arange(1.0, 4.0, 0.05):
ytotal = ytotal + sin(k * x - w(k) * t)
return ytotal
gdisplay(xtitle='x', ytitle='ytotal', title='Wave',
xmin = -10.0, ymin = -60.0, xmax =100.0, ymax = 60.0)
pulse = gcurve()
t = 0.0
# Make initial plot
for x in arange(-10.0, 100.0, 0.2):
pulse.plot(pos=(x, yt(x,t)))
npts = len(pulse.gcurve.pos)
while 1:
t = t + 0.5
for n in range(npts):
x = pulse.gcurve.pos[n][0]
pulse.gcurve.pos[n][1] = yt(x, t)
This would be rather simple for my students to follow (I hope!), but it
is painfully slow. I was able to speed it up considerably by doing
some obvious things:
1. Precomputing the dispersion relationship w(k) once and storing it
2. storing each set of points in a temporary array and updating the
display only once per time step (instead of every x value)
3. Calling math.sin() directly instead of the overloaded Numeric sin()
but it's still too slow. (Obviously, I can reduce the spatial and
temporal resolution, but that's beside the point here.)
To be specific, here's what I ended up with.
#!/usr/bin/python
# Show propagation of a wave pulse with dispersion.
from visual.graph import *
def w(k):
"Dispersion relation w(k), gives v ~ 1/sqrt(k)."
return sqrt(k)
# Range of wavenumbers 'k' to include.
ka = arange(1.0, 4.0, 0.05)
wk = w(ka)
nkmax = len(ka)
def yt(x,t):
"Total wave obtained by superimposing sinusoids with k from 1 to 4."
ytotal = 0.0
for nk in range(nkmax):
ytotal = ytotal + math.sin(ka[nk] * x - wk[nk] * t)
return ytotal
gdisplay(xtitle='x', ytitle='ytotal', title='Wave',
xmin = -10.0, ymin = -60.0, xmax =100.0, ymax = 60.0)
pulse = gcurve()
t = 0.0
for x in arange(-10.0, 100.0, 0.2):
pulse.plot(pos=(x, yt(x,t)))
# Make a temporary copy -- see below
temp = array(pulse.gcurve.pos)
npts = len(pulse.gcurve.pos)
while 1:
t = t + 0.5
# Compute the new values in a temp array that isn't displayed, then
# copy to the pulse.gcurve.pos array to update the display all at once.
for n in arange(npts):
temp[n][1] = yt(temp[n][0], t)
pulse.gcurve.pos = temp
Does anyone have any suggestions for how to significantly speed this up?
All the looping over arrays seems to me to be the sort of thing that I
ought to be able to get the Numeric extension to do for me efficiently,
but I confess I'm new to it and haven't seen how to do it yet.
Thanks for any advice,
Andrew Dougherty dou...@la...
Dept. of Physics
Lafayette College, Easton PA 18042
|
|
From: David S. <dsc...@vi...> - 2002-01-07 22:19:57
|
This happens to be an easy one.
Change the for loop in the na=EFve version from:
while 1:
t =3D t + 0.5
for n in range(npts):
x =3D pulse.gcurve.pos[n][0]
pulse.gcurve.pos[n][1] =3D yt(x, t)
to:
while 1:
t =3D t + 0.5
x =3D pulse.gcurve.pos[:,0]
pulse.gcurve.pos[:,1] =3D yt(x, t)
and you're done!
Dave
> -----Original Message-----
> From: vis...@li...=20
> [mailto:vis...@li...] On=20
> Behalf Of Andy Dougherty
> Sent: Monday, January 07, 2002 4:13 PM
> To: Vis...@li...
> Subject: [Visualpython-users] Speeding up a wave propagation=20
> simulation
>=20
>=20
> I've been considering using Vpython in a course to illustrate=20
> the propagation of wave packets with dispersion. The problem=20
> is it's just too slow. I'm hoping someone on this list with=20
> experience with the Numeric module might be able to help me=20
> speed it up. I think I'm spending way too much time=20
> explicitly looping over arrays.
>=20
> Sorry for the length of what follows, but it might be helpful to be
> specific:
>=20
> Here's my original naive attempt:
>=20
> #!/usr/bin/python
> # Show propagation of a wave pulse with dispersion.
> from visual.graph import *
>=20
> def w(k):
> "Dispersion relation w(k), gives v ~ 1/sqrt(k)."
> return sqrt(k)
>=20
> def yt(x,t):
> "Total wave obtained by superimposing sinusoids with k=20
> from 1 to 4."
> ytotal =3D 0.0
> for k in arange(1.0, 4.0, 0.05):
> ytotal =3D ytotal + sin(k * x - w(k) * t)
> return ytotal
>=20
> gdisplay(xtitle=3D'x', ytitle=3D'ytotal', title=3D'Wave',
> xmin =3D -10.0, ymin =3D -60.0, xmax =3D100.0, ymax =3D 60.0)
> pulse =3D gcurve()
>=20
> t =3D 0.0
> # Make initial plot
> for x in arange(-10.0, 100.0, 0.2):
> pulse.plot(pos=3D(x, yt(x,t)))
> npts =3D len(pulse.gcurve.pos)
> =20
> while 1:
> t =3D t + 0.5
> for n in range(npts):
> x =3D pulse.gcurve.pos[n][0]
> pulse.gcurve.pos[n][1] =3D yt(x, t)
>=20
> This would be rather simple for my students to follow (I=20
> hope!), but it is painfully slow. I was able to speed it up=20
> considerably by doing some obvious things:
> 1. Precomputing the dispersion relationship w(k) once=20
> and storing it
> 2. storing each set of points in a temporary array and=20
> updating the=20
> display only once per time step (instead of every x value)
> 3. Calling math.sin() directly instead of the overloaded=20
> Numeric sin() but it's still too slow. (Obviously, I can=20
> reduce the spatial and temporal resolution, but that's beside=20
> the point here.)
>=20
> To be specific, here's what I ended up with.
>=20
> #!/usr/bin/python
> # Show propagation of a wave pulse with dispersion.
> from visual.graph import *
>=20
> def w(k):
> "Dispersion relation w(k), gives v ~ 1/sqrt(k)."
> return sqrt(k)
>=20
> # Range of wavenumbers 'k' to include.
> ka =3D arange(1.0, 4.0, 0.05)
> wk =3D w(ka)
> nkmax =3D len(ka)
> def yt(x,t):
> "Total wave obtained by superimposing sinusoids with k=20
> from 1 to 4."
> ytotal =3D 0.0
> for nk in range(nkmax):
> ytotal =3D ytotal + math.sin(ka[nk] * x - wk[nk] * t)
> return ytotal
>=20
> gdisplay(xtitle=3D'x', ytitle=3D'ytotal', title=3D'Wave',
> xmin =3D -10.0, ymin =3D -60.0, xmax =3D100.0, ymax =3D 60.0)
> pulse =3D gcurve()
>=20
> t =3D 0.0
> for x in arange(-10.0, 100.0, 0.2):
> pulse.plot(pos=3D(x, yt(x,t)))
>=20
> # Make a temporary copy -- see below
> temp =3D array(pulse.gcurve.pos)
> npts =3D len(pulse.gcurve.pos)
>=20
> while 1:
> t =3D t + 0.5
> # Compute the new values in a temp array that isn't=20
> displayed, then
> # copy to the pulse.gcurve.pos array to update the=20
> display all at once.
> for n in arange(npts):
> temp[n][1] =3D yt(temp[n][0], t)
> pulse.gcurve.pos =3D temp
>=20
> Does anyone have any suggestions for how to significantly=20
> speed this up? All the looping over arrays seems to me to be=20
> the sort of thing that I ought to be able to get the Numeric=20
> extension to do for me efficiently, but I confess I'm new to=20
> it and haven't seen how to do it yet.
>=20
> Thanks for any advice,
>=20
> Andrew Dougherty dou...@la...
> Dept. of Physics
> Lafayette College, Easton PA 18042
>=20
>=20
> _______________________________________________
> Visualpython-users mailing list=20
> Vis...@li...
> https://lists.sourceforge.net/lists/listinfo/visualpython-users
>=20
|
|
From: Andy D. <dou...@la...> - 2002-01-08 14:22:37
|
On Mon, 7 Jan 2002, David Scherer wrote:
> This happens to be an easy one.
>=20
> Change the for loop in the na=EFve version from:
>=20
> for n in range(npts):
> x =3D pulse.gcurve.pos[n][0]
> pulse.gcurve.pos[n][1] =3D yt(x, t)
>=20
> to:
>=20
> x =3D pulse.gcurve.pos[:,0]
> pulse.gcurve.pos[:,1] =3D yt(x, t)
Ah, wonderful. Thanks, that indeed does it and was the sort of insight I
was looking for.
(Even looking back, I still don't see in the NumPy documentation just why
this works. It looks to me from the documentation that only "ufuncs"
operate on arrays element-by-element, not arbitrary user-defined
functions. Yet the actual behavior is exactly what I want in this case.)
Thanks again,
Andy Dougherty dou...@la...
Dept. of Physics
Lafayette College, Easton PA 18042
|
|
From: David S. <dsc...@vi...> - 2002-01-08 15:08:56
|
> (Even looking back, I still don't see in the NumPy
> documentation just why this works. It looks to me from the
> documentation that only "ufuncs" operate on arrays
> element-by-element, not arbitrary user-defined functions.
> Yet the actual behavior is exactly what I want in this case.)
Python will of course allow you to call a function with parameters of
any type, including arrays. So the call
pulse.gcurve.pos[:,1] = yt(x, t)
to the function
def yt(x,t):
ytotal = 0.0
for k in arange(1.0, 4.0, 0.05):
ytotal = ytotal + sin(k * x - w(k) * t)
return ytotal
is just like writing
ytotal = 0.0
for k in arange(1.0, 4.0, 0.05):
ytotal = ytotal + sin(k * x - w(k) * t)
pulse.gcurve.pos[:,1] = ytotal
The only operations that operate on the array x in this code are
arithmetic operators and sin, which is a ufunc. So it works fine. If
yt() attempted to do something with x that isn't supported for arrays,
you would get a TypeError at that point.
Incidentally, if you wanted to eliminate the "for k" loop as well, you
can do it with more Numeric code:
def yt(x,t):
k = arange(1.0, 4.0, 0.05)[:,NewAxis]
return sum( sin(k * x - w(k) * t) )
and, of course, you could calculate k and w(k) just once. But since
readability was one of your goals, I didn't recommend it. (If you want
to use this code, you will have to change the first call yt(x,t) in the
plotting loop to yt(x,t)[0], since ytotal comes back a rank-1 array
instead of a scalar)
Dave
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X