问:
How to find eigvals/eigvecs of huge symmetric tridiagonal matrix using multiple processes or threads?
In python, there is a function in scipy (scipy.linalg.eigh_tridiagonal) which is extremely efficient for an exact diagonalization of a symmetric tridiagonal matrix (I think it uses MRRR algorithm, not sure). But only a single core is assigned to solve the exact diagonalization. Is there a way of using remaining cores to speed up finding eigenvalues and eigenvectors?


问:
How to implement a generic monte carlo algorithm for n-dimensional integration?
A very visual picture for Monte Carlo integration is the approximation of $\pi$, by sampling in a square which contains a quarter of the unit circle.

We can extend this picture to 3 dimensions, by sampling a cube which contains 1/8 of a sphere.

But both cases are quite simple. I want to write a function that generalises Monte Carlo sampling for an arbitrary (yet integrable) function. Given a function $f(x_1, x_2, ..., x_n)$ how would I implement a generic algorithm? My thought process is this:

Treat $f$ as the "height" of my sampling hyperspace. If I randomly sample a point $P = (x_1, x_2, ..., x_n)$ and a random height $y$, then if $f(P) <= y$ my MC algorithm returns true, otherwise false. This picture holds in 2D and 3D dimensions, yet I'm not sure if it generalises to n-Dimensions.

Furthermore, the boundaries for a random $P$ are clearly given by the integration boundaries, which are provided. But how do I find a good boundary for $y$? I could try to find the maximum in the given interval, and then sample the hypercube.

This seems like a very visual (and thus likely ineffective) approach.

Is this a correct way to do it?

Is there a better way?

ContentLicense: CC BY-SA 4.0


问:
Automatic differentiation of integral
I would like to compute the integral

$$
\frac{\partial}{\partial{x}} \int_{-1}^1 \int_{-1}^1 F(x,y,\xi, \eta) \; d\xi \; d\eta
$$

or moving the derivative inside the integrals

$$
\int_{-1}^1 \int_{-1}^1 \frac{\partial}{\partial{x}} F(x,y,\xi, \eta) \; d\xi \; d\eta
$$

where $F(x,y,\xi, \eta)$ is a complicated rational function, the derivative of which is certain to result in expression swell.

I have come across only one example (a blog post containing an example) of automatic differentiation (AD) of an integral and I would like to see more, specially involving multiple integrals.

Title: Automatic differentiation of integral


问:
Iterative methods for underestimate of smallest eigenvalue for large sparse matrices
I recently read the paper "EUCLIDEAN-NORM ERROR BOUNDS FOR SYMMLQ AND CG" by Estrin et al. and there they use an underestimate (i.e. something in $(0,\lambda_{min}]$) of the smallest eigenvalue of the matrix in order to derive a practical upper bound for the L2 error of the iterations.

In my case I am dealing with very large and sparse matrices, and I would like to find such an underestimate for the smallest eigenvalue in order to be able to apply the above method. Is there some algorithm that is considered state of the art in that regard? What I am aware of is the inverse power iteration however it requires solving a sequence of problems with the system matrix which is not very fast, and I also do not know what the precision guarantees are or whether I am guaranteed to get an underestimate of the eigenvalue. I also looked into Saad's book, but there weren't very many details on this specific problem.


问:
Examples of kinetic modeling with optimization techniques in Python
We are looking for studies or tools that implemented kinetic modeling with parameter estimation differential evolution or similar optimization techniques in Python.

We are trying to understand what already exists in this area.
Thanks in advance!


问:
Numerical instability in the inverse Laplace transform
I have a problem with Laplace inversion and my function is not numerically stable for the Laplace inverse, but I do not understand the cause of this problem.
Here is my code and graph of this problem. Does anyone have any idea? I have also tried all 3 Dehoog Stehfest and Talbot, they all have a different type of instabillity

```<code>import mpmath as mp
import numpy as np
import matplotlib.pyplot as plt
def Laplace_Max(s):
    #mp.dps = 100
    A= s**(5/4)/ (
        0.00480931 +
        0.0244077*s**(1/4) +
        0.0129056*mp.exp(-35*s)*s**(1/4) +
        0.00329167*mp.exp(0.707997*s)*s**(1/4) * mp.gammainc(0.0, 35.708*s,mp.inf, regularized=True) -
        0.00530593*mp.gammainc(1.25, 35*s,mp.inf,   regularized=True)
    )
    return A
t_R_inv = np.linspace(1,1200,12000)
Max_R=np.zeros(len(t_R_inv))
for i in range(len(t_R_inv)):
    Max_R[i]=mp.invertlaplace(Laplace_Max, t_R_inv[i], method = 'stehfest', dps = 100, degree = 1000)
plt.plot(Max_R)
</code>```

Update
The inverse laplace of this function!

$\frac{1}{s}^{\frac{5}{4}}\cdot (0.00480931 + 0.02440766 \cdot s^{\frac{1}{4}} +
   0.0129056 \cdot e^{-35 \cdot s} \cdot s^{\frac{1}{4}}+
   0.003291670 \cdot e^{0.7079967\cdot s} \cdot s^{\frac{1}{4}}
    \cdot \Gamma[0, 35.7079967 \cdot s] -
   0.00530593 \cdot\Gamma[1.25, 35\cdot s])$


问:
Open boundary condition for 1d wave equation with variable wave speed using finite differences
I have implemented a finite difference solver for the 1d wave equation with variable wave speed:

$$ u_{tt} = c(x)u_{xx}, \hspace{10mm}c(x) = \dfrac{6 -x^2}{2} \hspace{5mm} $$
on $-2 \leq x \leq 2, t > 0$ with initial conditions: $$u(x,0) = \cos(\dfrac{\pi x}{0.2}), \hspace{5mm} u_t(x,0) = 0
$$ for $-0.1 \leq x \leq 0.1$.

I use second order finite differences in order to solve for each grid point in time:
$$\frac{u_j^{n+1} - 2u_j^n + u_j^{n-1}}{k^2} = c(x_j)\frac{u_{j+1}^n - 2u_j^n + u_{j-1}^n}{h^2}$$

Rearranging for $u_j^{n+1}$:
$$u_{j}^{n+1} = 2u_j^n - u_j^{n-1} + c(x_j)q^2(u_{j+1}^n - 2u_j^n + u_{j-1}^n)$$

Where $q = k/h$. The boundary condition to be enforced is the open condition, so  $qu_{x} = u_{t}$ at $x = -2$, and $qu_{x} = -u_{t}$ at $x = 2$. Using the first order finite differences, we can obtain an equation for the ghost points at the boundaries:

$$ q(u_{1}^{n} - u_{-1}^n) = u_0^{t+1} - u_0^{t-1}
$$
at $x = -2$.

Next, we can substitute this into the second order equation ($c(x_j)$ disappears since it is 1 at $\pm 2$) in order to end up with our iterative formula for the boundary at $x=-2$:

$$u_{0}^{n+1} = \dfrac{2u_0^n + (q - 1) u_0^{n-1} + 2r^2(u_{1}^n - u_0^n)}{q + 1}$$

The same can be done for the boundary at $x = 2$.

My python implementation is shown below.

```<code>def one_dimensional_wave_solver(x_steps, t_steps):
    # discrete points in space
    x = np.linspace(-2, 2, x_steps + 1)
    # distance of step
    dx = x[1] - x[0]
    # discrete points in time
    t = np.linspace(0, 10, t_steps + 1)
    # time steps
    dt = t[1] - t[0]

    # define cur_u and previous u
    cur_u = np.zeros(x_steps + 1)
    prev_u = np.zeros((x_steps + 1, t_steps + 2))

    # Define mesh constant
    c = 1
    mesh = (c * dt) / dx

    # Initial conditions
    for i in range(0, x_steps + 1):
        if (-0.2 <= x[i] and x[i] <= 0.2):
            prev_u[i, 0] = np.cos((np.pi * x[i]) / (0.2))
        else:
            prev_u[i, 0] = 0

    # Iterate through time steps
    for t in range(0, t_steps + 1):
        # For the first time step
        if t == 0:
            for i in range(1, x_steps):
                cur_u[i] = prev_u[i, t] + (cx(x[i]) * (mesh ** 2) / 2) * (prev_u[i - 1, t] - 2 * prev_u[i, t] + prev_u[i + 1, t])
            # Boundary conditions
            cur_u[0] = prev_u[0, t] +  mesh * (prev_u[1, t] - prev_u[0, t])
            cur_u[-1] = prev_u[-1, t] - mesh * (prev_u[-1, t] - prev_u[-2, t])

        else:
            # equation for computing u
            for i in range(1, x_steps):
                cur_u[i] = 2 * prev_u[i, t] - prev_u[i, t - 1] + cx(x[i]) * (mesh ** 2) * (prev_u[i - 1, t] - 2 * prev_u[i, t] + prev_u[i + 1, t])
            # Boundary conditions
            cur_u[0] = (2 * prev_u[0, t] + (mesh - 1) * prev_u[0, t - 1] + (2 * mesh ** 2) * (prev_u[1, t] - prev_u[0, t])) / (mesh + 1)
            cur_u[-1] = (2 * prev_u[-1, t] + (mesh - 1) * prev_u[-1, t - 1] + 2 * (mesh ** 2) * (prev_u[-2, t] - prev_u[-1, t])) / (mesh + 1)

        # Swap variables
        prev_u[:, t + 1] = cur_u.copy()

    # Return final u value
    return prev_u[:, -1], prev_u
</code>```

Using this allows me to get an almost fully transparent boundary. But there is still some slight reflection. Is there an issue with the first order equations for the boundaries?


问:
How "kinematic" rigid bodies are implemented in physics engines
In most physics engines there's this concept of "static" bodies, which act as rigid bodies with infinite mass. Then there are "kinematic" bodies that act as static bodies, but their motion is simulated by animation data or by user specifying position or velocity directly. And finally there are the "dynamic" bodies which are physics-driven rigid bodies with non-infinite mass.

I'm curious how these kinematic bodies fit into Lagrangian mechanics, and how the equality and inequality constraints are resolved for these. Say we want to simulate a hinge joint between rigid bodies. For two dynamic rigid bodies we can define a position constraint and apply constraint forces to satisfy the constraint. But how would we go about simulating this joint for a kinematic and a dynamic rigid body?

I can't find any information on this subject, and any insight on these types of kinematic bodies would be helpful.


问:
structured quad mesh with transition zone
I wonder which software can be used to generate the mesh like below. Anybody has an idea?


问:
How to estimate stability and stiffness of a system of coupled ODEs?
I'm running into issues with Python/Julia ODE solvers requiring prohibitively small timesteps to evolve a system of 4 coupled ODEs (the order of magnitude of the state variables and time unit span ~40-50 orders of magnitude). Interestingly if I use a simple explicit Euler approach then I can get past the stiffness at early times (where the general adaptive solvers get stuck) and roughly match theoretical expectations at later times.

Apart from the technical coding details, I have noticed that it is common to talk about the mathematical stability and stiffness of ODE systems. Is there a standard way to calculate the expected stability and stiffness given a set of coupled ODEs, initial conditions, integration timescales, input parameters, and maybe a solver algorithm?

For what it's worth, my system of 4 coupled ODEs looks like

$$\frac{dM_1}{dt} = g(t) - c(t,M_1,E)M_1(t) + w(t)M_3(t) - e(M_1,E)h(M_1,E)E(t)$$

$$\frac{dM_2}{dt} = c(t,M_1,E)M_1(t) - s(t)M_2(t) - w(t)M_3(t)$$

$$\frac{dM_3}{dt} = s(t)M_2(t)$$

$$\frac{dE}{dt} = g(t)v(t) - c(t,M_1,E)E(t) + b(t)w(t)M_3(t) - h(M_1,E)E(t)$$

where the various lettered functions are input parameters that depend on time and/or state variables. I know this isn't enough for doing stability/stiffness analysis but it's a start. I'd appreciate any pedagogical references, even for much simpler illustrative systems -- I vaguely remember this involves linear algebra and eigendecomposition but it's been so long since undergrad.


问:
CFD in Python: Mesh Generation
I am trying to set up a 3D CFD scheme for thermal and flow modelling in Python using the finite volume method. The first concern is to build the geometry and an accompanying mesh that is efficient for the purpose (resembling something like the 3rd picture below).

For Python, there are many libraries for geometry and mesh creation (the last two are for optimization and formatting respectively):

MeshPy

pygmsh

