问:
Order of local error when integrating ODE with discontinous derivatives
I'm working with ODEs, $$\dot{x} = f(x, t),$$ where the (higher) derivatives of the right-hand side have discontinuities. In particular, $f(x, t)$ is obtained by interpolation of discrete samples, and the number of continious derivatives thus depend on the order of interpolation.

I know from Theorem II.3.1 on page 157 of Hairer, Nørsett & Wanner that if a Runge-Kutta method is of order $p$, then the local error is bounded by $C h^{p+1}$ if all partial derivatives of $f(x,t)$ up to order $p$ exist and are continuous.

Hence, I know that if I for example use cubic spline interpolation to obtain $f(x, t)$, then only the first and second derivatives are continuous, and hence I cannot expect fourth-order accuracy from a fourth-order Runge-Kutta method. In practice, though, I do see fourth-order convergence when $f(x, t)$ is obtained by cubic spline interpolation, while if $f(x, t)$ is obtained with linear interpolation (where not even the first derivative is continuous) I get second-order convergence from fourth-order Runge-Kutta.

Hairer, Nørsett & Wanner only tells me that I cannot expect a method of order $p$ to work as advertised unless I have $p$ continous partial derivatives, but they don't tell me how large the error will be. I've done some numerical experiments, but I'm curious to see if there exists some more rigorous literature on this point.

My question is therefore if there exists any (reasonably accessible) theory that I can use to reason about this.


问:
Quadrature rules with the weight function $w(x) = |x|^\gamma$
I am interested in integrals of the form
$$
\int_{[0,1]^{d}} |x|^{\gamma}f(x)dx.
$$

$\gamma>0$ and $f$ has some singular behavior at $\vec{0}$. The weight function $|x|^\gamma$ is commonly used in weighted Sobolev spaces which include reduced regularity functions.

For $d=1$, a Gauss-Jacobi quadrature rule with $(\alpha,\beta) = (0,\gamma)$ can be used to approximate the above integral. For higher dimensions, I don't think this is the case as $|x|^\gamma \neq \prod_{i=1}^{d} x_{i}^{\gamma}$. The choice $(\alpha,\beta)=(0,0)$ (i.e. Gauss-Legendre) does not work well for me since $|x|^{\gamma}f(x)$ still contains a singularity (depending on the value of $\gamma$).

I'm wondering if anyone has literature suggestions for integrals of the above form, even just for 2D would be sufficient. I couldn't find anything in my search even though I thought this would be a weight function someone's looked at before.

ContentLicense: CC BY-SA 4.0


问:
Is it normal for a genetic algorithm to not have mating phase?
I know that is possible for a genetic algorithm to not have a mating phase because some genetic algorithms use a "mutation-only" approach, where new candidate solutions are generated by randomly mutating existing solutions without crossover.

But what about creating a genetic algorithm that does not include crossover and mutation? Is that a good approach?

I am thinking of creating a recommendation service that retrieves audio files from the database.
This service uses an algorithm that calculates the fitness score of each individual based on some properties of that file and after that, I sort the individuals based on the fitness score and get the best candidate.

But at this point isn't that just a "scoring" algorithm and not a genetic one? Or that's perfectly fine to call it like that?


问:
Does exact diagonalization of a matrix allow for efficient computation of a Lanczos basis?
Suppose that we are given a large, real-symmetric matrix $L$, which is simply too large to perform exact diagonalization on numerically. If we want to study its spectrum, one tool we can use is the Lanczos algorithm. This algorithm takes $L$ and an initial vector, call it $q_0$, and uses $q_0$ to produces a tridiagonal representation of $L$. A clear review of the algorithm is available on wikipedia. Call the orthonormal basis which spans this subspace the "Lanczos basis" with respect to $q_0$. Usually, this tridiagonal form and the Lanczos basis are just byproducts produced along the way towards computing eigenvalues of $L$.

I am in the strange position where my goal is to invert this procedure. For my problem, $L$ is a very large matrix that I happen to have good analytic control over, so I can compute its eigenvalues and eigenvectors by hand. My question is, is it possible to use this knowledge to efficiently construct the Lanczos basis and corresponding tridiagonalization of $L$ with respect to (a given but in principle arbitrary) $q_0$?

It's hard to imagine that knowledge of the complete spectrum doesn't give a computational advantage, but so far I have not been able to figure out how to use it. It's an unfortunate position to be in, because the Lanczos algorithm suffers from severe numerical stability issues and I would rather avoid running it; there should be a basis transformation which relates the Lanczos basis of $L$ with respect to $q_0$ to its eigenspace, but I don't know how to construct it. Any ideas from linear algebra practitioners?

Edit: It was asked in the comments (1) if the initial Lanczos vector is fixed or arbitrary and (2) what exactly is needed from the output of my computations. The initial Lanczos vector (call it $q_0$) I'll look at will depend on the precise application, but I think it's sufficiently generic to consider the case where $q_0$ is not an eigenvector of $L$.

As for the output that I need, what I am ultimately after are the off-diagonal elements of $L$ when it is put in tridiagonal form. That is, I don't need specific knowledge of the Lanczos basis per se, only the elements of $L$ in that basis. If that makes things simpler, that's great, I just don't see how we could get one without the other.

Edit 2: There has been some confusion over what constitutes a valid tridiagonal form; as pointed out in the comments, all diagonal matrices are tridiagonal. I have clarified my question above to emphasize that I am after the Lanczos basis with respect to a particular vector, $q_0$. This uniquely specifies a tridiagonal form of $L$, but unless $q_0$ happens to be an eigenvector of $L$, the eigenspace of $L$ is not a valid Lanczos basis of $L$ with respect to $q_0$.


问:
Finite volume method using Chebyshev polynomials
I want to solve the following set of coupled advection-diffusion equations:
$$
\frac{\partial f}{\partial t}=\nabla\cdot(\kappa\nabla f)+\nabla\cdot(\boldsymbol{u}f)+s_f(g),
$$
$$
\frac{\partial g}{\partial t}=\nabla\cdot(\kappa\nabla g)+\nabla\cdot(\boldsymbol{u}g)+s_g(f),
$$
in polar coordinates $(r,\theta)$. Using finite elements is natural in that case. However, to help for stability and computation time, I would like to have a grid equally spaced in $\cos\theta$ and the zeros of Chebyshev polynomials for $r$. I never tried using finite volume methods so I am a bit unsure of what I am doing. My question is this: what is the best way to go about this? Can I simply use Chebyshev interpolation and derivation for values on volume faces? Or can I solve the system for the Chebyshev coefficients instead? I plan to use a fully explicit time stepping scheme (probably Runge-Kutta second order)


问:
Discontinuous Galerkin for transport equation with non-constant advection
This question is mainly an inquiry about the usefulness of Discontinuous Galerkin (DG) for the time-independent transport equation of the form
$$\sigma u+\beta\cdot\nabla u =f,\;\;\;\text{on }\Omega\subset\mathbb{R}^2$$
when the advection field $\beta$ is a non-constant function of $x$ (or $u$). Boundary conditions are
$$ u=g\;\;\;\text{ on the 'inflow' boundary }\;\;\Gamma_-:=\{x\in\partial\Omega:\beta(x)\cdot\mathbf{n}(x)<0\}. $$
Now let $\mathcal{T}_h$ be a triangulation of $\Omega$ and $V_h=\{v\in L^2(\Omega):v_{|K}\in\mathbb{P}_d(K),\forall K\in\mathcal{T}_h\}$ be the usual DG space. Then the local formulation on each triangle $K\in\mathcal{T}_h$ is given by
$$\sigma\int_K u_hv +\int_K (\beta\cdot\nabla u_h)v-\int_K fv=\int_{\partial K^-} (\beta\cdot\mathbf{n}_K)(u_h-\xi_h)v,$$
where $\partial K^-$ is in the inflow boundary and
$$\xi_h=\begin{cases}0 &\text{ on }\;\;\partial K^-\cap \Gamma_-,\\ \beta\cdot{n}_K(x)\lim_{\epsilon\to 0^+}u_h(x-\epsilon\beta)&\text{ on }\;\;\partial K^-\setminus\Gamma_-. \end{cases} $$

In the case of when $\beta$ is constant, Lesaint and Raviart [1] showed that there is an enumeration of the triangles following the characteristic direction $\beta$ which renders the matrix of the discrete problem block triangular.
I understand the interest of DG over CG (Continuous Galerkin) here, where it is easier to invert the system matrix, thus more or less justifying the extra d.o.f (degrees of freedom).

However, I don't think this advantage remains when $\beta$ becomes non-constant and has loops (if $\beta$ is an electric field for instance). Finding a proper enumeration becomes harder and there is more global coupling. So why should I use DG with its extra d.o.f? Am I missing some fix that makes DG advantageous again?


问:
Deviation between Analytic DFT and FFT in Python
Within my work, I am trying to compare analytically retrieved power spectra with ones calculated from fft packages in python.

The problem I have, is that the analytic form of the peaks I derived does not match the fft result.

The minimal example I tested for that is the fourier transform of a cosine, in which I only consider positive frequencies. Analytically, I get via a geometric sum:

\begin{align} \mathcal{F}[\cos(2 \pi n dt)] &= \sum_{n=0}^{N-1} \frac{1}{2}(e^{-i2\pi (\nu -1)ndt}+e^{i2\pi(\nu+1)ndt}) \\
&= \frac{1}{2}\left( e^{i2\pi(\nu-1)(N-1)dt} \frac{\sin(\pi(\nu-1)T)}{\sin(\pi(\nu-1)dt)} + \dots\right).
\end{align}

As I am only interested in positive frequencies, I leave out the second term. Now when I compare this with numerical implementations, explicitly:

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

T = 10
Nsamp = int(T*8)
t = np.linspace(0,T,Nsamp)
nu = np.linspace(0,Nsamp/(2*T),int(Nsamp/2))

sig = np.cos(2*np.pi*t)

