I am trying to use Python and NumPy to use Euler’s formula e^i(π) represented as a matrix as e^A where

A = [0 -π] [π 0]

,and then apply it to the Maclaurin series for an exponential function e^x as

SUMMATION(n=0, infinity) x^n/n! = 1 x x^2/2! x^3/3! ...

So I am trying to compute an approximation matrix S^N 1 and print the matrix and it's four entries.

I have tried emulating euler's and maclaurin's series, which i think the final approximation matrix for this will be when N = 20, but currently my values do not add up. I am also trying to use np.linalg.norm to compute a 2 norm as well.

import math
import numpy as np

n = 0
A = np.eye(2)
A = math.pi * np.rot90(A)
A[0,1] = -A[0,1]
A
mac_series = 0
while n < 120:
  print(n)
  n  = 1
  mac_series = (A**n) / (math.factorial(n))

print("\n",mac_series)
np.linalg.norm(mac_series)

CodePudding user response：

The main problem here is that you are confusing A**3 with A@A@A.

Just look at case n=0.

A**0
#array([[1., 1.],
#       [1., 1.]])

I am pretty sure, you were expecting A⁰ to be identity (that is only that way that this thinking of x iy ⇔ np.array([[x,-y],[y,x]]) makes sense)

In numpy, you have np.linalg.matrix_power for that (or you could just accumulate power your self)

sum(np.linalg.matrix_power(A,i) / math.factorial(i) for i in range(20))

is

array([[-1.00000000e 00,  5.28918267e-10],
       [-5.28918267e-10, -1.00000000e 00]])

for example. Pretty sure that is what you were expecting (that is the matrix that represents real -1 using the same logic. And whole point of Euler identity is e^(iπ) = -1).

By comparison,

sum(A**i / math.factorial(i) for i in range(20))

returns

array([[ 1.        ,  0.04321392],
       [23.14069263,  1.        ]])

Which is just the maclaurin series computed for all four elements of the matrix. In other words, since your matrix is [[0,-π],[π,0]], you are evaluating using a MacLauring series [[e⁰, exp(-π)], [exp(π), e⁰]]. And it works. e⁰=1, obviously. exp(π) is 23.140692632779267, so we got a very good approximation in our result. And exp(-π) is the inverse, 0.04321391826377226. We also got a good approximation. So it works. Just not at all to do what you obviously intend to do: prove Euler identity's in matrix form; compute exp(iπ) not just exp(π).

Without matrix_power, and with a code closer to your initial code, you could

n=0
mac_series = 0
Apowern=np.eye(2) # A⁰=Id for now
while n < 20:
  print(n)
  mac_series  = Apowern / (math.factorial(n))
  Apowern = Apowern @ A # @ is the matrix multiplication operator
  n =1

Note that I've also moved n =1 which was misplaced in your code. You were stacking Aⁿ⁺¹/(n 1)! not Aⁿ/n! with your code (in other words, your sum misses the A⁰/0!=Id term).

With this, I get the expected result

>>> mac_series
array([[-1.00000000e 00,  5.28918724e-10],
       [-5.28918724e-10, -1.00000000e 00]])

Last problem, more subtle: you may have noticed that I do only 20 iterations, not 120. That is because after 20, you start to have a numerical problem. Apowern (or np.linalg.matrix_power(A,n), it is the same problem for both methods) becomes to big. Since it is divided by n! in the stacking, that doesn't prevent convergence. But it does prevent numeric convergence. And, in practice, after a while, numpy change the type of Apowern.

So, we should not have big matrix divided by big number, and try to iterate things that stay small enough. Like this for example

n=0
mac_series = 0
NthTerm=np.eye(2) # Aⁿ/n!. A⁰/0!=Id for now
while n < 120: # 120 is no longer a problem
  print(n)
  mac_series  = NthTerm 
  n  = 1
  NthTerm = (NthTerm @ A) / n # so if nthterm was 
  # Aⁿ/n!, now it becomes Aⁿ/n! @ A/(n 1) = Aⁿ⁺¹/(n 1)!

Result

>>> mac_series
array([[-1.00000000e 00, -2.34844612e-16],
       [ 2.34844612e-16, -1.00000000e 00]])

tl;dr

You have 4 problems

• The one already mentioned by Roy: you are not accumulating the Aⁿ/n!, just replacing them, and eventually keeping only the last. In other words, you need a = instead of =
• A**n is not Aⁿ. It is just A, with all the elements to the power n. Said otherwise [[x,-y],[y,x]]**n is not [[x,-y],[y,x]]ⁿ it is [[xⁿ,(-y)ⁿ],[yⁿ,xⁿ]]. So you'll end up computing [[e⁰, 1/e^π], [e^π, e⁰]] ≈ [[1, 0.0432], [23.14, 1]] which is irrelevant.
• n =1 is misplaced
• The numerical problem due to Aⁿ becoming huge (even if you intend to divide it by a even huger n!, so it does not theoretically/mathematically pose a problem, but numerically it does, since intermediate result is to big for computer)