PyDistMesh

mshr

pygalmsh

dmsh

optimesh

meshio

Somewhat overwhelmed, I tried the MeshPy library since its description mentions finite element simulations specifically. However, there does not appear to be a lot of documentation with the libraries. I was unable to construct a mesh that could be fed to a solver (such as fipy) - also inhibiting the required properties or types for each material; in order to specify for example the initial temperature, thermal conductivity, specific heat and density.

I would like to inquire for a working example of a grid generation that could be fed to a CFD solver, from which ever library is preferred, preferably such that different grid types can be created as seen on the image above, and such that the mesh generator is compatible with STEP format. As geometry a simple box, sphere or cylinder is abundant.

Title: CFD in Python: Mesh Generation


问:
Finite volume method for 1D heat equation in 1D
I wish to solve the following using the finite volume method:
$$\frac{\partial u}{\partial t}=\frac{D}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u}{\partial r}\right)+Q(t,r)$$
with the following boundary conditions:
$$\frac{\partial u}{\partial r}\Bigg|_{r=0}=0,\quad -D\frac{\partial u}{\partial r}\Bigg|_{r=R}=hT(t,R)$$
Integrating the PDE (multiplied with $r$) from $r_{j-\frac{1}{2}}$ to $r_{j+\frac{1}{2}}$ to get:
$$\frac{d}{dt}\int_{r_{j}-\frac{1}{2}}^{r_{j}+\frac{1}{2}}rudr=D\left[ r\frac{\partial u}{\partial r}\right]_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}+\int_{r_{j}-\frac{1}{2}}^{r_{j}+\frac{1}{2}}rQ(t,r)dr$$
I denote for general $f(t, r) $:
$$\tilde{f}_{j}(t)=\frac{1}{\delta r}\int_{r_{j}-\frac{1}{2}}^{r_{j}+\frac{1}{2}}rf(t,r)dr$$
which makes my equation:
$$\frac{d\tilde{u}_{j}}{dt}=\frac{D}{\delta r}\left[ r\frac{\partial u}{\partial r}\right]_{r_{j-\frac{1}{2}}}^{r_{j+\frac{1}{2}}}+\tilde{Q}_{j}$$
Now is this a valid approximation?
$$r\frac{\partial u}{\partial r}\Bigg|_{r=r_{j+\frac{1}{2}}}=r_{j+\frac{1}{2}}\cdot\frac{\tilde{u}_{j+1}-\tilde{u}_{j}}{\delta r}$$
and likewise for the other flux. I'm also having issues with implementing the inner boundary condition at $r=0$ which is confusing me a lot.


问:
Convergence of Evolutionary Algorithms
When it comes to Evolutionary Algorithms (e.g. Genetic Algorithm), I have often heard people make the following broad statement: "Evolutionary Algorithms Do Not Converge."

I was curious about know how true this statement is, or does it come with any caveats. For instance, I was interested in knowing if perhaps there might be some conditions in which Evolutionary Algorithms have been known to converge. Reading online, I found the two following links to recent research papers on this topic:

Average Convergence Rate of Evolutionary Algorithms in Continuous Optimization, Chen et al., 2021

It is proven that the ACR is linear for evolutionary programming using positive landscape-adaptive mutation, but sublinear for that using landscape-invariant or zero landscape-adaptive mutation.

Sufficient Conditions for Global Convergence of Differential Evolution Algorithm, Hu et al., 2013

Sufficient condition and a corollary for the convergence of DE (Differential Evolution - a type of Evolutionary Algorithm) to the global optima are derived by using the infinite product.

All this brings me to my question: There clearly seems to be some research and evidence suggesting that Evolutionary Algorithms can "sometimes" converge for some conditions - yet, we are constantly (informally) reminded that the "Evolutionary Algorithms Do Not Converge"; or in the instances where Evolutionary Algorithms do perform well, it's quickly said that "its very likely a Non-Evolutionary Algorithm could have performed equally well".

In the end, I am not sure what to believe, my knowledge of mathematics is insufficient to fully understand the literature and thoroughly be able to comment on this - Do we have any mathematical knowledge suggesting that in certain cases (e.g. certain types of Loss Functions) Evolutionary Algorithms have a "strong ability" to converge? Or do we only have the famous Schema Theorem - all we can be sure of, is that Evolutionary Algorithms are probabilistically guaranteed to provide better results relative to previous iterations, but no guarantee on convergence.

Can someone please comment on this?

Title: Convergence of Evolutionary Algorithms


问:
Numerically solving a 6th order non-linear differential equation in Matlab
I've posted yesterday a question about solving a non linear equation : it was not clear so I am reformulating my question.

I am trying to solve a high-order non linear differential equation presented right below with $F \in [-\eta_{0}, \eta_{0}]$

$\left(5F - 4\eta\frac{dF}{d\eta}\right) = 9 \frac{d}{d\eta}\left(F^{3} \frac{d^{5}F}{d\eta^{5}}\right)$

As boundary conditions, I have the following ones :

$\left\{
    \begin{array}{ll}
        F' = F'''=0 & \mbox{in } \space \eta = 0 \\
        F = F' = 0  & \mbox{in} \space \eta =\pm \space\eta_{0}
    \end{array}
\right.$

I am trying to apply a shooting method algorithm in order to solve this equation on $[0, \eta_{0}]$, I have converted it in a system of 1-st order differential equations as following :

$f_{i}' = \left[\begin{array}{@{}c@{}}
    y_{2}=\frac{dF}{d\eta} \\
    y_{3} = \frac{d^{2}F}{d\eta^{2}}\\
    y_{4} = \frac{d^{3}F}{d\eta^{3}} \\ y_{5} = \frac{d^{4}F}{d\eta^{4}} \\ y_{6} = \frac{d^{5}F}{d\eta^{5}}\\ y_{7} = \frac{d^{6}F}{d\eta^{6}}
    \end{array} \right]$ with : $y_{7} = \left[\left(\frac{5}{9y_{1}^{2}}\right) - \left(\frac{4xy_{2}}{9y_{1}^{3}}\right)-\left(\frac{3y_{2}y_{6}}{y_{1}}\right)\right]$

According to the boundary conditions we thus have : $y_{2}(0) = y_{4}(0) = 0$ and we have to guess the remaining ones : $y_{1}(0), y_{3}(0),y_{5}(0) \space \mbox{and} \space y_{6}(0)$

I have attached to my post, the algorithms that I wrote which are inspired by biblographical researchs : the main code is called "shoot4nl.m" and I have tried to guess some initial values in order to run the algorithm. Also, I have choses a value of $\eta_{0}$ very "high" in order to match the conditions...

Here are the algorithms that I am using : These are MATLAB scripts so I am sharing it below. The main function is the one called "shoot4nl".

```<code>function shoot4nl

% Shooting method for nonlinear 4th-order boundary
% value problem

global XSTART XSTOP H   % Make these params. global.
XSTART = 0; XSTOP = 50; % Range of integration.
H = 0.25;                % Step size.
freq = 1;               % Frequency of printout.
u = [1 -1 1 0];             % Trial values of u(1)
                        % and u(2), u(3) and u(4).

x = XSTART;
u = newtonRaphson2(@residual,u);
[xSol,ySol] = runKut5(@dEqs,x,inCond(u),XSTOP,H);
printSol(xSol,ySol,freq)

function F = dEqs(x,y)    % Differential equations.
F = zeros(1,6);
F(1) = y(2); F(2) = y(3); F(3) = y(4); F(4) = y(5); F(5) = y(6);
F(6) = ((5/(9*y(1)^2))-(4*x*y(2)/(9*y(1)^3))-(3*y(2)*y(6)/y(1)));

function y = inCond(u)    % Initial conditions; u(1),u(2),u(3) and u(4) are unknowns
y = [u(1) 0 u(2) 0 u(3) u(4)];

function r = residual(u)  % Bounday residuals.
global XSTART XSTOP H
r = zeros(length(u),1);
x = XSTART;
[xSol,ySol] = runKut5(@dEqs,x,inCond(u),XSTOP,H);
lastRow = size(ySol,1);
r(1) = ySol(lastRow,3);
r(2) = ySol(lastRow,4) - 1;
</code>```

In order to run it we have to define other functions.

runKut5 :

```<code>function [xSol,ySol] = runKut5(dEqs,x,y,xStop,h,eTol)

% 5th-order Runge-Kutta integration.
% USAGE: [xSol,ySol] = runKut5(dEqs,x,y,xStop,h,eTol)
% INPUT:
% dEqs  = handle of function that specifyies the 1-st order differential
% equations : F(x,y) = [dy1/dx dy2/dx dy3/dx ...]

%
% x,y
% xStop = terminal value of x.
% h     = trial value of increment of x.
% eTol  = per-step error tolerance (default = 1.0e-6).
% OUTPUT:
% xSol = x-values at which solution is computed.
% ySol = values of y corresponding to the x-values.

if size(y,1) > 1 ; y = y'; end % y must be row vector

if nargin < 6; eTol = 1.0e-6; end

n = length(y);
A = [0 1/5 3/10 3/5 1 7/8];
B=[    0          0       0          0           0
       1/5        0       0          0           0
       3/40      9/40     0          0           0
       3/10     -9/10    6/5         0           0
     -11/54      5/2   -70/27      35/27         0
    1631/55296 175/512 575/13824 44275/110592 253/4096];
C = [37/378 0 250/621 125/594 0 512/1771];
D = [2825/27648 0 18575/48384 13525/55296 277/14336 1/4];
% Initialize solution
xSol = zeros(2,1); ySol = zeros(2,n);
xSol(1) = x; ySol(1,:) = y;
stopper = 0; k = 1;
for p = 2:5000
    % Compute K’s from Eq. (7.18)
    K = zeros(6,n);
    K(1,:) = h*feval(dEqs,x,y);
    for i = 2:6
        BK = zeros(1,n);
        for j = 1:i-1
            BK = BK + B(i,j)*K(j,:);
        end
        K(i,:) = h*feval(dEqs, x + A(i)*h, y + BK);
    end
    % Compute change in y and per-step error from
    % Eqs.(7.19) & (7.20)
    dy = zeros(1,n); E = zeros(1,n);
    for i = 1:6
        dy = dy + C(i)*K(i,:);
        E = E + (C(i) - D(i))*K(i,:);
    end
    e = sqrt(sum(E.*E)/n);
    % If error within tolerance, accept results and
    % check for termination
    if e <= eTol
        y = y + dy; x = x + h;
        k = k + 1;
        xSol(k) = x; ySol(k,:) = y;
        if stopper == 1;
            break
        end
    end
    % Size of next integration step from Eq. (7.24)
    if e~=0; hNext = 0.9*h*(eTol/e)^0.2;
    else; hNext=h;
    end
    % Check if next step is the last one (works
    % with positive and negative h)
    if (h > 0) == (x + hNext >= xStop )
        hNext = xStop - x; stopper = 1;
    end
h = hNext;
end
</code>```

newtonRaphson2 :

```<code>function root = newtonRaphson2(func,x,tol)

% Newton-Raphson method of finding a root of simultaneous
% equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n.
% USAGE: root = newtonRaphson2(func,x,tol)
% INPUT:
% func = handle of function that returns[f1,f2,...,fn].
% x    = starting solution vector [x1,x2,...,xn].
% tol  = error tolerance (default is 1.0e4*eps).
% OUTPUT:
% root = solution vector.

if nargin == 2
    tol = 1.0e4*eps;
end

if size(x,1) == 1; x = x'; end   % x must be column vector
for i = 1:30
    [jac,f0] = jacobian(func,x);
    if sqrt(dot(f0,f0)/length(x)) < tol
        root = x; return
    end
    dx = jac\(-f0);
    x = x + dx;
    if sqrt(dot(dx,dx)/length(x)) < tol*max(abs(x),1.0)
        root = x; return
    end
end
error('Too many iterations')
function [jac,f0] = jacobian(func,x)
% Returns the Jacobian matrix and f(x).
h = 1.0e-4;
n = length(x);
jac = zeros(n);
f0 = feval(func,x);
for i =1:n
    temp = x(i);
    x(i) = temp + h;
    f1 = feval(func,x);
    x(i) = temp;
    jac(:,i) = (f1 - f0)/h;
end
</code>```

and printSol :

