问:
Jacobian of the electron and hole drift diffusion equations with respect to potential in semiconductors
These are the discretized drift diffusion equations as taken from the book "Analysis and Simulation of Semiconductor Devices". The electron continuity equation is:

$$((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^n.B((\varphi_{i+1,j}-\varphi_{i,j})/V_T).n_{i+1,j}  + ((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^n.B((\varphi_{i-1,j}-\varphi_{i,j})/V_T).n_{i-1,j}  - [(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^n.B((\varphi_{i,j}-\varphi_{i+1,j})/V_T)) + (((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}).D_{i-1/2,j}^n.B((\varphi_{i,j}-\varphi_{i-1,j})/V_T))  + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^n.B((\varphi_{i,j}-\varphi_{i,j+1})/V_T)) +  (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^n.B((\varphi_{i,j}-\varphi_{i,j-1})/V_T))].n_{i,j}  + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^n.B((\varphi_{i,j+1}-\varphi_{i,j})/V_T).n_{i,j+1}  + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^n.B((\varphi_{i,j-1}-\varphi_{i,j})/V_T).n_{i,j-1}  = R_{i,j}^{EFFECTIVE}.((\Delta_{i-1}^x+\Delta_i^x)/2).((\Delta_{j-1}^y+\Delta_j^y)/2)$$

The hole continuity equation is:

$$((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^p.B((\varphi_{i,j}-\varphi_{i+1,j})/V_T).p_{i+1,j}  + ((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^p.B((\varphi_{i,j}-\varphi_{i-1,j})/V_T).p_{i-1,j}  - [(((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_i^x).D_{i+1/2,j}^p.B((\varphi_{i+1,j}-\varphi_{i,j})/V_T)) + (((\Delta_{j-1}^y+\Delta_j^y)/2\Delta_{i-1}^x).D_{i-1/2,j}^p.B((\varphi_{i-1,j}-\varphi_{i,j})/V_T))  + (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^p.B((\varphi_{i,j+1}-\varphi_{i,j})/V_T)) +  (((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^p.B((\varphi_{i,j-1}-\varphi_{i,j})/V_T))].p_{i,j}  + \\ ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_j^y).D_{i,j+1/2}^p.B((\varphi_{i,j}-\varphi_{i,j+1})/V_T).p_{i,j+1}  + ((\Delta_{i-1}^x+\Delta_i^x)/2\Delta_{j-1}^y).D_{i,j-1/2}^p.B((\varphi_{i,j}-\varphi_{i,j-1})/V_T).p_{i,j-1}  = R_{i,j}^{EFFECTIVE}.((\Delta_{i-1}^x+\Delta_i^x)/2).((\Delta_{j-1}^y+\Delta_j^y)/2\ ,$$

where $\Delta_i$ and $\Delta_j$ are the mesh spacing for finite difference grid along X and Y directions. $D_n$,$D_p$ are electron and hole diffusion coefficients.$B(x) = x/(e^x-1)$ is Bernoulli's function and $V_T$ is the thermal voltage.

Also, $$R^{EFFECTIVE} = R^{SRH} + R^{AU} + R^{OPT} + R^{II}$$ where
$R^{EFFECTIVE}$ is the effective generation/recombination term. $R^{SRH}$, $R^{AU}$, $R^{OPT}$ and $R^{II}$ are SRH, Auger, Radiative and Impact Ionization Recombination terms.

Now, I am applying the coupled Newton's iterative scheme to solve the above equations where I need to differentiate the equations. When differentiating with respect to potential '$\varphi$', I noticed that $\varphi$ appears in three places namely,

High field Electron/Hole Mobility, $\mu_{n,p}^{high}=\mu_{n,p}^{low}⁄ [1 + (\mu_{n,p}^{LICN}.E_{n,p}/v_{n,p}^{sat})^{\beta_{n,p}}]^{1/\beta_{n,p}}$ and $D_{n,p} = \mu_{n,p}.V_T$;

Bernoulli's function, $B(x)$ where 'x' here means difference between potentials at neighboring grid points and

Impact Ionization $R^{II} = -\alpha_n.(|J_{i,j}^n|/q) - \alpha_p.(|J_{i,j}^p|/q)$, $\alpha_n = a_n.e^{-b_n/|E_{i,j}|}$, $\alpha_p = a_p.e^{-b_p/|E_{i,j}|}$ where $J^{n,p}$ and $E$ represent electron/hole currents and electric field respectively.
Current terms are made of potential terms as in,
$J_{i+1/2,j}^{nx}=D_{i+1/2,j}^n.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).n_{i,j} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).n_{i+1,j}]/\Delta_i^x  \rightarrow  (31)$, $J_{i,j+1/2}^{ny}=D_{i,j+1/2}^n.[B(\frac{\varphi_{i,j}-\varphi_{i,j+1}}{V_T}).n_{i,j} - B(\frac{\varphi_{i,j+1}-\varphi_{i,j}}{V_T}).n_{i,j+1}]/\Delta_j^y  \rightarrow  (32)$, $J_{i+1/2,j}^{px}=D_{i+1/2,j}^p.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).p_{i+1,j} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).p_{i,j}]/\Delta_i^x  \rightarrow  (33)$, $J_{i,j+1/2}^{py}=D_{i,j+1/2}^p.[B(\frac{\varphi_{i,j}-\varphi_{i+1,j}}{V_T}).p_{i,j+1} - B(\frac{\varphi_{i+1,j}-\varphi_{i,j}}{V_T}).p_{i,j}]/\Delta_j^y  \rightarrow  (34)$

So, basically I need to differentiate all of these to find the Jacobian with respect to $\varphi$. But, the problem arises when the difference between potentials become zero. This means that the entire electron/hole equation becomes a constant when the potential difference is zero. Hence, the Jacobian is also zero.

This would mean that the final Jacobian matrix would have entire columns and rows made out of zeros thus making the matrix singular and it is not possible to solve a singular matrix.

Is this the correct approach to solving using the Newton's method. What should I do when the potential difference is zero?

References

Selberherr, Siegfried. Analysis and simulation of semiconductor devices. Springer Science & Business Media, 2012.


问:
Update model parameter with new data, discarding old data
I have this dataset, and I am using y = (a * x^n) / (b + x^n) Hill function as the model, where a is the limit of the Hill curve, b is the point at which a/2 is reached (for n = 1) and n is the cooperativity or steepness of the curve.

Currently, I am storing all X,y values, computing the parameters from scipy.optimize.curve_fit, and plotting the curve. If new data points come along, I re-calculate the parameters with the old+new data.

Is there a way to update the parameters of the model without storing all of the previous old data points, once the initial parameters are obtained from the previous data points?

Example, I fit the curve to the first 1000 data points and have my parameters. Next, I discard some or all of the old data. Then, when I see the 1001st point I simply update my parameters and plot the curve again and so on for every new data point.

EDIT

My existing code is as follows (not super elegant).

```<code>import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def file_stream(file_name):
    with open(file_name, 'r') as in_file:
        for line in in_file:
            yield map(float, line.strip().split('\t'))

def hill_model(X, a, b, n):
    return [float((a * x**n)) / (b + x**n) for x in X]

def get_params(X_all, y_all, prev_par=None, fn=hill_model):
    if prev_par is None:
        a_init, b_init, n_init = y_all[-1], y_all[0], 1.0
    else:
        a_init, b_init, n_init = prev_par
    opt_par, opt_cov = curve_fit(fn, X_all, y_all, p0=[a_init, b_init, n_init])
    a_final, b_final, n_final = opt_par
    return a_final, b_final, n_final

def main():
    file_name     = 'data.tsv'
    file_streamer = file_stream(file_name)
    X_all, y_all  = [], []

    # Get some intial data from stream
    for _ in xrange(1000):
        X, y = file_streamer.next()
        X_all.append(X)
        y_all.append(y)
    plt.scatter(X_all, y_all)

    # Initialize params of model
    a, b, n = get_params(X_all, y_all)
    y_model = hill_model(X_all, a, b, n)
    plt.plot(X_all, y_model, 'r-')
    plt.show()

    # Rolling update
    seen_all = False # Helps stop when all data is fit
    while True:
        for _ in xrange(1000):
            try:
                X, y = file_streamer.next()
                X_all.append(X)
                y_all.append(y)
            except:
                seen_all = True
                break
        a, b, n = get_params(X_all, y_all, prev_par=[a, b, n], fn=hill_model)
        y_model = hill_model(X_all, a, b, n)
        plt.scatter(X_all, y_all)
        plt.plot(X_all, y_model, 'r-')
        plt.show()

        # Nothing more to update, return
        if seen_all:
            return

if __name__ == '__main__':
    main()
</code>```

The code currently reads in some X,y values, calculates the a, b, n parameters and when more X,y values are added, the code updates a, b, and n params. As you can see, I need to store previous X,y values, which I do not want. I want to update the parameters as new X,y values are seen and from the previous a, b, and n values only.

ContentLicense: CC BY-SA 3.0


问:
Computing the change of function at two close points without cancellation
I want to compute the difference $\Delta f(x_1,x_2) = f(x_1)-f(x_2)$ of a smooth function $f(x)$ at two points $x_1$ and $x_2$ which are close to each other. The magnitude of the expected result, $|\Delta f(x_1,x_2)|$, is often small compared to $|f(x_1)|$ and $|f(x_2)|$ and several significant digits are lost if I evaluate the expression $f(x_1)-f(x_2)$ straightforwardly. What is a recommended way to compute $\Delta f(x_1,x_2)$ accurately?

I have an analytical expression of $f(x)$. If a detail matters, it has a form like
$f(x) = \exp(-a x^2) \sin(x)$.

Perhaps one method would be to use the Taylor expansion of $f(x)$, but I am afraid that the convergence might not be very fast. Also, writing down the higher order derivatives is tedious. I'd like to know something easier and more efficient.


问:
Precision not improving by decreasing step-size in nonlinear Schrödinger
I tried to simulate soliton propagation by solving the nonlinear Schrödinger equation using the split-step Fourier method. The following is an example of the Matlab code copied from a textbook.

```<code>alpha = 0;
beta_2 = 1;
gamma = 1;

A = N*sech(T);   % initial pulse shape
L = max(size(A));
delta_omega = 1/L/delta_T*2*pi;
omega = (-L/2:1:L/2-1)*delta_omega;

solution = A;   % analytical solution (same as original it is periodic)

step_num_hist = round(10.^linspace(1,3.5,20));   % total number of steps
avg_err = zeros(3,max(size(step_num_hist)));

end_pt = 0.5*pi;   % the end of simulation, one soliton period

for step_num = step_num_hist
    h = end_pt/step_num;   % step size
    A_t = A;
    for n=1:step_num
        A_f = fftshift(fft(A_t));
        A_f = A_f.*exp(-alpha*(h/2)-1i*beta_2/2*omega.^2*(h/2));
        A_t = ifft(ifftshift(A_f));

        A_t = A_t.*exp(1i*gamma*(abs(A_t).^2*h));

        A_f = fftshift(fft(A_t));
        A_f = A_f.*exp(-alpha*(h/2)-1i*beta_2/2*omega.^2*(h/2));
        A_t = ifft(ifftshift(A_f));
    end
    comp_result = abs(A_t);
    avg_err(1,m) = sum(abs(comp_result.^2 - abs(solution).^2))/max(size(solution))/max(abs(solution).^2);   % avg relative intensity error, found in a paper
    m = m+1
end
</code>```

The code seems to produce something that makes sense, so I don't think it is wrong. The problem is that when the step number increases (step size h decreases) beyond a certain number (in my case around 1e-3), the error between the computed solution and the analytical solution doesn't decrease any further. This happens in both Matlab and Python. I also tried implementing a few higher-order methods that I found in the literature. They do converge faster, but the 'average relative intensity error' does not go below somewhere about 1e-3, while a paper I found which compares various algorithms can go down to like 1e-7 or even 1e-10 (total number of steps also of the order of 10e4).

So my question is why is that? Why does the error decrease with step size initially but then stops decreasing after the step size gets to a certain point? What can possibly go wrong?


问:
Minimize number of rectangles that cover all the points
I have a 2d distribution of moving points with known trajectories represented in a 640x480 image. Here is the initial state:

I have to find the minimum number of rectangles with fixed dimensions (35x42 pixels) that are able to cover all my input data, minimizing also the overlapping areas for each step.

So far I tried to solve the problem using the hierarchical clustering algorithm cutting the dendogram when observations in each flat cluster have no greater a cophenetic distance than $t$. This approach works pretty well when points are quite scattered but, little by little they get closer and the solution gets worse because the number of overlapping rectangles increases:

.

I was thinking to use an optimization approach but I don't know how to model the system. Do you have some idea how to do this? Or, do you know other approaches to face this problem? Any help is greatly appreciated!


问:
Parallel dense solve with submatrices from mesh refinement with Petsc
For a Bounday Element Method problem I require the solution of a system of linear equations with multiple right-hand sides. Though this is a dense system, I still want to do it via Petsc in parallel. Would the best way to do this be through something such as -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps?

Furthermore and importantly, I refine my grid in a regular way by adding $N-1$ points between the already present $N$ points. This leads to a new system of equations of size $N\left(N-1\right) \times N\left(N-1\right)$ in $N\left(N-1\right)$ unknowns.

The upper-left matrix $A_{11}$ in this new system is identical to the unrefined matrix (multiplied by 0.5). The other three blocks $A_{12}$, $A_{21}$ and $A_{22}$ are new:

$$\begin{bmatrix}A_{11}&A_{12}\\A_{21}&A_{22}\end{bmatrix}$$

The question is now whether knowledge about the unrefined matrix $A_{11}$ can be used to speed up calculation of the refined system?

I was thinking, for example, about using the LU decomposition of $A_{11}$ to calculate the LU decomposition of the entire matrix $A$ using the well-known formula's

$$
\left\{\begin{aligned}
L_{21} U_{11} &= A_{21} \\
L_{11} U_{12} &= A_{12} \\
L_{22} U_{22} &= A_{22} - L_{21} U_{12}
\end{aligned}\right.
$$

where $L_{21}$ and $U_{12}^T$ are $N-1\times N$ full matrices, $L_{11}$, $L_{22}$ are lower triangular $N\times N$ and $N-1 \times N-1$ matrices and $U_{11}$, $U_{22}$ are upper triangular with the same sizes.

Is there any way to do this in Petsc?

Or is there a better thing I can do?

For example, due to the fact that $A_{21}$, $A_{12}$ and $A_{22}$ result from mesh refinement of the original system, they all look a lot like $A_{11}$: Definining

$$ V = \left.\underbrace{\left[\begin{array}{cccc}
0.5 & 0.5 & 0 & \cdots & 0 & 0 \\
0 & 0.5 & 0.5 & \cdots & 0 & 0\\
\vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 0.5 & 0.5
\end{array}\right]}_{N}\right\} N-1
$$
an $N-1 \times N$ matrix, it can be seen that

$A_{21} \approx V A_{11} $, i.e. an average of the rows of $A_{11}$.

$A_{21} \approx A_{11} V^T $, i.e. an average of the columns of $A_{11}$.

$A_{22} \approx V A_{11} V^T $, i.e. an average of the rows and columns of $A_{11}$.

This would therefore mean that a good approximate for $A$ would be

$$ A \approx \left(\begin{array}{c} 1_N \\ V \end{array}\right) A_{11} \left(\begin{array}{c} 1_N \\ V \end{array}\right)^T , $$

from which it might be possible to construct some preconditioner easily using the LU decomposition from $A_{11}$:

$$
\left\{\begin{aligned}
L_{21} &\approx V L_{11} \\
U_{12} &\approx U_{11} V^T
\end{aligned}\right. ,
$$
from the equations above. However, this implies that $L_{22} \approx 0$ and $U_{22}\approx 0$, which I do not really understand well as it leads to non-invertible approximations for $U$ and $L$:

$$ L \approx \left(\begin{array}{cc} L_{11} & 0 \\ V L_{11} & 0 \end{array}\right) \ \text{and} \ U \approx \left(\begin{array}{cc} U_{11} & U_{11} V^T \\ 0 & 0 \end{array}\right) $$

How would this be used as a preconditioner?


问:
FEM Stiffness and Mass matrices for 2d cubic trial functions
I want to use FEM to solve a 2D eigenvalue problem for the biharmonic equation

\begin{align}
&K \psi= \lambda M \psi\\
&\psi=\nabla \psi\cdot \hat{n}=0\quad \text{on boundaries}
\end{align}

where $K$ is the stiffness matrix corresponding to

$$K_{ij} = \int_{A} \nabla^2 \psi_i  \nabla^2 \psi_j \mathrm{d}A$$

and $M$ is the mass matrix

$$M_{ij} = \int_{A} \psi_i \psi_j \mathrm{d}A$$

Because the equation is fourth order I am using cubic trial functions, both height and slope.

I know how to get the stiffness and mass matrices in 1D, but I am stumped by the 2d problem.

Does anyone know of a suitable reference?


问:
Calculate activation energy of inter molecular displacement on a nano particle
I'm simulating a nano particle and wish to calculate the activation energy of the reaction where an atom is displaced from one point on the lattice to a neighboring void( might be at higher energy than initial).
Calculating this activation energy using energy difference will require me to calculate the complete trajectory to determine the point with highest energy, which is computationally very expensive.
Is there any other tested method to find this activation energy without having to perform so many computations each time.
Note: the new site is a neighbor of the old site and I expect( not assume) the atom to roll over two of its adjacent atoms to reach this new site.
Also all the atoms are identical.</p>


问:
How to measure euclidean distance between points with vtkDelaunay3D package?
I'm working with python vtkDelaunay3D package for a special purpose. I have set some points but I do not know how can I measure euclidean distance between points with this package. Is there anyone who knows measuring euclidean distance between points with this package in python?</p>


问:
One Dimensional Schrodinger's Equation solution using Numerov Method
I have been trying to solve Time Independent Schrodinger's equation in one dimension using Numerov Method as discussed in this excellent lecture notes I found on net.

The Numerov method can solve an equation of the following kind:

$$\frac{{d^2}y}{dx^2}=-g(x) y(x) +s(x) $$

We can compare this with out Time Independent Schrodinger Equation :

$$\frac{{d^2}\psi}{dx^2}=-\frac{2m}{\bar h^2} \big(E-V) \psi$$

Where $ g(x)= -\frac{2m}{\bar h^2}(E-V) $

We can discretize the equation using Taylor expansion and write the iterative formula as :

$$\psi_{i+1}= \frac{\psi_i\big[2-\frac{5}{6}\cdot h^2 g_i]-\psi_{i-1}\big[1+\frac{h^2}{12}\cdot g_{i-1}]      }{(1+\frac{h^2}{12}\cdot g_i)}$$

Using the known Energy of Harmonic Oscillator I solved the equation and my solutions were working well.

I also used the shooting method to determine Energy value for a particular Eigen State and got a great result.

I was succesfully able to solve the schrodinger equation for Most of Even potential including Morse Potential and The Double Well Potential.

The Problem arrises if I try to take a non symmetric solution. Because I have taken the symmetry of my wave function into account and just solved for $\psi$ from $x=0 $ to $ x=xmax$, and then took output in a file by multiplying a parity to the $\psi$ based on weather my function is odd or not.

But When I solved the same problem from $x=-xmax $ to $x=xmax$ using same algorithm so that I could devise a general algorithm for Numerov Method. The output started getting bizarre. While the ground state solution is same , the odd state solutions are getting inverted somehow.

This is a part of my program for the first (working) Numerov Algorithm. (ignore the signchange counter I have used it to determine number of nodes for shooting method)

```<code>cout<<"Enter mesh number"<<endl;
 cin>>n;

 cout<<"Enter maximum value of x"<<endl;
 cin>>xmax;
 dx=xmax/n;

 if (nodes % 2==0) //Even Nodes, function is even function , i.e. y(-1) =      y(1) and f(-1) = f (1); Thus using algorithm!
        {
      y[0]=1.0;
      y[1]=( y[0]*( 12.0 - 10.0 *f[0] ) )/( 2*f[1] );
    }
  else //nodes are odd, i.e there is a node at x=0;
    {
      y[0]=0;
      y[1]=dx;           //arbitrary small value
    }

    signchange=0;

    //outward integration using algorithm and counting number of sign changes//

 for(int i=1; i <= n; i++)

    {
      y[ i+1 ]= ( ( 12.0 - 10.0*f[i] )*y[i] -( y[i-1]*f[i-1] ) ) / f[i+1];

     if (y[i] != copysign(y[i],y[i+1]))
       ++signchange;
    }
//For Output :

if (nodes % 2 == 0)
    p=1;
    else
    p=-1;
    for (int i = n; i > 0; --i)
    file<<setw(1)<<-1*x[i]<<setw(15)<<p*y[i]<<setw(15)<<y[i]*y[i]<<setw(15)<<V[i]<<endl;

    //-------x<0-------//
    for (int i = 0; i <= n; ++i)
    file<<setw(1)<<x[i]<<setw(15)<<y[i]<<setw(15)<<y[i]*y[i]<<setw(15)<<V[i]<<endl;
</code>```

And this is part of the second program which is giving me inverted output. Where I am iterating across the whole mesh

```<code>cout<<"Enter maximum value of x"<<endl;
cin>>xmax;

xmin=-1*xmax;

h=(2*xmax)/n;

y[0]=0;
    y[1]=h;

    /*calculating the wave function psi or y at all points for Energy e*/
        for(int i=1; i <= n; i++)

    {
      y[ i+1 ]= ( ( 12.0 - 10.0*f[i] )*y[i] -( y[i-1]*f[i-1] ) ) / f[i+1];

     if (y[i] != copysign(y[i],y[i+1]))
       ++signcount;
    }
</code>```

I would love to know why this is happening. It does not make sense to me. But i am also a beginner in Computational Physics. Also I would like to know what should be the best way to solve the equation for Non-Symmetric potentials, for example a Coulumb potential (inverse R).

EDIT: Since I am trying to solve Harmonic Oscillator problem the wave function should be bound on both sides. So i am taking $\psi_0$ as 0. And $\psi_1$ as some random small value which will be taken in to account when the wave function will be normalized later. This is for the case when I am solving the equatioin across the whole mesh, from $-xmax$ to $xmax$ .

But When I am solving the equation for symetric potential(i.e. from 0 to xmax), I am using a property of solution that I have already known, that if i am solving for odd Energy state, the wave funcion $\psi_{x=0}$ will always be zero. I have provided my declaration of initial value in code too. I hope this helps, in clarifying my problem a little.


问:
shallow water equation maccormack method
I am trying to make a code for 1D shallow water equation (nonlinear without source terms) using the MacCormack method for sinusoidal wave propagation. My issue is that the wave fluctuates and does not produce a smooth result. I think there might be an issue with my boundary conditions as I am giving boundary condition in linear terms:

\begin{align}
h=a\sin\left(\frac{2\pi n}{T}\right)\\
u=\left(\frac{g}{h_0}\right)^{1/2}h
\end{align}

where $n$ is temporal, $a$ is $0.1$, and $h_0$ is $1$ meter. $\Delta t$ is $0.01$ and $\Delta x$ is $0.1$. If I understood my problem correctly, please tell me how I can give boundary conditions for the nonlinear terms of shallow water equation.


问:
Numerical Precision in Matrix Inversion Routines
Let's say I have a matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$, and its associated inverse, $\mathbf{A}^{-1}$. The elements of $\mathbf{A}$ are given in IEEE single precision, i.e. 23-bit mantissa, and I use the lapack routines $\texttt{sgetrf}$ and $\texttt{sgetri}$ to compute $\mathbf{A}^{-1}$ through a single-precision LU-decomposition.

How am I to ascertain the error associated with the matrix product $\mathbf{A}\mathbf{A}^{-1}$. More concretely, I'd quite like to be able to approximate, a priori, what the maximum error, $\varepsilon = \max_i \Big|\big( \mathbf{AA}^{-1}\big)_{ii} - 1.0\Big|$, will be. Furthermore, is there a way to approximate the errors for other matrix inversion routines, such as Gauss-Jordan or a QR decomposition? I'm guessing I should use their computational scaling, e.g. the LU-decomposition scales as $\frac{4}{3}n^3$, but I'm unsure how to proceed beyond that.


问:
Does excitation type matter in a time-domain simulation/computation of a transfer function of an LTI system?
Let's say I am running a FDTD simulation of a wave-equation to determine a transfer function of an LTI system:

\begin{equation}
H(f) = \dfrac{Y(f)}{X(f)}\
\end{equation}

where $Y(f)$ and $X(f)$ are Fourier transforms of signals measured at some point in the FDTD grid for a grid with an "obstacle" inserted and an empty grid, respectively, due to an excitation $s(t)$. Assume that it's an acoustic wave simulation and that the "obstacle" is modelling a fully reflective object.

My question is: does the type of an excitation signal $s(t)$ matter in that kind of problem? For example, how broadband the signal is? And if does matter, then - how? My intuition tells me since in the end we're computing a relation of two signals, any propagation issues of the input signal will be the same in both simulations and therefore cancel each other in the above equation.


问:
Solution to 1D consolidation problem python implementation
A solution to the 1D consolidation problem is given by

$$\frac{\partial}{\partial t} p = c_{v} \frac{\partial^{2}}{\partial y^{2}}  p$$

where $p$ is the pore water pressure, $c_v$ is the consolidation coefficient, $t$ is the elapsed time and $y$ is the height of the column. Further background is provided by this document.

The solution is reproduced below:

$$p(y,t) = \sum_{n=1}^{n = \infty}\frac{2p_0}{n\pi}(1-\cos(n\pi))\sin(\frac{n \pi y}{2h}\exp(\frac{-n^2\pi^2c_v t}{4h^2}))$$

with initial conditions

$$p(y,0) = p_0$$

and boundary conditions

$$p(0,t) = 0$$
$$p(2h,t) = 0$$

To my knowledge that initial condition is a straight line equivalent to the initial loading applied on the soil layer. At $t = 0$ the graph of $y$ versus $p(y,t)$ should look like this for $p_0 = 5000000$

Here is my python implementation of the above solution.

```<code>from numpy import cos, exp,pi, linspace, array,sin
import matplotlib.pyplot as plt

#PRESSURE VERSUS HEIGHT

G = 5000000 # Pa = N/m2
E = 10000000# Pa = N/m2
K = E*G/3/(3*G-E)#Bulk modulus of soil skeleton Pa = N/m2
Kl = 2000000000 #Bulk modulus of liquid Pa = N/m2
n = 0.33 #porosity
gammal = 9810 #unit weight of liquid N/m3
k = 1.157E-17 #intrinsic permeability liquid m2
h = 0.5 #m

alpha = 1/(K+4*G/3) #m*s2/kg
beta = 1/Kl #m*s2/kg

cv = k/(gammal*(alpha+n*beta))#m^2/s

n = 200#number of terms to include in series
timevalues = array([0])
yvalues = linspace(0,1,100)
h = 0.5
p0 = 5000000

for t in timevalues:
    pvalues = []
    for y in yvalues:
        SUM = 0
        SUM2 = 0
        for j in range(1,n,1):
            SUM = SUM+(2*p0/j/pi)*(1-cos(j*pi))*sin(j*pi*y/2/h)*exp(-j**2*pi**2**(cv*t/h**2)/4)
        pvalues.append(SUM)
    plt.plot(pvalues, yvalues)
plt.grid()
plt.show()
</code>```

I am using arbitrary values for all the parameters and a special relation between liquid bulk modulus, soil bulk modulus and porosity to determine the consolidation coefficient. These things are not relevant however. What is relevant it the initial loading p0 = 5000000 and the timevalues array assignment of 0. I am creating this array and looping over it because later I wish to use other time values. For now however, using the initial time of $t = 0 s$ I am getting the nonsensical result as shown below:

I am plotting the values for $t = 0$. Why am I getting this weird sine-like function rather than a constant function? I am also quite sure the analytical result is correct because I have tested it with sympy. Please help thanks.


问:
How to choose seed in flood/fill algorithm
I was able to draw to the outer airtight boundary of hole. Now I would like to implement interior flood/fill algorithm. To get an internal seed, only one method comes to my mind which is testing neighbors of red cells for ray casting. However, I would like to know is there a faster method.


问:
Extracting Real Part of Twiddle Factors from `fftw_plan`
When calculating a type-II DCT with FFTW, you create an fftw_plan, via

   fftw_plan plan = fftw_plan_r2r_1d(size, in_ptr, out_ptr, FFTW_REDFT10, FFTW_ESTIMATE);

As I understand it, the plan precomputes the twiddle factors so that the plan can be efficiently reused for multiple calls to fftw_execute_r2r(plan, in_ptr, out_ptr).

Since I'm trying to calculate
\begin{align*}
c_{k} = \frac{1}{N} \sum_{j=0}^{N-1} f\left(\cos\left(\frac{(j+1/2)\pi}{N}\right)\right) \cos\left(\frac{(j+1/2)k\pi}{N}\right)
\end{align*}
I figure that the fftw_plan already has computed the argument to $f$, so I don't want to recompute it.

How do you extract the twiddle factors from the fftw_plan?


问:
Solve ill-posed linear system without transposing matrices?
I am attempting to use an iterative solver to solve $p$ in

$$
Jp = -r
$$

where $J$ is an $m\times m$ matrix ($m$ is in the order of $10^5$ and never explicitly stored). $J$ is a dense matrix representing a Jacobian for a certain function, and, in general, is not symmetric.

Unfortunately, this is an ill-posed problem and possibly infinite versions of $p$ exist.

I cannot compute the least-squares version

$$
J^TJp=J^Tr
$$

because I only can access the product of $Jp$ and not $J^Tp$ (a limitation of an automatic differentiation approach I am taking).

Is there a way I can compute $p$ in $Jp = -r$ according to some criterion (e.g. least-squares solution) without transposing
$J$?

If I use an iterative solver (e.g. GMRES, BiCGStab, etc...) with an
initial guess of $p=0$, what solution will the solver converge to?
The minimum-norm?

ContentLicense: CC BY-SA 4.0


问:
Optimization Problem: Minimize the absolute value of a set of points
I have the following optimization problem:

Given is a set of $n$ points $x_i, \ldots, x_n$ with $x_i \in \mathbb{C}, \|x_i\| < \infty$. Note that $n$ is usually not large, i.e. $n < 1000$. We want to find a point $c \in \mathbb{C}$ such that

$\max_i \| x_i + c \|$ is minimized

Currently, this problem is solved by means of an iterative descent algorithm:

Initialize $c^{(0)} = 0$

For $j = 0, \ldots, j_\mathrm
{iter}$

Find $\max_i \| x_i + c^{(j)} \|$

Set $c^{(j+1)} = c^{(j)} - \alpha \frac{x_i + c^{(j)}}{\| x_i + c^{(j)} \|}$
,where $\alpha$ is some stepwidth.

This algorithm works fine, but: I am looking for another (faster) solution. I wonder if it is possible to use any non-iterative solver. So far, I couldn't come up with any. So do you have any ideas how to tackle this?
Thanks


问:
Dealing with simultaneous collisions in N-body sim
I have written an 2d N-body simulation in Python which allows collisions between the bodies.

A body is modeled as a circle whose area is proportional to its mass. Each time the sim advances by one time-step, I identify which bodies have overlapping radii and I replace them by one particle (respecting conservation of momentum and treating the collision as inelastic).

This worked well for N<100 but for higher N I encountered (not systematically, but depending on the initial conditions of the sim) a weird error: the number of particles N actually increased (and quickly exceeded the maximum integer value Python can handle).

I believe what is occuring is, after the new positions of the bodies have been calculated, it happens that (at least) 4 bodies mutually overlap. Each of the 6 pairs creates a new body and but we also delete the 4 bodies,so N increased by os we have 6-4=2. But the 6 new bodies are spawned in roughly the same positions as the original 4 bodies so now we have even more overlap (N can now increase by potentially (6 choose 2)-6=9 bodies). Hence why N can grow very fast.

How can I avoid this error?

The only thing I can think of is to reduce the radii of the bodies to reduce the change of 4 particles mutually overlapping. But this does not definitively fix the problem.


问:
3D Finite-Difference using Kronecker Products
As described in this Wikipedia article, a discrete laplacian matrix can be made for a 3D regular grid using Kronecker products. I'd like to use the same methodology for $(n-1)\times n$ matrices of the form

$$
\begin{bmatrix}
    -1 & 1 \\
    & \ddots & \ddots \\
    & & -1 & 1 \\
\end{bmatrix}
$$

Given 3 of these matrices with sizes appropriate for each dimension, I'd like use Kronecker products to build a matrix for 3D problems. For example, given 3D dimensions of $m\times n \times k$, I'd expect a finite-difference matrix of dimension $(m-1)(n-1)(k-1) \times mnk$.

How do I build this desired matrix using Kronecker products?


问:
3-parameter optimization of a convex function
I am modeling the structure of a crystal and need to find what the "lowest energy structure" is. There are three free parameters (let's say, $a$, $b$, and $c$) to change, and the search space, $E(a,b,c)$ is convex. What is the best way to set this up so I vary $a$, $b$, and $c$ to find the combination that yields the minimum of $E(a,b,c)$ without running extraneous calculations?

What I am currently doing is picking a relatively wide range of $a$, $b$, and $c$ values and computing $E$ for each combination. Then I fit a paraboloid to series of 3D plots of $E$ vs. $a$ vs. $b$ at fixed $c$ value and find where the global minimum is located. This was fine when my calculations were quick, but I need to be more efficient now with how to improve how I am searching the parameter space. One can think of something simple like making a range of $a$, $b$, and $c$, and keep expanding this range until $E$ increases when changing in all 3 dimensions from a given point, for instance.

Any suggestions of an adaptive search algorithm for a problem like this?


问:
Deterministic method to compute "Process noise covariance matrix, Q" for a Kalman filter when parameter variations of the model is known apriori
I am implementing a Kalman filter (for a linear ODE system for now).

My model represents a physical device that has 6 "parameters", i.e. those values of the device do not evolve over time (within a time-domain simulation environment for now).

I have 100 such "slightly differing devices", i.e. whose parameter values differ from the nominal, but the distribution of all 6 of the model parameters are known apriori (Gaussian, with known mean (nominal values) and standard deviations).

I am assuming that there is no other source of process noise (i.e. no unmodeled states or other dynamics are present)

I am implementing a Kalman Filter with a single tuning parameter across all these 100 devices. I think that these deviations in model parameters (from the nominal) can be treated as process noise.
Is there a deterministic way, i.e. a recipe, to compute the "Process noise covariance matrix", typically denoted as $Q$ in Kalman Filter literature?

I am thinking that one can somehow exploit the affine properties of Guassian distributions to arrive at the error covariance matrix, although I am a bit rusty with the math.

Partially related, but different question asked by someone else:
https://dsp.stackexchange.com/questions/21796/question-about-q-matrix-noise-process-covariance-in-kalman-filter


问:
How to optimize for decay constant in exponential-like function?
I've got a data set of points $M_O .. M_N$ for time points $t_0 .. t_N$, where $N$ is approximately 10-20, and the spacing of time is not uniform (i.e., $t_{i+1}-t_i$ is not constant for all i). It is believed that this data should be able to be modeled by a differential equation of the form

$\frac{dM}{dt} = -k (M - M_{e})$

My goal is to solve for the constant $k$. If $M_e$ were a constant, then the solution of the function $M$ would of course be something like

$M = M_{e} + Ce^{-kt}$

and I could then fit this function to the data and determine $k$ with any regular minimization routine (say, using the RMSE).

However, $M_{e}$ is not a constant, but instead varies with $t$. I already have the values of $M_e$ for all $N$ points in my data set.

My question is, what is the best way to numerically optimize for the $k$ which fits the model best? Naively, I could approximate $\frac{dM}{dt}\approx \frac{\Delta M}{\Delta t} $, compute the differences for each segment, and choose something like the mean value of $k$ from each $\Delta t$. I'm curious as to whether this is the best way to do this, or if there's another more appropriate method. One possible weakness of what I'm describing is that averaging $k$ for all segments ignores the fact that some $\Delta t$ intervals are larger than others, and perhaps smaller $\Delta t$ should be systematically favored over larger time intervals. Finally, is there any good method to get a somewhat closed form solution for $M(t)$, or does this hinge on whether there's a simple closed-form solution to $M_e(t)$ (which is doubtful in my case because of its specific physical origin).


问:
Selecting n points to have a given mean and covariance
I have a set $N$ points in $d$-dimensional space, $\mathbf{x}_i \in \mathbb{R}^d ~~ \forall i=\{1,2,...,N\}$. How to choose $n$ points from the set such that it maximises
$$
f = |\mathbf{S}|^{-n/2} \exp\left[-\frac{1}{2} \mathrm{tr}\left(\mathbf{S}\boldsymbol{\Sigma}^{-1}\right) - \frac{n}{2} (\boldsymbol{\mu} - \mathbf{\bar{x}})^T \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu} - \mathbf{\bar{x}})\right],
$$
for given $\boldsymbol{\Sigma}$ and $\boldsymbol{\mu}$, where $\mathbf{S} = \sum_j (\mathbf{x}_j - \mathbf{\bar{x}})  (\mathbf{x}_j - \mathbf{\bar{x}})^T$ and $\mathbf{\bar{x}} = \sum_j\mathbf{x}_j/n$?

The solution I am thinking of is to try combinations of the points, but that will take exponential time to try all combinations. I wonder if there is a faster way, utilising the knowledge of function $f$? I don't need an exact set of point that has the global maximum of $f$. Random algorithms that achieve reasonable result in faster time are much appreciated.


问:
Why do the policy iteration and value function iteration not match in this exercise
I'm working through the book Approximate Dynamic Porgramming by Warren B. Powell. Currently I'm solving exercise 3.3 and 3.5. Let me state the exercise (I'm summarizing the text from the boook)

We have the following Markov Decision Problem. We have a two-stage Markov chain with one-step transition matrix
  $$
P=
\begin{pmatrix}
0.7 & 0.3  \\
0.05 & 0.95 \\
\end{pmatrix}
 $$
  The contribution from each transition from state $i\in\{1,2\}$ to state $j\in\{1,2\} $is given by
  $$
C = \begin{pmatrix}
10 & 30  \\
20 & 10 \\
\end{pmatrix}
 $$

that means going from state $1$ to state $2$ has a contribution (reward) of $30$.

Since the exercise are asking to solve the problem using value iteration (infinite horizon) and policy iteration let me state the two algorithms from the book (again I'm quoting)

Value Iteration

Set the initial value $v^0 = 0 $ for all state $s\in S$, fix a tolerance parameter
  $\epsilon$ and set $n=1$.

For each state $s$  compute
  $$ v^n(s) = \max_a (C(s,A) + \gamma\sum_{s'}P(s'|s,a) v^{n-1}(s'))$$

If $\|v^n-v^{n-1}\| <\epsilon(1-\gamma)/2\gamma $ let $ \pi^\epsilon$ be the resulting policy and $v^n$ the resulting value. Else set $n=n+1$ and go to step 2

and the policy iteration (again quoting the book)

Policy Iteration

Select a policy $\pi^0$ and set $n=1$.
  $\epsilon$ and set $n=1$.

Given the policy $\pi^{n-1}$ compute the transition matrix $P^n$ and the contribution vector $c^n$ where the element for state $s$ is given by $c^{n-1} = C(s, \pi^{n-1})$

let $v^{n-1}$ be the solution to ($I$ the identity matrix) $$(I - \gamma P^n) v = c^{n-1}$$

find a policy $\pi^n$ defined by $$ a^n(s) = \arg\max_a(C(a)+\gamma P^nv^n)$$ this requires to compute the action $a$ for each state $s$.

if $a^n(s)=a^{n-1}(s)$ for all $s$ then stop. Otherwise set $n=n+1$ and go to step $2$.

Now let me state the two exercises. In exercise $3.3$ we need to achieve the following two tasks:

Exercise 3.3
  Plot the value of being in state 1 as a function of the number of iterations if your initial estimate of value being in each state is $0$. Show the graph for $100$ iteration. Repeat the exercise with initial estimates of $100$ each and $100$ for state $1$ and $0$ for state. The author says (note that you are not choosing a decision so there is no maximization step)

I've used R to implement the value iteration:

```<code>P <- matrix(c(0.7, 0.3, 0.05, 0.95), nrow=2, byrow=T)
C <- matrix(c(10, 30, 20, 5),nrow=2, byrow = T)
discount <- 0.8
v <- matrix(rep(0,200), nrow=100)
for(i in 2:100){
  v[i,] <- (C[1,] + discount*sum(P[1,]*v[i-1,1]))
}
</code>```

where I've used $P[1, ]$ since we are only interested in the value of being in state $1$. Here are the graphs for the different initial states

both $0$ initial value

both $100$ initial value

and $100$ and $0$ initial value

FIRST QUESTION is the solution correct? I'm a bit unsure since the last initial value pair did not change the plot compared to the two above ones.

the second exercise is the following

Exercise 3.5
  Apply policy iteration to the above problem. Plot the average value function after each iteration alongside the average value function found using value iteration after each iteration (initial value $0$).

With this I'm struggling. It's not really clear to me how to use policy iteration if there is no policy to choose from as the author outlined above.

From the discussion with Kirill in the comments you can see that I've first tried to introduce a policy to solve this exercise which seems to be wrong. However, for completeness I leave the results in the old outdated part below.

OLD OUTDATED PART

Clearly there are policies of the form: being in state $i$ go either to state $j$ or stay at state $i$. We can model this via tuples $(i,j)$ meaning if I'm in state $1$ I go to state $i$ and if I'm in state $2$ I go to state $j$. I've also implemented the algorithm in R. Note that several steps of the above policy iteration are not needed to be recalculated since they are constant. What confuses me is that the author say above there are no decision but now we should solve it via policy iteration. Moreover, my resulting value is different which means something is wrong. Can someone help me to understand what's going wrong.

```<code>P <- matrix(c(0.7, 0.3, 0.05, 0.95), nrow=2, byrow=T)
C <- matrix(c(10, 30, 20, 5),nrow=2, byrow = T)
discount <- 0.8
v <- c(0,0)
current_policy = c(1,2)
id <- matrix(c(1,0,0,1), nrow=2, byrow=T)
inverse_factor = solve(id-discount*P)
new_policy = c(1,1)
while(any(new_policy != current_policy)){
  current_policy = new_policy
  c_pi = c(C[1,current_policy[1]],C[2,current_policy[2]])
  current_value = inverse_factor%*%c_pi
  a_1 = which.max(c(C[1,1]+discount*(P%*%current_value)[1],C[1,2]+discount*(P%*%current_value)[2]))
  a_2 = which.max(c(C[2,1]+discount*(P%*%current_value)[1],C[2,2]+discount*(P%*%current_value)[2]))
  new_policy = c(a_1,a_2)
  v <- rbind(v,t(current_value))
}
</code>```


问:
Von Neumann analysis for coupled PDE's
I am interested in understanding how to perform stability analysis for coupled (to keep things do-able, lets say linear) PDE's

In the case of a single PDE, i understand the logic behind the VN analysis: mode expansion and look at the amplification factor of any given mode.

However for coupled PDE's it seems like you won't be able to cancel things as nicely.

For example, basing off the notes here
https://en.wikipedia.org/wiki/Von_Neumann_stability_analysis

specifically starting at Eqn (2): If this represents the error for one function it will depend on the other via the coupling. It seems that the $e^{at}$ term will not cancel nicely to get Eqn (6). Presumably i can still restrict the analysis to one mode at a time assuming both PDEs are linear?

That is, rather than the $e^{at}$ terms canceling as they do in that simple example, in the coupled case i would be left terms containing the relative size i.e. $e^{(a-b)t}$ where $a$ and $b$ are the exponents of the mode expansion amplitudes of the errors for the two PDE functions respectively (whoa, a mouthful...)

It seems like the solution is to cast as a linear matrix problem but I get a bit stuck here too. Certainly one can write the time discretization in terms of an amplification factor/matrix $A$ acting $n$ times on some initial vector to give the new state, but how does one then look at the spatial discretization (contained inside the amp matrix $A$ in this context? Or do we just look for eigenvalues of that amplification matrix and then do finite-difference afterwards?

Any and all help is appreciated.

Title: Von Neumann analysis for coupled PDE's


问:
Assembled human ESTs for annotation
I am annotating protein-coding genes in a region of the human genome and am looking for a set of assembled ESTs I can use as evidence for constructing gene models. A quick search of NCBI's SRA shows 4,318 human RNA-seq data sets, so I understand there is no shortage of data. I am hoping, however, to save the time that would be required to search through this massive amount of data, select an appropriate subset, and assemble the ESTs myself. I'm sure this has been done many times before, and for this particular task it doesn't make sense to repeat this process.

Is there any sort of general (i.e. not tissue-specific) reference EST/transcriptome assembly available for Homo sapiens?

Title: Assembled human ESTs for annotation


问:
Multi-matrix orthogonal basis problem
Suppose we are given a set of symmetric, positive definite matrices $A_1,A_2,\ldots,A_k\in\mathbb{R}^{n\times n}$.  Is there any numerical method or reduction to a known problem (e.g. eigenvalue computation) that can solve the following problem?

$$
\begin{array}{rl}
\min_{v_1,\ldots,v_k\in\mathbb{R}^n} & \sum_{i=1}^kv_i^\top A_iv_i\\
\textrm{subject to }&v_i^\top v_j=\delta_{[i=j]}.
\end{array}
$$

Note: The problem has some resemblance to the "joint diagonalization" problem, e.g. in this paper.  The difference seems to be that in joint diagonalization you use the same basis to diagonalize multiple matrices.

Note 2: If the $A_i$'s do not have the same eigenvectors, then the $v_i$'s will not be exactly eigenvectors of the $A_i$'s.

Note 3:  My intuition is that this will give a set of mutually orthonormal vectors $v_1,\ldots,v_k\in\mathbb{R}^n$, with the property that vector $v_i$ is similar to an eigenvector of $A_i$. But this is not an eigenvalue problem.

Title: Multi-matrix orthogonal basis problem


问:
Numerical scheme to solve Maxwell equations with fixed potential boundaries?
We have a 2D electromagnetic field (in the sense that: $E=(E_x,E_y,0)$, $B=(0,0,B_z)$, and all derivatives with respect to $z$ are $0$), and we are considering a system made up of two walls at $x=-b$ and $x=b$ held at fixed potential $\phi(y)=e^{-y}$, with a third conducting wall at $y=0$ held at $\phi(x)=1$. The current in our system is $\vec{j} \propto \nabla \cdot E\left(B \times E+\sigma E \right)$. The initial conditions are that $E=-\nabla \phi$, where $\phi$ is the solution to the Laplace equation $\nabla^2 \phi = 0$, and $B=0$.

Yee's scheme is resulting in diverging $E_y$, so what numerical scheme would you suggest be used to solve Maxwell's equations in this case?


问:
Reconstruction of cells information when given mesh in polyMesh format
The polyMesh format used in OpenFOAM is an intelligent format for which 'cells' file, i.e. the file giving list of nodes that make a cell, is not needed. The way to reconstruct geometry data such as cell centroids and cell volume is by using divergence theorem an therefore only face data is needed.

If the OpenFOAM is used, the post-processing tool Paraview knows how to build a geometry for visualization given the OpenFOAM case.

However, if I want, with my own code that can read and process polyMesh format, to write a Paraview ".vtu" format (a format for unstructured mesh) I need to give cell data.

This is where I found a problem, and still can't find proper resolution.

A side note is that I use CSR format for sparse matrices, which I create using 'owner' and 'neighbour' polyMesh data. From CSR format I can easily find cell neighbour indices in the form of ja(ia(icell)),..,ja(ia(icell+1)-1), except the diagonal element, which is the particular cell. The 'ja', 'ia' are typical CSR sparse matrix format arrays.

To make it precise:
I'm looking for a list of nodes for each cell (or element), given the list of nodes for each cell face (element face), two lists of cells (elements) that share the face, a list for each side of the face, and of course the list of nodal coordinates. The list of nodes for each cell should be in correct order, consistent with faces list.


问:
General-purpose Numerical Calculus of Variations Problem
Calculus of variations problems are generally cast in in the following simple form: find $u(t)$ that satisfies some boundary conditions and minimises

$$
F[u] = \int_{t=0}^{t=t_f} f(u(t),u'(t),t) dt.
$$

I have run into a far more general problem. I have a functional $F[u]$ that is substantially more complicated. Given $u$ (or in this case the set of functions $u^i(\mathbf{x})$ where $\mathbf{x}\in\mathbb{R}^n$) I solve a pde which depends on these functions, and the solution is analysed to give the value of the functional $F$. Given a set of functions $u^i(\mathbf{x})$, I can numerically evaluate $F[u^i]$ very efficiently. However, I want to explore the function space to find the set of $u^i$ that maximises $F$. To complicate matters, I have further constraints on the $u^i$ I can choose: for instance, I require a particular path integral of $u^i$ to equal a certain value.

All the methods I have found in the calculus of variations seem to assume, in one way or another, that the functional $F$ is of the simple form outlined above. Do more general numerical methods exist? I'm looking for a method or algorithm that requires nothing more than the ability to numerically evaluate the functional for any given function, and from that will search the function space and converge on an optimal solution. The ability to place (ideally arbitrary) constrains on $u^i$ would be an added bonus.

Perhaps I ask for too much. I'd appreciate any help or advice.

Thanks


问:
Linear regression with inequality constraint in Java
I haven't been doing math in years and I'm facing the following problem. I'm trying to implement in Java a linear regression under a set of inequality constraints. Sorry in advance for all the approximations while trying to formulate this problem...

More precisely, I have a linear function $f_\Theta:R^n->R$ defined as such: $f_\Theta(X) = \theta_0 + \theta_1x_1+ \theta_0x_2+... \theta_nx_n$, with  $\Theta$ a (n+1) vector.

My aim is to $\min_\Theta \sum_{i}||f_\Theta(X_i)-Y_i||^2$ under the constraint that $f_\Theta(X_i)-Y_i < 0$ for a maximum of $i$ (I wish I could avoid to take into account "obvious" outliers in my observations)

So, bottom line, I'm looking for the linear regression that will act as the "best fit" providing a lower bound for most/all the observations.

How can I do that ? I've been looking for a Java Library that would implement such a thing. I known the Apache Common Math package and I also found this Flanagan library, but it does not seem to address this.

Any help greatly appreciated. Thanks

OwnerDisplayName: BlackLabrador


问:
Appropriately handling boundary conditions in a PDE eigenvalue problem
Suppose I have a nonlinear second order Cauchy PDE $\dfrac{\partial p(x,t)}{\partial t}=N(p(x,t))$, where $N:L^2(\mathbb{R})\rightarrow L^2(\mathbb{R})$, and a known fixed point $u(x)$. Mathematically, the only boundary condition known is that any solution $p$ decays to 0 as ${x\rightarrow\pm \infty}$.

From physical reasoning however, it is known that $u$ decays rapidly away from $x=0$, and hence, I am able to solve for $u$ numerically on a finite domain by prescribing zero Dirichlet boundary conditions on a finite domain around 0 (i.e, on $[-L,L]$, where $L$ is chosen judiciously).

Now, I want to linearize $N$ around $u$, and find eigenvalues of the resulting linear operator $L_{u}$. I have a dilemma as to what boundary conditions should I impose on the eigenfunctions to be able to solve this problem numerically on a finite domain.

Are there any standard ways to do this ? Is there a way to justify same boundary conditions as $u$.


问:
Minimum box size in order for an object in the center to not 'see' duplications of itself under periodic boundary conditions?
I am working on developing a Dissipative Particle Dynamics (DPD) model of a colloid in a bulk fluid. I essentially want to ensure that the perturbations of my meshed colloid are completely uncoupled from the duplicated copies of itself that are present as a result of using periodic boundary condition. This includes any hierarchical interactions. For example, if fluid particles are disturbed from interacting with  a duplicated version of the colloid, the propagation of these interactions should dissipate out enough such that by the time this propagation reaches the colloid of interest, that it essentially has a negligible effect on the behavior of the colloid. Unfortunately I myself am not certain about what 'negligible' means, but perhaps someone with more experience in the field could offer some advice on that point in particular.

As a starting point, I would assume the following:

$L-D > r_c$,

where $L$ is the box length, $D$ is the diameter of the colloid and $r_c$ is the largest cutoff for the bead-bead interactions forces. This ensures that the colloid will not interact with itself as a direct result of colloid-colloid interactions, but does not ensure the prevention of the hierarchical interactions I described earlier.


问:
Definition of CFL number in Arbitrary Lagrangian-Eulerian framework
In an Eulerian frame of reference, the CFL number is defined as $$\sigma=\frac{u \Delta t}{\Delta x}$$
with $u$ the magnitude of the fluid velocity. A restriction such as $\sigma<1$ for time integration then roughly means that the initial information travels at most over one cell-width during one time-step.

In an Arbitrary Lagrangian-Eulerian (ALE) framework the grid also moves and the convective velocity becomes $u-u_\mathrm{grid}$. Should the definition of the CFL number now be based on $u-u_\mathrm{grid}$ rather than $u$?


问:
Round-off errors in ratio of Gaussian functions
I would like to understand how division of round-off errors can lead to large errors???

I am using a method to measure the frequency dependent dissipation in simulations performed by Discontinuous Galerkin(Finite element based) Time Domain method on wave propagation problems. To this end, I excite my domain with a Gaussian pulse and measure the propagated signal $y_1$ and $y_2$ at 2 monitoring points $z_1$ and $z_2$ of the domain. The discrete Fourier transform(DFT) $F_1(f)$ and $F_2(f)$ of both signals are computed and the dissipation $\alpha_{diss}$ is evaluated through :

$\alpha_{diss}(f) = \frac{log(|F_2(f)|/|F_1(f)|)}{z_2-z_1}$

If there is no dissipation, the amplitude spectrums $|F_1|$ and $|F_2|$ should theoretically be the exact same Gaussian $F_{th}$. However, due to finite length signal, the DFT makes round-off errors(different for $y_1$ and $y_2$). As a result, there is an increasing overall error in the measure of $\alpha_{diss}(f)$ which prevents the method to be used far away from the mid-band frequency $f_0$ of the Gaussian spectrum.

The figure below shows this error increase of $|\alpha_{diss}(f)-0|$(in log scale) in red. The blue and cyan curves show the round-off errors $F_2-F_{th}$ and $F_1-F_{th}$ respectively. The mid-band frequency is $f_0=2$ and Gaussian line-width $\Delta f_0=0.5$. We can clearly see that beyond $3\Delta f_0$ the evaluation of $\alpha_{diss}(f)$ is not reliable anymore.

I really wonder why is this happening?


问:
Plotting yield surface for geotechnical constitutive model
I have found the following example which is produced as documentation for RS3 geotechnical FE software. This is a validation document that goes through an analytical verification of triaxial testing for modified cam clay material model. It relies heavily upon the p-q stress representation of the model to make its argument.

I find it to be quite pedagogical in the way it presents its content. It has enabled me to somewhat understand the following diagram:

For this example triaxial test the isotropic stresses are initialized at point $(p,q) =(p_i, q_i)=(200,0)$. The deviatoric stress $q$ is then applied until the sample 'fails'. This is the intersection between the green line and the red line (CSL).

This means that for a given axisymmetric FE analysis with following boundary conditions:

It is only possible to generate ONE p-q green dashed line.

I would like to know if there is a way to plot the ENTIRE surface. This would perhaps imply that I would have to force the sample to take several different stress paths and have them each intersect with the ellipse? Computational efficiency is not really the issue here, it is more the method. How do I force the sample to take different stress paths such that I can plot the entire yield surface? Is this even possible? Am I even asking the right question?


问:
Code for solving the heat equation on the semi-infinite rod
Cross posted in mathematica.SE.

Question : I want to test the solution which is given below is right by Matlab/Maple/Mathematica.

Please look the post in mathstackexhange

or

Please look below.

Q: Solve the following heat equation on the semi-infinite rod

$$u_t=ku_{xx}$$ where $x,t>0$ and
    $u_x(0,t) =0$ and $u(x,0)=\begin{cases}
      1, & 0 < x <2 \\
      0, &  2\leq x
   \end{cases}
$

Answer

We found the following answer
$$
u(x,t) = \frac{2}{\pi}\int_{0}^{\infty}e^{-s^2 t}\frac{1-\cos(2s)}{s}\cos(sx)ds.
$$

Code

But I am not sure the solution is right. I am not capable of to write a code in order to test it. Could you help me?


问:
Successive iteration method for solving eigenvalue ploblem
I have a question concerning the branch of successive iteration methods (Newton, Runge-Kutta). I definitely know (or can read in Wikipedia) the implementation of these methods. But I was wondering about the precision in special case.

I have a function
$$
D(\omega)=det(L_{m,n}(\omega)), \omega = \omega'-i\omega''
$$
and need to solve the equation
$$
D(\omega)=0
$$
For that task I want to use simple fixed-point iteration method
$$
\omega_{0}, \omega_{1} \mathrm{- given}\\
\omega_{n+1}=\alpha \omega_{n}+(1-\alpha)[\omega_{n}-\frac{D(\omega_{n})}{dD/d\omega|_{\omega_{n}}}]\\
dD/d\omega|_{\omega_{n}}=\frac{D(\omega_{n})-D(\omega_{n-1})}{\omega_{n}-\omega_{n-1}}
$$
So, could someone tell me if this method is well converging or not.

I am asking because previously I've used this method only for real-valued variables and now I have complex-valued. This may lead to component entanglement while calculating the derivative.


问:
CFD Mass exchange across a reactive front
I have a CFD solver at hand which can deal with 2 phases and so far is able to exchange momentum and energy across the interface, where both phases meet.
How this is done is well documented in the following paper "A conservative interface method for compressible flows" by Hu, Khoo and Adams, which publicly available.

In the paper, equation (10) is giving the general idea of how to update the conservative variables in so called "cut-cells", cell's which are cut by the interface (tracked by a level-set method).

In this update routine $ \hat{X} $ holds the average energy and momentum exchange. Now I want to extend the code such that it can deal with reactive fronts (like combustion or cavitation).

The arising Riemann problem is now solved by a reactive Riemann solver, which provides next to pressure and velocity across the contact discontinuity, also the mass flow through the reactive front.

So far all good, but now I am having troubles to adjust $\hat{X}$ such that it transports the mass (density) across the front.

I would be very grateful for any thoughts or hints.


问:
Hello word in FEniCS?
I am trying to start using FEniCS, but have a problem with the simple hello world examples given in the books. Could you please give me the simplest hello world example for solving the Fick's second law.

$$\frac{\partial c}{\partial t} = D \frac{\partial^2 c}{\partial x^2}$$

How can we apply the boundary condition to the simplest example of FEniCS for solving this PDE (e.g., drawing the concentration profiles at given times)?

Title: Hello word in FEniCS?


问:
Eigenvalue with largest imaginary part
Iterative eigensolvers such as ARPACK, give the option to find a subset of the eigenvalues which have the largest imaginary part. My question is how do these algorithms work.

As I understand it, generally such algorithms use Krylov spaces and shift-invert methods which are very good at finding eigenvalues either at the extremes of the spectrum or around a given value. This is essentially because repeated applications of a matrix projects onto the eigenvector with largest magnitude eigenvalue. However, since the imaginary part of an eigenvalue need not have any correlation with its real part, I can't see how it's able to `ignore' the real part.

If anyone has any idea how these algorithms work, or alternatively know a way to search for eigenvalues which is ignorant of the real part of the eigenvalue I'd be very grateful.

Note that I previously asked this question on the maths stackexchange but didn't receive any answers:

https://math.stackexchange.com/questions/2437587/eigenvalue-with-largest-imaginary-part

Title: Eigenvalue with largest imaginary part


问:
Algorithm design to filter on 5,000 stocks each of which has 4 months worth of data points
I want to filter on 5000 stocks, each of which has 4 month or more worth of data (>= 500 data points each). my filtering criteria will be based on 8 calculated values from the data points. for example, moving average, standard deviation and values alike.

my tools of choice are python pandas, numpy and CUDA C++; please take this as the default setting which can be expanded on, especially on the C++ side, I believe I can expand on Linear Algebra packs but I am not very familiar with them.

my current design is to:

read the data into pandas dataframe.

pass one stock at a time to C++ to compute all filtering criteria values.    here I plan to use CUDA to speed up the calculations for
  each criteria value.    e.g. running moving average calculation
  parallel on CUDA cores.

pass a list containing filter criteria values of each stock back to python.

sorting the stocks in terms of the filter criteria values; each criterion will weigh differently when I do the sorting.

I am still exploring, I know there will be better design/strategy with different tool set, but with these tools as given, is there a better strategy? Can I use Linear Algebra to speed this even further (potentially less codes)?


问:
Why does PETSc take unexpectedly long to set up its KSP solver with a custom preconditioner?
I am attempting to solve a large system, $\bf{Ax} = \bf{b}$ with the help of PETSc. Due to the size of the problem, I'm using a matrix-free approach, where $\bf{A}$ is just a shell. I'm also providing my own preconditioner (which is not a shell), and I'm using ilu(2) factorization on the preconditioner.

The issue: the setup stage of the solver (see relevant block of code below) is taking an extremely long time. I suspect that this would mostly be the ilu of the preconditioner that's taking a while. I know that it is expected that the ilu will take time, but I'm concerned for the following reason: when I try to solve the same problem using a direct solver (using MKL Lapacke to find the LU decomposition and then inversion of $\bf{A}$, outside of PETSc), the LU is 10x faster. I expect that PETSc's ilu should take a comparable time to the full LU factorization, but it is 10x slower, which seems very strange. (By the way, you may ask, why even use an iterative solver if I can do it with LU, but this example is not as large as the ones I actually want to run, at which point I won't be able to use a direct solver).

Here is the code snippet relevant to the issue:

```<code>MatCreateShell(comm, Nu, Nu, Nu, Nu, ctx, &A_shell);
MatShellSetOperation(A_shell, MATOP_MULT, (void(*)(void))usermult);

KSPCreate(comm, &solver);
KSPSetOperators(solver, A_shell, PreconditionerMatrix);
KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
KSPSetNormType(solver, KSP_NORM_UNPRECONDITIONED);

KSPSetFromOptions(solver);
KSPSetUp(solver);
</code>```

Things I know / have tried:

The condition number of the matrix can be as large as $10^7$, but I don't think that my issue has anything to do with that, because again, the time sink is in the setup. If this was the problem, it would also manifest itself when I do the full LU decomposition, but it doesn't.

I know that I have to provide a large enough fill factor guess for the post-ilu matrix, and I've set the option -pc_factor_fill to 3. After running the code with -info, I confirmed that this is sufficient to prevent any reallocation of memory. Interesting side note: when I do run it with -info, it reports the required fill factor fairly quickly. Does this mean that it actually does the ilu expectedly fast, but then gets stuck somewhere else? Am I barking up the wrong tree? Here's what it reports:

```<code>  [0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
  [0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
  [0] PCSetUp(): Setting up PC for first time
  [0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
  [0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
  [0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
  [0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
  [0] MatILUFactorSymbolic_SeqAIJ(): Reallocs 0 Fill ratio:given 3. needed 1.7385
  [0] MatILUFactorSymbolic_SeqAIJ(): Run with -[sub_]pc_factor_fill 1.7385 or use
  [0] MatILUFactorSymbolic_SeqAIJ(): PCFactorSetFill([sub]pc,1.7385);
  [0] MatILUFactorSymbolic_SeqAIJ(): for best performance.
  [0] MatSeqAIJCheckInode_FactorLU(): Found 2030 nodes of 6096. Limit used: 5. Using Inode routines
</code>```

Then it gets stuck for a really long time...so maybe I set the fill factor too large? I tried the same thing again with a fill factor of 2 instead, but it made no difference.

I have already made sure I'm not using the debugging installation of PETSc; when timing the code, I'm definitely using --with-debugging=0 in the PETSc configuration step.

I am not using any parallelization.

 Here is the output generated with -log_view:

```<code>Using Petsc Release Version 3.7.6, Apr, 24, 2017

                         Max       Max/Min        Avg      Total
Time (sec):           4.057e+02      1.00000   4.057e+02
Objects:              7.050e+02      1.00000   7.050e+02
Flops:                2.161e+11      1.00000   2.161e+11  2.161e+11
Flops/sec:            5.327e+08      1.00000   5.327e+08  5.327e+08
MPI Messages:         0.000e+00      0.00000   0.000e+00  0.000e+00
MPI Message Lengths:  0.000e+00      0.00000   0.000e+00  0.000e+00
MPI Reductions:       0.000e+00      0.00000

Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
                            e.g., VecAXPY() for real vectors of length N --> 2N flops
                            and VecAXPY() for complex vectors of length N --> 8N flops

Summary of Stages:   ----- Time ------  ----- Flops -----  --- Messages ---  -- Message Lengths --  -- Reductions --
                        Avg     %Total     Avg     %Total   counts   %Total     Avg         %Total   counts   %Total
 0:      Main Stage: 4.0567e+02 100.0%  2.1612e+11 100.0%  0.000e+00   0.0%  0.000e+00        0.0%  0.000e+00   0.0%

------------------------------------------------------------------------------------------------------------------------
See the 'Profiling' chapter of the users' manual for details on interpreting output.
Phase summary info:
   Count: number of times phase was executed
   Time and Flops: Max - maximum over all processors
                   Ratio - ratio of maximum to minimum over all processors
   Mess: number of messages sent
   Avg. len: average message length (bytes)
   Reduct: number of global reductions
   Global: entire computation
   Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().
      %T - percent time in this phase         %F - percent flops in this phase
      %M - percent messages in this phase     %L - percent message lengths in this phase
      %R - percent reductions in this phase
   Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)
------------------------------------------------------------------------------------------------------------------------
Event                Count      Time (sec)     Flops                             --- Global ---  --- Stage ---   Total
                   Max Ratio  Max     Ratio   Max  Ratio  Mess   Avg len Reduct  %T %F %M %L %R  %T %F %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------

--- Event Stage 0: Main Stage

MatMult              624 1.0 3.2336e+01 1.0 4.23e+10 1.0 0.0e+00 0.0e+00 0.0e+00  8 20  0  0  0   8 20  0  0  0  1307
MatMultAdd          2480 1.0 5.4122e-02 1.0 5.00e+07 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0   923
MatMultTranspose    3100 1.0 2.8153e+00 1.0 9.05e+09 1.0 0.0e+00 0.0e+00 0.0e+00  1  4  0  0  0   1  4  0  0  0  3215
MatSolve             608 1.0 4.6020e+01 1.0 7.41e+10 1.0 0.0e+00 0.0e+00 0.0e+00 11 34  0  0  0  11 34  0  0  0  1611
MatLUFactorNum         1 1.0 4.4004e+01 1.0 9.31e+10 1.0 0.0e+00 0.0e+00 0.0e+00 11 43  0  0  0  11 43  0  0  0  2115
MatILUFactorSym        1 1.0 3.4659e+01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  9  0  0  0  0   9  0  0  0  0     0
MatAssemblyBegin      27 1.0 1.5497e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatAssemblyEnd        27 1.0 6.8570e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatGetRow        10650098 1.0 6.6575e+00 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  2  0  0  0  0   2  0  0  0  0     0
MatGetRowIJ            1 1.0 9.5367e-07 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatGetOrdering         1 1.0 1.3018e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatZeroEntries        39 1.0 1.5883e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatMatMult             2 1.0 3.5977e+00 1.0 5.97e+09 1.0 0.0e+00 0.0e+00 0.0e+00  1  3  0  0  0   1  3  0  0  0  1660
MatMatMultSym          2 1.0 7.2353e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatMatMultNum          2 1.0 2.8741e+00 1.0 5.97e+09 1.0 0.0e+00 0.0e+00 0.0e+00  1  3  0  0  0   1  3  0  0  0  2078
VecMDot              302 1.0 4.2068e-02 1.0 2.06e+08 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  4907
VecNorm              621 1.0 1.4746e-02 1.0 3.03e+07 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  2054
VecScale             314 1.0 2.0843e-03 1.0 7.66e+06 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  3674
VecCopy              636 1.0 1.0492e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecSet              4096 1.0 4.1316e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecAXPY             5278 1.0 9.9347e-02 1.0 1.96e+08 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  1977
VecAYPX              306 1.0 4.3933e-03 1.0 7.46e+06 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  1698
VecMAXPY             608 1.0 5.9476e-02 1.0 4.16e+08 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  6993
VecAssemblyBegin    2493 1.0 1.5733e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecAssemblyEnd      2493 1.0 1.5340e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecNormalize         314 1.0 1.0375e-02 1.0 2.30e+07 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  2214
KSPGMRESOrthog       302 1.0 7.3104e-02 1.0 4.13e+08 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  5647
KSPSetUp               1 1.0 1.2398e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
KSPSolve               4 1.0 7.8504e+01 1.0 1.17e+11 1.0 0.0e+00 0.0e+00 0.0e+00 19 54  0  0  0  19 54  0  0  0  1491
PCSetUp                1 1.0 7.8663e+01 1.0 9.31e+10 1.0 0.0e+00 0.0e+00 0.0e+00 19 43  0  0  0  19 43  0  0  0  1183
PCApply              608 1.0 4.6022e+01 1.0 7.41e+10 1.0 0.0e+00 0.0e+00 0.0e+00 11 34  0  0  0  11 34  0  0  0  1611
------------------------------------------------------------------------------------------------------------------------

Memory usage is given in bytes:

Object Type          Creations   Destructions     Memory  Descendants' Mem.
Reports information only for process 0.

--- Event Stage 0: Main Stage

              Matrix    33             30   1306648980     0.
              Vector   663            663     65247472     0.
       Krylov Solver     1              1        35264     0.
      Preconditioner     1              1         1008     0.
              Viewer     2              0            0     0.
           Index Set     5              5        87624     0.
========================================================================================================================
Average time to get PetscTime(): 5.96046e-07
#PETSc Option Table entries:
-ksp_atol 1e-8
-ksp_converged_reason
-ksp_monitor
-ksp_monitor_true_residual
-ksp_rtol 1e-8
-log_view
-pc_factor_fill 3
-pc_factor_levels 2
-pc_type ilu
#End of PETSc Option Table entries
Compiled with FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 16 sizeof(PetscInt) 4
Configure options: PETSC_ARCH=arch-linux2-cxx-nodebug --with-scalar-type=complex --with-fortran-kernels=1 --with-clanguage=c++ --with-debugging=0 --with-cxx=g++ CXXOPTFLAGS=-O3 COPTFLAGS=O3 FOPTFLAGS=-O3 --download-openmpi --with-blaslapack-dir=/opt/intel/mkl
</code>```

Questions

I know that there may not be an easy / obvious solution to this, but at least I'd like insight on why this happens, and what PETSc is doing internally that takes this long

If there is no clear solution, what are some steps I can take to try and mitigate this or investigate this further?

Is this expected / normal, and should I stop worrying about it and just suck it up?

Sorry about the length of this question. Thanks for your time!


问:
Piecewise Constant Enrichment of the continuous galerkin method
I am interested to study crack propagation in a hyperelastic material in a variational setting. The crack surface exhibits a jump discontinuity. The function space for displacement field should consist of discontinuous functions.

The discontinuous galerkin finite element (DGFEM) approximation space consist of piecewise discontinuous functions which seems to a valid choice to solve crack problem. However, the required degree of freedom is extremely large for a high resolution discretization.

The recently proposed, enriched galerkin finite element method (EGFEM) enriched the continuous galerkin approximation space with a piecewise constant function. This allows discontinuities along the element face. For example, a one dimensional function $f(x)$ can be approximated in a unit interval as,

$ f(x) = f_1 \phi_1(x) + f_2 \phi_2(x) + f_3$

where $\phi_1(x)$ and $\phi_2(x)$ are basis functions. EGFEM seems to be a better choice to solve crack problem as required degree of freedom is quite less than for DGFEM.

EGFEM does introduce the discontinuity along the edges/faces, but in the deformed configuration the disjoint edges/faces would always be parallel. This is one of the limitation of EGFEM.

Before going to implement these methods, I would like to understand advantage and disadvantage of DGFEM and EGFEM approach to solve the hyperelastic energy minimization problem.

What all points I should consider to evaluate these methods ?


问:
One dimensional differential equation with resonance recombination (RHS)
I am trying to solve the following equation:
$$\dfrac{\partial }{\partial \epsilon}\left[B\left(\epsilon\right)f_e\left(1-f_e\right)+D\left(\epsilon\right)\dfrac{\partial f_e}{\partial \epsilon}\right] = R\left(f_e, f_h, E\right)$$

Till now, there is nothing special at all, and I can solve this (with a simple $R$).

But, the most important part of this equation is a form of recombination $R$:

$$R \approx \dfrac{1}{A-B\cdot f_e}$$

$A,B$ - are known constants and $(A-B\cdot f_e)$ should always be positive!

But numerical scheme which I use does not know anything about it, and it simply jumps through this boundary, which means no special recombination is included.

Typical value of $(A-B\cdot fe) \in [ 10^{-14},10^{-7}]$

Are there any methods to solve it?
The main problem is resonance characteristic of recombination $R$.

I trying to change $A$ on each iteration and wait for the convergence, after change boundary parameters (which is a parameter to change exactly) and solve again. I perform this procedure until $|A - A_\text{const}| = 1e-4$.

(This method works bad because of the relaxation coefficient within Newton iteration ~ 1e-6.)


问:
Computing Direct Scattering Transform
I'm working on the Nonlinear Schrodinger equation (NLSE) in 1d:
$$i\psi _t (t,x)+ \psi _{xx} + |\psi|^2\psi = 0 \, ,$$ for $t\geq 0$ and $x\in \mathbb{R}$.

This equation is integrable, and so similarly to the KdV equation, there is a way to solve it using the Direct and Inverse Scattering Transform.

What I'm looking for is an algorithm for computing the Direct Scattering Transform (DST), i.e., given $\psi(t,x)$, compute the discrete and continuous spectrum. The help could be in the form of a tutorial, a paper, pseudo-code or even an implementation.

Thank!

Title: Computing Direct Scattering Transform


问:
Mesh error quantification using error norm for a sharp logistic function
Consider the following logistic function in 1D for the domain $0$ to $1$.

$$ f(x) = \frac{1}{1 + \exp(-50(x-0.5))} $$

Let's say I resolve this on two meshes:

first with 11 equispaced points (0.0, 0.1, 0.2, ... 1.0)

second with 101 equispaced points (0.0, 0.01, 0.02, ... 1.0)

The plots for the two meshes look like this (blue is first, green is second):

How do I compute the error norm for the two meshes? I don't know if L2 norm is the best choice but let's say it is. In that case, the discrete version of the L2 norm is usually stated as follows:

$$ ||\mathrm{error}||_{2}
^{2} = h \sum_{i} |f_{i,\mathrm{ref}} - f_{i,h}|^2 $$

Where $f_{i,h}$ is function value for a given mesh at mesh point $i$, and $f_{i,\mathrm{ref}}$ is the value of the analytical function used as a reference.

However this will produce zero error norm for both cases, since the function evaluated at the mesh points is exact, i.e., it comes directly from the analytical expression. But we can clearly see in the plot that it's a wrong answer for the mesh with 11 points, i.e., error should be non-zero around the middle of the domain, resulting in a non-zero error norm. The zero-error answer is wrong for 101 points too but in the plot the error is not very visible (but it's still there, due to discretization/sampling). In other words, the second mesh has reduced the error compared to the first mesh. (what I described is beginner stuff; just putting things in context).

For correct error computation, a possible alternative is to compute the norm through quadrature. But before we even think about quadrature, we need to compute the error function itself. But how do we do that? Do we use a high resolution mesh (like 1001, or 10001 points) where we resolve $f_{i,\mathrm{ref}}$ as well as interpolate $f_{i,h}$ (say we use linear interpolation)? but this feels very sloppy/adhoc. Is there a better, more systematic way?

(I'm just trying to find out what are the "standard" ways of doing such error analysis for mesh resolution, such that I can show that the error in the second case goes down compared to that in the first case. Any tips for doing similar analysis in 2D on unstructured meshes would also be highly appreciated).


问:
Using C/C++ for Markov chain Monte Carlo (MCMC) methods
I'm working on optimizing the parameters of a mathematical model to fit experimental data, using an existing formula for the likelihood of observing the data given a set of parameter values. At the moment I use an algorithm (bound optimization by quadratic approximation, BOBYQA) to find the values that give the maximum likelihood, but I'd like to try using MCMC methods to integrate over the likelihood function, but haven't used them before. Can anyone recommend a good approach and/or existing library for implementing MCMC methods in C/C++?</p>


问:
Fourier transform spherical system
I need to take the Fourier transform of a 3D function $h(r)=h(|r|)$ so that I can invert a convolution problem.

What is the best way to do this with Python?

I know the the FT is equivalent to a sine transform (see below) but I can't find a Python library that appears to let me do this integral. One issue I have is that although my sample points for $h(r)$ are regularly spaced they do not start from $r=0$ due to how they've been generated, does this cause any issues?.

$$\tilde{h}(k)=\frac{2}{k}\int_0^{\infty}rh(r)\sin(2\pi kr)dr$$

Title: Fourier transform spherical system


问:
Neumann-Neumann boundary intersection

I am trying to solve a Vertex Centered, Finite Difference, Poisson equation $-\nabla^{2}u=f$ with Dirichlet boundaries on the left and bottom and Neumann boundaries on the top and right using Geometric Multigrid. An illustration is shown in 2-D (apologies for the hand drawn figure). Crosses indicate unknowns, red squares indicate ghost points for Neumann boundary, D indicates Dirichlet, 9-pt stencil is shown in green lines.  It is given that $\frac{\partial u}{\partial n}=0$ at Neumann boundaries.

If I use a 9-pt stencil in 2-D to approximate variable $u$ at various mesh points, what will be the value of the point at the Green questions mark?

Should I be using a 5-pt stencil at the Neumann-Neumann boundary intersection?

Will using a 5-pt stencil at the corner point not deteriorate the convergence rate in Multigrid?

If I use the 9-pt stencil at the top-left boundary, are we not taking a Dirichlet point which is outside the physical domain ?

Title: Neumann-Neumann boundary intersection


问:
Computing IFFT for only $k$ samples in a high dimensional signal
I have a very high dimensional signal, say $15$ dimensions. Across each dimension the width is $N$ points. So total number of points is $N^{15}$. I already have the FFT given to me. Only low frequency coeffs are non zero. The non zero coeffs are the ones inside of the hypercube of width 5. So only $5^{15}$ are non zero.

Now I want to take an IFFT to compute samples in signal domain.But I dont want to compute all samples of signalbut only a $k$ of them. $k$ is much smaller compared to $N^{15}$. The samplesin signal domain are by no means sparse, but its just that I want to compute only a small number $k$ of them. These $k$ samples are distributed arbitrarily and are not confined to any region.
How can I compute them efficiently?

I came across Sparse FFT, but here my signal is not sparse but just that i am interested only in a few samples in signal domain. So I am not sure I can use SFFT.

PS : I don't want to take a full IFFT due to memory constraints.


问:
Split-step-method for coupled equations
I have implemented a split-step-method for an equation of the shape
$$\partial_z E = i\partial_x^2E+ic|E|^2E$$
resulting in a split into the linear part
$$L=\partial_x^2$$
and the nonlinear part
$$N=c|E|^2E$$
giving me
$$E(x, z+\Delta z)=\exp(i\Delta z(N+L))E(x, t)$$
But now my $c$-variable is non-constant, but rather
$$\partial_tc=|E|^4$$
How do I couple that to my equation above in the split-step approach (either using an FFT or a DHT)? I can not include it as nonlinear term, after it does not depend on $E$, only on the absolute value.</p>

Title: Split-step-method for coupled equations


问:
Solving complicated coupled ODE using RK4/ODE45 in Matlab
I have the following coupled differential equations also known as Guiding Center Approximation. It is used to explain the position- and velocity change of particles (electrons and protons, N = 1000) moving in plasma.
I am trying to either solve this using the Runge-Kutta 4th order method or the built-in ODE45 in Matlab.

The equations are as follows:

\begin{align}
  {d u_{||} \over {dt}} &= {q \over m}  E_{||} - {\mu \over {\gamma m B}} (\mathbf{B} \cdot \nabla) B \\
  {{d \mathbf{R}} \over dt} &= {u_{||} \over {\gamma B}} \mathbf{B} + {{\mathbf{E} \times \mathbf{B}} \over B^2}
  + { \mu \over {\gamma q B^2}} \mathbf{B} \times \nabla B + {{ m u^2_{||}} \over {\gamma q B^4}} \mathbf{B} \times
  (\mathbf{B} \cdot \nabla) \mathbf{B}\\
  \mu & = {{m u^2_{\perp}} \over {2B}} = \mathrm{constant}
  \end{align}

where,

|| is the velocity component parallel to the magnetic field

$u$ is the particle velocity

$q$ is the charge

$m$ is the mass of the particle

$\mathbf{E}$ is the electric field vector

$\mu $ is the constant determined in the last equation above

$\mathbf{B}$ is the magnetic field vector

$\mathbf{R}$ is the position of the particle

$\gamma $ is the Lorentz factor given by $\gamma = (1-\frac{v^2}{c^2})^{-\frac{1}{2}}$

I am given the initial position and velocity of each particle. Furthermore I'm given $\mathbf{B}$ & $\mathbf{E}$. I'm also given the initial dt.
When $B$ is used without being bold it's either the norm of $\mathbf{B}$ or the component going in the direction of $\mathbf{B}$ and same goes for the electric field.

\begin{align}
\mathbf{B} & = {B_0 \over L} (x, y, -2 z)\\
\mathbf{E} & = E_0 L \left( { -y^2 \over {(x^2+y^2)^{1.5}}}, {xy \over
     {(x^2+y^2)^{1.5}}}, 0 \right)
\end{align}

To solve these coupled equations I tried using the RK4 method. What I'm unsure of is how to write my equations in a way of $f(t_n,y_n)$

where,

$k_1 = f(t_n,y_n) $

$k_2 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot k_1)$

$k_3 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot k_2)$

$k_4 = f(t_n + h, y_n + h \cdot k_3)$

$y_{n+1} = y_n + \frac{h}{6} \cdot (k_1 + 2\cdot k_2 +2\cdot k_3 + k_4 )$

$t_{n+1} = t_n + h$

And furthermore, how do I using Matlab, calculate these two equations at the same time (as I assume I must). Lastly, would it be easier or more accurate to use the builtin ODE45?


问:
Finding minimum of Wronskian determinant
If I have a partial differentail equation such as $\frac{\partial^2 u}{\partial t^2} + c^2 \frac{\partial^2 u}{\partial x^2}$, with boundary conditions $u(0,t) = 0$ and $u(1,t)=0$, I can solve this system using separation of variables: $u=X(x)T(t)$. This means I have an ODE: $X''(x) + (c^2-\lambda)X(x)=0$, with $\lambda$ as separation constant.

To look for solutions I set up the Wronskian and search for $\lambda$ values for which the Wronskian determinant is zero, meaning I found an eigenvalue. For 'easy'boundary conditions this can also be done manually, however I want to do it numerically.

Now I look over a range of $\lambda$, look for minima, and optimize around those minima.

Knowing how the Wronskian is set up, is there an optimization method or numerical trick I can perform to speed up the process?


问:
How to fast estimate derivates for calculating quantiles
I would like to know if there exists a package or how one can fast calculate the quantiles of a function within python, where the inverse of the function for calculating the quantile depends on the first and second derivative. More concretely.

Suppose we are given a $C^2([a, b])$ function $f$ and define the curvature $\kappa$ as

$$ \kappa(x) = \frac{|f''(x)|}{|1+(f'(x))^2|^{\frac{3}{2}}}$$

This is a positive function. I then normalize $\kappa$ s.t. $\int_a^b \kappa = 1$. For simplicity let's denote the normalized function with $g(x):=\frac{\kappa(x)}{\int_a^b \kappa}$.

Now I define the cumulative distribution function (CDF) $F$ as:

$$ F(x) = \int_a^x g(t) dt $$

we know that $F(a) = 0$ and $F(b) = 1$. The problem I want to solve is:

For given value $y$ of $F$ find $x$ s.t. $F(x) = y$, i.e. the inverse problem of $F$.

I'm interested in a solution within python (or any other similar language like Matlab, R). In my problem I only have data points of $f$ and need to approximate it. The following code snippet is an example for the $sin^2(x)$ curve.

```<code>import timeit
import numpy as np
from scipy.interpolate import UnivariateSpline
import scipy.integrate as integrate
from scipy.optimize import brentq

x = np.linspace(0, 6, 100)
interp = UnivariateSpline(x, np.sin(x)**2, k=4, s=0)

# calculate the normalization constant

def kappa(x):
    return abs(interp.derivative(2)(x))/abs((1+interp.derivative(1)(x)**2)**(1.5))

k = integrate.quad(kappa, 0, 6)[0]

def g(x):
    return kappa(x)/k

x = np.linspace(0, 6, 100)
interp2 = UnivariateSpline(x, g(x), k=4, s=0)

def F(x):
    return integrate.quad(interp2, 0, x)[0]

# get 20 test points
# not 0 and 1 are clear so we don't need
# to estimate these points.

N = 20
x2 = [(i - 1)/(N-1) for i in range(2, N)]

# first approach would be
start_time = timeit.default_timer()
y = [brentq(lambda u: F(u) - i, 0, 6) for i in x2]
elapsed = timeit.default_timer() - start_time
print(elapsed)

# just add the two default points

# y = [0] + y
# y = y.insert(len(y), 6)

def F2(x):
    return integrate.quad(g, 0, x)[0]

start_time = timeit.default_timer()
y = [brentq(lambda u: F2(u) - i, 0, 6) for i in x2]
elapsed = timeit.default_timer() - start_time
print(elapsed)
</code>```

Note, I've calculated the inverse with two different methods. First I use another interpolation to approximate $\kappa$. In the second version I run it via directly using $\kappa$ where it has to evaluate the first and second derivative. If I'm running this on my pc I get the elapsed times:

```<code>0.449614048004
7.88159608841
</code>```

What a huge different! Why is the second one so much faster? Now given what I want to do what is the best way to speed this up (even more than in case 1).

EDIT I've actually had a bug in the code. Since I've used python 2.7 I needed to import the division module from future. With this I get the following. Time is still a huge difference and values are very close for both system:

Time and values for first method:

```<code>2.38353395462
[0.1437681151139289, 0.35456065909107926, 1.1539903423627926, 1.3971602056486263, 1.5460856607363669, 1.686478225698852, 1.8751150788777717, 2.6398043923321652, 2.9355086219617466, 3.091753998559681, 3.230312179133905, 3.4029515812391478, 4.058035698689046, 4.470426379024948, 4.636764697040471, 4.774960068200856, 4.93583175345252, 5.273899782448736]
</code>```

Time and values for second method

```<code>57.0844209194
[0.1437734552194811, 0.354572374885639, 1.1530574841452361, 1.3967524608790778, 1.545758378815468, 1.6861161414005308, 1.8745107127405394, 2.637806584725155, 2.9348568464428215, 3.0912665934557833, 3.229804581185758, 3.4021836211797964, 4.0530805425833485, 4.469683619221007, 4.636255181639525, 4.774458238739285, 4.935134505525859, 5.271168509499432]
</code>```


问:
Choice of velocity grid - staggered or not?
I'm trying to understand when and why one would use a staggered vs. a colocated grid in problems that have velocities and scalars that they transport (e.e. density).

If scalars are defined cell-centered, then in the physics world you often hear that velocities should be defined between cells, because apart from being intuitive, the resulting numerical schemes are stable.

I'm certain to having had heard this statement in the context of finite difference codes.

But now I'm learning finite volume methods with the Riemann Solvers by E. Toro. There, in the math world, vectors like velocities are just 3-scalars in equation systems that are all solved on one and the same colocated grid. There is never any mention of staggering grids at all (at least not as far as I got in the book, around Chapter 16, 3rd Ed.).

So is the choice of defining velocities on staggered vs. colocated grids related to the numerical methods one is using, or rather if one is solving ODEs vs. PDEs, or is the answer something completely different?


问:
the augmented global stiffness matrix is not positive semi-definite using Lagrange Multipliers method within FEM
The augmented global stiffness matrix is not positive semi-definite when using Lagrange Multipliers method to enforce boundary constraints on a simple square domain of integral form:<img src="https://i.stack.imgur.com/GofAL.png" alt="enter image description here">
I am considering linear elastic material constitutive.
I checked the condition number which is order 10e3 indicating the augmented matrix is not much ill-conditioned. However, there are some negative eigenvalues for the matrix. What are possible reasons for the lack of semi-definite positiveness? It should be always positive semi-definite for the (linearized) global stiffness matrix whatever Lagrange multipliers method is employed. Is my Lagrange multipliers discretization/approximation space is not well set up or I lost something very important? Anyone can share some experience? </p>


问:
Uniform distribution in 3D space
Posted this at math stack exchange as well, but alas no replies! So, I have been trying to find ways of distributing particles of spherical or other shape in 3D space, e.g. rectangular space. Random() and checking for overlaps will only get you up to a certain volume fraction percentage, then it is nearly impossible to find a non-overlapping position. I know there is a way to "drive" particles in empty spaces by some sort of Monte Carlo distribution or using some kind of potential difference? A dirty example in pseudo code or Matlab would be awesome!</p>

Title: Uniform distribution in 3D space


问:
Do practice and theory differ substantially when implementing Neumann Boundary Conditions using a Mixed Method?
I have implemented a pretty straightforward finite element solver for the following Poisson equation.  For the purposes of this question we can assume the source term and the Dirichlet data both vanish.  I assume that $g\in H^{-1/2}(\partial \Omega_N) \cap L^2(\partial \Omega_N)$ (instead of just $g\in H^{-1/2}(\partial\Omega_N)$)

\begin{aligned}
\mathbf{q} + \mathbf{\nabla}u &= 0 & \mathrm{in\,} \Omega,\\
\mathbf{\nabla} \cdot \mathbf{q} &= 0& \mathrm{in\ } \Omega,\\
u& = 0 &\mathrm{on\ } \partial\Omega_D,\\
\mathbf{q}\cdot \mathbf{\eta} &= g & \mathrm{on\ } \partial\Omega_N.
\end{aligned}

I accounted for Neumann boundary conditions on all or some of the boundary by using Lagrange multipliers as follows.

Let the The domain $\Omega$ have a triangulation $T$ made up of simplexes $K$.
  We denote $$\mathbf{W}_h= \{\mathbf{p}\in \mathbf{H}^{\mathrm{div}}(\Omega)\colon \forall K\in T \quad \mathbf{p}\big|_K \in \mathrm{RT}^k(K))$$
$$V_h = \{v\in L^2(\Omega)\colon \forall K \quad v\big|_K \in P^k(K)\}$$
$$X_h = \{\mu \in L^2(\partial \Omega_N)\colon \forall K\in T\,\, \forall \mathrm{edge}\subset\partial K\cap \partial\Omega_N\quad \mu \big|_{\mathrm{edge}} \in P^k(\mathrm{edge})\}$$

Where $\mathrm{RT}^k(K)$ is the standard Raviart Thomas space and $P^k(K)$ is the space of polynomials of total degree no more than $k$.

The problem is to find the triple $(\mathbf{q}_h,u_h,\lambda_h) \in \mathbf{W}_h \times V_h \times X_h$ so that for all test functions $(\mathbf{p}, v, s) \in (W_h \times V_h \times X_h)$ we have
\begin{aligned}
0 &= \int_\Omega \mathbf{p}\cdot \mathbf{q}_h\, \mathrm{d} x - \int_\Omega (\mathbf{\nabla} \cdot \mathbf{p}) u_h\, \mathrm{d} x + \int_{\partial\Omega_N} \mathbf{p} \lambda_h \cdot \mathbf{\eta}\, \mathrm{d} x \\
0 &= \int_\Omega v (\mathbf{\nabla} \cdot \mathbf{q}_h)\, \mathrm{d} x\\
\int_{\partial\Omega_N}s g \, \mathrm{d}x&= \int_{\partial\Omega_N} s\mathbf{q}_h \cdot \mathbf{\eta}\, \mathrm{d} x
\end{aligned}

The scheme works (well defined, converges to true solution, etc.) and is (I think) close to standard practice. However, in all the references I have seen in the literature (On the Mixed Finite Element Method with Lagrange Multipliers Babuska and Gatica, 2003 Numerical Methods for Partial Differential Equation and section 4.4 of A simple Introduction to the Mixed Finite Element Method  by Gatica) require that the trial and test space for the Lagrange multiplier, $X_h$, have elements that are continuous across elements of the "boundary mesh".

Is this scheme close to standard practice?

Are there any references I could read that discuss this scheme?

What problems am I allowing by choosing a piecewise discontinuous space for the Lagrange multiplier?

A disadvantage of the scheme is that it requires extra regularity on $g$.  Are there more disadvantages?


问:
An algorithm to turn descending order into ascending order, but why correct?
I wonder why does this algorithm below work well? Many appreciates!

Task:

Suppose all elements in an array with length $n = 2^m$ are in descending order: $a_1 > a_2 > ... > a_n$, we need to turn them into ascending order.

Algorithm:

For each step from $1$ to $m$:

Swap each adjacent (odd index, even index) pair into ascending order, so that:
$a_1 \le a_2, a_3 \le a_4, ... , a_{n-1} \le a_n$

Move all elements with odd index to the front of those
with even index, so $a_1, a_2, a_3, a_4, ..., a_{n - 1}, a_n$ turns to : $a_1, a_3, ..., a_{n - 1}, a_2, a_4, ..., a_n$

Example:

Sort result after each step with $n = 8$:

[0] 8 7 6 5 4 3 2 1

[1] 7 5 3 1 8 6 4 2

[2] 5 1 6 2 7 3 8 4

[3] 1 2 3 4 5 6 7 8

Code:

```<code>#include <stdlib.h>

void swap(int *arr, int i, int j)
{
    int temp = arr[i];
    arr[i] = arr[j];
    arr[j] = temp;
}

void magic_sort(int *arr, int *result, int n)
{
    int half = n >> 1;
    int i;

    /* Swap adjacent (odd index, even index) pair into ascending order */
    for (i = 1; i <= n; i += 2)
    {
        if (arr[i] > arr[i + 1])swap(arr, i, i + 1);
    }

    /* Move all elements with odd index to front half */
    for (i = 1; i <= n; i += 2)
    {
        result[(i + 1) / 2] = arr[i];
    }

    /* Move all elements with even index to tail half */
    for (i = 2; i <= n; i += 2)
    {
        result[half + i / 2] = arr[i];
    }

    for (i = 1; i <= n; i++)
    {
        arr[i] = result[i];
    }
}

int main()
{
    int m = 3;
    int n = 1 << m;
    int *arr = (int *)malloc(sizeof(*arr) * (n + 1));
    int *result = (int *)malloc(sizeof(*arr) * (n + 1));
    int step, i;

    printf("[%d] ", 0);

    for (i = 1; i <= n; i++)
    {
        arr[i] = n + 1 - i;

        printf("%d ", arr[i]);
    }

    printf("\n");

    for (step = 0; step < m; step++)
    {
        magic_sort(arr, result, n);

        printf("[%d] ", step + 1);

        for (i = 1; i <= n; i++)
        {
            printf("%d ", arr[i]);
        }

        printf("\n");
    }

    free(arr);
    free(result);

    return 0;
}
</code>```