All Questions

Filter by
Sorted by
Tagged with
13
votes
4answers
508 views

In floating point arithmetic, why does numerical imprecision result from adding a small term to a difference of large terms?

I have been reading the book Computer Simulation of Liquids by Allen and Tildesley. Starting on page 71, the authors discuss the various algorithms that are used to integrate Newton's equations of ...
13
votes
4answers
2k views

Is C slower than Fortran on the spectral norm shootout (using gcc, intel and other compilers)?

The conclusion here: How much better are Fortran compilers really? is that gfortran and gcc are as fast for simple code. So I wanted try something more complicated. I took the spectral norm shootout ...
13
votes
3answers
2k views

Mathematically, why does mass matrix / load vector lumping work?

I know that people often replace consistent mass matrices with lumped diagonal matrices. In the past, I've also implemented a code where the load vector is assembled in a lumped fashion rather than ...
13
votes
4answers
415 views

Testing suites for numerical applications in C++?

Recently, I've been pushing my group to include more testing when writing their code. There were several major bugs that took a lot longer to catch than probably speaking was necessary, because we ...
13
votes
3answers
9k views

Understanding the “rate of convergence” for iterative methods

According to Wikipedia the rate of convergence is expressed as a specific ratio of vector norms. I'm trying to understand the difference between "linear" and "quadratic" rates, at different points of ...
13
votes
2answers
5k views

What is the purpose of the test function in Finite Element Analysis?

In the wave equation: $$c^2 \nabla \cdot \nabla u(x,t) - \frac{\partial^2 u(x,t)}{\partial t^2} = f(x,t)$$ Why do we first multiply by a test function $v(x,t)$ before integrating?
13
votes
3answers
3k views

Is the Thomas algorithm the fastest way to solve a symmetric diagonally dominant sparse tridiagonal linear system