```<code>function printSol(xSol,ySol,freq)

% Prints xSol and ySoln arrays in tabular format.
% USAGE: printSol(xSol,ySol,freq)
% freq = printout frequency (prints every freq-th
%        line of xSol and ySol).

[m,n] = size(ySol);
if freq == 0;freq = m; end
head = '     x';
for i = 1:n
    head = strcat(head,'     y',num2str(i));
end
fprintf(head); fprintf('\n')
for i = 1:freq:m
    fprintf('%14.4e',xSol(i),ySol(i,:)); fprintf('\n')
end
if i  ~= m; fprintf('%14.4e',xSol(m),ySol(m,:)); end
</code>```

Unfortunately, when I run it : I have an error stating that one of the matrices is singular : I assume that there is somewhere a division by 0 that occurs...

```<code>Warning: Matrix is singular to working precision.
> In newtonRaphson2 (line 23)
In shoot4nl (line 14)

Error using newtonRaphson2
Too many iterations

Error in shoot4nl (line 14)
u = newtonRaphson2(@residual,u);
</code>```

Could someone help me with his mathematical expertise in order to solve this problem ?

Thank you in advance,

Best regards.


问:
Implementation of integration schemes for ordinary differential equations in Python and peformance comparison
I look for a book/manual where I can find implementations of different integration schemes for ordinary differential equations (like 4-th order Runge-Kutta) in Python with Numba.

To be more specific, let me provide some details and background. In my research I deal with a system of coupled non-linear differential equations that looks like
$$\frac{dx_i}{dt}=F(g,x_i(t),x_j(t)),$$
where $F$ is non-linear function, $g$ is the parameter and there are $N$ equations. My goal is to obtain the solutions $x_i=x_i(t)$ on the interval $[0,T]$ for all the parameter values from the interval $[g_1,g_2]$ as fast as possible and easily perform basic mathematical manipulations with values $x_i(t=T)$.

Now background:

Currently I use Wolfram Mathematica built-in NDSolve and NDSolveValue functions and use parallelization (Parallelize, ParallelTable). I see that the performance is not good enough for values $N>1000$

Recently I have known that it seems possible to improve performance in Mathematica: one can implement integration scheme by hand and next use Compile (it is quite similar to byte-code compilation, as I understand) and/or use C++ compilation. For me, this program is desirable but I cast some doubts.

These doubts come from the fact that I can use my GPU with CUDA (it is quite "small" card, NVIDIA GTX 1050) and try to improve performance. Next, I am not sure that built-in functions and any manipulations with compilation in Wolfram Mathematica will be faster than another programming language.

From other programming languages, the appropriate choice for me is Python. I know the following opportunities in Python:

odeint from scipy

torchdiffeq package

diffeqpy that operates with Julia

numba in order to wrap numpy-code and make it more fast

So, my ultimate question is:

What is the fastest tool for the stated problem: torchdiffeq or diffeqpy or combination numpy+numba (where integration scheme is written by hand)?

Does this fastest tool allow me to run my code at cluster/server (in future, now I have only my laptop)?


