33 votes

C, Julia, Python, Maxima, Mathematica, ChatGPT and numerical errors

This function, as any similar hill, tent or hat functions, implements the expansion-and-folding scheme. This means that in the most benign interpretation you lose in each iteration one bit from the ...
Lutz Lehmann's user avatar
  • 5,984
32 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-...
Federico Poloni's user avatar
30 votes

C, Julia, Python, Maxima, Mathematica, ChatGPT and numerical errors

I thought it would be interesting to see how the number of bits of precision affects the error. I wrote the following using arb, a library for interval arithmetic (and in particular a library that ...
davidlowryduda's user avatar
29 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 ...
Chris Rackauckas's user avatar
29 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 ...
Jeff Hammond's user avatar
  • 2,116
24 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 ...
Brian Borchers's user avatar
24 votes

More stable algorithm to calculate `sqrt(a^2 + b^2) - abs(a)` in MatLab

A manipulation that may help is the following. Assume for simplicity $a>0$. We have the identity $$b^2 = (a^2+b^2)-a^2 = (\sqrt{a^2+b^2}-a)(\sqrt{a^2+b^2}+a),$$ hence $$ \sqrt{a^2+b^2}-a = \frac{b^...
Federico Poloni's user avatar
18 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 ...
Chris Rackauckas's user avatar
18 votes
Accepted

More stable algorithm to calculate `sqrt(a^2 + b^2) - abs(a)` in MatLab

You can break down the domain of your function into three distinct cases: $|a|\gg |b|$: In this case, $\sqrt{a^2+b^2} \approx |a|$ and a naive application of the formula will likely result in poor ...
Wolfgang Bangerth's user avatar
17 votes
Accepted

Why aren't Krylov subspace methods popular in the Machine Learning community compared to Gradient Descent?

On a basic level, I don't buy the argument that you have to "solve a linear system for many machine learning algorithms". Much more, you usually have to optimize a non-linear equation which ...
davidhigh's user avatar
  • 3,127
16 votes

C, Julia, Python, Maxima, Mathematica, ChatGPT and numerical errors

Other answers have explained the results from Python, Julia, Maxima, C, and Mathematica. I'll explain the result you got from ChatGPT. Generally, the way that ChatGPT works is that you ask a question, ...
Tanner Swett's user avatar
14 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.$$
coolguy1000000's user avatar
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; ...
J. M.'s user avatar
  • 3,155
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 ...
Lutz Lehmann's user avatar
  • 5,984
13 votes
Accepted

How can I numerically integrate the Kepler problem?

I'm going to assume for the moment that your code is correctly implemented and that this problem isn't a bug. Believe it or not, gradual increase of the energy is the expected behavior of most simple ...
Daniel Shapero's user avatar
13 votes

C, Julia, Python, Maxima, Mathematica, ChatGPT and numerical errors

Perhaps the reason you are stunned is that you don't understand that what you are asking to calculate is impossible to ever calculate exactly. So if you want the answer to some precision, you need to ...
Secto Kia's user avatar
  • 239
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 ...
whpowell96's user avatar
  • 2,259
12 votes
Accepted

Real-world applications of eigendecomposition?

You can compute analytic functions of matrices using the eigendecomposition (or more generally by using the Jordan normal form in case the matrix is defective), you cannot do so with the singular ...
lightxbulb's user avatar
  • 1,994
12 votes

C, Julia, Python, Maxima, Mathematica, ChatGPT and numerical errors

Unfortunately this comes down to the math library. IEEE754 does not mandate exactly how the cosine should be computed (technically one doesn't even need to include a function to compute the cosine to ...
Federico Poloni's user avatar
12 votes

C, Julia, Python, Maxima, Mathematica, ChatGPT and numerical errors

We could use exact arithmetic to compute the correct value. A range of exact numbers $e$ correspond to the same $a=5.0$ representation in float32 IEEE754, plot over this range: You can see that after ...
Yaroslav Bulatov's user avatar
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 ...
Pseudonym's user avatar
  • 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 ...
njuffa's user avatar
  • 1,865
10 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 ...
Carl Christian's user avatar
10 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 ...
spektr's user avatar
  • 4,228
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 ...
Wolfgang Bangerth's user avatar
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 ...
Steven Roberts's user avatar
10 votes

Accurate computation of logbinomial(x,y)

The proposed computation of $\ln{x \choose y}$ via the logarithm of the beta function does not give rise to catastrophic cancellation other than for cases where $y=0$, or $y =x$, which are easily ...
njuffa's user avatar
  • 1,865
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)-...
Kirill's user avatar
  • 11.4k
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 ...
Reid.Atcheson's user avatar
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}{\...
Kirill's user avatar
  • 11.4k

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