Have a look at Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods (Barrett et al.). You can find it here. Here's why I'm recommending this over other references: the ...

Compute the largest-magnitude eigenvalue $\lambda_{max}$ of $A$ (with, say, eigs('lm')). Then compute the largest magnitude (negative) eigenvalue $\hat{\lambda}_{max}$ of $M = A - \lambda_{max}I$ (...

Further to Chris' answer: Yes, weather (or the equations describing it) is extremely sensitive to the initial conditions. The fact that the weather system contains phenomena at pretty much all time ...

Given $A \in {\bf S}^n$ (a positive definite matrix) with eigenvalues $\lambda_1 \leq \lambda_2 \leq \ldots \leq \lambda_n$, then: $\displaystyle f_k(A)=\sum_{i=1}^{k} \lambda_i$ is concave. Why? $$... View answer 7 votes There exist special techniques for updating the eigen-decomposition of time-dependent covariance matrices. Given a "prior" eigenvalue decomposition (say at some initial time t^0), these recursive ... View answer 7 votes Here's another answer in the same vein as the above. You should look into the Armadillo C++ linear algebra library. Pros: The function syntax is high-level (similar to that of MATLAB). So no ... View answer Accepted answer 6 votes No, there is no specialized ARPACK routine for complex Hermitian matrices. The ARPACK authors recommend using the znaupd routine for both Hermitian and non-Hermitian problems: https://www.caam.rice.... View answer Accepted answer 6 votes While preconditioned CG with incomplete Cholesky (ICC) is reasonably straightforward to formulate mathematically, writing an efficient implementation is a non-trivial matter. Here's some of the ... View answer 6 votes I would like to hear comments from users that have some practical models (e.g. black-box hyperparameter optimization) which are still needed to be solved acceptably - whether this method works or not ... View answer Accepted answer 6 votes The measurements of this field are often spotty and missing chunks; why interpolate to get a continuous field of dubious fidelity if that can be avoided? You're perfectly right - most of the time, ... View answer Accepted answer 6 votes For domains that are logically square or cubic (like your quadrilateral), you can use the tensor product (dimension-by-dimension) approach. That is, generate the 2D Gauss-Lobatto point matrix as the ... View answer Accepted answer 5 votes This 2004 paper by Roy et. al in Int. J. Numer. Meth. Fluids 2004; 44:599–620 (DOI: 10.1002/d.660) should contain exactly what you're looking for: ftp://ftp.demec.ufpr.br/CFD/bibliografia/... View answer Accepted answer 4 votes This topic has been discussed at some length on Cross Validated (aka stats.stackexchange) and Reddit: Why is Newton's method not widely used in machine learning? (see in particular Nick Alger's ... View answer 4 votes Have a look at Pyclaw. This library has been around for quite a while and is fairly robust. It offers: Implementations of several Godunov-type methods and Riemann solvers in 1/2/3D. Adaptive mesh ... View answer Accepted answer 4 votes Some remarks: Notation: [t_1,t_2]: the time interval; \Omega: the spatial domain; \bar{u}^i(x): the known tumour profile at t_i; \left\| \cdot \right\|_\Omega: a suitably chosen norm on \... View answer 4 votes How about doing an SVD to get the null space of the 1 \times d matrix A = [\bar{v}] with \bar{v} := v^T/\|v\|? The null function in MATLAB does exactly this: d = 4; v = rand(d,1); vbar = v.'/... View answer Accepted answer 4 votes OK - so your original equation is$$Ax = b$$Say you've come up with a good preconditioner for A, call this M. Also, say you have pre-computed an LU-decomposition for this M, i.e.$$M = L_m ...

Short answer: yes (in exact arithmetic). You'll have to use the centered difference formula evaluated at $x \pm \frac{1}{2}\delta x$, like this: $$u_x = \frac{u(x + \frac{1}{2}\delta x) - u(x - \... View answer Accepted answer 4 votes Here's a short Matlab benchmark. I've used a sparse tridiagonal matrix. function benchmark_qr maxit = 25; time = zeros(maxit,2); for iter = 1:maxit mdim = 500 * iter; A = gallery('... View answer Accepted answer 4 votes If your matrices are large, why not use a library like ARPACK? The shift-and-invert mode of ARPACK will help you calculate the eigenvalues close to \sigma. There are interfaces to ARPACK for most ... View answer 3 votes An addition to the other two (excellent) answers: Understanding the source code of something like OpenBLAS can be a daunting task. As an alternative starting point, you could read through How to ... View answer Accepted answer 3 votes Your problem reads Given N points in, say, two dimensions, find an axes-parallel rectangle of minimal area that encloses at least 95% of the points. This is essentially equivalent to the ... View answer Accepted answer 3 votes I've summarized the comments thread of your original question into an answer. Here are a few things that you can try: Increase the number of your Arnoldi vectors (NCV) generated at each iteration. ... View answer Accepted answer 2 votes You want to solve$$ \arg \min_{D(x),k(x)} \mathcal{G}(u;D,k) := \frac{1}{2} \int_{\Omega} \left(u(t=T) - {\bar{u}}(T) \right)^2 d\Omega + \frac{1}{2}\int_{\Omega} D^2 d\Omega + \frac{1}{2}\int_{\...

If I understood you correctly, you are using the IDA solver from the LLNL SUNDIALS suite. If you have access to the source code that computes the residual of your DAE system, you may then consider ...

Odespy needs from you a 1st order system of ODEs, namely something of the form \begin{align} \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y},t), \tag{1} \end{align} where, for your particular problem, ...

One remark (seeing that the two answers given thus far come with quite a few comments, I've decided to post this as a separate answer instead): After you've solved the problem with your ODE ...

You can avoid the explicit computation of $A^{-1}$ in the Newton iteration for the sign of $A$. Refer to chapter 5 from this book by Higham, more specifically equation 5.22.

Starting with $b=Ax$ the forward differentiation gives $\dot b=\dot A x+A\dot x$, or, with $\dot y = \dot A x$: $$\dot x = -A^{-1} \dot y + A^{-1} \dot b$$ Let's write this in terms of the local ...
It's good to be aware of Floater-Hormann interpolants when you have to (or want to) work with equidistant points $\{x_i\}_{i=1\ldots n}$. Given the integer $d$ with $0 \le d \le n$, let $p_i$ be the ...