问:
Are there radial basis functions $\phi_n(x)$ which can be factored as $f_n(x) x^{\alpha}$ for $0<\alpha<1$ and $n>0$?
I am looking for a set of radial basis functions $\phi_n(x)$ on $\mathbb{R}^+=[0;\infty[$ which satisfy some othogonality condition

$$ \int x \phi_n(x) \phi_{n'}^*(x) dx =\delta_{nn'}$$

and "behave" like some (fractional) polynomial  $x^{\alpha}$, by which I mean they can be analytically factored as $\phi_n(x)=f_n(x) x^{\alpha}$ for $0<\alpha<1$ and $n>0$ and where $f_n$ is an "uncomplicated" function, i.e. without singularities on $\mathbb{R}^+$ and at least a few times differentiable. Of course, this would mean that $\phi_n(0)=0$ so I am allowing one constant term $\phi_0(x)=\mathrm{const}$.

The idea here is that given the expansion $g(x)=\sum_n c_n \phi_n(x)$ for some $g$ which is compatible with the basis set and for which $c_0=0$, we have $g(x)x^{-\alpha}=\sum_n c_n f_n(x)$ whithout having numerical issues near $x=0$.


问:
optimization problem with L2-norm constraint
I am currently trying to solve a regression problem, which leads me to an optimization problem. Say that we have measured data ($\hat{S}(\omega)\in \mathbb{C}^{N\times N}$), and each entry of this matrix is a $d$-sized vector containing frequency-dependent measurements. We can perform the regression by minimizing

\begin{equation}
\min_A ||S - \hat{S}||_2^2, \quad
\end{equation}

where $ S\in\mathbb{C}^{N\times N}$, and $S_{ij} = MA_{ij}$, where $M\in\mathbb{C}^{d\times k}$, and $A_{ij}\in\mathbb{C}^{k} \ \forall (i,j) \in (1,2, ...,N) \times (1,2,...,N)$. Note that each entry of the model has its respective vector of coefficients $A_{ij}$ that soughts to model the respective $\hat{S}_{ij}$. This is a fairly simple problem so far, just take the derivative of the objective function w.r.t. A and set it to zero. However, the complication arises when the constraint

\begin{equation}
 || S||_2^2 \leq 1
\end{equation}
must be met. I'm clueless as to how to implement this constraint in this problem. Per my understanding, this is a convex optimization problem owed to the L2-norm squared nature of the objective function. However, I'm not entirely sure how to proceed given the aforementioned constraint.


问:
Numerical scheme for the heat equation on the icosahedral hexagonal grid
I have a predefined grid(like this) that is spawned from a regular icosahedradron. It consists of many hexagons and 12 pentagons (corresponding to icosahedradron vertexes). I can tweak the granularity of a grid, but I did not find a way to port this method on a sphere.

I want to solve the heat equation on such a grid:

$\frac{\partial u}{\partial t} - \Delta u = f(x, y, t)$

with some simple border conditions for $u$

I failed to find good resources with numerical schemes and proof of their stability. There is an obvious scheme:

Hexagon:

$\frac{u_i^{n + 1} - u_i^n}{\Delta t} = \frac{1}{6}\cdot \frac{ - 6u_i^n + \sum\limits^{6}_{j = 1}u_j}{d^2} + f^n_i$

Pentagon:

$\frac{u_i^{n + 1} - u_i^n}{\Delta t} = \frac{1}{5}\cdot \frac{ - 5u_i^n + \sum\limits^{5}_{j = 1}u_j}{d^2} + f^n_i$

where $u_j$ are neighbours of $u_i$ and $d$ is a half of center-to-center distnance.

This scheme's error is probably $O(\Delta t) + O(\Delta h^2)$ as it is proofed here behind a paywall for a regular hexagonal grid, not for spherical one with the pentagon though. I bet there are some instabilities that are occurring near pentagons. Maybe they will go away by making the scheme implicit, I don't know.

Another problem is how to store such a grid, I know there is a way to store such a grid if it is on a plane, but I could not port it to a sphere on my own.

And another question that I have is how to extend the hexagonal grid from only the surface to some depth, making it model a spherical layer.

There are some resources for a hexagonal grid on a plane, but it seems like nothing for one on a sphere.

To summarise:

Where to find materials to get good schemas for the heat equation on this grid?

How to store it?

How to extend it to model spherical layer?


问:
scipy.optimize.minimize fails to converge but result is OK
I am trying to optimize a non-linear least squares problem with scipy.optimize.minimize. I have simplified my actual problem down to the case where I am just computing the top 'principal components' like in a PCA analysis and the method reports failure to converge (as it does for my full model). However, the result matches reasonably closely the expected result that you get through SVD decomposition as usually done in PCA. Note that my full problem has no closed-form solution, unlike this simplified one.

My questions are: can I fix the convergence issue? If not, can I ignore the non-convergence status (in the more complicated models)?

In short, I want to find the loadings of the top r PCA components of a matrix, i.e. a k x r orthogonal matrix L such that XL maximizes its variance (or X L L^T is as close to X as possible). I parametrize the orthogonal matrices with a skew-symmetric matrix C and use the matrix exponential map to convert it to L.

I have done the following 'obvious' checks:

NaNs/Infs - I don't see how any could occur and I don't observe any in the outputs.

Needing to rescale the loss function to be not huge/not small. Loss is around 10-20 in the range I am dealing with. The Jacobian is pretty small but non-zero at the final iteration (values around 1e-3 to 1e-5).

Non-differentiable or poorly behaved loss function. It's squared sum residuals and it looks OK to me. The matrix exponential should be well-behaved (finite) in the region where we are in, with the matrix norms all about 1 or lower.

Alternative optimization algorithms. None that I've tried have worked, including BFGS, CG, Newton-CG, and Nalder-Mead (which I think would be too slow for my real data anyway).

Below is a reproducible example on some small simulated dataset. I am using JAX to get the Jacobians (similarly, using jax.scipy.optimize.minimize gives convergence problems too, returning with status=3 indicating that the zoom step of BFGS failed). Disabling JAX jit doesn't change the results.

import jax
import jax.scipy.optimize
from jax import numpy as jnp
import scipy.optimize
import numpy

# Simulate data with one large PC in the first two components
N = 30
scores = numpy.linspace(-1,1, N)
simd = numpy.concatenate([
    [10*scores],
    [10*scores],
    numpy.random.normal(size=(1,N))*5,
    numpy.random.normal(size=(5,N))*0.01
], axis=0).T

# Optimal solution, via SVD
def low_rank_loadings(X, r):
    ''' give best rank r loadings to approximation X '''
    u,d,vt = numpy.linalg.svd(X, full_matrices=False)
    #return u[:,:r] @ numpy.diag(d[:r]) @ vt[:r,:]
    return vt[:r,:].T

# Loss function
@jax.jit
def eval(C, L0, X):
    def func(i, ssr):
        x = X[[i],:]
        L = jax.scipy.linalg.expm(C) @ L0
        return ssr + jnp.linalg.norm(x - x @ L @ L.T)**2
    ssr = jax.lax.fori_loop(
        0,
        len(X),
        func,
        init_val = jnp.asarray([0.]),
    )
    return ssr / len(X)
def extract(vars, L0):
    # Pull out our C matrix from a flattened vector
    k, r = L0.shape
    N = k*r - (r*(r+1)//2) # Num free vars per matrix
    C_ = vars[0:N]
    # Convert to rectangular lower triangular matrices
    idxs = jnp.tril_indices(n=k, k=-1, m=r)
    Ctri = jnp.zeros((k,r)).at[idxs].set(C_)
    # Convert to the right coordinates to match L0
    C = Ctri @ L0.T - L0 @ Ctri.T
    return C

def jax_pca(X, r):
    ''' X = data matrix, r = rank PCA to compute '''
    n,k = X.shape
    L0 = jnp.eye(k, r) # Start with projection to first r components; a bad PCA
    X = jnp.asarray(X)
    def f(vars):
        C = extract(vars, L0)
        return eval(C, L0, X)[0]
    N = k*r - (r*(r+1)//2) # Num free vars per matrix
    res = scipy.optimize.minimize(
        f,
        x0 = numpy.zeros(N),
        jac = jax.jacrev(f),
        method = "BFGS",
    )
    return (L0, extract(res.x, L0), res)
L0, C, res = jax_pca(simd, 1)
print(f"Success? {res.success} -  {res.message}")
loadings = jax.scipy.linalg.expm(C) @ L0
optimal_loadings = low_rank_loadings(simd, 1)
print("Computed:")
print(loadings)
print("Optimal:")
print(optimal_loadings)
```
<p>Which gives output:</p>
```Success? False -  Desired error not necessarily achieved due to precision loss.
Computed:
[[ 6.9389027e-01]
 [ 6.9399828e-01]
 [-1.9204833e-01]
 [-3.3489120e-04]
 [ 3.8082333e-05]
 [ 1.3443755e-04]
 [ 1.8192278e-04]
 [ 2.3052096e-04]]
Optimal:
[[ 6.93945195e-01]
 [ 6.93945195e-01]
 [-1.92041436e-01]
 [-3.36661681e-04]
 [ 3.91581503e-05]
 [ 1.34102010e-04]
 [ 1.81875549e-04]
 [ 2.31636442e-04]]
```


问:
Complexity of Branch-and-cut algorithm in terms of "Big O"
How can I compute the Big O complexity of the Branch and cut algorithm?
I am solving an integer linear program using MOSEK that includes $M$ binary variables, but I do not know how to calculate the exact complexity.
Is it just O($2^M$)? or are there problems solved at each branching node that have extra complexity to be computed?


问:
Numerical algorithms made stable by unums which are unstable on IEEE floats
For unums, there is good evidence (see figure 5) that accuracy is better than IEEE floats. (Note: I use the term "unum" broadly to refer to any of the various iterations and revisions to the format.)

However, generally we regard "stability as more of a constraint than accuracy". Are there numerical algorithms which can be made stable using some unum format that is unstable using IEEE floating point?


问:
cuSolver Sparse solver fails with CUSOLVER_STATUS_INTERNAL_ERROR
I'm endeavouring to use cuSolver to solve a sparse linear system.

But when I run my code, it fails with CUSOLVER_STATUS_INTERNAL_ERROR.

What's going wrong? How does one debug this error?

My MWE follows and uses two data files (chosen somewhat arbitrarily):

mat65k.mtx - Manually remove the lines beginning with %

rhs65025.txt

Code

// Compile with: nvcc src/cuda_solver.cu -lcusolver -lcusparse

#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <vector>

#include <cuda.h>
#include <cuda_runtime.h>
#include <cusolverSp.h>
#include <cusparse.h>
#include <eigen3/Eigen/Sparse>

#define CUSOLVER_CHECK(err) (HandlecusolverError(err, __FILE__, __LINE__))
#define CUDA_CHECK(err) (HandleError(err, __FILE__, __LINE__))

static void HandlecusolverError(cusolverStatus_t err, const char* file, int line) {
  if (err != CUSOLVER_STATUS_SUCCESS) {
    fprintf(stderr, "ERROR: %d in %s at line %d, (error-code %d)\n", err, file, line, err);
    fflush(stdout);
    exit(-1);
  }
}

static void HandleError(cudaError_t err, const char* file, int line) {
  if (err != cudaSuccess) {
    fprintf(stderr, "ERROR: %s in %s at line %d (error-code %d)\n", cudaGetErrorString(err), file, line, err);
    fflush(stdout);
    exit(-1);
  }
}

#define CUSPARSE_CHECK(func)                                                         \
  {                                                                                \
    cusparseStatus_t status = (func);                                                \
    if (status != CUSPARSE_STATUS_SUCCESS) {                                         \
      std::cerr << "Got cusparse error"                                              \
                << " at line: " << __LINE__ << " of file " << __FILE__ << std::endl; \
    }                                                                                \
  }

struct SparseSolverCUDA {
 private:
  cusolverSpHandle_t cusolver_handle = nullptr;
  csrqrInfo_t info                   = nullptr;
  cusparseMatDescr_t matrix_desc     = nullptr;
  cudaStream_t stream                = nullptr;
  const int32_t width                = 0;

  // We wish to solve the sparse linear system Ax=b where A is in CSR form
  double* d_csr_val_A    = nullptr;  // Values of the A matrix
  int32_t* d_csr_row_ptr = nullptr;  // Rows of the A matrix
  int32_t* d_csr_col_idx = nullptr;  // Columns of the A matrix
  double* d_b            = nullptr;  // Values of the b vector
  double* d_x            = nullptr;  // Values of the x vector

 public:
  using eigen_spmat_t = Eigen::SparseMatrix<double, Eigen::RowMajor, int32_t>;

  SparseSolverCUDA(const eigen_spmat_t& eigen_spmat) : width(width) {
    if (eigen_spmat.rows() != eigen_spmat.cols()) {
      throw std::runtime_error("SparseSolverCUDA matrix must be symmetric.");
    }

    // step 2: create cusolver handle, qr info and matrix descriptor
    CUSOLVER_CHECK(cusolverSpCreate(&cusolver_handle));

    CUDA_CHECK(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking));
    CUSOLVER_CHECK(cusolverSpSetStream(cusolver_handle, stream));

    CUSPARSE_CHECK(cusparseCreateMatDescr(&matrix_desc));

    // Options:
    //   CUSPARSE_MATRIX_TYPE_GENERAL
    //   CUSPARSE_MATRIX_TYPE_SYMMETRIC
    //   CUSPARSE_MATRIX_TYPE_HERMITIAN
    //   CUSPARSE_MATRIX_TYPE_TRIANGULAR
    CUSPARSE_CHECK(cusparseSetMatType(matrix_desc, CUSPARSE_MATRIX_TYPE_GENERAL));
    // CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE
    CUSPARSE_CHECK(cusparseSetMatIndexBase(matrix_desc, CUSPARSE_INDEX_BASE_ZERO));

    CUSOLVER_CHECK(cusolverSpCreateCsrqrInfo(&info));

    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&d_csr_val_A), eigen_spmat.nonZeros() * sizeof(double)));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&d_csr_row_ptr), eigen_spmat.outerSize() * sizeof(int32_t)));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&d_csr_col_idx), eigen_spmat.innerSize() * sizeof(int32_t)));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&d_b), eigen_spmat.rows() * sizeof(double)));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&d_x), eigen_spmat.rows() * sizeof(double)));
  }

  ~SparseSolverCUDA() {
    CUDA_CHECK(cudaFree(d_csr_row_ptr));
    CUDA_CHECK(cudaFree(d_csr_col_idx));
    CUDA_CHECK(cudaFree(d_csr_val_A));
    CUDA_CHECK(cudaFree(d_b));
    CUDA_CHECK(cudaFree(d_x));
    CUDA_CHECK(cudaStreamDestroy(stream));
    CUSPARSE_CHECK(cusparseDestroyMatDescr(matrix_desc));
    CUSOLVER_CHECK(cusolverSpDestroy(cusolver_handle));
  }

  void operator()(const eigen_spmat_t& A, const Eigen::VectorXd& b, std::vector<double>& x_out) const {
    // step 3: copy Aj and bj to device

    assert(A.rows() == b.size());
    assert(A.rows() == x_out.size());

    std::cerr << "type = " << typeid(A.innerIndexPtr()).name() << std::endl;

    CUDA_CHECK(
        cudaMemcpyAsync(d_csr_val_A, A.valuePtr(), A.nonZeros() * sizeof(double), cudaMemcpyHostToDevice, stream));
    CUDA_CHECK(cudaMemcpyAsync(
        d_csr_col_idx, A.innerIndexPtr(), A.innerSize() * sizeof(int32_t), cudaMemcpyHostToDevice, stream));
    CUDA_CHECK(cudaMemcpyAsync(
        d_csr_row_ptr, A.outerIndexPtr(), A.outerSize() * sizeof(int32_t), cudaMemcpyHostToDevice, stream));

    CUDA_CHECK(cudaMemcpyAsync(d_b, b.data(), A.rows() * sizeof(double), cudaMemcpyHostToDevice, stream));

    std::cerr << "A.rows() = " << A.rows() << std::endl;
    std::cerr << "A.nonZeros() = " << A.nonZeros() << std::endl;

    CUDA_CHECK(cudaStreamSynchronize(stream));

    int32_t singularity;
    const auto sparse_ret_code = cusolverSpDcsrlsvqr(
        cusolver_handle,
        A.rows(),
        A.nonZeros(),
        matrix_desc,
        d_csr_val_A,
        d_csr_row_ptr,
        d_csr_col_idx,
        d_b,
        0.0001,  // tol
        0,       // No ordering if reorder=0. Otherwise, symrcm, symamd, or csrmetisnd is used to reduce zero fill-in
        d_x,
        &singularity);

    if (sparse_ret_code != CUSOLVER_STATUS_SUCCESS) {
      std::string error_code_str;
      switch (sparse_ret_code) {
        case CUSOLVER_STATUS_SUCCESS:
          error_code_str = "The operation completed successfully.";
          break;
        case CUSOLVER_STATUS_NOT_INITIALIZED:
          error_code_str = "The library was not initialized.";
          break;
        case CUSOLVER_STATUS_ALLOC_FAILED:
          error_code_str = "The resources could not be allocated.";
          break;
        case CUSOLVER_STATUS_INVALID_VALUE:
          error_code_str =
              "Invalid parameters were passed (m,nnz<=0, base index is not 0 or 1, reorder is not 0,1,2,3)";
          break;
        case CUSOLVER_STATUS_ARCH_MISMATCH:
          error_code_str = "The device only supports compute capability 2.0 and above.";
          break;
        case CUSOLVER_STATUS_INTERNAL_ERROR:
          error_code_str = "An internal operation failed.";
          break;
        case CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED:
          error_code_str = "The matrix type is not supported.";
          break;
      }

      throw std::runtime_error("Sparse solver failed with " + std::to_string(sparse_ret_code) + " " + error_code_str);
    }

    std::cerr << "Singularity = " << singularity << std::endl;

    CUDA_CHECK(cudaMemcpyAsync(x_out.data(), d_x, A.rows() * sizeof(double), cudaMemcpyDeviceToHost, stream));

    CUDA_CHECK(cudaStreamSynchronize(stream));
  }
};

int main() {
  CUDA_CHECK(cudaSetDevice(0));

  int width, height, nnz;
  std::ifstream fin("mat65k.mtx.txt");
  fin >> width >> height >> nnz;

  std::cerr << "Populating matrix A" << std::endl;

  std::cerr << "nnz = " << nnz << std::endl;

  std::vector<Eigen::Triplet<double>> triples;
  triples.reserve(nnz);
  for (int i = 0; i < nnz; i++) {
    int row, col;
    double val;
    fin >> row >> col >> val;
    triples.emplace_back(row, col, val);
  }

  Eigen::SparseMatrix<double, Eigen::RowMajor> mat(height, width);
  mat.setFromTriplets(triples.begin(), triples.end());

  std::cerr << "rows = " << mat.rows() << std::endl;
  std::cerr << "cols = " << mat.cols() << std::endl;
  std::cerr << "nnz = " << mat.nonZeros() << std::endl;
  std::cerr << "mat.innerSize() = " << mat.innerSize() << std::endl;
  std::cerr << "mat.outerSize() = " << mat.outerSize() << std::endl;

  fin.close();

  std::cerr << "Populating RHS" << std::endl;
  fin.open("rhs65025.txt");
  Eigen::VectorXd b(width);
  for (int i = 0; i < height; i++) {
    double val;
    fin >> val;
    b(i) = val;
  }

  std::cerr << "Compressing matrix A" << std::endl;
  mat.makeCompressed();

  auto ssc = SparseSolverCUDA(mat);

  std::vector<double> x(width);
  ssc(mat, b, x);

  return 0;
}
```</span>
</code>```


问:
Numpy random generator seed choice
I am new to understanding pseudorandom number generators and sometimes feel daunted by them. For some time I have used the following approach in my work:

In one file, I'd have

```<code>import numpy as np
np.random.seed(1)
my_random_process()
</code>```

In another file, I'd have

```<code>import numpy as np
np.random.seed(2)
my_random_process()
</code>```

and so on. In the nth and final file, I'd have

```<code>import numpy as np
np.random.seed(n)
my_random_process()
</code>```

In my head, each of these processes were "independent" in the sense that if I averaged my outcomes over all the file's outputs, I would have an estimate for the mean of my outputs with an error decreasing as $c/\sqrt{n}$ for a constant $c$. My aim was merely to find the average of the outputs of my_random_process, which depends on np.random functions.

However, I've been reading the documentation and I fear that the method above actually has correlations. In particular, numpy's parallel random sampling documentation says that np.random uses MT19337 for which "two adjacent 32-bit integer seeds (i.e. 12345 and 12346) would produce very similar streams". This makes me very sad, because correlations would inflate the value of the constant $c$.

On the other hand, numpy has other methods of generating random numbers, one of which I'm not too familiar with, np.random.default_rng. On another page of the random sampling documentation, there's the sentence "Seeds can be passed to any of the BitGenerators. The provided value is mixed via SeedSequence to spread a possible sequence of seeds across a wider range of initialization states for the BitGenerator."

Does this mean that if I do

```<code>import numpy as np
rng = np.random.default_rng(1)
my_new_random_process(rng)
</code>```

and

```<code>import numpy as np
rng = np.random.default_rng(2)
my_new_random_process(rng)
</code>```

and so on, where I've replaced np.random with rng.random (and so on) in my_random_process, will the random streams be largely uncorrelated? Or does one still need to be a little more thoughtful to avoid correlations? This is not a large problem for me, since I just want to calculate means, but if a small change in my workflow could help reduce the number of iterations I have to do, I'd be happy.

Title: Numpy random generator seed choice


问:
How can I improve this matrix product calculation in OpenCL?
I am trying to compute a matrix-matrix product of N stacked complex double N x N matrices. For simplicity, I assume N = 512. I have written code in C++ parallelized with OMP and using OpenBLAS for the matrix-matrix product. This is the code:

```<code>#include <iostream>
#include <complex>
#include <chrono>
#include <openblas/cblas.h>

typedef std::complex<double> cx;

int main() {
    const int NUM_REPS = 10; // Number of speed test repetitions

    const int N = 512;
    const int N3 = N * N * N;
    const int N2 = N * N; // Distance between matrices

    cx *A = new cx[N3];
    cx *B = new cx[N3];
    cx *C = new cx[N3];

    cx alpha = 1.0;
    cx beta = 0.0;

    for (int i = 0; i < NUM_REPS; i++) {

        #pragma omp parallel for
        for (int j = 0; j < N3; j++) {
            A[j] = cx(1, 1);
            B[j] = cx(0, 1);
        }

        auto start = std::chrono::steady_clock::now();

        #pragma omp parallel for
        for (int k = 0; k < N; k++) {
            cblas_zgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans, N, N, N,
                        &alpha,
                        &(A[k * N2]), N,
                        &(B[k * N2]), N,
                        &beta,
                        &(C[k * N2]), N);
        }

        auto end = std::chrono::steady_clock::now();

        std::cout << "Delta [ms] = "
                  << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
                  << std::endl;

        cx checksum = 0.0;
        for (int j = 0; j < N3; j++)
            checksum += C[j];
        std::cout << checksum << std::endl;
    }

    return 0;

}
</code>```

Note that in the example, the first matrix should be transposed before the multiplication and the memory layout is row major.

The "checksum" for the output I receive is (6.87195e+10,6.87195e+10). On my laptop, each iteration in i takes about 2.381s.

Since this feels a little bit long, I was wondering whether this can be improved by using the GPU for the computation. The main problem here is that a 512³ matrix stack does not fit into the GPU memory all at once since it takes 2 GB of memory per stack. My GPU (Nvidia GeForce MX450) has exactly 2 GB of memory, so that will not work. The only solution would be to copy the stack in chunks into the GPU memory, perform the multiplication of this sub-stack and copy the results back into host memory which means that it becomes an "out of core" operation, see also for example this answer, which describes splitting a large matrix into blocks and multiplying blockwise on the GPU.

For this, I am using the CLBlast library, more specifically the GemmStridedBatched method. Here is the full code:

```<code>#include <iostream>
#include <complex>
#include <chrono>
#include <clblast.h>

#include "cl.hpp"

typedef std::complex<double> cx;

int main()
{
    const int NUM_REPS = 10; // Number of speed test repetitions

    const int N = 512;
    const int N3 = N * N * N;
    const int N2 = N * N; // Distance between matrices
    const int block_len = 128; // Length of each block, must divide N

    const auto platform_id = 0;
    const auto device_id = 0; // This selects the Nvidia GPU

    cx *A = new cx[N3];
    cx *B = new cx[N3];
    cx *C = new cx[N3];

    cx alpha = 1.0;
    cx beta = 0.0;

    // Get OpenCL platform and device
    std::vector<cl::Platform> platforms = std::vector<cl::Platform>();
    cl::Platform::get(&platforms);
    if (platforms.size() == 0 || platform_id >= platforms.size()) {
        std::cout << "No platform." << std::endl;
        return -1;
    }
    cl::Platform platform = platforms[platform_id];
    std::vector<cl::Device> devices = std::vector<cl::Device>();
    platform.getDevices(CL_DEVICE_TYPE_ALL, &devices);
    if (devices.size() == 0 || device_id >= devices.size()) {
        std::cout << "No device." << std::endl;
        return -1;
    }
    cl::Device device = devices[device_id];

    auto device_as_vector = std::vector<cl::Device>{device};
    auto context = cl::Context(device_as_vector);
    auto queue = cl::CommandQueue(context, device);
    cl_event event;

    auto device_a2 = cl::Buffer(context, CL_MEM_READ_ONLY, N2 * block_len * sizeof(cx));
    auto device_b = cl::Buffer(context, CL_MEM_READ_ONLY, N2 * block_len * sizeof(cx));
    auto device_c = cl::Buffer(context, CL_MEM_WRITE_ONLY, N2 * block_len * sizeof(cx));

    // Repeat the speed test NUM_REPS times
    for (int i = 0; i < NUM_REPS; i++) {
        for (int j = 0; j < N3; j++) {
            A[j] = cx(1, 1);
            B[j] = cx(0, 1);
        }

        auto start = std::chrono::steady_clock::now();

        // Matrix-matrix product over blocks of length block_len
        // in GPU memory
        for (int k = 0; k < N; k += block_len) {
            queue.enqueueWriteBuffer(device_a2, CL_TRUE, 0, N2 * block_len * sizeof(cx), &(A[k * N2]));
            queue.enqueueWriteBuffer(device_b, CL_TRUE, 0, N2 * block_len * sizeof(cx), &(B[k * N2]));

            auto queue_plain = queue();
            auto status = clblast::GemmStridedBatched(
                clblast::Layout::kRowMajor,
                clblast::Transpose::kConjugate, clblast::Transpose::kNo,
                N, N, N,
                alpha,
                device_a2(), 0, N, N2,
                device_b(), 0, N, N2,
                beta,
                device_c(), 0, N, N2,
                block_len,
                &queue_plain, &event);

            // Record the execution time
            if (status == clblast::StatusCode::kSuccess) {
                clWaitForEvents(1, &event);
                clReleaseEvent(event);
            }
            else {
                std::cout << "(1) FAILURE: " << (int)status << std::endl;
                return -1;
            }

            queue.enqueueReadBuffer(device_c, CL_TRUE, 0, N2 * block_len * sizeof(cx), &(C[k * N2]));
        }

        auto end = std::chrono::steady_clock::now();

        std::cout << "Delta [ms] = "
                << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count()
                << std::endl;

        cx checksum = 0.0;
        for (int j = 0; j < N3; j++)
            checksum += C[j];
        std::cout << checksum << std::endl;
    }

    return 0;

}
</code>```

Here, block_len defines the number of matrices from the stack which lives in the GPU memory at the same time. This number must divide N so as not to "miss" any matrices. (In a more general implementation this can of course be adapted to allow arbitrary stack sizes.)

The cl.hpp header can be found here. I have used the OpenGL code from this example from the CLBlast library as a basis for this code.

My problem is, that this version actually takes  10.838s per iteration in i. That is unexpected and I honestly don't know if I can then make the code actually faster than 2s per iteration using OpenBLAS. However, since I am new to OpenCL there may be some way to actually improve the speed of my code. This would be my first question. More generally, is there maybe another way to improve the speed of this stacked matrix-matrix product? I might have missed something important.


问:
How to accelerate the computing of implicit finite difference method for heat conduction between two solids
Edit on May 3rd: I have found the problem. Because the difference of between $k_1$ and $k_2$ is huge, a very small time step need to be chosen so that the right green part can "feel" the heat from the left. Otherwise the temperature profile in the green region would barely change. However, this requres a lot of computation resources especially when I non-dimensionalize the problem. In the following figure, the blue line indicate the condiction that $k_1=k_2$, whereas the red line is $k1<<k2$. In obtaining the blue line, I chose the non-dimensional time step $\Delta\tau$ to be 6e-11 to pass the grid sensitivity test, but as you can imagine this took a lot of time. So is there any technic to speed up the calculation? Can I compute the two sides separately?

I am solving a transient heat condcution problem involving 2 solids using implicit finite difference method.

At the interface, I soppose the finite difference equation at the interface should be
$$ \frac{k_1(T_{i-1}-T_i)}{\Delta x}=\frac{k_2(T_{i}-T_{i+1})}{\Delta x}$$
or
$$\frac{1}{2}(\rho_1Cp_1+\rho_2Cp_2)\frac{T^{n+1}_i-T^n_i}{\Delta t}=k_1\frac{T^{n+1}_{i-1}-T^{n+1}_i}{\Delta x^2}+k_2\frac{T^{n+1}_{i+1}-T^{n+1}_i}{\Delta x^2}$$
The thernmal conducvivities of the 2 solids are 0.5 and 400, and $\rho_1Cp_1$ and $\rho_2Cp_2$ are assumed to be equal. But when I solved the problem using the following code (in Matlab), the result is unexpected. That is, the heat transfer performance of the 2 solids with $k_1=0.5$ and $k_2=400$ is worse than $k_1=k_2=0.5$.

So I am wondering if there is any special criteria need to be applied if the difference between $k_1$ and $k_2$ is huge?

```<code>clear all
L=1;
d=0.5;%model dimension
dx=d/40;%grid size

GN=d/dx+1;%interface position
N=L/dx+1; %total grid number

rho1=900;
cp1=2000;
k1=0.5;
rho2=900;
cp2=2000;
k2=400;
a1=k1/(rho1*cp1);
a2=k2/(rho2*cp2);%thermomechanical properties

Time=100000; %total time
dt=1;%time step
alpha=dx^2/dt;

T0=273.15+25;%initial temperature
Tb=273.15+100;%boundary temperature
Ntime=round(Time/dt);%time steps
theta=zeros(Ntime,N);%temperature field history
T=T0*ones(1,N);%initial temperature
T(:,1)=Tb;%LHS B.C.
T2=T;
R=1;

for p=1:1:Ntime
    p
theta(p,:)=T;
while (R>10e-5)
  T3=T2;
  for i=2:1:GN-1
      T2(i)=(a1*(T2(i-1)+T2(i+1))+alpha*T(i))/(alpha+2*a1);
  end
  T2(GN)=(k1*T2(GN-1)+k2*T2(GN+1))/(k1+k2);
  for i=GN+1:1:N-1
      T2(i)=(a2*(T2(i-1)+T2(i+1))+alpha*T(i))/(alpha+2*a2);
  end
  T2(N)=(a2*(T2(N-1)+T2(N-1))+alpha*T(N))/(alpha+2*a2);

  R=max(abs(T3-T2)) ;

end
T=T2;
R=1;
end
</code>```


问:
The mathematical meaning of a zero gradient pressure boundary condition in the Navier-Stokes equations
I would like to solve the Navier-Stokes equations for the unsteady problem of the flow around a circular cylinder. I would like to understand how to write mathematically the boundary condition for the pressure on the surface of the cylinder (essentially the circle). The domain of the problem is shown in the image below with the cylinder having a diameter D and centered in the origin.

I usually solve this problem using the CFD solver OpenFOAM, therefore, I will describe the boundary conditions for this problem which are common among OpenFOAM users. At the inlet the velocity $u$ is $(1,0)$ and zero gradient pressure is imposed. On the sides the so-called slip condition is imposed for velocity and zero gradient is used for pressure, while on the outlet we have zero gradient for velocity and zero pressure. In OpenFOAM, zero gradient pressure is often used for the obstacles in the domain such as the cylinder in this case and a zero velocity is assumed. Whenever a zero gradient condition is mentioned, what is meant is zero gradient in the direction perpendicular to the flow (for example on the outlet only the $x$-derivative is equal to zero for the velocity components $U$ and $V$ that is $U_x = V_x = 0$).

For the case of the cylinder boundary, is it correct to assume that the scalar product between the pressure gradient $\nabla p = (p_x,p_y)$ and the normal vector at the boundary $\boldsymbol{n} = \frac{(x,y)}{0.5D}$ has to be zero, in other words:

$$
\nabla p \cdot \boldsymbol{n} = x p_x + y p_y = 0
$$

Is that the correct formulation of the pressure boundary condition ? and is there another physical condition for the pressure at the cylinder boundary ?


问:
Is there a mathematical algorithm that can determine whether a problem is parallelizable?
Is there a mathematical algorithm that can determine whether a problem is parallelizable?

Possibly similar to how the definition for "computable" works out:

https://en.wikipedia.org/wiki/Computable_function

So a parallelizable program is a program that ...?


问:
Lattice Boltzmann methods beyond BGK-style collision operators
I am trying to educate myself on the lattice Boltzmann method to see if it may be useful for simulating some problems in plasma physics that I am interested in. One thing that strikes me is the (seemingly) universal adoption of BGK-style collision operators. Is the method intractable with more sophisticated collision operators like full Boltzmann (for neutral gases) or Fokker-Planck (for charged particles)? I understand that BGK-style operators are useful pedagogically for learning about the method, but does this also represent the state of the art?</p>


问:
FEM applied to heat equation and incompatible conditions
Consider the problem $$u_t - \Delta u = f \text{ on } \Omega\times (0,T)\\u=0  \text{ on } \partial \Omega\times (0,T) \\ u(x,0)=g(x) \text{ on } \Omega$$
with $g$ NOT vanishing on the boundary. If I use a FEM discretization in space, what is the best a priori error (considering both space and time) that I can hope for? Under which time discretization? Is there anything else involving finite elements (e.g. space-time FEM) that would improve the situation?

Indeed, all the nice error estimates I came across with, assume compatibilty of initial and boundary condition...

Thank you in advance.

Cross posted also on Mathoverflow: here


问:
Reference request for finite elements theory
Consider a domain $\Omega \subseteq \mathbb{R}^{2,3}$ which is non convex and with $C^2$ boundary.  Could you recommend a good reference where it is explained how:

without needing isoparametric elements

with a triangulation of simplices

one can obtain:

some space of linear finite elements $S^h_0$ which vanish on $\partial \Omega$ (the reference should clearly explain how to define the hat functions outside of $\Omega$ and what properties the triangulation should satisfy...)

a suitable definition of the nodal interpolator $I_h:H^2(\Omega)\cap H^1_0(\Omega)\rightarrow S^h_0$

such that the approximation property $||I_hv-v||_{L^2}+h||\nabla(I_hv-v)||_{L^2}\leq C h^s , ||v||_{H^s}$ holds for $s=1,2$?


问:
How can I reduce the artifact in "Thin Plate Spline" interpolation?

At the Top "right", there is the 2D-density plot of the recorded data
(actual), fewer in number. Recorded data has been sampled a on the 8
arms of a regular octagon. These 8 arms are placed along x and y axes
and along the diagonals. On the top-left is the interpolated data
(thin-plate) of the actual data, i.e. larger in number. Compared to
the Right (actual) data-set we can observe that "Interpolated" data is
showing some kind of "Quadrupole" kind of (which also looks like wings
of a  "Butterfly") artifact which is not desired and expected. I was
expecting kind circular kind of "Interpolated" density plot as can be
shown at the feet of the bottom figures.

I have tried both my own written-code and in-built python module to perform "Thin-Plate-spline" interpolation and both the results are showing same "Quadrupole" kind of artifact in the interpolated data (see the fig. at top-left corner). Codes are attached below with both in-bult scipy.interpolate() module and my own code as well. In my own code I have solved the equation as described below,

https://elonen.iki.fi/code/tpsdemo/

My question is  why is "Quadrupole" or "Butterfly" sort of artifact is
arising instead of a rather circular or octagonal pattren (which is
seen in actual data, Top-Right corner Figure) in "Thin Plate Spline"
interpolation ?

How shall I be able to reduce that ? For example,The top-left figure
is the actual data and top-right one is the corresponding
"Interpolated" data which is different kind of pattern, right ?. This
implies that I might be making some errors in the process of
interpolation, but where am I committing the mistake is not clear to
me.

Now,I have tried 3 different codes to have the interpolated data. And,
in each case I found the artifact to be present there in the
interpolated data.

2 of these 3 codes have used "in-built"

scipy.interpolate()-thin_plate spline whereas I also ran my own code
of "Thin-Plate-Spline". I was hoping if there is any mistake in my
codes and thus I posted my codes along with the question. I was hoping
that it might aid reviewer to understand where does the onset of
"Quadrupole" or "Butterfly" kind of "artifact" lie. I would like to
request if any one could kindly response and guide me about what has
to be done. Kindly, please have your indispensable expertise.

CODE - 1 => Density-Plot

```<code>import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

###########################################################
#np.random.seed(1977)
#x, y, z = np.random.random((3, 10))
###########################################################
data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,:])/100
y = (data[1,:])/100
E1 = data1[0,:]
###########################################################

interp = scipy.interpolate.Rbf(x, y, E1, function='thin_plate')

#yi, xi = np.mgrid[0:1:100j, 0:1:100j]

###########################################################
X = np.linspace(-300, 300, 1200)
Y = np.linspace(-300, 300, 1200)
xi, yi = np.meshgrid(X, Y)
###########################################################

Ei = interp(xi, yi)

###########################################################
plt.plot(x, y, 'ko')
plt.imshow(Ei, extent=[0, 1, 1, 0], cmap='gist_earth')
plt.colorbar()
###########################################################
plt.show()
</code>```

CODE - 2 => 3D-Surface Plot

```<code>import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

#############################################################
from mpl_toolkits.mplot3d import axes3d  #Axes3d
import matplotlib.cm as cm
from pylab import rcParams
rcParams ['figure.figsize'] = 5, 5
#########################################################################################################

fig = plt.figure()
axes = fig.gca(projection='3d') #gca = get current axis

###########################################################
#np.random.seed(1977)
#x, y, z = np.random.random((3, 10))
###########################################################
data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,:])/100
y = (data[1,:])/100
E1 = data1[0,:]
###########################################################

