问:
Calculate the Bloch wave
The eigenvalue problem
$$\frac{d^2u}{dx^2}+2i k\frac{du}{dx}-[k^2-6\sin(x)^2]u(x)=-\mu u(x)$$

gives the first five eigenvalues with $k=0$ or $k=1$ which are $2.06$, $2.26$, $5.16$, $6.81$, and $7.74$.

And those eigenvalues I have calculated, but the corresponding eigenvector that I got did not agree with Fig.2 in the paper Z. Shi and J. Yang, "Solitary waves bifurcated from Bloch-band edges in two-dimensional periodic media," Phys. Rev. E, vol. 75, no. 5, pp. 056602, May 2007.

In this paper the author tells us the vector at the edge is all real, but I calculate it to be complex. Why is there an inconsistency with the paper?

This question on StackOverflow.

```<code>PRE75.056602

Clear["Global`*"]
Lx = \[Pi];
k0x = 2 \[Pi]/Lx;
V = 6 Sin[x]^2;
kx = 0;
{eig0, funs} = NDEigensystem[{-D[u[x],{x,2}]-2I kx D[u[x],x]+(kx^2+V)u[x],
                              u[0]== u[\[Pi]]}, u[x], {x,0,\[Pi]}, 3,
                              Method->{"VectorNormalization"->True,
                              "PDEDiscretization"->{"FiniteElement",
                              {"MeshOptions"->{"MaxCellMeasure"->0.01}}}}]

{2.06318, 6.81429, 7.74678},

kx=1;
{eig2,fun2} = NDEigensystem[{-D[u[x],{x,2}]-2I kx D[u[x],x]+(kx^2+V)u[x],
                            u[0]== u[\[Pi]]}, u[x], {x,0,\[Pi]}, 3,
                            Method->{"PDEDiscretization"->{"FiniteElement",
                                     {"MeshOptions"->{"MaxCellMeasure"->0.01}}}}]

{2.26673, 5.16594, 12.0926}
</code>```

Title: Calculate the Bloch wave


问:
How to implement adaptive step size Runge-Kutta Cash-Karp?
Trying to implement an adaptive step size Runge-Kutta Cash-Karp but failing with this error:

```<code>home/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:15: RuntimeWarning: divide by zero encountered in double_scalars from ipykernel import kernelapp as app
</code>```

The ODE i am trying to solve (and using in the example below, transformed from higher order to system of 1st order ODEs) is the following:

$$u''(x)+u(x)+u(x)^3 =0 ~~{\rm with}~~ u(0)=u'(0)=1 $$

Here is the implementation I am using:

```<code>def rkck(f, x, y, h, tol):
    #xn = x + h
    err = 2 * tol
    while (err > tol):
        xn = x + h
        k1 = h*f(x,y)
        k2 = h*f(x+(1/5)*h,y+((1/5)*k1))
        k3 = h*f(x+(3/10)*h,y+((3/40)*k1)+((9/40)*k2))
        k4 = h*f(x+(3/5)*h,y+((3/10)*k1)-((9/10)*k2)+((6/5)*k3))
        k5 = h*f(x+(1/1)*h,y-((11/54)*k1)+((5/2)*k2)-((70/27)*k3)+((35/27)*k4))
        k6 = h*f(x+(7/8)*h,y+((1631/55296)*k1)+((175/512)*k2)+((575/13824)*k3)+((44275/110592)*k4)+((253/4096)*k5))
        yn4 = y + ((37/378)*k1)+((250/621)*k3)+((125/594)*k4)+((512/1771)*k6)
        yn5 = y + ((2825/27648)*k1)+((18575/48384)*k3)+((13525/55296)*k4)+((277/14336)*k5)+((1/4)*k6)
        err = yn4[-1]-yn5[-1]
        h = 0.8 * h * (tol/err)**(1/float(5))
        yn = yn4
    return xn, yn

def integrate_sStepControl(f, t0, y0, tend, h, tol):
    T = [t0]
    Y = [y0]
    t = t0
    y = y0
    while (t < tend):
        h = min(h, tend-t)
        t, y = rkck(f, t, y, h, tol)
        T.append(t)
        Y.append(y)
    return np.array(T), np.array(Y)

def f_1(t,y):
    return np.array([ y[1], -y[0]-(y[0])**3 ])

Y0_f1 = np.array([1.0,1.0])

# Execution
h = 0.05
tv, yv = integrate_sStepControl(f=f_1, t0=0.0, y0=Y0_f1, tend=100.0, h=h, tol=1.0E-05)
print("[ %20.15f, %20.15f]"%(yv[-1,0], yv[-1,1]) )
plt.plot(tv, yv)
</code>```

Getting the error above, but it gets plotted. Don't know what I am doing wrong here :-/

ContentLicense: CC BY-SA 4.0


问:
Why does the correlation function of this stochastic differential equation starts at different points?
I am working with the following differential equation:

The equation is $$x=\beta +\sqrt{2D} \xi(t)$$ where $\xi(t)$ is a white noise term, with a reflecting wall boundary conditions.

After solving the equation numerically with Euler-Maruyama method, I used the Wiener-Khinchin theorem  (Using FFT) to find a numerical solution for the correlation function.
My problem is, get every time a new starting point for the correlation function even for exactly the same parameters, no matter how I change the parameters, and I am sure of my code because I have tried it on Orstein-Uhlenbeck and the results were right.

What is the problem of the varying starting points, or what did I misunderstand about correlation functions that is leading to this  ?? (Please see the code below, and the plots of the correlation function for the same given parameters)

This is the code for simulating the SDE with Euler method and simulating the correlation function

```<code> N = 8192
dt = 0.3
T = N*dt
time = np.arange(N)*dt
gamma = 0.1
D=1
v_factor = math.sqrt(2*D*dt)
v_data = []

v=1

###Eurler method to solve the SDE
for t in range(50000):  ##This part just to make sure we get the stationary process
        F = random.gauss(0,1)
        v = v - gamma*dt + v_factor*F
        if v<0:
            v=-v
for t in range(N):
        F = random.gauss(0,1)
        v = v - gamma*dt + v_factor*F
        if v<0:
            v=-v
        v_data.append(v)
# Compute numerically the autocorrelation via a Fourier transform

fft_cor = scipy.signal.fftconvolve(v_data, v_data[::-1])[N-1:]
fft_cor /= (N - np.arange(N))
plt.plot(time, fft_cor, 'k-', lw=2,label="Fourier-Transform")

# Compute the autocorrelation via the Wiener-Khinchin theorem

psd = np.fft.fft(v_data)*np.conj(np.fft.fft(v_data))/N
C = np.fft.ifft(psd).real

plt.plot(time, C,label='Wiener-Khinchin theorem')
</code>```


问:
The final Boundary Condition is Unknown, Is Backward Euler is still valid to be implemented?
I am working on conductive polymer modeling and supposed to do one-dimensional diffusion model in the thickness of the polymer, however, due to the small value of thickness in micro, when I use the forward approach it leads to instability, no matter how I decrease the number of nodes.

I have to use Backward Euler to counter this issue, however, due to the nature of our system, we can only know the initial boundary condition and not the final value, how I can solve this issue.

The question: Is there any alternative techniques I can use in numerical analysis while the final boundary condition is unknown?

The mathematical modeling:
For modelling the system, I am using RC transmission line, the mathematical description here: http://www.swarthmore.edu/NatSci/echeeve1/Ref/trans/Infinite.html


问:
Cell-segmentation from overstaining images Ask
I am currently trying to segment cells from digital pathology images.

The method I use is an algorithm based on color distance. This works for most of the cases, however, when dealing with the images which cells are overstaining, then the segmentation results is pretty bad. Is there any golden-standard method in treating these poorly stained with highly adhered-cell images?


问:
How to integrate the contents of a vector using an adaptive quadrature routine
I have a function which requires the return type to be a container. The problem is that I need to integrate the contents of the container as efficiently as possible and was hoping to use adaptive Gauss-Kronrod integration or something equally efficient (or better).

I was hoping to use the GNU Scientific Library quadrature routine qags, but it returns the result with type double.
As things stand I think I'm probably going to have to re-write sections of the GSL quadrature routines to convert the return type to std vector, but that would be a rather lengthy and possibly error-prone detour. I was hoping somebody might recommend a better solution!

Below is an example of the kind of function I'm trying to integrate over, presently employing the basic trapezoidal rule but indicating where I would instead prefer to implement the GSL routine gaq. whilst it is much simpler than my actual problem, it demonstrates that each element that is put in the vector is calculated from the previous result, hence the requirement of the container.

```<code>#include <iostream>
#include <vector>
#include <gsl/gsl_integration.h>
using namespace std;

vector<double> f(double E, int N) {
    vector<double> result;
    result.reserve(N);
    double x = E;
    for (int it=0; it < N; ++it){
        result.push_back(x);
        x *= x;
    }
    return result;
}

vector<double> f_qag(double E, void * params) {
    int N = *(int *) params;
    vector<double> result;
    result.reserve(N);
    double x = E;
    for (int it=0; it < N; ++it){
        result.push_back(x);
        x *= x;
    }
    return result;
}

int main(){
    vector<double> result;
    vector<double> integrate;
    int N = 20;
    result.reserve(N);
    integrate.reserve(N);
    for (int i = 0; i < N; i++)
        result[i] = 0.;
    double end = 1.0;
    double start = -1.0;

    // I would like to use qag here
    /* gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000); */
    /* vector<double> error; */
    /* gsl_function F; */
    /* F.function = &f_qag; */
    /* F.params = &N; */
    /* gsl_integration_qag (&F, start, end, 0, 1e-7, 1000, 1, w, &result, &error); */

    // instead of resorting to the trapezoidal rule here
    double E;
    const int n = 1000;
    double factor = (end - start)/(n*1.);
    for (int k=0; k<n+1; k++) {
        E = start + k*factor;
        integrate = f(E, N);
        for (int i = 0; i < N; i++){
            if ((k==0)||(k==n))
                result[i] += 0.5*factor*integrate[i];
            else
                result[i] += factor*integrate[i];
        }
    }

    for (int i = 0; i < N; i++)
        cout<<i<<" "<<result[i]<<endl;

    return 0;
}
</code>```

I intend to check for convergence against the pseudo-result conv;

```<code>double conv = 0.;
for (int i = 0; i < N; i++)
    conv += abs(integrate[i]);
