August 8th, 2014, 1:30 pm
Let's say you are trying to solve du/dt = A u with A the described sparse array and u(0) given? (If your problem is not in this form, I would try to reformulate it so it is.)Then, in Mathematica, which I use, this would be quite tractable with quite large A, say size N x N = 10^4 x 10^4 you mention.You just give it to the built-in ODE system solver.Alternatively, as you mention, a spectral approach might work. However, there's no need to compute an inverse.Again in Mathematica or your favorite CAS, just ask for the first n (orthonormal) eigenvectors [$]v_n[$]/eigenvalues [$]\lambda_n[$] of A,where n is small (< 10). Then, see what kind of convergence and run-time you get with(*) [$] u(T) = \sum_{j=1}^n e^{\lambda_j T} (v_j, u(0)) \, v_j[$], where [$](a,b) = a \cdot b[$], the usual vector inner product. [The solution (*) with [$]n = N[$] is exact, and follows from the spectral theorem for (in this case), real, symmetric matrices]. You will want the first (or last) n eigenvalues, with the appropriate ordering so these leading terms are dominant.Sometimes, this kind of thing convergences extremely rapidly; it depends on T and the eigenvalue spacing.Indeed, the relative error in (*) is [$]O(e^{(\lambda_{n+1} - \lambda_1) T})[$], so the larger T, the better. So, to my mind, a conventional CAS approach should handle this kind of thing, making parallelization(or indeed writing any significant amount of code in any system at all) unnecessary ..
Last edited by
Alan on August 7th, 2014, 10:00 pm, edited 1 time in total.