I am wondering if the Thomas algorithm is the fastest way (provably?) to solve a symmetric diagonally dominate sparse tridiagonal system in terms of algorithmic complexity (not looking for ...
13
votes
5answers
2k views

C++ or Python for a development of CFD library

What would you say would be the advantages/disadvantages of two approaches to coding a general (finite volume, fem, dg) library for Computational Continuum Mechanics? This is how I see things right ...
13
votes
1answer
4k views

Strong vs. weak solutions of PDEs

The strong form of a PDE requires that the unknown solution belongs in $H^2$. But the weak form requires only that the unknown solution belongs in $H^1$. How do you reconcile this?
13
votes
3answers
1k views

Is there any benefit to compiling LAPACK from source versus installing the prebuilt package from Ubuntu?

I know that ATLAS is able to optimize itself for the machine it is compiled on and thus maximum benefits are found by compiling from source. Is there any benefit to compiling LAPACK from source? It ...
13
votes
2answers
2k views

Does the “cofactor technique” for inverting a matrix have any practical significance?

The title is the question. This technique involves using the "matrix of cofactors", or "adjugate matrix", and gives explicit formulae for the components of the inverse of a square matrix. It is not ...
13
votes
3answers
1k views

SVD for finding the largest eigenvalue of a 50x50 matrix — am I wasting significant amounts of time?

I've got a program that computes the largest eigenvalue of many real symmetric 50x50 matrices by performing singular-value decompositions on all of them. The SVD is a bottleneck in the program. Are ...
13
votes
4answers
346 views

How to create a random 3D domain representing a plant's root structure?

I would like to model laminar flow of water from roots to the stem of a plant. At the very end of the roots, the tubes vary from millimeter to centimeter scale in diameter and length. As we get closer ...
13
votes
4answers
6k views

FLOP counting for library functions

When evaluating the number of FLOPs in a simple function, one can often just go down the expression tallying basic arithmetic operators. However, in the case of mathematical statements involving even ...
13
votes
1answer
926 views

What are the guidelines for conducting computational experiments?

Physics, biology, chemistry etc. have different sets of rules for making experiments: what events are considered relevant, how to avoid contamination of samples, how to create and fix a reproducing ...
13
votes
4answers
986 views

Estimating hardware error probability

Say I run a supercomputer computation on 100k cores for 4 hours on http://www.nersc.gov/users/computational-systems/edison/configuration, exchanging about 4 PB of data over the network and performing ...
13
votes
1answer
4k views

Pressure as a Lagrange Multiplier

In the incompressible Navier-Stokes equations, \begin{align*} \rho\left(\mathbf{u}_t + (\mathbf{u} \cdot \nabla)\mathbf{u}\right) &= - \nabla p + \mu\Delta\mathbf{u} + \mathbf{f}\\ \nabla\cdot\...
13
votes
2answers
12k views

Confusion about Armijo rule

I have this confusion about Armijo rule used in line search. I was reading back tracking line search but didn't get what this Armijo rule is all about. Can anyone elaborate what Armijo rule is? The ...
13
votes
6answers
5k views

What is a common file/data format for a mesh (for FEM)?

I'm developing an FEM simulation. For early testing, I will use simple self-written mesher and visualisation of the mesh graph. But I want to prepare my program to use data generated by an existing ...
13
votes
5answers
1k views

On Finding Open Source Projects To Contribute To

This question has been asked a billion times on Stackoverflow however, the focus has always been on Non-Numerical Coding. I am looking for a project to contribute to within the confines of Numerical ...
13
votes
4answers
336 views

Rapidly determining whether or not a dense matrix is of low rank

In a software project that I'm working on, certain computations are vastly easier for dense low-rank matrices. Some problem instances involve dense low-rank matrices, but they're given to me in full, ...
13
votes
2answers
5k views

Why does Matlab's integral outperform integrate.quad in Scipy?

I am experiencing some frustration over the way matlab handles numerical integration vs. Scipy. I observe the following differences in my test code below: Matlab's version runs on average 24 times ...
13
votes
3answers
2k views

Are HDF5 files suitable for git revision control?

I am not familiar with the file format used in HDF5, but I am wondering if HDF5 files are suitable for revision control with git (or for example Mercurial or Subversion)? I guess what I mean is: are ...
13
votes
1answer
1k views

Can the advection equation with variable velocity be conservative?

I am trying to understand the advection equation with variable velocity coefficient a bit better. In particular I don't understand how the equation can be conservative. The advection equation, $$ \...
13
votes
1answer
7k views

How does density functional theory scale with system size?

Theoretically, how does the time to do a density functional theory (DFT) calculation scale with the number of electrons? I'm interested in "typical" DFT implementations such as VASP, ABINIT, etc., not ...
13
votes
2answers
699 views

Design patterns in writing numerical software in C++

I'm looking for resources on design patterns and principles for numerical software, potentially but not necessarily with a focus on object-oriented approaches to numerical codes. I am aware of the ...
13
votes
3answers
4k views

Single versus double floating-point precision

Single precision floating point numbers take up half the memory and on modern machines (even on GPUs it seems) operations can be done with them at almost twice the speed compared to double precision. ...
13
votes
1answer
1k views

Why is leapfrog integration symplectic and RK4 not, if the latter is more accurate?

In a system where energy theoretically should be conserved, the most accurate simulation would conserve energy (as well as giving accurate positions, velocities and etc). RK4 is more accurate than ...
13
votes
1answer
383 views

CFD: Does order of time stepping scheme affect steady state solution? If yes why?

I am trying to solve Ideal MHD equation using semi discrete methods, ENO spatial reconstructions and TVD RK time stepping. I am getting different steady state solutions with different temporal order. ...
13
votes
5answers
2k views

Calculation of the sparsity structure for finite element matrices

Question: What methods are available to accurately and efficiently calculate the sparsity structure of a finite element matrix? Info: I'm working on a Poisson Pressure Equation solver, using Galerkin'...
13
votes
1answer
7k views

Understanding the Wolfe Conditions for an Inexact line search

According to Nocedal & Wright's Book Numerical Optimization (2006), the Wolfe's conditions for an inexact line search are, for a descent direction $p$, Sufficient Decrease: $f(x+\alpha p)\le f(x)+...
13
votes
3answers
5k views

Memory usage in fortran when using an array of derived type with pointer

In this sample program I'm doing the same thing (at least I think so) in two different ways. I'm running this on my Linux pc and monitoring the memory usage with top. Using gfortran I find that in the ...
13
votes
2answers
1k views

Alternatives to von neumann stability analysis for finite difference methods

I'm working on solving the coupled one-dimensional poroelasticity equations (biot's model), given as: $$-(\lambda+ 2\mu) \frac{\partial^2 u}{\partial x^2} + \frac{\partial p}{\partial x} = 0$$ $$\...
13
votes
1answer
1k views

What are possible methods to solve compressible Euler equations

I would like to write my own solver for compressible Euler equations, and most importantly I want it to work robustly in all situations. I would like it to be FE based (DG is ok). What are the ...
13
votes
2answers
2k views

