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

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/...


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, ...


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

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

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

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. ...


6

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 ...


5

Yes, the problem with mixed boundary conditions is well posed. What's not clear to me is this: Why do you approximate the derivative via the two-sides approximation? Shouldn't it be enough to just the following? $$ \frac{f_6 - f_{6-\Delta x}}{\Delta x} = 0. $$ In your animations, what is the size of $\Delta x$? Your curve looks very smooth, which can be ...


5

The problem that you are seeing is a well known problem of these basis functions: $$\phi_{m,n}(x,y) = x^{m}y^{n}$$ I quote from here: The reason why the coefficient matrix is nearly singular and ill-conditioned is that our basis functions $\Psi_{i}(x) = x^{i}$ are nearly linearly dependent for large i. That is, $x^{i}$ and $x^{i+1}$ are very close for i ...


5

Would a decomposition of the form $A = XX^T$ suffice? This would be enough, e.g., if the end goal is sampling from the Gaussian distribution with this given covariance. If so, you can use the following formula, which is quite similar to your approximation: $$X = D^{1/2} + \frac{\sqrt{u^T D^{-1} u+1}-1}{u^T D^{-1} u} u u^T D^{-1/2}$$ This follows from ...


5

You may have seen depictions of ODE dynamics for the scalar equation $x'(t)=f(t,x(t))$ in the form of a $t,x$ plot in which one plots little arrows in the entire $t,x$ plane. A trajectory starting from a particular $t,x_0$ value then follows the arrows it encounters along its way. Now think of an attractor, i.e., a curve in $t,x$ space to which all of these ...


5

The second formula is just the velocity verlet, and it's correct but if you adapt time steps then it's not symplectic. In a separate answer I describe in quite detail that symplecticness is a global property of the integration, it's not a stepwise property. Because of this, local error estimates and local changes of $\Delta t$ do not necessarily preserve the ...


5

As pointed out in the comments, the cost of evaluating $f$ is critical, and in most practical cases will be the dominant cost. Lets suppose it takes $C$ operations to evaluate $f$. For nontrivial functions, $C$ will be at least $\mathcal{O}(N)$ just from using the $N$ arguments. Also pointed out in the comments, the finite difference formula you have does ...


5

There is no need for numerical computation here. First, $T(q)$ is a well-known function, the logarithmic integral. Repeated integration by parts gives an asymptotic expansion $$\mathrm{Li}(q) = \frac{q}{\log q}\sum_{k=0}^{K-1} \frac{k!}{\log^k q} + O\left(\frac{q}{\log^{K+1}q}\right).$$ There's also a fairly rapidly convergent representation due to ...


5

The accepted answer provides an overview. I'll add a few more details about support in NVIDIA processors. The support I'm describing here is 16 bit, IEEE 754 compliant, floating point arithmetic support, including add, multiply, multiply-add, and conversions to/from other formats. Maxwell (circa 2015) The earliest IEEE 754 FP16 ("binary16" or "half ...


4

Crank-Nicolson is a very good classical approach for parabolic PDE like the heat transfer PDE to which it was originally applied. It is relatively easy to understand and implement so it is often presented in basic courses on numerical methods for PDE. pdepe is also very well-suited to this class of PDE (the second "p" in pdepe stands for parabolic). It has ...


4

Parareal is fairly effective for parabolic problems, but doesn't work well for hyperbolic problems. So it wouldn't be effective for the advection part, but could be effective for the diffusion part. Making it work well for hyperbolic problems is an area of ongoing research. Note that the parallel efficiency of parareal is quite poor, so if you haven't yet ...


4

It depends what you consider as an analytical solution. If you consider a Fourier series as an analytical solution, you have this: $$i \hbar \frac{\partial \Psi (\mathbf{r},t)}{\partial t} = -\frac{\hbar^{2}}{2m} \nabla^{2} \Psi(\mathbf{r},t)$$ Let's say you can use separation of variables and assume this for your wave function as: $$\Psi(\mathbf{r},t) = ...


4

A necessary condition for $\lfloor a/b\rfloor = \lfloor c/d\rfloor$ is $|ad-bc| < db$. This is not sufficient however, so if this test yields true (or if one of the multiplications overflows), you still must do the integer division. Depending on the probability distribution of your numbers a,b,c,d and the runtime complexity of multiplication compared to ...


4

You are almost there, just put $t$ under the square root, and it will become $ \displaystyle\int_0^{1/1095} {\frac{dt}{\sqrt{(1+644.153t)(4.17 \cdot 10^{-5} + 0.145 t)}}}, $ which eliminates the singularity at $t=0$


4

You write $S(q)$ and $T(q)$ as integrals, but it is easier to think of them as solutions of ODEs: $$ S'(q) = \sin^2\left(\frac{π\Gamma(q)}{2q}\right) $$ with initial conditions $$ S(2) = 0, $$ and similarly for $T(q)$. You can then use any of the common ODE integrators in matlab, mathematica, maple, ..., to solve and plot the solutions so that you can ...


4

You can find out if your hardware supports half-precision via: $ lscpu | grep Flags | grep f16c Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf ...


3

Obviously from the picture it seems the variation of concentration in the junction itself might be considerable that lead the authors to consider a spherical control volume at the junction. When you ignore the volume of the junction and consider it as a point and write the mass balance equation as an algebraic equation of balance of incoming and outcoming ...


3

If you use diffeqpy you can use the commands adaptive=false,dt=... to specify fixed time stepping. The following is for using the Dormand-Prince RK45 method with fixed time stepping on the Lorenz equation: from diffeqpy import de import matplotlib.pyplot as plt def f(u,p,t): x, y, z = u sigma, rho, beta = p return [sigma * (y - x), x * (rho - z)...


3

The following paper gives a good overview of well-posedness for the compressible Navier Stokes equations. https://link.springer.com/content/pdf/10.1007%2Fs10915-016-0303-9.pdf The basic building block is the energy method, where you linearize the equations, then left multiply with the primitive variables, and then integrate. This gives you a norm which you ...


3

Estimating the error: yes, this is standard material in rounding error analysis. See for instance (3.13) in Higham's Accuracy and Stability of Numerical Algorithms: if $C=AB$ and $\hat{C}$ is its computed version, then $$ |C-\hat{C}| \leq \gamma_n |A|\,|B|, $$ where the absolute values are taken componentwise and $\gamma_n = \frac{nu}{1-nu} = nu + O(u^2)$, ...


3

A search for the specific coefficients listed led me to the method ROS3PRL from J. Sieber, Konvergenzanalyse und Numerische Tests für die Prothero–Robinson–Gleichung (Master thesis), TU Darmstadt, 2014. I can't seem to find this thesis online, but the method is mentioned in the following which may be of interest. Rang, Joachim. "Improved traditional ...


3

After thinking about this some more, I can answer this one myself! I don't think the complex plane makes the log-sum-exp trick appreciably different, at least in Cartesian coordinates. In particular, if $z=u+iv$ then $e^z=e^{u+iv}=e^u (\cos v + i\sin v).$ Notice the $v$ part has magnitude 1 by construction, so overflow or underflow is principally caused ...


3

What you're looking for is an a-posteriori error estimate for your mesh study. Normally, these quantitative measures are line/area/volume- and/or time- averaged nodal or element quantities (e.g. avg. temperature of some body integrated over time). In structural analysis it is also customary to use the energy norm as a quantitative measure. A line/area or ...


Only top voted, non community-wiki answers of a minimum length are eligible