</code>```


问:
Cover a polygon with least amount of parallelograms
I am solving the task that is as follows:

Input: a polygon. Can be any kind of polygon without self intersections. Can be a non-convex and with holes inside.

Goal: to cover it with 2 (at least) or more parallelograms that together are equal to or contain the whole polygon. The following criteria should be met:

There is no point of the polygon that lies beyond parallelograms.

The number of the parallelograms should be the least possible. I.e. we want to find largest parallelograms that cover the polygon.

These parallelograms must not intersect within the polygon.

Question: Is this task solved and if it is solved - how? I am looking for a direction where to start with and what related algorithms/theory to learn.


问:
ILUTP in sparse.linalg.spilu?
In Matlab, an ILU with threshold and pivoting (ILUTP) can be passed by default as:

```<code>setup.type = 'ilutp';
[L, U] = ilu(A, setup);
</code>```

Looking for an equivalent in scipy, ex:

```<code>sparse.linalg.spilu(A, options = ??)
</code>```

Tried to look up how matlab does the pivot and threshold but didn't seem to find much information.

Title: ILUTP in sparse.linalg.spilu?


问:
Extracting raw data from a graph
I am supposed to do black-box modeling, I need the input and output data for training, however, the authors didn't agree to share the raw data with me, I can only use their publication where these attached images are given below. Of course, I do understand the curve, however, I don't know how to create raw data from figures, is it could be achieved, if so what could be the easiest approach to have raw data from figures.

Title: Extracting raw data from a graph


问:
Calculating depth mask from different lighting
I have a object which is static, the camera is static and light source is moving. How can the depth mask be calculated ?

Concept is to use - calculate height from shadow length

Lets imagine a have more images likes this.

And light source is moving, which is sun, and from time it's position, distance, orientation can be calculated.

Yellow lines are sun rays.

Those circles are brightest pixel previously(direction of raw) of a very dark pixel. considering them highest i'n my knowledge . And those green lines as shadow. How can I calculated the depth mask of this terrain.

I'm ignoring the fact that, object/terrain can be colorful - no reflection/absorption.


问:
Singular Spectrum Analysis Explanation
I need you to help me understand the Singular Spectrum Analysis algorithm. I already read a lot of articles about the subject but they never answered my questions like what is the mathematical reason for embedding the time series into a trajectory matrix and why applying the SVD gives us access to such trend and periodic and noise functions.

In fact, a lot of people compare SSA to PCA for time series but one could easily explain PCA by saying that we want to find the most relevant direction for explaining the variance of the dataset and thus we are aiming to maximize the variance of the projection of the data on this particular direction we are looking for leading to this well-known optimization and eigenvector and diagonalization problem. But I am absolutely unable to find such an explanation for the SSA algorithm even by trying to explain it myself.

So if someone could help me with that I would really appreciate, in fact, it is really important for me to understand deeply how the things work in order to, for instance, understand what window length to choose or how the eigenvalues are related to the importance of the principal components.

Title: Singular Spectrum Analysis Explanation


问:
Double mach reflection at a inclined wedge
I am running into a strange problem when solving the 2D compressible Euler equations on a inclined wedge. To elaborate, my top boundary condition seems to emitting some type of instability.  I have tried using Neumann BCs for the u,v velocities, but the issue still persists.

I am using a 2nd-order finite volume method on a curvilinear mesh (structured grid).  My flux scheme is the Roe scheme with the Harten entropy fix.  I have attached a movie that illustrates the problem.

I was curious if someone has any thoughts.

https://www.dropbox.com/s/sev3o7olhqqx7h7/animation.mp4?dl=0


问:
Solving a non-convex optimization problem using fmincon
I am trying to solve a non-convex optimization problem using fmincon().
At each iteration, I am iteratively looking for the optimum value and when the termination criterion is satisfied, I keep this value as an initial value for fmincon() and run the code for the second iteration.

The problem is the code is running fine for the first to third iteration and then I get the following error

Error using sqpInterface Finite difference derivatives at initial
  point contain Inf or NaN values. Fmincon cannot continue.

Error in fmincon (line 823)
      [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = sqpInterface(funfcn,X,full(A),full(B),full(Aeq),full(Beq), .

How can I get rid of this error?


问:
Stability of SVD, Eigendecompositions, and pseudoinverse procedures in modern LAPACK routines
I have proposed an optimisation algorithm which I claim has improved upon the previous algorithm in a number of ways.

One of these claims is that my proposed solution requires no explicit SVD and psuedoinverse calculation which I claim is good, since in my case perturbations are important, and I believe these routines are traditionally unstable to perturbations.

My question is given the implementations in modern linear algebra computational libraries (such as LAPACK v3.0), how strong is this claim?

To summarise:

My method per iteration:

Only one solution of a generalized eigenvalue problem

Old method per iteration:

One solution of a generalized eigenvalue problem, AND,

One explicit SVD calculation for a Jacobian, $J = USV^{\intercal}$

Two explicit calculations of a pseudoinverse based on the Jacobian, $J$. That is, (1) $J^{\dagger}$ and (2)  $V^{\dagger}$

Given that my algorithm will have input perturbations, and considering the stability of modern linear algebra routines, is my claim a strong selling point?


问:
Changing geometry scale breaks simulation
I'm trying to find the capacitance matrix for a small array of metal boxes in air using Comsol 5.4. The geometry presented here is a simple geometry that recreates the problem I'm experiencing.

Simulation Settings: What Works

For now I set my geometry length unit to "um". To create my block array I first created a 40 x 40 x 20 unit block and then made a 3 x 3 array of these with 50 spacing in x and y. Meshing this with the default (Free Triangular) mesh gives me the following.

For my simulation I am using the "Electrostatics, Boundary Elements (esbe)" physics. I have my simulation domain selection set to "All voids".  I have my materials set to "Air" for "All Voids" and did not define a material for the blocks. I left all the default esbe boundary conditions (Charge Conversion, Zero Charge 1, Zero Charge 2, Initial Values). But then overrode nearly all of these by making all 6 boundaries of each block a voltage terminal.

This simulation works. I'm able to run my stationary source sweep and extract the Maxwell capacitance matrix and it gives me numbers that look reasonable.

What is Broken

I switched the geometry length unit to "nm". The resulting mesh looks identical to this one, but with "nm" as the unit instead of "um". Everything else remains unchanged. I tried running this simulation and got the following error:

```<code>There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term1.V0_ode.
at coordinates:  (0,0,0), ...
There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term2.V0_ode.
at coordinates:  (0,0,0), ...
There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term3.V0_ode.
at coordinates:  (0,0,0), ...
There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term4.V0_ode.
at coordinates:  (0,0,0), ...
There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term5.V0_ode.
at coordinates:  (0,0,0), ...
There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term6.V0_ode.
at coordinates:  (0,0,0), ...
There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term7.V0_ode.
at coordinates:  (0,0,0), ...
There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term8.V0_ode.
at coordinates:  (0,0,0), ...
There are 1 void equations (empty rows in matrix) for the variable comp1.esbe.term9.V0_ode.
at coordinates:  (0,0,0), ...
and similarly for the degrees of freedom (empty columns in matrix).
Returned solution is not converged.
Not all parameter steps returned.
</code>```

I have been able to manually extract the capacitance matrix from this geometry on the nm scale by manually setting terminal voltages and running a stationary simulation for each terminal. This is tedious, time consuming, and shouldn't be necessary. It is my understanding that a semi-manual terminal sweep should be possible by following the instructions here, but I have been unable to try this as "PortName" never shows up in my Parametric Sweep as a valid parameter (yes, I enabled "Activate Manual Terminal Sweep"). I'm using esbe instead of the es physics used in the tutorial. Maybe there is something else I have to do?

What can I do to make the stationary source sweep simulation work properly at the nm scale? Or if that is impossible, the parametric terminal sweep.


问:
How to determine the order of convergence of the Euler-Maruyama method?
This question is originally posted in Quant.StackExchange but has been unanswered for some time so I ask in here.

To make this simple let us consider the Geometric Brownian Motions (GBM).

My questions:

1. How can I show that the Euler-Maruyama Method is convergent using GBM?

2. How can I determine the order of convergence?

Background and the Problems I am facing:

According to Numerical Analysis theory An SDE solver has order $p$ if the expected value of the error is of $p$th order in the time step size

Now letting $S_t$ be a GBM with $S_0=1,$ $\mu=0.1$ and $\sigma=0.15$ then $E[S_{10}] =e$

When I run the solver $10,000$ times for different sizes $dt$, I would expect that the difference between my sample mean and true mean, $e$, will decrease as $dt$ gets smaller. However, when I run these simulations this is what I get:

This does not indicate that the error converges to zero as $dt$ goes to zero?

Why is that ..... Is this because of the randomness? I have actually tried for $1 \,000$, $10\,000$, and $20\,000$ simulations. In general, $10\,000$ simulations should be enough to remove much of variance.

Please suggest how I can determine/show the order of convergence.

If anyone is interested in the code

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

T = 10
mu = 0.1
sigma = 0.15
S0 = 1
ns = 10000

solution = S0*np.exp( mu*T )

dt_ = np.array([0.1,0.05,0.01,0.005])
err = np.zeros( len(dt_ ));
for j in range ( len(dt_ )):
    dt = dt_[j]
    Sn = np.zeros( (ns) )
    for i in range(ns):
        N = int(round(T/dt))
        t = np.linspace(0, T, N)
        ex= np.linspace(0, T, N)
        W = np.random.standard_normal(size = N)
        W = np.cumsum(W)*np.sqrt(dt) ### standard brownian motion ###
        X = (mu-0.5*sigma**2)*t + sigma*W
        S = S0*np.exp(X) ### geometric brownian motion ###
        Sn[i]= S[-1]

    mn          = np.mean(Sn)
    print(mn)
    err[j]      = abs( mn - solution)

plt.clf();
plt.loglog(dt_,err, color ="black", label = "Error (abs)");plt.xlabel("dt",fontsize = 20); plt.ylabel("Error (abs)",fontsize = 20)
plt.loglog(dt_,err, 'o', color ="black" , label = "x")
plt.title("loglog-plot",fontsize = 30);
plt.loglog(dt_,0.05*dt_**0.5, linestyle = ":")
</code>```

The code is mostly just copy-pasted from StackOverflow


问:
Fast approximate solver for vehicle routing problem
I need to solve capacitated asymmetrical vehicle routing problem with time windows on ~30k points. Time limits for calculations are 2 hours. I've tried using Clarke and Wright savings algorithm, it is very fast but too inaccurate, and I've tried using LKH-3 util for solving my task. It finds almost optimal solution, but unfortunately it takes up to several days to calculate. Is there some intermediate option with best solution quality on 2 hour runs?

Currently I'm using python 3 for solving this task, although any language is acceptable. The graph is the graph of garbage producers in modern cities and around 20% of the points have time bounds, constraining time it can be visited e.g. from 16:00 to 17:00. The edges of the graph is an estimated time, that a vehicle will spend traveling between nodes.


问:
Density functional theory: Total energies and forces
In DFT, forces are calculated using the Hellman-Feynman theorem, such that:
$$\frac{\text{d} E_\lambda}{\text{d} \lambda} = \left \langle \psi_\lambda \left|\frac{\text{d} \hat{H_\lambda}}{\text{d}\lambda} \right |\psi_\lambda\right \rangle$$

The primary question I have is: Assuming planewave DFT using the PAW method (Vasp, GPAW), is it strictly necessary to converge total energies to obtain correct forces both in the large force (atoms closer to each other) and small force (atoms at nearly equilibrium spacing) regimes?

I've run a set of tests on systems that are relevant to me, where I've noticed the following:

For a planewave cutoffs of 600 eV - 800 eV using a particular PAW dataset (total energies converge at 1500 eV), the energies are clearly not converged, but the forces are converged to ~0.02 eV/Angs in the high force regime (~ 4 eV/Angs), and have the same direction (determined by normalizing both force vectors and projecting one onto another). In the low force regime (~0.1 eV/Angs and below), the difference in the two forces remains approximately proportional, as above.

The final geometry of my structure using these two cutoffs was also identical to a maximum of ~0.005 Angs of differences in atomic positions.

I realize that formally, the total energy would have to be converged for the Hellman-Feynman forces to be converged, but are there considerations of the implementation using planewaves, or using the PAW method that make forces converge before total energies? Such as the energy not converging, but derivative wrt parameter converging much faster?


问:
Numerically Approximating the Jacobian and Comparing the Eigenvalues With Analytical Form
I am trying to study the stability of numerical discretization schemes using the Jacobian matrix of the residues with respect to the vector of conserved variables.

For a simple diffusion equation discretized with central difference scheme, the RHS can be written as:

$RHS=\nu\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2}$

This equation is interesting since the Jacobian matrix has an analytical form which can be written, in banded matrix notation as:

$J=\nu/dx^2 (+1,-2,+1)$

Basically, what I am trying to do is to compare the eigenvalues of the numerical approximated Jacobian with the eigenvalues of the analytical solution.

However, I am having some troubles because the results are not making sense. The numerically approximated Jacobian is presenting completly different eigenvalues when compared with the exact Jacobian. I am suspecting that the problem is in the way I am trying to compute the Frechet derivatives to compute the Jacobian.  The code fragment below shows how I am doing so:

```<code># Residue function of a simple centered three point scheme.

def frhs(um1,uo1,up1,dx,nu):
    return nu*( (up1 - 2.0 * uo1 + um1)/dx**2.0 )

# Computes the derivative of the residues with respect to the solution vector.

eps = 0.00001

drhs_du = np.zeros((nx-1,nx-1))  # In order to take the eigenvalues, this shall be a matrix.

e = np.zeros(nx-1)

# This loop computes the jacobian matrix according to http://www.netlib.org/math/docpdf/ch08-04.pdf

for i in range(1,nx-1):
    for j in range(1,nx-1):

        e[j] = 1.0

        drhs_du[i,j] = ( frhs(un[i-1]+eps*e[j], un[i]+eps*e[j], un[i+1]+eps*e[j],dx,nu) - frhs(un[i-1], un[i],un[i+1],dx,nu) )/eps

        e[j] = 0.0
</code>```

The complete code can be found in:
https://github.com/lmarmotta/disc_stability/blob/master/1dDiffusionEq.py

Am I doing anything wrong ?

How would you guys approximate the Jacobian numerically ?  And Are there any other algorithms to approximate Jacobian matrices for these type of cases ?

Thank You !


