26
Your question is asking about model Verification. You can find numerous resources on methods and standards by searching for Verification and Validation (Roache 1997, 2002, 2004, Oberkampf & Trucano 2002, Salari & Knupp 2000, Babuska & Oden 2004), as well as the broader topic of Uncertainty Quantification. Rather than elaborate on methods, I would ...
26
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. ...
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 comment ...
19
No such standards exist, as reliable error estimates often cost much more than the approximate calculations.
Basically there are four kinds of error estimates:
(i) Theoretical analyses proving that a numerical method is numerically stable. This doesn't really give an error bar as the analysis only guarantees that the error made is not worse than a ...
18
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, ...
13
Sort of. There are theoretical error bounds that have been derived by numerical analysts that are usually overestimates, and may not be as useful in practice, because they may involve information that is difficult to obtain for problems in practice. A good example would be the bounds on numerical errors in the solution of ordinary equations, which you can ...
13
The formula
$$\mathrm{logsum}(x,y)=\max(x,y)+\mathrm{log1p}(\exp(-\operatorname{abs}(x-y))$$
should be numerically stable. It generalizes to a numerically stable computation of
$$\log \sum_i e^{x_i} = \xi+ \log\sum_i e^{x_i-\xi},~~~\xi=\max_i x_i$$
In case the logsum is very close to zero and you want high relative accuracy, you can probably use
$$\...
13
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.
...
10
This is the stochastic root-finding problem, as in The stochastic root-finding problem: Overview, solutions, and open questions.
10
I have a suspicion that your Python expression for the right-hand side is doing integer division, not floating-point division. As a result, ((l1-l2)**2)/(d**2) is being evaluated as zero, and the term inside the square root is one.
In fact, you forgot the reciprocal in your right-hand expression as well, but that's not the first problem...
In MATLAB:
>&...
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 ...
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 ...
8
$BV$ spaces are exactly what you want in many inverse problem. The point is that in many inverse problems you try to determine a function $q(x)$ that describes the internal properties of a body -- say, the water content (MRI), the density (X-ray) or elastic coefficients (ultrasound). A good approximation is that this function $q(x)$ is constant (or, at least,...
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 ...
7
In addition to the other answers, there are a few additional points to consider.
Numerical discretization errors, or at least the order of the schemes, can be determined analytically. The discussion of these errors may be omitted from papers if they use a commonly known scheme.
Grid refinement studies where the same problem, usually something simple, is ...
7
All you can compare in such cases are the statistics of your solution: averages, higher moments, heat fluxes across the boundary, and other integral quantities. Take a look at one of the many papers discussing turbulence models for the Navier-Stokes equations, for example: they're full to the brim with plots of power spectra, entalpies, entropies, ...
7
As Michael C. Grant pointed out, the problem is that the division operation is integer division, not floating point division. In versions of Python before 3, dividing two integers using / rounds down to an integer. This behavior can be changed by adding
from __future__ import division
at the top of your script. Note that this declaration changes the ...
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
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
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
This is a way of writing the flow of a vector field. In other words, if you have an ordinary differential equation (ODE) given by $\dot{x}(t) = v(x(t))$ with initial value $x(0) = x_0$, you could formally write its solution as $x(t) = \mathrm{e}^{t v} x_0$. The formal solution coincides with the actual solution when $v$ is a linear mapping (think of the case ...
6
FEM is a subset of XFEM. XFEM is a methodology for enriching finite-element spaces to handle problems with discontinuities (such as fracture). With classical FEM, attaining similar accuracy typically requires complicated conformal meshing and adaptive refinement, where as XFEM does it with a single mesh, moving that geometric complexity into the elements (...
6
Reformulating (1) and (2), the system reads
$$
\begin{pmatrix}
\mathcal{L}&\\
&\mathcal{L}
\end{pmatrix}
\begin{pmatrix}
W\\
h_2
\end{pmatrix}
=
\kappa
\begin{pmatrix}
-\mathcal{L} &\mathcal{M}\\
\mathcal{M} & -\mathcal{L}
\end{pmatrix}
\begin{pmatrix}
W\\
h_2
\end{pmatrix}
$$
where the calligraphic symbols are composed of some linear ...
Only top voted, non community-wiki answers of a minimum length are eligible
