0
$\begingroup$

I think this is a numeric problem but would like to see where it occurs. I create a complex rational polynomial as follows and make it a function of f

 bf = ButterworthFilterModel[{"Lowpass", 2, 2 \[Pi] 500.}];
  a = bf[I 2 \[Pi] f][[1, 1]]

(*
 9.8696*10^6/(((2221.44 - 2221.44 I) + 
   2 I f \[Pi]) ((2221.44 + 2221.44 I) + 2 I f \[Pi]))
*)

Now I use two methods for calculating the modulus squared.

b1 = ComplexExpand[Re[a Conjugate[a]]] // Simplify
b2 = ComplexExpand[Abs[a]^2] // Simplify

Why do I get two diferent answers?

(*

(0. + 0. I) + 3.90625*10^21/(6.25*10^10 + 1. f^4)^2 - (
 3.58535 f^2)/(6.25*10^10 + 1. f^4)^2 + (
 6.25*10^10 f^4)/(6.25*10^10 + 1. f^4)^2

(0. + 0. I) + 6.25*10^10/(6.25*10^10 + 1. f^4)
*)

I think this is to do with numerics because if I repeat with exact numbers

bf = ButterworthFilterModel[{"Lowpass", 2, 2 \[Pi] 500}];
a = bf[I 2 \[Pi] f][[1, 1]];
b1 = ComplexExpand[Re[a Conjugate[a]]] // Simplify
b2 = ComplexExpand[Abs[a]^2] // Simplify

I get two results the same.

(* 
62500000000/(62500000000 + f^4)

62500000000/(62500000000 + f^4)
*)

The problem is that I can't necessarily work with exact numbers. So why do I get different answers? Is the problem in ComplexExpand or Simplify? Is there a preferential approach (I note that Abs[]^2 gives simpler answers)? Thanks

$\endgroup$
  • $\begingroup$ Rationalize[ approx, 0] will make an exact constant from the approximate constant approx. $\endgroup$ – John Doty Oct 29 '19 at 13:02
  • $\begingroup$ @JohnDoty Thanks. I will look into this but I am not sure I can do this with more complex rational polynomials for my cases $\endgroup$ – Hugh Oct 29 '19 at 13:05
  • 1
    $\begingroup$ As you go to more complicated polynomials, numerical instability problems worsen. Manipulating roots of polynomials can be problematic, too. But for those, Mathematica offers exact representation as Root objects, which you may use without stability concerns. For numerics, you'll be in the more difficult territory of "which polynomial form is best for this problem?" It's never power series... $\endgroup$ – John Doty Oct 29 '19 at 13:41
  • $\begingroup$ While the structure of the two results are different, they are numerically equivalent. Maximize[{Abs[b1 - b2], 0 < f < 10^9}, f] evaluates to {2.22045*10^-16, {f -> 0.0525412}}. $\endgroup$ – Bob Hanlon Oct 29 '19 at 14:01
  • $\begingroup$ @BobHanlon Interesting result. If you apply Simplify and FullSimplify to the difference of the answers you also get different answers. I would never have thought two different expressions could be so similar. I guess this is why Pade approximations are so good. $\endgroup$ – Hugh Oct 29 '19 at 14:09
3
$\begingroup$

You may be interested what is the "right" answer. This can be easily found

bf = ButterworthFilterModel[{"Lowpass", 2, k}];
a = bf[I 2 π f][[1, 1]];
b1 = ComplexExpand[Re[a Conjugate[a]]] // Simplify
b2 = ComplexExpand[Abs[a]^2] // Simplify

Leading to $$\frac{k^4}{16 \pi ^4 f^4+k^4}$$ in both cases. I suggest to do numerics on this stage.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you. This is useful. I should have stressed that this is a minimum working example. For more complicated filters, e.g. band pass of high degree the ability to work out the roots becomes more difficult and a numerical method is necessary. Still good to see the "right" answer. $\endgroup$ – Hugh Oct 29 '19 at 19:57
1
$\begingroup$
bf = ButterworthFilterModel[{"Lowpass", 2, 2 π 500.}];

Convert a to exact values

a = Rationalize[(bf[I 2 π f][[1, 1]]), 0] // Simplify;

b1 = ComplexExpand[Re[a Conjugate[a]]] // Simplify;

b2 = ComplexExpand[Abs[a]^2] // Simplify;

b1 and b2 are identical

b1 === b2 // FullSimplify

(* True *)

Further, with exact values it is not necessary to use Re in the definition of b1

b1 = ComplexExpand[a Conjugate[a]] // Simplify;

b1 === b2 // FullSimplify

(* True *)
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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