Enhancement: The same code in Python

[StSav012]

The following code is a re-write of your wonderful piece of code in Python. The differences are the following:

• It assumes the initial data are biased by a constant.
• It scales to any number of exponential functions. Change N for that.

The required Python version should start from 3.9, although I used only 3.11 to run the code.

The additional packages required to be installed are numpy and scipy. As far as I have found out, the code should work with the packages some 2 years old (i.e., numpy>=1.21.0 and scipy>=1.6.0), so basically with any version available now.

from numpy import arange, column_stack, exp, eye, ones
from numpy.linalg import eigvals, pinv
from numpy.typing import NDArray
from scipy.integrate import cumulative_trapezoid

# get data
dx: float = 0.02
x: NDArray[float] = arange(dx, 1.5 + dx, dx)
y: NDArray[float] = -1 + 5 * exp(0.5 * x) + 4 * exp(-3 * x) + 2 * exp(-2 * x) - 3 * exp(0.15 * x)

# number of exponents in question
N: int = 4

# calculate integrals
iy: list[NDArray[float]] = [cumulative_trapezoid(y, x, initial=0)]
for n in range(1, N):
    iy.append(cumulative_trapezoid(iy[-1], x, initial=0))

# get exponents' lambdas
Y: NDArray[float] = column_stack((*iy, *[x**n for n in range(N, 0, -1)], ones(x.shape)))
A: NDArray[float] = pinv(Y) @ y

lambdas: NDArray[complex] = eigvals([A[:N], *eye(N - 1, N)])

# get exponents' multipliers
X: NDArray[complex] = column_stack((ones(x.shape), *[exp(lambdas[n] * x) for n in range(N)]))
P: NDArray[complex] = pinv(X) @ y

# print out the results, sorting by the magnitude of the multiplier
order_by_multiplier: NDArray[int] = abs(P[1:]).argsort()[::-1]
print(f'y(x) =\n    {P[0]:.4f}')
for i in order_by_multiplier:
    print(f'    {P[i + 1]:+.4f} × exp({lambdas[i]:.4f} x)')

The result I get is this:

y(x) =
    -0.9900
    +5.0018 × exp(0.5000 x)
    +4.0052 × exp(-2.9991 x)
    -3.0109 × exp(0.1500 x)
    +1.9939 × exp(-1.9997 x)

[juangburgos]

Awesome, if you want to add a PR to have it alongside the matlab code, I'd be happy to merge it. It'd have the same license as the original one, though you might want to put you name in a comment at the top of the file if you want.