2

(I'm not sure whether I should post this problem on this site or on the math site. Please feel free to migrate this post if necessary.)

My problem at hand is that given a value of k I'd like to numerically compute a rational function of nonlinear polynomials in k which looks like the following: (sorry I don't know how to typeset equations here...) where {a_0, ..., a_N; b_0, ..., b_N} are complex constants, {u_0, ..., u_N, v_0, ..., v_N} are real constants and i is the imaginary number. I learned from Numerical Recipes that there are whole bunch of ways to compute polynomials quickly, in the meanwhile keeping the rounding error small enough, if all coefficients were constant. But I do not think those ideas are useful in my case since the exponential prefactors also depend on k.

Currently I calculate it in a brute force way in C with complex.h (this is just a pseudo code):

double complex function(double k)
{
    return (a_0+a_1*cexp(I*u_1*k)*k+a_2*cexp(I*u_2*k)*k*k+...)/(b_0+b_1*cexp(I*v_1*k)*k+v_2*cexp(I*v_2*k)*k*k+...);
}

However when the number of calls of function increases (because this is just a part of my real calculation), it is very slow and inaccurate (only 6 valid digits). I appreciate any comments and/or suggestions.

7
  • 4
    Horner's method? a0 + k * (a_1*cexp(I*u_1*k) + k * (a_2*cexp(I*u_2*k) + k * ... )...))
    – zch
    Sep 19, 2013 at 23:08
  • What are the expected ranges of k and N?
    – rob mayoff
    Sep 19, 2013 at 23:08
  • @zch, yes I've been suggested Horner's rule. But I thought it's faster simply because one could store all constant coefficients in an array and multiply them using for loop, isn't it? Will test it later though. I just wanna make sure whether there's any similar situation (hopefully an algorithm :P) that has been analyzed thoroughly but I don't know yet. Thanks anyway!
    – Leo Fang
    Sep 20, 2013 at 1:12
  • @rob, k is in a wide range [-1000,1000] while N is from 1 to 20 according to the cases at hand.
    – Leo Fang
    Sep 20, 2013 at 1:14
  • 1
    @LeoFang it can be faster anyway. Without Horner you multiply *k about 200 times (for N ~= 20), with Horner it's about 20.
    – zch
    Sep 20, 2013 at 10:45

1 Answer 1

0

I trust that this isn't a homework assignment! Normally the trick is to use a loop add the next coefficient to the running sum, and multiply by k. However, in your case, I think the "e" term in the coefficient is going to overwhelm any savings by factoring out k. You can still do it, but the savings will probably be small.

Is u_i a constant? Depending on how many times you need to run this formula, maybe you could premultiply u_i * k (unless k changes each run). It's been so many decades since I took a Numerical Analysis course that I have only vague recollections of the tricks of the trade. Let's see... is e^(i*u_i*k) the same as (e^(i*u_i))^k? I don't remember the rules on imaginary numbers, or whether you'll save anything since you've got a real^real (assuming k is real) anyway (internally done using e^power).

If you're getting only 6 digits, that suggests that your math, and maybe your library, is working in single precision (32 bit) reals. Check your library and check your declarations that you are using at least double precision (64 bit) reals everywhere.

3
  • Although I haven't done any test yet, I would agree with your guess that factoring out k saves not too much (which is basically Horner's rule proposed by @zch in the comment). u_i's and v_i's are all constant, but k is going to be changed each run. e^(i*u_i*k) is indeed the same as (e^(i*u_i))^k, but k is a real number, not just an integer, so I don't think computing (e^(i*u_i))^k, which is done in C by using cpow(cexp(I*u_i),k), is a good idea.
    – Leo Fang
    Sep 20, 2013 at 1:23
  • BTW, I checked that complex.h does work with double precision (and even beyond double) providing the functions computed are not too complicated, so I think the only reasonable guess is that those nonlinearities in my function increase the rounding error dramatically.
    – Leo Fang
    Sep 20, 2013 at 1:24
  • Ah, so you're getting 15 or more digits, but only 6 of accuracy? I misinterpreted your earlier statement. Yeah, it sounds like even with double or quad math, the numerical methods need some tweaking. Beyond "watch out for subtracting similar-sized numbers, or adding/subtracting numbers vastly different in magnitude", I don't remember much of NA. Sorry.
    – Phil Perry
    Oct 30, 2013 at 19:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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