interp = scipy.interpolate.Rbf(x, y, E1, function='thin_plate')

#yi, xi = np.mgrid[0:1:100j, 0:1:100j]

###########################################################
X = np.linspace(-300, 300, 1200)
Y = np.linspace(-300, 300, 1200)
xi, yi = np.meshgrid(X, Y)
###########################################################

Ei = interp(xi, yi)

###########################################################
#p = axes.plot_surface(xi,yi,Ei, rstride=4, cstride=4, cmap=cm.RdBu, linewidth=0, antialiased=False)
p = axes.plot_surface(xi, yi, Ei,cmap='viridis', edgecolor='none')
fig.colorbar(p, shrink=0.5)
plt.tight_layout();
fig.savefig("surface.pdf")
###########################################################
plt.show()
</code>```

CODE - 3 => My own code - 3D

```<code>import numpy as np
import math as mt

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d  #Axes3d
import matplotlib.cm as cm
from pylab import rcParams
rcParams ['figure.figsize'] = 5, 5

#########################################################################################################

fig = plt.figure()
axes = fig.gca(projection='3d') #gca = get current axis

#########################################################################################################
print(" ")
################################################################

#data = np.loadtxt('data1.dat',unpack=True)
#x = data[0,:]
#y = data[1,:]
#E1 = data[2,:]

