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 ...
14
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
$...
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
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
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 ...
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^{-...
7
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 ...
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 = ...
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
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 \...
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 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 ...
5
1) If you're just looking to solve the PDEs without any other optimization, then my answer would be "none of them". Algorithms that discretize partial differential equations and then solve them as algebraic equations are massively parallelizable. It is possible to solve a PDE over a billion point mesh. Algorithms for nonlinear programming have made great ...
5
Since your problem is small, you're probably best off trying fsolve or root. Both of these are interfaces to MINPACK and call HYBRD or HYBRJ. Since calculating a Jacobian matrix for your system shouldn't be hard (either do it by hand, or use your favorite computer algebra system, like SymPy, Sage, Maple, or Mathematica), you should supply a Jacobian matrix. ...
5
I'm not sure if this is the answer, but consider the case as $\alpha \rightarrow \infty$. In this case you're trying to solve
$\iint_{\Omega} \cosh(P) dxdy = 0$
And since cosh looks like this:
You can see that integrating it over any domain will never give you zero. If the diffusion of $P$ ($\nabla^2 P$) is not large enough to counteract this limit then ...
5
I second the idea of using Eigen, which is pretty efficient, but also very simple to include.
If you need a lot more performance, you could try to use PETSc or Trilinos. They are very powerful libraries to store and solve sparse systems, they allow for a large number of iterative or direct solvers and are compatible with MPI for added performance. However, ...
5
I would say that the implementation + verification + unit testing would take you more than just 3 weeks. Although, if you are planning to invest that time, you might add that capabilities to scipy.sparse or scikits-sparse.
Regarding symmetric sparse matrices, you can check Pysparse. It has Sparse Skyline format (SSS) that is used for symmetric matrices. You ...
Only top voted, non community-wiki answers of a minimum length are eligible