问:
Metis: how to use and tutorial recommendation
I am new to METIS and trying to use it in my fortran code. I read the manual online. But still, I am not sure about how to implement it my code. I tried the test cases in the graphs directory. For example, I partitioned the 2d mesh metis.mesh by using mpmetis -gtype=nodal metis.mesh 4. And it gave me two output files:metis.mesh.npart.4 and metis.mesh.epart.4. I opened these files, they are point values that make no sense to me. How can I visualize the partition results? And is there any good tutorial about metis that can help me to start? Thank you!


问:
Neumann boundary conditions on arbitrary surface for finite difference diffusion
I am facing the following problem, formulated in practical terms: I have a region $\Omega$ in two or three dimensions, represented as a binary mask, and an initial density $u_0$ within that region that I'd like to evolve following a diffusion equation:

$$\partial_t u(\vec{x},t)=\nabla(D\nabla u)+f(u)$$

where D can be a quadratic matrix. I'm struggling with how to ensure Neumann boundary conditions of the form $(D\nabla u)\cdot\vec{n}|_{\partial\Omega}=0$, i.e. the direction of diffusion should be parallel to the surface along the entire surface. All I could find (and actually understand...) was in 1D, which I could translate to rectangular regions and scalar $D$, but now I'm trying to move to arbitrary regions. My first guess was to get the surface normals (e.g. from the gradient of a distance transform) and then project diffusion vectors along the surface onto a plane orthogonal to the surface normal, but I have no idea whether that's a good way to do it.

Are Neumann boundary conditions supposed to imply conservation of mass? Because I can't see how that's going to work in a discrete scheme. If it's in any way relevant, I'm using a simple forward model in time and central differences in space. Do I need an implicit time scheme instead? Any help is much appreciated!


问:
Strange spectral (finite) element results of a solid plate
I tried to implement the spectral element method1 as proposed in [1]. I simulated an aluminum plate and a vertical concentrated force was applied at the middle of the top left quarter of the plate2. I should get a wave with circular crest, coming from the point where I exerted the force. But apparently the wavefront got inverted as illustrated in the following picture.

(This is the z displacement at a particular instant.)

<s> There must be something wrong in my code </s>. What could have happened? <s>
 Something wrong during assembling of the stiffness and mass matrix</s>? This link  contains my python source code in pdf format and the video that illustrates the propagation of the wave. I will try rewrite my code with ordinary plate element first.

1Here, SEM is just finite element method with higher order polynomial shape functions. Some other works use sinusoidal shape functions.

2I didn't abuse symmetry since I will add something to the plate later.

Update: I tried using ordinary plate element instead, but the results are the same given that the mesh is fine enough.

Edit 2: I think my code is actually correct, given that the Abaqus simulation also show such circular pattern on the top, bottom, left and right of the main crest (they are more obscure though). But I am still puzzled about why those circular pattern appear. I guess it has to do with the concentrated load and the effect of such spurious pattern would be weakened only if I refine the mesh (h-refinement). But from my observation, p-refinement does not help.


问:
Finite dimensional optimization problem over dynamical system
I am interested in solving numerically the following mathematical problem
Consider an ode of the form
$$
\dot q(t) = f(q(t),t_1,\ldots, t_N),\qquad t\in [0,T],
$$
where $q\in \mathbb{R}^n$ is the state and $t_i,\, i=1,\ldots N$ are real parameters such that $0\le t_1 \le \ldots \le t_N\le T$.

Moreover $f$ is defined from the heaviside function: $H(t-t_i)$ (so I don't want to use auto differentiation)

I would like to solve:
$$
\min_{t_1,\ldots,t_N}\; \lVert q(T)\rVert
$$

Here is what I have done so far. I am using the scipy.optimize.minimize() function together with the 'SLSQP' method and I am integrating the ode using the 'dopri5' method from the scipy ode integrator.

Here some "pseudo code":

```<code>n = 2 #dimension of q
N = 10 # numbers of ti
qf = [] #vector containing the state to compute the objective function
tis0 = ... # containing initialization of t1,t2, ...

#init integrator
solver = ode(dyn).set_integrator('dopri5')

#function to compute the dynamcis
def dyn(t,z,n,N):
    #compute dq/dt
    dqdt = [0. for i in range(0,n+N)]
    dqdt[...] = ... #heaviside function applied to (t-ti) appearing here
    dqdt[n:-2] = 0. #dynamics for the (constants) ti
    return dqdt

#compute the state q
def solout(t,q):
    qf.append(q)

#function to compute the objective
def objfun(tis):
    solver.set_solout(solout)
    y0[n:n+N] = tis
    solver.set_initial_value(y0,0.).set_f_params(n,N)
    solver.integrate(T)
    return norm(qf[-1])

#argument for the optimization routines : bounds, constraints...
bnds = tuple([(0.,T) for i in range(0,N)])
cons = ({'type':'ineq', 'fun': lambda x: x[0] < x[1]},...)
#call to the optimization routine
res = minimize(objfun,tis0,method='SLSQP',bounds=bnds,constraints=cons)
</code>```

I am not convinced that this is a suitable approach for this problem. In fact I am observing that the algorithm is converging but the optimal solution is very close to the initialization vector tis0 ... so tis' variables don't move. Maybe I shouldn't define the tis as constants and don't put them in the dynamics ?

Also I wonder whether SLSQP is the right method to tackle this problem.

Any thoughts even not related to this solution is pretty welcome.

EDIT:
the dynamics is given by:
\begin{array}{l}
&\dot x(t) = y(t)\, \gamma(z(t)) -x(t)\, \beta(z(t)) \\
&\dot y(t) = x(t) - y(t) + 1
\end{array}
where $\gamma(z(t))=\frac{z(t)}{z(t)+1}$, $\beta(z(t))=\frac{1}{1+\gamma(z(t))}$
and
$$
z(t) = \sum_{i=1}^N R_i(t_i)\,(t-t_i)\, \exp(-(t-t_i))\, H(t-t_i)
$$
where $R_i(t_i) = 1+\exp(-(t_i-t_{i-1}))$ and with the conviention $t_0=-\infty$.


问:
Finding the polynomial for the solution of an ODE
I’m stuck trying to solve part (b) and (c) of the below problem, but part (b) is the one of main concern here as I think (c) should follow easily once (b) is completed.

I don’t know where to start with part (b). My first idea was to try to plug in $y’=-\lambda y$ in K1 instead of $f(t_n,y_n)$. I didn’t get far with that approach as the second step has the expressions $t_n+\frac{1}{3}h$ and $y_n+\frac{1}{3}K1$ in the arguments of $f$, and I don’t know what would relate to our $y$ here. I assume that we would get some expression with different powers of $\chi=h\lambda$ so that we could write this as a polynomial of $\chi$. I feel like I’m completely lost here, and looking in Applied Numerical Methods by Steven Chapra I don’t find much of use. I kindly ask for any help possible, and gladly provide clarifications. Thanks in advance!


问:
Stokes Equation fails to converge for an ellipse
This might be because of the mesh, but the following code blows up for all values of b not 1. Does anybody have any experience working with the ellipse mesh in Fenics?

import matplotlib
matplotlib.use('Agg')
from matplotlib import ticker
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *

def epsilon(u):
  return grad(u)+nabla_grad(u)

b=1
embryo = Ellipse(Point(0.0, 0.0), 1, b)
mesh = generate_mesh(embryo, 32)

# Define function spaces
P2 = VectorElement('CG', triangle, 5)
P1 = FiniteElement('CG', triangle, 3)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)
g = Constant(0.0)
mu = Constant(1.0)
force = Constant((0.0, 0.0))

# Specify Boundary Conditions
boundary = 'on_boundary'
flow_profile = ('-sin(atan2(x[1]*b,x[0]))*(a0+a1*sin(2*atan2(x[1]*b,x[0])))','cos(atan2(x[1]*b,x[0]))*(a0+a1*sin(2*atan2(x[1]*b,x[0])))')
bcu = DirichletBC(W.sub(0), Expression(flow_profile, degree=5, a0= -0.4202, a1 = 0.7653, b=b), boundary)
bc = [bcu]

# Define trial and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

a1 = inner( mu*epsilon(u)+p*Identity(2), nabla_grad(v))*dx + div(u)*q*dx
L1 = dot(force,v)*dx + g*q*dx
print('Preliminaries Done')

# Solve system
U = Function(W)
solve(a1==L1, U, bc)
print('Solving Done')

# Get sub-functions
u, p = U.split()

folder='./'
fig=plot(u, title='Velocity')
plt.colorbar(fig)
plt.savefig(folder+'vel.png', dpi=1000)
plt.close()
print('Plotting Done')
```


问:
Computing Trajectory Equations of Kerr Geodesics
I want to numerically solve the trajectory equations of a Kerr geodesic given by wikipedia in Matlab.

The trajectories look like:

I implemented the equations and solved it with the standard Runge-Kutta method. Here it is my code:

```<code>% Initial parameter (r, theta, phi, t)
r       = x(1);
theta   = x(2);
phi     = x(3);
t       = x(4);

% Constants of motion
mu = const(1);
E = const(2);
L = const(3);
Q = const(4);

% Abbreviations
Sigma = r^2 + a^2 * cos(theta)^2;
Delta = r^2 - 2 * M * r + a^2;
K = Q + (L - a*E)^2;
X = a^2*(E^2 + mu) - L^2 - Q;
R = (E^2 + mu)*r^4 - 2*M*mu*r^3 + X*r^2 + 2*M*K*r -a^2*Q;
O = Q + cos(theta)^2*(a^2*(E^2 + mu) - L^2/sin(theta)^2);

% ODE from O Neill Chapter 4.2.2 Geodesics
dx(1) = sqrt(R / Sigma^2);
dx(2) = sqrt(O / Sigma^2);
dx(3) = ((L - a*sin(theta)^2*E) / sin(theta)^2 + a*((r^2 + a^2)*E - L*a)/(Delta) ) / Sigma^2;
dx(4) = (a*(L - a*E*sin(theta)^2) + (r^2 + a^2)*((r^2 + a^2)*E - L*a)/(Delta)) / Sigma^2;
</code>```

I am a little confused by the R- and Q-term because there is a square root. How can I handle the signs of these in the calculations?

I know there exists an alternative form of the equations with the components of momenta, but I do not want to use this approach.


问:
Dealing with spurious oscillations in particle tracking methods
I work on modelling high intensity discharge xenon-filled lamps.
The model governing the discharge is quite complex and sadly includes fluid dynamics.
After some time, I managed to implement a finite-difference particle-tracking lagrangian fluid dynamics for axisymmetric 1D discharge:

$$
\frac{\partial}{\partial t} \left(\frac{1}{\rho}\right) = \frac{\partial(rv)}{\partial s}
\\
\frac{\partial v}{\partial t} = - r \frac{\partial p}{\partial s}
\\
\frac{\partial r}{\partial t} = v
\\
\frac{\partial}{\partial t} \left( \mathcal{E} + \frac{v^2}{2}\right) = - \frac{\partial}{\partial s}(prv) - \frac{\partial(rW)}{\partial s} + \frac{f}{\rho}
\\
W = -\lambda \frac{\partial T}{\partial r}
\\
p = p(\rho, T)
\\
\frac{\partial\mathcal{E}}{\partial T} = c_v(\rho, T)
$$

Computational cells are organised so that lagrangian particle speeds and eulerian coordinates are defined on cell-interfaces and energy equation parameters (temperature, density, pressure, thermal capacity, thermal conductivity) are defined on cell centers (image below). Domain boundaries coincide with cell interfaces.
.

Right-side boundary conditions are dirichlet for speed $$v(r)|_{r=R} = 0$$, keeping the scheme homogenous.
For energy balance equation it gets a bit trickier as there is no boundary node to specify dirichlet condition at, so I intruduce fictious node in the depth of the wall material and interpolate boundary condition at cell interface:

$$
\frac{\hat{T}_{last} + \hat{T}_{last + 1}}{2} = T_{boundary}
\\->\\
\hat{T}_{last + 1} = 2T_{boundary} - \hat{T}_{last}
$$

For computational efficiency reason I solve the equations separately. There's iterational process on energy equation; first I solve energy equation with implicit O(t^2)O(h^2) scheme, then I advance gas dynamics part on interpolated temperature profile with explicit scheme with time step several orders smaller to refine coefficients for energy equation. This is dictated by the extremely high cost to recalculate energy source function as there's radiative transfer over complex spectra involved.

Model works great for pulses with current of about 500A and duration greater than 50e-6 seconds, however once I go into kA ranges waves rebounding from the lamp wall tend to spuriously oscillate and destroy the model upon converging in the center ().

I tried some less radical methods of dealing with these oscillations like introducing artificial viscosity, however I found that my implementation of them greatly reduced stability due to loss of conservation.

With great dread I look forward to rewriting everything from scratch in terms of TVD scheme or with nested grids or with something equally time-consuming, and my question is: what are some less intrusive approaches to eliminating spurious oscillations that I overlooked?

Edit 1: added more information about energy balance equation, schema, and boundary conditions used.


问:
Simulating Anderson model, have problem with momentum representation (MATLAB)
I want to change from real-space representation to momentum-space representation I have a Hamilton-operator (Anderson-model), and I calculated some kind of entropy of its eigenstates (this is working, I see what I want). Next I want to change to momentum repr. using FFT and now my eigenstates in momentum space are not normalized. E.g. if I calculate the sum of the square of the eigenstates, it has to be 1, but it does not work.

I tried to sum the square of the eigenstates and normalized with them, but It does not work (I show the code without any failed trying).

```<code>N=100; %dim of matrix
Nx=15; %number of points
%because of log scale
xmin = -3.0;
xmax =  3.0;
dx = (xmax - xmin)/(Nx-1);
x = zeros(1,Nx); %x axis pre
ss = zeros(1,Nx); %entropy pre
spp=zeros(1,Nx); %entropy in Fourier space pre
eps=1.0e-6;

