Most Popular
1500 questions
7
votes
6 answers
Python implementations of Gillespie's direct method
I'm looking for a decent implementation of Gillespie's Direct Method in Python, as if I code the algorithm myself I'm nigh positive I'll do it inefficiently.
Anyone have a favorite?
Fomite
- 2,383
- 1
- 18
- 25
7
votes
0 answers
"Geometry of ill-conditioning" for least-squares problems
It is an idea that dates back to Demmel, 1987 that the condition number of a problem is often related to the distance to the closest ill-posed problems. In Section 3 of the above paper, the author cites as examples linear systems/matrix inversion…
Federico Poloni
- 11,344
- 1
- 31
- 59
7
votes
1 answer
Polynomial Fitting from Chebyshev Coefficients
I have been reading Numerical recipes about how to create a power series approximation to a function once you have a Chebyshev approximation to the function. However it is still very unclear to me how one can go from the coefficients $c_k$ in…
tau1777
- 451
- 2
- 5
7
votes
3 answers
Does a symmetric positive definite matrix also have $\mathbf{A} = \mathbf{L}^T\mathbf{L}$ (where $\mathbf{L}$ is a lower triangular matrix)?
As we know, for a symmetric positive definite (SPD) matrix $\mathbf{A}$, there is a theorem about the Cholesky factorization $\mathbf{A}= \mathbf{L}\mathbf{L}^T$, where $\mathbf{L}$ is a lower triangular matrix.
I am a little curious whether the…
Happy
- 961
- 4
- 11
7
votes
2 answers
How to solve calculus of variations problems numerically?
For example, how to solve the well-known isoperimetric problem (i.e., to enclose the largest area with a fixed-length curve)?
We can simplify things a bit and fix the two ends of the curve at $[a,0]$, $[b,0]$, then the problems is
$$\text{maximize}…
Taozi
- 277
- 2
- 9
7
votes
1 answer
Computing square root of diag(u)-uu'?
I need an efficient way to take square root of a matrix which is a sum of diagonal matrix and rank-1 matrix.
More specifically it's the following matrix
$$A=D-uu'=\text{diag}(u)-uu'$$
Where entries of $u$ are non-negative and $\sum_i^d u_i=1$. This…
Yaroslav Bulatov
- 2,655
- 11
- 23
7
votes
3 answers
Convergence of fixed point iterations of a non-linear matrix system
I'm working on modeling two phase immiscible flow in a porous medium. When I setup the system of equations, I obtain a non-linear system of equations that can be expressed in the form:
$A(x)x=b$
where the matrix $A(x)$ is a function of the…
Paul
- 12,045
- 7
- 56
- 129
7
votes
0 answers
Implementation of Lanczos method that returns tridiagonal matrix
The Lanczos method can be used to obtain extremal eigenpairs of sparse symmetric or hermitian matrices. I know there are several implementations of the Lanczos method (as well as Arnoldi, Davidson, and others) out there, but I am specifically…
delete000
- 171
- 3
7
votes
1 answer
Numerically estimating expected value of f(x) when x is normally distributed
I need to estimate
$$
\mathbb{E}_x[f_i(x)] = \int_{\mathbb{R}^n} f_i(x) p(x) dx
$$
for many functions $f_i(x)$, where $p(x)$ is the density of a normal distribution. The evaluation of all the functions $f_i(x)$ is expensive.
I was looking into…
Rosh
- 73
- 3
7
votes
1 answer
Numerical calculation of Integral of Si(x)/x
I'm interested in evaluating
\begin{equation}
\int_0^x \frac{Si(t)}{t}\;dt
\end{equation}
Where
\begin{equation}
Si(x) = \int_0^x \frac{\sin t}{t}\;dt
\end{equation}
I've found a nice method for evaluating $Si(x)$ using Pade and Chebyshev-Pade…
Michael Anderson
- 185
- 2
- 9
7
votes
0 answers
Is there any catch on using `zgemm3m` vs regular `zgemm`?
I've just (to my embarrassment) encountered a BLAS-like extension of a matrix-matrix product subroutine gemm in Intel MKL: gemm3m. This subroutine (particular versions: cgemm3m and zgemm3m) allows performing matrix-matrix multiplication for…
Anton Menshov
- 8,672
- 7
- 38
- 94
7
votes
1 answer
Appropriate space for weak solutions to an elliptical pde with mixed inhomogeneous boundary conditions
I'm working with the following mixed inhomogeneous boundary value problem:
$\nabla(\kappa\nabla u)=f$ in $\Omega$
with $\partial\Omega = \Omega_1 \bigcup\Omega_2$ such that
$u=g$ on $\partial\Omega_1$
$\kappa\nabla u\cdot n = h$ on $\Omega_2$ …
Paul
- 12,045
- 7
- 56
- 129
7
votes
3 answers
Spectral Element vs Finite Element
I am trying to understand the difference between SEM and FEM. If I go by this paper, spectral element methods are a subset of FEM methods and the only difference lies in the choice of basis functions. If this is the case, are there any advantages in…
efso
- 73
- 1
- 3
7
votes
3 answers
Choosing subset of vectors to approximate a subspace
Suppose I have a high-dimensional vector space $X$, a subspace $V \subset X$, and a collection of $n$ vectors $\{x_i\}_{i=1}^n \subset X$.
My question is: How can I choose a small collection $k < n$ of the vectors $x_i$ so that the span of this…
Nick Alger
- 3,143
- 15
- 25
7
votes
2 answers
transitive floating point comparison with (absolute) tolerance
I want to compare two floating point numbers for equality relative to a known absolute tolerance. However, this is inside an algorithm I wrote quite some time ago, and I believe the logic of that algorithm would get corrupted if the equality…
Thomas Klimpel
- 2,141
- 15
- 35