################################################################

data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,16:])/100
y = (data[1,16:])/100
E1 = data1[0,16:]

################################################################
################################################################
n1 = len(x)
#print(" n1 : ", n1)
r = np.stack((x, y), axis=-1)
print(r)

N = n1 + 3

o = np.zeros(3)

#################################################################

#################################################################

V1 = np.hstack((E1, o))
V = V1.reshape(N, 1)

##########################################################################################################

C = np.ones(n1)

P = np.stack((C, x, y), axis=-1)
P_T = P.transpose()  # Transpose of P

#print(" P : ")
#print(P)

#print(" P_T : ")
#print(P_T)

#########################################################################################################

O1 = np.zeros(9)
O = O1.reshape((3,3))

#print(" O : ")
#print(O)

#########################################################################################################

K = np.zeros((n1,n1))

for i in range(n1):
    for j in range(n1):

        if(i == j):
            K[i][j] = 0.0
        else:
            R = r[i] - r[j]
            s = np.sqrt( (R[1]*R[1]) + (R[0]*R[0]) )
#           print(s)
            K[i][j] =  (s*s) * (np.log(s))

#print(" U[i][j]")
#print(K)

##########################################################################################################

#L1 = np.block([K, P])
L1 = np.hstack((K, P))
#L2 = np.block([P_T, O])
L2 = np.hstack((P_T, O))
#L = np.block([[L1],[L2]])
L = np.vstack((L1,L2))
#print(" L : ")
#print(L)

LI = np.linalg.inv(L)

###########################################################################################################

W = np.dot(LI, V)   # Weights

print(" ")
#print(" W : ")
#print(W)

###########################################################################################################

def f(x, y, r, N, W):

    f = W[N-3] + (W[N-2])*x + (W[N-1])*y

    for i in range(N-3):

        r1 = np.sqrt(((r[i][0] - x)**2) + ((r[i][1] - y)**2))
#       r1 = (((r[i][0] - x)**2) + ((r[i][1] - y)**2))
        if(r1==0.0):
            U1 = 50     #IMPORATANT
        else:
            U1 = (r1*r1) * (np.log(r1))

        f = f + W[i]*U1

    if(f < 0.0):
        f = 0.0

    return f

#f9 = f(35, 35, r, N, W)  # value of f9 is in good approximation with the value achieved from simulation
#print(f9)

## upto this point CODE is working well

print("_______________________________________________________________________________________")
print(" ")

#F = f(xv, yv, x, y, N, W) # This is not working

x1 = np.linspace(-200, 200, 1000)
y1 = np.linspace(-200, 200, 1000)
xv, yv = np.meshgrid(x1,y1)

f_xy = np.zeros((1000, 1000))

for i in range(1000):
    for j in range(1000):
        f_xy[i][j] = f(x1[i], y1[j], r, N, W)

#############################################################################################################################

p = axes.plot_surface(xv,yv,f_xy, rstride=4, cstride=4, cmap=cm.RdBu, linewidth=0, antialiased=False)
##p = axes.plot_surface(xv, yv, f_xy,cmap='viridis', edgecolor='none')
fig.colorbar(p, shrink=0.5)
plt.tight_layout();
fig.savefig("surface.pdf")
plt.show()
</code>```


问:
Closed (Robin) boundaries in advection-diffusion equation with FDM
I am solving the equation
$$
\frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x} \left( D \frac{\partial \phi}{\partial x} + v\phi \right)
$$
using finite differences. I want to include Robin boundaries, that is, boundaries that are closed to both advection and diffusion. I read this excellent question on it:

Conservation of a physical quantity when using Neumann boundary conditions applied to the advection-diffusion equation

However, the actual solution to the problem is not given there.
Robin boundary conditions on the left edge of the domain should read
$$
\Big(D \frac{\partial \phi}{\partial x} + v\phi \Big)|_{x=0} = j_\mathrm{robin} = 0
$$
I have tried discretizing this as
$$
D \frac{\phi_2-\phi_0}{\Delta x} + v\phi_1 = j_\mathrm{robin}
$$
With a fictious cell with index 0 outside of the domain.

I can reorder this to get
$$
\phi_2 - (j_\mathrm{robin} - v \phi_1) \frac{2 \Delta x}{D} = \phi_0
$$
The discretized equation for the left cell - with the usual 2nd order central scheme for the diffusion term and a central scheme for the advective scheme -  reads

$$
\frac{\phi_1^{(t+1)}-\phi_1^{(t)}}{\Delta t} = \frac{D}{\Delta x^2} \Big(\phi_0 - 2 \phi_1 + \phi_2 \Big) + \frac{v}{2\Delta x}(\phi_2 - \phi_0)
$$
I have omitted the upper indizes for the terms on the right hand side - I use an higher-order explicit scheme.
We can eliminate the values for the index 0, yielding

$$
\frac{\phi_1^{(t+1)}-\phi_1^{(t)}}{\Delta t} = \frac{2D}{\Delta x^2} (\phi_2 - \phi_1) + (j_\mathrm{robin} - v\phi_1)(\frac{v}{D}-\frac{2}{\Delta x})
$$

I have tried this, basically with a Gaussian traveling into the domain boundary, and while there is a notable difference between this and the Neumann boundary (that would not be closed to advection) I nevertheless lose all mass.
My question is basically whether or not this is the right approach to begin with.
Thank you!


问:
Advance a Interpolation
Note; No special knowledge of Pykrige is needed to answer the question, as I already mention examples in the question!

Hi I would like to use Universal Kriging in my code. For this I have data that is structured as follows:

```<code>      Latitude   Longitude  Altitude          H2        O18        Date    Year  month       dates  a_diffO       O18a
0    45.320000  -75.670000     114.0  -77.500000 -11.110000  2004-09-15  2004.0    9.0  2004-09-15   -0.228 -10.882000
1    48.828100    9.200000     314.0  -31.350000  -4.880000  2004-09-15  2004.0    9.0  2004-09-15   -0.628  -4.252000
2    51.930000  -10.250000       9.0  -18.800000  -3.160000  2004-09-15  2004.0    9.0  2004-09-15   -0.018  -3.142000
3    48.248611   16.356389     198.0  -45.000000  -6.920000  2004-09-15  2004.0    9.0  2004-09-15   -0.396  -6.524000
4    50.338100    7.600000      85.0  -19.200000  -3.190000  2004-09-15  2004.0    9.0  2004-09-15   -0.170  -3.020000