for ix=1:Nx
    %log scale
    x(ix) = xmin + (ix-1)*dx;
    xx = 10.0^x(ix);

    average_s=0;
    average_spp=0;
    %anderson modell
    W=xx;
    r=rand(1,N)*W-(W/2);
    A=diag(ones(1,N-1),1)+diag(ones(1,N-1),-1)+diag(r);
    %diagonalization
    [V,D]=eig(A);
    %PROBLEM HERE:
    %Fourier transformation
    P=fft(V)/(sqrt(2*pi)*N);
    P=abs(P);
    for j=1:N
        four_sum=0; square_sum=0; entropy=0;
        four_sum_p=0; square_sum_p=0; entropyp=0;
        for i=1:N
            %Fou
            probp=(P(i,j)).^2;
            square_sum_p=square_sum_p+probp;
            if probp>eps
                entropyp=entropyp-probp*log(probp);
            end;
            four_sum_p=four_sum_p+probp.^2;

            %Real
            prob=V(i,j).^2;
            square_sum=square_sum+prob;
            if prob>eps
                entropy=entropy-prob*log(prob);
            end;
            four_sum=four_sum+prob.^2;
        end
        qp=square_sum_p.^2/(four_sum_p);
        average_spp=average_spp+entropyp-log(qp);

        q=square_sum.^2/(four_sum);
        average_s=average_s+entropy-log(q);
    end
    ss(ix)=average_s/N;
    spp(ix)=average_spp/N;
end
plot(x,ss,x,spp);
</code>```

The structural entropy in real space (ss vector) has the correct form, but in momentum space (spp) after FFT is not look like what I want, and it is not normalized.


问:
Fit exponential convergence
I'm working with a numerical algorithm whose output $y$ asymptotically approaches a certain unknown value $a$.

I expect an exponential convergence, i.e. the data $y$ given by my algorithm should be distributed as

$$
   y=a+b e^{-c x}
$$

where $x$ represents the number of iterations of the algorithm. For computational purposes, it's good to terminate the algorithm when $x$ is still relatively small, say when $x=x^*$. In this case, however $y(x^*)$ hasn't fully reached its asymptotic value $a$.

I kindly ask you to provide a Matlab code which takes some tens of pairs $(x,y)$ and, retrieves the value $a$ appearing in the aforementioned equation.

Title: Fit exponential convergence


问:
Compute cell values and isosurface constant from density values of particles
I am trying to reconstruct the surface for a fluid simulation based on a list of particles using the Marching Cubes algorithm. From different resources, such as http://paulbourke.net/geometry/polygonise/, or https://cseweb.ucsd.edu/classes/sp16/cse169-a/slides/CSE169_15.pdf

I see that I have to compute an isosurface based on a 3D scalar field. My question is, how can I set density values to the corners of the cubes for each cell of the grid and how should I compute the isosurface constant. Is that constant same for every cell, or unique for each one?

OwnerDisplayName: Petros Yannopoulos


问:
Numerically computing deflection due to thermal expansion
Using linear elasticity formulation, I am attempting to numerically compute the displacement due to thermal expansion. This is done for a 3-D isotropic material.

The governing equations are simply:

$$
\nabla \cdot \sigma = 0 \\
\sigma = \sigma_{\text{thermal}} - \sigma_{\text{mechanical}} \\
\sigma_{\text{thermal}} = C'\epsilon_{\text{thermal}} \\
\epsilon_{\text{thermal}} = \begin{bmatrix}
\alpha \Delta T \\
\alpha \Delta T \\
\alpha \Delta T \\
0 \\
0 \\
0 \\
\end{bmatrix} \\
C' = \text{4th order stiffness tensor reduced to a 2nd order 6x6 matrix} \\
\sigma_{\text{mechanical}} = C'\epsilon_{\text{mechanical}} \\
\epsilon_{\text{mechanical}} =
\begin{bmatrix}
\frac{du}{dx} \\
\frac{dv}{dy} \\
\frac{dw}{dz} \\
\frac{dv}{dz} + \frac{dw}{dy} \\
\frac{du}{dz} + \frac{dw}{dx} \\
\frac{du}{dy} + \frac{dv}{dx}
\end{bmatrix}
$$

From this, we can formulate an expression for the displacements, $u$,$v$,$w$. After writing these equations, I realized that it appears that the stiffness tensor has no bearing on the deflection generated from a thermal load. Mathematically, because $C'$ is invertible, we are left with the equation:

$$
\nabla \cdot \epsilon_{\text{mechanical}} = \nabla \cdot
 \epsilon_{\text{thermal}} .
$$
I am being a little lazy here with the notation where I previously
used $\epsilon$ to represent a $6\times1$ vector, here it is a $3\times 3$ matrix.

I was able to verify this in ANSYS by varying the structural properties of the material and it demonstrated no change in the deflection from the thermal loading.

But physically, why is the deflection due to thermal expansion independent of the structural properties?


问:
Testing the SUPG method and other methods for hyperbolic equations
I am interesting in integrating the simple equation
$$
\frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla \phi = 0
$$
with a Dirichlet boundary condition at the influx boundary ($\mathbf{u} \cdot \mathbf{n} < 0$). The domain is a simple unit square. I pretend to use SUPG and I want to validate the method using the method of manufactured solutions. I came up with a simple solution

$$\phi = \begin{cases} cos(5x - t) ~ \text{for} ~ y > 0.5 \\
0 ~ \text{else}
\end{cases}$$

with a advective field $\mathbf{u} = (1.0/5.0, 0)$. This makes the source term of the equation equal to zero and I adjust the Dirichlet BC appropriately

Discretizing this equation with just Lagrange elements, no stabilization methods added, and integrating in time with Crank-Nicholson results in nice solutions with no oscillations that I would expect from using the Lagrange elements on a hyperbolic equation with a discontinuous function being advected. Why is it so? Are there other manufactured solutions that will actually trigger oscillations?


问:
Optimization (best input variables search) for a non-smooth non-linear unknown function
I am trying to optimize a system that monitors and advises a user multiple times over a certain period of time depending on changing outside factors. The systems behavior can be altered by 5 constrained input variables, i.e. weights of a cost function. These parameters are fixed for one run.

At the end of the execution, the systems performance can be summarized by 3 output parameters (1 real, 2 integer) that are unrelated to the input variables.

The goal is to find those input values that minimize the output parameters.

(I have looked at classic optimization problems, but since there is no objective function that can be defined and the outputs are non-smooth / non-linear many don't seem to be applicable. From what I found "derivative-free optimization" might work, but I am completely new to the topic)

The system's code is written in JAVA and one execution takes about 10min.

Any suggestion on what solver / search algorithm could be used for this problem (ideally with a JAVA library available)?


问:
Trouble with creating correct element matrices in Finite Element Analysis for a cantilever beam
I'm trying to solve for displacements of a cantilever beam numerically with FEA. The beam is modeled as a 3D-solid made up of a set of 8-noded hexahedral elements, which are in their undeformed state, constructed as cubes. Every element has the same volume before solving for the displacements, where each side is 1 cm. See picture below:

I apply a uniformly distributed load at the 4 outermost nodes on the top of the beam in the positive z-direction (global coordinates). And the result I get is this horribly deformed structure!

Obviously there is something faulty with my solving, the beam should simply bend without these awful artifacts. I will now try to explain the things I have considered so far, bare with me:

Node numbering

For each element, I have numbered the 8 nodes (global node numbers) counter-clockwise starting from the bottom left doing a full rotation of the "bottom face" and then the same for the "top face", (see attached image below). With this numbering schema I've managed to get a positive Jacobian for every element stiffness matrix. It does become negative, as expected if the numbering is instead done clock-wise. However the Jacobian is very small, in fact it is almost constant for every element stiffness matrix, being for example 1.250000000000003e-07, only differing slightly with respect to the last decimal.

Eigenvalues of element matrices

When testing for the eigenvalues of a stiffness matrix many of the element stiffness matrices have a total of 6 eigenvalues that are zero, corresponding to the rigid body displacements stated in an answer to a similar post here. In that post the answerer states that:

If you get more than six, your stiffness matrix doesn't have the
  necessary rank and is therefore defective.

My element matrices are 24x24, as each node has 3 DOF's. Should there be a total of 6 eigenvalues that are zero if I have modeled everything only with 3 degrees of freedom? I reckon that the additional 3 values are in respect to the rotational degrees of freedom as well, even though I have only modeled my problem with the 3 translational DOF's.

I'm creating the element matrices with a total 2 Gauss points in each of the 3 directions: since the shapes functions are linear the number of Gauss points should be enough.

```<code>for i = 1:numgp
    for j = 1:numgp
        for k = 1:numgp
            xi = gp(i);
            eta = gp(j);
            zeta = gp(k);
            %In here I compute Jacobian, strain matrix
            mat = mat + w(i)*w(j)*w(k)*B'*c_mat*B*J;
        end
    end
end
</code>```

and compute the derivatives in natural coordinates for each of the 8 shape functions: for example the derivative of the first shape function:

```<code>dN1dXI = -(1/8)*(eta-1)*(zeta-1);
dN1dETA = -(1/8)*(xi-1)*(zeta-1);
dN1dZETA = -(1/8)*(xi-1)*(eta-1);
</code>```

And construct the gradient matrix as follows:

```<code>%Gradient matrix of shape functions N1-N8
GN = [
    dN1dXI, dN2dXI, dN3dXI, dN4dXI, dN5dXI, dN6dXI, dN7dXI, dN8dXI;
    dN1dETA, dN2dETA, dN3dETA, dN4dETA, dN5dETA, dN6dETA, dN7dETA, dN8dETA;
    dN1dZETA, dN2dZETA, dN3dZETA, dN4dZETA, dN5dZETA, dN6dZETA, dN7dZETA, dN8dZETA;
];

%Jacobian matrix
%C is a 8x3 matrix with the global (cartesian) coordinates
%the columns x y and z where the i:th row is the i:th node in the element.
Jmat = GN*C;
%Jacobian
J = det(Jmat)

Bmat = Jmat\GN;

B1x = Bmat(1, 1); B2x = Bmat(1, 2); B3x = Bmat(1, 3); B4x = Bmat(1, 4);
B5x = Bmat(1, 5); B6x = Bmat(1, 6); B7x = Bmat(1, 7); B8x = Bmat(1, 8);

B1y = Bmat(2, 1); B2y = Bmat(2,2); B3y = Bmat(2,3); B4y = Bmat(2, 4);
B5y = Bmat(2, 5); B6y = Bmat(2, 6); B7y = Bmat(2,7); B8y = Bmat(2,8);

B1z = Bmat(3, 1); B2z = Bmat(3, 2); B3z = Bmat(3, 3); B4z = Bmat(3, 4);
B5z = Bmat(3, 5); B6z = Bmat(3, 6); B7z = Bmat(3, 7); B8z = Bmat(3,8);

B1 = [B1x 0 0; 0 B1y 0; 0 0 B1z; 0 B1z B1y; B1z 0 B1x; B1y B1x 0];
B2 = [B2x 0 0; 0 B2y 0; 0 0 B2z; 0 B2z B2y; B2z 0 B2x; B2y B2x 0];
%..and so on for all the 8 nods

