I’m trying to compute the step response of a linear system with four degrees of freedom in Julia.
This can be done reasonably efficiently in Python using NumPy broadcasting:
import numpy as np
def step_analytical(A, B, F, n):
# Eigendecomposition: A = V D V_inv
eigvals, V = np.linalg.eig(A)
V_inv = np.linalg.inv(V)
# Precompute constant term
BF_transformed = V_inv @ (B * F) # (4, 1)
# Time steps (n,)
t = np.arange(1, n + 1)
# Compute exp(D * t) for all t: shape (n, 4)
exp_diag = np.exp(np.outer(t, eigvals)) # (n, 4)
# Compute (exp - 1) / lambda
diag_terms = (exp_diag - 1) / eigvals # (n, 4)
# Scale V columns by diag_terms
scaled_V = V[None, :, :] * diag_terms[:, None, :] # (n, 4, 4)
# Multiply by BF in transformed space
x = (scaled_V @ BF_transformed).squeeze(-1) # (n, 4)
return x
# Dynamics matrix A
A = np.array([[-0.19475053, 0.09852727, 0. , 0. ],
[ 0.01789924, -0.04116155, 0.02326231, 0. ],
[ 0. , 0.01585978, -0.01863287, 0.00277309],
[ 0. , 0. , 0.00061425, -0.00061425]])
# Input matrix B
B = np.array([[0.12675659],
[0. ],
[0. ],
[0. ]])
# Constant forcing term F
F = np.float64(5.156003293238043)
# Number of time steps n
n = 5000
Here is my attempt at a Julia implementation:
using LinearAlgebra
function step_reference(A, B, F, n)
# Eigendecomposition: A = V D V_inv
λ, V = eigen(A)
V_inv = inv(V)
# Precompute constant term
BF_transformed = V_inv * (B * F) # (4, 1)
# Time steps (n,)
t = 1:n
# Compute exp(D * t) for all t
exp_diag = exp.(λ * t') # (4, n)
# Compute (exp - 1) / λ
diag_terms = (exp_diag .- 1) ./ λ # (4, n)
# Scale V columns by diag_terms: shape (4, 4, n)
V_scaled = reshape(V', 4, 4, 1) .* reshape(diag_terms, 4, 1, n)
# Multiply by BF in transformed space
x = zeros(eltype(λ), 4, n)
for j in 1:n
x[:, j] = V_scaled[:, :, j]' * BF_transformed
end
return x
end
# Dynamics matrix A
A = [
-0.19475053 0.09852727 0.0 0.0;
0.01789924 -0.04116155 0.02326231 0.0;
0.0 0.01585978 -0.01863287 0.00277309;
0.0 0.0 0.00061425 -0.00061425
]
# Input matrix B
B = [
0.12675659;
0.0;
0.0;
0.0
]
# Constant forcing term F
F = 5.156003293238043
# Number of time steps n
n = 5000
I’ve reimplemented most of the function using broadcasting but I’m stuck at the final step, where I resort to a loop.
Unlike in NumPy, Julia doesn’t seem to support batched matrix multiplication without additional packages, which leads me to believe that what I’m trying to do is not idiomatic Julia. I notice that my Julia code here is also running slower than NumPy, possibly due to the large number of allocations.
Any hints on how to improve this would be much appreciated.