</code>```

You can find my data here:https://wetransfer.com/downloads/9c02e4fc1c2da765d5ee9137e6d7df4920220618071144/8f450e

I want to interpolate the data (Latitude, Longitude, Altitude and O18) with Universal Kriging and use the height as a drift function.

So far I have programmed this here but I am not getting anywhere, e.g. I don't know how to effectively use the height as a drift function and the information from the Pykrige documentation is of limited help:

```<code>from traceback import print_tb
from typing_extensions import Self
import numpy as np
from pykrige.uk import UniversalKriging
from pykrige.kriging_tools import write_asc_grid
import pykrige.kriging_tools as kt
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Path, PathPatch
import pandas as pd
from osgeo import gdal

def load_data():
    df = pd.read_csv(r"File")
    return(df)

def get_data(df):
    return {
        "lons": df['Longitude'],
        "lats": df['Latitude'],
        "values": df['O18'],
    }

def extend_data(data):
    return {
        "lons": np.concatenate([np.array([lon-360 for lon in data["lons"]]), data["lons"], np.array([lon+360 for lon in data["lons"]])]),
        "lats": np.concatenate([data["lats"], data["lats"], data["lats"]]),
        "values":  np.concatenate([data["values"], data["values"], data["values"]]),
    }

def generate_grid(data, basemap, delta=1):
    grid = {
        'lon': np.arange(-180, 180, delta),
        'lat': np.arange(np.amin(data["lats"]), np.amax(data["lats"]), delta)
    }
    grid["x"], grid["y"] = np.meshgrid(grid["lon"], grid["lat"])
    grid["x"], grid["y"] = basemap(grid["x"], grid["y"])
    return grid

def interpolate(data, grid):
    uk = UniversalKriging(
        data["lons"],
        data["lats"],
        data["values"],
        variogram_model='exponential',
        verbose=True,
        drift_terms=["point_log"],
    )
    return uk.execute("grid", grid["lon"], grid["lat"])

def prepare_map_plot():
    figure, axes = plt.subplots(figsize=(10,10))
    basemap = Basemap(projection='robin', lon_0=0, lat_0=0, resolution='h',area_thresh=1000,ax=axes)
    return figure, axes, basemap

def plot_mesh_data(interpolation, grid, basemap):
    colormesh = basemap.contourf(grid["x"], grid["y"], interpolation,32, cmap='RdBu_r') #plot the data on the map. plt.cm.RdYlBu_r
    color_bar = basemap.colorbar(colormesh,location='bottom',pad="10%")

df = load_data()
base_data = get_data(df)
figure, axes, basemap = prepare_map_plot()
grid = generate_grid(base_data, basemap, 90)
extended_data = extend_data(base_data)
interpolation, interpolation_error = interpolate(extended_data, grid)
plot_mesh_data(interpolation, grid,basemap)
plt.show()