B = [B1 B2 B3 B4 B5 B6 B7 B8];
</code>```

This will however not result in numerically symmetric element stiffness matrices. They are almost symmetric but differ slightly in the values causing them not to be exactly symmetric.

Boundary conditions

Before solving for the displacements I simply remove the first 48 rows and columns, in my assembled (global stiffness matrix), that correspond to the Dirichlet boundary conditions of the cantilever beam, that is, there are a total of 16 nodes (each having 3 translational DOF's) that should have zero displacements.

I appreciate any help I can get, I could include more code snippets but as you might understand there is a whole lot of code so I'll keep it rather stripped down as for now.


问:
Best optimizer for unconnstrained non-convex nonlinear least-square optimization problem?
I am looking for a very good optimizer to the following problem:

$$\min_{P,\Theta}\lVert APD(\Theta)P^{-1} -B \rVert_F$$
where $A,B \in \mathbb{R}^{n\times m}$, $P \in \mathbb{R}^{m\times m}$, $D\in \mathbb{R}^{m\times m}$.

For my problem, typically $n = O(10^3)$, $m = O(10^2)$.

where $D(\Theta)$ is a matrix of quadratic function of some variable $\Theta$. The dimension of $\Theta $ is exactly $m$.

I have some heuristic, the best I can find is

Adam optimizer from tensorflow for first 3000 epoches

then use L-BFGS to fine tune local minimum.

I am looking for any suggestions!


问:
How to construct an effective preconditioner for this particular problem
A quick introduction to my problem

I am currently developing a method for simulation of water waves in three dimensions based on potential flow theory. The computational bottleneck of the method is to solve the equation
\begin{equation}
f(x,y) = \bigg( \cos\Big( \eta(x,y) \nabla \ \Big) - \sin \Big( \eta(x,y) \nabla \Big) \cdot \tan \Big( h \nabla \Big) \bigg) g(x,y)
\tag{1}
\label{equation1}
\end{equation}
for the function $g(x,y)$. Here $f(x,y)$ and $\eta(x,y)$ are known functions, $h$ is a positive constant, $\nabla = (\partial_x, \partial_y)$ and the operators  $\cos\Big( \eta(x,y) \nabla \ \Big)$,  $\sin \Big( \eta(x,y) \nabla \Big)$ and $\tan \Big( h \nabla \Big) $ are defined in terms of their Taylor series, for example
\begin{equation}
\cos \Big( \eta(x,y) \nabla \Big) = \sum_{m = 0}^{M} \frac{(-1)^m}{(2m)!} \eta(x,y)^{2m} \nabla^{2m}.
\end{equation}
I solve this equation by truncating the $\cos$ and $\sin$ operators at order $M$ and using the Fourier collocation method in which all functions are approximated by truncated Fourier series, for example
\begin{equation}
f(x,y) = \sum_{n_x = -N_x}^{Nx-1} \sum_{n_y = -N_y}^{Ny-1} \widehat{f}_{n_x, n_y} \exp \Big( i \, \mathbf{k}_{n_x, n_y} \cdot (x,y) \Big).
\end{equation}
Truncating the infinite operators at order $M$ gives an approximate equation for $g(x,y)$ which reads
\begin{equation}
\begin{aligned}
(1) \, \, \, f(x,y) = & \sum_{m = 0}^{\left \lfloor{M/2}\right \rfloor } \frac{1}{(2m)!} \eta(x,y)^{2m} \sum_{n_x = -N_x}^{Nx-1} \sum_{n_y = -N_y}^{Ny-1} | \mathbf{k}_{n_x, n_y} |^{2m} \widehat{g}_{n_x, n_y} \exp \Big( i \, \mathbf{k}_{n_x, n_y} \cdot (x,y) \Big) \\
& + \sum_{m = 0}^{\left \lfloor{(M-1)/2}\right \rfloor } \frac{1}{(2m+1)!} \eta(x,y)^{2m+1} \sum_{n_x = -N_x}^{Nx-1} \sum_{n_y = -N_y}^{Ny-1} | \mathbf{k}_{n_x, n_y} |^{2m+1} \tanh\Big( |\mathbf{k}_{n_x, n_y}| h \Big) \\
& \times \widehat{g}_{n_x, n_y} \exp \Big( i \, \mathbf{k}_{n_x, n_y} \cdot (x,y) \Big).
\end{aligned}
\end{equation}
Since the right hand side of this equation is easily evaluated when using the fast Fourier transform and because the matrix representation of the operator turning the values of $g(x,y)$ into the values of $f(x,y)$ is a full matrix, I solve this equation iteratively using the GMRES method. Unfortunately, this typically requires many iterations.

My question

I would like to ask for ideas and suggestions for how to construct a preconditioner for equation $\eqref{equation1}$ above.

What have I tried so far

So far I have tested the following ideas, but unfortunately none of them have worked:

Approximating the operator $\cos\Big( \eta(x,y) \nabla \ \Big) - \sin \Big( \eta(x,y) \nabla \Big) \cdot \tan \Big( h \nabla \Big)$ by keeping terms up to second order. This gives the operator $1 - (\eta(x,y)^2/2 + h \eta(x,y)) \nabla^2$  which I have discretized using the simplest possible finite difference method.

Not truncating the infinite operators gives rise to the equation
\begin{equation}
\begin{aligned}
f(x,y) = & \sum_{n_x = -N_x}^{Nx-1} \sum_{n_y = -N_y}^{Ny-1} \bigg( \cosh \Big( |\mathbf{k}_{n_x, n_y}| \eta(x,y) \Big) + \sinh \Big( |\mathbf{k}_{n_x, n_y}| \eta(x,y) \Big) \\
& \tanh\Big( |\mathbf{k}_{n_x, n_y}| h \Big) \bigg) \widehat{g}_{n_x, n_y} \exp \Big( i \, \mathbf{k}_{n_x, n_y} \cdot (x,y) \Big).
\end{aligned}
\end{equation}
I have tried replacing $\cosh \Big( |\mathbf{k}_{n_x, n_y}| \eta(x,y) \Big)$ and $\sinh \Big( |\mathbf{k}_{n_x, n_y}| \eta(x,y) \Big)$ by their average values over the computational domain in order to get something independent of $x$ and $y$ inside the summation such that the equation can easily be inverted using the fast Fourier transform.

Constructing the matrix representation of the operator to be inverted and ilu-decomposing it with different drop tolerances. For large drop tolerances this preconditioner does not help much, while for small drop tolerances this preconditioner is effective but this is because essentially the whole matrix is decomposed.

Performing a literature search. Among the different preconditioning strategies for potential flow water waves, I have not been able to find a strategy applicable to the equation I am trying to solve.


问:
Super Computer cluster service for students
I would like to ask whether there is a free of charge supercomputer for students, in my current institution they don't have a supercomputer, I am running MD simulations, I would be grateful for suggestions.</p>


问:
Constraint solver vs Bayesian optimizer for fast discontinuous processes
I have a complex domain-specific process that accepts inputs:

10-500 inputs, where each input is of type:

enum: choice between multiple string or numeric values

int: integers

float: floating point numbers

boolean: can be considered a enum(0,1)

It outputs several dozen outputs that I have clubbed together with weights to form an output singular composite loss function.

The process:

Is non-linear with an unknown polynomial degree

Is discontinuous

Can be evaluated quickly in parallel (<10 s per evaluation)

The third point above seems to make pure Bayesian Optimization a suboptimal choice for this problem, as per skopt:

If you do not have these constraints [f is expensive to evaluate], then there is certainly a better optimization algorithm than Bayesian optimization.

Should I be looking at constraint solvers, instead, since the inputs are well defined (if discontinuous)?

What is the difference between a constraint solver and a Bayesian optimizer in this case?

Details:

The real-world process is actually extremely time-consuming, but I have generated a decision-tree Regression model that evaluates predictions rather quickly.


问:
How to use Wolfe-Powell step-size control in quasi-Newton method?
I'm trying to find the minimum of a function using the quasi-Newton method with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.

But I want to change the following implementation, so that:

1) Wolfe-Powell step-size control is used to determine the minimum of a function using the quasi-Newton method.

2) Levenberg-Marquardt method is used to guarantee that the Hk (Hessian) stays positive definite.

How can I do that? What part of the code do I have to change? I don't get it.

Here is the working implementation of quasi-Newton with BFGS:

```<code>import numpy as np
import scipy.optimize, scipy.linalg
from scipy.optimize import minimize

##
## Finite difference quotient for vectors
##
def finite_differences(f, x, h = 1e-5):

    g = np.array([0.0]*len(x))
    x_tmp = np.array([0.0]*len(x))

    for i in range(len(x)):
        np.copyto(x_tmp, x)
        x_tmp[i] += h
        g[i] = ( f(x) - f(x_tmp) ) / h

    return g

##
## limit some quotient
##
def ensure_quotient(q, eps):
    if q >= 0:
        q = max(q, eps)
    else:
        q = min(q, -eps)
    return q

##
## DFP Update of the Inverse
##  M_{k+1} = M_k + \alpha_k \frac{d_k d_k^T}{\delta g_k^T d_k} - \frac{p_k \delta g_k^T}{p_k^T\delta g_k^T} M_k
##
def Update_Inverse_DFP(M_k, alpha_k, d_k, dg_k):
    ##
    ## get difference of search directions for conjugating
    ##
    p_k = M_k.dot(dg_k)

    ##
    ## get quotients and limit them so that there is no devision by zero
    ##
    q_1 = np.inner(d_k, dg_k)
    q_1 = ensure_quotient(q_1, 1e-8)
    q_2 = np.inner(p_k, dg_k)
    q_2 = ensure_quotient(q_2, 1e-8)

    ##
    ## calculate update
    ##
    M_new = M_k + alpha_k * np.outer(d_k, d_k) / q_1 - np.outer(p_k, dg_k).dot(M_k) / q_2

    return M_new

##
## BFGS Update of the Inverse
##
def Update_Inverse_BFGS(M_k, dx_k, dg_k):

    ##
    ## get quotients and limit them so that there is no devision by zero
    ##
    q = np.inner(dx_k, dg_k)
    q = ensure_quotient(q, 1e-98)

    ##
    ## get difference of search directions for conjugating
    ##
    U = np.outer(dx_k, dg_k)

    ##
    ## Identity
    ##
    assert U.shape[0] == U.shape[1], "Matrix dimension mixed up!"
    I = np.identity(np.shape(U)[0])

    ##
    ## calculate update
    ##
    M_new = ( I - U/q ).dot(M_k).dot( I - U.T/q ) + np.outer(dx_k, dx_k) / q

    return M_new

##
## perform a step of a quasi newton's method
##
def Quasi_Newton_Step(f, x_k, g_k, M_k, h = 1e-6, grad=lambda x: finite_differences(f, x), Update_Inverse = Update_Inverse_BFGS):

    ##
    ## determine direction
    ## d_k = - M_k g_k
    ##
    d_k = -M_k.dot(g_k)

    ##
    ## do golden rationline search to find an optimal alpha
    ## \alpha_k = argmin_\alpha \{  f(x_k + \alpha d_k)  \}
    ##
    line = lambda alpha: f( x_k + alpha * d_k )

    alpha_k = scipy.optimize.golden( line )

    ##
    ## get update x
    ## x_{k+1} = x_k + \alpha_k d_k
    ##
    x_new = x_k + alpha_k * d_k

    ##
    ## get gradient at new position
    ## g_{k+1} = \nabla f(x_{k+1})
    ##
    g_new = grad(x_new)

    ##
    ## get difference of gradients
    ##
    dg_k = g_new - g_k

    ##
    ## get difference of values
    ##
    dx_k = x_new - x_k

    ##
    ## get new M matrix
    ##
    M_new = Update_Inverse(M_k, dx_k, dg_k)

    return([x_new, g_new, M_new])

##
## Break condition for newton method
##
def Break_Condition(f, x_new, x_old, eps):
    return sum(abs(x_new-x_old))**2 < eps and abs(f(x_new) - f(x_old)) < eps

