问:
How to approach geographic data interpolation by distance?
let's say I have a set of geographic locations (lat, lng) resulting from a query. Those locations have some kind of internal ranking, my set is sorted by this number in a descending order.

Now I'm trying to interpolate this data by eliminating all locations in a certain distance from each other. For example, starting with the highest ranking, point A must not have any adjacent points closer than 300 meters, so point B must be in distance of at least 300 meters away.

First thing came to mind was creating some kind of chess table which contains all distances and odd the low ones out. This may work for 5, 10, maybe 100 points but obviously won't scale well. I expect an average of 10000 points to be processed regularly. So I would have to query and post-process 9999 times for each point existing (I know this could be optimized minimally but you get the idea).

Second idea is processing the data like a tree, starting with the highest ranking point(s) and working my way down. Like creating an incremental search 'map', each 'allowed' point adding a circle of 300m radius as 'disallowed' area used as template for the next tier. I have used the haversine formula before and understand how to achieve this but I'm not sure performance-wise and fear eliminating too many points from my set.

So my question is, how to approach this with a goal of minimizing cost? Is there maybe an algorithm targeting this problem I don't know about yet? Help appreciated.


问:
Evolutionary dynamics in vascularised tumors, PDE-ODE coupled system
I have to solve the following PDE-ODE system