</code>```

I now only use universal kriging and create these images:

I get the expected error: ValueError: Must specify location(s) and strength(s) of point drift terms.

I just know that I have to create a grid with the height, but I don't know how and I don't know how to make the drift dependent on the altitude. The altitude formula is:

where 100 m represents 100 m hight difference.

The interesting thing is that there is this website with examples: however, I am too inexperienced in coding to understand the examples and to transfer them to my example: https://python.hotexamples.com/examples/core/-/calc_cR/python-calc_cr-function-examples.html

I've been trying to solve these problems for 3 weeks now, but I'm really getting nowhere.

The question, how can I create an external drift for the UK interpolation that works as a function of altitude?

Title: Advance a Interpolation


问:
schrodinger eq time propagation with dissipation using split step operator
I am looking in ways to include energy dissipation while propagating a coherent wavepacket in a 1d TDSE. for example I use the split step method: exp[Δt(D+V)]≈exp[ΔtV/2]exp[ΔtD]exp[ΔtV/2], and subtract some epsilon value from the kinetic energy term every step. What other ways can I dissipate energy from the system so it'll relax to its equilibrium position?</p>


问:
Can you describe the Galerkin numerical method to solve the wave equation?
How would you describe the Galerkin method to solving the 3D wave equation
$$u_{tt}= c^2\Delta u$$
to someone who wants to implement it immediately?

More precisely, we want to solve the Cauchy problem
$$u_{tt}= c^2\Delta u , \\  u_t(x,y,z,0)= g(x,y,z) , \\ u(x,y,z,0)= f(x,y,z)$$
where $f, g$ are given compactly supported functions.

Say in a bounded region of the space we are interested in, say a ball of radius $R$, note that the wave propagates towards the boundary of the region, the sphere, thus we can't impose the condition at the boundary.

How would you implement this method ?
I am looking for a description of the finite elements, how to decompose this region? (you can use any other shape, not necessarily a ball, take a cube for example), then describing the numerical algorithm necessary to find the values of the solution $u$ at certain element or node $(x,y,z,T)$ at time $T$


问:
Reformulate a problem with concave objective function into a QP
I would like to convert this problem into a QP (Quadratic program).

$$\text{Maximize } \sum_{k=1}^{K}\sum_{n=1}^{N}log2(1+p_{kn}b_{kn})\\
\text{subject to } \sum_{k=1}^{K}\sum_{n=1}^{N}p_{kn}\leq P_{0}\\
p_{kn}\geq 0, \forall k,n$$

The variables are $p_{kn}, \forall k,n$
Since the optimization variables are continuous and are inside the logarithm function, I was unable to convert it into a QP. Can you help me to resolve this issue please?


问:
Book recommendation request: A good C++ based software development book for a computational mathematician
I'm hoping that this type of question is on-topic, and I have put in a little hunting to make sure it isn't a flat out repeat. If it is either of the above, I apologize in advance.

I'm an early career computational mathematician who is at least minimally capable with C++, and pretty proficient with something like MATLAB, so the act of writing code itself is not generally a major issue for me. I most often use libraries developed by others, but even when coding doesn't hold me back I'm constantly finding myself losing hour after hour to things such as: Not using the correct commands when I'm trying to download something from GitHub, not using the correct file structure, not giving path to my compiler correctly, not configuring dependent libraries properly, not building a library properly, not building my own code properly and so on.

It is almost certainly a situation of I don't know what I don't know, but I have the feeling a lot of the pain comes from skills not directly related to coding, but a light form of software development. I would like to learn about CMake, using compliers properly, how to download and link a library properly, and similar things. When a given software is well-documented and provides me with an example CMakeLists.txt file then I can sort it out. If it does not, then I might as well not even bother downloading the library because I won't figure it out.

Thus I'm asking for a book that would teach me those things from a standpoint of what I need to know and not too, too much else (we all want to get the code working and do the more fun stuff, right?) Bonus points if it covers git. I feel it would be a good experience for me to try to contribute to some of these libraries as I become capable. It doesn't have to be a book necessarily, but a mostly self-contained resource would be great. Trying to fit several pieces of the solution together from multiple sources, is what I have been doing and has lead me here.

Thanks for the help in advance.

Edit/Addendum: Maybe another piece of the puzzle here is that I lack knowledge of bash commands and probably also have a lot to learn about Ubuntu. Again, I'm hoping I've given enough information to give an idea of the type of things I feel hold me back. So that even though I don't know what I don't know, you all will be able to provide some helpful suggestions.


问:
Describing a vector field
I am currently making use of a 2D vector field that I would like to be able to describe and characterize. I thought of a few metrics that I could use, but they seem limited for both a qualitative and quantitative analysis.

I believe that the divergence and curl of the field would be interesting. The vector field could also be mapped with a 2D polynomial function as well as harmonics using 2D Fourier Transform.

Do you happen to know other tools to describe or decompose a vector field?

Title: Describing a vector field


问:
Recursion relations for integrating Gaussian functions
I'm trying to implement a numerical method used in quantum chemistry from scratch. I'm using this paper as a reference. It's also available on Sci-Hub. So, the method requires calculating integrals of each spatial dimension (x, y, z) using recursive relations but I got completely lost on how it can be implemented in the code. Can someone, please, guide me?

Here, the authors suggest starting from term (11). The question that I'm having here is what do 0 arguments mean? Is this the starting member of a functional sequence?
After calculating that, how is it related to (15) and (16)? How do I use (11) to calculate (15) and (16)? Given both (15) and (16) correspond to 2nd and 4th arguments equal 0. I understand that (15) is used to "step up" in index 1 and (16) is used to step up in index 3, but how is it supposed to be implemented? It also seems to be a fair question, what exactly do i, j, k and l mean in this function if not the basis function number? What do values n and m mean, as well?

Eventually, the algorithm leads to calculating (19) and (18) with (18) being the final value of the term. But how does it result from (15), (16) and (11)?

I would be very grateful if one provides me some pseudocode or Python/C++ functions. The method has been implemented many times already, however it's usually either Fortran or very high-level C++ which I can't understand. Because I can't understand it I'm making it from scratch. It has to be a compact and a portable code.


问:
Can someone explain why RK4 is less accurate for very small timesteps?
I am currently working on a project where I have used an RK4 integrator to attempt to solve the three-body problem. An interesting result that I found, was that decreasing the size of the time steps only improved accuracy until a certain point where the accuracy started to decrease. Is this normal for RK4 integrators or have I done something wrong? if this is normal, can anyone explain why this behaviour occurs?


问:
Unit comparison between two molecular dynamics simulations with different time step
I need some clarification regarding the time step used in the molecular dynamics simulations. In the first simulation the time step was 1 nanosecond and in the second simulation the time step in 10 ns. Now I want to compare the results of these two simulations. For example take mean-square displacement.
$$
MSD(1ns) = \frac{nm^{2}}{ns} = 10^{-9} m^{2}/s \\
MSD(10ns) = \frac{nm^{2}}{10 ns} = 10^{-8 } m^{2}/s
$$
What I did above is correct or I should not include 10ns in the denominator of the second MSD?


问:
Getting gradient from voxel, is it possible and how?
The R package rmarchingcubes implements a marching cubes algorithm. Given a voxel contructed from a function defining an isosurface, it returns a 3d mesh with some vertex normals. When I plot this surface, it is incredibly smooth. If instead of the rmarchingcubes normals, I use the classical normals obtained by averaging the surrounding faces, the plot is not as smooth. If I use the gradient of the function to get the normals, it is smooth.

So it is clear that rmarchingcubes does not compute the classical vertex normals. But when I run it, I don't give the function, I only give the voxel. So I'm wondering whether it is possible to approximate the gradient from the voxel? And if it is possible, how?


问:
Does an alias method exist for sampling a discrete distribution that is slighlty modified at each iteration?
I have the following problem. I must sample from a discrete distribution that is changing at each sort.

Let me explain, with a "vivid" description, I draw a color ball from a bag. The probability associated with each color is different (for this part I could choose an alias method, for example, the Vose's algorithm). BUT, after I choose a paint can, I paint the ball with this color and I put it back in the bag. So, If I use a classical alias method, I have to initialize prob and alias each time and it takes time. So, I wonder if an algorithm derived from the alias method exists for this kind of problem. As only two probabilities are changed at each iteration, I could imagine that the initialization part could be done in a very efficient way.


问:
Stretched Elastic Sheet with Horizontal Cut in Interior
Given a rectangular finite elastic sheet $ABCD$ containing an arbitrary horizontal cut (or discontinuity) $EF$ in its interior. The edge $AB$ of the sheet is tethered but $CD$ is stretched by $\Delta$ such the horizontal cut is pulled down and "stretched open" as shown below.

I would like to determine the function $x = F(X)$ that maps any point $X$ in the unstretched configuration to a point $x$ in the stretched configuration, assuming the stretching of the sheet does not tear the horizontal cut any further. I am not interested in the stresses on the sheet only its deformation. I would also be satisfied with an approximate "strength of materials" type solution.


问:
Finding the weak form of a PDE with a tensor argument
I am trying to solve for the order parameter ($A$) in the Ginzburg Landau equations. I had asked on the math SE site but was recommended to ask here.
We are trying to solve the following equation, (with the implied sum over $j$, for each component identified by $\mu$, $i$), $K_i$, $\beta_i$, and $\alpha (T)$ are known parameters (constants wrt position),
$$
K_1 \partial^2_j A_{\mu i} + K_{23} \partial_i (\partial_j A_{\mu j}) = 2 \beta_1 Tr(AA^T) A_{\mu i} + 2 \beta_2 Tr(AA^\dagger)A_{\mu i} + 2 \beta_3[AA^TA]_{\mu i} + 2 \beta_4[AA^\dagger A]_{\mu i} + 2 \beta_5[AA^TA]_{\mu i} + \alpha(T) Tr(AA^\dagger) = (\text{rhs})_{\mu i},
$$
assuming,
$$
A =
\begin{pmatrix}
A_{uu} & A_{uv} & A_{uw} \\
A_{vu} & A_{vv} & A_{vw} \\
A_{wu} & A_{wv} & A_{ww}
\end{pmatrix}.
$$

I have a code in C++ that implements the FDM with relaxation, but we have found that our mixed derivative approximations have large enough error that the code doesn't always converge to a solution. We are looking at the option of using others' FEM solvers (like FreeFEM or MOOSE), but I'm having a hard time getting the weak form of our set of equations. I am following the description here, and obtained
$$
0 = K_1 \oint_{\Gamma} {\psi \vec{\nabla}A_{\mu i} \cdot \hat{n}}
- K_1 \int_{\Omega} {\vec{\nabla} \psi \cdot \vec{\nabla}A_{\mu i}}
- \int_{\Omega} {\psi (\text{rhs})_{\mu i}}
+ K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})},
$$
but I don't know how to get the last term. I saw this post and thought that maybe I'd have to explicitly write out all 9 equations and use 9 test functions?

My main question is: can this system of equations be reasonably written for a FEM solver? If so, how?

I played around with the equations a little more and figured I could write them generally as,
$$
(\text{rhs})_{\mu i} = K_1 \vec{\nabla}^2 A_{\mu i} + K_{23} \partial_i \Big[ \vec{\nabla} \cdot \vec{A}_{\mu;r} \Big],
$$
where $\vec{A}_{\mu;r}$ is the vector formed from the $\mu$-th row of $A$ (is there a better way to represent that?). (A note here: the rows and columns of the order parameter do not necessarily correspond to spatial variables, and for our ultimate goal will be in cylindrical coordinates; so it may not be a good idea to represent this as a vector.) Thus,
$$
0 = K_1 \oint_{\Gamma} {\psi \vec{\nabla}A_{\mu i} \cdot \hat{n}}
- K_1 \int_{\Omega} {\vec{\nabla} \psi \cdot \vec{\nabla}A_{\mu i}}
- \int_{\Omega} {\psi (\text{rhs})_{\mu i}}
+ K_{23}\int_{\Omega} {\psi \partial_i \Big[ \vec{\nabla} \cdot \vec{A}_{\mu;r} \Big]},
$$
which might be easier (or more obvious how) to integrate by parts? So, would I still have to use 9 test functions? We expect the solution function for each element to be different, so I might still need to have all 9 separately...

Trying a slightly different approach: taking the last term, $K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})}$, explicitly writing out the sum over $j$, and using the product rule, we can say,
$$
K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})} \to K_{23} \sum_{j=u,v,w}{ \Bigg[ \int_{\Omega}{\partial_i (\psi \partial_j A_{\mu j})} - \int_{\Omega}{(\partial_i \psi)(\partial_j A_{\mu j})} \Bigg] }
$$

The first term here ($\int_{\Omega}{\partial_i (\psi \partial_j A_{\mu j})}$) looks like the divergence term that we previously rewrote as a boundary integral...except that there are no vectors. Is this now in the right form so that I can use it in a FEM solver?

After applying the product rule, I found that I could select a different order of derivatives (because these derivative operators commute, right?) and instead obtain,
$$
\to K_{23} \Bigg[ \int_{\Omega}{ \vec{\nabla} \cdot (\psi \partial_i \vec{A}_{\mu;r})} - \sum_{j=u,v,w}{ \int_{\Omega}{(\partial_j \psi)(\partial_i A_{\mu j})} \Bigg] }
$$
Now the first term could be rewritten using the divergence theorem? (Again, assuming that $\vec{A}_{\mu;r}$ can be treated as a "normal" vector.)


问:
How to maximize effectiveness of 2 check bits overseeing 6 data bits?
Having a single parity bit gives a Hamming distance of 2, so only one single bit can be corrupted to reliably detect an error.

The next step I know of would be a Hamming code, which uses at least 3 parity bits.

However, is there a method which uses two parity bits and performs better than just a single parity bit?

I don't know of any error checking method which uses 2 bits. They either use 1, or at least 3 bits.

A naive approach would be to use one parity bit for the first half of the byte, and the other one for the second half of the byte. However, this won't detect two consecutive corrupted bits (unless it's right in the middle of the byte)

Another approach would be to have a regular parity bit, and the second bit would be a parity bit for every other bit. This would correctly detect two consecutive errors, but in case of randomly distributed errors it wouldn't be any better than using a single parity bit.

Is there a better approach? Is my assumption correct, that in the worst case two parity bits are no better than only one?


问:
Lowest eigenvalues of a matrix with divergent entries
In high energy physics we oftentimes encounter the following problem. For a given parametrized matrix $\{H_{ij}(\Lambda)\}$, we know that in the limit $\Lambda\to\infty$ some of its entries become large:

$$
\text{For some $(i,j)$:  }H_{ij}(\Lambda)\to\infty\text{ as }\Lambda\to\infty.
$$

We also know that in the limit $\Lambda\to\infty$ the lowest eigenvalues in the spectrum of $\{H_{ij}(\Lambda)\}$ — which we are most interested in — tend to converge to particular values (while the larger eigenvalues may potentially diverge).

Have we had access to an "infinite-precision" computer, we could just diagonalize the matrix for increasing values of $\Lambda$ extrapolate to $\Lambda\to\infty$. In reality, this is impractical.

I am wondering if there exist any well-known ways to deal with such a situation, at least for some specific cases (such as divergent elements staying on the diagonal and matrix being sparse).


问:
Does loss landscape provide a means to compare optimization algorithms?
I'm comparing two optimization algorithms for deep neural networks. To that end, I train a network with the same data and starting rom the same initial weights using two different optimization algorithms. Then, I visualize the optimization trajectory within the loss landscape using this method.

Besides the qualitative observation that how trajectory moves w.r.t. the loss level-sets, are there any quantitive measures to compare the two methods?


问:
positivity preservation for mixed finite element discretization
I'm interested in mixed discretizations of the diffusion (and related) equations:

$$\begin{align}
\frac{\partial h}{\partial t} + \nabla\cdot \mathbf q & = f \\
k^{-1}\mathbf q + \nabla h & = 0
\end{align}$$

where $h$ is a scalar field representing thickness or hydraulic head and $\mathbf q$ is a vector field representing flux or velocity.
The primal form of this problem is obtained by eliminating $\mathbf q$, giving a second-order PDE for $h$.
When discretizing the primal form with degree-1 continuous Galerkin basis functions, one can prove that, when the triangulation is nice, the resulting matrices are M-matrices and thus the problem has a discrete maximum principle.
See for example the Kuzmin book.

Is the maximum principle still true with some mixed discretizations?
For example, the MINI element $CG_1 \times (CG_1 + B_3)^2$ is inf-sup stable for mixed Poisson; does it have a discrete maximum principle for the thickness?
What about DG $\times$ Raviart-Thomas (or some other hipster element) for the flux?
The Godunov theorem tells us that there can't be a discrete maximum principle for higher-order basis functions, but there are still low-order bases for which it isn't necessarily an obstruction.
Does the analysis change appreciably if the constitutive relation is nonlinear?


问:
Adding wall friction to a 3D model - Stokes equation
Maybe that is not the best place for the following question, if such, let me know and sorry.

I would like to model the flow of expanding bentonite in a fracture following the given reference. In order to take into account the swelling pressure, they modified the Stokes equation (momentum equation) by adding the swelling pressure explicitly (equation 4-6 of reference):
\begin{equation}
\frac{\partial \bf{u}}{\partial t} -\nabla \cdot \bf{\tau} + \nabla (p+p_{sw}) = 0 \\
\end{equation}

After that they add the static friction at the wall (equation 4-7 in reference):
\begin{equation}
\frac{\partial \bf{u}}{\partial t} -\nabla \cdot \bf{\tau} + \nabla (p+p_{sw}(1-2\mu)) = 0 \\
\end{equation}

the friction coefficient is \begin{equation} \mu \end{equation}
I wonder if such assumption in a 3D model would be correct. I assume that if I solve the system with the last equation, the nodes that are not close to the wall are also going to take into account the friction coefficient. I think it has more sense to add the Friction force into the wall as a boundary condition. I thought to use a navier slip condition.

Does it make sense what I am saying, or I am wrong.

Thanks

Best Regards

Dani

Ref:

Pont, A., Coene, E., & Idiart, A. (2020). Bentonite erosion project. Svensk kärnbränslehantering AB [Swedish Nuclear Fuel and Waste Management Company]}


问:
How to solve the following SDP with cvxpy in Python?
The SDP problem is
$$
\min_{Z \in S^{n},Y \in S^{m}} {\rm trace}(Z) +{\rm trace}(Y)\\
{\rm s.t.} \begin{bmatrix}
Y & X\\
X^T & Z
\end{bmatrix} \succeq 0\\
X \in C
$$
Where $C$ is a convex set.But I have some difficulty in transforming it to the standard form:
$$
\min_{X}\operatorname{trace}(CX) \\
\text{s.t.} \ \operatorname{trace}(A_{i}X)=b_i,i=1,...,p \\
X\succeq 0
$$
If let $\overline X= \begin{bmatrix}
Y & X\\
X^T & Z
\end{bmatrix}$ ,then the problem transform to
$$
\min_{\overline X} {\rm trace}(\overline X)\\
{\rm s.t.} \overline X \succeq 0\\
\overline{X}_{m+i,j} =\overline{X}_{j,m+i},i=1,...,n,j=1,...,m
$$
But $\overline{X}_{m+i,j} =\overline{X}_{j,m+i}$ seems can not writen as ${\rm trace}(A_i\overline{X})=b$,where $A_i$ symmetric matrice
This question comes form pape: A Rank Minimization Heuristic with Application to
Minimum Order System Approximation


问:
Unconstrained convex optimization: correlation between dimensionality and Lipschitz constant
The author of the SIAM News article "Optimization Theory and Perspectives on the Field of Machine Learning" mentions:

... For unconstrained convex optimization, GD (gradient descent) convergence rates do not depend on the space’s dimensionality. In other words, if we measure a notion of how close we are to the global minimum, convergence rates depend on the function’s properties and the number of descent steps; they are independent of whether a problem has 10 variables or 10 million. Readers might find this dimension-independence surprising, though perhaps there is a “catch” that practical problems may have a correlation between dimensionality and their Lipschitz constant.

I was wondering if there are any research papers/reports/studies that support the claim (in bold) that practical optimization problems have a correlation between dimensionality and their Lipschitz constant.

This may be a general study that performs an analysis on a large set of various problems or a narrow application with a dimensionality-Lipschitz constant study. It would be particularly helpful if the study is related to computational electromagnetics, but I am also interested in other domains. A more general optimization theory reference that would justify this [quite intuitive] claim (as optimization is not my main area of expertise) would be also quite welcome.


问:
How to classify ODE equilibria that are stable but slowly changing in value with time?
I'm numerically solving a system of coupled ODEs where time is the independent variable. At each time, I can solve for the equilibrium values of my state variables where their respective derivatives go to zero and classify the equilibria as stable or unstable. For one of my variables, the equilibrium value is stable at all times but gradually changes with time. My actual solution also tends to stay near this equilibrium -- i.e., my solution variable also gradually changes in value to keep up with the stable equilibrium.

My question is: what is the formal classification name for this kind of equilibrium? The best I can come up with is a wordy description: the equilibrium is "locally stable" (at any given time, a perturbation away from it would bring the variable toward the equilibrium) but at the same time also "globally unstable" since the equilibrium value is not constant in time but rather evolving. I should say that after evolving the system for a long time, the equilibrium does become constant with time so in the end it is also "asymptotically stable". Just wondering if I am using these equilibrium/stability classification terms correctly.


问:
Reference for preconditioning nonlinear conjugate gradient with LU decomposition of jacobian
When solving non-linear systems of equations, @Arnold Neumaier suggested in How can I precondition a non-linear problem before linearization? to use an approximate LU decomposition of an approximate Jacobian at the starting point to precondition nonlinear conjugate-gradient to minimize $\|F(x)\|_2^2$. Can anyone point out a reference for this scheme, ideally with convergence results?