28
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-Fike Theorem.)
If $A$ is non-symmetric, then the eigenvalues of $A+E$ could be just anywhere in the complex plane. No bound can be formulated a priori. (...
27
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 symplectic manifold if it comes from one, and so if the system comes from a Hamlitonian system, then it will solve on some perturbed Hamiltonian trajectory. Some ...
22
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 advantage but the compute is done with single precision after converting to and from the IEEE half precision format.
https://software.intel.com/content/www/us/en/...
20
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 sufficiently distinguishing between two interconnected but different steps: creating a mathematical model of a physical process and solving it numerically. Let me ...
20
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 analysis (as distinguished from discrete
mathematics).
Numerical analysts are typically interested in proving mathematical results about their algorithms, ...
15
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.$$
13
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; usually, those have a slightly wider applicability for the same amount of computational effort.
I'll be discussing $\exp x-1$ for the rest of the answer, since $10^x-1=\exp(x\log 10)-1$.
To ...
13
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 method is two, the implemented order is one. This is due to the implementation of the first derivative as a one-sided difference quotient. Changing that to the ...
12
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 of an attractor, all points in some neighborhood of the attractor will converge to the attractor under the flow of the dynamical system as $t\to\infty$. However, ...
11
The Finite Element Method (FEM) is the parent method which has inspired many, many other methods and methods which are actually FEM but pretend not to be.
In the finite element method, "shape functions" are used to provide an approximation space so that the solution can be represented by a vector. In the classical FEM, these shape functions are polynomials.
...
11
This is the stochastic root-finding problem, as in The stochastic root-finding problem: Overview, solutions, and open questions.
11
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 page here.
Broad but shallow knowledge of all of the above is absolutely necessary for the team as a whole. Deep knowledge can be acquired as needed, or hired as ...
10
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 floating-point evaluation:
$$
{u_{0} = 2,\\
u_{1} = -4,\\
u_{n} = 111 - \frac{1130}{u_{n-1}} + \frac{3000}{u_{n-1}u_{n-2}},\>\>\>\>n \ge 2.}
$$
...
10
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 Taylor expansion plus a remainder term.
One can consider the Bramble-Hilbert lemma as a variation of Taylor's theorem, but it has different applications and is ...
10
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 timestep restrictions that are related to spatial refinement and choice of time integration method. Due to the diffusion term in 13a, the timestep must satisfy an ...
9
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)-x| \leq 2^{-23}|x|$ when $x>y>1$ and $x+y\neq\infty_{32}$ in 32-bit floating point arithmetic:
λ> import Data.SBV
λ> :set -XScopedTypeVariables
λ&...
9
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 changes in the input. Specifically, we have
$$ \underset{\epsilon \rightarrow 0_+}{\lim}\sup \left\{ \frac{1}{\epsilon} \left|\frac{s(x+\Delta x) - s(x)}{s(x)} \...
9
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 nonlinear function. In this case the Jacobian just equals $A$.
In other words: this is a little circular because it requires you to solve the system $Ax=b$ in the ...
9
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}{\log 2} - \frac{\log\log 2}{\log 2} + \sum_{n\geq 1}a_n x^n, $$
where $a_n$ are Taylor series coefficients of the l.h.s. after the log-singularity is subtracted. ...
8
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 very small, then these time scales are very different.
Stiff ODEs are difficult to solve because in order to solve them accurately, you have to resolve the ...
8
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 can modify the time derivatives accordingly. This means that you should have the differential equation $\frac{dy}{d\tau} = -f(-\tau, y)$ with $y(\tau = 0) = 0$ as ...
8
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 simulations. In between pure theory and expensive real-world experiments, we can now run simulations!
When it comes to these simulations, you may observe two types of ...
8
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 IEEE-754 half (10 bit mantissa, 5 bit exponent, 1 bit sign) but also bfloat16 (7 bit mantissa, 8 bit exponent, 1 bit sign) which favors dynamic range over ...
8
Amazing question which has a long answer, but I will try to be concise. In the context of Krylov subspace methods for general matrices, the eigenvalues of a non-symmetric matrix mean very little. In “Any nonincreasing convergence curve is possible for GMRES”, Greenbaum et al. show that any nonincreasing convergence curve is possible for GMRES independent of ...
7
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 observations give you a field $u_{\text{exp}}$; you'd like to find a value of $\mu$ for which
$\frac{1}{2}\iint(u - u_\text{exp})^2dx\hspace{2pt}dt$
is a minimum, ...
7
Both Mike's answer and Jed's one describe well the XFEM/FEM dichotomy and correctly point out that the most important area of application is 3D Fracture Mechanics, where you have a crack, i.e. a displacement discontinuity across a surface inside your domain.
Cracks are hard to model in classical FEM for two reasons:
The mesh has to be congruent across the ...
7
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" when you are learning these things in class, but when researchers choose them, they do so drawing on the combined experience of the field, as published in ...
7
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}\begin{bmatrix}x_{1}(t)\\
x_{2}(t)
\end{bmatrix}=0 \tag1
\end{equation}
where the natural frequency is $\omega \triangleq \sqrt{1-\eta^2/4}$. It is easy to verify ...
7
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 in three dimensions for a
vector function $F$ can be written
$$
\int_V \nabla \cdot {\bf F} dv = \int_S {\bf F} \cdot ds
$$
That is, an integral over the ...
7
Performing $k$ steps of GMRES uses $O(n k^2)$ time and $O(n k)$ memory. In other words, the algorithm gets more and more expensive with each additional iteration. In theory, the algorithm terminates in $k \le n$ steps (ignoring round-off), thereby yielding a worst-case cubic $O(n^3)$ time and quadratic $O(n^2)$ memory figure.
In practice, we manually ...
Only top voted, non community-wiki answers of a minimum length are eligible