25

Accuracy Trefethen and Schreiber wrote an excellent paper, Average-case Stability of Gaussian Elimination, which discusses the accuracy side of your question. Here are a few of its conclusions: "For QR factorization with or without column pivoting, the average maximal element of the residual matrix is $O(n^{1/2})$, whereas for Gaussian elimination it is $O(...


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

There are lots of modern finite element references, but I will just comment on a few books that I think are practical and relevant to applications, plus one containing more comprehensive analysis. Wriggers Nonlinear Finite Element Methods (2008) is a good general reference, but will be most relevant to those concerned with applications in structural ...


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, ...


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 ...


15

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,...


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

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

This answer partly responds to JackPoulson's comment (because it is long), and partly answers the question. Interval arithmetic is a computational procedure to give rigorous bounds on calculated quantities, only in the sense that the interval extension of a real-valued function over an interval encloses the image of that function over the same interval. ...


11

I'm surprised no one mentioned linear least squares problems, which occur frequently in scientific computing. If you want to use Gaussian elimination, you have to form and solve the normal equations, which look like: $$A^{T}Ax = A^{T}b,$$ where $A$ is a matrix of data points corresponding to observations of independent variables, $x$ is a vector of ...


11

For the second question, as a reader of Claes Johnson's book myself, I would say that you didn't miss much as a beginner in the finite element method, that book is pretty well-rounded with every aspect of the FEM except for the implementation. However, lots of developments have been made since the book published 20 years ago, like other people already ...


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

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-...


11

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$ ...


10

I'd probably start with Gil Strang's Introduction to Linear Algebra. It's best to get a solid foundation of the subject without proofs before moving on to a rigorous introduction, like learning calculus before studying real analysis. After you study Strang's book, if you're still interested in learning more about the rigor behind linear algebra, you could ...


10

Here is a 97-line example of solving a simple multivariate PDE using finite difference methods, contributed by Prof. David Ketcheson, from the py4sci repository I maintain. For more complicated problems where you need to handle shocks or conservation in a finite-volume discretization, I recommend looking at pyclaw, a software package that I help develop. ""...


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

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

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

Original answer (May 2012) As far as I am concerned, the Blaze library has not been publicly released. A link to the software as well as the license for its use should have both been in the paper. If you're interested in a modern, freely available, numerical linear algebra library that heavily leverages expression templates, I recommend Eigen. Update (...


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 ...


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