Virtual Polywell

drmike · Post by **drmike** » Mon Mar 10, 2008 12:52 am

Yeah, that's kinda been my experience. For really learning things, it's a good thing to do.

I posted the full version of the latest Virtual Polywell Model. Over the next week or so I'll start writing the code that implements it. Other than the electron source, it should pretty much just be cut and paste from previous stuff, so hopefully it'll be debugging more than writing.

It'll be interesting!

scareduck · Post by **scareduck** » Mon Mar 10, 2008 4:57 pm

It just gets weirder and weirder. I found a chapter from a proposed book on numerical methods by a couple professors at the Technische Universitaet Carolo-Wilhelmina (in English, no less) that gives analytic roots of the Legendre polynomials through n=4, and the results I'm getting are the same. So far, so good, but I'm utterly puzzled by the lack of agreement by the numbers that appear in the GSL -- still. One interesting clue: every other root seems to agree between the code I'm running and the GSL Kronrod constants.

I just ordered the Rabinowitz book, even though it's a little dated, because it seems to have been well-regarded by various other sources. drmike, if you can recommend any books on the subject, that would be most welcome, too.

Edit: Apparently -- and I didn't notice this when I ordered it -- the second edition was published just last year. Good stuff.

drmike · Post by **drmike** » Tue Mar 11, 2008 2:00 pm

Sorry, the Kronrod stuff is not in my list of things I've got texts for. Elliptic curve math for crypto is a different story

scareduck · Post by **scareduck** » Wed Mar 12, 2008 5:12 am

Slow, slower, slowest ... while waiting for the Rabinowitz book to arrive, I managed to dig up the original Kronrod treatise on the subject. The first nine pages of his thin little book (most of it's numerical tables) is dedicated to proofs. I'll be reading it soon ...

And I must say ... I received the 3rd edition of Numerical Recipes, and was sorely disappointed to see the same treatment on Kronrod as was in the second edition. Oh, well.

drmike · Post by **drmike** » Wed Mar 12, 2008 1:36 pm

But you never would have heard of it otherwise. So in the sense that NR gives keywords on what to go look for, it's a good place to start.

My kids wanted me to read one of their fun books, so I've been putting off coding. It's a great break to read about monsters and gods sending email, but puzzles of reality are more fun for me. The book is "The Sea of Monsters" and the main characters are greek heros (human children of gods) born in the modern world. Hopefully I'll be back to coding soon.

scareduck · Post by **scareduck** » Wed Mar 12, 2008 3:21 pm

Okay, now I have a better understanding of what's going on. Gauss-Legendre quadrature uses exclusively Legendre polynomials L(n,x), which looks in pseudocode like this:

Code: Select all

for(k=0;k<=n;k++)
  s+=w[k]*L(n,x);

where w[k] is a weight vector based on n. Kronrod interposed another unknown but orthogonal polynomial P(n+1,x) into the mix, so the equation became

Code: Select all

for(k=0;k<=n;k++)
  s+=W[k]*L(n,x)*P(n+1,x)*x^k;

Kronrod didn't know what P(n+1,x) would look like, and didn't prove that there would always be one available, but he did manage to generate them up to n=40. The roots of P(n+1,x) would always be interspersed between the roots of the Gaussian Legendre polynomial. Patterson showed that for any n, P(n+1,x) would always exist, and furthermore that it could be written in terms of the Legendre polynomials itself. He also found a numerically stable way to calculate the weights using a recurrence function, something Kronrod had failed to do; Kronrod ended up carrying 65 decimal digits to calculate only 12-14 digits of precision.

Far more detail than I'm sure anyone cared to read, but there it is.

drmike · Post by **drmike** » Thu Mar 13, 2008 2:17 am

You'd think we would be using K(n+1, x) instead of P()

Math details are always interesting. Eventually they come back to haunt you, so having seen it some place once is going to be useful. Thanks!

scareduck · Post by **scareduck** » Thu Mar 13, 2008 4:30 am

Actually, I'm using Kronrod's own notation there. Others (Patterson in particular) do use K(n+1,x).

I got the Davis and Rabinowitz book today, and they spend several pages treating Gauss-Kronrod in some detail. Needless to say, this was cause for no small amount of celebration in these parts.

drmike · Post by **drmike** » Thu Mar 13, 2008 4:33 pm

scareduck wrote:Actually, I'm using Kronrod's own notation there. Others (Patterson in particular) do use K(n+1,x).

I got the Davis and Rabinowitz book today, and they spend several pages treating Gauss-Kronrod in some detail. Needless to say, this was cause for no small amount of celebration in these parts.

Sounds great! It'll take me a while to convert my math into code, so not much celebration here yet. Don't have too much fun.

scareduck · Post by **scareduck** » Sat Mar 15, 2008 5:26 am

The truth is, I'm a lightweight at this stuff. Real math majors kick sand on me at the beach. Fortunately, I have persistence and Google to help:

http://www.netlib.org/toms/672

Well-formed C it's not, but it's translatable.

drmike · Post by **drmike** » Sat Mar 15, 2008 2:12 pm

Just use f77, that compiles fortran directly. No need to translate. I think there have been upgrades to fortran in the '90's too, but that code all says "FORTRAN-77" so I don't think you need to worry about it.

Hard to believe fortran is still used, but no point in debugging old code that works!

scareduck · Post by **scareduck** » Sat Mar 15, 2008 6:47 pm

drmike wrote:Just use f77, that compiles fortran directly. No need to translate. I think there have been upgrades to fortran in the '90's too, but that code all says "FORTRAN-77" so I don't think you need to worry about it.

Hard to believe fortran is still used, but no point in debugging old code that works!

Well --

1) No matter what happens, I'll need to use extra bits to get better precision on the abscissae and weights for final use with the qd libs.
2) Ergo, I'll need to use the GMP lib.
3) I don't hate FORTRAN, but I'd just as soon avoid it. And if I have to use the GMP lib, that means translating two things backwards instead of just one. Blecch.

drmike · Post by **drmike** » Sat Mar 15, 2008 8:27 pm

Yeah, that makes sense. I think there's an automatic translator because f77 itself is just an ftoc type translator. A quick look and I found f2c which says:

Source for f2c, a Fortran 77 to C translator jointly developed by
folks from Bell Labs, Bellcore, and Carnegie Mellon, is now freely
available.

F2c was derived from the original UNIX operating system's f77(1),
and the generated C follows f77's calling conventions; on some machines, the
resulting object files are interchangeable with (and behave
indistinguishably from) objects compiled by f77.

The date on that is 1990!

scareduck · Post by **scareduck** » Sun Mar 16, 2008 6:31 am

Turns out some folks have been writing a fully C++-ized scientific library called O2scl. Haven't looked at it closely yet, but it could come in handy.

drmike · Post by **drmike** » Sun Mar 16, 2008 2:19 pm

Yup, that looks interesting. I note it says "pre-beta", but it probably is a good place to start. All the routines you want are probably done, and it is just a matter of double checking them.

I'm in the process of writing code for multiple particle tracking. It's interesting because I can account for current lost to MaGrid, external wall or electron source all independently. I haven't quite figured out how to follow the data yet, but I think showing the particle electric field as an independent entity may be the simplest thing I can do. Somehow I have to show the particles too, but I'm not sure yet how to do that.