Compute all eigenvalues of a very big and very sparse adjacency matrix

I have two graphs with nearly n~100000 nodes each. In both graphs, each node is connected to exactly 3 other nodes so the adjacency matrix is symmetric and very sparse. The hard part is I need all ...
13
votes
1answer
927 views

How is Krylov-accelerated Multigrid (using MG as a preconditioner) motivated?

Multigrid (MG) may be used to solve a linear system $Ax=b$ by constructing an initial guess $x_0$ and repeating the following for $i=0,1..$ until convergence: Compute the residual $r_i = b-Ax_i$ ...
13
votes
3answers
302 views

Computing slightly oscillatory series to high precision?

Suppose I have the following interesting function: $$ f(x) = \sum_{k\geq1} \frac{\cos k x}{k^2(2-\cos kx)}. $$ It has some unpleasant properties, like its derivative not being continous at rational ...
13
votes
2answers
947 views

Impose the compatibility conditions for mixed finite elements method in Stokes equation

$\newcommand{\v}[1]{\boldsymbol{#1}}$ Suppose we have following Stokes flow model equation: $$ \tag{1} \left\{ \begin{aligned} -\mathrm{div}(\nu \nabla \v{u}) + \nabla p &= \v{f} \\ \mathrm{div} \...
13
votes
3answers
835 views

What are the basic principles behind generating a moving mesh?

I am interested in implementing an moving mesh for an advection-diffusion problem. Adaptive Moving Mesh Methods gives a good example of how to do this for Burger's equation in 1D using finite-...
13
votes
1answer
1k views

Why is pinning a point to remove a null space bad?

A Poisson equation with all Neumann boundary conditions has a single constant dimensional null space. When solving via a Krylov method, the null space can be removed either by subtracting the mean of ...
13
votes
2answers
697 views

Why can ill-conditioned linear systems be solved precisely?

According to the answer here, large condition number (for linear system solving) decreases the guaranteed number of correct digits in the floating point solution. Higher order differentiation matrices ...
13
votes
2answers
9k views

Periodic boundary condition for the heat equation in ]0,1[

Let us consider a smooth initial condition and the heat equation in one dimension : $$ \partial_t u = \partial_{xx} u$$ in the open interval $]0,1[$, and let us assume that we want to solve it ...
13
votes
1answer
354 views

How to construct well-balanced finite volume and discontinuous Galerkin methods for hyperbolic PDEs with source terms?

Source terms, such as those due to bathymetry in the shallow water equations, need to be integrated in a special way in order to preserve physical steady states. Is there a general way to construct ...
13
votes
3answers
649 views

Fitting Implicit Surfaces to Oriented Point Sets

I have a question regarding quadric fit to a set of points and corresponding normals (or equivalently, tangents). Fitting quadric surfaces to point data is well explored. Some works are as follows: ...
13
votes
1answer
4k views

PDE solvers for Drift-diffusion and related models

I'm trying to simulate basic semiconductor models for pedagogical purposes--starting from the Drift-diffusion model. Although I don't want to use an off-the-shelf semiconductor simulator--I'll be ...
13
votes
2answers
521 views

Verification in Eigenvalue problems

Let us start with a problem of the form $$(\mathcal{L} + k^2) u=0$$ with a set of given boundary conditions (Dirichlet, Neumann, Robin, Periodic, Bloch-Periodic). This corresponds with finding the ...
13
votes
2answers
484 views

Which time-integration methods should we use for hyperbolic PDEs?

If we employ the Method of Lines for discretization (separate time and space discretization) of hyperbolic PDEs we obtain after spatial discretization by our favorite numerical method (fx. Finite ...
13
votes
3answers
1k views

Best practice for storing hierarchical simulation data

TL,DR What is the accepted best practice in scientific computing circles for storing large quantities of hierarchically structured data? For example, SQL does not play nicely with large sparse ...
13
votes
1answer
1k views

Can an approximated Jacobian with finite differences cause instability in the Newton method?

I have implemented a backward-Euler solver in python 3 (using numpy). For my own convenience and as an exercise, I also wrote a small function that computes a finite difference approximation of the ...
13
votes
1answer
355 views

Is there a tool out there that can generate interval extensions of Fortran (or C) functions by parsing Fortran (or C) code?

Case studies in my PhD thesis require that I have interval extensions of Fortran subroutines in CHEMKIN-II (apologies for the link; it's the best one I could find for a package no longer distributed ...

15 30 50 per page
1
5 6
7
8 9
178