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


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


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

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


7

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


6

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


6

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


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

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

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

Depending on how large $n$ can get and how many evaluation points $x$ you wish to use, this summation problem is well-suited to the use of fast multipole methods (FMMs); for instance, see the black-box FMM, which only requires you to tell it what kernel function you want to use. In your case, it's a simple Gaussian kernel.


5

The short answer is that you need $$\phi_{-1} = \phi_0$$ $$\phi_N = \phi_{N-1}$$ to impose $\nabla\phi=0$. A quick check by making the following change if idx == -1: idx = 0 elif idx == N: idx = N-1 in the code, you have posted shows that the average $\phi$ remains constant up to 14 decimal places. To see why this is the correct boundary condition ...


5

Ill conditioning is a property of the system of equations rather than the algorithm used to solve the system of equations. Using a bad algorithm can certainly make the situation worse, but you're already in trouble when you try to solve an ill-conditioned system of equations with A coefficients or right-hand side b with even tiny errors even if you use ...


5

Q1: No. Here's a counter-example: >> A = eye(4)*1e-300 A = 1.0e-300 * 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 >> rank(A) ans = 4 >> >> rank([A, ones(4, 1)]) ans = 1 >> so, if the added column ...


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


4

The set of $x$ such that $|x_{1}| \geq \sum_{i=2}^{n} | x_{i} |$ is a non-convex set, but it's relatively easy to look for solutions where $x_{1} \geq \sum_{i=2}^{n} | x_{i} | $ and then look for solutions where $-x_{1} \geq \sum_{i=1}^{n} | x_{i} |$. The two subproblems can be solved as linear programming problems feasibility problems. If you want the ...


3

If you have a solution $\psi$ to the stationary Schroedinger equation $$ H\psi(x) = E \psi(x) $$ then the time dependent Schroedinger equation $$ i\hbar \frac{\partial}{\partial t}\Psi(x,t) = H\Psi(x,t) $$ has the solution $\Psi(x,t) = e^{-iE t/\hbar}\psi(x)$. (In other words, every pure eigenstate remains a pure eigenstate -- only the phase rotates ...


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

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


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

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

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


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