TransWikia.com

Calculating residue of a rational function

Computational Science Asked on October 23, 2021

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)= limlimits_{zto 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?

One Answer

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.

Answered by cdalitz on October 23, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP