19
If what you want is to solve
$\min \| Ax - b \|_{2}^{2} + \lambda^{2} \| x \|_{2}^{2}$
subject to
$x \geq 0$,
then this is easily implemented. Construct a matrix
$C=\left[
\begin{array}{c}
A \\
\lambda I
\end{array}
\right]$
and a vector
$d=\left[
\begin{array}{c}
b \\
0
\end{array}
\right]$.
Then use your nonnegative least squares solver on
$...
18
There are two issues that you are likely to be encountering.
Ill-conditioning
First, the problem is ill-conditioned, but if you only provide a residual, Newton-Krylov is throwing away half your significant digits by finite differencing the residual to get the action of the Jacobian:
$$ J[x] y \approx \frac{F(x+\epsilon y) - F(x)}{\epsilon}$$
If you ...
11
The best approach is to use an ODE solver that is guaranteed to conserve the norm of the initial condition, i.e., for which $\|y_n\| = \|y_0\|$ for all $n\in\mathbb{N}$. Such solvers exist, and are called geometric integrators, since they preserve geometric properties of the exact solution (in this case, that energy is conserved, i.e., $\frac{d}{dt}\|y(t)\| =...
9
You can generate a noise sequence with whatever noise spectrum you want (including $1/f$, also known as pink noise) by generating the noise coefficients in spectral space. The magnitudes of the coefficients should be chosen to give the desired spectrum and the phases should be chosen randomly. You then simply perform an inverse Fourier transform to give the ...
9
Any general-purpose ODE solver should be able to handle this linear coupled system of ODE very easily, for example:
scipy.integrate.ode
CVODE from the Sundials solver suite; it appears to have Python bindings here, and perhaps there are others.
This kind of thing is typically discussed in any textbook on numerical methods.
In general, computing the matrix ...
8
The problem is almost definitely with how QUADPACK (which is the backend used by scipy.integrate.quad) handles numerically small integrands. Essentially the integrand is so small (at $x=0$ it is $6.58\times 10^{-12}$), that it is comparable to the absolute error tolerance. I'm not sure quite why it says it's "probably divergent", but that doesn't matter. I ...
8
The feature that you demand is called event location in Matlab ODE solvers pack, or rootfinding in SUNDIALS solvers suite terminology. Essentially this feature allows to stop integration exactly at the point where some vector function of free and dependent variables has a root. Namely, for system
$$
\frac{d\boldsymbol y}{dt} = f(t, \boldsymbol y),
\quad \...
8
When you use $r=5$, the initial condition is
$$ x(-10) \approx e^{-50} \approx 1.9\times 10^{-22}. $$
This is much smaller than the machine epsilon, $2\times 10^{-16}$, and it is very likely that LSODA just concludes that the solution is identically zero.
To check this idea, I replaced $t_0$ with $-5$ instead of $-10$, so that $x(-5)\approx 1.4\times 10^{-...
8
The matrix B (M in the documentation) needs to positive definite according to the documentation: "If sigma is None, M is positive definite", this is in addition to the first requirement "M must represent a real, symmetric matrix if A is real" which your B follows. The eigenvalues of your current matrix B are -1, 1 and 6. So matrix B is ...
7
The reason the resulting geodesic curve was deviating was because the calculated Christoffel symbol of second kind was incorrect. Using the correct Christoffel symbol :
C = Matrix([[(0, -tan(v)), (0,0)],[(sin(v)*cos(v),0),(0, 0)]])
results in the proper output (as displayed in the reference) :
Now, I suppose I have to figure out why the calculated ...
7
Took me quite a while to figure this out and as usual it becomes obvious after you find the culprit.
After checking the problematic cases reported in David S. Watkins. A case where balancing is harmful. Electron. Trans. Numer. Anal, 23:1–4, 2006
and also the discussion here (both being cited in arXiv:1401.5766v1), it turns out that matlab uses the ...
7
One of your problems is the system of units that you are using. Just changing the units improves the results
import numpy as np
import scipy.integrate as integrate
eigenvalue = [0.9, 1.3]
fermi = 1.0
T = 300
kB = 8.6173303e-5
def fermi_integral(E, fermi, T):
return 1 / (1 + np.exp((E - fermi) / (kB * T)))
for i in range(len(eigenvalue)):
result = ...
7
The code proposed by the OP can indeed made be more efficient, mainly by noting the fact that to form the sequence $A^i B$, with $i=0\,\dots,N$ you do not have to compute $A^i$ at each step, but you can exploit the fact that $A^i B = A\,(A^{i-1}B)$, reusing the result of the previous step.
My proposed implementation is
import numpy as np
N = 3
M = 1
A = ...
7
I believe you can accomplish what you want efficiently using the recursive LU algorithm. In brief, recursive LU on a $M \times N$ matrix $A$ proceeds by partitioning the matrix into 4 blocks:
\begin{align}
\pmatrix{A_{11} & A_{12} \\ A_{21} & A_{22}} &= \pmatrix{L_{11} & 0 \\ L_{21} & L_{22}} \pmatrix{U_{11} & U_{12} \\ 0 & U_{22}}...
7
I think you are overestimating the overhead of computing L. There are zero extra operations needed; the only additional cost is writing to RAM some numbers that you have already computed anyway. The algorithms commonly used (in Lapack, for instance) to compute U also compute L along the way, and you'd save 0 flops by omitting it.
For instance, if you think ...
7
When one says an algorithm is of order $O(n)$, that may mean that the complexity is given by: $c + b*n$. With every new element you add you increase in runtime (effectively). What mathematically minded people often forget is that these statements do not include how large the constants are. That of course carries over to $O(n²)$ and such. I can not answer ...
6
Looking at the information of nympy.linalg.solve for dense matrices, it seems that they are calling LAPACK subroutine gesv, which perform the LU factorization of your matrix (without checking if the matrix is already lower triangular) and then solves the system. So the answer is NO.
Otherwise, it makes sense. If you do not have an easy (cheap) way to verify ...
6
No. The numpy.linalg.solve method uses LAPACK's DGESV, which is a general linear equation solver driver. If you know that your matrix is triangular, you should use a driver specialized for that matrix structure.
scipy.linalg.solve does something similar.
MATLAB detects triangularity in a solve if you use the backslash operator; see this page for pseudocode....
6
FILTLAN is a C++ library for computing interior eigenvalues of sparse symmetric matrices. The fact that there is a whole package devoted to just this should tell you that it's a pretty hard problem. Finding the largest or smallest few eigenvalues of a symmetric matrix can be done by shifting/inverting and using the Lanczos algorithm, but the middle of the ...
6
From Ablowitz and Zeppetella we know that the analytical solution reads:
\begin{equation}
u(x,t)=\frac{1}{1+e^{{(-\frac{5}{6} t+ \frac{\sqrt{6}x}{6}})^2}}
\end{equation}
Usually, analytical solutions require the imposition of boundary conditions. Are there any used in the derivation of this expression? If so, you must use the same boundary conditions ...
6
First of all, your function $x\sin(\frac{1}{x})$ is singular in $x=0$. You might want to add an if clause like this:
def f(x):
if abs(x) < 1e-10:
res = x
else:
res = x*sin(1/x)
but this does hurt speed (masked arrays would be better).
The reason why your code doesn't work is because
scipy.quad only accepts a single value as ...
6
Here's a few options that are relatively easy to work with:
One node - multiprocessing is the most straightforward thing to do. multiprocessing.map works well for an embarrasingly parallel problem.
GPU - pyCUDA allows you to do some GPU level programming with just python. I haven't toyed with it that much.
MPI - mpi4py- Exposes MPI interface at the python ...
6
The reverse Cuthill-McKee algorithm produces a reordering that applies to both the rows and columns. This is because it works by considering matrices as graphs of (undirected) connected nodes. According to the function's documentation in SciPy, the output array is the permuted row/column indices, so you can simply do the following
perm = ...
6
The function $q(e)$ satisfies a first order linear ODE
$$ \frac{\mathrm{d}q}{\mathrm{d}e} = \frac{111 e^4+876 e^2+288}{(e^2-1)
(121 e^2+304)} q(e), $$
which can be solved very easily by using an integrating factor:
$$ q(e) =
C_1 \exp\left(
\tfrac{111}{121} e-3\operatorname{arctanh}e+\tfrac{10440}{1331 \sqrt{19}} \operatorname{arctan}\left(\tfrac{11}{4
\sqrt{...
6
If the only non-zero entries of $A_{ij}$ have $j$ in $\{i - 1, i, i + 1\}$, then $A$ is a banded matrix with bandwidth 1.
More generally, you can talk about matrices of bandwidth $k$ where $k$ is any integer.
For example, if you were using a higher-order finite difference discretization that used more points to calculate a derivative, you'd get higher-...
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
Do you have more suggestions to improve the precision and accuracy of the method ? What are the cause of such errors from the mathematical point of view ? For instance if we have different time scale in the dynamics ...
You need to use a reversible ODE solver method if you want to do this. I actually recently showed in a blog post that there are many cases ...
6
The interpolated polynomial does not have roots. Considering that the behavior outside the interpolation region holds is termed extrapolation.
You can explicitly use the polynomial, given by (as I explained in this post)
$$f(x) \approx N_1(x) u_1 + N_2(x) u_2 + |J|(N_3(x) u'_1 + N_4(x) u'_2)\quad \forall x\in [a, b]$$
with $|J| = (b - a)/2$ the Jacobian ...
6
The Savitzky-Golay filter uses a constant delta (the spacing of the samples,) and the default value of the delta in the filter implementation is 1, according to https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.signal.savgol_filter.html. Your data set has irregular deltas, not 1, so the result from the Savgol filter is incorrect.
5
In general, numerical methods do not conserve quantities that may be conserved in a continuous system. The discretized system is typically non-conservative. A smaller time step (equivalently a smaller error tolerance) should help, but if you have a quantity that must be precisely conserved, you will have to use a method that is designed to conserve the ...
Only top voted, non community-wiki answers of a minimum length are eligible