##
## Quasi Newton algorithm
##
## TODO: This implementation is very memory consuming because it keeps track of all iterates
## TODO: The precision is inherently dependent on the differencing step with h. What is the best value for that? eps?
def Quasi_Newton(f, x_0, max_iter = 100, h = 1e-6, eps = 1e-8, method = "BFGS"):

    if method == "BFGS":
        Update_Inverse = Update_Inverse_BFGS
    elif method == "DFP":
        Update_Inverse = Update_Inverse_DFP
    else:
        print('ERROR! Unknown Method: {}'.format(method))
        return ["-", "-"], 0
    ##
    ## store all intermediate values in lists
    ##
    x, g, M = [], [], []

    ##
    ## starting configuration
    ##
    g_0 = finite_differences(f, x_0)
    M_0 = np.identity(len(g_0))

    ##
    ## append to lists
    ##
    x.append(x_0)
    g.append(g_0)
    M.append(M_0)

    ##
    ## start iterations
    ##
    for k in range(max_iter):

        ##
        ## perform one quasi newton update step
        ##

        [ x_new, g_new, M_new ] = Quasi_Newton_Step(f, x[-1], g[-1], M[-1], h = h, Update_Inverse = Update_Inverse)

        ##
        ## append results of step to lists
        ##
        x.append(x_new)
        g.append(g_new)
        M.append(M_new)

        ##
        ## break condition
        ##
        if Break_Condition(f, x[-1], x[-2], eps):
            print('Converged after {} iterations'.format(k+1))
            break
    if not Break_Condition(f, x[-1], x[-2], eps):
        print('Diverged after {}!'.format(max_iter))

    return x[-1], f(x[-1])
</code>```

TEST the implementation:

```<code>def himmelblau(x, y):
    return (x**2 + y - 11)**2 + (x + y**2 - 7)**2

f = lambda x: himmelblau(x[0], x[1])

# Test with Himmelblau function

x_0 = np.array([6.0,0.0])

print("Starting Point:", x_0)
x_opt, f_opt = Quasi_Newton(himmel_blau, x_0, max_iter = 100, h = 1e-6, eps = 1e-6 )
print('Optimal value {} found at: {}'.format(f_opt, list(x_opt)))
</code>```


问:
Solved : Damped spring-mass system, wrong position, correct speed and acceleration
I am modulating a spring-mass system with gravitation and aero drag,  with python programming. The spring is hanging vertically and attached a weight. The user then selects a length to drag it down and the code will generate the graphs

I solved the differential equation this way.

\begin{align}
\sum F_x &= 0\\
\sum F_y &= -ky + mg - F_\text{drag}\\
\mathbf{F} &= m \mathbf{a}\\
m a_y &= -ky + mg + F_\text{drag}\\
m\frac{\mathrm{d}v}{\mathrm{d}t} &=  -ky + mg - F_\text{drag}\\
m \mathrm{d}v &= - ky \mathrm{d}t + mg \mathrm{d}t - F_\text{drag}\mathrm{d}t\\
m\Delta v &= -ky\Delta t + mg\Delta t - F_\text{drag}\Delta t\\
\Delta t &= 0.001 \text{ s}\\
\Delta v &= \frac{\sum F_y}{m}\\
v^{i + 1} &= v^{i} + \Delta v\\
\Delta s &= v \Delta t\\
x &= x + \Delta s
\end{align}

I solved this in python with this code:

```<code>from pylab import *
from matplotlib import pyplot as plt # kilde: https://matplotlib.org/users/pyplot_tutorial.html

# housekeeping
g = 9.81  # gravitation
m = float(input("how much do you want the mass to be? "))
l = 0.02     # aerodrag constant
k = float(input("insert springconstant:"))
#startvalues
t0 = 0.0  # start time
v0 = 0.0   # start speed
dv0 = 0.0#start axcelleration
xstart = float(input("how long do you want do drag it down? "))
x0 = (-m*g/k) - xstart  # start position

time = list() # list for time
time.append(t0) # insert start value for time
speed= list()
speed.append(v0)
position = list()
position.append(x0 + mg/k)
aksellerasjon = list() # accelleration
aksellerasjon.append(dv0)
dt= 0.001
x = x0
v = v0
t = t0

# main algorithm
while t < 20.0:

    aerodrag= l*v
    F = (-m*g - k*(x) - aerodrag)  #-m*g because gravitation is negative
    dv = F*dt/m
    v+= dv
    dx=  v*dt
    x += (dx)
    t += dt
    aksellerasjon.append(dv)
    time.append(t)
    speed.append(v)
    position.append(x)

# Plotting
plt.figure(figsize=(6,10))
subplot(3,1,1)
title("Akselerasjon")
xlabel("tid/s")
ylabel("akselerasjon m/s^2")
plot(tid, aksellerasjon, "b-")

subplot(3, 1, 2)
plot(tid, fart, "r-")
title("Hastighet")
xlabel("tid/s")
ylabel("fart m/s")

subplot(3,1,3)
plot(tid, posisjon, "g-")
title("Strekning")
xlabel("tid/s")
ylabel("posisjon/ m")
plt.tight_layout()

show()
</code>```

I have translated the comments and parts of the code from my language, so sorry for my bad English

The speed and acceleration are correct, but the position is wrong. It's too low, although the amplitude of the wave seams correctly.

I solved the differential equation in GeoGebra and the position should look like this.

What am I doing wrong?

After getting good help in the comments, i added + mg/k in the append

```<code>posisjon.append(x + (m*g)/k)
</code>```

The position now converges to zero, thanks!


问:
How to solve the Poisson equation with KINK aligned with mesh facet
I have a problem that solving the Poisson equation with kink ( discontinuous gradient but solution is continuous ) in the analytical solution, I want to solve this problem with FEM. To approximate this problem with continuous Galerkin, very fine resolution is required to capture the sharp gradient jump. With the attempt to solve it with less computation, I go through the discussion in the community and found some very useful inofrmation like

Why is my second order accurate method only converging at first order when the coefficients are rough?

and

Example of a continuous function that is difficult to approximate with polynomials

I alsoe tried out solving it with discontinuous Galerkin or add some penalty stabilization term on the standard Galerkin form, but still cannot get an idea approximation.

So my question is, what is the best method that I can rely on to solve this specific problem, that in this Poisson equation, the kink of the solution always aligned with the mesh facet.

Thanks a lot for your comments and help!


问:
Derivative of Whittaker-Shannon interpolant
Last time we looked at how to improve the accuracy of Whittaker-Shannon interpolation, where user njuffa demonstrated that judicious use of sin_pi could greatly improve the accuracy of computing
$$
f(t) = \sum_{j = 0}^{n-1} y_{j} \frac{\sin\left(\pi( \frac{t-t_0}{h} -j)\right)}{\pi\left(\frac{t-t_0}{h}-j\right)}
$$
Now I am struggling with accurate evaluation of
$$
f'(t) = \frac{\pi}{h} \sum_{j = 0}^{n-1} y_{j} \mathrm{sinc}'\left( \pi \frac{t-t_0}{h} -\pi j\right)
$$
My code for evaluating this derivative is here and for brevity the relevant portion is presented here:

   Real prime(Real t) const {
        using boost::math::constants::pi;
        using std::isfinite;
        using std::floor;

        Real x = (t - m_t0)/m_h;
        // The integer branch is not a problem:
        if (ceil(x) == x) {
            Real s = 0;
            long j = static_cast<long>(x);
            long n = m_y.size();
            for (long i = 0; i < n; ++i)
            {
                if (j - i != 0)
                {
                    s += m_y[i]/(j-i);
                }
                // else derivative of sinc at zero is zero.
            }
            if (j & 1) {
                s /= -m_h;
            } else {
                s /= m_h;
            }
            return s;
        }
        Real z = x;
        auto it = m_y.begin();
        Real cospix = boost::math::cos_pi(x);
        Real sinpix = boost::math::sin_pi(x);

        Real s = 0;
        auto end = m_y.end();
        while(it != end)
        {
            s += (*it++)*(pi<Real>()*z*cospix - sinpix)/(z*z);
            z -= 1;
        }

        return s/(pi<Real>()*m_h);
    }
```

I'm looking to achieve errors of roughly 10 ULPs, but this is giving errors of roughly ~1000 ULPs, which is a sizable fraction of the bit budget. I believe the source of the error comes from the (pi<Real>()*z*cospix - sinpix)/(z*z) term, but I haven't been able to rearrange it to get better accuracy.

Is there a way to get this more accurate?


问:
Techniques to optimise the integral of a function of known analytical form
I need to compute repeatedly a function that depends on an integral. The integral is not solvable analytically, but it depends on the argument of the function parametrically, like this:

$$
f(x) = \int_0^1g(x, u)du
$$

The function $g(x,u)$ is of known and fixed analytical shape, infinitely differentiable. My question is, is there any specialised algorithm that I can use to optimise this integral as much as possible, exploiting what I know of the function? I tried to see if it could be expressed as a Chebyshev series but no luck.


问:
how to estimate the number of people on a street within an hour?
I tried to use opencv to analyze a video filmed on the street.
But the problem is the performance is not enough. </p>

<p>I think  the number of people must follow the poisson distribution.
So I want to know that is there any statistical method to work out this question? </p>


问:
maximum eigenvalue across subsamples
I have an $N$-dimensional vector of data, say $X_{t}$, with $1 \leq t \leq T$.

Of this vector $X_{t}$, I want to consider vectors, say $X_{t}^{b}$, which are $m$-dimensional combinations of elements of the original vector $X_{t}$. In total, I have ${N}\choose{m}$ of such combinations.

Note that $N$ can be a large number; I need to allow for $N \rightarrow \infty$. Also, I will choose $m$ in such a way that $m$ is "smaller" than $N$, e.g. $m=O(N^{1/2})$.

I want to compute the $k$-th eigenvalue ($k$ is user-defined) of the second moment matrix of each of these $X_{t}^{b}$s, i.e.

