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

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


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

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

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

In my opinion, an example could be the calculation of the minimal solution of a three term recurrence relation (TTRR) $$y_{n+1} +a_{n}y_{n}+b_{n}y_{n-1} = 0, \quad n=1,2,3,\ldots$$. For example, the Bessel $J_{n}(x)$ function for a fixed $x$ satisfies the TTRR $$y_{n+1}-\frac{2n}{x}y_{n}+y_{n-1} = 0$$ Suppose you know $y_{0} = J_{0}(1)$ and $y_{1} = J_{1}(...


6

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


6

After you compute $Q$ and $D$, form $D'=\max(D,0)$, and compute $A'=QD'Q^\top$, the algorithms involved in multiplying those matrices do not promise that $A'$ will be exactly $QD'Q^\top$. Most commonly, they are backward stable, and promise that the actual floating-point output will be $(Q+\delta Q)(D'+\delta D')(Q+\delta Q)^\top$, for some small ...


6

What you want seems inherently impossible, and that’s not due to restrictions of Python. The only way we can arrive at a situation where we only need to apply a single quadrature is to get analytically get rid of all dependencies of $T$ in the integral. To this end, the best we can do is to apply the substitution $x=Ty$ to your integral: $$ I(T) = \int_0^∞ \...


6

As suggested in the comments, the eigendecomposition $\mathbf J = \mathbf Q \mathbf D \mathbf Q^{T}$ can be used to generate the matrix $\mathbf J^{-1/2}$, just take the eigenvectors from $\mathbf J$ (denoted $\mathbf Q$) and the inverse square root of the eigenvalues (denoted $\mathbf D^{-1/2}$) .. this is a "scalar"/"entrywise" operation because $\mathbf ...


5

What you are looking for is called "bootstrapping". It is a common problem of all multistep ODE integrators and is discussed in many books on the topic. Among your options are to use a lower-order method with smaller time step, or to use a one-step method of higher order for the first few steps (e.g., a Runge-Kutta method).


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

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

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


4

In BDF schemes for $\dot y = f$, one uses $$ f(t_n)=\dot y(t_n) $$ and tries to approximate $\dot y(t_n)\approx \sum_{j=0}^k\alpha_k y_{n-j}$ by the current value $y_n$ (that is to be computed) and the $k$ previously computed approximations. In the presented approach, in $(5)$, $y$ is approximated as a polynomial $p$ in $t$ fitted to $y_{n-j}$, so that the ...


4

I did some more investigating about the accuracy of the standard recurrence relation method vs. the newer Bunck algorithm. It seems that in fact the Bunck algorithm is generally more accurate for all $n$ and $x$. I made a python script using the mpmath arbitrary precision library to calculate the Hermite functions to lots of decimal places, and then compared ...


4

The methods are unrelated. The finite element method is a way to convert a partial differential equation into a finite dimensional problem ("discretization") so it can be solved on a computer. The Conjugate Gradient method is a method to solve linear systems such as those that arise from the finite element method.


4

The difficulty is relative to something, in this case it is relative to diffusion dominated problems. Diffusion dominated aren't "easy" either, they have their own set of problems. I'll start with some favorable qualities of diffusion problems, and then mention why they are not present for advection problems: Discretizations are often symmetric and ...


4

In Gauss-Seidel, you are using this $A=L+U$ splitting implicitly. So, you never form $L$, $U$, or $L^{-1}$ explicitly. Which is extremely good, since forming an additional matrix (not even talking about a calculation of an explicit inverse) is a huge burden. Instead, since $L$ is lower triangular, you can change the explicit inverse of $L$, by performing a ...


4

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


3

Check Algebraic Multigrid Solvers for Complex-Valued Matrices by MacLachlan and Osterlee. On the implementation side of things, PyAMG supports complex-valued matrices. I've used it before and it works well.


3

I think that there is a bit of confusion with change of variables in some of the previous answers as well as some errors. The integral of a log function is not the log of the integral. I think in general is difficult to write out the the integral of a function knowing the integral of its log. If anyone knows how to do that I would be interested. In the ...


3

Higham's algorithm seems very useful for online algorithms or algorithms in which you have limited storage capacity such as edge processing. In these situations it is probably always worth the cost. But, to address somewhat your question as posed, I implemented an SSE2 version which I think captures your question: #include <chrono> #include <...


3

Suppose you have $m$ numbers and the vector length $l$ divides $n$, i.e., $m=ln$. The you can apply Higham's scheme to each group of $n$ numbers. This vectorizes nicely, especially if the numbers are interleaved, i.e., all the first elements are contiguous in memory, followed by all the second elements, etc. When you have computed the average of each group ...


3

The common way of tracking such properties is indeed to add an additional ODE to the system, here: $$\frac{df}{dt}=a(t)+b(t)+x(t)+z(t)$$ with initial condition $$f(0)=0.$$ If this is not possible, then a smart choice of output times e.g. quadrature points in t for an integration afterwards is recommend.


3

The motivation for arranging the computation in this fashion was almost certainly the issue of wobbling precision when computing with the hexadecimal floating-point format of the IBM System/360. There, the mantissa is normalized such that $num = \pm 16^{e} \cdot f$, with $\frac{1}{16} \le f \lt 1$. This means that up to three leading bits of the mantissa can ...


3

The square root is a relatively tough mathematical operation. Thus, it is reasonable to assume that in pre-IEEE-754 dominated world, the argument range could significantly influence the accuracy/performance of SQRT. Especially, as @njuffa mentioned in the comment with ...wobbling precision of IBM System/360 hexadecimal floating point format However, I ...


3

This algorithm suffers from similar numerical stability problems as the symmetric Lanczos tridiagonalization algorithm, see here In exact arithmetic, after $k$ steps the following holds: $U_k^T U_k =I$, $V_k^T V_k =I$, and $U_k^T A V_k = B_k$, where $B_k$ is the bidiagonal matrix with diagonal elements $\alpha_i$ and superdiagonal elements $\beta_i$. ...


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