27
votes
Accepted
Conserving Energy in Physics Simulation with imperfect Numerical Solver
There are a few ways to conserve energy during ODE integration.
Method 1: Symplectic Integration
The cheapest way that is to use a symplectic integrator. A symplectic integrator solves the ODE on a ...
- 11.5k
27
votes
Accepted
Why do we usually not want the eigenvalues of non-symmetric matrices?
Stability under perturbations
Let $E$ be a perturbation such that $\|E\| \leq \varepsilon$.
If $A$ is symmetric, then the eigenvalues of $A+E$ are at a distance $\varepsilon$ from those of $A$. (Bauer-...
- 8,553
25
votes
Accepted
Is half precision supported by modern architecture?
Intel support for IEEE float16 storage format
Intel supports IEEE half as a storage type in processors since Ivy Bridge (2013). Storage type means you can get a memory/cache capacity/bandwidth ...
- 1,991
23
votes
Accepted
Scientific computing vs numerical analysis
Wikipedia gives a good definition
Numerical analysis is the study of algorithms that use numerical
approximation (as opposed to general symbolic manipulations) for the
problems of mathematical ...
- 17.6k
20
votes
Accepted
Use of machine learning in computational fluid dynamics
It's a long-running joke that CFD stands for "colorful fluid dynamics". Nevertheless, it is used -- and useful -- in a wide range of applications.
I believe your discontent stems from not ...
- 11.9k
17
votes
Accepted
Why is RK45 used as the "default" method for non-stiff ODEs rather than a multistep one?
First, let's establish that they are a good choice. The SciMLBenchmarks are probably the most comprehensive that there are as of right now for modern methods. This uses the vast number of methods ...
- 11.5k
15
votes
Accepted
Compute powers close to zero
Use the Taylor series expansion of $a^x$ about $x=0$ and evaluate a small number of terms for $a=10$:
$$a^x - 1 \approx x\log(a) + \frac{1}{2}x^2\log^2(a)+\cdots.$$
- 974
13
votes
Accepted
How to write integration tests for numeric simulation software?
I don't think you can avoid using a tolerance for floating-point comparisons. Error due to round-off, discretization, and so on using floating-point numbers is unavoidable.
What I typically do to ...
- 29.8k
13
votes
Compute powers close to zero
Personally, in the absence of a special function like expm1 (see also this scicomp.SE question), I'd use a Padé approximant instead of a Taylor/Maclaurin series; ...
- 3,075
13
votes
Accepted
Why is the central difference method dispersing my solution?
I'll write the equation short as $$\ddot x(t)+c\dot x(t)=a(t,x(t))$$ to separate the "easy" linear parts from the non-linear and forcing terms.
On the first method
The claimed order of the ...
- 3,511
12
votes
Runge-Kutta in the presence of an attractor
The problem you are encountering is likely not a consequence of your choice of algorithm, but in fact a consequence of the resulting dynamical system after applying time reversal. Per the definition ...
- 1,097
11
votes
Accepted
How important is learning hardware/architecture for scientific computing?
I haven't worked in quantum chemistry specifically, but I've worked in other areas where high performance is a correctness requirement (along with scientific accuracy), so I think we're on the same ...
- 351
10
votes
Accepted
Analytical convergent sequence and numerical divergent sequence
Jean-Michel Muller, et. al., "Handbook of Floating-Point Arithmetic 2nd ed.", Birkhäuser 2018, gives the following example due to Muller, specifically constructed to deliver incorrect results with ...
- 1,030
10
votes
What are the most important theorems in computational science?
You'll get everyone to give different answers to this question, and maybe that's alright. Here are some of my favorite ones:
Taylor's theorem that a function (of sufficient smoothness) equals its ...
- 50.2k
10
votes
How to solve a second order differential equation (diffusion) with boundary conditions using Python
I have found that I must keep the value of dt near dx or the results become unstable
This behavior you have noticed is known as the Courant–Friedrichs–Lewy (CFL) condition. Indeed, there are ...
- 934
10
votes
Why in scientific papers convergence of finite difference and finite volume schemes is tested using multiple norms ($l_1$, $l_2$ and $l_{\infty}$)?
In finite dimensional spaces (say, in $\mathbb R^n$), all norms are equivalent and as a consequence, if something converges with a specific rate in one norm, it also converges with the same rate in ...
- 50.2k
9
votes
How to write integration tests for numeric simulation software?
Geoff has already given an excellent overview, but I wanted to provide another real world look on it. In the deal.II project (http://www.dealii.org/) we run some 7,000 tests with every change in the ...
- 50.2k
9
votes
Accepted
Is it possible to proof a-b+b = a for all double floating-point numbers?
You can sometimes prove such results (or get counterexamples) using an SMT solver such as Z3 that supports floating point arithmetic. Here is a proof of a version of your theorem that says $|((x+y)-y)-...
- 11.4k
9
votes
Accepted
Poorly conditioned, easily evaluated sum for unit testing
The condition number of sum $s(x) = \sum_{j=1}^n x_j$ is given by
$$ \kappa(x) = \frac{\sum_{j=1}^n |x_j|}{|\sum_{j=1}^n x_j|} = \frac{s(|x|)}{|s(x)|}$$ and reflects the sums sensitivity to small ...
- 1,351
9
votes
Accepted
Solve linear system with Newton-Raphson method
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 ...
- 3,203
9
votes
Accepted
Accurate computation of Gauss-Kuzmin entropy
It's fairly easy to evaluate, to do this expand the logs in Taylor series in $x=(k+1)^{-2}$:
$$ \log_2(1-x) = \frac{-1}{\log 2}\sum_{m\geq1}\frac{x^{m}}{m}$$
$$ \log_2(-\log_2(1-x)) = \frac{\log x}{\...
- 11.4k
8
votes
Scientific computing vs numerical analysis
As someone who moved from Engineering to Scientific Computing during Grad school as an incidental need of the kind of work I was doing here are my two cents:
Numerical analysis would focus on the ...
- 259
8
votes
ODE $x''(t)+\eta x'(t)+x(t)=0$ with the $\eta$ extremely small
The problem you have is what is called "stiff": it has two time scales, namely one for the oscillation (which is ${\cal O}(1)$) and one for the damping (which is ${\cal O}(\eta^{-1})$). If $\eta$ is ...
- 50.2k
8
votes
Accepted
Going back in time in an initial value problem
This is technically still an IVP if you do an appropriate change of variables. Given your time is between $t \in [t^*, 0]$, make a new time variable $\tau = -t$ so that $\tau \in [0, -t^*]$ and you ...
- 3,673
8
votes
Accepted
What kind of a researcher am I?
Up until a couple of decades ago, science was based on two large pillars. Those were theory and actual physical experiments. It is an exciting time to see a third pillar arise with numerical ...
- 1,936
8
votes
Is half precision supported by modern architecture?
In my opinion, not very uniformly. Low precision arithmetic seems to have gained some traction in machine learning, but there's varying definitions for what people mean by low precision. There's the ...
- 4,296
7
votes
Use of machine learning in computational fluid dynamics
I think you are mixing a couple different ideas that are causing confusion. Yes, there are a wide variety of ways to discretize a given problem. Choosing an appropriate way may look like "voodoo" ...
- 4,567
7
votes
Accepted
Inverse advection-diffusion problem, solving for a drift coefficient with experimental data?
There is another approach called the adjoint method, which is commonly used in inverse problems for PDE and which is quite easy to generalize to other problems. This is going to be long.
Your ...
- 8,087
7
votes
ODE $x''(t)+\eta x'(t)+x(t)=0$ with the $\eta$ extremely small
Let us write the equation in state-space form
\begin{equation}
\frac{d}{dt}\begin{bmatrix}x_{1}(t)\\
x_{2}(t)
\end{bmatrix}+\begin{bmatrix}\eta/2 & -\omega\\
\omega & \eta/2
\end{bmatrix}\...
- 2,375
7
votes
Accepted
How do I integrate a function defined over an arbitrary area?
Instead of directly integrating over the area, it is often more convenient
to use the
divergence theorem to replace the area integral with an integral over the
boundary edges.
The divergence theorem ...
- 5,744
Only top scored, non community-wiki answers of a minimum length are eligible