$$\lambda_{k}^{b} \left( \frac{1}{T} \sum_{t=1}^{T} X_{t}^{b} (X_{t}^{b})' \right).$$

Then, finally, I want to compute

$$ \max_{1 \leq b \leq B} \lambda_{k}^{b}\, ,$$

where $B$ is defined as ${N}\choose{m}$. In principle, I may need to compute also other measures such as the average of $\lambda_{k}^{b}$.

I have two questions:

Is this problem NP-hard? Is there a reference to back up either statement (i.e. "it is" or "it is not")?

Is there some heuristic to make the problem less computationally burdensome? Is there any way to prove, for a given heuristic, that the solution found by the heuristic is "very close" to $ \max_{1 \leq b \leq B} \lambda_{k}^{b}$ - e.g. by showing that the solution found by the heuristic and $ \max_{1 \leq b \leq B} \lambda_{k}^{b}$ are equal almost surely, or something similar?

Many many thanks for your help.

Title: maximum eigenvalue across subsamples


问:
A preconditioner for self-consistent iteration
I tried to derive a preconditioner for self-consistent iteration similar
to section IX in arXiv:0804.2583.
For simplicity, consider here only
one orbital (one or two electrons) systems.

Suppose that the iteration starts with function $f_1$ and write
$f_1 = \psi_0 + \Delta \psi_0$, where $\psi_0$ is the normalized
exact solution of the Schrödinger equation. Suppose also that
$H$ is real and Hermitian, $H \psi_0 = E \psi_0$, and $\psi_0$ and
$\Delta \psi_0$ are real. Assume that $\Delta \psi_0$ is small.

Define $\Lambda := \left< f_1 \vert H f_1 \right>$ and
$g := H f_1 - \Lambda f_1$. We have
$$\Lambda = E + 2 \left< H \psi_0 \vert \Delta \psi_0 \right>
+ \left< \Delta \psi_0 \vert H \Delta \psi_0 \right>$$
and we approximate
$$
\Lambda \approx E + 2 \left< H \psi_0 \vert \Delta \psi_0 \right> .
$$
We have
$$
g = E \psi_0 + H \Delta \psi_0 - \Lambda \psi_0 - \Lambda \Delta \psi_0
$$
from which it follows that
$$
(H - \Lambda) \Delta \psi_0 = (\Lambda - E) \psi_0 + g
\approx 2 \left< H \psi_0 \vert \Delta \psi_0 \right> \psi_0 + g
\tag{1}
$$
We have
\begin{eqnarray}
2 \left< H \psi_0 \vert \Delta \psi_0 \right> \psi_0
& = & 2( f_1 - \Delta \psi_0) \left< H(f_1 - \Delta \psi_0) \vert \Delta
\psi_0 \right> \nonumber \\
& = & 2( f_1 - \Delta \psi_0)
\left( \left< H f_1 \vert \Delta \psi_0 \right> -
\left< H \Delta \psi_0 \vert \Delta \psi_0 \right> \right) \nonumber
\end{eqnarray}
and we approximate
$$
2 \left< H \psi_0 \vert \Delta \psi_0 \right> \psi_0
\approx
2 \left< H f_1 \vert \Delta \psi_0 \right> f_1
$$
We define linear function $G$ by
$$
G(x) := 2 \left< H f_1 \vert x \right> f_1 .
$$
Now it follows from equation (1) that
$$
(H - \Lambda - G) \Delta \psi_0 \approx g
$$
Consequently we approximate $\Delta \psi_0$ by
$$
\tilde{g} := (H - \Lambda - G)^{-1} g
$$
and set
$$
f_2 := \frac{f_1 - \tilde{g}}{\Vert f_1 - \tilde{g} \Vert}
$$
for the next iteration.
Unless we have an one electron system we also compute the new Hartree
potential and the new Hamiltonian operator for the next iteration.

I tried this method for the hydrogen atom and the hydrogen molecule.
For hydrogen atom I computed the approximate wavefunction with Arnoldi
method and then tried to enhance the result with this method. The
iteration converged well until both the gradient norm and quantity
$\Vert H \psi - E \psi \Vert$ were about $4 \cdot 10^{-15}$ and then
stopped converging.
For hydrogen molecule I computed the initial orbitals with Arnoldi method
setting Hartree potential to zero. The SC iteration converged well
until the gradient norm was about $10^{-4}$ and
$\Vert H \psi - E \psi \Vert \approx 10^{-6}$. After that the SC
iteration stopped converging.

Q: Is the derivation of the preconditioner correct?

Q: Is there some reason why the method does not work when the iteration gets close to a self-consistent solution?


问:
Kinetic preconditioning
Publication arXiv:0804.2583 describes a method for doing self-consistent iteration without having to diagonalize the Hamiltonian operator at every step.

IX. PRECONDITIONING

As already mentioned, direct minimisation of the total energy is used to find the converged wavefunctions. The gradient $g_i$ of the total energy with respect to the $i\text{-th}$ wavefunction $\left| \Psi_i \right>$ is given by
$$
\left| g_i \right>
= H \left| \Psi_i \right>
- \sum_{j}{\Lambda_{ij} \left| \Psi_j \right>} \,,
\tag{26}
$$
where $\Lambda_{ij} = \left< \psi_j \middle| H \middle| \psi_i \right>$ are the Lagrange multipliers enforcing the orthogonality constraints. Convergence is achieved when the average norm of the residue $\sideset{}{^{1/2}}{\langle \overline{g_i | g_i} \rangle}$ is below an user-defined numerical tolerance.

Given the gradient direction at each step, several algorithms can be used to improve convergence. In our
method we use either preconditioned steepest-descent algorithm or preconditioned DIIS method (27; 28). These
methods work very well to improve the convergence for non-zero gap systems if a good preconditioner is available.

The preconditioning gradient $\left| \tilde{g}_{i} \right>$ which approximately points in the direction of the minimum is obtained by solving the linear system of equations obtained by discretizing the equation
$$
\left(\frac{1}{2} \nabla^2 - \epsilon_i\right) \tilde{g}_i\left(\mathbf{r}\right)
= g_i \left(\mathbf{r}\right)
\,.
\tag{27}
$$
The values $\epsilon_i$ are approximate eigenvalues obtained by a subspace diagonalization in a minimal basis of atomic pseudopotential orbitals during the generation of the input guess. For isolated systems, the values of the $\epsilon_i$ for the occupied states are always negative, therefore the operator of Eq. (27) is positive definite.

Eq. (27) is solved by a preconditioned conjugate gradient (CG) method. The preconditioning is done by using
the diagonal elements of the matrix representing the operator $\frac{1}{2} \nabla^2 − \epsilon_i$ in a scaling function-wavelet basis. In the initial step we use $\ell$ resolution levels of wavelets where $\ell$ is typically 4. To do this we have to enlarge the domain where the scaling function part of the gradient is defined to a grid that is a multiple of $2^\ell .$ This means that the preconditioned gradient $\tilde{g}_i$ will also exist in a domain that is larger than the domain of the wavefunction $\Psi_i .$ Nevertheless this approach is useful since it allows us to obtain rapidly a preconditioned gradient that has the correct overall shape. In the following iterations of the conjugate gradient we use only one wavelet level in addition to the scaling functions for preconditioning. In this way we can do the preconditioning exactly in the domain of basis functions that are used to represent the wavefunctions (Eq. 5). A typical number of CG iterations necessary to
obtain a meaningful preconditioned gradient is 5.

–"Daubechies wavelets as a basis set for density functional pseudopotential calculations" (2008)

Question: How are the formulas (26) and (27) derived?

Title: Kinetic preconditioning


问:
Reverse automatic differentiation and integration
In Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control and more Sanz Serna writes:

It is well known that the reverse mode of differentiation implies
  an integration of the adjoint equations.

Any source that shows the fact explicitly?
If I want to differentiate $ z = f(x, y) $ w.r.t. $(x,y)$, what would be the adjoint equation that should be integrated?


问:
How to get the derivatives of the determinant and inverse of 2nd-order tensor wrt itself in SymPy?
I have a second-order tensor for which I need to compute the derivatives of its determinant and inverse w.r.t. itself. The equations are as follows:

$$\frac{\partial \, det(\mathbf{F})}{\partial F_{ij}} = det(\mathbf{F}) \, F^{-1}_{ji}$$

$$\frac{\partial \, F_{kl}^{-1}}{\partial F_{ij}} = - F^{-1}_{ki} \, F^{-1}_{jl}$$

These are the equations (49) and (60) in the matrix cookbook.

I have been working with SymPy for about a week. I can get the derivative of F wrt itself but can't figure out how to get the derivatives of its determinant and the inverse. For reference, I am posting the code and the output along with the error message.

Python code

```<code>from sympy import *

i = tensor.Idx('i',3)
j = tensor.Idx('j',3)
k = tensor.Idx('k',3)
l = tensor.Idx('l',3)

F = MatrixSymbol('F', 3, 3)

print("Derivative of F wrt F")
print("---------------------")
print(diff(F[k,l], F[i,j]))
print("\n")

J = det(F)
print("Derivative of det(F) wrt F")
print("--------------------------")
print(diff(J, F[i,j]))
print("\n")

print("Derivative of inv(F) wrt F")
print("--------------------------")
Finv = Inverse(F)
print(Finv[i,j])
print(diff(Finv[k,l], F[i,j]))
print("\n")
</code>```

Output

```<code>Derivative of F wrt F
---------------------------
KroneckerDelta(i, k)*KroneckerDelta(j, l)

Derivative of det(F) wrt F
-------------------------------
Derivative(Determinant(F), F[i, j])

Derivative of inv(F) wrt F
-------------------------------
Traceback (most recent call last):
  File "matdiffinverse.py", line 25, in <module>
    print(Finv[i,j])
  File "/usr/lib/python3/dist-   packages/sympy/matrices/expressions/matexpr.py", line 248, in __getitem__
return self._entry(i, j)
  File "/usr/lib/python3/dist-   packages/sympy/matrices/expressions/matpow.py", line 46, in _entry
    raise NotImplementedError(("(%d, %d) entry" % (int(i), int(j)))
  File "/usr/lib/python3/dist-packages/sympy/core/expr.py", line 207, in __int__
    raise TypeError("can't convert symbols to int")
TypeError: can't convert symbols to int
</code>```


问:
Is it Grid/Cluster/Cloud Computing? How are those terms defined?
There are three very connected and widely used terms:

Grid and grid computing

Cluster and cluster computing

Cloud and cloud computing

In many situations, it is not obvious which term to use, as I fail to find an easy to use, laconic, precise set of criteria which would separate those terms.

The reference I know about is

N. Sadashiv and S. N. Dilip Kumar, "Cluster, grid and cloud computing: A detailed comparison," in 6th Int. Conf. Comp. Sci. Education, Singapore, Aug. 2011, pp. 477–482. Also available here.

This reference (pretty well cited) provides detailed definitions and a comparison table; however, I find it hard to formulate "laymen's" definitions and criteria based on it. Probably, I am missing some very well-known publication that is used much more frequent.

I am looking for such simple definitions/criteria and reference(s) to ensure it's not an opinion-based question.


问:
Difference between contraction rate and convergence rate
I am slightly confused about the concept of contraction rate. For me it sounds equivalent to a convergence rate. Could somebody clear up the difference for me?
I have those two definitions of convergence and contraction rate from the book Numerical Optimal Control (M. Diehl):


问:
How to solve $y(x) y'''(x)=f(x)$
I have a PDE of the form $\partial_t y(x,t)+\partial_x(y(x) y'''(x)-f(x))=0$, where $f(x)=\cos(x)$.

Suppose a stable equilibrium exists, and I want to find the steady-state solution $y(x) y'''(x)=f(x)$ and I expect the physical solution to be periodic (as the source term is also periodic). How do I solve this equation (or 3-order nonlinear ODE in general) numerically (or analytically) in Scipy/Matlab/Mathematica?

For the boundary conditions, for the original problem I can just specify $y(x,0)$ and evolve the system, and the steady-state solution can be obtained for  large $t$. But if I want to solve the steady-state ODE directly, do I need to figure out three equations for the boundary conditions for it to be solvable?

Title: How to solve $y(x) y'''(x)=f(x)$


问:
Duplicate Surface Meshes in different Gmsh files
My objective is to get a conformal surface mesh between two geometries meshed independently in different files.
By example, I can start with a simple cube in one .geo file and mesh this geometry.

The script for this test1.geo file is as follows:

```<code>        // points - Block 1
       L = 1; H =1; Fin = 1/5; e =1; // L = length, H = width, e = elevation (z axis), Fin = mesh-size

    Point(1) = {0,0,0, Fin};
    Point(2) = {L/2,0,0, Fin};
    Point(3) = {L/2,H,0, Fin};
    Point(4) = {0,H,0, Fin};
    Point(5) = {0,0,e, Fin};
    Point(6) = {L/2,0,e, Fin};
    Point(7) = {L/2,H,e, Fin};
    Point(8) = {0,H,e, Fin};

Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,5};
Line(9) = {1,5};
Line(10) = {2,6};
Line(11) = {3,7};
Line(12) = {4,8};

// surfaces Block 1
Line Loop(1) = {-1,-2,-3,-4};
Plane Surface(1) = {1};
Line Loop(2) = {5,6,7,8};
Plane Surface(2) = {2};
Line Loop(3) = {1,10,-5,-9};
Plane Surface(3) = {3};
Line Loop(4) = {2,11,-6,-10};
Plane Surface(4) = {4};
Line Loop(5) = {3,12,-7,-11};
Plane Surface(5) = {5};
Line Loop(6) = {4,9,-8,-12};
Plane Surface(6) = {6};

// Volume Block 1
Surface Loop(1) = {4,1,3,2,5,6};
Volume(1) = {1};
</code>```

For various reasons, I need to write explicitly define all the lines and surfaces.

A part of the surface 4 has a physical interface with the second .geo file. i.e.

test2.geo

```<code>L = 1; H= 1; e = 1; Fin = 1/50;

Point(13) = {L/2,0,0.1*e, Fin};
Point(14) = {0.9*L,0,0.1*e, Fin};
Point(15) = {0.9*L,H,0.1*e, Fin};
Point(16) = {L/2,H,0.1*e, Fin};
Point(17) = {L/2,0,0.9*e, Fin};
Point(18) = {0.9*L,0,0.9*e, Fin};
Point(19) = {0.9*L,H,0.9*e, Fin};
Point(20) = {L/2,H,0.9*e, Fin};

Line(13) = {13,14};
Line(14) = {14,15};
Line(15) = {15,16};
Line(16) = {16,13};
Line(17) = {17,18};
Line(18) = {18,19};
Line(19) = {19,20};
Line(20) = {20,17};
Line(21) = {13,17};
Line(22) = {14,18};
Line(23) = {15,19};
Line(24) = {16,20};

Line Loop(13) = {-13,-14,-15,-16};
Plane Surface(13) = {13};
Line Loop(14) = {17,18,19,20};
Plane Surface(14) = {14};
Line Loop(15) = {13,22,-17,-21};
Plane Surface(15) = {15};
Line Loop(16) = {14,23,-18,-22};
Plane Surface(16) = {16};
Line Loop(17) = {15,24,-19,-23};
Plane Surface(17) = {17};
Line Loop(18) = {16,21,-20,-24};
Plane Surface(18) = {18};

// Volume Block 2
Surface Loop(13) = {16,13,15,14,17,18};
Volume(13) = {13};
</code>```

I want a conformal mesh on the mating physical interface between the .geo files
which is shown by the red arrow.

The previous posts I read had a complete surface mating so it was possible to duplicate the surface use the Geometry.CopyMeshingMethod = 1 option, however this is not the case that I have.
any suggestions, comments could be helpful!
Thanks a lot in advance,


