Reid.Atcheson
  • Member for 10 years
  • Last seen this week
How mature is the "Julia" scientific computing language project?
23 votes

I believe Julia is worth learning. I have used it to produce a few research finite element codes, and produce them very quickly. I have been over all very pleased with my experience. Julia has ...

View answer
Why do equi-spaced points behave badly?
22 votes

This is a really interesting question, and there are a lot of possible explanations. If we are attempting to use a polynomial interpolation, then note that polynomial satisfy the following annoying ...

View answer
What is the principle behind the convergence of Krylov subspace methods for solving linear systems of equations?
Accepted answer
22 votes

In general, all Krylov methods essentially seek a polynomial that is small when evaluated on the spectrum of the matrix. In particular, the $n$th residual of a Krylov method (with zero initial guess) ...

View answer
What kinds of problems lend themselves well to GPU computing?
22 votes

Problems which have a high arithmetic intensity and regular memory access patterns are typically easy(ier) to implement on GPUs, and perform well on them. The basic difficulty in having high ...

View answer
What is the purpose of using integration by parts in deriving a weak form for FEM discretization?
18 votes

Nothing stops you from doing that technically, but when you integrate by parts you get more flexibility with the solution space in that they need not have $H^2$ regularity (required for the non I.B.P ...

View answer
Do there exist low-storage Runge–Kutta methods with an order larger than four?
10 votes

This area has been fairly well researched, you may check e.g. Ketcheson's review of such methods: https://doi.org/10.1016/j.jcp.2009.11.006 which does contain some low-storage Runge-Kutta methods ...

View answer
Functional Programming and Scientific Computing
10 votes

I have maybe a unique perspective on this because I am a HPC practitioner with a scientific computation background as well as a functional programming language user. I don't want to equate HPC with ...

View answer
Solve linear system with Newton-Raphson method
Accepted answer
9 votes

Yes you can do this, and it will converge in one iteration regardless of the starting value. This is because each step of Newton's method involves solving a linear system with the Jacobian of the ...

View answer
Problems where Conjugate gradient works much better than GMRES
8 votes

One thing that CG has in its favor is that it's not minimizing the discrete $l^2$ norm for its residual polynomials (what GMRES does). It's minimizing a matrix-induced norm instead, and very often ...

View answer
How to efficiently invert $K \otimes M+I_T\otimes \Sigma$?
Accepted answer
7 votes

Generally there isn't a way to compute the inverse of a sum of Kronecker products. However, suppose there is a factor in common, let's say $I_T$ here and your sum is $$ A = K \otimes I_T + I_T \...

View answer
Optimal ODE method for fixed number of RHS evaluations
6 votes

There aren't many results in this direction because it is more difficult than just fixing accuracy, since stability considerations can often require you to pick time-steps that are smaller than you ...

View answer
How does gmres method iteration behave for this non-diagonalizable matrix?
5 votes

Estimates for GMRES convergence based on eigenvalue distribution often implicitly assume that the matrix is normal. Sometimes the convergence rate is still provable in an asymptotic sense in the non-...

View answer
Missing something fundamental about condition number estimation
Accepted answer
5 votes

This algorithm is most useful for two situations, which are related to each other in practice: You don't know the matrix entries explicitly, but instead can only compute matrix-vector products with ...

View answer
Choose a subset of $m$ columns that maximize $|A^T A|$?
5 votes

A related task to this is to find a subset of column vectors that are maximally linearly independent. Linear independence isn't exactly the same thing as asking for a large determinant, but if we can ...

View answer
Why does sparse linear algebra have a low arithmetic intensity?
5 votes

This really depends on the operations you are including in your question. If you took the sparse equivalent of any level 1 BLAS or level 2 BLAS algorithm, then yes they are memory bound (not compute ...

View answer
Building Gaussian-type quadrature schemes with Zernike polynomials
5 votes

Unfortunately the 1D reasoning will not directly carry into 2D and beyond. Depending on the domain for your integral, you will have considerable latitude in the node choice - but of course the ...

View answer
In what regime do the continuous and discontinuous Galerkin method become unstable for advection-diffusion systems?
Accepted answer
5 votes

In the advection dominated case the problem can develop physics which is invisible to a computational mesh that is too coarse (say by having elements which contain many wavelengths). This leads to ...

View answer
What's the most efficient way to compute the eigenvector of a dense matrix corresponding to the eigenvalue of largest magnitude?
5 votes

Power iteration is the simplest, but as mentioned above it would likely converge very slowly if the matrix is very non-normal. You get a "hump" phenomenon where the sequence appears to diverge for ...

View answer
Why the numerical solution of advection-dominant problem is challenging
4 votes

The difficulty is relative to something, in this case it is relative to diffusion dominated problems. Diffusion dominated aren't "easy" either, they have their own set of problems. I'll start with ...

View answer
Is it a good idea to use vector<vector<double>> to form a matrix class for high performance scientific computing code?
4 votes

I don't recommend it, but not because of performance issues. It will be a little less performant than a traditional matrix, which are usually allocated as a big chunk of contiguous data that is ...

View answer
What type of matrices is approximate inverse preconditioner $||I - AM||_F$ well suited for?
4 votes

I'm not an expert on this topic, but I have read a few papers concerning this in the hopes of finding a useful preconditioner for my problem. I think though that understanding of this topic is rather ...

View answer
Convergence/stagnation of BiCGStab(l)
4 votes

You might be familiar with the following paper already: http://link.springer.com/chapter/10.1007%2F978-3-642-22061-6_10 Problems which are highly indefinite and oscillatory are very difficult to ...

View answer
What is the meaning of "preasymptotic" and "superconvergent"?
4 votes

The meanings of those terms depend on context. Superconvergence is usually used to mean you are converging faster than the "optimal" rate, and occasionally this sort of weirdly fast convergence can ...

View answer
Reference request: Rigorous analysis of algorithms for PDE and ODE
4 votes

In addition to Jed's great recommendations (I can personally vouch for Brenner+Scott as a great intro finite elements book), an excellent book for the numerical solution of ODEs is Butcher: http://...

View answer
Is it possible to use Octave to learn MATLAB programming?
4 votes

Yes you absolutely can, I did. Much of the power out of MATLAB however does come with some of its hugely easy to use toolboxes and builtins which may or may not have equivalents in Octave. Also be ...

View answer
Runge-Kutta and Reusing Datapoints
4 votes

In general explicit Runge-Kutta methods of order $N$ require at least $N$ function evaluations, and there is absolutely no way to avoid this. Past $N=4$ they require more than $N$ function evaluations....

View answer
Is it possible to use BLAS if I have a function rather than a matrix?
4 votes

In a word: No. The reason MKL and other optimized libraries are so good is they are optimized to balance calculation versus memory use (i.e. your matrix is fully stored in RAM explicitly). If your ...

View answer
Sufficient conditions to ensure convergence of the conjugate gradient method
4 votes

I think diagonal dominance + symmetry gives positive semidefinite (definite if the dominance is strict). Since you know it's invertible that gives positive definite. Edit: The diagonal entries must ...

View answer
Does a new proposed method for solving Ax=b must beat matlab command 'A\b' be a successful method?
Accepted answer
3 votes

It's almost impossible to say whether a direct solver will outperform an iterative solver or vice-versa without knowing more specific information about the sparse matrix. The key problem with direct ...

View answer
Role of the numerical flux in DG-FEM
3 votes

Very loosely speaking there are two things most discretization techniques need in order to converge to the actual solution of your PDE as you increase their approximation quality, regardless if you're ...

View answer