Numerical instability can be a common issue when calculating the matrix exponential, especially for large matrices or when dealing with complex numbers.
In your code, the high condition number suggests that the matrix is ill-conditioned, which can lead to such instability. The condition number can be evaluated using MATLAB's “cond” function. You can read more about “cond” function here: https://www.mathworks.com/help/matlab/ref/cond.html
To mitigate numerical instability, you can consider the following strategies:
- Matrix Balancing: You can use matrix balancing to improve the conditioning of the matrix before computing the exponential. MATLAB's “balance” function can help with this. You can read more about it here: https://www.mathworks.com/help/matlab/ref/balance.html
- Regularization: Consider adding a small regularization term to the matrix to improve its condition number. This can be done by adding a small multiple of the identity matrix, “epsilon * eye(d)”, to H.
- Alternative Matrix Exponential Methods: You can also use “Krylov subspace” techniques or other specialized libraries designed for stable matrix exponential calculations. You can read about it here: https://www.mathworks.com/matlabcentral/fileexchange/45170-jacobian-free-newton-krylov-jfnk-method
Here’s my modified version of your code involving regularization and matrix balancing to reduce condition number and bring numerical stability.
d_values = 10:110; % Range of d values to scan
t = 0.6; % Time value
first_element = zeros(length(d_values),1); % Preallocate array to store the first element of vec1
for idx = 1:length(d_values)
d = d_values(idx);
A = diag(sqrt(1:d-1), 1);
Ad = A';
A3 = A^3;
Ad3 = A3';
H = (Ad3 + A3);
% Regularization (small identity matrix addition)
epsilon = 1e-10; % Regularization parameter, adjust as needed
H_reg = H + epsilon * eye(d);
% Balance the matrix
[H_bal, D] = balance(H_reg);
% Check condition number
cond_number = cond(H_bal);
fprintf('Condition number for d = %d: %e\n', d, cond_number);
vec0 = zeros(d,1);
vec0(1) = 1;
% Compute the matrix exponential
try
U_tsq = expm(-1i * t * H_bal);
catch ME
warning('Matrix exponential computation failed for d = %d: %s', d, ME.message);
continue;
end
vec1 = U_tsq * vec0;
first_element(idx) = vec1(1);
end
% Plot d vs the first element of vec1
plot(d_values, abs(first_element), '-o');
xlabel('matrix size');
ylim([0,1]);
ylabel('Abs of the first element of vec1');
grid on;
Hope it helps!