问:
Isosurface extraction when you know the derivatives
I have a complicated equation that defines a shape in 3D, and I would like to generate a surface mesh. The shape is defined by an isosurface, i.e. the function is positive inside the shape and negative outside it.

The standard basic technique for this is the 'marching cubes' algorithm. However, it has a few disadvantages for me:

It generates meshes that are not particularly good, in ways that are well known.

At least with the most naïve implementation, I have to evaluate my function at every single point in a 3D grid, even though many of these are far away from the surface. I would rather have an algorithm that makes fewer evaluations of the function.

Marching cubes doesn't take advantage of any information about derivatives, which I can easily provide since I know the function.

Are there known, 'recommended' algorithms for this kind of case? There are many, many variations on marching cubes, but it's pretty hard to get an overview of their comparative advantages, and many of them seem only to address the first point above, not the second two.

This is for a little hobby project, so ideally I'm looking for something that's comparatively easy to implement.


问:
Regularized least squares with QR factorization
Consider the regularized least squares problem
$$
\min_x || b - A x ||^2 + \lambda^2 ||x||^2
$$
which is equivalent to
$$
\min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{A \\ \lambda I} x \right|\right|^2
$$
Suppose I have factored $A$ with a $QR$ factorization, so that $A = QR$ and therefore
\begin{align}
\min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{A \\ \lambda I} x \right|\right|^2 &=
\min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{QR \\ \lambda I} x \right|\right|^2 \\
&= \min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{Q & 0 \\ 0 & I} \pmatrix{R \\ \lambda I} x \right|\right|^2 \\
&= \min_x \left|\left| \pmatrix{Q^T b \\ 0} - \pmatrix{R \\ \lambda I} x \right|\right|^2
\end{align}
Now, I want to solve this regularized problem for many different $\lambda$ values, so my question is how to do this efficiently once I have already computed $R$ and $Q^T b$.

My thinking so far is to do another $QR$ factorization on the matrix
$$
\pmatrix{R \\ \lambda I} = W S
$$
The matrix is almost in upper triangular form except for the $\lambda I$ term, so I initially thought one could use Givens rotations to zero out the $\lambda I$ part, but it seems this will not work. For example, in the $N = 2$ case, we would have
$$
\pmatrix{ R_{11} & R_{12} \\ 0 & R_{22} \\ \lambda & 0 \\ 0 & \lambda}
$$
A first Givens transformation would zero out the $\lambda$ in the first column, however it would then destroy the sparse structure of the $\lambda I$ block. It would produce something like this:
$$
\pmatrix{ \tilde{R}_{11} & \tilde{R}_{12} \\ 0 & R_{22} \\ 0 & X \\ 0 & \lambda}
$$
where the $X$ term used to be zero but isn't anymore.

So a Givens approach seems inefficient since it does not preserve sparse structure. This seems to indicate that an optimized Householder approach is better since Householder transformations do preserve sparse structure.

Just wondering if anyone has any clever insights before I try coding up a Householder approach.

EDIT:

The book "Numerical Methods in Matrix Computations" by Ake Bjorck has an Exercise 2.6.2 in which he claims it is possible to perform a QR factorization of the matrix
$$
\pmatrix{B \\ \mu D}
$$
where $B$ is upper triangular and $D$ is diagonal, using Givens rotations in $O(11 n)$ flops. The only hint given (earlier in the chapter) is a reference to the following paper,

Eldén, L. (1984). An efficient algorithm for the regularization of ill-conditioned least squares problems with triangular Toeplitz matrix. SIAM journal on scientific and statistical computing, 5(1), 229-236.

This paper develops a Givens-based approach to QR factor the matrix
$$
\pmatrix{K \\ \mu L}
$$
where both $K$ and $L$ are upper triangular Toeplitz matrices. His algorithm is a straightforward application of Givens rotations, which work in this case because both $K$ and $L$ have the same non-zero pattern. When applying a Givens transformation to vectors $x$ and $y$, the result $(x',y')$ has a non-zero pattern which is the union of the patterns of $x$ and $y$.

Elden's algorithm for upper triangular Toeplitz matrices is $O(n^2)$ but I don't see an immediate way to adapt his method to the case of an upper triangular matrix and diagonal matrix.

Unfortunately Bjorck's exercise 2.6.2 doesn't give any further hints.


问:
Why I am not able to generate Levy distribution using inverse transform method on MATLAB?
I am trying to generate random numbers using Levy distribution in MATLAB. I am using the inverse transform method to generate random numbers and then plot its histogram for pdf estimation. Then I analytically plot the PDF of Levy distribution.

Problem: The estimated PDF and analytical PDF are not matching. I am putting my code here for reference where I have taken location parameter=5 and scale parameter=5.

```<code>clear all;
close all;
clc;
for loopr=1:1:10^8
X11(loopr)=5+5/(2*(erfcinv(rand))^2);
end
bin=[0:0.1:100];
l=0;
for loop1=5:0.1:100
l=l+1;
pdf12(l)=sqrt(5/(2*pi))*exp(-5/(2*(loop1-5)))/(loop1-5)^1.5;
end
loop1=5:0.1:100;
plot(loop1, pdf12)
hold on;
histogram(X11, bin, 'normalization', 'pdf');
hold off;
</code>```

The figure I am getting is:


问:
difference between PLSA and LDA
What is the difference between Probabilistic Latent Semantic Analysis(pLSA) and Latent Dirichlet Allocation (LDA)?</p>

Title: difference between PLSA and LDA


问:
How to perform local sensitivity analysis for partial differential equations
I am looking for a way to do local sensitivity analysis for PDEs, preferably in Python.

I get the impression that discretizing the equation then treating it as an ODE could work; however, would that be computationally tractable?

Wouldn't there be an efficiency gain in keeping the PDE structure rather than reverting to a (hard) ODE? Especially, allowing adaptive spatial grids (as in this paper that does adjoint sensitivity with adaptive mesh) seems like it would be a huge boost. Is this intuition correct? If yes, is there any prior existing implementation of PDE-specific sensitivity analysis, or should I tweak an ODE one or write from scratch?


问:
How to read data from simply formatted text file (c++)?
Mesh information like points, faces and cells is to be stored into separate files:

e.g. for points file:

# points data: x y z

N_Points 100

x1  y1  z1

x2  y2  z2

...

cell file:

# cells data: cell type and points index - clockwise (6: tri; 9:quad)

N_Cells 20

9 3 1 2 0

9 2 5 6 1

6 2 3 5

...

here # is comment, and N_Points  is the size of vector used to store these information.

I used to store these information directly to ascII files with just numbers, in which each line has same format.

For AscII file with comments, name-value pair and numbers, it is more readable, but I have no idea how to read different type of lines from txt file properly. I appreciate any tips to get this simple job done.


问:
Algorithm to join hexahedra and obtain outline volume
I would like to join several hexahedra and obtain an outline volume.

First, I started with 2D implementation. In 2D, there are non-intersecting quadrangles which always touch each other as shown in the figure. I used CGAL to join polygons and obtain the outline which is also a polygon.

Now I want to extend my code to 3D. In this case, instead of quadrangles, there  are hexahedra which also touch each other and do not intersect. CGAL, in this case, is not my first choice because it was a pain to achieve my goal in 2D when I used CGAL.

I would like to know if there is an intuitive algorithm to solve this problem without reading bunch of articles.


问:
OpenFOAM and CFDEM on GPU
I have a simulation project written for OpenFOAM and CFDEM and would like to find an alternative to run it on GPU since raising the number of cores already provided a promising speed up and extrapolation shows that GPU could be the way to go.

For OpenFOAM exists a project called "RapidCFD" which brings OpenFOAM computations to the GPU but it isn't compatible with CFDEM/LIGGGHTS and even porting it doesn't seem very promising referencing the answer of one of the developer (https://github.com/Atizar/RapidCFD-dev/issues/64).

My suspicion is that will not be easy or efficient to use the RapidCFD base code for a CFDEM solver.

LIGGGHTS is not GPU enabled, as far as I know.

The CFDEM libraries are not setup for GPU calculations.

RapidCFD is based on OpenFOAM 2.3.x, and presently CFDEM is based on    OpenFOAM 5.x (maybe the current version of CFDEM still works
  with OF
  2.3.x, or for sure you can use an earlier version of CFDEM, to get around this concern).

But the main idea of RapidCFD is to minimize data transfer between CPU
  memory and the GPU. Without porting the CFDEM libraries (or LIGGGHTS,
  too?) for GPU, I am not sure this will be as speedy as some of the
  other results people have discussed with RapidCFD, since data will
  need to be transferred between CPU and GPU (this is slow).

Finally, I don't think that all calculations lend themselves to
  speedup on the GPU, even if you were to port CFDEM libraries and/or
  LIGGGHTS to GPU. One would need to examine if the type of calculations
  done by CFDEM and/or LIGGGHTS are matrix type calculations that lend
  themselves to GPU speedup. If not, then after all the work is done,
  the speedup may not be so impressive.

Another project I found which could probably replace LIGGGHTS is EDEM (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6267572/) and due to the following publication (https://www.edemsimulation.com/news/new-cfd-coupling-interface-for-edem-enables-particle-fluid-simulation-with-any-cfd-software/) it should be possible to couple it with any CFD software.

As I have no experience in the area of CFD and DEM Particle simulations I would like to ask you if you know alternatives and possible expectations for a speed up using those tools. Does it make sense to switch to a commercial tool like EDEM and combine it with an GPU enabled CFD Solver.

Title: OpenFOAM and CFDEM on GPU


问:
Best way of porting code from the GPU to MPI-nodes
I have a program, structured in two parts, $A$ and $B$. Both parts are capable of running as standalone units, and written in C++. $A$ is written for cluster systems, running entirely on CPU-nodes, connected via MPI (both on the same CPU for multiple cores and for different nodes), and $B$ is targeted to run on a single GPU. $A$ is doing FEM-calculations, while $B$ is only doing matrix-vector-multiplications, with a static dense matrix and varying vectors (matrix-size is typically 5kx5k-15kx15k complex double elements). The result of the matrix-vector-multiplication must be accessible from all threads in $A$.

In order to keep the memory load on the GPU as low as possible my strategy until now is to create the static matrix only in the first MPI thread in $A$, but the result matrix in all threads (giving me a large data chunk in GPU memory for the first thread, and small data chunks for all other threads) by creating as many instances of $B$ as there are threads in $A$, but only generating the matrix for one instance. After doing the multiplications in thread 0, the result is distributed using MPI to all the other threads. $A$ does calculations based on that result, and afterwards restarts $B$ with different parameters.

Now I have to port $B$ to the same nodes as $A$, without a GPU. I was looking at PBLAS-functions (PZGEMV) or the PLASMA-library, and at the Trilinos-library (which I already use in $A$, thus integration in $B$ should be easy), but there I am not sure if dense distributed matrix-vector multiplications are available. What would be the best strategy for porting, and to make the resulting program as efficient as possible?


问:
Sequence planning with 3 machines
together!

First of all, I have to mention that because of my background as an Industrial Engineer, I have limited abilities in mathematics, but am disciplined enough to expand myself from competencies.
I am working in a production plant in the production planning. The manufacturing process is divided into 6 steps, but I only control the first. That means I am responsible for the order start. Our planning is currently based on Excel. Per product code, we check each week what the need is for the next two weeks. The possible capacities are then distributed according to the requirements.

There is a great challenge here.
We have three machines, which tend to be able to produce only a certain product portfolio. This means:
Machine 1 can produce 40% of the products.
Machine 2 can produce 50% of the products.
Machine 3 can produce 60% of the products.
This also means that we can not produce every product code every week and have to prioritize accordingly. This is done so far only after our feeling. Here we would like to create a rule.

Actually, the product portfolio or the variance results from two parameters:
The thickness of the product itself, as well as the size of the hole. The machines are drills. They grasp the product with pliers and drill it.
Not every machine can reach any parameter. That means there are the first ruins here.
The products are already ordered in classes, which are assigned according to the thickness. Here we can roughly say what needs to be served per class every two months.

How can an optimal sequence planning be determined for three machines, each of which can only produce a certain product portfolio, and the individual needs differ per product class.

I am grateful for every tip, approach, concept. Like to name just a few keywords or literature, so I can deepen myself.
So far I could only identify the simplex method and the Johnson algorithm as a potential approach. If you need further information, please let me know!

greeting
Joshua

Title: Sequence planning with 3 machines