19

So there is a ton to say about this, and we will actually be putting a paper out that tries to summarize it a bit, but let me narrow it down to something that can be put into a quick StackOverflow post. I will make one statement really early and keep repeating it: you cannot untangle the efficiency of a method from the efficiency of a software. The details ...


18

See https://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (sorry, I would have put that in a comment but I've registered just to post this so I can't post comments yet). But since I'm writing it as an answer, I'll also write the method: $$E=\frac{m_{00}+m_{11}}{2}; F=\frac{m_{00}-m_{11}}{2}; G=\...


17

As Jed Brown mentioned, the connection between gradient descent in nonlinear optimization and time stepping of dynamical systems is rediscovered with some frequency (understandably, since it's a very satisfying connection to the mathematical mind since it links two seemingly different fields). However, it rarely turns out to be a useful connection, ...


17

If you're doing celestial mechanics over long time scales, using a classical Runge-Kutta integrator will not preserve energy. In that case, using a symplectic integrator would probably be better. Boost.odeint also implements a 4th-order symplectic Runge-Kutta scheme that would work better for long time intervals. GSL does not implement any symplectic methods,...


16

Discretization of the continuum PDEs usually commits much more error than the finite precision. I find that about 90% of the people requesting higher precision have just been lazy with problem formulation and are trying to solve a problem using poor scaling, bad discretizations, or bad continuum modeling. The remaining 10% may have justifiably ill-...


15

This issue has become much more nuanced as the changes in architectures has shifted the HPC landscape. As Wolfgang Bangerth mentions one current longstanding view, I'll split my answer into basic definitiions and further details. Basic Definition A node refers to the physical box, i.e. cpu sockets with north/south switches connecting memory systems and ...


14

You can't beat an explicit formula. You can write down the formulas for the solution $x=A^{-1}b$ on a piece of paper. Let the compiler optimize things for you. Any other method will almost inevitably have if statements or for loops (e.g., for iterative methods) that will make your code slower than any straight line code.


13

Steve McConnell's book Code Complete, 2nd edition has an extensive bibliography discussing these issues from more of the standpoint of software developers than computational scientists. The book is starting to become a little dated, in that it's approaching a decade old, so it doesn't cover more recent testing methodologies like behavior-driven development. ...


13

While I haven't seen the exact formulation that you have written down here, I keep seeing talks in which people "rediscover" a connection to integrating some transient system, and proceed to write down an algorithm that is algebraically-equilavent to one form or another of an existing gradient descent or Newton-like method, and fail to cite anyone else. I ...


12

For the sake of notation, let's suppose that $f: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}$ (i.e., it's a vector-valued function that takes a vector as input and outputs a vector of the same size). There are two concerns: computational cost and numerical accuracy. Calculating the derivative $\mathrm{D}f(x)$ (the Jacobian matrix, $J(x)$, or $(\nabla f(x))^{T}...


12

Two examples of libraries that use modern C++ constructs: Both the eigen and armadillo libraries (linear algebra) use several modern C++ constructs. For instance, they use both expression templates to simplify arithmetic expressions and can sometimes eliminate some temporaries: http://eigen.tuxfamily.org http://arma.sourceforge.net/ http://hpac.rwth-...


12

Yes, this is the standard Aubin-Nitsche (or duality) trick. The idea is to use the fact that $L^2$ is its own dual space to write the $L^2$-norm as an operator norm $$\|u\|_{L^2} = \sup_{\phi\in L^2\setminus\{0\}} \frac{(u,\phi)}{\|\phi\|_{L^2}}.$$ We thus have to estimate $(u-u_h,\phi)$ for arbitrary $\phi\in L^2$. To do that, we "lift" $u-u_h$ to $H^1_0$ ...


12

Though it is a relatively rare situation when you actually have to calculate an inverse of the matrix, not all techniques were created equally. I would use the term badly-conditioned instead of ill-conditioned. For badly conditioned matrices, you might opt in the SVD-route to calculate the inverse: $$ A=U\Sigma V^H \implies A^{-1}=V\Sigma^{-1}U^H. $$ If ...


11

OPT++ is used internally by Dakota (Sandia), which is much more than an optimization library and is released under the LGPL. You should also take a look at TAO (ANL), released under a BSD-like license. An introduction to both OPT++ and TAO can be found here. Other alternatives are MOOCHO, NOMAD, and HOPSPACK, which are, as far as I know, also LGPL licensed. ...


11

I posted some benchmarks here: http://www.cecm.sfu.ca/~rpearcea/mgb.html These are for total degree orders. To solve systems you typically need to do more work. Timings are for a typical midrange desktop as of 2015 (Haswell Core i5 quad core). The fastest system on one core is Magma, which uses floating point arithmetic and SSE/AVX. Magma is the ...


10

I would suggest taking a look at Deal.II. It uses the STL, it's own iterators, shared pointers, etc. The various linear solvers can use the various matrices because of how it was designed. I haven't come across any use of move semantics, but that doesn't mean they aren't there. Here is a link.


10

For parabolic/elliptic PDE's, I highly recommend Beatrice Riviere's book: Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation. For hyperbolic PDE's and general (i.e. nonlinear) conservation laws, I recommend Hesthaven & Warburton's book: Nodal discontinuous Galerkin methods: algorithms, analysis, and ...


9

Usually one doubles the initial step until the Goldstein condition is violated or (in a feasible point method) the boundary is reached. Then one has a bracket. (If no such step exists, the objective function is unbounded below.) One can also use less conservative extrapolation procedures, but these require good tuning to be robust enough in a general purpose ...


9

In the 15 years that we have provided FEM software in the form of the deal.II project (http://www.dealii.org/), I don't think that we've ever had a genuine request to solve PDEs to higher accuracy than double precision. The reason is as Jed suggests in the other answer: The error one makes discretizing the PDE is much larger than the 16 digits of accuracy ...


9

You won't find one reference systematically covering the analysis of all the important methods for PDE. The field of discretization techniques for PDE is at least an order of magnitude larger than either topic you mentioned above. For any methods involving implicit solves, studying discretizations without also considering solution methods (e.g., associated ...


9

From skimming the table of contents to the book you listed, I'd say that computational books of that type for physics (or in my case, engineering, since that is my background) tend to sacrifice depth and quality of explanation for breadth. The best book of this ilk I can think of is probably Strang's Computational Science and Engineering, because he's a ...


9

@Pedro Gimeno "I doubt it can be any more robust than that." Challenge accepted. I noticed the usual approach is to use trig functions like atan2. Intuitively, there shouldn't be a need to use trig functions. Indeed, all the results end up as sines and cosines of arctans--which can be simplified to algebraic functions. It took quite a while, but I managed ...


9

Both GNU Scientific Library (GSL) (C) and Boost Odeint (C++) feature 8th order Runge-Kutta methods. Both are opensource, and under linux and mac they should be directly available from the package manager. Under windows, it will probably be easier for you to use Boost rather than GSL. GSL is published under the GPL license, and Boost Odeint under the ...


9

Since the matrix is so close to the identity, the following Neumann series will converge very rapidly: $$A^{-1} = \sum_{k=0}^\infty (I-A)^k$$ Depending on the accuracy required it might even be good enough to truncate after 2 terms: $$A^{-1} \approx I + (I - A) = 2I - A.$$ This might be slightly faster than a direct formula (as suggested in Wolfgang ...


8

Since the size of each type of atom is fixed, for a given level of accuracy the asymptotic cost is dominated by far field electrostatic interactions. These are $O(n)$ using multigrid and $O(n \log n)$ using FFTs. Thus the optimal complexity is $O(n)$ per time step as a function of the number of atoms simulated. The time step is also asymptotically ...


8

The best high-level overview that I know of is Trefethen and Bau. If I had to boil it down to a list, it would be (somewhat in pedagogical order): Dense $QR$ factorization Dense symmetric/Hermitian Eigenvalue Decomposition (EVD) Dense Singular Value Decomposition (SVD) The Conjugate Gradient Method (CG) Generalize Minimum Residual method (GMRES) Sparse ...


8

My personal favorite for linear structural mechanics and dynamics has not been mentioned yet: Finite Element Procedures, from KJ Bathe. If you have a structural engineering background this book is the best introduction to the FEM that I have seen. It discusses the formulation of structural elements in depth, the inf-sup condition, error estimation, and ...


8

Kate, you asked for a reference request as a high school student interested in getting started in computational science, which I think I can tackle fairly specifically. As long as you are comfortable (or getting comfortable) with Calculus, I think there are two great self-interest textbooks for you to go, depending on your interest and access to Python and ...


8

The GSL has a 2-by-2 SVD solver underlying the QR decomposition part of the main SVD algorithm for gsl_linalg_SV_decomp. See the svdstep.c file and look for the svd2 function. The function has a few special cases, isn't exactly trivial, and looks to be doing several things to be numerically careful (e.g., using hypot to avoid overflows).


8

I think you start by looking at something like FEniCS. Marie Rognes has a presentation with code examples and a paper discussing the theory and implementation. libMesh is supposed to be able to do something similar for 2-manifolds in 3-space, and so is deal.II, judging from this manuscript. Developers of deal.II and FEniCS answer questions on SciComp, and ...


Only top voted, non community-wiki answers of a minimum length are eligible