ft1 = np.abs(sc.fft.fft(sig))[:Nsamp//2]
ft2 = np.abs(np.fft.fft(sig))[:Nsamp//2]
ft3 = 1/2*np.abs(np.sin(np.pi*(nu-1)*T)/np.sin(np.pi*(nu-1)*T/Nsamp))
</code>```

I get noticeable deviations from between the two, as shown here:

What could be the reason for this?


问:
What is the upper bound of the size of the neighborhood of a vertex in a twin-free $K_{1,5}$-free interval graph?
Suppose that I have an interval graph that does not contian $K_{1,5}$ (the star on 6 vertices) as an induced subgraph and that I only care about twin-free graphs, i.e., no two vertices have the same neighborhood. What is the maximal degree of such a graph? Is it bounded by a constant?

As a remainder, interval graphs are chordal and without asteroidal triples.

Is it roughly 30?


问:
Can we get the exact solution of large-scale quadratic programming problems (quadratic objective, linear inequality constraints) using KKT condition?
Crossposted at MathOverflow

Consider a quadratic programming problem with the following format:
$$
\text{min} Q(x) = c^Tx+\frac{1}{2}x^TDx \\
$$
$$
\text{s.t.} Ax\leq b, \\
x\geq 0
$$
where $D$ is a $n\times n$ positive definite matrix and $A$ is a $m\times n$ matrix. We know when $n$ and $m$ are small, for example, when both of them are 2, we can solve the problem and get the exact solution using the KKT condition. However, when $n$ and $m$ are large (meaning that we have a large-scale quadratic programming problem. I think we are also able to get the exact solution using the KKT condition.

However, to the best of my knowledge, I cannot find any reference for this. The book 'quadratic programming' said that this type of question can be solved in polynomial time (page 3), for example, the interior point method. But as far as I know, the interior point method only can get the near-optimal solution making the objective function within the acceptable tolerance (like $10^{-7}$).
We can also refer to this post about using the interior point method to solve the quadratic programming method.

My question: is it possible for us to get the exact solution to the aforementioned question when $n$ and $m$ are large? How expensive it is if we want to get the exact solution (global optimum) of the quadratic programming problem I mentioned before? And how it relates to the number of $m$? Is it impossible to calculate it even with computers? I will really appreciate it if you can offer some references for this problem.
I also posted a question in "mathoverflow", but it seems that there are no answers there, and I received a suggestion to post my question here.


问:
Is the Hessian of the strain energy of a hyperelastic material positive definite in general
Is the spatial second derivative of the strain energy of a hyperelastic material positive definite in general?

If this is not a general property of hyperelastic materials are there techniques for approximating positive definiteness? Or are some hyperelastic material models positive definite?

I am doing some large-scale nonlinear optimization on hyperelastic materials and I want to be able to leverage fast linear solvers such as LU factorization and Cholesky decomposition both of which require positive definite systems.

Update

I found a text that appears to suggest that the principle of minimum potential energy corresponds to the requirement that the Hessian of the potential energy is positive definite. Is this principle applicable in this case?


问:
Help with code for FDM Couette flow of two fluids on python
I'm trying to solve numerically with a python code a very simple BVP that represent the velocity field of two fluids between two plates where the driving force of the momentum transport is only for the movement either upper plate or bottom plate.
$V$ is the velocity (not in module) of bottom and $U$ is the velocity of upper plate. $m1$ and $m2$ are the viscosity of fluids. Here the equations:

$$\frac{d^2 v(y)} {dy^2} = 0 $$

with:

$v(0)=V$

$v(\varepsilon H)=u(\varepsilon H)$

And:

$$\frac{d^2 u(y)} {dy^2} = 0 $$

with:

$ v'(\varepsilon H) m1=m2 u'(\varepsilon H) $

$u(H)=U$

To solve the problem I tried to use the centered scheme of FDM at #2 order of accuracy (there isn't for first order acc.), and for the Neumann condition I tried to use for $dv/dy$ the backward scheme and for %du/dy$ the forward scheme.
The problem is the interface, I cant' find a numerical solution because the problem doesn't converge.

This is my python code, i hope that there are an error:

```<code>import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

# Parameters
e1 = 1 #viscosities
e2 = 5
H = 2
h = 0.001
U = 14
V = 0
ep = 0.3 # ratio between thickness of fluid 1 and total thickness

#nodes
N = int(H/h+1)
a = e1/e2
M = int(ep*H/h+1) #interface node

def eq_field(z):

  v = z[0:M+1] #da 0 a M
  u = z[M+1:] #da M+1 fino alla fine
  f = np.zeros(N+1)

  f[0]= -V + v[0]

  for i in range(1,M):
    f[i]=v[i-1]-2*v[i]+v[i+1]

  f[M] = v[-1]-u[0]
  f[M+1] = v[-1]-e2/e1*(u[1]-u[0])-v[-2]
  f[M+2] = e2/e1*(u[1]-u[0])+v[-2]-2*u[1]+u[2]

  for i in range(1,N-M-3):
    f[i+M]=u[i-1]-2*u[i]+u[i+1]

  f[-1]=u[-1]-U

  return f

p1 = np.linspace(0,U/2,M+1)
p2 = np.linspace(U/2,U+1,N-M)
p = np.concatenate((p1,p2))

field = fsolve(eq_field,p)
print(field[int((ep*H)/h+1)])
print(field)
x = ep*H*np.ones(len(field))
y = np.linspace(0,H,len(field))
plt.plot(field,x)
plt.plot(field,y)
plt.show()'
</code>```

I appreciate every tip.

OwnerDisplayName: Mario Riccio


问:
Efficient cutting of mesh edges
I am looking for efficient algorithms to cut a mesh along edges. I have a (half-edge) mesh and a list of inner edges that I want to cut, such that both are new boundary edges. At each vertex there can be an arbitrary number of cuts, either starting at the vertex or crossing the vertex. If you are cutting incrementally, you can also assume that each cut starts at the vertex. I do not want to cut the mesh into several connected components, but only to splice it along a chain of edges, which often will not result in a new connected component.

In my current approach I mark each vertex that belongs to one or more cut edges and then iterate over the adjacent edges and if they are marked as a cut edge I add the edge and the triangles seen before the edge to a list. Then I create a new vertex for each entry in the list, remove the triangles that belong to the entry (i.e. the triangles seen between the cut edge and the previous cut edge) and re-insert them using the new vertex.

I think this approach may be quite inefficient, but I haven't found much literature on how to cut meshes, as most papers that turn up when looking for algorithms explain how to find a good cut, but not how to actually cut a mesh along these edges.

Splitting a mesh into connected components when the cut edges connect two boundary vertices is also rather simple, as you can do a search from a seed and add neighbouring triangles if the edge between the triangles is not a cut edge, and then rebuild the mesh by copying the vertex set and adding only the triangles of the connected component found. However, this algorithm does not work if you want to add cuts that do not create a new connected component.

Are there any known algorithms for splitting mesh edges in an efficient way, i.e., splitting an inner edge into two boundary edges? If it is difficult with a half-edge data structure, I would also consider using other data structures that are more efficient for such operations.

Title: Efficient cutting of mesh edges


问:
How to find volume of a depression in a bone
This is known as an articular pillar. It lies on the outer surface of the bone. We have 3D Slicer openCV, but we are unable to find it as its a little irregular in shape and also its 3D in nature. Kindly help in suggesting how we may find it.

OpenCV couldnt be used because its to calculate volume only of 2D structures. 3D Slicer was used to get this image from the image of whole body's CT Scan.

in this image, you can see that lower right image, the left space's volume is what we wanted. And it looks like a void/cavity/depression in the top image.

Hope the question is understandable


问:
Matrix for Marker and Cell grid?
I have an assignment question that reads:

Show that the combined matrix for the Marker and Cell (MAC) grid
for velocity and pressure for the steady Stokes equations is symmetric.
You can consider the case of a 2D square domain with homogeneous
Dirichlet boundary conditions for the velocity

I am trying to understand what this actually means. I have only used MAC grids in the context of enforcing compresibility conditions for fluid simulations. In that context the MAC grid is more a data structure than it is a mathematical object.

I am trying to understand what the assignment means by the "combined matrix". The prof for my course is on a symposium atm so I can't ask them.

Title: Matrix for Marker and Cell grid?


问:
Need help implementing finite difference Beam Propagation Method to Solve 2-D Helmholtz equation
I am trying to implement beam propagtion method in a two-dimensional lattice to solve Helmholtz equation by following the scheme given this paper. I am using Matlab for implementation.

The expected results are:

For ordered Lattice:

For disordered lattice:

There is localization in the disordered lattice due to Anderson Localization.

Now here's my attempt at reproducing these results:
First, I define lattice parameters and constants.

```<code>    clear all;
    close All;

    %constants---------------------------------------------%
    lambda = 1e-6; % wavelength in micrometer
    k_0 = 2*pi/lambda;
    n_r = 1.46; %constant refractive index

    %------------------------------------------------------%
    N_x = 100; %steps along x

    x1 = -350e-6;
    x2 = 350e-6;
    delta_x = (x2-x1)/(N_x);
    x = linspace(x1,x2-delta_x,N_x); %steps along x

    z1 = 0;
    z2 = 7e-3;
    delta_z = 0.2e-6;
    N_z = int32((z2-z1)/delta_z);%number of steps along z
    z = linspace(z1,z2-delta_z,N_z);
</code>```

Now there since this paper aims to reduce the algebra of finite difference method into simple matrix operations, the $\delta_x^2$ is simply given by the $D_x$ matrix which is tridiagonal and has -2 on its primary diagonal and 1 on both the sub diagonals. In Matlab, I have written

```<code>D_x = diag(-2*ones(1,N_x)) + diag(1*ones(1,N_x-1),1) + diag(1*ones(1,N_x-1),-1);
</code>```

There is an $S$ operator. I am not sure of its functionality other than that it is used to define the $P$ matrix used in the step calculation $\phi(z+\Delta z) = PQP \phi(z)$. The $S$ operator is defined as: $$S = k_0^2 n_r^2 I + \frac{1}{\Delta_2} \left( D_x-\frac{1}{12} D_x^2+\frac{1}{90} D_x^2-\frac{1}{560} D_x^2 \right)$$, which I defined the code as:

```<code>S = (k_0^2)*(n_r^2)*eye(N_x)+(1/(delta_x^2))*(D_x-((1/12)*D_x^2)+((1/90)*D_x^3)-((1/560)*D_x^4));
</code>```

The elements of $P$ operator is defined as:

```<code>P_11 = cos(sqrtm(S)*delta_z/2);
P_12 = (sin(sqrtm(S)*delta_z/2))/sqrtm(S);
P_21 = -sqrtm(S)*sin(sqrtm(S)*delta_z/2)/sqrtm(S);
P_22 = cos(sqrtm(S)*delta_z/2);

</code>```

I am confused about defining the R  matrix part, which corresponds to the lattics structure. In the paper, it is mentioned to specify the refractive index by specifying the field at different nodes $x_j$, at which the refractive index is defined as $n_j^2(z) = n^2(x_j,z)$. Then this array is assigned to the diagonal of $R$ matrix such that $R_{jj}(z) = n^2(x_j,z)$. This matrix is used in defining the $Q$ matrix which is given by $$
\begin{bmatrix}
 I & 0 \\
 -k_0^2 R(z) \Delta_z & I\\
\end{bmatrix}
$$

I define $R$ matrix in my code as follows:

```<code>R = zeros(N_x);
delta_n = 0.003 %difference between the high and low refractive index

flag = 0;
c = 0; % this is strength of randomness in the lattice refractive index
sqwv = zeros(N_x,1); %this matrix is pretty much the <span class="math-container">$x_j$</span> matrix which is later assigned to the diagonal of the R matrix
del_c = 0;
del_mat = [];

if c == 0
    for i = 1:N_x
        del_mat(i) = del_c;
        if flag == 0
            sqwv(i) = 1.46^2-n_r^2+del_c;
        elseif flag == 1
            sqwv(i) = 1.463^2-n_r^2+del_c;

        end
        if rem(i,1) == 0
            del_c = c*delta_n*rand(); %disorderness parameter for the refractive index when I want to generate a disorder in lattice
            if flag == 0
                flag = 1;
            elseif flag == 1
                flag = 0;
            end
        end

    end

    figure();
    plot(x,sqwv);

    for i = 1:N_x
        R(i,i) = sqwv(i);
    end
end
</code>```

For $c = 0$, I get $R$ matrix(lattice) of form:

Now I define my input signal which is a gaussian beam, and it's derivative:

```<code>psi = exp(-(x/W0).^2);
psi = transpose(psi); %transposed the input gaussian beam for consistent matrix multplication later in the code

psi_der = diag(zeros(N_x));

mat_new(1,:) = psi./max(psi); % mat_new matrix is two by two matrix that stores the output. I have assigned the first row to be the input signal
</code>```

Now I run the main loop across $N_z$ i.e. the number of lattice points in z-direction

```<code>for i = 1:N_z-1

    if c~=0
        for i = 1:N_x
            del_mat(i) = del_c;
            if flag == 0
                sqwv(i) = 1.46^2-n_r^2+del_c;
            elseif flag == 1
                sqwv(i) = 1.463^2-n_r^2+del_c;

            end
            if rem(i,1) == 0
                del_c = c*delta_n*rand();
                if flag == 0
                    flag = 1;
                elseif flag == 1
                    flag = 0;
                end
            end

        end
        R = diag(sqwv);

        eta = P_11*psi+P_12*psi_der;
        eta_der = P_21*psi+P_22*psi_der;

        QP_11 = eta*delta_z;
        QP_21 = -(k_0^2)*R*(normalize(eta_der))*delta_z;

        psi_new = P_11*QP_11 + P_21*QP_21;
        psi_der_new = P_21*QP_11+P_22*QP_21;

        psi = psi_new;
        psi_der = psi_der_new;
        mat_new(i+1,:) = psi_new;
        %mat_new(i,:) = psi_new;
    end
    if c==0
        eta = P_11*psi+P_12*psi_der;
        eta_der = P_21*psi+P_22*psi_der;

        QP_11 = eta*delta_z;
        QP_21 = -(k_0^2)*R*(normalize(eta_der))*delta_z;

        psi_new = P_11*QP_11 + P_21*QP_21;
        psi_der_new = P_21*QP_11+P_22*QP_21;

        psi = psi_new;
        psi_der = psi_der_new;
        mat_new(i+1,:) = psi_new;

    end

end
</code>```

Notice that there is conditional for c~=0 and $c=0$. That is done so that when I require a disordered lattice (i.e. c~= 0), $R$ matrix is calculated for every iteration to introduce disorder throughout the system. Otherwise, $R$ matrix calculated above is used. The multplications in the loop are as follows:
$$
P\Phi(z) = \begin{bmatrix}
            P_{11} \psi + P_{12} \frac{\partial \psi}{\partial z} \\
            P_{21} \psi + P_{22} \frac{\partial \psi}{\partial z}
       \end{bmatrix} = \begin{bmatrix}
                            \xi \\
                            \frac{d \xi}{d z}
                        \end{bmatrix}
$$

then,
$$
Q(z)\{P\Phi(z)\{ = \begin{bmatrix}
                   \xi \Delta z \\
                   -k_0^2 R(z) \frac{d\xi}{dz} \Delta z
                \end{bmatrix}
$$
Then, this matrix is multiplied with $P$ matrix again and the next step of propagtion is obtained ($\phi(z+\Delta z) = PQP \phi(z)$). Then I assign it to the mat_new matrix every iteration and then finally mesh plot the mat_new matrix.

```<code>mat_new = (abs(mat_new).^2)./max(abs(mat_new).^2);%to only plot the absolute values
figure();
mesh(mat_new);
xlabel("x-axis");
ylabel("z-axis");
zlabel("Intensity");
colormap default;
colorbar;
</code>```

For c=0 i.e. ordered lattice, I am getting a result which is somewhat okay:

I am not sure why the signal is spread out so much in the beginning (the parallel yellow lines) even though the input signal is a gaussian beam and such a spread is unexpected.

For c not eqaul to zero, there is the expected localization but it is random, and broken. It should be concentrated at the centre, but there extra lines in the lattice besides that.

Please help me figure out what am I doing wrong. It would be helpful to know how to implement the $R$ matrix properly. Or how to generate a disordered lattice.


问:
Numerical method for space fractional derivative in 1 dimension
I am very new to the subject of fractional derivatives which arise while characterizing the anomalous transport of passive scalar in turbulence.

I have found an equation of the following form, to describe such anomalous transport, which I have to solve numerically,
$$\partial_{t}f(x,t)=kD^{\alpha}_{|x|}f(x,t)+g(x,t)$$
where $k$ is some constant, $D^{\alpha}_{|x|}$ is symmetric Riesz fractional derivative of order $\alpha$ and $g(x,t)$ is a source/sink term.
The numerical algorithm that I can find are too mathematically abstract for me,
therefore it would be of great help if anyone could sketch out the basic ideas behind constructing numerical algorithms for such equations.


问:
About Convolution Theorem
```<code>import numpy as np
import scipy.signal as ss
from scipy.fft import fft, ifft
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
np.convolve(x,y) --> gave me result as array([ 4, 13, 28, 27, 18])
</code>```

Now, I wanted to get the same answer using FFT and the code is below.

```<code>a = fft(x)
b = fft(y)
d = ifft(x*y)
d --> gives array([31.+0.j, 31.+0.j, 28.-0.j])
</code>```

Would appreciate if anyone could provide the right code or let me know why it is so?

Title: About Convolution Theorem


问:
Adding stability to MPM simulation?
I am writing a 2D implementation of MLS-MPM, I have fluids working perfetly fine, solids technically work as well, at low time steps.

This is the fluid simulation at a large time step:
https://imgur.com/a/fEuzN9z

And this is a pciture of the solids at a low time step after a large amount of simulation:

I am running into the following problem, the MLS-MPM method uses a deformation matrix for particles who are meant to remain solid. At large time steps the inaccuracies cause this matrix to become degenerate and my simulation explodes.

This is a rust implementation of the critical part where the problem happens:

let volume = particle.mass / density;

            let j: f32 = properties.deformation_gradient.determinant();

            assert!(j > 0.0, "{}", properties.deformation_gradient);
            let volume_scaled = volume * j;

            let f_t = properties.deformation_gradient.transpose();
            let f_inv_t = f_t.try_inverse().unwrap_or(Mat2::identity());
            let f_minus_f_inv_t = properties.deformation_gradient - f_inv_t;

            let p_term_0 = properties.elastic_mu * f_minus_f_inv_t;
            let p_term_1 = properties.elastic_lambda * j.ln() * f_inv_t;
            let p_combined = p_term_0 + p_term_1;

            let stress = (1.0 / j) * p_combined * f_t;
            let eq_16_term_0 = -volume_scaled * 4.0 * stress * self.simulation_speed;
```

Which I have taken from here who got it from this tutorial who modified it from the original paper's author.

The short of it is, it seems that large timesteps make the deformation matrix flip from a positive determinant to a negative determinant. SInce this code uses logarithms, that negative determinant is a no-no, as it destroys everything.

I was wondering if there is a way for me to get numerical stability at larger time steps. I am fine sacrificing accuracy as this is for a toy implementation.

Title: Adding stability to MPM simulation?


问:
When is it worth it to use the dual simplex method with warm starts?
Let's say I have a linear program that incurs a series of slight changes to it, so I want to warm-start it. I've read various things that recommend using the dual simplex algorithm over the primal depending on the nature of the change, e.g.

If the objective coefficients are changed or columns are added, warm start the primal.

If the RHS of the constraints is changed or a constraint is added, warm start the dual.

These guidelines essentially recommend solving the LP that maintains feasibility. I had a couple of concerns though:

Warm starts usually imply supplying a basis — will the corresponding dual basis simply be the variables not included in the primal basis?

Let's say the changes being made are mostly constraint RHSs changing. The criteria above would suggest to use the dual. Is this still the case if the number of constraints is small compared to the number of variables?


问:
Where to find datasets to research the optimal parameter search space for HPC systems?
I am a student researcher who is new to the HPC domain.

I have recently taken a project where I am working on optimizing the parameter search space (both application level and hardware level parameters) for running HPC applications. In order to do this, I am trying to implement the Multi-arm bandit approach to find the optimal parameters by taking each parameter as an arm of the said algorithm and the performance metric(Accuracy, energy, execution time, etc.) as the reward for that particular arm.

But, I am facing difficulty in how implementing this approach experimentally. I was thinking about getting a dataset that can map a parameter (or a certain group of parameters) to the said performance metric. I have seen some previous works which deal with this kind of problem finding the best optimal parameters (Example: https://www.cs.umd.edu/~bhatele/pubs/pdf/2018/ics2018.pdf). But, these works do not share the dataset that they worked with. Does anyone know where I can get access to the dataset which is befitting to my particular problem?


问:
Accessing nodal and degree of freedom organisation in the solution vector in Fenics and Firedrake
Context

Let suppose a solution to a finite element problem (associated to a mesh) stored in a vector $X$  of size $(N \cdot d)$ corresponding to $N$ nodes in the mesh and $d$ degrees of freedom for each node.

This vector could be organised this way :

$X = (X_1^1, \ldots, X_1^d, X_2^1, \ldots, X_2^d,\mathbf{\cdots}, X_N^1,\ldots,X_N^d)$

i.e for every node, we concatenate every dofs evolution, then we iterate for every node...

However, as Fenics or Firedrake optimise their computation cost using sparse mass and stifness matrices, the final organisation of the solution might not be structured as simply as it is in this example.

Needs

For a project, I need to isolate the evolution of the first dof of each nodes (in my case, it will permite plotting the temporal displacements of a plate).

Question

How do you access the organisation of nodes and dofs in the solution vector, so I can isolate only the first dof evolutions, with Fenics and/or Firedrake ?


问:
Gradient based optimization involving a black box function where gradients (not the objective) are approximated using a surrogate model
Can we use gradients of surrogates of a black box function and the actual function evaluation for optimization involving an expensive black box simulation? Is there any merits to it?</p>


问:
What is meant exactly by the domains in this context and how can I simulate their size?
I was working on a problem from "Computational Physics using python: Chapter 17." The problem is concerned with simulating the thermodynamics behind the spin of electrons and the resulting change in energy and magnetization. The last part of the problem asks for making a graph for the average domain size versus temperature. I'm familiar with magnetic domains, but I'm not really sure how to simulate them.

My initial thinking was domains here is similar to measuring how many neighbouring electrons have the same spin. So, I though about writing a code that checks every electron in the configuration and counts the average size of domains by checking whether they are similar or not. I tried it but didn't work out as I want. Does anyone have any suggestions?

My 1D Ising model:

```<code>%matplotlib inline
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt

# Defining the initial spin of all N particles of the system
def initialstate(N):
    state = 2*np.random.randint(2, size=(1,N))-1
    return state

# Defining the monte carlo random flipping of the spins
def mcmove(config, beta):
    for i in range(N):
        b = np.random.randint(0, N) # what spin we try to flip (a random spin)
        s =  config[0,b] # initial spin of the spin we will try to flip
        nb = config[0,(b+1)%N] + config[0,(b-1)%N]
        cost = 2*s*nb
        if cost < 0:
            s *= -1 # flip the spin
        elif rand() < np.exp(-cost*beta):
            s *= -1 # flip the spin
        #otherwise we don't flip the spin
        config[0,b] = s # new state (with s either flipped or not)
    return config

# Calculating the energy of the system
def calcEnergy(config):
    energy = 0
    for i in range(N):
        S = config[0,i]
        nb = config[0,(i+1)%N] + config[0,(i-1)%N]
        #print nb
        energy += -nb*S
    return energy/2.

# Calculating the magnetization of the system
def calcMag(config):
    mag = np.sum(np.sum(config))
    return mag

nt      = 2**6        # number of temperature points
N       = 2**4        # size of the lattice, N
eqSteps = 2**10       # number of MC sweeps for equilibration
mcSteps = 2**10       # number of MC sweeps for calculation

n1, n2  = 1.0/(mcSteps*N*1), 1.0/(mcSteps*mcSteps*N*1)

T=np.linspace(.2, 2, nt)

Energy       = np.zeros(nt);   Magnetization  = np.zeros(nt)
SpecificHeat = np.zeros(nt);   Susceptibility = np.zeros(nt)

for m in range(len(T)):

    config = initialstate(N)
    iT=1.0/T[m]; iT2=iT*iT; # beta = 1/(k_B*T)

    for l in range(eqSteps):         # equilibrate
        mcmove(config, iT)           # Monte Carlo moves

    #Calculate quantities post equilibriation
    E1 = M1 = E2 = M2 = 0
    for i in range(mcSteps):
        mcmove(config, iT)
        Ene = calcEnergy(config)     # calculate the energy
        Mag = calcMag(config)        # calculate the magnetisation

        E1 = E1 + Ene                 #Cumulative energy
        M1 = M1 + Mag
        M2 = M2 + Mag*Mag
        E2 = E2 + Ene*Ene

        Energy[m]         = n1*E1     #Weighing for MC averages.
        Magnetization[m]  = n1*M1
        SpecificHeat[m]   = (n1*E2 - n2*E1*E1)*iT2
        Susceptibility[m] = (n1*M2 - n2*M1*M1)*iT2

f = plt.figure(figsize=(15, 8)); # plot the calculated values

sp =  f.add_subplot(1, 2, 1 );
plt.plot(T, Energy, 'o', color="#348ABD", label='simulation');
plt.plot(T, -np.tanh(1.0/T), color="#A60628", lw=2,label='analytical')
plt.xlabel("Temperature (T)", fontsize=20);
plt.ylabel("Energy / spin ", fontsize=20); plt.legend(fontsize=16)

sp =  f.add_subplot(1, 2, 2 );
plt.plot(T, Magnetization, 'o', color="#348ABD", label='simulation');
#plt.plot(T, (1.0/T**2)*(np.cosh(1.0/T))**(-2), color="#A60628", lw=2,label='analytical')
plt.xlabel("Temperature (T)", fontsize=20); #plt.ylim([-10, -2])
plt.ylabel("Magnetization / spin ", fontsize=20);  plt.legend(fontsize=16)

</code>```

A snippet of the code I tried but always gave me a size of zero:

```<code>domain = 0
for i in range(len(T)):
  for i in range(N):
     s  = config[0, i]
     if s == s+1:
       domain += 1
</code>```


问:
Encountering blow-up when solving the one-way heat equation using Lax-Wendroff
This is my first time attempting to implement a finite difference method for a PDE in Python, and I am having a bit of trouble. The PDE I am trying to solve is as follows:
$$
 \begin{cases}
        v_t - 3v_x = 0, & (x,t) \in [0,1] \times [0,1] \\
        v(x,0) = 3 - 6x, & 0 \leq x \leq 1    \end{cases}
 $$

The analytical solution is $v(x,t) = 3 - 6x - 18t$.  I am attempting to use the so-called Lax-Wendroff method:

$$
     u_k^{n+1} = u_k^{n} - \frac{R}{2}(u_{k+1}^n - u_{k-1}^n) + \frac{R^2}{2}(u_{k+1}^n - 2u_k + u_{k-1}^n), \quad k=1,2,\ldots,M-1,
$$

where $R := -3 \frac{\Delta t}{\Delta x}$ and $M := 1/\Delta x$  For the boundary conditions, we have
$$ u_k^0 = 3 - 6k \Delta x, \quad k = 0,1,2,\ldots,M $$

As for the spatial boundary values, I tried the following one-sided difference approximations:
\begin{align*}
    u_0^n = u_0^{n-1} + \frac{u_0^{n-1} - u_0^{n-2}}{\Delta t}, \qquad u_M^n = u_M^{n-1} + \frac{u_M^{n-1} - u_M^{n-2}}{\Delta t}, \quad n = 0,1,\ldots,N
\end{align*}

where $N := 1/\Delta t$. Unfortunately, my solution quickly starts to blow up. I suspect this is has to do with my handling of the boundary conditions... Here is my code:

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

#In this program we solve v_t + a*v_x = 0

a = -3
dx = 0.01   #spatial discretization
dt = 0.01   #time discretization

M = int(1/dx)   #Number of x sub-intervals; k = 0,1,2,...,M
N = int(1/dt)   #Number of t sub-intervals; n = 0,1,2,...,N

R = (a*dt)/(dx)

#The x coordinates of the gridpoints
x_pts = np.linspace(0, 1, M + 1, endpoint = True)

#Define the initial condition: v(x,0) = f(x)
def v_0(x):
    return 3-6*x

u = np.zeros((M+1,N+1))  #Create M+1 x N+1 array to store the solutions

#Populate the solution array with the initial condition values
u[:,0] = v_0(x_pts)

#%%

#Each *column* of the array is the solution at time step n

for n in range(1,N+1):     #n = 1,2,...,N

    #Handle the boundary (x = 0)
    u[0,n] = u[0,n-1] + (u[0,n-1] - u[0,n-2])/dt

    #Handle the boundary (x = 1)
    u[M,n] = u[M,n-1] + (u[M,n-1] - u[M,n-2])/dt

    #Lax-Wendroff scheme:
    for k in range(1,M): #k = 1,2,...,M-1
        u[k,n] = u[k,n-1] - 0.5*R*(u[k+1,n-1] - u[k-1,n-1]) + \
                 +  0.5*(R**2)*(u[k+1,n-1] - 2*u[k,n-1] + u[k-1,n-1])
</code>```

The plots at times $t = 0, 0.01, 0.02$ show a very rapid blow up that starts at the boundaries:

What can I do to get this scheme to work? I have tried altering the discretization parameters, but the problem doesn't go away. In theory this method should be unconditionally stable, at least according to section 5.3 of JW Thomas' numerical PDEs book...Any help or suggestions would be greatly appreciated.

Update 4/16/2023:

As @Lutz Lehman indicated in the comments, the difference quotient (u[0,n-1] - u[0,n-2])/dt is not defined when $n = 1$; by default Python used a value on the order of $10^{100}$ in place of $u[0,-1]$, which would explain the blow up! The following changes appear to have solved the problem:

I handled the boundary conditions as follows:

```<code> #At x = 0, use *forward* difference approx for u_x
    u[0,n] = u[0,n-1] - R*(u[1,n-1] - u[0,n-1])

    #At x = 1, use *backward* difference approx for u_x
    u[M,n] = u[M,n-1] - R*(u[M,n-1] - u[M-1,n-1])
</code>```

I realized that according to the stability theory in Thomas, we require $|R| \leq 1$ for the method to be stable. With my initial choice of discretization parameters I had $|R| = 3$, so I set dx = 0.05 and
dt = 0.01 so that I now have $|R| = 0.6$.

The solution plots now appear correct!


问:
Resource to learn assembly code
I'm a PhD student in mechanical engineering and I have to perform a lot of simulations for my project. In my lab we use several well-known libraries, from FEM to machine learning. As we're doing performance critical simulations, I've heard the advisor saying that being able to read assembly code is certainly a very promising way to see what the compiler is doing and hence optimise some relevant parts of a code.

Of course, I will never need to write directly assembly instructions, but during my PhD I think it'd be nice if I can read assembly code with GDB or objdump like other computational scientists in my group are doing every now and then.

So, coming to my question: is the some good reference to recommend, especially targeted to people that just need to understand it?

Title: Resource to learn assembly code


问:
Refluxing step on Finite difference AMR
Hi I am a computer scientist working on MHD code for astrophysics simulation. We use a finite difference scheme where we first solve the spatial derivatives and with them solve the right hand side and take timestep forward.

Now I am interested in looking into adding AMR into our code. However most AMR material concerns Finite Volume Method or at least it seems so. In most cases this is not a problem but I am having difficulties understanding how to implement refluxing step in Finite Difference case, or is that even possible.
To give a concrete example let's say we want to solve simple equation:
u_t = f_x, in 1d where u_t means time derivative and f_x spatial derivative.
Now f_x can be calculated with f_x(x) = (f(x+h)-f(x-h))/2h.

Okay we can solve all of these equations on the coarse and fine level, using coarse level values as boundary value problems for the fine level and then averaging the fine level values and overwriting the coarse level with them where the fine and coarse level overlap.

However there usually is presented the refluxing step where one should calculate the flux on the coarse/fine-boundary with the coarse values and fine values. Then the difference of the fine flux and the coarse flux should be spread evenly on the boundary to ensure conservation.

The problem I am having that what is supposed to be the flux this case? In cell-centered Finite Volume methods the flux is easier to implicitly see on the faces of the control volumes. Does refluxing in finite difference context mean that we have to recalculate the f as f(x) =(f_x(x+h0.5)-f_x(x-h0.5))/h where h is the grid distance on the coarse level and we have the finer f_x(x+h0.5) from the finer level.

Again not really sure what would be the corresponding operation to refluxing in this case. In case someone understood my not best formulated question all answers are appreaciated :)

Title: Refluxing step on Finite difference AMR


问:
How can i study stability for a new method that solves second degree non lineara differential equations?
I developed a new method to solve this equations $\frac{d{y}^{2}}{dx^{2}}=g(x,y,\frac{dy}{dx})$ for a general g(x,y,z), wich their solution is a function f(x)=y. Obviously you need $y(x_0)=y_0$ and $y'(x_0)=y_1$.

The thing is i do not have a clue how to study stability of solutions, or to see if converges rapidly to the solution, or it's range of aplication. Can someone orientate me? the method first aproximates the function, and then iterates with the data available, wich is the second derivative the first and the value of the function, i don't use taylor aproximations. I programmed it and works perfectly fine.


问:
A question from boundedness property of LMT
I have the following PDE
$$
u_t = k \Delta u + \alpha u H(u-c)
$$
I am trying to show the boundedness property to apply LMT. I get confused with the estimates of third term.

The space I have
$$V = \{v\in H^1(\Omega) : v(0,y)=v(1,y)~ and~ \frac{\partial v}{\partial y}(x,0)=0=\frac{\partial v}{\partial y}(x,1)\}$$
and the function $H$ is a linear function and it is $0 \leq H(x) \leq 1$.

The third term in the bilinear form is

$$
|(\alpha u^{n+1} H(u^n-q_c),w)_{L^2(\Omega)}| \leq \alpha \|u^{n+1}\|_{L^2(\Omega)} \|w\|_{L^2(\Omega)} \leq \alpha \|u^{n+1}\|_{H^1(\Omega)} \|w\|_{H^1(\Omega)}
$$

Is third estimate correct here? Can I split the $L^2$ norm of $u^{n+1}H(u^n-q_c)$ integrands?


问:
General dimensional solution to fast sweeping quadratic equation
I am reading Hongkai Zhao's paper The Fast Sweeping method. I have implemented the method in 2D and now I want to move onto 3D.

However section 2.6 confuses me.

The article says:

But I have no idea how that is actually supposed to work. i.e. On a purely computational level, what steps do I actually follow to solve that quadratic equation?


问:
Compatibility condition for Poisson equation in cylindrical symmetry
I'm trying to implement multigrid approach for a Poisson equation $\frac{1}{r}\frac{\partial}{\partial r}\left( r \frac{\partial H}{\partial r} \right) = f$ with all Neumann boundary conditions. Jacobi iteration with underrelaxation is used for smoothening ($\omega = 2/3$ as in Briggs Multigrid tutorial). After few iterations on the coarsest grid, residual ($R_i = \frac{1}{r_i}\frac{1}{\Delta r}\left( r_{i+1/2} \frac{H_{i+1} - H_{i}}{\Delta r}  - r_{i-1/2} \frac{H_{i} - H_{i-1}}{\Delta r} \right) - f_i$) and the residual norm ($\sum |R_i|$) converges to some value, which is not zero. In the case of cartesian coordinates I used discrete compatibility condition $\sum f \Delta x = 0$ to ensure solution uniqueness in the case of all Neumann conditions. Should I modify this condition some how to ensure that the proсedure converges to the solution with zero residual?


问:
Rank-one updates for symmetric matrix eigen-system
Are there existing implementations for rank-one updating of symmetric matrices eigensystems?
This is the mathematical statement of the problem. Let $S=QDQ^T$
$$S + vv^T = QDQ^T +vv^T = Q_{new}D_{new}Q_{new}^T$$
This paper describes the algorithm for efficiently computing update, but the closest thing I found was SLAED7 in LAPACK, which assumes S is tridiagonal symmetric.


问:
How best to iterate over three dimensional grid?
My problem is to iterate through $[0.95,1] \times [0.95,1] \times[0.95,1] \in \mathbb{R}^{3}$. Is there a generic way to do this? I was considering rolling it all out and effectively vectorising the problem.

I would like to use Eigen but I think using a vector of MatrixXd elements seems not wise.

Any tips or help would be appreciated!


问:
Estimating/Tuning a Coefficient from a Quadratic Programming Objective Function so Optimal Solution Reflects Data
I am working on a problem that is a modified version of a two-knapsack knapsack problem. I am able to find the optimal choices by using Gurobi. However, I would like to estimate a coefficient that is added linearly to my objective function so that the results most closely match my data. I will describe the problem below, but the thoughts I have had so far:

Brute force/use grid search to try several parameter values by reoptimizing each time. This, however, would be computationally too expensive to be feasible.

Use Gurobi's multiple scenarios option. However, this suffers from a similar problem, as the coefficient can take any real value (although it realistically will be between roughly [-10,10]). This option also feels too computationally expensive/would lack precision unless I could try different scenarios based on the results of previous scenarios (ie to focus on trying values closest to the results with the best outcomes.

Setting up a bilevel optimization problem by minimizing the sum of squared errors, where the values come from the first-level optimization problem. Then use a Kuhn-Tucker/Lagrangean transformation. However, I have been having trouble with this approach.

Any other options that people can think of would be incredibly helpful!

The problem involves many schools needing to be placed into one of two knapsacks, $X$ or $Y$. By being placed in either knapsack, they will get value from data, $w_{i}$, and will get no value by being out of the knapsacks. Knapsack $Y$ must have an average weight (average of the w's of the schools in the knapsack) above 0.7, and knapsack $X$ must have an average weight below 0.7. Additionally, the value a school gives from being in knapsack $Y$ is 0.7. However, the value from being in knapsack $X$ is the average weight (again, the average of all $w_{i}$s in the knapsack). This problem I can solve in Gurobi.

However, I want to assume that a school has a value from being in either knapsack, $\theta$. I am interested in the minimum $\theta$ required for it to be optimal for a school to be placed in a knapsack (or maximum such that it won't be in a knapsack). I want to estimate theta such that my model's prediction if a school is in any knapsack most closely matches the reality from the data. I would be satisfied estimating a $\theta$ for each school or one shared $\theta$ for the entire problem.

Below is my formulation of the single-level and bilevel problems. $A$ is the average weight of the $X$ knapsack, $X_{i}$ and $Y_{i}$ represent being placed in a knapsack, and for each school, their sum has to be less than or equal to 1 (can only go in at most one knapsack). $P_{i}$ represents "participating" or being placed into either knapsack. $P_{i}$ is observed in my data, and I am trying to minimize the squared error from comparing my model's prediction of going into a knapsack ($X_{i} + Y_{i}$) to the data, $P_{i}$. As $P_{i}$, $X_{i}$, and $Y_{i}$ are binary, the first equation in the bilevel problem represents the sum of squared errors.

Does anyone have any thoughts on how I could go about estimating theta to match my data? Either a method to reformulate the problem, a trick Gurobi (or similar software) can do, tips for the Kuhn-Tucker reformulation, or something different would be greatly appreciated.

Note, this is a simplified version of my problem, but it contains the core details. I will also include screenshots from the more complicated version of this problem if you are curious/able to provide support on the real version.

Original Optimization Problem:
\begin{equation}
    \max_{X_{i},Y_{i}} \sum_{i=1}^{n} A \cdot X_{i} + 0.7 \cdot Y_{i} + (X_{i} + Y_{i})\theta
\end{equation}
such that:
\begin{equation}
    A = \frac{\sum_{j=1}^{n} w_{j}X_{j}}{\sum_{j=1}^{n} X_{j}} \leq 0.7
\end{equation}
\begin{equation}
    \frac{\sum_{j=1}^{n} w_{j}Y_{j}}{\sum_{j=1}^{n} Y_{j}} \geq 0.7
\end{equation}
\begin{equation}
    \forall i\quad X_{i} + Y_{i} \leq 1
\end{equation}
$\forall i \quad X_{i}$ and $Y_{i}$ are Binary Variables, $\theta \in \mathbb{R}$ and $\forall j \quad w_{j} \in [0, 1]$

Bilevel Problem for Estimating $\theta$:

\begin{equation}
     \min_{\theta} \sum_{i=1}^{n} P_{i}(1-X_{i}-Y_{i}) +  (1-P_{i})(X_{i}+Y_{i})
\end{equation}
such that:
\begin{equation}
   [X_{i},Y_{i}]_{i=1}^{n} = arg\max_{X_{i},Y_{i}} \sum_{i=1}^{n} A \cdot X_{i} + 0.7 \cdot Y_{i} + (X_{i} + Y_{i})\theta
\end{equation}
\begin{equation}
    A = \frac{\sum_{j=1}^{n} w_{j}X_{j}}{\sum_{j=1}^{n} X_{j}} \leq 0.7
\end{equation}
\begin{equation}
    \frac{\sum_{j=1}^{n} w_{j}Y_{j}}{\sum_{j=1}^{n} Y_{j}} \geq 0.7
\end{equation}
\begin{equation}
    \forall i\quad X_{i} + Y_{i} \leq 1
\end{equation}

Real Bilevel Problem

My attempt at a KKT Reformulation, to be used in Gurobi (with the M and $\pi$ terms for dealing with the bilinearity of the complementary slackness conditions


问:
MPM implicit integration

Setup

I am having a hard time translating the formulas in the 2016 MPM course from siggraph into algorithmic computations.

In particular we have equations 200 and 201

$$h(v^{n+1}) = Mv^{n+1} - \Delta t f(x^n + \Delta t v^{n+1}) - Mv^n = 0$$
$$v^{i+1} = v^i - \bigg(\frac{\partial h}{\partial v}(v^i) \bigg)^{-1}h(v^i)$$

I understand the first is implicit Euler and the second is newton's method to find a suitable approximation for the next stages velocity.

$M$ is a diagonal matrix containing the mass of grid cell $i$, $v^{k}$ are the velocities at the center of each grid cell. $h$ is the conservation of momentum equation.

General question

What I want to do is translate these two paragraphs into pseudocode of what I must tell a computer to do. i.e. I know I need to make a sparse linear system of equations which I must then solve and the solution of that sparse system will be the new velocities at grid cell centers for the next stage. But that is much too general still, I am not clear as to what specific values should be inserted into the sparse matrix.

Specific questions

What is $\frac{\partial h}{\partial v}$? i.e. do I take derivatives with respect to each instance of v regardless of the superscript? Or do I only take it with respect to $v^{n+1}$. Notation here isn;t clear to me since I cannot symbol match directly.

What is $x^n$? It would seem to me that $x$ here represents position, but grid cell positions don't change should that not be a constant?

Can you avoid solving the full linear system directly and instead solve for things locally? Like not make the full matrix but rather query your cell neighbourhood and solve for your next velocity purely from your immediate environment.

Title: MPM implicit integration


问:
Solving a minimization problem without flattening inputs?
I am reading this paper on improving time steps for solving simulation problems:
https://www.math.ucla.edu/~jteran/papers/GSSJT15.pdf

The authors developed this energy  function:

$$E(x) = \frac{1}{2\Delta t ^ 2}(x-\hat x)^TM(x - \hat x) +  \phi$$

In this particular formulation, $x$ is the global, flattened, vector of particle positions, $M$ a lumped mass amtrix, $\hat x$ the global flattened vector of positions in the current time step and $\phi$ the flattened forces acting on all points.

The author's then solve for the minimizer by doing a hessian weighted gradient descent on this energy function.

This works, but has 2 things I don't like.

You have to flatten your vectors from each particle into a large global state vector,

and this does not parallelize too well.

Is it possible to retain the main contribution of the paper (large time steps)  but solve localized problems for each input particle instead of a global equation?

It really doesn't seem to me like you need global state for this as particles which are physically far away from each other should not impact the solution.


问:
Absence of discontinuity in Specific Heat in liquid-gas transition (Based on the Ising Model)
I'm trying to do a model for the transition liquid-gas based on the Ising model and the metropolis algorithm, instead of using values of spins, I'm saying that a cell is occupied by a particle or not. My specific heat doesn't show a discontinuity and I'm struggling to find the problem. I have tried to increase the size of my lattice, the numbers of sweeps and also to decrease the steps that I do to change temperature.The Montecarlo function works as followed :

1.We choose randomly a cell, if it's occupied we choose one of its neighbors

2.If the neighbor is occupied we do nothing

3.If it's not we calculate the energy of the lattice and then we made the particle jump to the second cell

4.We then calculate the energy in this configuration

5.We do the metropolis algorithm as for the Ising model : if conditions are not met we reverse back to the 1st configuration and if they are we add dE to the energy.

To calculate the energy I'm using the same expression as for the Ising Model, which means that if a particles as its 4 nearest neighbors occupied the energy of this site will be 4. And we are going to count like this the number of nearest neighbors site occupied for every site.

To calculate the heat capacity I'm trying to use this equation :
$C=\big[\langle E^2\rangle-\langle E\rangle^2\big]/k_BT^2$
I have also tried to zoom in around 1~3 and the heat capacity only decreases and show 0 discontinuity

Here are my plot for the energy and heat capacity vs temperature with parameter that are in the code :

Here is my code :

```<code>#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define L 40            /*Lattice size*/
#define STEPS 10000     /*Numbers of step to obtain equilibrium*/
#define mol 1000        /*Numbers of particles in the system*/
#define J 1             /*Value of the interaction*/
#define T_MIN 0.1       /*Minimal temperature*/
#define T_MAX 7         /*Maximal temperature*/
#define T_STEP 0.08
#define kB 1

double energy;
int N[L][L];

void display();
void init_params();
double calc_energy();
void MonteCarlo(double temperature);
void save();
 /*This is the main function that will do the calculation for the energy and heat_capacity vs temperature*/
int main() {
    FILE *file1, *file2 ;
    file1 = fopen("energy.dat","w");
    file2 = fopen("heat_capacity.dat","w");

    double temperature;
    for (temperature = T_MIN; temperature <= T_MAX; temperature += T_STEP) {
/*at each temperature we initialize the lattice*/
        init_params();
/*calculate the energy and energy squared and initialize the heat capacity*/
        energy = calc_energy();
        double en = energy;
        double energy_squared = energy * energy;
        double heat_capacity = 0.0;
/* here we are going to do a lot of steps of the function Monte Carlo to obtain equilibrium and to be able to calculate a mean of energy*/
        for (int step = 0; step < STEPS; step ++){
            MonteCarlo(temperature);
            en += energy;
            energy_squared += energy * energy;

        }
        en /=(STEPS*L*L);
        energy_squared /= (STEPS*L*L) ;
        heat_capacity = (energy_squared - en*en)/(kB * temperature*temperature);
/*Saving up the quantities into file to graph them*/
        fprintf(file1,"%lf %lf \n",temperature,en);
        fprintf(file2,"%lf %lf \n",temperature,heat_capacity);

    }
}
/*Initialization of the lattice*/
void init_params() {
    time_t *tm=NULL;
    srand48(time(tm));
    int k=0;
    int a,b;
/*First of all we made it completely empty*/
    for(int i = 0; i < L; i++) {
        for(int j = 0; j < L; j++) {
            N[i][j] = 0;
        }
    }
/* and then we randomly fill it with the number of particles choose before*/
    while (k<mol) {
        a = floor(drand48()*L);
        b = floor(drand48()*L);
        if(N[a][b]==0) {
            N[a][b]= 1;
            k++;
        }
     }
}
/* This is the function that calculate the energy of the lattice with cyclic boundaries, its the same expression as the Ising Model*/
double calc_energy() {
    double energy = 0.0;
    int i,j;
    for(i = 0; i < L; i++) {
        for(j = 0; j < L; j++) {
            int up = (i == 0 ? L-1 : i-1);
            int down = (i == L-1 ? 0 : i+1);
            int left = (j == 0 ? L-1 : j-1);
            int right = (j == L-1 ? 0 : j+1);
            energy += -J*N[i][j]*(N[up][j]+N[down][j]+N[i][left]+N[i][right]);
        }
    }

    return energy/2.0;
}
/*This is the function that will do the MonteCarlo and Metropolis algorithm*/
void MonteCarlo(double temperature) {
    int n,ix,iy,N_init;
    double E1, E2, dE;
    int a,b;
    double p,r;
  /*We are choosing one random site on the lattice*/
    ix = floor(drand48()*L);
    iy = floor(drand48()*L);
    N_init = N[ix][iy];
/*If it's occupied we can do the rest of the code*/
    if (N_init == 1) {
/*Here we are choosing one of its neighbors randomly*/
        n = floor(drand48()*4);
        if (n == 0) {
            a = (ix + L - 1) % L;
            b = iy;
        }
        else if (n == 1) {
            a = ix;
            b = (iy + 1) % L;
        }
        else if (n == 2) {
            a = (ix + 1) % L;
            b = iy;
        }
        else {
            a= ix;
            b = (iy + L - 1) % L;
        }
/*Calculation of the energy of the lattice before a swap*/
        E1 = calc_energy();
/*If the neighbor site is empty we are going to jump the particle from the first site to it*/
        if (N[a][b] == 0){
            N[ix][iy] = 0;
            N[a][b] = 1;
/*We calculate the energy of the new lattice after the jump and make the metropolis algorithm like for the ising model.*/
            E2 = calc_energy();
            dE = E2 - E1;
            p = exp(-dE/(kB*temperature));
            r = drand48();
/*If the condition are not met we return to the first lattice and the energy of the system is E1*/
            if ((dE>0) && (p<r)){
                        N[ix][iy] = 1;
                        N[a][b] = 0;
                        energy = E1;
            }
/*If conditions are met we keep the lattice as it is and have the energy E2 for our system*/
            else {energy = E2;}

        }

    }
}
</code>```


问:
RFEM: Equilibrium over plate section
I'm calculating a simple plate arrangement in RFEM:

The plate has dimensions of 1m X 1m, thickness of 50mm. Elastic modulus of the material is 205000.0 N/mm2, Poisson's ratio is 0.3.

The plate is loaded in the middle, with 10kN load. Self-weight has been neglected. The plate is simply supported at 4 corners:

According to RFEM, the support reactions at the corners is 2.5kN, which at the 4 corners adds to 10kN as expected.

I took a section cut 200mm from one edge (shown in the first picture) and studied the results.

The internal moment m_y looks like this:

The section average is reported as 1kNm/m, meaning that the total moment across the 1m section is 1kNm. A simple check calculation can be done: since the distance from this section to the edge is 200mm, and we have total of 5kN at the two corners, we have 5kN * 0.2m = 1kNm. This matches the section diagram result and the section is at equilibrium.

The shear force in z-direction looks like this:

The total shear is reported as 3.805kN. This does not seem to make the section at equilibrium with the supporting forces. Should the total section shear be equal to the reaction forces, which are a total of 5kN?

Would anyone know if I am misunderstanding something about plate theory in general, or is my mistake something in the software itself?

According to RFEM documentation, the shear force V_y is calculated by integrating the shear stresses in y-face of an element, into z-direction, so the internal force in question should be correct.

Title: RFEM: Equilibrium over plate section


问:
What is the current state of preconditioning highly heterogeneous equations? Where do I start?
Consider the simple standard Laplacian in 2D or 3D:
$$\nabla (\alpha \nabla u) =f $$

with $\alpha$ being a scalar. $\alpha$ can take values that can vary largely between 1e-8 and 1 throughout the domain. After discritizing the equation above with FEM one arrives at a linear system which will encounter convergence difficulties if a standard iterative solver is used.

Problems like these arise an highly heterogeneous media and require special treatment for the iterative solver. What are the current techniques out there? How do they compare in terms of performance, scalability, and ease of application? Where can I start?


问:
Energy potential involvign pressure and viscoscity?
Making a really long explanation short. I am looking for an energy density function for an incompressible isotropic fluid that involves both viscoscity and pressure.

The reason is, there is a way to take really large time steps when simulating something using MPM, but it only works if one has an energy density function.

For example, for a neo hookean constitutive model the energy density function can be found in equation 46 of the MPM 2016 siggraph course:

$$\Psi(F) = \frac{\mu}{2}\bigg(tr(F^TF - d)\bigg) - \mu \log(J) + \frac{\lambda}{2}\log^2(J)$$.

In the same notes a stress tensor is derived through the use of equations 48 and 38.

In these notes the stress is computed as follows:

```<code>// end goal, constitutive equation for isotropic fluid:
// stress = -pressure * I + viscosity * (velocity_gradient + velocity_gradient_transposed)

float2x2 stress = math.float2x2(
    -pressure, 0,
    0, -pressure
);

// velocity gradient - MLS-MPM eq. 17, where derivative of quadratic polynomial is linear
float2x2 dudv = p.C;

// build strain from the velocity gradient
float2x2 strain = dudv;
strain[0][1] = strain[1][0] = trace(dudv);

float2x2 viscosity_term = dynamic_viscosity * strain;
stress += viscosity_term;
</code>```

Which comes from the wikipedia description of the stress tensor for a newtonian fluid.

This is fine for the forward Euler step the author of the notes is using in his demos, but I need an energy density function. GIven that for the neo-hookean model there is an energy density function associated with the stress tensor, I was hoping someone could point me to an energy density function for a Newtonian Fluid.


问:
How to model data that is repeated over time within the same group, from which outcome variable is categorical and explanatory variables are linked?
I have 1 dependent variable that is categorical. And I have 3 explanatory variables, the two continuous variables are likely interrelated and the 1 categorical variable is not interlinked.
The sample group has been measured at two different times and has mostly similar participants.
Now I want to find out how a change in dependent variable can be explained by changes in explanatory variables.
What is the best statistical model to use?

I would like to make a graph out of it, in which the y-axis = categorical outcome, x-axis = percentual change of explanatory variables (-100% till 100%) and the lines are based on two dots from the average explanatory variables per year. I was thinking about using a GEE model. Would this be a good model?


问:
What architecture of CNN to solving a image regression problem (case study: solving Poisson equation)?
I've been working on solving Poisson problem using CNN model (you can ignore the Poisson problem part if you not familiar and jump to the image processing/CNN part). More specifically, I am solving electric potential problem ($\phi$) by CNN algorithm with particle distribution ($\rho$) as the input in Cartesian coordinates, similar to this but with no permittivity.

The training input ($\rho$) is $20\times20$ pixels matrix and the data label ($\phi$) is also $20\times20$ matrix. Either input and label is a 5000 matrix (5000, 20, 20, 1). The input-label data set was obtained from Gauss-Seidel method.

The goal is to obtain electric potential ($\phi$) from input data particle distribution ($\rho$) with CNN model.

With this regression case, what is the best CNN architecture to use?


问:
Video frame homography alignment
My project requires me to track objects such as cars via a drone.

I have a live video stream in which I would like to detect and track objects, the view doesn't change that much as the drone is nearly stationary moving only about 50 meters left and right at an altitude of 500m. For the proper tracking (and in general better video quality), I would need to stabilise the video using homography. I thought of these steps:

Align the images to the 1st frame with OpenCV ORB preferable so that the roads and building would be at the same exact location as in the previous frames but moving objects would be not accounted for. This approach is also a bit slow but I can just throw more power at it so it's not that crucial.

Track the objects using a tracking algorithm.

Unfortunately, when I tried this approach these are the results I got which where unsatisfactory:

Original stock video for testing (the actual drone won't move that much)

Frame 0 and frame 96 matches:

Frame 96 after alignment:

As you can see by the black bar on the right the alignment is completely off, I suspect the moving cars to be the problem. I thought of alining the frame x to the one before it but then the alignment error will build up over time.

Any ideas how I could improve my solution? Maybe some algorithm with knowledge of several past frames?

Small update:
My take at first was to eliminate the missalignment problem with a image flow detector like Lucas-Kanade. By detecting the major direction points of movement i could filter out the points that moved in a different direction to the rest. Then using the rest of the points (good points) to find a homography matrix to align the images. However that still wasnt good enough to position the cars. Now I think that I can align the image relative to the satelite imagery. That would eliminate the drifting problem but also introduces many others.

Title: Video frame homography alignment


问:
Why do we use modified pressure in incompressible multiphase solvers with gravity?
The context of my question is two-phase incompressible solvers such as interFoam in OpenFOAM, but I have seen this trick used elsewhere. To handle gravity these solvers subtract the hydrostatic pressure forming a modified pressure:

$p^* = p - \rho \boldsymbol{x}\cdot\boldsymbol{g}$

and then apply Chorin's projection to Navier-Stokes equation, deriving a Poisson equation. The effects of source terms will appear on the right hand side of the latter.

What is a formal justification for this trick and does it work always or under some circumstances only? The problems should be mathematically identical, but is clearly different numerically. Comparing the resulting Poisson equations a modified pressure will contain a divergence of $\boldsymbol{g} \cdot \boldsymbol{x} \nabla\rho$ and the other a divergence of $\rho\boldsymbol{g}$ with the remaining terms staying the same.

An explanation that I always bought into is that with a modified pressure and with clear separation between incompressible phases a zero solution for pressure is now valid in parts of the domain. Therefore, a zero vector is a reasonable initial guess. With the full pressure, we are asking the linear solvers to find the hydrostatic pressure so it has to work harder. Similarly, BCs can be also set to zero. However, a second derivative of $\rho$ can be quite large around the interface so I think in some circumstances modified pressure may be detrimental and possibly leading to instabilities?


问:
Why weighted harmonic mean for pchip slopes leads to monotone interpolator
In Fritsch and Carlson's paper on monotone interpolation, they identify numerous conditions under which a cubic Hermite interpolator will be monotone. For example: On the subinterval $[t_i, t_{i+1}]$ where $p(t_i) = y_i$ and $p(t_{i+1}) = y_{i+1}$, if we choose endpoint derivatives $d_i = d_{i+1} = 0$, then the interpolator is monotone. There are other conditions. Let $h_i := t_{i+1}-t_i$ and $\Delta_i := (y_{i+1} - y_{i})/h_i$. The condition
$$
d_i + d_{i+1} - 2\Delta_i \le 0
$$
is sufficient to ensure monotonicity. If $y_i' + y_{i+1}' - 2\Delta_i >  0$ then we could either have $x^{*} := \frac{1}{3} \frac{2d_{i} + d_{i+1} - 3\Delta_i}{d_{i} + d_{i+1} - 2 \Delta_i}$ not in the interval $[0, 1]$, or if it is in the interval $[0,1]$, then $\mathrm{sgn}(p(x^{*})) = \mathrm{sgn}(\Delta_i)$ will again guarantee monotonicity.

This is somewhat complicated, but there is quite a bit of freedom on a single interval. However, $d_i$ is constrained by data to both its right and left. It's clear that if $\Delta_{i-1}$ and $\Delta_i$ do not have the same sign, that $d_i = 0$ is the only solution. However, when $\Delta_{i-1}$ and $\Delta_i$ have the same sign, we need to find a choice of $d_i$ that satisfies the constraints listed above on both the intervals $[t_{i-1}, t_i]$ and $[t_i, t_{i+1}]$. Matlab solves this problem by setting $d_i$ to the "weighted harmonic mean"
\begin{align}
\frac{w_1 + w_2}{d_i} = \frac{w_1}{\Delta_{i-1}} + \frac{w_2}{\Delta_i}, w_1 :=2h_k + h_{k-1}, w_2 := h_k + 2h_{k-1}
\end{align}
I am unfamiliar with the properties of weighted harmonic means, and hence I cannot see why this choice for the derivatives satisfied the constraints for monotonicity listed in Fritsch and Carlson's paper.

How do we prove that the weighted harmonic mean leads to a monotone interpolator?

Partial solution: Using @MPIchael's suggestion, I was at least able to split of some cases. For example, if $d_{i+1}$ must be zero, then
$$
d_i - 2\Delta_i \le \max(\Delta_i, \Delta_{i-1}) - 2\Delta_i.
$$
And if $\max(\Delta_i, \Delta_{i-1}) = \Delta_i$, we satisfy $d_i + d_{i+1} - 2\Delta_i \le 0$. Otherwise, use $d_i \le \sqrt{\Delta_i \Delta_{i-1}}$, whereupon
$$
d_i - 2\Delta_i \le \Delta_i(\sqrt{\Delta_{i-1}/\Delta_i }-2)
$$
which is non-positive whenever $\Delta_{i-1} \le 4\Delta_i$.
I suspect to get the more general case, I need the inequality $H_{w}(\Delta_i, \Delta_{i-1}) \le 2\min(\Delta_i, \Delta_{i-1})$, which I know to be true for the standard harmonic mean, but I can't find a proof for the weighted harmonic mean.


问:
Resultant Lift Direction Confusion
Quite a small scale question but I have gotten myself confused. Inside Ansys Fluent I am simulating a wing at different angle of attacks, instead of rotating the airfoil I have rotated the inlet velocity, into its X and Y components as such that the flow will enter from a slight angle below the airfoi like so.

which would make X direction COS(AoA) and Y Direction SIN(AoA)
now heres where I've gotten myself stuck. when calculating lift force which is perpendicular to the flow, this requires a resultant force triangle to the left.

Now I cant decide whether it should be X = -SIN(AoA) and Y = COS(AoA) to resolve for lift

Title: Resultant Lift Direction Confusion


问:
MPM course, implicit updat clarifification
In the 2016 SIGGRAPH notes on MPM the authors explain how to do an implicit update in section 11.4 throguh by this formula:

It is not clear to me what the $j$ is over in this case (for eq 202), is it over the velocities of all grid nodes? Or is it just over the velocities of the nodes in a local neighborhood or over something else entirely?


问:
I encountered some issues simulating the shape transition of silicon clusters using LAMMPS
==

I have been trying to write an input file for lammps to simulate the shape transition of a nanoparticle of silicon using multiple potentials to obtain the best structure, which I will compare with a density-functional-theory simulation results. However, I have encountered an issue with the structure of the nanoparticle being distorted from the beginning when using the fix ID nve command, and I don't know how to fix it.

I've tried to change the fix nve command with the fix nvt command or just leaving the velocity command running alone, but the resulting structure is not as expected when compared to the density-functional-theory results.

Her is the input file that i have been using:

```<code>variable        dt equal 0.001
units           metal

variable t equal 300
variable t_start equal 300
variable t_final equal 1

# ----------------  Loop start -------------------------------------
variable i loop 12 26
variable cluster equal <span class="math-container">$i
variable epair_cluster equal epair/v_cluster
variable num_steps equal ($</span>{t}-${t_final})*1000

#-----------------------ATOM DEFINITION----------------------------
read_data       data/si_${i}.data
mass            1 28.085  ##Silicon

#------------------------FORCE FIELDS------------------------------
#pair_style      sw
#pair_coeff      * * Si.sw Si

#pair_style edip
#pair_coeff * * Si.edip Si

pair_style tersoff
pair_coeff * * Si.tersoff Si

#pair_style     tersoff/mod/c
#pair_coeff * * Si.tersoff.modc Si

neighbor        2 bin
neigh_modify    every 10 delay 0 check no

thermo_style custom step temp v_epair_cluster etotal v_cluster
thermo          1000

#dump            0 all xyz 100 min/si_min_$i.xyz
#dump_modify     0 flush yes
#dump_modify     0 sort id element Si
#minimize 1.0e-4 1.0e-6 100 1000
#undump          0

reset_timestep  0
timestep        $(dt)

#------------------------------ Setting the velocities  -------------------------
fix         1 all nve
velocity        all create  ${t_start} 4928459 rot yes dist gaussian

dump            5 all xyz 1000 si_<span class="math-container">$i_$</span>t.xyz
dump_modify     5 append yes

run     250000

unfix       1
#------------------------------- annealing -------------------------
fix             1 all nvt temp <span class="math-container">${t_start} $</span>{t} 0.1

dump            1 all xyz 1000 si_<span class="math-container">$i_$</span>t.xyz
dump_modify     1 append yes

run             ${num_steps}

#------------------------------- Relaxation  -------------------------
fix             2 all nvt temp  <span class="math-container">${t} $</span>{t} 0.1

dump            2 all xyz 1000 si_<span class="math-container">$i_$</span>t.xyz
dump_modify     2 append yes

run             250000

#-------------------------------  couling -------------------------
fix             4 all temp/berendsen <span class="math-container">${t} $</span>{t_final} 0.1

dump           3 all xyz 1000 si_<span class="math-container">$i_$</span>t.xyz
dump_modify    3 append yes

run             ${num_steps}

#------------------------------- Relaxation -----------------------
fix     5 all nvt temp <span class="math-container">${t_final} $</span>{t_final} 0.1

dump           4 all xyz 1000 si_<span class="math-container">$i_$</span>t.xyz
dump_modify    4 append yes

run            250000

 clear
 next i
 jump annealing_si_loop.lmp

print         "Simulation complete"
</code>```


问:
How to accelerate a convolution (laplace kernel) with FFT
I have the following computation I'm trying to program and accelerate with the FFT.

$$
\phi(x) = \sum_{y \in Y} K(x, y) q(x), \> \> \forall x \in X
$$

Where $X$ and $Y$ are sets of Cartesian coordinates in 2D. and $K(x, y) = \frac{1}{2\pi} \ln(\sqrt{(x_1 - y_1)^2 + (x_2 - y_2)^2})$.

I want to accelerate this convolution using the FFT, however I'm a little confused as to exactly how to formulate it.

In my case $x$ and $y$ denote points discretising the surfaces of two squares (for those who know it, this arises in the Kernel Independent formulation of the FMM algorithm).

I know from theory that the fourier transfrom can be applied to a kernel of this kind with the following reasoning, to accelerate the convolution

$$
q(x) * K(x-y) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} q(y) K(x-y)dy
$$
$$
q(x) * K(x-y) = \int_{-\infty}^{\infty} q(y) \frac{1}{\sqrt{2\pi}}\left ( \int_{-\infty}^{\infty} \hat{K}(\omega)e^{-i\omega(x-y)} d\omega\right ) dy
$$

$$
q(x) * K(x-y) = \int_{-\infty}^{\infty} \hat{K}(\omega) \left ( \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty}q(y) e^{i\omega y} dy \right ) dx
$$

$$
q(x) * K(x-y) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} \hat{K}(\omega) \hat{q}(\omega) d
$$

and that therefore,

$$
F[f * g](x) = F[K(x)]F[q(x)]
$$

In practice I don't think I'm formulating the calculation correctly.

I'm calculating the kernel as follows,

def K(x, y):
    diff = (x[0]-y[0])**2+(x[1]-y[1])**2
    return 1/2*np.pi * np.log(np.sqrt(diff))
```
<p>and for two sets of target/source points</p>
<pre class="lang-py prettyprint-override"><code>yc = np.array([1, 1])
xc = np.array([3, 1])
y = np.array([
    [0, 0],
    [1, 1],
    [1, 0],
    [0, 1]
]) + yc

x =  np.array([
    [0, 0],
    [1, 1],
    [1, 0],
    [0, 1]
]) + xc
</code>```

I'm calculating the fourier transforms as,

# Fourier transform of kernel values wrt to origin
origin = np.array([0, 0, 0])
y_s = np.zeros(len(y))
for i, source in enumerate(y):
    y_s[i] = K(origin, source)

y_k = fft(y_s)

# Fourier transform of signal
q_s = np.zeros(len(y))
for i, source in enumerate(y):
    q_s[i] = 1.0
q_k = fft(q_s)
```
<p>And performing the convolution/inverse transform.</p>
<pre class="lang-py prettyprint-override"><code>conv_fft = (y_k * q_k)
conv_fft = ifft(conv_fft)
</code>```

I'm not getting the same answer as with direct calculation, but don't understand why.

This is my code for direct convolution

conv_x_direct = np.zeros(len(x))
for i, yi in enumerate(y):
    conv_sum = 0
    for j, xi in enumerate(x):
        conv_sum += K(yi, xi)
    conv_x_direct[i] = conv_sum
```


问:
Example: Velocity Verlet reduced accuracy
Velocity Verlet is often held to far more accurate than forward Euler while being no more expensive. Technically, this requires some degree of regularity on the potential. But, is there a convincing example of a potential for which Velocity Verlet will be no more accurate than forward Euler?</p>


问:
What are K neighbor global scale and range global scale?
I am asking these terms as used in the following link. I understand k-neighbor but not k neighbor global scale.
When I used CGAL it for the 10 points 0 to 9 on x-axis(That is (0,0,0),(1,0,0),....(9,0,0)), it reported, Global K scale = 6;
Global range scale was 3.

But I don't understand why.
https://doc.cgal.org/latest/Point_set_processing_3/index.html#Point_set_processing_3Example_scale_estimation_global


问:
Algorithm for Probabilities in a "Random" Integers Runs test
Suppose you have a set of N random integers. ( aka bytes 0..255 )

You then consider each subsequent neighbor as "ascending" or "descending", counting the number of consecutive neighbors who are also of the same kind. ( i.e. 1234 is a "run" of 3 ascending neighbors; 123431 is a ascending run of 3 followed by a descending run of 2 ).  As an assumption, the problem considers both ascending and descending to be functionally symmetrical/same.  We can only consider runs of one direction.

In the set of N members, what is the probability of any Run of length L?   Such that one could test a series of bytes generated by a PRNG and gauge how closely the series is "Random" by how well its Runs fit the theorethical model.  ( This is only one test, and with the usual disclaimers about "determining randomness" in place. )

i.e.  In the set of N, how many runs of length 1 should there be?   Is there a formula for Runs-of-Length given N ( and the range of values of bytes ) ?

I feel that there should be a "simple" formula/algorithm, but I cannot figure it out and I have not found any similar information with searches.

Other similar tests with Runs use empirical results from a "golden" PRNG, but I feel that is unwise and really just testing similarity to some other PRNG ( which may have a serious unknown flaw/bias ).

An EXAMPLE data set, a specific series generated the following "rising" runs.

$N$ = $67108864$

```<code>[Length]=   num of occurs
[0]=    0
[1]=    13926757
[2]=    6144339
[3]=    1779212
[4]=    391946
[5]=    71071
[6]=    10716
[7]=    1390
[8]=    164
[9]=    16
[10]=   3
[11]=   0
[12]=   0
</code>```

Total number of Runs = $22325614$

As you can see, the distribution appears to be some sort of bell perhaps.  However, I am unable to figure out what kind and what the PDF is.  All series I generate from various "random" sources seem to yield similar results, but not exact ( of course ).

Thanks in advance!


问:
Numerical algorithm to compute second order solution of non-linear advection with source
I would like to solve the following equation,
$$\frac{\partial f}{\partial t}+\frac{\partial}{\partial x}\left({\mathcal A}f\right)=S(x,t)$$
where $f$ is a function of $x$ and $t$, ${\mathcal A}$ is a function of $f$ and $S$ is a source independent of $f$ but depends on $x$ and $t$.
My problem is I can find second order solution in $x$ of only advection equation using flux limiter but such solutions are not second order in time $t$.
Moreover, with source such equation introduces stiffness and often give negative solutions.
Is there any algorithm which can take care of these issues and provide second order in time and $x$.
Any help would be highly appreciated.


问:
Interpreting iterative smoothers and solvers as krylov preconditioners
Various literature and library implementations like petsc use preconditioners based on simple smoothers that themselves could be used the solve the systems directly. e.g. say I have a function

```<code>void do_jacobi_smooth(const matrix A, vector x, const vector b);
</code>```

which updates $x$ for the system $Ax=b$. Say I have a system where repeated calls to this eventually converge.

What I'm confused by in the literature of preconditioned methods, is that e.g. from the wiki page https://en.wikipedia.org/wiki/Conjugate_gradient_method for conjugate gradient you see the use of a preconditioner $M^{-1}$ written as

$$
z_0 := M^{-1}r_0
$$

so that $M^{-1}$ is some kind of (here linear) operator.

My question is what is the form of this operator, as implemented by something like the do_jacobi_smooth(A,x,b) above? In other words, say I have an unpreconditioned conjugate gradient implementation, and a jacobi smoother implementation (or a multigrid solver), how do I connect these parts to get a preconditioned CG?

Is it as simple as applying one iteration with a 0 initial guess?


问:
When is Lanczos tridiagonalization accurate?
Suppose that we are given a random, symmetric matrix $A$, and a random vector $q$. For concreteness, assume the dimensions of $q$ and $A$ are both $1,000$. I would like to use the Lanczos algorithm to compute the tridiagonal matrix, $T$, which is uniquely defined by $A$ and $q$. I don't necessarily need access to the entirety of $T$, but I want to understand any systematic errors which appear in the algorithm due to finite-precision arithmetic.

So far, I have computed the first $100$ rows and columns of $T$ by two methods. The first uses high-precision arithmetic:  I instantiate $A$ and $q$ with a huge number of digits of accuracy, and carry out the Lanczos algorithm as though numerical precision is effectively infinite. This should generate essentially exact results, and I have checked that orthogonality and normalization of Ritz vectors is maintained in this process. So far, so good.

My second approach uses standard machine precision, and in order to counter the loss of orthogonality among Ritz vectors, I explicitly re-orthogonalize each vector with respect to previously encountered vectors using the Gram-Schmidt procedure. I thought this would yield the same result as method (1), but in general it does not: I have found that after the first 50 or so iterations (this depends slightly on the initial matrix/vector), there is a significant divergence between the results of these two methods. I don't understand why this is the case; my understanding was that the loss of accuracy in the Lanczos algorithm was primarily attributable to the loss of orthogonality among these vectors.

So, my question: when is the tridiagonal matrix $T$, generated by standard precision arithmetic, reliable? Why isn't complete re-orthogonalization of the generated vectors sufficient?


问:
How to reconcile difference between Arpack manual and example?
The ARPACK manual for znaupd on pp. 128–129 says:

Mode 2:[...]
 OP = inv[M]*A and B = M.
 [...]
 IDO = 1: compute Y = OP * Z and
Z = B * X where
 IPNTR(1) is the pointer into WORKD for X,
 IPNTR(2) is
the pointer into WORKD for Y,
 IPNTR(3) is the pointer into WORKD for
Z.

I interpret that as:

When IDO=1, compute the matrix=vector product $$ Z = MX $$ and store
the result Z at IPNTR(3) then compute another
matrix-vector product $$ rhs = AZ $$ then solve the linear system
$$MY = rhs $$ and store the result Y at IPNTR(2).

However, the example zndrv3.f does this:

When IDO=1, compute the matrix-vector product $$ rhs = AX $$ Then
solve the linear system $$ MY = rhs $$ and store the result Y
at IPNTR(2).

Which seems contradictory since it lacks one of the multiplications by $M$ and doesn't store Z.

So I'm not sure what I should do. I tried to copy the example but it didn't work (seems like wrong solutions). There are other possible reasons it didn't work but this is one of them I want to confirm. There's a similar inconsistency for dsaupd and doing it like the examples works OK there.


问:
Theory behind the MAT_HIGH_EXPLOSIVE_BURN in LS-DYNA
I'm looking for the reference of the high explosive material model in LS-DYNA, which is represented by keyword MAT_HIGH_EXPLOSIVE_BURN. The theory manual has an entry for this material model but it is quite vague and does not point to a good reference. Does anyone know a good literature behind it?</p>


问:
Monte Carlo simulation of classical Heisenberg model doesn't represent Curie curve
I've created a JavaScript file to execute and log average energy and magnetization values of 2D lattice classical Heisenberg model.

I run the simulation with parameters,

L=32 (Number of atomic spins on each side of the lattice)

N=L*L (total number of spins)

J=1 (exchange constant)

T=0, 0.25, ..., 4.75, 5.00 (temperature values ranging from 0 to 5 with 0.25 steps)

k_B = 1.380649e-23 (Boltzmann constant, i'm not sure if I should use Boltzmann constant or not)

maxSteps=1000*N (maximum Monte Carlo steps per each temperature value, in this case 1,024,000 steps)

Question #1:

When I run the simuation with the code below I get random looking average magnetization values which doesn't look any similar to Curie curve where Average Magnetization approaches to 0 as temperate increases.

The sample values I got from runnig my code:

```<code>0.00 - 112.3611554761474,  0.25 - 23.5203320525913,
0.50 - 53.703480153558935, 0.75 - 92.0844269044712,
1.00 - 125.14550180873817, 1.25 - 74.21278454080122,
1.50 - 71.35901940700411,  1.75 - 34.8608142613192,
2.00 - 31.958062690297456, 2.25 - 87.57994096782724,
2.50 - 72.47576840153171,  2.75 - 68.64552406668619,
3.00 - 150.81073588464645, 3.25 - 124.51410153666252,
3.50 - 54.3372397486911,   3.75 - 20.10328424666724,
4.00 - 42.026588080636444, 4.25 - 52.42965924964965,
4.50 - 76.06067331506105,  4.75 - 17.969373459157673,
5.00 - 30.795866110306843
</code>```

Why might this be happening? I tried to run simulation for larger lattice (100x100 with 10,000,000 steps per each temperature value and temperatures ranging from 0 to 100 with 1 steps but got the similar results)

Question #2:

I previously run Ising model and got the same results for average magnetization-temperature graph with random scattering of points representing $M(T)-T$ values. I used the Boltzmann distribution to calculate the acceptance probability to decide whether or not to accept the spin flip (Math.random() < Math.exp(-dE/(k_B*T))). I found another way to calculate the acceptance probabilities in a book related to this subject and it was using the following function to calculate it:

// Transition probability calcualtion
function transitionProbability(t, neighboringSpins) {
    t=(t+1)/4;
    let e4 = Math.exp(-4/t);
    let e8 = Math.pow(e4, 2);
    let c = 0;

    if(neighboringSpins == 0) {
        c = 0;
    }
    if(neighboringSpins == 2) {
        c = e4;
    }
    if(neighboringSpins == 4) {
        c = e8;
    }
    if(neighboringSpins == -2) {
        c = 1/e4;
    }
    if(neighboringSpins == -4) {
        c = 1/e8;
    }

    return c;
}
```
<p>In this code variable  <code>neighboringSpins</code> represents the summed spin value of neighbors of the spin to be flipped which may be -4, -2, 0, 2, 4 since the spins may have -1 or +1 values (<span class="math-container">$\sigma_i=\{+1, -1\}$</span>). And this function calculates the acceptance probability for each of these cases rather than using Boltzmann  distribution. When I used this function to calculate acceptance probabilities the results for <span class="math-container">$M(T)-T$</span> plot was similar to Curie curve. But classical Heisenberg model, by definition, doesn't have spins with exact values like +1 or -1, rather the spins are unit vectors with random directions in 3D space (<span class="math-container">$\vec{\sigma}_i\in R^3, |\vec{\sigma}_i|=1$</span>). So why do I not get the same results for using Boltzmann distribution and mentioned acceptance probability calculation function? Why do I get better results with the mentioned function? What can I do to get similar results for classical Heisenberg model's Monte Carlo simulation which it's code given below?</p>
<p>The code for simulation:</p>
<pre class="lang-javascript prettyprint-override"><code>var L = 32;
var N = L*L;
var J = 1;
var T = 0;
const k_B = 1.380649e-23;

var system = Array.from(new Array(L), () => new Array(L));
var energyValues = new Array();
var magnetizationValues = new Array();

var maxSteps = 1000*N;
var stepsCounter = 0;

// Vector operations
function dotProduct(vectorA, vectorB) {
    return (vectorA[0]*vectorB[0])+(vectorA[1]*vectorB[1])+(vectorA[2]*vectorB[2]);
};

function addVectors(...vectors) {
    return vectors.reduce(
        (resultVector, vectorToAdd) => [
            resultVector[0]+vectorToAdd[0],
            resultVector[1]+vectorToAdd[1],
            resultVector[2]+vectorToAdd[2]
        ],
        [0, 0, 0]
    );
};

function flipVector(vector) {
    return vector.map(
        component => -component
    );
};

function magnitude(vector) {
    return Math.sqrt(
        Math.pow(vector[0], 2)
        +Math.pow(vector[1], 2)
        +Math.pow(vector[2], 2)
    );
};

// Initialize 2D lattice with random spins with magnitude of 1 unit
function initializeSystem() {
    for(let i = 0; i < L; i++) {
        for(let j = 0; j < L; j++) {
            let theta = Math.random()*(2*Math.PI);
            let phi = Math.acos(2*Math.random()-1);

            system[i][j] = [
                Math.sin(theta)*Math.cos(phi),
                Math.sin(theta)*Math.sin(phi),
                Math.cos(theta)
            ];
        }
    }

    // Calculate initial energy and magnetization
    let { energy, magnetization } = calculateEnergyAndMagnetization();

    energyValues.push(energy);
    magnetizationValues.push(magnetization);
};

function calculateEnergyAndMagnetization() {
    let energy = 0;
    let magnetization = [0, 0, 0];

    for(let i = 0; i < L; i++) {
        for(let j = 0; j < L; j++) {
            energy += -J*dotProduct(system[i][j], addVectors(system[(i-1+L)%L][j], system[i][(j+1)%L]));
            magnetization = addVectors(magnetization, system[i][j]);
        }
    }

    return {
        energy,
        magnetization
    };
};

function monteCarloStep() {
    let i = parseInt(L*Math.random());
    let j = parseInt(L*Math.random());

    let neighboringSpins = addVectors(
        system[(i-1+L)%L][j],
        system[i][(j+1)%L],
        system[(i+1)%L][j],
        system[i][(j-1+L)%L]
    );
    let dE = 2*J*dotProduct(system[i][j], neighboringSpins);
    let dM = flipVector(system[i][j]).map(component => 2*component);

    if(dE <= 0) {
        system[i][j] = flipVector(system[i][j]);
        energyValues.push(energyValues[stepsCounter]+dE);
        magnetizationValues.push(addVectors(magnetizationValues[stepsCounter], dM));
    } else {
        if(Math.random() < Math.exp(-dE/(k_B*T))) {
            system[i][j] = flipVector(system[i][j]);
            energyValues.push(energyValues[stepsCounter]+dE);
            magnetizationValues.push(addVectors(magnetizationValues[stepsCounter], dM));
        } else {
            energyValues.push(energyValues[stepsCounter]);
            magnetizationValues.push(magnetizationValues[stepsCounter]);
        }
    }
    stepsCounter++;
};

for(let t = 0; t < 21; t++) {
    T = t/4;
    stepsCounter = 0;
    energyValues = new Array();
    magnetizationValues = new Array();
    initializeSystem();

    while(stepsCounter < maxSteps) {
        monteCarloStep();
    }

    let avgEnergy = energyValues.reduce((sum, value) => sum + value, 0)/maxSteps;
    let avgMagnetization = magnetizationValues
        .reduce((sum, value) => addVectors(sum, value), [0, 0, 0])
        .map(component => component/maxSteps);

    console.log(`Temperature: ${T}`);
    console.log(`Average Energy: ${avgEnergy}`);
    console.log(`Average Magnetization: ${magnitude(avgMagnetization)}`);
}
</code>```


问:
Python code of explicit method of a nonlinear a BVP
I am trying to have a Python code for the following nonlinear BVP:
$$\frac{\partial N}{\partial t}=\frac{\partial^2 N}{\partial x^2}+N(1-N)-\sigma N$$  $$N(0,x)=\sin(2\pi x)$$
$$N(t,0)=0 \hspace{3mm}N(t,L)=0 \hspace{1.5mm};\hspace{2mm}t>0,0<x<L.$$

Without the non-linear part I discretize with explicit FTCS method.
$$\frac{N_{i,j+1}-N_{i,j}}{k} = c^2\frac{N_{i+1,j}-2N_{i,j}+N_{i-1,j}}{h^2}$$
$$N_{i,j+1} = \lambda[N_{i+1,j}-2N_{i,j}+_{i-1,j}]+N_{i,j}$$
$$\text{where}\hspace{3mm} \lambda = \frac{c^2k}{h^2} $$

But the problem is that maybe I can't discretize the nonlinear part. That's why my Python code does not respond well and I am new to Python. So, can anybody help me with a Python algorithm?


问:
Rayleigh-Ritz under LOBPCG : Nested iteration?
I am trying to understand the basic LOBPCG algorithm, as used in popular python libraries,
primarily from the suggested resource on wiki.

LOBPCG involves a subroutine to compute the coefficients of the orthonormal basis.
This computation is equivalent to another eigenvalue problem; it is equivalent to finding the optimal step size.
The Rayleigh-Ritz solves this eigenvalue problem, by projecting the given matrix on a smaller subspace.

I am not able to identify how exactly this projected eigenvalue problem is solved ??
Specifically, Algo2, step 4; "who is finally solving" ?

As suggested here, pg5., Rayleigh-Ritz utilize standard methods.

What are these standard methods ? QR,Power-iteration ?

Which of them are popularly implemented in python modules, such as scipy, torch etc. ?

Considering this method as an iterative algorithm, does it mean that LOBPCG have nested iteration ?