All Questions

Filter by
Sorted by
Tagged with
23
votes
3answers
2k views

Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation?

Is there an $O(n^3+n^2 k)$ method to solve $k$ linear systems of the form $(D_i + A) x_i = b_i$ where $A$ is a fixed SPD matrix and $D_i$ are positive diagonal matrices? For example, if each $D_i$ is ...
23
votes
2answers
10k views

Libraries for solving sparse linear systems

There are a number of different libraries out there that solve a sparse linear system of equations, however I'm finding it difficult to figure out what the differences are. As far as I can tell there ...
23
votes
10answers
16k views

Fast, lightweight C++ tensor library for dimension-agnostic code

I am looking for a C++ tensor library that supports dimension-agnostic code. Specifically, I need to perform operations along each dimension (up to 3), e.g. calculating a weighted sum. The dimensions ...
23
votes
4answers
3k views

Algorithms for (adaptive?) function plotting

I am looking for algorithms to draw standard 2d-graphs for functions that may or may not have singularities. The purpose is to write a "Mini-CAS", so I have no a priori knowledge of the types of ...
23
votes
2answers
4k views

A good finite difference for the continuity equation

What would be a good finite difference discretization for the following equation: $\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho u\right)=0$? We can take the 1D case: $\frac{\partial \...
23
votes
3answers
5k views

Using unconventional programming languages for scientific computation [closed]

Note: the following post may include controversial opinions, so please note that they are only my opinions, and not intended to offend anyone. I'm being programming in some form or the other since ...
23
votes
3answers
19k views

Intel Fortran Compiler: tips on optimization at compilation

I will start with my personal experience in our lab. Back in the ifort 9 and 10 days, we used to be quite aggressive with the optimizations, compiling with -O3 and processor specific flags (-xW -xSSE4....
22
votes
6answers
16k views

How can I numerically differentiate an unevenly sampled function?

Standard finite difference formulas are usable to numerically compute a derivative under the expectation that you have function values $f(x_k)$ at evenly spaced points, so that $h \equiv x_{k+1} - x_k$...
22
votes
3answers
15k views

Recommendation for Finite Difference Method in Scientific Python

For a project I am working on (in hyperbolic PDEs) I would like to get some rough handle on the behavior by looking at some numerics. I am, however, not a very good programmer. Can you recommend ...
22
votes
2answers
2k views

What simple methods are there for adaptively sampling a 2D function?

I have a two-dimensional function $f(x,y)$ whose values I would like to sample. The function is very expensive to compute and it has a complex shape, so I need to find a way to get the most ...
21
votes
9answers
39k views

Basic explanation of shape function

I just started studying FEM in a more structured basis compared to what I used to do during my undergraduate courses. I am doing this because, despite the fact that I can use the "FEM" in commercial (...
21
votes
3answers
8k views

Why should non-convexity be a problem in optimization?

I was very surprised when I started to read something about non-convex optimization in general and I saw statements like this: Many practical problems of importance are non-convex, and most non-...
21
votes
4answers
8k views

F2Py with allocatable and assumed shape arrays

I would like to use f2py with modern Fortran. In particular I'm trying to get the following basic example to work. This is the smallest useful example I could ...
21
votes
4answers
70k views

How to determine the amount of FLOPs my computer is capable of

I would like to determine the theoretical number of FLOPs (Floating Point Operations) that my computer can do. Can someone please help me with this. (I would like to compare my computer to some ...
21
votes
3answers
2k views

Parallel I/O options, in particular parallel HDF5

I have an application that can be trivially parallelized but its performance is to a large extent I/O bound. The application reads a single input array stored in a file that is typically 2-5 GB in ...
21
votes
1answer
2k views

How does the performance of Python/Numpy array operations scale with increasing array dimensions?

How do Python/Numpy arrays scale with increasing array dimensions? This is based on some behaviour I noticed while benchmarking Python code for this question: How to express this complicated ...
21
votes
1answer
3k views

Diagonal update of a symmetric positive definite matrix

$A$ is a $n \times n$ symmetric positive definite (SPD) sparse matrix. $G$ is a sparse diagonal matrix. $n$ is large ($n$ >10000) and the number of nonzeros in the $G$ is usually 100 ~ 1000. $A$ has ...
20
votes
4answers
9k views

Why are higher-order Runge–Kutta methods not used more often?

I was just curious as to why high-order (i.e. greater than 4) Runge–Kutta methods are almost never discussed/employed (at least to my knowledge). I understand it requires greater computational time ...
20
votes
3answers
3k views

What guidelines should I use when searching for good preconditioning methods for a specific problem?

For the solution of large linear systems $Ax=b$ using iterative methods, it is often of interest to introduce preconditioning, e.g. solve instead $M^{-1}(Ax=b)$, where $M$ is here used for left-...
20
votes
4answers
5k views

What is the best way to find discontinuities of a black-box function?

It was suggested that this might be a better place for this question than Mathematics Stack Exchange where I asked it before. Suppose one has a black-box function which can be evaluated anywhere (...
20
votes
6answers
3k views

Analyzing Numerical Error in C++ Function

Suppose that I have a function that takes as input several floating-point values (single or double), does some computation, and produces output floating-point values (also single or double). I am ...
20
votes
4answers
3k views

Is half precision supported by modern architecture?

I am new to computer science and I was wondering whether half precision is supported by modern architecture in the same way as single or double precision is. I thought the 2008 revision of IEEE-754 ...
20
votes
1answer
2k views

BDF vs implicit Runge Kutta time stepping

Are there any reasons for why one should choose high order implicit Runge Kutta (IMRK) over BDF time stepping? BDF seems much easier to me as $q$ stage IMRK needs $q$ linear solves per time step. ...
20
votes
6answers
14k views

How to get started with OpenFOAM for CFD

I'm looking at using OpenFOAM for solving basic internal flows in CFD. What is the best way to get started, and could anyone please point me to a good online reference to go to with any questions I ...
20
votes
2answers
2k views

Algorithms for a many-to-many generalized assignment problem

I can't seem to find any literature on algorithms which can be used to solve a many-to-many generalized assignment problem (GAP), i.e. models where not only can more tasks be assigned to one agent, ...
19
votes
6answers
44k views

Python vs FORTRAN

Which one is better: FORTRAN or Python? And I guess that in both cases you need Gnuplot, am I right? I'm working on a Windows machine at the moment. I'd like to use it to get numerical solutions for ...
19
votes
6answers
1k views

How do I write dimensionally agnostic code?

I often find myself writing very similar code for one, two, and three dimensional versions of a given operation/algorithm. Maintaining all of these versions can become tedious. Simple code ...
19
votes
3answers
14k views

Binary vs. ASCII file size

I need to write some data from a computation, that will be read later by Paraview (.vtu or vtk file). When it comes to file size , should I go for the ASCII format or the Binary format ?
19
votes
3answers
30k views

Euclidean distance in Octave

I would like to know if there is a quick way to compute the Euclidean distance of two vectors in Octave. It seems that there is no special function for that, so should I just use the formula with <...
19
votes
4answers
2k views

Numeric Quadrature with Derivatives

Most numerical methods for quadrature treat the integrand as a black-box function. What if we have more information? In particular, what benefit, if any, can we derive from knowing the first few ...
19
votes
3answers
836 views

Is it well known that some optimization problems are equivalent to time-stepping?

Given a desired state $y_0$ and a regularization parameter $\beta \in \mathbb R$, consider the problem of finding a state $y$ and a control $u$ to minimize a functional \begin{equation} \frac{1}{2} \...
19
votes
2answers
12k views

What is pseudo time-stepping?

While reading some literature on PDE solvers I came across the term pseudo time-stepping today. It seems to be a common term, however I failed to find a good definition or an introductionary article ...
19
votes
1answer
8k views

What is the general idea of Nitsche's method in numerical analysis?

I know that the Nitsche's method is a very attractive methods since it allows to take into account Dirichlet type boundary conditions or contact with friction boundary conditions in a weak way without ...
19
votes
5answers
3k views

Parallel Scientific Computation Software Development Language?

I want to develop a parallel scientific computation software from scratch. I want some thoughts on which language to start. The program involves reading/writing data to txt files and doing heavy ...
19
votes
6answers
1k views

What is the best way to do reproducible research if you need proprietary libraries?

Reproducible research in computation aims to make the code needed to generate the results in a computational paper available to other researchers so that they can run this code to reproduce the ...
19
votes
3answers
859 views

Strategies for unit testing and test-driven development

I'm a huge advocate of test-driven development in scientific computing. It's utility in practice is just staggering, and really alleviates the classic troubles that code developers know. However, ...
19
votes
2answers
2k views

How to determine if a numerical solution to a PDE is converging to a continuum solution?

The Lax equivalence theorem states that consistency and stability of a numerical scheme for a linear initial value problem is a necessary and sufficient condition for convergence. But for nonlinear ...
19
votes
4answers
4k views

Is Fortuna or Mersenne Twister preferable as an algorithmic RNG?

A recent answer mentioned the use of Fortuna or Mersenne Twister Random Number Generators (RNGs) to seed a Monte Carlo simulation. I hadn't heard of Fortuna before so I looked it up - looks like it is ...
19
votes
4answers
3k views

For which statistical methods are GPUs faster than CPUs?

I have just installed a Nvidia GT660 graphic card on my desktop and, after some struggle, I manage to interface it with R. I have been playing with several R packages that use GPUs, especially ...
19
votes
1answer
650 views

Difficulty with Spectral Method using Chebyshev Polynomials

I am having a bit of difficulty in trying to understand a paper. The paper uses spectral method to solve for an eigenvalue that comes from a system of coupled ODEs. I will write out only one equation ...
18
votes
5answers
6k views

State of the Mac OS in Scientific Computing and HPC

Back towards the dawn of OS X, there seemed to be a great deal of hubbub, at least in the Mac world (I was nowhere near scientific computing at the time) about the Mac OS as a platform for scientific ...
18
votes
3answers
2k views

Log-log parallel scaling/efficiency plots

A lot of my own work revolves around making algorithms scale better, and one of the preferred ways of showing parallel scaling and/or parallel efficiency is to plot the performance of an algorithm/...
18
votes
6answers
6k views

How does one determine the point group of a molecule?

You've managed to finally find out how the atoms are spatially arranged on your newly-discovered molecular entity. Through, say, spectroscopic means, you are now in possession of a bunch of atom ...
18
votes
4answers
4k views

Is there a general-purpose library for structured grid adaptive mesh refinement?

Adaptive mesh refinement (AMR) is a common technique for dealing with the problem of widely varying spatial scales in the numerical solution of PDEs. What general-purpose libraries exist for AMR on ...
18
votes
3answers
4k views

C++ Best practices for dealing with many constants, variables in scientific codes

I am developing a code to simulate fluid flow with biological substances present in the flow. This involves the standard Navier-Stokes equations coupled to some additional biological models. There are ...
18
votes
5answers
738 views

20% performance penalty for a nice software design

I'm writing a small library for sparse matrix computations as a way to teach myself to make the best use of object-oriented programming. I've worked really hard on having a nice object model, where ...
18
votes
4answers
607 views

Which methods can ensure that physical quantities remain positive throughout a PDE simulation?

Physical quantities like pressure, density, energy, temperature, and concentration should always be positive, but numerical methods sometimes compute negative values during the solution process. This ...
18
votes
4answers
5k views

The definition of stiff ODE system

Consider an IVP for ODE system $y'=f(x,y)$, $y(x_0)=y_0$. Most commonly this problem is considered stiff when Jacobi matrix $\frac{\partial f}{\partial y}(x_0,y_0)$ has both eigenvalues with very ...
18
votes
2answers
8k views

Null-space of a rectangular dense matrix

Given a dense matrix $$A \in R^{m \times n}, m >> n; max(m) \approx 100000 $$ what is the best way to find its null-space basis within some tolerance $\epsilon$? Based on that basis can I then ...
18
votes
1answer
2k views

How can wavelets be applied to PDE?

I would like to learn how wavelet methods can be applied to PDE, but unfortunately I do not know a good resource to learn about this topic. It seems that many introductions to wavelets focus on ...

15 30 50 per page