$$ \displaystyle{\partial_{t} n = \bigl[a(s) - b(s)(y - h(s))^{2} -  d\int_{\mathbb{R}}  n \, dy \; + \; \beta \, \partial^2_{yy} n \;}
\\\\
\displaystyle{s' = -\int_{\mathbb{R}} \bigl[\eta_{s}(g(y,s))_{+}n \bigl] \, dy - \lambda_{s}s + S0 \;} $$
where $n = n(t,y)$ is the distribution of tumor cells depending on time and phenotype state $y$, $s = s(t)$ the oxygen concentration with intial oxygen concentration $S_{0} = 6.3996 \times 10^{-7}$ and $(.)_{+}$ the positive part. The phenotypic state is $y \in \mathbf{R}$, the time $ t \in [0,+\infty]$ . The function $a(s),b(s),h(s)$ depending on the oxygen concentration are so defined:
$$ a(s) := \Bigg( \gamma_{s} \frac{s}{\alpha_{s} + s} + \frac{\varphi^{2}}{\varphi + \gamma_{s}\dfrac{s}{\alpha_{s} + s}} \Bigg) $$
$$ b(s) := \bigg(\varphi + \gamma_{s}\dfrac{s}{\alpha_{s} + s} \bigg) $$
$$ h(s) := \frac{\varphi}{\varphi + \gamma_{s}\dfrac{s}{\alpha_{s} + s}} $$
the function $g(y,s)$ depending on phenotype and oxygen concentration as
$$g(y,s) := \gamma_{s} \frac{s}{\alpha_{s} + s}(1 - y^{2})$$

The parameters are the following
$$\alpha_{s} \ \text{Michaelis-Menten constant of oxygen} \ 1.5 \times 10^{-7}$$
$$d \ \text{Rate of death due to competition for space} \ 2 \times 10^{-13}$$
$$\varphi_{} \ \text{Maximum rate of cell division via anaerobi pathway}  \ 1 \times 10^{-5}$$
$$ \gamma_{s} \ \text{Maximum rate of cell division via aerobic pathways} \ 1 \times 10^{-4} \ $$
$$ \eta_{s} \ \text{Conversion factor for cell consumption of oxygen} \ 2 \times 10^{-11} \  $$

My question is how to solve this problem numerically (I'm using Matlab). To start, I used a uniform grid for the domain phenotype $y∈[−2L,2L]$ with $L= 2,Ny= 1200$ points (https://biologydirect.biomedcentral.com/articles/10.1186/s13062-016-0143-4).  After that I also discretized time $t$ between the starting time $t=0$ and the final time $t_{f} = 5$ in units of $10^{4}s$ with a time step $dt= 1s$ (and $N_{t} $time step), according with the CFL condiction.

The method for calculating numerical solutions is based on a time splitting scheme between the conservative part and the reaction term.  I approximate the diffusion term through finite difference approximation (three point stencil) and we used an explicit finite difference scheme for the reaction term. The scheme that I implemented is the following
$$n^{k+1}_{i} = n^{k}_{i} +  dt \, \bigl[ a(s^{k}) - b(s^{k})(y_{i} - h(s^{k}))^{2} - d \, dy \, \sum_{j=1}^{N_{y}}(n^{k}_{j})\bigl] \,n^{k}_{i} + \frac{\beta dt}{dy^{2}} \biggl[ n^{k}_{i+1} - 2n^{k}_{i} + n^{k}_{i-1}\biggl] $$
$$ s^{k+1} = s^{k} -  dt \, dy \,  \sum_{j=1}^{Ny}(\bigl[\eta_{s}(g(y_{j},s^{k}))_{+}n^{k}_{j} \bigl]) - dt \, \lambda_{s} \, s^{k} + dt \, S0 $$
$$k = 1,...,N_{t}-1, \ i = 2,...,N_{y}-1.$$
But unfortunately the solution isn't correct and oxygen concentration has very abnormal behavior

Here is the documentation on which this work is based
https://arxiv.org/abs/1910.08566
This is an easy case in which there is not a spatial dependence.

Here i upload the code:

```<code>    clc
    clear all
    close all

    set(0,'DefaultAxesFontName', 'Times New Roman')
    set(0,'DefaultAxesFontSize', 18)
    set(0,'defaultaxeslinewidth',1)
    set(0,'defaultpatchlinewidth',1)
    set(0,'defaultlinelinewidth',4)
    set(0,'defaultTextInterpreter','latex')

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%% To solve
    %%%%% n_t = beta*n_yy + (a(s) - b(s)*(y - h(s))^2 - chi*rho)*n
    %%%%% s_t = ...
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Parameters
    betay = 1e-6; % rate of spontaneous phenot. variations
    chi = 2*1e-13; d = chi; % rate of death for competition
    phi = 1e-5;
    alfac = 2*1e-6;
    alfas = 1.5*1e-7;
    muc = 4.4*1e-6;
    mus = 2*1e-5;
    gammac = 1.8*1e-4;
    gammas = 1e-4;
    etac = 4*1e-11;
    etas = 2*1e-11;
    lambdac = 2.3*1e-4;
    lambdas = 2.78*1e-6;

    %% Time, Space and Phenotype
    Ly = 4;
    Ny = 1200;
    t0 = 0;
    tf = 5*1e4;

    %% Discretization
    ysp = linspace(-Ly,Ly,Ny)';
    dy = ysp(end-1) - ysp(end-2);
    dt = 1;
    tsp = t0:dt:tf;
    Nt = length(tsp);

    %% Oxygen concentrarion at time 0, Oxyeng into system
    S0 = 6.3996*1e-7;

    %% Def : a,b,h;
    a = @(s) gammas*(s./(alfas + s))  + + (phi).^2./(phi + gammas*s./(alfas + s));
    b = @(s) phi + gammas*(s./(alfas + s));
    h = @(s) (phi)./(phi + gammas*s./(alfas + s));
    g = @(s) (gammas*s./(alfas + s)).*(1 - ysp.^2);

    %% Numerical approach, Initialization
    Rho0 = 1e8; Mu0 = 0.5; Vu0 = 1; % Data from reference
    n0 = Rho0*sqrt(Vu0/(2*pi))*exp(-(Vu0/2)*(ysp - Mu0).^2).*ones(1,1)';
    R = zeros(Ny,1);    % initialization
    nOld = n0;          % initialization
    nNew = n0;          % initialization
    sOld = S0;
    sNew = S0;
    sMemory = S0*ones(Nt,1); % To memorize s(t)
    %% I-M-V
    rho = ones(Nt,1).*Rho0;
    mean = ones(Nt,1).*Mu0;
    var = ones(Nt,1).*Vu0;
    initial = [1/Vu0,Mu0,Rho0];

    for kk = 1:Nt-1 %%% TIME
        tk = kk*dt; % time counter
        %% Explicit-Eu
        %%%% Reaction and Integral for n
        R = a(sOld) - b(sOld).*(ysp - h(sOld)).^2 - d*dy*sum(nOld);
        %%%% New n
        nNew(2:Ny-1) = nOld(2:Ny-1) + dt*R(2:Ny-1).*nOld(2:Ny-1) + ...
                + (betay*dt/(dy^2))*(nOld(3:Ny) - 2*nOld(2:Ny-1) + nOld(1:Ny-2));
        %%%% Neumann boundary cond
        nNew(1) = nNew(2);
        nNew(Ny) = nNew(Ny-1);
        %%%% rs
        rs = etas*max(0,g(sOld));
        sNew = sOld - dt*(dy*sum(rs.*nOld)) - dt*lambdas*sOld + dt*S0;

        %% Int - Mean - Variance
        rho(kk+1) = dy*sum(nNew);
        mean(kk+1) = dy*sum(ysp.*nNew)./rho(kk+1);
        var(kk+1) = dy*sum(ysp.^2.*nNew)./rho(kk+1) - mean(kk+1).^2;

        %% Updating
        nOld = nNew;
        sOld = sNew;
        sMemory(kk+1) = sNew;
    end

    %% I-M-V with ODE Plot
    figure()
    subplot(1,3,1)
    plot(tsp,rho,'--')
    % axis([t0 tf min(min(rho))-.1e7 max(max(rho))+.1e7])
    axis square
    xlabel('time')
    title('<span class="math-container">$\rho(t)$</span>')
    subplot(1,3,2)
    plot(tsp,mean,'--')
    % axis([t0 tf min(min(media))-.1 max(max(media))+.1])
    axis square
    xlabel('time')
    title('<span class="math-container">$\mu(t)$</span>')
    subplot(1,3,3)
    plot(tsp,var,'--')
    % axis([t0 tf min(min(var))-.1 max(max(var))+.1])
    axis square
    xlabel('time')
    title('<span class="math-container">$\sigma^{2}(t)$</span>')
    s = sprintf('I-M-V of solution n(y,s) with tf = %d, dt = %f',tf,dt);
    suptitle(s)

    %% s,a,b,h,g
    figure()
    subplot(1,5,1)
    plot(tsp,sMemory)
    title('S(t)')
    xlabel('time')
    axis square
    subplot(1,5,2)
    plot(tsp,a(sMemory))
    title('a(S(t))')
    xlabel('time')
    axis square
    subplot(1,5,3)
    plot(tsp,b(sMemory))
    title('b(S(t))')
    xlabel('time')
    axis square
    subplot(1,5,4)
    plot(tsp,h(sMemory))
    title('h(S(t))')
    xlabel('time')
    axis square
    subplot(1,5,5)
    G = (gammas*sMemory./(alfas + sMemory)).*(1 - ysp'.^2);
    plot(tsp,G(:,1))
    title('g(S(t))')
    xlabel('time')
    axis square
</code>```

ContentLicense: CC BY-SA 4.0


问:
Why GA convergence curves continue as two parallel lines?
I'm working on a optimization problem and using GA algorithm (in MATLAB, ga function).

As you know MATLAB plots GA result with two curves, one for the best values and other to show the mean values and when this two curves touch each others it means algorithm has been converged.

On my cause these two curves don't get along and continue as two parallel lines until the end of max Generations number and finally a premature convergence,What makes it happen?

I have tested with other parameters value but same result. On pre-processing I used Missing value,Normalized and smoothing and tested data set belongs to UCI.

GA Parameters:

```<code>MutationFcn      : @mutationadaptfeasible: 0.03
CrossoverFcn     : @crossoverheuristic : 1.2
maxGenerations   : 200;
populationSize   : 180;
Chromosome length: 33
Test Instances   : 71
</code>```


问:
polylog implementation of fully-dynamic graph connected-components?
I have read about papers in the last 20 years that have solved this problem. Many are mentioned in http://jamiemorgenstern.com/teaching/s18-6550/notes/notes-lec4-dgc.pdf

Unfortunately the only implementations I have found are also 20 years old and no longer compile due to changing LEDA and C++ (from the ACM DL).

Can anybody suggest any implementations that are open source, poly logarithmic time complexity, and support fully dynamic connected components? By that I mean the capability to add and remove edges without recalculating everything. Polylog seems to be the best possible so that is what I am hoping to find.

Some example papers that sounded promising:

https://u.cs.biu.ac.il/~liamr/p723-holm.pdf

https://people.csail.mit.edu/rahul/papers/dyncon-jea2001.pdf

https://www.semanticscholar.org/paper/An-Empirical-Study-of-Dynamic-Graph-Algorithms-Alberts-Cattaneo/d3efbed24e526b42da5d001036519c9b0fa79856


问:
CFD and finite volume method: Dirichlet boundary conditions for the Euler equations
Please point me to an answer if one already exists, but after some searching, I still can't find the answer to what seems like a very simple question. There are plenty of references out there for solving the Euler equations with FVM. I am using the Toro (1997) book, Riemann Solvers and Numerical Methods for Fluid Dynamics, as a reference on Riemann solvers and the Godunov method.

I am attempting to solve the 2-dimensional Euler equations on a fixed, uniform grid of cells using the Godunov method and HLLC flux. Having implemented one of the Sod shock tube tests in the book, I feel confident that I haven't made a blatant error in my Riemann solver. The book does not address Dirichlet boundary conditions, however, which is exactly the boundary condition I'd like to implement.

Equations

I have the 2-dimensional Euler equations with a source term:

$$ \frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho \mathbf{u}
 \\ \rho E \end{pmatrix} + \nabla \cdot \begin{pmatrix} \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbf{I} \\ \left(\rho E + p\right) \mathbf{u} \end{pmatrix} = \mathbf{S} $$

where I have used bold notation for vectors, $\rho$ is the fluid density, $\mathbf{u} = \left(u, v\right)$ is its velocity, $p$ is the pressure, and $\rho E = \frac{p}{\gamma - 1} + \frac{1}{2} \rho \left(u^2 + v^2\right)$, where $\gamma$ is the (constant) ratio of specific heats.

Code

I tried to write my finite volume code in a general way so that I could easily test it and apply it to different kinds of problems. This means that each subproblem is expected to provide a function with a specific interface for each of the north, south, east, and west boundaries. To prescribe the Dirichlet boundary condition I want, there is a layer of ghost cells surrounding the mesh where I prescribe the fixed value $\mathbf{g}_D$. Then, I compute the flux as I normally would, specifying the fixed state in the appropriate slot for the Riemann solver. The resulting flux is then returned to the driver function, which just accumulates all the flux, weighted by the inverse grid spacing, into a big vector of derivatives.

Problem

The prescribed value on the ghost cells is easy to preserve by messing with how the source term is added, etc. Assuming the initial conditions satisfy $\mathbf{g}_D$, I can just set the total flux and source term to zero on the boundary so that the derivative there is zero. This made sense when I did it, but I'm no longer convinced it's correct, since I now get rather serious spikes in velocity at the boundaries of the domain where the value was forced to be constant (image below). The velocity spikes eventually work their way into every other field.

Question

Am I on the right track by using ghost cells with a prescribed value? Clearly, the exact method I've described doesn't work as expected. Is there a better way to accomplish this?


问:
How to check if my stiffness matrix is correct
I built the stiffness matrix for the Poisson equation on a 2-dimensional domain with the shape of "almost" an octagon, using pyramid basis functions. I used almost to intend the fact that I have an "irregular" shape.

The mesh has been given to me by a point matrix $P$ and a connectivity matrix $T$. In $P$ I have the points of the domain, while in the matrix $T$ I have in each entry the vertices of the triangles. On page 47 of Larson & Bengzon book (section "Data Storage Structures") is explained better what I mean with a simple example.

Anyway, I implemented the assembly procedure as described on page 85 and I

obtain the right result for simple squared domains. I also tried to use the code for the $[0,1]$-interval with grid size $h$ and things are good.

Of course here the point matrix is made by a linspace from $0$ to $1$, while $T$ is a matrix like $$T =\begin{bmatrix} 1& 2&3  \\ 2&3&4\end{bmatrix} $$ and I end up with the usual Laplacian matrix (modulo a factor of $h$). So it's right.

Now the main question:

How can I check that the matrix I obtain for the "irregular" domain is correct?

I already checked it's symmetric, and since there were Homogeneous Dirichlet boundary conditions I followed the approach of Wolfgang Bangerth explained here at slide 14 (approach 2a). After this, I computed the condition number $K$ of the matrix and it is $K \approx 6$.


问:
How to use discrete cosine and discrete sine transforms in fftw
I work on fluid-related simulations. I have used FFT for fluid simulation. I want to use discrete cosine transform (DCT) and discrete sine transform (DST) to transform my velocity field to wavenumbers.

I am using REDFT00 and RODFT00 as DCT and DST respectively (so the inverse are themselves).

I understand that the normalization is $2(n-1)$ and $2(n+1)$, respectively.

for DCT and DST, I take

$u(x) = \cos(3x); n=64$

But when I invert it, values after $(n/2)+1$ are all incorrect.

What is the right way to do?

UPDATE: I was using an old library which was heavily modified(FFT worked fine though). I used the original library and everything works fine. But now I have difficulty in using the library for my purposes(Don't know how to take care of the wavenumbers) which I have posted as a separate question. So I guess I don't need the answer. It was a defective library.


问:
fastest way to compute many small dot products
I have two n-by-3 blocks contiguous in memory ("n vectors of length 3") and I'd like to compute the dot product between each of the rows as fast as possible. In numpy, using einsum is the fastest variant to the best of my knowledge (einsum("ij,ij->i", a, b)).

I've compared this to the transposed variant (also contiguous in memory), i.e., two 3-by-n blocks. Here, the each of the components (x, y, z) form a big block in memory rather than n contiguous 3-blocks. I found that computation here can be about twice as fast:

I'm out of memory at $2^{25}$ vectors so I'd be interested in results beyond that point.

Any idea on why the transposed variant is so much faster at that size range?

Code to reproduce the plot:

import numpy
import perfplot

def setup(n):
    a = numpy.random.rand(n, 3)
    b = numpy.random.rand(n, 3)
    aT = numpy.ascontiguousarray(a.T)
    bT = numpy.ascontiguousarray(b.T)
    return (a, b), (aT, bT)

perfplot.save(
    "rel.png",
    setup=setup,
    n_range=[2 ** k for k in range(1, 25)],
    kernels=[
        lambda data: numpy.einsum("ij, ij->i", data[0][0], data[0][1]),
        lambda data: numpy.einsum("ij, ij->j", data[1][0], data[1][1]),
    ],
    labels=["einsum", "einsum.T"],
    xlabel="len(a), len(b)",
    relative_to=1,
    flops=lambda n: n,
)
```


问:
Fast approximate evaluation of Fourier-Legendre series
Suppose I know that a function from $[0,\pi] \to \mathbb{R}$ may be written as

$$
\sum_{k=0}^\infty A_l \frac{2l+1}{4\pi} P_l(r)
$$

where $A_l$ all are known. Is there a way in which I may very quickly approximate the value of the function? My naive implementation is very slow and also suffers problems with numerical issues stemming from adding the very small numbers for certain $l$ and $r$ to large numbers.

A package or something in whatever language, Python, Julia, Fortran, etc, anything goes really, would be really appreciated.


问:
Subspaces for Iterative methods
In the original paper of Conjugate Gradients, the authors mention that if we pick the canonical basis $\{e_1,e_2,\ldots,e_n\}$, to obtain A-orthonormal vectors, we end up with the Gaussian elimination method (section $12$). My question is the following:

Apart from this canonical basis and of-course the Krylov subspace, are there other attractive subspaces to look at?

For instance, intuitively, it looks like if $A$ is well-conditioned, picking the subspace generated by $\{b,A^Tb,(A^T)^2b,(A^T)^3b,\ldots\}$ might be attractive option. If $A$ is the most well-conditioned matrix, the method is guaranteed to converge in $1$ iteration. Has there been some analysis done on other possible subspaces (apart from Krylov and the canonical ones)?

Title: Subspaces for Iterative methods


问:
Calculate integrals using numpy.fft
Good evening,

I would like to understand why I do not get the correct result:
I assume that I know my function on discrete data points and expand it as a discrete Fourier transform: $\text{sin}(x)=\sum_{k=-m}^m A_ke^{i2\pi kx}$

This I use to calculate an integral:
$$\int_{-R}^{R}\text{sin}(x+y)\text{d}x = \sum_{k=-m}^m A_ke^{i2\pi ky}\int_{-R}^R e^{i2\pi kx} \text{d}x = \sum_{k=-m}^m\frac{A_k}{i2\pi k}\left( e^{i2\pi kR}-e^{-i2\pi kR}  \right) e^{i2\pi ky}$$

In other words, I calculate $A_k$ from FFT(sin(x)), multiply each $A_k$ with the term $\frac{1}{i2\pi k}\left( e^{i2\pi kR}-e^{-i2\pi kR}  \right) $ and go back into real space with inverse FFT of these modified $A_k$.

Yet I do not get the correct result (the analyitc solution is $-\text{cos}(R+y)+\text{cos}(-R+y)$ ). What am I doing wrong?

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, fftfreq, ifft

pi2 = 2*np.pi

L = pi2
N = 1001   ## should be odd!!
R = 0.2

x = np.linspace(0, L-L/N, N)
f = np.sin(x)
FOU_f = np.zeros(N, dtype=complex)
phi = np.zeros(N, dtype=complex)
FOU_phi = np.zeros(N, dtype=complex)

FOU_f = fft(f, N)    ## function to fourier space

for k in range(1, N):
    FOU_phi[k] = ( FOU_f[k]/(1j*pi2*k) ) * ( np.exp(1j*pi2*k*R) - np.exp(-1j*pi2*k*R) )

phi = ifft(FOU_phi, N)  ## back to real space

analytic_sol = -np.cos(x+R) + np.cos(-R+x)

plt.plot(x, phi, label="test")
plt.plot(x, analytic_sol, label="analytic")
plt.legend()
plt.show()
```

Title: Calculate integrals using numpy.fft


问:
How do I calculate the amplitude after a 2D r2c transform using FFTW?

#include<stdio.h>
#include<fftw3.h>
#include<math.h>
#include <stdlib.h>

int main(){
    int N=2000;
    int i,j;
    FILE *filepointer;
    filepointer=fopen("2DDFT_spacetime.plt","w");
    //double in[N][N];
    double *in;
    fftw_complex *out;
    fftw_plan p;
    double fx=13.0;
    double fz=9.0;
    double x[N];
    double xstart=0.0;
    double xend=5.0/fx;
    double z[N];
    double zstart=0.0;
    double zend=5.0/fz;
    double dx=(xend-xstart)/(N-1);
    double dz=(zend-zstart)/(N-1);
    x[0]=xstart;
    z[0]=zstart;
    in = (double*) malloc(sizeof(double) * N * N); //allocates input array
    out = fftw_alloc_complex(N*((int)floor(N/2)+1)); //wrapper function ;allocates output array
    p = fftw_plan_dft_r2c_2d(N,N, in, out, FFTW_MEASURE);
    for(i=1;i<N;i++) {
        x[i]=x[i-1]+dx;
    }
    for(i=1;i<N;i++) {
        z[i]=z[i-1]+dz;
    }
    for(i=0;i<N;i++) {
        for(j=0;j<N;j++) {
            in[i*N+j]=cos(2*M_PI*fx*x[i]+2*M_PI*fz*z[j]);
        }
    }

    fftw_execute(p);

    fprintf(filepointer,"TITLE =\"FFTW\"\nVARIABLES=\"Wavenumber-x\"\n\"Wavenumber-z\"\n\"Amplitude\"\nZONE T=\"Amplitude\"\n I=%d, J=%d, K=1, ZONETYPE=Ordered\n DATAPACKING=POINT\n DT=(SINGLE SINGLE SINGLE)\n",N,(int)floor(N/2)+1);
    for(j=0;j<(int)floor(N/2)+1;j++) {
        for(i=0;i<N;i++) {
            fprintf(filepointer," %.9E %.9E %.9E\n",i/(xend-xstart),j/(zend-zstart),sqrt(pow(out[i*((int)floor(N/2)+1)+j][0],2)+pow(out[i*((int)floor(N/2)+1)+j][1],2)));
        }
    }
    fftw_destroy_plan(p);
    free(in);
    fftw_free(out);
    fclose(filepointer);
    return(1);
}
```

I am using the above program to initialize a 2D field of 2000x2000 elements of data (say in a x-z plane) and am using FFTW to transform the data from the real plane to the wavenumber plane, which, as per my best knowledge, has 2000x1001 elements as per FFTW documentation. The data is given by the function:

$$
in=\cos(2\pi f_x+2\pi f_z)
$$

where $f_x$ is the frequency in x-direction, $f_z$ is the frequency in z-direction. I expect the maximum amplitude at my chosen wavenumbers to be 1.0. However, I get the number 1999959 as the maximum amplitude. I know that I have to normalize this to get the maximum amplitude, but how do I calculate the number that I have to divide the maximum amplitude with to arrive at Amplitude=1.0? If the number is 2000x1001, does this mean that I cannot obtain the amplitude exactly as 1.0 due to some numerical error?


问:
FFT convolution works only with certain domain length
in my quest to understand how I can use FFT to compute integrals (see my other question click, still no answer there), I came across the fact that a convolution of two functions can be calculated by multiplying the Fourier transforms of the functions and taking the inverse transform of that product.

Now suppose I want to calculate
$$\int_{0}^R x \text{sin}(y-x)dx. $$
I can write this as the convolution of $\text{sin}(x)$ with $x\cdot \text{Ind}_{x<R}$ with $\text{Ind}$ as an indicator function that gives the proper integral boundaries.

My $x$ vector is given on N data points, equally spaced between 0 and some L, and to these N points I also know $\text{sin}(x)$.

The strange thing is, if $L=2\pi$, the method seems to work. But with another L, say $L=1$ it does not work. What is the reason for this, and what can I do?

I have read about the method of zero-padding where I add N-1 zeros to both function vectors. Would that help? And if so, I want to have a vector of the original size N after the ifft (which will give me a vector of 2N-1), so that I have the convolution value to each of my N $x$ values. Is there a way to do this?

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, fftfreq, ifft

pi2 = 2 * np.pi

N = 2001   ## should be odd!!
R = 0.2
L = 1

x = np.linspace(0, L-L/N, N)
f = np.sin(x)
g = x.copy()
g[g>R]=0      ## indicator function

# f=np.pad(f, (0, N-1), 'constant')   ## zero padding??
# g=np.pad(g, (0, N-1), 'constant')

F = fft(f)
G = fft(g)
C = (L/N)*ifft(F*G)

analytic_sol = -np.sin(R-x)+R*np.cos(R-x)-np.sin(x)   ## analytic result

plt.plot(x, C, label="test")
plt.plot(x, analytic_sol, label="analytic",linestyle="--")
plt.legend()
plt.show()
```

Notice how the solution only becomes valid at x=0.2=R, i.e. when the indicator function becomes zero. What happens here?


问:
How to use RODFT00 and REDFT00
I have some difficulty in implementing RODFT00 and REDFT00. I want to use them for fluid simulations.

I would really appreciate some help with my following problems:

I have read the documentation
(http://www.fftw.org/fftw3_doc/Real-even_002fodd-DFTs-_0028cosine_002fsine-transforms_0029.html#Real-even_002fodd-DFTs-_0028cosine_002fsine-transforms_0029)
and have been unable to understand

The r2r kinds for the various REDFT and RODFT types supported by FFTW, along with the boundary conditions at both ends of the input array (n real numbers in[j=0..n-1]), are:

FFTW_REDFT00 (DCT-I): even around j=0 and even around j=n-1.
FFTW_RODFT00 (DST-I): odd around j=-1 and odd around j=n.

Note that these symmetries apply to the “logical” array being transformed; there are no constraints on your physical input data.

What does it explains about the implementation of cosine and sine transforms?

Since I want to obtain the differentiation using the transforms. How should I arrange the wavenumbers to multiply with the matrix in both the cases??

Title: How to use RODFT00 and REDFT00


问:
How to move particles in a variational Monte Carlo simulation?
I'm attempting to implement some code that moves particles and then calculated the acceptance/rejection of the move, but I'm stuck in a rut. Here is my Python code:

```<code>def VMC(WF,numSteps):
     seed=5   ### my random number seed
     EnergyList=[]
     R=numpy.zeros((2,3),float)
     movesAttempted=0.0
     movesAccepted=0.0
     print "num-steps: ",numSteps
     for i in range(numSteps):

            OldPos= R.copy()
            oldWfc= WF.WaveFunction(R)
            for ptcl in xrange(0,len(R)):
             R[ptcl] = numpy.add( R[ptcl], (numpy.random.rand(3) - 0.5)*3.0 )
            newWfc= WF.WaveFunction(R)

            ratio = (newWfc**2/oldWfc**2)

            rander = numpy.random.rand(1)
            if ratio > rander:
             Eloc = WF.LocalEnergy(R)
             EnergyList.append(Eloc)
             movesAttempted+=1
             movesAccepted+=1
            else:
             movesAttempted+=1
             R = OldPos
     print "movesAcepted: ",movesAccepted
         print "  movesAttempted: ", movesAttempted
     print "Acceptance Ratio", movesAccepted/movesAttempted
     return EnergyList
</code>```

Here is a pastebin of more of the code: http://pastebin.com/7nTLDesy
I'm  getting a ridiculous value of the acceptance ratio, which makes me think I'm not getting the concept of how the electrons should be moved at all.

Specifically: When I run this code I get an acceptance ratio of 0.478, while for VMC calculations, moving a particle in a box of length 1.5 (in each direction) around its current position (moving in 3d) the acceptance ratio is supposed to be   0.323.

My problem is mainly to know how atoms are supposed to be moved in a VMC calculation and whether how I'm doing it is o.k.

I've seen an implementation in Fortran (metropolis), and I was wondering how different mine is and how I can correct it to give better results.

OwnerDisplayName: mike

ContentLicense: CC BY-SA 3.0


问:
Fitting a multivariate PDE (using Java)
I'm doing simulations of 2 coupled PDE's with Comsol Multiphysics. I want to fit some data (using the Application method, whose language is Java) to those simulations.

In order to answer my question you don't need to know what is Comsol Multiphysics or what is an Application Builder don't worry :)

So basically the equations look like : $\phi_t=\nabla.(\lambda\nabla\mu(\phi))+\Gamma(\phi,n)\tau_1$ and $n_t=D\nabla^2+G(\phi,n)\tau_2$.

And I have 3 parameters to fit :$\lambda,\tau_1,\tau_2$.

Here is a curve of the data with a curve of a simulation :

So I heard about many methods to fit univariate functions like the Bisection method, or the Brent-Decker, but I was wondering about how to apply them to a multivariate function, and more than this to my case with the simulation.

I was also wondering if there was a tool in Java to do it.

If I omitted some things are not clear in what I wrote please ask :)

Thanks in advance

Title: Fitting a multivariate PDE (using Java)


问:
Numerical integration with singularity term
In https://www.johndcook.com/blog/2012/02/21/care-and-treatment-of-singularities, the author explains the subtraction method to get rid of singularities when performing numerical integration.

The example he gives is the integral $\int_{0}^{1}\frac{1}{\sqrt{\sin x}}dx$.

To accurately compute this integral numerically, one first finds the approximation of integrand near $x=0$, which is $\frac{1}{\sqrt{x}}$ for this case.

Then, the original integral can be written like $\int_{0}^{1}(\frac{1}{\sqrt{\sin x}}-\frac{1}{\sqrt{x}})dx+\int_{0}^{1}\frac{1}{\sqrt{x}}dx$.

The latter term can be computed analytically by hand, so the author mentions that the work only left is computing the former one with numerical integration method such as trapezoidal, Simpson, etc.

I understand that the singularity is eliminated from the former integrand but I could not understand how one can actually compute it with numerical method.

For example, let's say I want to compute it with Simpson's method:

$\int_{0}^{2h}f(x)dx=\frac{h}{3}[f_{0}+4f_{1}+f_{2}]+\mathcal{O}(h^{5})$.

Then, I must evaluate $f(0)$ with the computer to fully compute the former integral. However, how can one compute $\frac{1}{\sqrt{sin(0)}}-\frac{1}{\sqrt{0}}$ with the computer? Shouldn't it give DivisionByZero error?

I want to understand how can this problem solved or am I misunderstand something?


问:
Using MATLAB to simulate the Ising Model
I am using MATLAB to simulate a 1D Ising Chain. I am running into an issue where when trying to find heat capacity, my system has a tremendous amount of noise. I'll post my code and an image of the heat capacity (as well as it smoothed 1000 times).

I have tried to smooth the energies first, I've tried having the "equilibrium energy" be the average of the last 5, 10, 100, 1000 timesteps (runs in the code), and nothing has worked. The max should be around 450, but smoothing will change the values.

My heat capacity definition is
$$C = \dfrac{\partial E}{\partial T}$$
But I have also tried using this definition to no avail:
$$C = kT \sigma_E^2$$

Any advice is welcome.

Edit 1

I changed the code so that the energy is smoothed before it is differentiated, and I smoothed it with different inputs. It somewhat resolved the problem.

%% Ising1D.m is a script designed to simluate a 1d Ising Chain at different Hfields and Temps
% The Metropolis Algorithm will be used
function [k, T, J, N, E, C, M, X] = Ising1D
%% Initializing the chain and other constants
N = 1000; % Number of spins
J = 1;
k = 1;
H = linspace(0,10,100); % Magnetic Energy
T = linspace(.1,10,100); % Temperature
flip = -1; % If I want to flip a spin, I'll multiply it by this.

%% Applying Metropolis Algorithm
time = 30*N; % Number of times I run the algorithm / mcs / Monte Carlo Steps
E = zeros(length(T), length(H), time); % A vector consisting of energy of system at different points in algorithm.
M = zeros(length(T), length(H)); % Magnetization
randSpin = randi(N,time,1); % An array of randomly chosen spins
randNum = rand(time,1); % An array of random numbers between 0 and 1.
for Tindex = 1:length(T)
    for Hindex = 1:length(H)
        temporary = rand(1,N); % Building my chain of spins randomly
        ups = temporary >= .5; % Up spins
        downs = -(temporary < .5); % Down spins
        chain = ups + downs; % Random initial set of spins.
        site = 1:(N-1);
        E(Tindex, Hindex, 1) = -J.*sum(chain(site).*chain(site+1)) - Hindex.*(sum(chain)); % Initialize energy
        for run = 2:time
            if randSpin(run) == 1 % calculate change in energy
                dE = 2*J*chain(randSpin(run))*chain(randSpin(run)+1) + 2*Hindex*chain(randSpin(run));
            elseif randSpin(run) == N
                dE = 2*J*chain(randSpin(run)-1)*chain(randSpin(run)) + 2*Hindex*chain(randSpin(run));
            else
                dE = 2*J*chain(randSpin(run))*(chain(randSpin(run)-1) + chain(randSpin(run)+1)) +  + 2*Hindex*chain(randSpin(run));
            end
            if ( dE < 0 )
                E(Tindex, Hindex, run) = E(Tindex, Hindex, run-1) + dE; % Update Energy.
                chain(randSpin(run)) = flip*chain(randSpin(run)); % Flip the spin for good.
            else
                p = exp(-dE/(k*T(Tindex)));
                if (randNum(run) <= p) % Accept the change
                    E(Tindex, Hindex, run) = E(Tindex, Hindex, run-1) + dE;
                    chain(randSpin(run)) = flip*chain(randSpin(run));
                else % Reject the change
                    E(Tindex, Hindex, run) = E(Tindex, Hindex, run-1);
                end
            end
        end
        M(Tindex, Hindex) = sum(chain); % M = N*avg(mu) = N*(ups+downs)/n = sum(chain)
    end
end

%% Smooth the data
E(:, 1, end) = smooth(E(:, 1, end), 0.15, 'rloess');
M(:,1) = smooth(M(:,1), .5, 'rloess');

%% Calculate and plot C
C = diff(E(:, 1, end))./diff(T)';
C = smooth( C/(N*k), .75, 'rloess'); % Scaled heat capacity
figure('Color', 'w');
subplot(1,2,1);
plot(linspace(0,8,length(C)), C); % The random energies at high temps make this very noisy. Smooth reduces noise via moving average.
title('Heat Capacity [C] vs Temperature [T]');
xlabel('T');
ylabel('C');
legend('H = 0');

%% Calculate X
X = diff(M,1,2)./diff(H);
subplot(1,2,2);
hold on;
plot(H(1:end-1), smooth(X(1,:), 1000)); plot(H(1:end-1), smooth(X(20,:), 1000)); plot(H(1:end-1), smooth(X(40,:), 1000)); plot(H(1:end-1), smooth(X(60,:), 1000)); plot(H(1:end-1), smooth(X(80,:), 1000)); plot(H(1:end-1), smooth(X(100,:), 1000));% Smoothing reduces noise but alters the data, namely by making the values smaller.
title('\chi vs Magnetic Field Strength [H]');
xlabel('H');
ylabel('\chi');
legend('T = 0.1', 'T = 2', 'T = 4', 'T = 6', 'T = 8', 'T =10');

end
```

Energy vs Temp

Not Smoothed Heat Capacity

Smoothed Heat Capacity

Expected behavior is like that of $x\times$Gaussian.


问:
Question regarding 1D implementation of the DG method
I'm pretty new to the DG method and have been writing a 1D code to help me understand the coding aspect.  With respect a reference, I've been following these notes https://www3.nd.edu/~zxu2/acms60790S15/DG-general-approach.pdf.  Where my confusion is comes from how to define the computational cells.  In the notes the author defines

$$
x_{i} = 0.5(x_{i-1/2} + x_{i+1/2})\\
I_{i} =  [x_{i-1/2},x_{i+1/2}].
$$

From a practical point of view, I am defining my mesh as follows

```<code>do i=1,ne
   x_r(i) = a + (i+0.5)*dx
   x_l(i) = a + (i-0.5)*dx
   x_c(i) = 0.5*(x_r(i) + x_l(i))
end do
</code>```

As a note, a refers to the left side of the domain, i.e. [a,b]. $x_{i+1/2}$ is x_r and $x_{i-1/2}$ is x_l. So my question is whether or not that is how I am supposed to define the left and right grid points.  I know that in structured finite volume methods, thats how you calculate the grid points at the half points (faces).  Thanks!


问:
What kind of problem or matrices are suitable for multigrid method?
For Poisson or Convection-diffusion equation as follows:
$$
-\Delta u=f,\qquad u|_\Omega = g.
$$
or
$$
-\Delta u +\vec{w}.\nabla u=f,\qquad u|_\Omega = g.
$$
using FDM or FEM discretization, we can get
$$
Ax=b.
$$
I have known that Multigrid method is a very efficient method for these equations and I have done some numerical examples using AMG code bulit-in some known software, e.g., ifiss. Indeed, I do not know why AMG works so efficient (because the principles seem so complex that I cannot understand, then I just use it), from which I indeed can get an optimal reslut, i.e., iteration step independent on mesh size $h$, and CPU time is almost linear with the unknowns of the system of linear equations.

My questions is as follows:

why amg works so efficient for Poisson and C-D equations? (can you give me some simple and clear explations without detailed principles?)

What other else matrices are suitable for amg method? so far, I have just known that matrices from Poisson and C-D equations are best fit for AMG? anyone else?


问:
What is the fastest algorithm for computing log determinant of a PSD matrix? (All possible PSD matrices)
I am diving into some literature to understand which is the best algorithm for computing the log-determinant of a PSD matrix. More generally, I am interested in a list of resources to read, which cover all the possible cases of matrix (i.e. sparse, dense, low rank, etc..) with all the probability of failure and error (relative, absolute, etc..).

So far I have found the following two papers:

Large-scale Log-determinant Computationthrough Stochastic Chebyshev Expansions

Fast Estimation of $tr(f(A))$ via Stochastic Lanczos Quadrature

A Randomized Algorithm for Approximating the Log Determinant of a Symmetric Positive Definite Matrix

We were just wondering if there are other algorithms with better asymptotics, different techniques, or other paper that we should be aware of.

Similar question on Computational Science are:

Calculating the log-determinant of a large sparse matrix


问:
Is a complete bacteria simulation with an exascale supercomputer possible?
Will it be possible to simulate a complete (at least simple) bacteria atom by atom on an exascale supercomputer? or is it possible already today with the largest systems?

Here, I've read that scientists can simulate an organelle of bacteria.
For that, they modeled a system consisting of 100M atoms. How to know which computational capacity was required for that?

Paper: "Atoms to Phenotypes: Molecular Design Principles of Cellular Energy Metabolism"

On Quora, there's an estimation that an E.coli organism would have about 9x1010 atoms (its genome size is about 4 million base pairs).

The smallest known bacterium so far has got about 112K nucleotides in its genome.

In 2014, there has been an atom by atom simulation of the polio virus capsid on the K computer with 10 petaflops. To include enough water molecules to fill and surround the capsid, Okazaki and colleagues needed to model the dynamics of nearly 6.5 million atoms, which they did for a simulated 200 ns.

A 2019 paper claims to simulate a non-biological environment with a total of 20 trillion atoms ... at up to 1.33 PFLOPS, trillion is here 1012. The used tool in this case models fixed molecules, compare to more advanced and complex ab initio simulation.


问:
Using nondimensionalization to solve an ode in MATLAB
I am trying to solve an ode that uses some extremely large numbers and some extremely small numbers, namely

$$
e = 1.6\times 10^{-19}\\
E = 10^6\\
\tau = 6\times 10^{-24}\\
m = 9.1\times 10^{-31}\\
c = 3\times10^8\\
$$

I introduced some dimensionless constants, one in particular being
$$
\varepsilon = \frac{\tau e}{cm} E \approx 4\times10^{-15}
$$

MATLAB's ode45, ode15s, and ode23s cannot handle these numbers. Am I allowed to just set $\varepsilon = 10$? That number produces a pleasant graph with  a nice solution, but I am not sure how to go from an $\varepsilon = 10$ world back to the original world.


问:
Truncated power series algebra implementation
1) I am looking for references for an efficient implementation and usage of TPSA. What sources exist besides Berz's 1989 original paper and the incomplete chapter in Dragt's book?

2) Are there reliable TPSA packages that can be used in Julia or Python?


问:
Can't plot correctly precession of perihelion of Mercury in MATLAB using ode45 or ode23
I was trying to plot precession of perihelion of Mercury using matlab. For this I am following a book Computational Physics by Nicholas J. Giordano and Hisao Nakanishi 2nd Edition. In that book authors used following force law predicted by general relativity:

$$
F \approx \frac{G M_m M_s}{r^2}\left(1+\frac{\alpha}{r^2}\right)
$$

Where $\alpha$, they state to be expressed by speed of light, mass of sun, eccentricity of the orbit.

To plot the orbit I have used that equation and ode45, ode23 functions in matlab. My code is as follows:

% Initial Conditions
% Initial position in AU
xm0 = [0.47 0];
% Initial velocity in AU/Year
vm0 = [0 8.2];

inivec = [xm0;vm0]

% tspan for 1000 years
tspan = [0, 1000];

[t, y] = ode45(@sunmer, tspan, inivec);

% x positions of Mercury
xmer = y(:,[1:2]);
% y positions of Mercury
vmer = y(:,[3:4]);
plot(xmer(:,1),xmer(:,2))
hold on

% Position of Sun in plot
plot(0,0,'o')
hold off
axis([-0.5, 0.5, -0.5, 0.5])
```

<p><strong>function for ode45 or ode23</strong></p>

<pre class="lang-matlab prettyprint-override"><code>function xpr = sunmer(t,x)
alpha = 0.0008;
xmer = x([1:2]);
vmer = x([3:4]);
r = norm(xmer);
rcubed = r^3;
rsqr = r^2;

% amer = -((G*Ms*x)/r^3)*(1+(alpha/r^2))
% G*Ms = 4*pi^2 in AU^3/Yr^2
amer = -(((4*(pi)^2))/rcubed)*xmer*(1+(alpha/rsqr));
xpr = [vmer;amer];
end
</code>```

I am getting following plot when using ode45:

If I use ode23 I get following figure:

These are not the plots I was expecting. I can't understand what is wrong in my code. Is it my code or the ODE solvers?


问:
Fast algorithm for computing lower mode shapes and natural frequencies in MATLAB using sparse stiffness and mass matrices
I am looking for a fast algorithm for computing eigenvalues and eigenvectors from sparse stiffness and mass matrices in MATLAB. The eig(K, M) doesn't work with sparse matrices.


问:
Methods to approximate obective function gradients from point cloud
Problem statement:
Assume that I have an objective function $f(x)$ which takes as input a  $D$-dimensional vector $x\in\mathbb{R}^D$, and that $f(x)$ is sufficiently smooth. Assume further that I have a set of $N$ samples $\{x_1,x_2,...,x_N\}$ for which I know the output of the objective functions $\{f(x_1),f(x_2),...,f(x_N)\}$. It is my task to minimize the objective function $f(x)$ by moving the samples towards its minimum.

Complication:
Now assume that I cannot evaluate the partial derivatives of the objective function $\nabla f(x)$ analytically, and that evaluations of $f(x)$ are very expensive and $D$ is very large, so that obtaining $\nabla f(x)$ through local finite differences is not an option. Assume further that the points are distributed arbitrarily.

Question:
I am looking for methods which would permit me to approximate the partial derivatives of the objective function $\{\nabla_{x_1}f(x_1),\nabla_{x_2}f(x_2),...,\nabla_{x_N}f(x_N)\}$ at each sample point $x_i, i=1,...,N$, using only the evaluations of the objective function at the sample points. Do you know of any methods which could be used towards this end?


问:
Tight binding model calculation with Extended Huckel Approximation
I've been reading Albright's Orbital Interactions in Chemistry. In the chapter on solids, he provided a general approach to find the band structure of a solid state system

Now if we are to model a one dimensional hydrogen chain, let $|1s(r)\rangle$ denote the $1s$ orbital centered at $r$. The band energy at each $k$ point will just be $e(k)=\frac{H_{11}(k)}{S_{11}(k)}$ according to the secular determinant. I then use a first neighbor interaction simplification, so $H_{11}(k)=\langle 1s(r)|H_{eff}|1s(r)\rangle+2\cdot\exp(ika)\langle 1s(r)|H_{eff}|1s(r-a)\rangle$ and $S_{11}=1+2\cdot\exp(ika)\langle1s(r)|1s(r-a)\rangle$.

It's worth noticing that assuming orthogonality of neighboring orbitals will set $S_{11}=1$, then $\langle 1s(r)|H_{eff}|1s(r-a) \rangle$ and $\langle 1s(r)|H_{eff}|1s(r)\rangle$ are just analogous to $\beta$ and $\alpha$ parameters in the simple Huckel model. But for the sake of realism I’m considering overlaps(the textbooks is also doing that) That is why the plots below show an asymmetric split with regard to hydrogen 1s ground state energy.

So far so good. I can even get a quick numerical estimate of $H_{11}$ using the extended Huckel approximation because then I just need to know the overlap integral between neighboring $1s$ orbitals. However, when I plot the result in Matlab it's different from the plot in the textbook:

Especially, the energy at the edge of the Brillouin zone is too positive. Anybody knows why that might happen? Thanks!

PS: I have to cross post this from Chemistry stackexchange because I'm feeling like this is more of a computational problem I have...


问:
Necessary information a topological optimisation solver needs to collect from a pre-processed CAD model
I am developing a solver that gets a CAD model as entry and does the topological optimisation calculation on it. My solver is inspired by the open-source codes presented in the literature.

Since it is just a solver I need to do the pre-processing and post-processing using a third-party code. I create my geometry and mesh using Abaqus and my solver retrieves the information in the .inp file in order to launch the calculation.

Currently, the information that my solver retrieves are as follows:

The coordinate of the nodes

element types

Material type

B.C., etc.

I am wondering what indispensable information any finite element code needs to retrieve from the CAD model and what information is dependent on the type of calculation and can be ignored from one simulation to the other?


问:
How to avoid density getting "deleted" in two way rigid body coupling with LBM CFD?
I've been reading this paper recently, which talks about using Lattice Boltzmann methods and two way coupling.  Specifically, it outlines fluid solid coupling, and solid fluid coupling, and how simply using both results in two way coupling.  The force from the fluid on the ridgid body makes sense, what doesn't make sense is how movements from the ridgid body can actually effect the fluid.

It only explains boundary conditions, it never explains what happens to fluid when you physically move the rigid body.

For example, what if during an update step, my box moves from a set of voxels A, with centroid at (0,0,0), to a set of voxels B, with centroid (10,0,0)?  When the body leaves an area, that makes sense, you are left with zero density voxels which will be filled in automatically by LBM. What happens on the opposite end makes no sense to me. When we move through fluid in real life that fluid should be displaced.  I do not see how this is accomplished in the paper.

To me it looks like they just... ignore and delete this density.  That does not make sense to me, and would be a huge source of error.  The closest thing I can find where they structurally consider even thinking about this is in two parts:

Section 5.3 pg 40 (simple boundary reflection):

This technique is limited in two ways: First, moving boundaries do not
  accelerate the fluid, and second, the accuracy is limited to grid
  spacing. This poses an additional issue, since it creates a
  “staircase-effect”, which shows non-realistic characteristics (Figure
  5.4).

and here:

Section 5.3 pg 45 (approach using fluid solid fraction, ε):

A disadvantage of this approach is that solid cells are still handled
  as fluid cells. This means that the collision and streaming steps are
  still applied, even when the collision result is multiplied by 0 and
  the streaming step just bounces the result between two neighboring
  nodes. Analytically this does not cause problems, but incoming fluid
  molecules that cannot escape the solid cause the density in border
  cells to rise. When the float parameter overflows, the GPU treats the
  value as infinity, which is not rectified by multiplying by 0.
  Additionally, this can cause problems when a solid moves away from a
  fluid cell holding high density values. By capping the density value,
  this can be avoided. Holdych [Holdych, 2003] proposes a solution to
  this problem for every cell whose ε > 0.95 that applies the assumption
  that the cell density difference between two neighboring cells is
  negligible. A fluid fraction $\bar{ε}$ is defined as $$ε = \sum_{i} (1 − ε_i)$$ (5.9)
  where $ε_i$ is the solid fraction ε of the neighboring cell in direction
  $e_i$ . The fluid density for these cells is then calculated depending on
  $\bar{ε}$:   $$ρ =
\begin{cases}
\frac{1}{\bar{ε}} \sum_{i}(1 − ε_i)ρ_i,  & \text{if $\bar{ε}$ > 0.01} \\
0, & \text{otherwise}
\end{cases}$$ where $ρ_i$ is the
  density value of the cell in direction $e_i$ .

And while I believe I understand what is going on here, there still doesn't appear to be any solution to the "moving boundary doesn't move fluid" problem, yet the whole thing comes off as if this has already been solved!

What am I missing here?


问:
Large-scale optimization of nonlinear equations
I'm looking to find a computationally efficient solution to a large system of nonlinear equations. I'm trying to maximize the following function:

$$ f(\vec{x}) = \sum_i^N C_i (x_i-A_i)x_i^{\epsilon_{ii}} \prod_{j\not = i}^N x_j^{\epsilon_{ij}}$$

Subject to the constraints:
$$ C_i > 0  \quad\forall\: i$$
$$ x_i > A_i > 0  \quad\forall\: i$$
$$ {\epsilon}_{ii} < -1  \quad\forall\: i$$
$$ {\epsilon}_{ij} \ge 0  \quad\forall\: i,j \quad\text{s.t.}\; i \not= j$$

I've been able to write a few solvers that can solve this, but I'm struggling to find/write a solver that can scale to $N=1000000$ (one million). Unfortunately, the $\epsilon_{ij}$ matrix is dense.

Things I've tried:

scipy.optimize BFGS algorithm: takes 8 hours when N = 2000

a fixed point solver: takes 3 hours for N = 10000, but seems to get stuck in local maxima

I would really appreciate advice on what to investigate next. I was gonna start going down the path of either (a) hand deriving the Jacobian and Hessian and taking another stab at using a scipy.optimize solver or (b) trying to coerce this to look more like a geometric programming problem and using gpkit. I'm pretty new to optimization so I could very well be missing something obvious - please don't hesitate to point out anything I've overlooked!

If it helps, I do have distributed computing resources available (spark, dask, kubernetes, etc).


问:
Lagrangian Duality

Possible Duplicate:

  Linear programming boundedness

Consider the following LP:

$\max$  $\sum_{i=1}^N b_i \pi_i$

s.t. $\;\;$   $\pi_i-\pi_j\leq 1 $ $\quad$   for each $(i,j) \in \tilde{A}$

$\sum_{i=1}^N a_s^i \pi_i \leq 1$ for each $s \in S$

This LP is always feasible, since zero is a feasible solution. Now, in order to find the optimum value, we relax the LP by multiplying second set of constraints by nonnegative Lagrangian multipliers $\lambda_s, s\in S$ and add them to the objective
function. Thereafter, we get the following relaxed dual LP:

$L(\lambda)$: $\max$ $\sum_{i=1}^N (b_i+\sum_{s\in S}\lambda_s a_s^i)\pi_i - \sum_{s\in S}\lambda_s$

s.t.  $\pi_i-\pi_j\leq 1$ $\quad$ for each $(i,j) \in \tilde{A}$

Let us consider the Lagrangian dual LP. This LP can be formulated
as follows:

LD: $\min$ $L(\lambda)$

s.t. $\lambda \geq 0$

Now, the question is as follows. Both LP and LD problems are feasible (e.g. zero vecor). Moreover, if the LD is bounded so is the LP. What cab be said if the LD is unbounded. Consider that in my case LP is always feasible.

Title: Lagrangian Duality


问:
Immersed boundary method in FEniCS?
I have looked at the FEniCS tutorials and documentation but I cannot find any mention to the possibility of implementing an immersed boundary method (IBM) for fluid dynamics.

In particular, I want to develop a solver for the incompressible Navier-Stokes equations of internal flows with elastic walls. I am thinking of using a cartesian mesh with an IBM to avoid remeshing.

Is this possible in FEniCS? If yes, is there any code example or paper to demonstrate the procedure? If no, how tedious would it be to dig into the source code and implement it?

Title: Immersed boundary method in FEniCS?


问:
How to reproduce the numerical examples in Prof. Saad's Book about Krylov subspace methods?
After reading Prof. Saad' Book, "Iterative methods for Sparse Linear Systems, 2nd version", I want to do the numerical examples about the Krylov subspace methods not only to reproduce the results in his book but also to guarantee me really understand the idea of Krylov subspace methods, but I have some trouble.

First, in Page 98, section 3.7, he uses 5 matrices to test all his methods, the first 3 matrices are named F2DA,F2DB,F3D are from a software "SPARSKIT" using Fortran and under UNIX system, which is impossible for me to generate the 3 matrices, (I just can write matlab code under windows). so, I want to ask that whether anyone has done the numerical examples with these 3 matrices? Where can I download these 3 matrices? Or how to generate these 3 matrices?

Second, the last 2 matrices are from matrix market, named ORS and FID . So, I can only get these 2 matrices and use these 2 matrices to do numerical examples. (I feel that if an author uses the matrix from matrix market, it is so convient for us to get easily. ;). So I use the matrix ORS to do the first 2 experiments about full orthogonalization method (FOM) and GMRES on Page 168 and 180 respectively as follows:

Page  168 as follows:
  Example 6.1. Table 6.1 shows the results of applying FOM(10) with no preconditioning
  to three of the test problems described in Section 3.7. The column labeled
  Iters shows the total actual number of matrix-vector multiplications (matvecs) required
  to converge. The stopping criterion used is that the 2-norm of the residual be reduced by a factor of 107 relative to the 2-norm of the initial residual. A maximum
  of 300 matvecs are allowed. Kflops is the total number of floating point operations
  performed, in thousands. Residual and Error represent the two-norm of the residual
  and error vectors, respectively. Note that the method did not succeed in solving the
  third problem.

Table 6.1 A test run of FOM with no preconditioning
\begin{array}{ccccc}
\hline
Matrix&Iters&Kflops&Residual&Error\\
F2DA&109&4442&0.36E-03&0.67E-04\\
F3D&66&11664&0.87E-03&0.35E-03\\
ORS&300&13558&0.26E+00&0.71E-04\\
\hline
\end{array}
Page 180, Table 6.2 A test run of GMRES with
 no preconditioning

Example 6.2. Table 6.2 shows the results of applying the GMRES algorithm with
  no preconditioning to three of the test problems described in Section 3.7. See Example
  6.1 for the meaning of the column headers in the table. In this test, the dimension
  of the Krylov subspace is m = 10. Observe that the problem ORS, which could not
  be solved by FOM(10), is now solved in 205 steps.

\begin{array}{c|c|c|c|c}
\hline
Matrix&Iters&Kflops&Residual&Error\\
F2DA&95&3841&0.32E-02&0.11E-03\\
F3D&67&11862&0.37E-03&0.28E-03\\
ORS&205&9221&0.33E+00&0.68E-04\\
\hline
\end{array}

My question is that when I write the FOM(10) method in matlab 2018b, I really cannot obtain the same resluts (Example 6.1) with that of his book (and I am sure my code is correct). In addition, I use the matlab built-in function gmres.m to test the result (example 6.2). I also cannot obtain the same results with his book. As we know, in FOM and GMRES method, one iteration step will require 1 matrix-vector, so the number of matrix-vector multiplication can be obtained by the iteration step. I want to know that how to obtain the results in his book i.e., What should I do to obtain the same results with his book? or indeed, I may not get the same results with his book if using matlab? or Indeed, I do not have to obtain the same results as long as my matlab code can converge? Thanks very much.
below is my FOM(restart) matlab code.

% main.m
clc;clear;close all;
%   load the sparse matrix 'ORS' from matrix market and convert it to matlab format
load orsirr_1.mtx
i = orsirr_1(2:end,1);
j = orsirr_1(2:end,2);
elem= orsirr_1(2:end,3);
A = sparse(i,j,elem);
%   initizlization
n = length(A);
x_star = ones(n,1);%    set the exact solution is [1,...,1]'
b = A*x_star;%  generate the right hand side of Ax=b
tol = 1e-7;
maxit = 30;%    maximum the outer iteration
restart = 10;%  restarted steps
x0 = zeros(n,1);%   initial guess vector

%%  FOM(restart)
fprintf('\n************  FOM(restart) method :  ****************\n')
[x,flag,relres,iter,resvec] = myFOM_restart(A,b,tol,maxit,x0,restart);
Iters=iter, Residual =resvec(end), Error = norm(x_star-x)

%%  gmres(restart)
fprintf('\n************ gmres(restart)method:   ****************\n')
[x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,[],[],x0);
Iters=(iter(1)-1)*restart+iter(2);
Iters, Residual =resvec(end), Error = norm(x_star-x)
```

<p>FOM(restart) function is </p>

<pre class="lang-matlab prettyprint-override"><code>function [x,flag,relres,iter,resvec] = myFOM_restart(A,b,tol,maxit,x0,restart)
%   Full Orthogonalization Method (FOM(m)) using modified Gram-Schmidt with
%   restarted       solving Ax=b
%   input
%           A       real nonsingular matrix in n X n
%           b       right hand side
%           tol     tolerance of relative residual norm: ||r_m||/||r_0||<tol
%           maxit   outer maxmimu iterations
%           restart restarted steps
%           x0      initial guess vector usually taken zero
%   output
%           x       approximate solution
%           flag    convegence :0
%                   unconvergence 1
%           relres  relative residual norm ||r_m||/||r_0||
%           iter    iteration steps
%           resvec  all the residual vector norm  [||r0||,...||rm||]
%   written  on 2019 12.14
%   reference  P167 Algorithm 6.4.1 from
%   Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia,
%                                           PA, second edition, 2003.

%%  initialize space to guarantee efficience
m = restart;
n = length(b);
H = zeros(m+1,m);
resvec = zeros(m+1,1);
V = zeros(n,m+1);
r0 = b-A*x0;
beta = norm(r0);
V(:,1) = r0/beta;
resvec(1) = beta;
flag = 1;

%%  begin to iteration
for k=1:maxit%  loop of maximum outer iteration steps
    %   modified Gram-Schmidt orthogonalization
    for j = 1:m
        w = A*V(:,j);
        for i=1:j
            H(i,j) = w'*V(:,i);
            w = w-H(i,j)*V(:,i);
        end
        H(j+1,j) = norm(w);
        if H(j+1,j) ==0
            flag = 0;
            fprintf('lucky breakdown!!!!!!!\n')
            break;
        end
        V(:,j+1) = w/H(j+1,j);

        %   compute the new residual vector
        e1 = zeros(j,1);e1(1)=1;
        ej = zeros(j,1);ej(j)=1;
        y = H(1:j,1:j)\(beta*e1);%  ym
        resvec(j+1,1) = norm(-H(j+1,j)*ej'*y*V(:,j+1));%    norm of residual r_j

        %   check convergence
        relres = resvec(j+1,1)/norm(r0);
        if relres < tol
            flag = 0;
            break;%
        end
    end%end of restart

    if flag==0
        break;
    end

    %   update new approximate solution
    x = x0+V(:,1:j)*y;
    x0 = x;
end %   end of outer maximum iteration

iter = j+(k-1)*restart;%    total iteration steps
if flag==0
    x = x0+V(:,1:j)*y;
    fprintf('** FOM(m) converges in %d steps to relative residual %5.2e*****\n',iter,relres);
else
%    fprintf('** FOM(m) unconverges in %d steps with relative residual %5.2e  *\n',iter,relres);
end
resvec(iter+2:end) = [];
end
</code>```

Below is my experiments results:

************  FOM(restart) method :  ****************
Iters =
   300
Residual =
   6.9544e+03
Error =
  187.4011

************ gmres(restart)method:   ****************
Iters =
   300
Residual =
  216.7682
Error =
   24.0811
```


问:
How to simulate a PDF using data samples like this?
I know two methods to simulate a PDF from random data samples using MATLAB :

1) Using a histogram where I use this command
histogram(data,'Normalization','pdf'), it gives PDF like bins.

2) Another one is using kernel density estimation

From (1) I get a PDF in bin format and from another (2), I get a continuous(smooth) PDF that I can match with my analytical derived PDF. I have seen in some scientific/engineering papers they plot simulated PDFs like the one given in the figure below:

where red is the simulated PDF and blue is analytical.
I want to know how one gets shaky simulated PDF like this and do they offer any benefits. May be represent data much clearly, have very small bins?


问:
Nonlinear least squares optimized Jacobian calculation
I have a nonlinear least squares problem, in which I am trying to minimize residuals which can be divided into four classes:
$$
\min_x ||\epsilon(x)||^2 + ||\xi(x)||^2 + ||\delta(x)||^2 + ||s(x)||^2
$$
The parameter vector $x$ can be naturally divided into four classes of parameters, $x = \left( x_1, x_2, x_3, x_4 \right)$. Some of the residuals are linear in the parameters, and others are not. The full Jacobian matrix is:
$$
J = \pmatrix{
\mathcal{E}_1 & \mathcal{E}_2 & 0 & \mathcal{E}_4(x) \\
\Xi_1 & \Xi_2 & \Xi_3 & 0 \\
\Delta_1 & \Delta_2 & 0 & 0 \\
S_1(x) & S_2(x) & S_3(x) & 0
}
$$
So you can see some blocks in the Jacobian don't depend on the parameters $x$ (these correspond to the residuals which are linear in $x$). So for example, $\mathcal{E}_i = \partial \epsilon / \partial x_i$.

Now since its a nonlinear least squares problem, I need to iterate, and solve the following linear least squares problem at each iteration:
$$
\min_p
\left| \left|
\pmatrix{ \epsilon(x) \\ \xi(x) \\ \delta(x) \\ s(x) } - J p
\right| \right|^2
$$
and then set $x := x + p$.

I want to take advantage of the sparse structure of $J$, and also the parts of $J$ which don't depend on $x$ to make the calculation as fast as possible. Another issue, is that $J$ is too large to fit in memory at once, so I am using the normal equations approach. $J$ has roughly 10 million rows and 20,000 columns.

With the normal equations approach, I can separate $J^T J$ into a part which can be precomputed prior to the iteration, and a part which needs to be computed on each iteration:
$$
J^T J =
\pmatrix{
\mathcal{E}_1^T \mathcal{E}_1 + \Xi_1^T \Xi_1 + \Delta_1^T \Delta_1 &&& \\
\mathcal{E}_2^T \mathcal{E}_1 + \Xi_2^T \Xi_1 + \Delta_2^T \Delta_1 & \mathcal{E}_2^T \mathcal{E}_2 + \Xi_2^T \Xi_2 + \Delta_2^T \Delta_2 && \\
\Xi_3^T \Xi_1 & \Xi_3^T \Xi_2 & \Xi_3^T \Xi_3 & \\
0 & 0 & 0 & 0
}
+ \pmatrix{
S_1(x)^T S_1(x) &&& \\
S_2(x)^T S_1(x) & S_2(x)^T S_2(x) && \\
S_3(x)^T S_1(x) & S_3(x)^T S_2(x) & S_3(x)^T S_3(x) & \\
\mathcal{E}_4(x)^T \mathcal{E}_1 & \mathcal{E}_4(x)^T \mathcal{E}_2 & 0 & \mathcal{E}_4(x)^T \mathcal{E}_4(x)
}
$$
So the first term doesn't depend on $x$, is small enough to fit in memory, and can be precomputed. The second term does depend on $x$ and must be computed on each iteration. I wrote only the lower triangles since these matrices are symmetric. Once I have the full $J^T J$, I use a Cholesky decomposition to compute the step $p$.

My question is: is the above approach the best way to implement something like this?

In the future, I would like to add more parameters to my model (i.e. $x_5, x_6, \dots$) which would increase the number of columns in my Jacobian. But then the $J^T J$ matrix will get more and more complicated, and it is already becoming difficult to compute all the different terms in the matrix (without making mistakes!). Does anyone have any suggestions for a strategy to compute all of these terms efficiently and simply? Is there a better method I should be using, such as TSQR?


问:
Stably solve transport equation with source term
I am trying to solve a transport equation of the form for the variable $\psi(t,r)$

\begin{equation}
\partial_t\psi-\alpha(r)\partial_r\psi-\beta(r)^2\psi-f(t,r)=0
,
\end{equation}

where I am solving for the function $f(t,r)$ from another PDE (a wave equation-so the amplitude of $f(t,r)$ oscillates in time and space). My problem is: I know from other arguments that $\psi$ should $not$ grow exponentially (essentially the term $\beta^2\psi-f$ should oscillate in time, so while sometimes $\psi$ will grow exponentially, it will also shrink exponentially at a later time so there is no long term secular growth in $\psi$). I am having trouble stably integrating this equation though, as if I get a small numerical error that makes $\psi$ too large, then $\beta^2\psi-f$ stops oscillating and becomes positive definite and I get runaway growth in $\psi$.

My question is: are there field redefinitions/solution techniques that could lead to more stable integration of this transport equation? Trying to, e.g. define $p\equiv \partial_t\psi-\alpha\partial_r\psi$ does not cure the basic problem.


问:
Multibody Systems modeling disadvantages
Multibody Systems modeling is a very systematic approach usually results in large sparse Jacbian matrix. I am working to model a system consisting of 11 bodies and 63 constraint equations as soon as i finished the assembly of the system`s Jacobian matrix I got a very ill conditioned matrix even I reviewed the model several times but still have this problem. Later I discovered that it is a common problem using such this method.
I am now in mid of my research there is no time to go step back so i tried to deal with this matrix and compared results with other methods there were great differences shown in graphs below . I need an urgent help how to overcome this problem as I said no time to go step backward.
I also tried to read about the saparsity there are dozens of research talking about that with a very complicated mathematics far away from my research field.
So Any suggestions please to solve this problem.


问:
How to compute the computational cost and storage of the Full Orthogonalization Method?
About the analysis of Full Orthogonalization Method (FOM) in Prof. Saad's book, wrote as follows:

Algorithm 6.4 (FOM):
\begin{array}{l}
 r_0=b-Ax_0,\beta=\|r_0\|_2,v_1 = r_0/\beta\\
 Define \quad H_m = \{h_{ij}\}_{i,j=1,\dots,m}, Set\quad H_m=0;\\
 For \quad {j=1,\dots,m}\\
 {\quad w_j = Av_j\\
  \quad\quad For\quad{i=1,\dots,j}\\
   \quad\quad{h_{ij}=(w_j,v_i)\\
   w_j=w_j-h_{ij}v_i\\
  }\\
\qquad end\\h_{j+1,j}=\|w_j\|_2\\
 v_{j+1}=w_j/h_{j+1,j}\\
 }\\
 end \\y_m = H^{-1}_m(\beta e_1),x_m=x_0+V_my_m.\\
\end{array}

This is the FOM method from Krylov subspace for solving $Ax=b$, where $A\in \mathbb{R}^{n\times n}$ is nonsingular.
Considering the computational cost and storage it is wrote as

A rough estimate of the cost of each step of the algorithm is determined as follows.
  If Nz(A) is the number of nonzero elements of A, then m steps of the Arnoldi
  procedure will require m matrix-vector products at the cost of 2m × Nz(A). Each
  of the Gram-Schmidt steps costs approximately 4 × j × n operations, which brings
  the total over the m steps to approximately $2m^2n$. Thus, on the average, a step of
  FOM costs approximately
  2Nz(A) + 2mn.
  Regarding storage, m vectors of length n are required to save the basis Vm. Additional
  vectors must be used to keep the current solution and right-hand side, and a
  scratch vector for the matrix-vector product. In addition, the Hessenberg matrix Hm
  must be saved. The total is therefore roughly
  $$(m + 3)n +\frac{m^2}{2}.$$
  In most situations m is small relative to n, so this cost is dominated by the first term.

I want to ask that How to estimate the computational cost and storage, in general?


问:
Nodal and Element Equilibrium in FEM Solution
I am learning FEM from KJ Bathe's textbook and it is mentioned that for a general FEM solution, nodal and element equilibrium is satisfied. The explanation takes steps that I don't understand.

At any node, the sum of the element nodal point forces is in equilibrium with the externally applied nodal loads (which include all effects due to body forces, surface tractions, initial stresses, concentrated loads, inertia and damping forces, and reactions).

Most of these terms, I don't think of as corresponding to nodes but being within the element.

How are reactions included as externally applied nodal loads? I was thinking of reactions as being the element nodal point forces due to stress.

He also mentions that the element nodal forces are equivalent to the element internal stresses. That is not totally clear to me. I am guessing most of this comes from "lumping" but I haven't come across an explanation of why we can do that yet.

Any explanation would be greatly appreciated!


问:
Numerically solving a partial differential equation in python with Runge Kutta 4
I'm supposed to solve the following partial differential equation in python using Runge-Kutta 4 method in time.

$$
\frac{\partial}{\partial t}v(y,t)=Lv(t,y)
$$

where $L$ is the following linear operator

$$
L=\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)^{-1}\left\{-i\alpha(1-y^2)\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)-2i\alpha+\frac{1}{R}\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)^2\right\}
$$

The spatial dimension is discretized along the $N = 101$ grid points from $y = -1$ to $y = 1$.

The boundary conditions are v[0] = 0 and v[N-1] = 0 and

$$
\left.\frac{\partial v}{\partial y}\right|_{y=-1}=\left.\frac{\partial v}{\partial y}\right|_{y=1}=0
$$

I have to discretize the second- and fourth-derivative operators using central finite difference and create a function for them. Using those functions, I can then create the $L$ operator.

Finally, using the initial condition $v(y) = 0.02(1+\cos(\pi y))$ time step $dt= 0.01$ and parameter values $R=500$ and $\alpha=0.3$, I have to find $v$ after 10s.

I have written the following code where the $L$ operator is put in a matrix form. The first few iterations seem to work giving somewhat reasonable results but then the values explode and I get an error. I'm not sure where the problem is, I might be misunderstanding how to implement the problem.

Another point I'm unsure about is whether the inverse of the $\left(\frac{\partial^2}{\partial y^2}-\alpha^2\right)$ operator can be understood as the inverse of the matrix representation of that operator or if that's not the case.

Can anybody give me some pointers as to how to solve this?

```<code>import numpy
import scipy
from scipy import linalg
from matplotlib import pyplot
%matplotlib inline

N=101
L=1.0

#creating the grid points
y=numpy.linspace(-L+0.04,L-0.04,N-4)
dy=(2*L)/(N-1)

v=numpy.empty(N-4)

# Operator of the second derivative acting on v with respect to y
def D2_v(N,dy):

    # Setup the diagonal of the operator.
    D = numpy.diag((-2.0/ (dy**2)) * numpy.ones(N-4))

    # Setup the upper diagonal of the operator.
    U = numpy.diag((1.0/ (dy**2)) * numpy.ones(N - 5), k=1)
    # Setup the lower diagonal of the operator.
    L = numpy.diag((1.0/ (dy**2)) * numpy.ones(N - 5), k=-1)
    # Assemble the operator.
    D2_v = D + U + L
    # Setup the Neumann condition
    D[-1, -1] = (-7.0/4.0) / (dy**2)
    D[0, 0] = (-7.0/4.0) / (dy**2)

    return D2_v

# Operator of the fourth derivative acting on v with respect to y
def D4_v(N,dy):

    # Setup the diagonal of the operator.
    D = numpy.diag((6.0/ (dy**4)) * numpy.ones(N-4))

    # Setup the upper diagonal of the operator.
    U = numpy.diag((-4.0/ (dy**4)) * numpy.ones(N - 5), k=1)
    U2 = numpy.diag((1.0/ (dy**4)) * numpy.ones(N - 6), k=2)

    # Setup the lower diagonal of the operator.
    L = numpy.diag((-4.0/ (dy**4)) * numpy.ones(N - 5), k=-1)
    L2 = numpy.diag((1.0/ (dy**4)) * numpy.ones(N - 6), k=-2)

    # Assemble the operator.
    D4 = D + U + L + U2 + L2

    # Setup the Neumann condition
    D4[-1, -1] = (5.0) / (dy**4)
    D4[0, 0] = (5.0) / (dy**4)
    D4[1, 0] = (-7.0/4.0) / (dy**4)
    D4[-2, -1] = (-7.0/4.0) / (dy**4)

    return D4

D2=D2_v(N,dy)
D4=D4_v(N,dy)

# Operator L
def L_v(N,y,dy,R,alpha):

    Ma = alpha*numpy.identity(N-4)
    My = numpy.diag(y)
    I = numpy.identity(N-4)
    Lv = numpy.linalg.inv(D2-Ma**2)*(-1j*Ma*(I-My**2)*(D2-Ma**2)+1j*
     (-2.0)*Ma+(1.0/R)*(D4-D2*(Ma**2)-(Ma**2)*D2+(Ma**4)))

    return Lv

#parameters
alpha=0.3
R=500.0
dt=0.01

L1=L_v(N,y,dy,R,alpha)

#initial value for v
for i in range(N-4):
    v[i]=0.02*(1.0+numpy.cos(numpy.pi*y[i]))

#Runge Kutta implementation for 1000 time steps
for m in range(1000):

    k1 = dt*numpy.matmul(L1,v)
    k2 = dt*numpy.matmul(L1,v + 0.5*k1)
    k3 = dt*numpy.matmul(L1,v + 0.5*k2)
    k4 = dt*numpy.matmul(L1,v + k3)

    v = v + (1.0/6.0)*(k1+2*k2+2*k3+k4)
</code>```


问:
Computation of Troullier-Martins pseudowavefunctions
The computation of Troullier-Martins pseudo-wavefunctions has been
described in [1].
The pseudo-wavefunction $R^{\textrm{PP}}_l$ is defined by
$$
R^{\textrm{PP}}_l(r) =
\left\{
\begin{array}{ll}
R^{\textrm{AE}}_l(r)  & r \geq r_{\textrm{cl}} \\
f_l(r)  & r \leq r_{\textrm{cl}}
\end{array}
\right.
$$
Here
$$
f_l(r) = r^l \exp(p(r))
$$

$$
p(r) = c_0 + c_2 r^2 + c_3 r^3 + c_4 r^4 .
$$
where $c_i$ are unknown coefficients ($c_1=0$).

The functions $R^{\textrm{AE}}_l$ are known and
the coefficients $c_i$ are defined by:
$$\int_0^{r_{\textrm{cl}}} \vert f_l(r) \vert^2 r^2 dr =
\int_0^{r_{\textrm{cl}}} \vert R^{\textrm{AE}}_l(r) \vert^2 r^2 dr
$$
$$
f_l(r_{\textrm{cl}}) = R^{\textrm{AE}}_l(r_{\textrm{cl}})
$$
$$
f'_l(r_{\textrm{cl}}) = (R^{\textrm{AE}})'_l(r_{\textrm{cl}})
$$
$$
f''_l(r_{\textrm{cl}}) = (R^{\textrm{AE}})''_l(r_{\textrm{cl}})
$$
How are these equations solved? What algorithm should be used?

References:

Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 1991, 43 (3), 1993–2006. DOI: 10.1103/PhysRevB.43.1993.


问:
Fitting a plane with the Prewitt gradient operator
Prewitt gradient operator

Show that the Prewitt gradient operator can be obtained by fitting the least-squares plane through the 3 × 3 neighborhood of the intensity function.

Hint: Fit a plane to the nine points (x + $\delta$x,y + $\delta$y,I[x + $\delta$x,y + $\delta$y]) where $\delta$x and $\delta$y range through −1, 0, +1. Then, having the planar model z = ax + by + c that best fits the intensity surface (i.e. minimizes the least-squares error in the z-direction), show that the two Prewitt masks actually compute a and b.

Solution: We denote the nine points by ($x_j$,$y_j$,$z_j$), $j$ = 1,...,9. The least-squares error for the planar model is:
$$
\epsilon(a,b,c) = \sum_{j= 1}^n(ax_j, + by_j + c + z_j)^2\tag{1}
$$
By differentation we get
\begin{align*}
\nabla\epsilon &= \begin{bmatrix}
\sum_{j= 1}^n(ax_j, + by_j + c + z_j)x_j\\ \sum_{j= 1}^n(ax_j, + by_j + c + z_j)y_j\\
\sum_{j= 1}^n(ax_j, + by_j + c + z_j)
\end{bmatrix}\\\tag{2}
\end{align*}
by setting the gradient to zero
$$
\implies=\begin{bmatrix}
\sum_{j = 1}^nx_j^2&\sum_{j = 1}^nx_jy_j&\sum_{j = 1}^nx_j\\
\sum_{j = 1}^nx_jy_j&\sum_{j = 1}^ny_j^2&\sum_{j = 1}^ny_j
\\
\sum_{j = 1}^nx_j&\sum_{j = 1}^ny_j&n
\end{bmatrix}
\begin{bmatrix}
a\\
b\\
c
\end{bmatrix}
= \begin{bmatrix}
\sum_{j = 1}^nx_jz_j\\
\sum_{j = 1}^ny_jz_j\\
\sum_{j = 1}^nz_j
\end{bmatrix}\tag{3}
$$
For our 3 $\times$ 3 neighborhood we may fix the coordinate origin at the center of the mask. In this case $x_j$ and $y_j$ range through −1, 0, +1 and the normal equations have the form
$$\begin{bmatrix}
\color{blue}6 & 0 &0 \\
0 & \color{blue}6 & \color{red}9\\
0 & 0 & \color{blue}9
\end{bmatrix}
\begin{bmatrix}
a\\
b\\
c
\end{bmatrix}
= \begin{bmatrix}
\sum_{j = 1}^nx_jz_j\\
\sum_{j = 1}^ny_jz_j\\
\sum_{j = 1}^nz_j
\end{bmatrix}\tag{4}
$$
Therefore the optimal plane is obtained by
$$
\begin{bmatrix}
a\\
b\\
c
\end{bmatrix}
= \begin{bmatrix}
\frac{1}{6}\sum_{j = 1}^nx_jz_j\\
\frac{1}{6}\sum_{j = 1}^ny_jz_j\\
\frac{1}{9}\sum_{j = 1}^nz_j
\end{bmatrix}\tag{5}
$$
where the first two components are obtained by using the Prewitt masks:
$$
M_x = \begin{bmatrix}
-1 & 0 & 1\\
-1 & 0 & 1\\
-1 & 0 & 1\\
\end{bmatrix}, \quad
M_y = \begin{bmatrix}
1 & 1 & 1\\
0 & 0 & 0\\
-1 & -1 & -1\\
\end{bmatrix}\tag{6}
$$

My questions

In (4), where do the $\color{blue}6$, $\color{blue}6$ and $\color{blue}9$

In (4), shouldn't the $\color{red}9$ be a $\color{red}0$ instead?

How can we get the insight between (5) and (6) (the relation with Prewitt mask)?


问:
Computing convolution of two characteristic function over a 1D Cartesian mesh
I am trying to compute the convolution of two characteristic functions over a  Cartesian mesh.  First, I define my Cartesian mesh of the interval $[0,1]$ as follows

$$
x_{i} = i \Delta x, i = 0, 1, 2\\
\Delta x = \frac{1}{N},
$$
with $N$ being the number of grid points.  Secondly, we can compute the convolution of two characteristic functions analytically, which gives

$$
f(x) = \mathcal{X}_{(-1/2,1/2)}(x) \star \mathcal{X}_{(-1/2,1/2)}(x) = \begin{cases}
      0 & x > 1/2\;\text{or}\;x < 1/2 \\
      x + 1 & x \in [-1/2,0] \\
      1-x & x \in [0,1/2]
   \end{cases}
$$

When computing $f(x)$ over a mesh, would I calculate the characteristic function like this?
$$
f(x_{i}) = \mathcal{X}_{(-1/2,1/2)}(x_{i}) \star \mathcal{X}_{(-1/2,1/2)}(x_{i}) = \begin{cases}
      0 & x_{i} > 1/2\;\text{or}\;x_{i} < 1/2 \\
      x_{i} + 1 & x_{i} \in [-1/2,0] \\
      1-x_{i} & x_{i} \in [0,1/2]
   \end{cases}
$$

As a side note, I know that the convolution of these two functions provide a B-Spline.


问:
Processing time steps in chunks with Fortran
My PDE simulation program written in Fortran has to make about 2 million variable time steps. But with each time step it slows down more and more, so that if it initially makes 1000 time steps per second, then when it reaches time step 600k and above, it makes 1000 time steps in as many as several minutes. And the slow down continues the further the simulation runs.

As far as I have established and understood, the problem is related to bufferization. There is enough memory and it does not get filled completely as the simulation progresses, but it does slow down quite a bit as the time steps progress. It looks like some caches get filled so that more and more latency occurs.

The solution would be to free up memory after every certain number of time steps, thus running the simulation in chunks. However, I don't yet know how to do that in Fortran. Does one have to use something like call Malloc() to dynamically allocate each arrays in Fortran, and then use call Free() to deallocate memory after a certain number of time steps has been done? Or is there a better approach toward this problem?
Would appreciate some advice.


问:
How to account for a corner node with zero-flux condition at an extrapolated distance
I am trying to implement a numerical solver and am having troubles dealing with boundary conditions, especially in the corners.

I have a 2D mesh, and on the left I have a Dirichlet condition, on the bottom I have a Neumann condition, and on the top and on the right I have a zero-flux vacuum boundary condition.

That is:

On the left: $u = B$, a constant value

On the bottom: $\frac{du}{dn} = 0$, symmetry of the problem

On the right and top: $\frac{1}{u} \frac{du}{dn} = \frac{-1}{d_n}$, where $d_n$ is an extrapolation length. So we're saying that the flux will be technically zero at that distance from the boundary

I can set up my finite difference equations for those mostly, with my Laplacian depending on $u_{i,k-1}$, $u_{i-1,k}$, $u_{i,k}$, $u_{i+1,k}$ and $u_{i,k+1}$, with my boundary conditions looking like:

Bottom left corner (i=0, k=0):

$u_{i,k} = B$

Bottom right corner (i=0, k=L):

$u_{1,k}$ = $u_{-1,k}$ [ghost point]

$\frac{u_{i,k+1} - u_{i,k-1}}{\Delta z + d} = -\frac{u_{i,k}}{d}$   $\rightarrow$ $u_{i,k} = \frac{d}{d+\Delta z} u_{i,k-1}$

And I replace everything in my Laplacian so that it only depends on $u_{i+1,k}$ and $u_{i,k-1}$

Rest of bottom points (i=0, k between 0 and L):

Use Neumann to replace the ghost point $u_{i-1,k}$ with $u_{i+1,k}$

Top left corner (i=R, k=0)

$u_{i,k} = B$

Top points (i=R, k between 0 and L):

$\frac{u_{i-1,k} - u_{i+1,k}}{\Delta r + d} = -\frac{u_{i,k}}{d}$       $\rightarrow$ $u_{i,k} = \frac{d}{d+\Delta r} u_{i-1,k}$

Right side points (i between 0 and R, k=L):

$\frac{u_{i,k+1} - u_{i,k-1}}{\Delta z + d} = -\frac{u_{i,k}}{d}$   $\rightarrow$ $u_{i,k} = \frac{d}{d+\Delta z} u_{i,k-1}$

And I replace everything in my Laplacian so that it only depends on $u_{i-1,k}$, $u_{i+1,k}$ and $u_{i,k-1}$

And I replace everything in my Laplacian so that it depends on u(i,k-1), $u_{i,k+1}$ and $u_{i-1,k}$

Top right corner (i=R, k=L):

This one has me stuck. Basically, I have two "concurrent" conditions applying there.

$\frac{u_{i-1,k} - u_{i+1,k}}{\Delta r + d} = -\frac{u_{i,k}}{d}$       $\rightarrow$ $u_{i,k} = \frac{d}{d+\Delta r} u_{i-1,k}$

$\frac{u_{i,k+1} - u_{i,k-1}}{\Delta z + d} = -\frac{u_{i,k}}{d}$    $\rightarrow$ $u_{i,k} = \frac{d}{d+\Delta z} u_{i,k-1}$

And I do not know how to apply them.

Any ideas? Or does something look off somehow?

Hopefully this is the right Stack Exchange place, if not please direct me to the correct one!


问:
Radial Hartree and exchange-correlation potentials
According to [1] the one-dimensional Kohn-Sham equation is given by
$$
\left( -\frac{1}{2} \frac{d^2}{dr^2} + \frac{l(l+1)}{2r^2} + V[\rho;r]\right)
rR_{nl}(r) = \varepsilon_{nl} rR_{nl}(r)
$$
where
$$
V[\rho;r] = -\frac{Z}{r} + V_\textrm{H}[\rho;r]+V^\textrm{LDA}_\textrm{xc}(\rho(r))
$$
How are the radial Hartree potential $V_\textrm{H}[\rho;r]$ and radial exchange-correlation potential $V^\textrm{LDA}_\textrm{xc}(\rho(r))$ computed?

[1] Troullier and Martins. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B vol. 43, no. 3, pp. 1993-2006 (1991).


Title: Unforeseen behaviour of fix() in Octave and Scilab

Body:
I encountered some odd results when using the function fix() in Octave and Scilab. The following is input and output of the Scilab console but the exact same thing happens in Octave.

I start with the following matrix:

```<code>A=[1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
   0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0;
   0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0;
   0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0;
   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1;
   1,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0;
   0,1,0,0,1,0,0,0,1,0,0,0,1,0,0,0,0,0,0;
   0,0,1,0,0,1,0,0,0,1,0,0,0,1,0,0,1,0,0;
   0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,1,0;
   0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1;
   0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0;
   0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0;
   1,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,1;
   0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0]
</code>```

I want to take it to the reduced row echelon form:

```<code>-->rref(A)
 ans  =

   1.   0.   0.   0.   0.   0.   0.   0.   0.   1.   1.   0.   0.   1.   2.   1.   0.   1.          1.
   0.   1.   0.   0.   0.   0.   0.   0.   0.  -1.   0.   0.   0.  -1.  -1.   0.   0.   0.          0.
   0.   0.   1.   0.   0.   0.   0.   0.   0.   0.  -1.   0.   0.   0.  -1.  -1.   0.  -1.         -1.
   0.   0.   0.   1.   0.   0.   0.   0.   0.  -1.  -1.   0.   0.   0.  -1.   0.   0.  -1.388D-16   0.
   0.   0.   0.   0.   1.   0.   0.   0.   0.   0.  -1.   0.   0.  -1.  -1.  -1.   0.  -1.          0.
   0.   0.   0.   0.   0.   1.   0.   0.   0.   1.   1.   0.   0.   1.   1.   1.   0.   0.          0.
   0.   0.   0.   0.   0.   0.   1.   0.   0.   0.   1.   0.   0.   0.   1.   0.   0.   1.          0.
   0.   0.   0.   0.   0.   0.   0.   1.   0.   0.   0.   0.   0.  -1.  -1.  -1.   0.  -1.         -1.
   0.   0.   0.   0.   0.   0.   0.   0.   1.   1.   1.   0.   0.   1.   1.   0.   0.   1.          0.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   0.   0.   0.   1.   0.   0.          1.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   1.   1.   1.   0.   0.          0.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   1.          1.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.          0.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.          0.
</code>```

So far so good, I guess, but when I now use fix() on the last obtained matrix to get integer entries, I get zeros at (6,10) and (6,11), where I expected ones:

```<code>-->fix(rref(A))
 ans  =

   1.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   0.   0.   1.   2.   1.   0.   1.   1.
   0.   1.   0.   0.   0.   0.   0.   0.   0.  -1.   0.   0.   0.  -1.  -1.   0.   0.   0.   0.
   0.   0.   1.   0.   0.   0.   0.   0.   0.   0.  -1.   0.   0.   0.  -1.  -1.   0.  -1.  -1.
   0.   0.   0.   1.   0.   0.   0.   0.   0.  -1.  -1.   0.   0.   0.  -1.   0.   0.   0.   0.
   0.   0.   0.   0.   1.   0.   0.   0.   0.   0.  -1.   0.   0.  -1.  -1.  -1.   0.  -1.   0.
   0.   0.   0.   0.   0.   1.   0.   0.   0.   0.   0.   0.   0.   1.   1.   1.   0.   0.   0.
   0.   0.   0.   0.   0.   0.   1.   0.   0.   0.   1.   0.   0.   0.   1.   0.   0.   1.   0.
   0.   0.   0.   0.   0.   0.   0.   1.   0.   0.   0.   0.   0.  -1.  -1.  -1.   0.   0.  -1.
   0.   0.   0.   0.   0.   0.   0.   0.   1.   0.   0.   0.   0.   1.   1.   0.   0.   0.   0.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   0.   0.   0.   1.   0.   0.   1.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   1.   1.   1.   0.   0.   0.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   1.   1.   1.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
</code>```

What happened there?


问:
Effective fitness formula in the Moran process on a game
I was recently reading some literature about evolutionary game theory and I got particularly interested in Moran process linked to prisoner's dilemma as a model of evolution of two sub-populations.

Please take a look at this paper:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5701760/

I understand the model rather well yet I am missing the point of eq.(4) and (5). Why does the fitness of a single cell $f$ takes this particular form of the payoff $\pi$:

$f = 1 - w + w\pi$

I understand that $w$ is supposed to be a 'selection pressure' but why do we have the exact equation as above? why not just $w\pi$ or $1-w\pi$? This is troubling me especially it is the same equation for both strategies i.e the selection would act in the same way on both subpopulations...


问:
Comparison of diffusion time - theoretical value vs computed
This is a follow up to my previous post

I've been trying to compare the diffusion time obtained from theoretical
derivation(answered in my previous post) and what is obtained computationally, for a 1D diffusion equation.

Theoretically, the diffusion time is given by

$$t_D = \frac{l^2}{2D}$$
Calculating $t_D$ , l = 5 nm and D = 500 $nm^2/min$ gives  $t_D$ = 0.025 min.

The following is the MATLAB code that simulates 1D diffusion system using pdepe solver.

```<code>function sol=so()
format short
global D nnode init_co find_index
m = 0;
xend = 5;
D = 500;
x = 0:1:xend;
find_index  = 0:1:xend;
t = 0:0.00001:0.5;
init_co = 1*ones(length(x),1);
nnode = length(x);
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t)
plot(t,sol)
function [g,f,s] = pdefun(x,t,c,DcDx)
g = 1;
f = D*DcDx;
s = 0;
end

function c0 = icfun(x)
c0 = init_co(find(find_index==x));
end

function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
% Dirichlet at left (concentration at left boundary = 2, Neumann at right(dC/dx = 0)
    pl = cl - 3;
    ql = 0;
    pr = 0;
    qr = 1;
end
end
</code>```

Result:

From the plot of C vs time, all the curves reach a steady-state at 0.1 min.

Could someone explain why there is a difference in what is computed theoretically and computationally? Is it appropriate to make this comparison?


问:
What is the reason for this finite-difference high errors on non-uniform grid?

tl;dr

Using a Taylor-matched method to find coefficients for the discretized equation $ \mathbf{A} \vec{f}'' = \mathbf{B} \vec{f} $, a Fortran code has been implemented to find the second derivative and tested for $f=\sin$, with bad results on non-uniform grids. Which is the cause for that?

Methodology

I am new to non-uniform grids.
Using the concepts described in the following papers:

Lele, Compact finite difference schemes with spectral-like resolution, 1992

Gamet, Compact finite difference schemes on non-uniform meshes. Application to direct numerical simulations of compressible flows, 1999

I started implementing a Fortran module to solve differential equations on non-uniform grids. The idea behind these method is to discretize the (second-order, in my case) derivative with the following matricial equation: $ \mathbf{A} \vec{f}'' = \mathbf{B} \vec{f} $.
The matrices are defined as multi-diagonal, where the number of secondary diagonals reflects the formal error of the method. For example, in Gamet the general form $\alpha_i f''_{i-1} + f''_i + \beta_i f''_{i+1} = A_i f_{i+1} + B_i f_{i-1} + C_i f _{i+2} + D_i f_{i-2} + E_i f_i $ is used and the matrix $\mathbf{A}$ will be tridiagonal, while $\mathbf{B}$ will be pentadiagonal (in Lele, a similar form is defined, but with different symbols and concerning just uniform grids).
$f_i$ and $f''_i$ can be expanded with a Taylor series on the grid points; as an example,

$f''_{i+n} = f_i'' + nhf_i^{(3)} + \frac{(nh)^2}{2!}f_i^{(4)} + \mathcal{O}(h^3f^{(6)})$

In order to determine the value of the coefficients, these expansions can be substituted in the above equation, then providing a system of equations matching terms of the same order (see below for how to solve this system).

The method is closed at the boundaries imposing the following forms:

$ i=1 \Rightarrow f''_1 + \alpha f''_2 = Af_1 + Bf_2 + Cf_3 + Df_4 $

$ i=2 \Rightarrow \alpha f_1'' + f_2'' + \beta f_3'' = A f_1 + Bf_2 + Cf_3 + Df_4 $

Strictly speaking, it follows that the matrices can have a higher number of diagonals at the boundary points.

In order to find the second derivative, the matrix $\mathbf{A}$ is inverted and $ \vec{f}'' = \mathbf{A}^{-1} \mathbf{B} \vec{f} $. To do so, the Lapack library is used, as shown in http://fortranwiki.org/fortran/show/Matrix+inversion.

Grid

Different grids have been considered. The Fortran code which defines them is here reported (corresponing lines should be uncommented):

```<code>do ig=1,ng
      !! uniform
      !grid(ig,1) = a + (b-a)*dble(ig)/ng
      !grid(ig,2) = (b-a)/ng
      !! Shukla
      !grid(ig,1) = (b+a)/2. + (b-a)*asin(-alpha*cos(PI_MATH*ig/ng))/(2.*asin(alpha))
      !grid(ig,2) = (b-a)*( asin(-alpha*cos(PI_MATH*ig/ng))/(2.*asin(alpha)) &
      !  - asin(-alpha*cos(PI_MATH*(ig-1)/ng))/(2.*asin(alpha)) )
      !! random
      !grid(ig,1) = a + (b-a)*(dble(ig)-1+rand())/ng
      !if(ig.eq.1) then
      !  grid(ig,2) = (b-a)/dble(ng)
      !else
      !  grid(ig,2) = grid(ig,1) - grid(ig-1,1)
      !endif
end do
</code>```

In particular, Shukla refers to the grid defined in the article by Shukla, Very high-order compact finite difference schemes on non-uniform grids for incompressible Navier–Stokes equations, 2007. The important thing to notice here is that this scheme has a parameter, alpha, which varies from 1 (uniform grid) to $\alpha\rightarrow 0$ (Chebyschev grid).

Testing

As a first test, the second derivative of a sine function has been evaluated and compared with the analytical second derivative:

$ x\in [0,1];\; f(x) = \sin(2\pi x); f''(x) = -(2\pi)^2\sin(2\pi x) $

Using a uniform grid, the results are fine for Lele

For Gamet, using a uniform grid and a Shukla grid with alpha = 1 the results are approximately as good as with Lele (the two grids are of course the same; I compared them in order to be sure about the results):

Here is a comparison of the errors

Problems arise when the grid is imposed as non-uniform: for alpha = 0.9

and for alpha = $10^{-10} \sim 0$

Comparison of errors for different values of alpha:

It appears that the problem is due to the boundaries, but it is not clear to me why it happens. I double-checked the code in search of typos (but having solved the system with Python they should be ruled out) and I observe the same problems with a different method (cfr. Shukla 2007). I would like to know if this behaviour is to be expected and what I could do about that.

Solving the system for the coefficients

The cited papers provide sets of coefficients for different orders of the formal error. However, some situations could require to develop a custom method, for example setting some coefficients to zero, so here is a possible way to solve the system, using Sympy, a library for Python. The second argument of sym.solve is the list of variables for which to solve the system.

If a lower order is desired, it is sufficient to remove one of them (which can be then fixed to an arbitrary value, or in order to match a given method, e.g. for adapting Gamet's method to Lele's one) and comment a corresponding number of equations in the system (starting of course from the bottom one).

If some constraint is desired, a corresponding equation can be substituted in the system (again, starting from the bottom one).

import sympy as sym
a,b,c,d,e,alpha,beta = sym.symbols("a b c d e alpha beta")
# hmm = h_{i-2}; h0 = h_i; hp = h_{i+1}; etc.
hmm,hm,h0,hp,hpp   = sym.symbols("hmm hm h0 hp hpp")
solution = sym.solve((
    a+b+c+d+e, # order -1
    hp*a-h0*b+(hpp+hp)*c-(h0+hm)*d, # order 0
    hp**2*a+h0**2*b+(hpp+hp)**2*c+(h0+hm)**2*d-factorial(2)/factorial(0)*(1+alpha+beta), # order 1
    hp**3*a-h0**3*b+(hpp+hp)**3*c-(h0+hm)**3*d-factorial(3)/factorial(1)*(hp*beta-h0*alpha), # order 2
    hp**4*a+h0**4*b+(hpp+hp)**4*c+(h0+hm)**4*d-factorial(4)/factorial(2)*(hp**2*beta+h0**2*alpha), # order 3
    hp**5*a-h0**5*b+(hpp+hp)**5*c-(h0*hm)**5*d-factorial(5)/factorial(3)*(hp**3*beta-h0**3*alpha), # order 4
    hp**6*a+h0**6*b+(hpp+hp)**6*c+(h0+hm)**6*d-factorial(6)/factorial(4)*(hp**4*beta+h0**4*alpha), # order 5
), (a,b,c,d,e,alpha,beta))
print("Gamet bulk:")
for sol in solution:
    print(sol)
    print(solution[sol])
    print()
```


问:
Heat diffusion simulation in a 3D piston using FENICS
I'm trying to simulate the heat diffusion in a 3D piston. I marked the boundaries on GMSH.

I have used a Dirichlet BC of 300 on the top face of piston. But the results look abnormal. There is a change in temperature along the edge on the top face, which should not be case.

Here is code snippet

```<code>from __future__ import print_function
from fenics import *
import numpy as np

mesh = Mesh("piston-3d.xml");

V = FunctionSpace(mesh, 'P', 1)

subdomains = MeshFunction("size_t", mesh, "piston-3d_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "piston-3d_facet_region.xml")

#Define neumann bcs
g_1= Constant(300.0)
g_2 = Constant(100.0)
g_3 = Constant(30.0)
g_4 = Constant(0.0) // On the cross-sectional faces, Neumann
g_5 = Constant(0.0) // On the cross-sectional faces, Neumann

bc1 = DirichletBC(V, g_1, boundaries, 1) // Top face of piston
bc2 = DirichletBC(V, g_2, boundaries, 2) // Liner face of the piston
bc3 = DirichletBC(V, g_3, boundaries, 3) // Inner side of piston and bottom

bcs = [bc1, bc2, bc3]

#define measures
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

u_D = Constant(100.0)

# Define initial value
u_n = interpolate(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = 0.0

a = dot(grad(u), grad(v))*dx

L = f*v*dx - (g_4*v*ds(4) + g_5*v*ds(5))

# Time-stepping
u = Function(V)

# Compute solution
solve(a == L, u, bcs)

vtkfile = File("solution_steady.pvd")
vtkfile << u
</code>```

My question is why is the temperature of Top face not remaining constant (300)? Is it because of the boundary sharing?


问:
Numerical methods. MDF (ILU) implementation
I am trying to implement Minimum Discarded Fill (MDF) Ordering algorithm for incomplete matrix factorization. The algorithm description is here on page 60 Preconditioning Techniques for a Newton–Krylov Algorithm
for the Compressible Navier–Stokes Equations.

There is a thing that I don't understand. For initialization I should use Algorithm 7 (page 62)

Compute the discard value for all nodes vj in the graph of A0 using Algorithm 7

But this algorithm uses the value that is calculated later. We have only variables with k-index

Compute the discard value, discard(vk) = ||Fk||F using the following iterations

And there is the k+1-index variable

lev(k+1)

I can’t understand whether it is a recursion or something regarding iterative methods... if it will be easy for anyone to understand this MDF algorithm, please tell me how to implement this moment.


问:
Levinson Recursion for Non Square Toeplitz Matrices
Given a rectangular Toeplitz Matrix $ H $, how could one solve:

$$ y = H x $$

For instance, $ H $ can be Linear Convolution Matrix of the filter $ h $:

$$ H = \begin{bmatrix}
{h}_{1} & 0 & 0 & \ldots &  & 0 \\
{h}_{2} & {h}_{1} & 0 & \ldots &  & 0 \\
{h}_{m} & \ldots & {h}_{1} & 0 & \ldots & 0 \\
0 & {h}_{m} & \ldots & {h}_{1} & 0 & \ldots & 0 \\
\vdots &  & \ddots &  &  & \vdots \\
0 & \ldots & {h}_{m} & {h}_{m - 1} & \ldots & {h}_{1} \\
0 & \ldots & 0 & {h}_{m} & \ldots & {h}_{2} \\
\vdots &  & \ddots &  &  & \vdots \\
0 & \ldots &  & 0 & \ldots & {h}_{m} \\
\end{bmatrix} $$

Which is clearly not square.

In the general case it would be generated by:

numRows = 10; %<! Or any other number
numCols = 20; %<! Or any other number
vR = randn(numRows, 1);
vC = randn(numCols, 1);
mH = toeplitz(vC, vR);
```

Since the Matrix is not square, I'm after the least squares solution:

$$ \arg \min_{x} {\left\| H x - y \right\|}_{2}^{2} $$

I have implemented the Levinson Recursion for the square case yet I'd like to know if there is an extension to it.

At the moment the trick I use for this specific case it to partition the matrix and data for a subset of square Toeplitz Matrices and use the algorithm I implemented.

I wonder if there is a better way.


问:
Integer partition algorithms
I am familiar with and have written MathCad algorithms for the partition functions (,),which gives the number of ways of partitioning  into   parts, (,), which gives the number of ways of partitioning  into   distinct parts, and (,,), which gives the number of ways of partitioning  into exactly  parts of size at most . I seek a recursive algorithm for (,,), which gives the number of ways of partitioning   into exactly  distinct parts of size at most . If there is an existing python algorithm, all the better. Thanks.

Title: Integer partition algorithms