0

I have a function $$ f(z) = \frac{1}{(z-z_1)(z-z_2)(z-z_3)(z-z_4)} $$ All of $\{z_1,z_2,z_3,z_4\}$ are simple poles. The residues for this function are given as $$ \text{Res}(f(z),z_i)= \lim\limits_{z\to zi} \frac{(z-z_i)}{(z-z_1)(z-z_2)(z-z_3)(z-z_4)} $$ For example to find $\text{Res}(f(z),z_1)$ one first cancels the $(z-z_1)$ numerator and denominators and then takes the limit $z \to z_1$. so the result is $$ \text{Res}(f(z),z_1) = \frac{1}{(z_1-z_2)(z_1-z_3)(z_1-z_4)} $$ Similarly one can find $\text{Res}(f(z),z_2)$, $\text{Res}(f(z),z_3)$, $\text{Res}(f(z),z_4)$.

I want to implement this residue finding algorithm in a function. In Cpp I tried to implement this like this

double z1 = 1.0;
double z2 = 2.0;
double z3 = 3.0;
double z4 = 4.0;

auto res = [&](double z){
return [&](double zi){
    return (z-zi)/((z-z1)(z-z2)(z-z3)*(z-z4));
}(z);
};

This returns -nan when I compute res(z1) as the function becomes of $\frac{0}{0}$ form. I wanted to define a function that will first get rid of the common factor in the numerator and the denominator and then puts the value $z_1$ in the function. For simple enough functions with simple poles, this should be enough to find the residue.

How to do this in Cpp?

Galilean
  • 151
  • 6
  • How do you feel about switch case statements? – Abdullah Ali Sivas Jul 20 '20 at 18:55
  • @AbdullahAliSivas yes! that is an option I've been doing till now, Hard coding all the expressions of the residues. – Galilean Jul 21 '20 at 01:07
  • That is pretty much the way to go. Doing symbolic computation is expensive and complicated, so if you have that solution already stick with it – Abdullah Ali Sivas Jul 21 '20 at 23:42
  • What is your input? The vector of the $z_i$s? And what is your output? If that function res is the proper signature, how is it supposed to behave for $zi \not \in {z_1,z_2,z_3,z_4}$? – Federico Poloni Jul 23 '20 at 18:52
  • @FedericoPoloni for any other point which is not a pole, the residue function Will just output the number computed by the function. – Galilean Jul 23 '20 at 20:08

1 Answers1

1

Instead of hard coding all cases with a switch clause, you can parametrize the function by its poles:

double residue(size_t i, const std::vector<double> &poles) {
  double res = 1.0;
  for (size_t j=0; j < poles.size(); j++) {
     if (j != i) {
        res *= 1 / (poles[i] - poles[j]);
     }
  }
  return res;
}

As a side note, I wonder whether de l'Hospital's rule might be helpful in case of less simple functions.

cdalitz
  • 481
  • 3
  • 7