问:
Generalized Hermite Function as eigenfunction of a differential operator
I'm going through this paper.

The article defines function function $\phi_n^\mu(x)$ that is orthonormal on $L^2$ with measure $dm = dx$:
\begin{equation}
 \phi^\mu_n =\left(\frac{\gamma_\mu(n)}{\Gamma(\mu+\frac{1}{2})}\right)^\frac{1}{2}\frac{1}{\sqrt{2^n}n!}e^{-\frac{x^2}{2}}H^\mu_n(x)x^\mu;  \mu,n\in{N},
\end{equation}

Also, it defines differential operator (one-dimensional case, $d=1$, $\mu_1=\mu$)

\begin{equation} \label{eq:hcal}
 \mathcal{H_\mu} = \frac{1}{2}\left(
   \sum^d_{k=1}\left[
    \mathcal{D}^2_{\mu_k}-\frac{2\mu_k}{x} \left(
     \frac{\partial}{\partial x}+\mu_k-1
    \right)
   \right] -|x|^2
 \right),
\end{equation}
where \begin{equation} \label{eq:dunkl}
 \mathcal{D}_\mu\phi(x)=\phi'(x) + \frac{\mu}{x}\left(
  \phi(x)-\phi(-x)
 \right),
\end{equation}

It is said (in the paper) that $\phi$ is the eigenfunction of $\mathcal{H}^\mu$ operator, i.e. $-2(\mathcal{H}^\mu \phi_n^\mu) (x) = 2(n+\mu+1/2)\phi_n^\mu(x)$ for one-dimensional case.

I have the functions written in Python: dunkl2 for $\mathcal{D}_\mu$, and calh for $\mathcal{H}_\mu$ operator:

```<code>def dunkl2(n,m,func,x):
return misc.derivative(func,x,0.001,2) + 2*m*(misc.derivative(func,x,0.001))/x - m*(func(x)-func(-x))/x**2

def calh(n,m,func,x):
return -dunkl2(n,m,lambda x: func(n,m,x),x) + 2*m/x*misc.derivative(lambda x: func(n,m,x),x,0.001) + (2*m/x*(m-1)+x**2)*func(n,m,x)

for i in range(0,10):
print("n =", i, "; m =", 10-i, "; x =", 1, "; (Cal.H * f)(x) / f(x) =", calh(i,10-i,phi,1)/phi(i,10-i,1)
</code>```

I test my Python functions by simply dividing both sides of the differential equation by $\phi_n^\mu(x)$ and checking if the result is $2(n+\mu+1/2)$, and this gives me the set of numbers below. But, as you can see, there's something wrong, at least because the resulting coefficient depends not only on $n$ and $\mu$, but on $x$ also.

Where is my mistake? I can assure you that this phi I use is indeed orthonormal.

Data:

```<code>n = 0 ; m = 9 ; x = 1 ; (Cal.H * f)(x) / f(x) = 108.999910168
n = 1 ; m = 8 ; x = 1 ; (Cal.H * f)(x) / f(x) = 74.9999101676
n = 2 ; m = 7 ; x = 1 ; (Cal.H * f)(x) / f(x) = 75.0000086295
n = 3 ; m = 6 ; x = 1 ; (Cal.H * f)(x) / f(x) = 49.0000086295
n = 4 ; m = 5 ; x = 1 ; (Cal.H * f)(x) / f(x) = 49.0000089668
n = 5 ; m = 4 ; x = 1 ; (Cal.H * f)(x) / f(x) = 31.0000089668
n = 6 ; m = 3 ; x = 1 ; (Cal.H * f)(x) / f(x) = 30.9999828753
n = 7 ; m = 2 ; x = 1 ; (Cal.H * f)(x) / f(x) = 20.9999828754
n = 8 ; m = 1 ; x = 1 ; (Cal.H * f)(x) / f(x) = 20.9999722442
n = 9 ; m = 0 ; x = 1 ; (Cal.H * f)(x) / f(x) = 18.9999722445

n = 0 ; m = 9 ; x = 2 ; (Cal.H * f)(x) / f(x) = 77.5000026668
n = 1 ; m = 8 ; x = 2 ; (Cal.H * f)(x) / f(x) = 61.0000026667
n = 2 ; m = 7 ; x = 2 ; (Cal.H * f)(x) / f(x) = 53.9999972207
n = 3 ; m = 6 ; x = 2 ; (Cal.H * f)(x) / f(x) = 41.4999972198
n = 4 ; m = 5 ; x = 2 ; (Cal.H * f)(x) / f(x) = 36.5000043004
n = 5 ; m = 4 ; x = 2 ; (Cal.H * f)(x) / f(x) = 28.0000043004
n = 6 ; m = 3 ; x = 2 ; (Cal.H * f)(x) / f(x) = 24.9999862338
n = 7 ; m = 2 ; x = 2 ; (Cal.H * f)(x) / f(x) = 20.4999862334
n = 8 ; m = 1 ; x = 2 ; (Cal.H * f)(x) / f(x) = 19.4999826454
n = 9 ; m = 0 ; x = 2 ; (Cal.H * f)(x) / f(x) = 18.9999826456
</code>```


问:
Clustering similar words and mapping to numbers in Python
I'm familiar with kmeans to cluster data points, but not text.. So I have one column of words in cvs format, which I want to cluster those which have similar word or more, and then mapping those cluster to numbers as index, those index numbers needed to be added as a second column. I know there is scipy packahes and word2vec in python, but this is the first time for me dealing with clustering text.. Any idea on how to do this??</p>

ContentLicense: CC BY-SA 3.0


问:
How come the use of delay differential equations in model parameter estimation better than ordinary differential equations?
in systems biology why is the use of delay differential equations better than ordinary differential equations i.e. compartment models in delay modelling? is there a data independence angle in models parameter estimation?</p>


问:
Second and Higher Order Order Corrector in Spectral Deferred Correction
I am trying to work out a second order or higher order correction for the method of Spectral deferred Correction (SDC). Specifically using as a corrector a second order or third order multi-step.

In SDC, the timestep $[t_{n},t_{n+1}] $ is subdivided in m quadrature nodes $[t_{n,1},t_{n,2},...,t_{n,m-1},t_{n,m}] $.

The ODE problem is written in its integral form:
$$ \begin{equation}
\label{eq:ODEPicard}
\begin{array}{ccl}
\phi(t) & = & \phi_{0} + \displaystyle\int_{0}^{t} F(\tau,\phi(\tau))d\tau
\end{array}
\end{equation}$$

And later rewritten in function of the error.

$ \begin{equation}
\begin{array}{ccl}
\label{eq:CorrectionEmeasure}
\delta(t) & = & \phi_{0} + \displaystyle\int_{0}^{t} [F(\tau,\tilde{\phi}+\delta(t)) - F(\tau,\tilde{\phi})]d\tau +   E(t,\tilde{\phi})
\end{array}
\end{equation}$

where $\tilde{\phi} $ denotes the approximate solution. This previous eq is discretized and based on the definitions of the error $\delta(t)$ and the residual $E(t,\tilde{\phi})$ can be discretize by a r-multi-step method as:

$ \begin{equation}
\label{eq:CorrectionMultiStepDiff1}
\begin{array}{ccl}
\displaystyle\sum^{r}_{j=0}{\{\alpha_{j}\tilde{\phi}^{[k]}_{i-j} + \displaystyle\sum_{w=0}^{j} \alpha_{w} I^{i-j}_{i-j-1}(F(\tilde{\phi}^{[k-1]})) \}}
-  \{F(t_{i-1},\tilde{\phi}_{i-1}^{[k]} ) - F(t_{i-1},\tilde{\phi}_{i-1}^{[k-1]})\} = 0
\end{array}
\end{equation} $ \
Being $I^{i+1}_{i}(F)$ an interpolatory quadrature rule
$\begin{equation}
\label{eq:QuadratureInterpolatory}
  \int^{t_{i+1}}_{t_{i}} F(\tau,\phi(\tau))d\tau \approx I^{i+1}_{i}(F) = \int^{t_{i+1}}_{t_{i}} p(\tau)\omega(\tau)d\tau = \displaystyle\sum_{j=1}^{n} \omega_{j}F(t_{j})
\end{equation}$

Well, this as every multistep method needs an initialization in the first step. What I have done is to resolve that initialization by RK2 or midpoint rule. But furthermore, the first subinterval $[t_{n,1},t_{n,2}]$ require values of the previous r-1 subnodes values. This is a problem, because as it is pointed out in Layton [1]. $\delta(t)$ is continous in the interval $[t_{n},t_{n+1}]$ but may be discontinous in the previous one. I do not totally understand why this happen, so if someone could explain me this I would be so grateful, I think it is because of the quadrature term but I don't still see it. Well, but she say that this produces a lack of smoothness in the solution and the order or correction is only 1 per step instead of r for an r-step.  And later in that paper she proposes some solutions for this, (page 5), but it is not clear what she exactly mean by "fixed starting values" and "variable starting values".

What I have tried to keep the order of accuracy of the corrector (specially I have focused in Second order Adams Bashforth)

$
\displaystyle \beta_{1} = \frac{\Delta t_{i-1}^2}{2\Delta t_{i-2}} \\
\displaystyle \beta_{2} = \frac{2\Delta t_{i-1}\Delta t_{i-2} + \Delta t_{i-1}^2 }{2\Delta t_{i-2}}
$

$ \begin{equation}
 \phi_{i}^{[k+1]} = \phi_{i-1}^{[k+1]} + \beta_{1}[F(t_{i-1},\phi_{i-1}^{[k+1]}) - F(t_{i-1},\phi_{i-1}^{[k]})] + \beta_{2}[F(t_{i-2},\phi_{i2}^{[k+1]}) - F(t_{i-2},\phi_{i-2}^{[k]})] +I^{i}_{i-1}(\phi^{[k]})
\end{equation}
$

And the first I tried it was the use of midpoint rule in the first subinterval, but as it is pointed out in IDC [2] this lack of smoothness (this concept is also redefined in this paper), it has a degree of smoothness one, in the discrete sense of the numerical set. And later I have tried to use other initialization but also lacking of smoothness.

So my question is does anybody has an idea how to obtain a second or higher order smoothness in the discrete sense based on a combination of Multistep and other method (maybe multistage). Or how to initialize at each subinterval the multistep method not to suffer order reduction due to the lack of smoothness.

P.S. If someone could clarify me also the comment that Layton does in the page 5 of [1] about continuity of the error it would be great too!


问:
How to make the following data separable for the classification into three classes?
The figure below shows the PCA projections of inputs which are 14 meteorological features, (i.e. wind, temperature, humidity, pressure, and so on.)
I would like to use any technique to make it more separable than this, The ISOMAP method is also used instead of PCA, but it gave a non-separable distribution as well! Any suggestions would be highly appreciated.


问:
What are some important pros and contras considering using Trilinos?
I'm currently developing a large simulation that will run on many nodes, and, as far as possible, on GPUs. I'm now considering to build it using libraries from the Trilinos project. I would for sure need Kokkos, Tpetra and some linear solver libraries. I also need to deal with a mesh and with construction of FEM equations.

So far, I have the following pros and cons:

Pro:

I get a large number of libraries.

My solution will be somewhat standardized and therefore more interoperable with other libraries, software, developers, ...

There's an user community, however, of course, much smaller than for many other programming topics.

Contra:

Installation is tedious. I'm now on day 6...

According to the hands on tutorial, Tpetra seems to be not very well established yet.

Compiling for CUDA also seems to be tricky and not very mature.

The risk that Trilinos does not support something I need. It would probably be very difficult to understand and extend it by myself. Particularly GPU-stuff.

Same risk for bugs.

I'm looking for some more arguments pro and contra Trilinos, or some advices considering when you would use it and when not. I know that it might be a bit opinion based, but I would be very glad if I you could share some experience and knowledge with me.


问:
Issues with self-consistent Poisson-Schrodinger solver
I'm currently in the process of writing a self-consistent Schrodinger-Poisson solver for a device heterstructure (High Mobility Electron Transistor). The algorithm is based off of this journal(1). I am having a few difficulties however. It seems like I can either have Dirichlet BCs or Neumann BCs at the substrate (i.e. the right end of the structure). Below is a snapshot of a Dirichlet and Neumann BC, respectively where $ x_n=0 $ for the first, and $\partial x/\partial n = 0$ for the second plot. Note that the boundary in question $\partial \Omega$ is at the right end of the plot.

I'm really interested in having the boundary $\partial\Omega $ be like the plot below.

 (1) A self‐consistent solution of Schrödinger–Poisson equations using a nonuniform mesh. Tan, I‐H. and Snider, G. L. and Chang, L. D. and Hu, E. L., Journal of Applied Physics, 68, 4071-4076 (1990)

Below outlines the discretization scheme that I am using.

My plan is to simulate an AlGaAs/GaAs HEMT stack that looks almost
identical to the one in the journal that I referenced above (The
stack is in Fig. 3). There's only a few Angstroms of difference in
some of my layers. I begin with a simple trial potential $$
V_{trial}(x) = -q\phi(x) + \Delta E_c(x) $$ Initially, this is a
simple line that begins at 0.7 eV with a step function of 0.2 eV
that spans the AlGaAs region.

I then solve the Schrodinger equation using the finite difference scheme.

The function Schrodinger($V_{trial}$) gets called first.

```<code>function [prob_Density, eig_Vals] = schrodinger(phi)
p = parameters;
n = p.grid_points;
t0 = p.hbar^2 / (2*p.m_AlGaN*p.dz^2);

A = zeros(n,n);

% Define first and last rows
A(1,1) = 2*t0; A(1,2) = -2*t0; A(n-1,n) = -2*t0; A(n,n) = 2*t0;

for i=2:(n-1)  % Create Hamiltonian
    A(i,i) = 2*t0 + phi(i);
    A(i+1,i) = -t0; A(i,i+1) = -t0;
    A(i-1,i) = -t0; A(i,i-1) = -t0;
end

[eig_states, eig_Vals] = eig(A);  % Solve for eigensolutions

prob_Density = eig_states.*conj(eig_states);  % Calculate probability function

end  % End of function
</code>```

Using the result "prob_Density" from step 2, I then solve for the current distributions using $$n(x) = \sum_{k=1}^{m} \psi_k^* (x) \psi_k (x) n_k$$ where $$n_k = \frac{m^*}{\pi \hbar^2} \int_{E_k}^{\infty} \frac{1}{1+e^{(E-E_F)/kT}} dE $$

Next, I use the $n(x)$ and solve Poisson's equation. I once again use the finite difference scheme with the Poisson's equation. My Dirichlet boundary condition for x=0 is defined at the Schottky barrier height. My second boundary is of the Neumann type, where the derivative is zero.

The code for the "Charge Density" (i.e. $n(x)$) and Poisson solver are written below, respectively.

```<code>%#### FIND CHARGE DISTRIBUTION ########
function charge_density = n_x(prob_Density, eig_Vals, Ef, final_Eigstate)

p = parameters;

n = p.grid_points;
n_array = zeros(n,1);
const = (p.m_AlGaAs ./ (pi*p.hbar^2));

fermi_dirac = @(E) 1./(1+exp((E-Ef)./p.Vt)); % Unit eV

for eig_state = 1:final_Eigstate

    Ek = eig_Vals(eig_state,eig_state); % First,second,third .. eigenenergy
    nk = const .* integral(fermi_dirac, Ek, inf)*p.q;   %  [m]^-2
    n_x_arr = prob_Density(:,eig_state).*nk;    % [m]^-1 * [m]^-2 -> [m]^-3
    n_array = n_array + n_x_arr;
end

charge_density = n_array; % Electron volumetric density [m]^-3

end % End program

%######### START POISSON SOLVER ###########
function phi = poisson(n_x)

first_boundary = 0.7; % eV
second_boundary = 0;

p=parameters;
n = p.grid_points;

% Create Poisson matrix elements
A = zeros(n,n);

for i = 2:n-1
   A(i,i) = -2;
   A(i+1,i) = 1;
   A(i-1,i) = 1;
   A(i,i+1) = 1;
   A(i,i-1) = 1;
end

A(1,1) = -2;
A(1,2) = 1;
A(n,n) = -1; % -1 a result of Neumann BC
A(n,n-1) = 1;

F = zeros(n,1);

%% RHS of matrix equation (the 'B' n x 1 term) AU=B, U=A\B
h = p.dz;

for i = 2:n-1
   f = p.q*(p.ND(i) - n_x(i))./p.eps_hemt(i);
   F(i) = f*h^2;
end

f1 = p.q*(p.ND(1) - n_x(1))./p.eps_hemt(1);
fn = p.q*(p.ND(n) - n_x(n))./p.eps_hemt(n);

F(1) = f1*h^2 - first_boundary;

F(n) = fn*h^2 - second_boundary;

phi = (A\F);                      %  Voltage profile

end
%##### END POISSON SOLVER
</code>```

Finally, I use the result from the Poisson solver and return to Step 2. I once again run the function Schrodinger($V_{new}$) and iterate through this process until a self consistent solution is found. This method brings me to the very first plot you see above, and I'm aiming to simulate the third plot.


问:
Once and for all: Which FEM plattform should I use for a very large multiphysics simulation?
I'm fooling around with the decision on how to build a multiphysics simulation for too long now (also several questions in this forum). First, I thought it would be possible/necessary to write most of it from scratch and by combining various libraries. (I did this in Matlab and Python to get familiar with the model, so I basically have an idea how to do it.) Now I learned more about more integrated platforms, for example Trilinos and PETSc. Then I learn about MOOSE that pretty much claims to be the jack of all trades. And probably a dozen other platforms.

For me as an unexperienced PhD student, it seems to be almost impossible to do an educated decision about which platform to use. So I really, really need the help of the community.

Requirements:

Multiphysics and nonlinear equations

Various regions with different physics and internal boundaries to each other

Very large simulations on many computing nodes

In these days it's probably mandatory to leverage GPUs.

Geometry will be generated programmatically and it would be very convenient to have an integrated mesher that can generate meshes from parametrically defined geometries (signed distance functions for example).

Can be integrated into own software

Liberal license, like BSD or (L)GPL

If you started such a project right now, what FEM platform would you count on?


问:
Collocated Grid Navier Stokes Solver
I want to solve Navier Stokes equations on a collocated grid. Earlier, I was using a MacCormick scheme based solver where I discretized predictor step in forward differences and corrector step in backward differences as mentioned in this paper  (An explicit finite-difference scheme for simulation of moving particles, http://dl.acm.org/citation.cfm?id=1127013). This  works, but the convergence is very slow at high-grid resolutions.

In this regard, can anyone suggest me some explicit NS solvers that are not very complicated and could be converged in relatively shorter times. Also, it should be mentioned here that I am writing NS equations in artificial compressibility framework.

Title: Collocated Grid Navier Stokes Solver


问:
Manufacturing a solution for non-smooth coefficients in elliptic problems
This question is a continuation of this answer (I can't comment) If we were going to manufacture a solution for a problem with discontinuous coefficients, I understand that the solution should have discontinuous gradient in the interface normal between the coefficients discontinuity, but this would not be arbitrary for problems where the interface is strangely defined. How could we go about this? Is there an automated process to manufacture a solution with discontinuous gradient in the interface?

Also, when going from the strong form to the weak form, shouldn't we have a boundary term in the interface?
$$
\int_{\Omega_1} k_1\nabla u_1 \cdot \nabla v_1 \; d\Omega + \int_{\Omega_2} k_2\nabla u_2 \cdot \nabla v_2 \; d\Omega = \\
\int_\Omega f v \;d\Omega + \int_{\partial \Omega_1}k_1 \nabla u_1 \cdot n_1\cdot v_1\; d\partial \Omega +\int_{\partial\Omega_2}k_2 \nabla u_2 \cdot n_2 \cdot v_2 \;d\partial \Omega \;\;\forall v
$$
This term is not going to vanish unless our solution has a gradient such that the flux will be zero. Say that we have two materials $k_1$ and $k_2$ separated by $\Gamma = \partial \Omega_1 \cap \partial \Omega_2$


问:
Optimal ordering in Jacobi SVD algorithm
In Jacobi SVD algorithm as given here every pair of columns of the matrix is orthogonalized until convergence. I want to know that how does the order of selection of the pair of columns affect the convergence of the algorithm.

I am working on implementing this algorithm. In my program even though I am covering every pair of columns but still with a certain order of selection of columns the algorithm doesn't converge even after thousands of iteration. But with ordering like cyclic ordering, the algorithm converges in 10-12 iterations.

The Jacobi rotation of one pair interferes with rotation of another pair. I want to know which order is the best for converging the algorithm in the minimum number of steps.


问:
Solving system of constrained linear and non-linear equations in MATLAB
Solving system of constrained linear and non-linear equations in MATLAB

I'm solving a FEM problem in MATLAB with use of the direct stiffness method. The problem is now formulated as a system of nn linear equations:

$KU=F$ or $CF=U$ with $C = K^{-1}$.

the problem is formulated in 3-Dd and has over 200 nodes. The model is used as a deformation estimator, the force is always the input and the deformation output. So there are a lot of linear equations. Now this problem needs to be extended, 2 types of constraints need to be implemented. Some of the nodes have to be constrained so that they only can move in a circular motion around a given point. And some of the nodes have to be bound with a maximum amount of deformation.

in total there will hundreds of equations of the form equations of the form:

$a_1 X(1) + a_2 X(2) + \cdots - F(1) = 0$

about 10 equality constraints of a nonlinear form like (could also be a sinusoid):

$X(1)^2 + X(2)^2 = L$

and 10 inequality constraints of the form:

$X(1) < 30\pi/180$

What is the best method to solve a nonlinear system of equations with constraints of this type simultaneously? Is it possible to use fsolve? Or fmincon? I've been looking into solvers but I can't find a solid conclusion on what solver can handle all these equation types at once and what would be most efficient in terms of computation time.


问:
Solving a 3D (almost radial) convolution with FFT
I have a 3D integral that is almost a radial convolution of the form
$$ \int d^{3}k'h(\mathbf{k'})g(|\mathbf{k-k'}|) $$
and I am looking for a fast and efficient algorithm (e.g. FFT) to solve it numerically. In the following I give you a brief description of the functions I am using and the problems I encountered.

The functions and the radial case

Consider the Fourier transform of a spherical top hat function:
$$ \tilde{W}_{i}(k,R_{i},R_{i+1})=\frac{1}{V_{i}}\int W_{i}(r)e^{i\mathbf{k\cdot r}}d^{3}r\\=\frac{4\pi}{V_{i}k}\int_{R_{i}}^{R_{i+1}}\sin(kr)rdr\equiv\tilde{W}_{i}(k)$$
where $R_{i}$ are some radiuses. If $f(k)$ is a general function (it comes from an interpolation, so there is no analytical formula) the convolution integral
$$ h_{ij}(k)=\int d^{3}k'f(|\mathbf{k'}|)\tilde{W_{i}}\left(|\mathbf{k}-\mathbf{k'}|\right)\tilde{W}_{j}\left(|\mathbf{k}-\mathbf{k'}|\right)$$
is radial, thus it can be solved with a one-dimensional FFT algorithm.

The $k$-derivative and the convolution $d_{ij}^{(2)}(\mathbf{k})$

Considering now the $k$-derivative of one of the two window function, some terms in the integral are no longer radial:
$$\int d^{3}k'f(|\mathbf{k'}|)\tilde{W_{i}}\left(|\mathbf{k}-\mathbf{k'}|\right)\frac{\partial\tilde{W}_{j}\left(|\mathbf{k}-\mathbf{k'}|\right)}{\partial k} = \\ = \int d^{3}k'f(k')\frac{\tilde{W_{i}}\left(|\mathbf{k}-\mathbf{k'}|\right)}{|\mathbf{k}-\mathbf{k'}|}\frac{\partial\tilde{W}_{j}\left(|\mathbf{k}-\mathbf{k'}|\right)}{\partial|\mathbf{k}-\mathbf{k'}|}\left(k-\frac{\mathbf{k}\cdot\mathbf{k}'}{k}\right) \\
 \equiv kd_{ij}^{(1)}(k)-\frac{\mathbf{k}}{k}\cdot d_{ij}^{(2)}(\mathbf{k}).$$
 Particularly the integral $d_{ij}^{(2)}$   is still a convolution but it isn't radial:
$$d_{ij}^{(2)}(\mathbf{k})=\int d^{3}k'f(|\mathbf{k'}|)\mathbf{k'}\frac{\tilde{W_{i}}\left(|\mathbf{k}-\mathbf{k'}|\right)}{|\mathbf{k}-\mathbf{k'}|}\frac{\partial\tilde{W}_{j}\left(|\mathbf{k}-\mathbf{k'}|\right)}{\partial|\mathbf{k}-\mathbf{k'}|}$$

Problems

Now, solving the integral by using common integration procedures is not really doable, because the functions $\tilde{W}_{i}(k)$ oscillate really rapidly so the numerical integral takes a lot of time and the result is never really reliable and it depends a lot on the chosen precision.

On the other hand a three dimensional FFT algorithm is really memory consuming. In fact $\tilde{W_{i}}(k)$ is a rapidly oscillating function and it needs a lot of points to be well sampled (around 5000/7000 points for each dimension) thus in three dimensions the matrices become really big. But for example I don't know anything about possible methods to parallelize the process.

I also checked other possibilities like Hankel transforms, Bessel functions for the angular part or Legendre polynomials, but without any success. If it was possible to rewrite it in the form of a 2-dimensional convolution, perhaps it would also be enough.


问:
Estimating the maximum absolute value (magnitude) of the Laplacian for a given function?
As motivation, consider a function which is smooth and continuous but for some reason it is very expensive to perform routine calculations of finding the Laplacian on it (maybe because it is over a large domain, etc.).

Given an original function, say, in one dimension, is it possible to estimate (or slightly overestimate) the maximum absolute value (i.e., the magnitude) of the Laplacian of that function anywhere over its given domain (without actually computing the Laplacian at all points and just finding the maximum)?

To start discussion, perhaps this can be done my taking a subset of points along the known function and using them as a basis for the estimate.


问:
Finite differences for incompressible viscous fluid equations
I am working with the equations for incompressible viscous fluid:

$$
\partial_t \vec{\omega} + (\vec{u}\cdot\nabla)\vec{\omega} = \nu\nabla^2\vec{\omega}
$$
$$
\nabla^2 \vec{\psi} = -\vec{\omega}
$$
$$
\vec{u}=\nabla \times \vec{\psi}
$$

with an usual notation:

$\omega$ ... vorticity

$u$ ... velocity

$\nu$ ... kinematic viscosity

$\psi$ ... stream function

I need to discretize the system (by means of finite differences). Since this is very common tasks, the schemes must be somewhere available, but I was not lucky in googling this time. Could you please share a solution?


问:
Numerical solution of non-linear advection equation other than inviscid burgers
I am solving a non-linear advection equation of the form $u_t + f(u)_x = 0$ where $f(u)$ is a complicated function of $u$. I am solving this equation using a first order fully implicit scheme (backward Euler in time and upwind for space) with second order linearization using Newton's method for non-linear convergence. In order to analyze the results (whether they make sense physically) I wanted to get an idea of the general characteristics that they may exhibit (like shocks and rarefactions) as a function of the nature of $f(u)$. However, in the literature, every treatment of non-linear advection that I find ultimately ends up dealing with inviscid Burgers equation. While I understand the rationale, I would like to look at some other examples. Specifically, my flux function can be convex or concave i.e., $f^{\prime\prime}(u)$ can be $< 0$ or $> 0$ and I wanted to look at cases where the non-linear wave-speeds $f^\prime(u)$ are something other than $u$ itself. Any relevant reference is much appreciated.</p>


问:
Stiffness emerges as number of ODEs increases
I want to solve a system of ordinary differential equations with Matlab. I need this to solve a mechanical engineering related problem.
If $n$ is the number of degrees of freedom of my mechanical system, then the number of variables (and also of equations, obviously) is $n^2$.
I can easily numerically integrate the ODE for example for $n=6$. I can tell that the problem is stiff because ode45 takes a very long time, while ode15s takes seconds and gives accurate results.
For a 13 degrees of freedom system (i.e. 169 variables) instead the computations are extremely slow and that's because the step size is tiny, at the beginning is $\approx 10^{-18}$, after one minute it becomes $\approx 10^{-15}$ and it gets bigger and bigger but at some point I have to stop the computation and the integration hasn't reached $10^{-10}$ (whereas with $n=6$ I easily integrate from $0$ to $100$).
I was wondering if this is something I could expect because of the 169 variables. What can I do in this cases? I tried all the available solvers for stiff problems but it's not any better.

EDIT: I understand that the 6 degrees of system and the 13 degrees of freedom system are just different models so I can't make a comparison, but from the mechanical point of view the bigger one is not much more complex.
I'm very naive for what concerns numerical solution of ODEs and I expected that doubling the size of the problem (thus solving 4x equations) would have led to (roughly) 4x computation time. Of course it cannot be as simple as this, but I'm surprised by the difference between the two models. The computation time was "seconds" for the small system and it has become "neverending" for the big one!

EDIT: I eventually discovered what my problem was! As I said, the ODE was related to a mechanical engineering problem. The issue was that the variables represented actual physical quantities the were in different units and very different orders of magnitude! I just scaled the problem by means of appropriate scalar factors and the computation time decreased significantly! Obviously, you have to scale the solution back to the right units, but it's not a big deal. For additional information I recommend this paper http://dx.doi.org/10.1137/S0895479803434914


问:
Optimally conditioned 3-tensor factorization
I have a 3-tensor $A = A_{ijk}$ with each dimension between 9 and 25 (roughly), and an integer $n > 0$.  I would like the factor this tensor as

$$A_{ijk} = \sum_{0 \le \alpha \lt n} B_{\alpha i} C_{\alpha j} D_{\alpha k}$$

such that that $B$, $C$, and $D$ are optimally conditioned possibly rectangular matrices.  I'm not sure whether optimally should mean max, sum, or product of $\operatorname{cond}(B)$, $\operatorname{cond}(C)$, $\operatorname{cond}(D)$.

Are there algorithms to do this?  Software?  I expect it's an NP-complete problem, but $A$ is fixed and I'm willing to spend quite a lot of compute time.

Clarification: As noted in the comments, I'm specifically interested in the overcomplete case $n > \operatorname{rank}(A)$, where more terms than necessary are used in order to improve the conditioning.


问:
Implementing pressure inlet boundary condition
I am interested in implementing pressure inlet boundary conditions for the 2D compressible Euler equations.

My equation of state is an ideal gas ($p=\rho R T$) which is thermally prefect but not calorically perfect, i.e. $c_p\equiv c_p(T)$ is a non-constant functions of $T$ alone.

The boundary condition specifies total pressure $p_t$, total density $\rho_t$, and angle of attack $\alpha$ on the boundary. I would like to know how to compute a consistent state at the boundary assuming an isentropic and adiabatic process.

I see from the ANSYS Fluent manual, that in case of a perfect gas one does the following:

Static pressure on the boundary $p_b$ is set equal to the static pressure in the adjacent cell $p_i$

Mach number on the boundary is computed from the isentropic relation
$$ p_t = p_b\left(1+ \frac{\gamma-1}{2}M^2\right)^{\gamma/(\gamma-1)}$$

Static temperature on the boundary $T_b$ is computed from the isentropic relation
$$ T_t = T_b\left(1+\frac{\gamma-1}{2}M^2\right)$$

Any references on implementing this kind of boundary condition in case of a ideal gas which is solely thermally perfect is highly appreciated.

My first approach was

set the mach number on the boundary equal to that of the adjacent cell

solve the adiabatic process equation
$$ h_t = h(T) + \frac{1}{2}M^2 c(T)^2$$
for $T$ to determine the temperature on the boundary

Solve $s(T,p) = s(T_t,p_t)$ for $p$ to determine the pressure on the boundary

However, this approach leads to a flow which differs significantly from what I see in ANSYS. In fact, the computed mach numbers are much lower.


问:
Corner Transport Upwind for Linear Advection in Arbitrary Velocity Field
I need to implement a 3D version of the Corner Transport Upwind (CTU) finite volume method (in python); and so I've been reading Leveque, "Finite Volume Methods for Hyperbolic Problems" which I think is quite good. Shown below is the critical section regarding how to generalize CTU for linear advection to arbitrary velocity fields (non-constant in space).

I don't understand something about the flux updates $G_{i, j-1/2}$ and $G_{i, j+1/2}$. Suppose $v_{i, j-1/2}$ is negative and $v_{i, j+1/2}$ is positive but they have the same magnitude. Then there is transverse flux across both $G_{i, j-1/2}$ and $G_{i, j+1/2}$, and hence those flux updates will be non-zero. Further, $G_{i, j-1/2}$ will be positive (the sign flips due to the sign of the velocity right?) and $G_{i, j+1/2}$ will be negative.

These cells are then used in the flux differencing formula for cell $Q_{ij}$ as: $G_{i, j+1/2} - G_{i, j-1/2}$. Because the two terms had opposite sign before, the minus sign causes them to add together, rather than cancel.

That doesn't seem right to me. With respect to cell $Q_{ij}$, I thought these transverse terms were meant to account for parts of the wave that don't flux into the cell, but in this case the donor-cell wave traveling from $Q_{i-1, j}$ into $Q_{ij}$ is stretching on the top and bottom equally, so shouldn't the flux into $Q_{ij}$ be the same as if the transverse fluxes weren't there at all, as if $v_{i, j-1/2}$ and $v_{i, j+1/2}$ were both zero? That is, the donor cell flux into $Q_{ij}$ sweeps out the same area as if $v_{i, j-1/2}$ and $v_{i, j+1/2}$ were both zero.


问:
PDE solver for drift-diffusion with generation
I am trying to solve the drift-diffusion model with generation,

$$\frac{\partial N_e}{\partial t} = \alpha(x) \left| \Gamma_e (x,t) \right| - \frac{\partial \Gamma_e}{\partial x} (x,t)$$
$$\frac{\partial N_i}{\partial t} = \alpha(x) \left| \Gamma_e (x,t) \right| - \frac{\partial \Gamma_i}{\partial x} (x,t)$$

where

$$\Gamma_e (x,t) = -\mu_e(x) E(x) N_e(x,t) - D_e(x) \frac{\partial N_e}{\partial x} (x,t)$$
$$\Gamma_i (x,t) =  \mu_i(x) E(x) N_i(x,t) - D_i(x) \frac{\partial N_i}{\partial x} (x,t)$$

I was looking at upwinding schemes but am not sure if I can use this method when there is some diffusion? If so how would I handle the $0^{\mathrm{th}}$ and $2^{\mathrm{nd}}$ order terms?

I am trying to use finite differencing.


问:
Conservative Short Sales Portfolio Optimization with MatLab code help
I'm trying to solve the optimization problem below with MatLab, but I'm unsure of how to modify the constraints in the quadprog function (or how to add the constraints with the Portfolio object).

In the formulation below, omega is the weight of the asset in the portfolio (which we are trying to find).
The capital sigma is the covariance matrix of asset returns.
mu is the expected return.
delta is the maximum allowed deviation of the expected return of an asset i, from an estimated value mu(i).
gamma is a bound for the amount of shorting allowed.
l and h are lower and upper limits for the weights.

The problem is trying to find the optimal portfolio with weights between l and h y considering only the worst outcome for the amount of short positions that exceed gamma.


问:
BC's for intermediate velocities in Implicit Fractional Step Methods
Lately, I was reading some seminal papers on Fractional Step Algorithms and I found this one:

Kim, D., Choi, H. A Second-Order Time-Accurate Finite Volume Method for Unsteady
Incompressible Flow on Hybrid Unstructured Grids. Journal of Computational Physics, Volume 162, Issue 2, pp. 411-428 (2000)

And I am confused on as to how to use/define the boundary conditions  for the auxiliary variables  in the method proposed. Specifically, in this paper the    use of a second auxiliary velocity is proposed, $\hat{u_i}$, and a certain    boundary condition is proposed for it in equation (19). Anyway, the    proposed integration method described by equations (12)-(15) seem not to require  $\hat{u_i}$   on the boundaries; the rotational part of the momentum    equation, equation (12), is solved for $\delta \hat{u}$ which, I think, it    can be extrapolated from the cell-center to the faces (in specific, the boundary faces); and equation (13) only requires the surface integration of the left-hand side which only involves the pressure. My question is: Is the latter assertion (the extrapolation of $\delta \hat{u}$) correct? If not, How can I couple the    boundary condition $\hat{u_i} = u_i^{n+1}+O(\Delta t^2)$ with equation (12) and equation    (14)? How can be implemented?


问:
Quasi Newton taking very small steps
I have implemented a Quasi-Newton method based on the Hessian approximation. I am noticing that the algorithm takes too many iterations to converge, even though it does converge. What I am not able to understand is the really small step size ($l_{\infty} \text{norm} \le 10^{-5}$) for many of these iterations. The problem is ill-conditioned and I read that Quasi-Newton can degenerate to steepest descent for ill-conditioned problems. Is that what is happening? If so, is there any intuition about the step size being so small? I have tried the damping modification given in section 18.3 from Nocedal and Wright. However, it is not improving the step length, since, the damping criterion itself becomes very small. Is there some other principled way to ensure larger step size? Will larger step size help with faster convergence?</p>

Title: Quasi Newton taking very small steps


问:
Efficiently creating an adjacency matrix from a lattice in numpy
I have an $n$ by $m$ numpy array representing a rectangular lattice $L$, where each site contains a one or a zero, representing two different materials. I'm modelling heat flow across this lattice. The idea is that the conductivity between two lattice sites depends on the material in both sites.

To solve the heat equation I need to create from this a sparse $nm\times nm$ matrix $C$, where each of the $nm$ rows and columns corresponds to a cell in my original matrix, and where the $(i,j)^\text{th}$ element represents the thermal conductance between sites $i$ and $j$ on the lattice.

My question is whether there is an efficient way to implement this in numpy. Currently I have nested Python loops that iterate over every element of $L$ and add corresponding entries to $C$, but for large lattices this obviously gets rather slow. Is there a way to do this using numpy primitives rather than Python loops?

For clarity, my current code looks like this:

```<code>A = sp.lil_matrix((w*h,w*h))

def coords_to_index(x, y):
    return x*h+y

def setlink(p1x, p1y, p2x, p2y):
    i = coords_to_index(p1x,p1y)
    j = coords_to_index(p2x,p2y)
    ci = grid[p1x,p1y]
    cj = grid[p2x,p2y]
    if ci==cj:
        d = 1
    else:
        d = 1e-3 # conductance is much smaller between different materials
    A[i,j] = d
    return d

for x in range(w):
    for y in range(h):
        d = 0.
        if x>0:
            d += setlink(x,y,x-1,y)
        else:
            d += 1
        if x<w-1:
            d += setlink(x,y,x+1,y)
        else:
            d += 1
        if y>0:
            d += setlink(x,y,x,y-1)
        else:
            d += 1
        if y<h-1:
            d += setlink(x,y,x,y+1)
        else:
            d += 1
        i = coords_to_index(x,y)
        A[i,i] = -d

A = A.tocsr()
</code>```


问:
Stiff ODEs coupled with PDEs (computational efficiency)
I am simulating in COMSOL a system of 3 coupled PDEs (parabolic & elliptic) along with 10 stiff ODEs. In order to have the system working, I am downsizing the time step size too much to achieve convergence. I am lucky enough that I am doing 1D simulations which don't take too much time, but this drives me to think about more complex 3D cases in which the distretized system will result in hugh DOF to be solved, with the very low step size this takes too much time.

So, is there computationally efficient method to solve such system (e.g. downsizing step size for ODEs but having larger ones for PDEs. Just imagining)


问:
Choose of basis set
I am engaged with the syntheses and DFT calculations of coordination complexes. I wonder how can I set up the basis set for carried out using B3LYP method and mixed basis sets of LanL2DZ for Pd and 6-31+G(d) for all other atoms, in the use of DFT theory by Gaussian 09W. can any one please explain for me?  </p>

Title: Choose of basis set


问:
Anisotropic cover with n-cuboids
I'm working on an algorithm for which I would like to cover an $n$-dimensional unit cube by a set of $n$-cuboids (i.e., $n$-dimensional rectangles).  The size and orientation of these cuboids is determined solely by the position of their centrepoint -- so given a point in the cube, one has a set of coordinate axies based at that point, a set $a_1,\ldots,a_n$ of positive numbers, and the cuboid for that point is $-a_1 \leq x_1 \leq a_1,\ldots,-a_n \leq x_n \leq a_n$ with respect to these coordinate axes.  The cuboids may overlap, but I'd like the overlap to be small if possible (for reasons off efficiency).  The unit cube must be contained in the union of the cuboids. The dimenison $n$ is small and odd, say $n=3,5,\ldots,11$. Any idea about how one could construct such a covering?  Computational cost is a minor issue here, a solution needing serious resources is acceptable (probably inevitable).

Motivation: This problem comes from a subdivide-and-reject method for a global optimisation problem with a non-convex objective.  I have a bound on a quantity of interest whose domain of validity is an $n$-ellipsoid at each point, but the shape and orientation of the ellipsoid varies over the cube.  The idea would be to put a cuboid inside the ellipsoid (so the bound is valid across it), this cube might then be rejected or subdivided uniformly, the subdivided cuboids would then have an a smaller covering ellipsoid, an improved bound, and so on. This is the reason that the cuboids must cover the cube, so that the union of these bounds covers the entire cube.

A key point is that the anisotropy is "localy almost constant", once we get to a small enough size, the ellipsoids covering the cuboids are practically the same shape as the ones covering the subdivided cuboids.  One could express this in terms of a bound away from zero of the Jacobian of the anisotropy tensor or such-like, suffice it to say that the anisotropy is "well behaved".

Title: Anisotropic cover with n-cuboids


问:
Sparse matrix format and sparse-matrix sparse-matrix multiplication
I'm having some performance problems with my code dealing with the multiplication of big sparse matrices (stiffness and aerodynamic influence coefficient matrices).
Mainly I have to multiply such matrices frequently in my code with each other. Basically, one should of course not multiply sparse matrices but rather use matrix-vector multiplications. However, in my case, I can hardly prevent matrix-matrix multiplications.

Now...
I get sparse matrices in CSC, CSR and COO formats as an input. (Not sensible as well, but I cannot influence this interface -.-).
Are there rules of the thumb that help me decide which format the matrices in my multiplication should best have to end up with a high(er) performance code?

I know Intel's mkl_dcsrcsc (multiplies CSR-matrix with CSC-matrix) for example, cannot find a mkl_dcsccsr (multiplies CSC-matrix with CSR-matrix) on the other hand. Is it, therefore, sensible for $A \cdot B$ to convert $A_\text{CSC}$ from CSC-format into CSR-format $A_\text{CSR}$ and $B_\text{CSR}$ into $B_\text{CSC}$? Might it even be sensible to take a long way and convert $A_\text{CSC}$ to $A_\text{CSR}$ and $B_\text{CSR}$ to $B_\text{COO}$ (in order to enable mkl_dcsrcoo)?

How can I find out which format I should best use for my matrices and which functions/routines I should best use for my matrix-matrix multiplications?

ContentLicense: CC BY-SA 4.0


问:
Inner Products vs Discretizations of Functions
Let $f:(0,1)\rightarrow [0,1]$ be a $\mathcal{C}^{\infty}$ function.
We say that a function $D_r^f:\{1,...,r\}\rightarrow [0,1]$ is an $r$-discretization of $f$ if $$D_r^f(j) = \frac{1}{|a(j)-b(j)|}\cdot\int_{(a(j),b(j))} f(x)\;dx$$ where $a(j)=\frac{(j-1)}{r}$ and $b(j) = \frac{j}{r}$. In other words, we divide the interval $(0,1)$ in $r$ consecutive sub-intervals  and $D^f_r(j)$ is the average of $f$ on the $j$-th sub-interval $(a(j),b(j))$.

For any two $\mathcal{C}^{\infty}$ functions $f,g:(0,1)\rightarrow [0,1]$, we let we consider the following inner product. $$\langle f,g\rangle = \int_{(0,1)} f(x)\cdot g(x)\; dx$$
On the other hand we define the inner product on $r$-discretizations as follows.
$$\langle D_r^f,D_r^g\rangle = \sum_{j=1}^r D_r^f(j)\cdot D_r^g(j).$$

Questions:

Given an $\varepsilon \in (0,1)$, how large must $r$ be so that $|\langle D_r^f,D_r^g \rangle - \langle f,g\rangle| \leq \varepsilon$?

Is it enough to assume that $f,g$ are $\mathcal{C}^{\infty}$ functions in order to be able to prove a bound on $r$ in terms of $\varepsilon$? If not, what are the weakest conditions $f$ and $g$ must satisfy?

Obs: This question was asked at mathoverflow some days ago (with no answers). Now I realize that scicomp may be the most adequate place to ask it.


问:
How to calculate multi-class prediction probabilities using One vs All and One vs One classifiers
To be specific, let's consider Support Vector Machines.

In the 1 vs all case, one has $n$ binary classifiers where $n$ is the number of classes. In each classifier one has two labels 0, 1 and for a given test data set, one can predict $\rm Prob(0)$ and $\rm Prob(1)$.

But how to compute prediction probability $\rm Prob(i)$ of a given test data point in class $i=1,2,\cdots,n$.

In sklearn, it seems that ${\rm Prob}(i) = {\rm Prob}(1;i)/\sum_{j=1}^n {\rm Prob}(1;j)$. How to interpret it?

In the 1 vs 1 case, one has $n(n-1)$ binary classifiers. Then how to compute prediction probability of classes from binary classifiers?


问:
Weak form for elastoplastic wave propagation
I am trying to simulate elastoplastic seismic wave propagation using Fenics Solid Mechanics Application.

The app. provides some quasi-static demos to show elastoplastic behaviour in a cube/ beam/ square. The weak form they implement is taken from page 199 of this paper written by the developers themselves, and is given as:

All the terms are defined in a paragraph below the equation in the paper. The equation is very similar to the weak form obtained for linear elasticity. It must be noted that in computational plasticity, return mapping algorithms (Simo and Hughes 1998) are popular and hence the application uses a 'closest point projection' return mapping algorithm along with Newton method. The stress (sigma) in the linear form 'L' is the residual stress (or stress from previous iteration). In the LHS (or bilinear form 'a') the variable 'C' is the consistent elastoplastic tangent moduli. Other variables are Body force 'f'; traction 'h'; test function 'v'; and trial function 'u'.

So far so good. Now, for elastic linear wave propagation the weak form implemented is:

Here, the test function 'w' is the same as the test function 'v' in equation (6) and (7). The only new term is the 'intertia term', where $/rho$ is the mass density of material and $/ddot{u}$ is the acceleration.

Considering the above 2 equations, is it fair to say that I can implement an elastoplastic wave using the following weak form?

Here $C_{ep}$ is the same as $C$ in equation (6).


问:
Reaction-diffusion equations
I'm simulating a biological phenomena with reaction diffusion equations. There are multiple diffusing materials and there are some complex relations about consumption and production of such materials. What is the best way to simulate these interactions? Is finite elements method appropriate? and is there any tool, which solves the PDE (with any method) and give me distribution of chemicals at each time step so that I apply my computations on them?</p>

Title: Reaction-diffusion equations


问:
How to add a Ricker Wavelet (Mexican Hat) to a 2D/ 3D fem mesh?
I have a 2D square mesh and a 3D beam shaped mesh and I want to propagate a seismic wave in them. I am trying to simulate them using Open source FEM codes (fenics). I have left the top surface to be free because that is where we would see a response (ground motion). I have also set the other boundaries except the bottom one to zero Dirichlet condition. These will cause the waves hitting them to reflect back. Usually, these would be absorbing boundary conditions or PML but for now I keep them reflective. I believe that the input seismic excitation will be set as a Dirichlet condition on the bottom boundary (2D/ 3D). In literature, a Ricker Wavelet pulse is considered suitable to represent a seismic input.

My questions are:

Is the input source given as a Dirichlet condition on bottom boundary or is it given as a point source ricker wavelet somewhere within the mesh close to the bottom boundary?

If it is a point source ricker wavelet; how would I implement it in fem mesh using c++ or python at the point where I want it? Does anyone have a piece of code that would help? Is it added as the source term '$f$' in the '$f∗v∗dx$' part of the weak form of wave equation (where $v$ is the test function) or is it projected as an initial condition on the mesh??

Thank you very much for your guidance!


问:
Finding boundary intersection points with Cartesian grid?
I would like to find the intersection points e.g. $M_x$ and $M_y$ as in the attached figure. The boundary (solid line) is defined by the Lagrangian points.
I am working in C++ (basic - medium proficiency).

In this regard, I need a computer algorithm or code to get the intersection points both in x- and y-in a separate arrays.


问:
Size of the support
I am working on moving least square approximation. To well-define of this method, the following matrix of size $n \times m$ must be of rank $m$.
For $m=3$
$${\bf A }=\begin{pmatrix}
1&x_1&y_1\\
1&x_2&y_2\\
\vdots & \vdots & \vdots\\
1&x_n&y_n
\end{pmatrix},$$
where $x_i,y_i$ are point's coordinates lie in the supprt and $n$ is the  number of them.
Also, for $m=6$
$${\bf A }=\begin{pmatrix}
1&x_1&y_1&x_1^2&x_1y_1&y_1^2\\
1&x_2&y_2&x_2^2&x_2y_2&y_2^2\\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
1&x_n&y_n&x_n^2&x_n y_n&y_n^2
\end{pmatrix}.$$
If for the first one, there are only $3$ points in the support and these points are the vertices of a triangle, can we sure that $rank(A)=3?.$
Does anyone what the best strategy for finding support size is? For uniform distributed points usually we can take $r=1.1*\delta$  for $m=3$ and $r=2.1*\delta$ for $m=6$, where $\delta$ is the distance between points. But for a mesh obtained by the adaptive technique or non-uniform this is an important thing that should be choose carefully not too large and not too small.</p>

Title: Size of the support


问:
Dynamic Successive Over/under Relaxation (SOR) with several variables
I am solving a partial differential algebraic equation (PDAE) system which has the following dependent variables:

$f=f(X,T)$ and $g=g(T)$, along with a few others

My current method for coupling is to use fixed-point iteration with constant under-relaxation factors for the two coupled variables $f$ and $g$. The issue is that convergence (determined from the norm of the residuals) is too slow. To resolve this I have tried to implement dynamic under-relaxation using Aitken's method. However, this doesn't seem to offer any speed-up - it actually often requires more iterations to converge than using constant relaxation factors. This is the case even if I apply upper and lower bounds on the relaxation factors calculated by Aitken's method. Reducing the starting relaxation factors to small values (~0.01) also doesn't seem to help much.

My initial thoughts are that the Aitken's method is not accounting for the coupling between the two variables - the relaxations factors $\omega_f$ and $\omega_g$ are calculated with no consideration for the variation of each other. I have tried having a look at the literature, but I can only find the application of Aitken's method for dynamic relaxation for a single variable.

Any suggestions for improving the performance of Aitken's method in this case would be much appreciated


问:
Weighted Frobenius norm in BFGS
In what sense is the weighted Frobenius norm "adimensional"/"scale-invariant" for any symmetric positive definite weight matrix $W$? If we plug in a positive diagonal matrix into $W$ wee see that $||A||_W$ is not always $||TA||_W$ for a one-to-one transformation $T$.

See page 196 of Nocedal's book on numerical optimization and page 60 of Fletcher's "Practical methods of numerical optimization" for reference.

Title: Weighted Frobenius norm in BFGS


问:
Error in MD simulation with Lennard-Jones potential in C, resulting in a lack of conservation of energy
I'm working on a MD simulation using the Lennard-Jones potential and the Velocity Verlet algorithm. The code compiles and runs, but the energy for each time step is not conserved; it decreases constantly when run. The magnitude of the potential energy is much larger than that of the kinetic energy, but I haven't been able to detect where the error is.

Thank you for any help you can offer!

```<code>#include <stdio.H>
#include <math.H>
#include <stdlib.H>
struct particle
{
    double x[3],v[3],a[3];
};

double sigma, epsilon, m, dt;
#define N 1000
struct particle particles[N];

double rfinder(int id_a, int id_b)
{
    int i;
    double rtot,dist,dist2;
    dist2 = 0;
    for(i=0;i<3;i++)
    {
        dist = particles[id_a].x[i] - particles[id_b].x[i];
        dist2 += dist * dist;
    }
    rtot = sqrt(dist2);
    return rtot;
}

double calcacceleration()
{
    int b,c,i;
    double dist,r,r7,r13,sigma6,sigma12, r12, sor6, sor12, pe;
    pe = 0;
    for(b=0;b<N;b++)
    {
        for(i=0;i<3;i++)
        {
            particles[b].a[i] = 0;
        }
    }

    for(b=0;b<N;b++)
    {
        for(c=0;c<N;c++)
        {
            if(b!=c)
            {
                double dist2;
                dist2 = 0;
                for(i=0;i<3;i++)
                {
                    dist = particles[b].x[i] - particles[c].x[i];
                    dist2 += dist * dist;
                }
                if (dist2<8*8)
                {
                    r = sqrt(dist2);
                    //printf("%d %d %f\n",b,c,r);
                    double rinv = 1.0/r;
                    double r2 = rinv*rinv;
                    double r4 = r2*r2;
                    double r6 = r4 * r2;
                    r7 = r4*r2*rinv;
                    r13 = r7*r4*r2;
                    r12 = r6 * r6;
                    sigma6 = sigma*sigma*sigma*sigma*sigma*sigma;
                    sigma12 = sigma6*sigma6;
                    sor6 = sigma6*r6;
                    sor12 = sigma12*r12;
                    for(i=0;i<3;i++)
                    {
                        dist = particles[b].x[i] - particles[c].x[i];
                        particles[b].a[i] += -(0.5 * (4 * epsilon *((-12 * sigma12 * r13) +(6 * sigma6 * r7))) * ((dist/r))/m);
                        particles[c].a[i] +=  (0.5 * (4 * epsilon *((-12 * sigma12 * r13) +(6 * sigma6 * r7))) * ((dist/r))/m);
                        pe += 0.5*(4*epsilon*(sor12 - sor6));
                    }
                }
            }
        }
    }
    return pe;
}

void position_mover()
{
    double x, x0, v;
    int p, i;
    for(p=0;p<N;p++)
    {
        for(i=0;i<3;i++)
        {
            particles[p].x[i] = particles[p].x[i] + (particles[p].v[i] * dt);
        }
    }
    return;
}

double velocity_mover()
{
    double ke;
    int p, i;
    ke = 0;
    for(p=0;p<N;p++)
    {
        for(i=0;i<3;i++)
        {
            particles[p].v[i] = particles[p].v[i] + (0.5 * particles[p].a[i] * dt);
            ke += 0.5 * m * particles[p].v[i] * particles[p].v[i];
        }
    }
    return ke;
}

void output_to_file()
{
int p;
FILE * fp;
fp = fopen("data.xyz", "a");
fprintf(fp,"%d \n\n",N);
for(p=0;p<N;p++)
{
    fprintf(fp,"H %f %f %f\n",particles[p].x[0],particles[p].x[1],particles[p].x[2]);
}
fclose(fp);
}

int main()
{
int p, i, j, k;
double ft,t,pe, ke,te;
sigma = 1;
epsilon = 1;
FILE * fp;
fp = fopen("data.xyz", "w");
fclose(fp);
printf("What timestep would you like?\n");
scanf("%lf", &dt);
printf("In seconds, what is the final time?\n");
scanf("%lf", &ft);
printf("What mass would you like?\n");
scanf("%lf", &m);
p = 0;
//sets initial positions of atoms to be evenly spaced in a 10 * 10 * 10 cube
for(i=0;i<10;i++)
{
    for(j=0;j<10;j++)
    {
        for(k=0;k<10;k++)
        {
            particles[p].x[0] = i*2;
            particles[p].x[1] = j*2;
            particles[p].x[2] = k*2;
            p++;
        }
    }
}

/*sets initial velocities and accelerations of atoms to 0 to get rid of any
garbage in memory */
for(i=0;i<3;i++)
{
    for(p=0;p<N;p++)
    {
        particles[p].v[i] = 0;
        particles[p].a[i] = 0;
    }
}

for(t=0;t<ft;t+=dt)
{
output_to_file();
calcacceleration();
velocity_mover();
position_mover();
pe = calcacceleration();
ke = velocity_mover();
te = ke + pe;
printf("Time: %f\n PE: %f\n KE: %f\n TE: %f\n",t,pe, ke, te);
}
printf("Done!\n");
return 0;
}
</code>```


问:
temperature set point maintenance with Langevin thermostat + shake algorithm
During the past ~4.5 years of my PhD, I have written an MD code, based largely on LAMMPS and NAMD, however reparallelized to suit my computational scenario better. One important thing to mention possibly for this discussion is that the code runs in SPDP mode to take advantage of the wide vector registers on xeon phi coprocessors. I am currently working with acetaminophen in water, however I cannot seem to rationalize the behavior I am observing with the temperature vs time data. I am using GAFF force field for acetaminophen and tip3p for water. The langevin thermostat is being used for temperature maintenance, and is being ran with a relaxation time of 2 ps. The shake algorithm is being used to constrain the motion of all bonded hydrogens. Currently, the system is a 2.75nm acetaminopehn nanocrystal embedded in a previously equilibrated box of acetaminopen + tip3p water. These are fairly large simulation I am trying to run, and consist of about 2600 acetaminophen molecules in about 30K water molecules for about 150K total atoms.  I have noticed that the system equilibrates at ~321 K, even though the set point is 320 K. I have figured out through testing that the set point offset is somehow related to the shake algorithm being applied to the solid nanocrystal. Currently,  I have the following test cases. I have more plots to show, but the website wont let me post more than 2.

solid liquid shake : this is the case where shake is applied to all hydrogens in the system (solid acetaminophen nanocrystal + liquid acetaminophen +. The system was ran in the NPT ensemble with langevin thermostat (2 ps relaxation time) and berendsen barostat (1 ps relaxation time). The timestep was 1.0fs. As you can see from the first image, the system equilibrates at ~321 K, and is clearly above the 320 K set point.

fluid-system shake: To test if there is a bug in my shake code that appears only when applied to acetaminophen, I ran an equilibration of a fluid mixture of acetaminophen and water. This system has no embedded solid sphere, and was generated using the insert-molecules and solvate command from gromacs. The following image shows the system equilibrating at 320 K correctly. Not sure what is so specific about the embedded sphere case thats making it bug out.

one final thing of note is that I have tried varying the friction constant of the langevin thermostat. Increasing the friction accelerations temperature convergence to the steady state value, but does not affect the steady state value itself (i.e. 321 K)

My specific questions are (1) Do others share any concern regarding the system equilibrating at a slightly higher (~0.8-1.0 K) temperature than the set point. I.E. Is it a general phenomenon for the shake algorithm to possible induce temperature set point offset? Is this kind of offset even considered strange? (2) Does anyone see something I am missing? I have been pounding my head on a desk for about a week now. Do others have thoughts as to why the shake algorithm, when applied to my embedded solid nanocrystal, causes the system to miss the set point, but not when the shake algorithm is applied to an acetaminophen + water fluid mixture? Any help is greatly appreciated. I am the only one in my research group who does MD, and cant ask these kind of very detailed questions to anyone else.


问:
explicit scheme stability restriction
Looking at the plain heat equation $u_t=u_{xx}$ the explicit scheme for it would look like the following iteration:
$$u_{m,n+1}=\rho u_{m-1,n}+(1-2\rho)u_{m,n}+\rho u_{m+1,n}$$
I noticed this equation resembles some probability equation: basically If I look at it as:
$$u_{m,n+1}=p_1 u_{m-1,n}+p_2u_{m,n}+p_3u_{m+1,n}$$
I notice $p_1+p_2+p_3=1$ and if I restrict those numbers to be between $[0,1]$ if as in case of probabilities it matches the stability criteria. Is there some sort of connection between probability and FDM here? It looks to me as $u_{m,n+1}$ is some sort of expectations but I can't make my argument complete so I would appreciate some help on that.</p>

Title: explicit scheme stability restriction


问:
FEM libraries with trace spaces
To implement hybridizable discontinuous Galerkin methods, one needs finite element spaces defined on the skeleton of the mesh. deal.II has support for HDG through FE_FaceQ class which provides discontinuous Lagrange functions supported on the skeleton of the mesh. I want to know if there are other libraries which provide such FE spaces for "embedded DG" (EDG) and "interior embedded DG" (IEDG) methods. EDG methods use continuous trace spaces while IEDG methods use trace space continuous on interior edges but discontinuous on boundary edges.</p>

Title: FEM libraries with trace spaces


问:
Spectral/hp-finite elements for 4th order PDEs
Does anyone know of references discussing the solution of 4th order PDEs by way of spectral/hp-finite element methods?  Specifically, I'm interested in the extension of the spectral element method, which uses high-order Lagrange polynomial interpolants (at a suitable distribution of nodes) as shape functions.  In this case, the resulting finite elements are C0 continuous, which is sufficient for 2nd order PDEs.

For fourth order PDEs, C1 continuous finite elements are required. Hence, cubic (or quintic) Hermite interpolation is typically used to construct the shape functions, as they produce finite elements with this property.  I'm interested in how one might best construct higher order "spectral" elements that are C1 continuous and have the fast convergence properties associated with the spectral element method.

I'm aware that an alternative approach would be to split the 4th order PDE into two 2nd order PDEs.  However, I think this would require solving for the 2nd derivative of the unknown, which may become infinite at points on the boundary in my case.  Thus, I'm interested in methods where I need only solve for the unknown and its first derivative to avoid this trouble.


问:
Oscillations in Chorin's method due to the BC
I am pretty new to the CFD and I wanted to start with Chorin's projection. The starting problem is just a free jet flowing in the investigated area. I got terrible oscillations almost immediately and I think it is due to the wrong treated boundary conditions for intermediate velocity $u^*$. There is just a noob mistake somewhere in the process:

Calculation of the $u^*$ using:

$$
u^* = u^{n}+\Delta t(u^n\cdot\nabla)u^n
$$

Solving the Poisson equation with the source term $\frac{\Delta t}{\rho}\nabla \cdot u^*$ and then updating the velocity:

$$
u^{n+1} = u^* + \frac{\Delta t}{\rho}\nabla p^{n+1}
$$

Where in the process should the BC for $u^*$ step in?

How can I use a boundary condition $u^*=u^{n+1}$ if I don't have the $u^{n+1}$ yet?

I need zero Neumann BC for velocity and zero Dirichlet BC for pressure.


问:
What is some of the best software design you've seen in a numerical code?
I don't think I've ever seen any scientific computational software that I would consider "good" in terms of software design. I don't mean in terms of functionality, but in terms of good high-level design, coding standards, readability, etc. Most things I've encountered have some degree of bad code smells.

What are some examples of really good computational science software engineering? Why do you like it? Preferably open source so we can see it, any language.

To minimize the subjective nature of the question, and because there's no one "right way" to do things, please give concrete examples of things like software design patterns, agile development practices, regression testing, and similar things implemented in scientific codes.

I like this post titled "The Exceptional Beauty of Doom3's Source Code." I'd love to see a similar example in the scientific world.


问:
Mesh partitioner with user-defined overlap
I am looking for a mesh partitioner, where I can specify overlap, for example h = 3.
I have looked into metis, but I wasn't able to find such a functionality.

Is there any other package which provides this functionality?


问:
Computing 3-term Connection Coefficients for Wavelets
I am trying to calculate the three-term connection coefficients

$$
Λ_{l,m}^{d_1,d_2,d_3} = ∫_{-∞}^∞ φ^{(d_1)}(x) φ^{(d_2)}_l(x) φ^{(d_3)}_m(x) dx
$$

for Daubechies wavelets numerically using Python. This report (p. 20) gives a nice overview of how to calculate them.

I wrote a (well documented) implementation in python and for the two-term connection coefficients had no problems but when I compare my 3-term results with the ones of Resnikoff and Wells (1998) I find lots of errors (only the (l=m) coefficients appear to be correct).

I am starring at the code for weeks now, but I cannot find any errors. I would be very pleased if somebody would give me a hint where the error might be.


问:
Difference between Chebyshev first and second degree iterative methods
Consider linear equation $Au = f$.

We want to solve it with iterative method (assuming $A$ is good).
First order iterative method is:
$$
u^{k+1} = u^k - \alpha_{k+1}(Au^k - f),
$$
The second degree method is:
$$
u^{k+1} = u^k - \alpha_{k+1}(Au^k - f) - \beta(u^k - u^{k-1}).
$$

For both methods we can define iteration parameters $\alpha_k$, $\beta_k$ via minimax problem which solution is Chebyshev polynomials.

This is good, but it seems to me, that convergence speed is the same for both cases and is
$$
\|\varepsilon^{k}\| \leq \frac{2\sigma^k}{1+\sigma^{2k}}\|\varepsilon^{0}\|,
$$

where $u - u^k = \varepsilon^{k}$ approximation error after the $k$-th iteration and $\sigma$ is constant dependent on operator spectrum.

The only idea I have, that first order iterations optimal for chosen $k$ for which coefficients are computed, while second-order iteration is optimal on every step.

I would greatly appreciate any work on this to clear-up those details.


问:
linearly interpolate and determine gradients for data on non-uniform grid
I have measurements of a quantity on a 3d grid. My measurements are distributed on four x-y planes similar to what is shown in the image below. The measurements roughly follow a Cartesian grid but have varying spacing and gaps. The measurements are more sparse in the z-direction but I expect the quantity to vary by an order of magnitude more in the x and y direction.

The whole domain is not convex but for now I will concentrate only on the subset of data points where the measurement planes overlap.

I am interested in linearly interpolating the quantity in 3D (in the x-y region where the measurement planes overlap), integrating the quantities over all space, and determining gradients in each direction.

I am aware of trilinear interpolation along a cartesian grid (e.g. scipy.interpolate.RegularGridInterpolator) and barycentric interpolation with Delaunay triangulation (e.g. scipy.interpolate.LinearNDInterpolator). Is there an advantage to choosing one method over another?
I also found a gradient calculation in matplotlib.tri.TriInterpolator.gradient. Is this a good choice?


问:
RCWA for Non-Periodic Structure?
I have a quick question regarding Rigorous Coupled Wave Analysis (RCWA) for numerically solving problems (in optics). I am trying to simulate a device which is non-periodic in x, y directions (but is longitudinally homogenous in z-direction). Is this possible? If so, would you be able to point me toward some readings or materials (haven't been able to find any)? Thanks in advance!</p>

Title: RCWA for Non-Periodic Structure?


问:
Basis over edges of a mesh
In 2D, the hat functions are the usual function basis for the finite dimentional space $\{v\in \mathcal{C}:v|_T\in P_1(T)\,\forall T\}$: set of continuously functions that are polynomials less than or equal to one over each element (triangle) $T$ of the mesh.

But, which are the function bases for $\{v:v|_e\in P_k(e)\,\forall e\}$? (space of functions that are polynomials less than or equal to one on edge $e$ of the mesh).

In other words, how can I define a base of polynomials less than or equal to one over an single edge $e$ of the mesh?

Title: Basis over edges of a mesh


问:
Eigenvalue decomposition of the sum: $AA^T$ + diag($u$)
Suppose $A\in\mathbb{R}^{n\times c}$,$u\in\mathbb{R}^n$,$n\gg c$.
The time complexity of eigenvalue decomposing directly for matrix $AA^T+\text{diag}(u)$ is $O(n^3)$. And it is easy to avoid $O(n^3)$ for matrix $AA^T+I$. So can we avoid $O(n^3)$ for matrix $AA^T+\text{diag}(u)$? Thanks.

Regards.

Jay.


问:
Time integration for elastodynamics
I'd like to solve the elastodynamics equation for a problem with fast scales from an impact (although no shocks).

$\rho\ddot{\mathbf{u}} - \nabla \cdot \sigma(\mathbf{u}) = \mathbf{f}$

In structural dynamics, this equation is typically solved with the Newmark scheme with Lagrange elements for the spatial domain discretization. The Newmark scheme is implicit, but I believe given the fast scales, it will be cheaper to solve with an explicit scheme. To avoid inverting the mass matrix from the second-order derivative term, I have seen that a common technique is to lump the mass matrix into a diagonal matrix. What's the accuracy of this technique?

Other methods involve using DG for the spatial discretization, which gives a block diagonal matrix that can be easily inverted, could this be cheaper and more accurate than the above approach? What would be the right way to approach this problem? Thanks.

Title: Time integration for elastodynamics


问:
Square error estimate for adaptive mesh refinement
In a particular implementation(Finite volume advection using upwind) of adaptive mesh refinement the error square estimate for a cell C is given as

$$
\sum_{i = x,y,z} vol * \frac{1}{12} * h^{2} * (\frac{ds}{di})^{2}
$$

where $vol $ is the volume of the cell, $h $ is the size of the cell and ${ds}/{di} $ is the derivative of the solution $s $ in each of the direction $i$.

What I still fail to understand that where is this error expression coming from and is the error expression specific for the scheme being implemented which is first order accurate in space or, can it be applied for any other scheme too. I tried to find the truncation error in space for upwind scheme but the expression doesn't matches the expression given above. What i also tried, is to use the same expression for error estimate for a second order accurate central differencing scheme for diffusion scheme( nomal 5 point stencil in 2D or 7 point stencil in 3D) and it seems to work for that though I face feeble issues.

Any other better ideas or pointers for error estimation to dynamically refine the mesh is also welcome.


问:
Gauss Seidel moving mesh AMR hamilton Jacobi
I was trying to implement a moving mesh algorithm using Gauss Seidel, so:

I have a pde like this (Hamilton - Jacobi eq)

$\begin{cases}\phi_t+H(\phi_x)=0 & [-1,1]\times [0,T]\\\phi(x,0)=\phi_0 & t=0\end{cases}$

and I do this each step:

```<code>a) Use a standard solver to approximate solution (one step more in time)
b) Do 5 Gauss Seidel iteration to adapt grid
c) interpolate solution on the new grid
d) repeat until t=T
</code>```

Now I have some problems with (b), maybe my weight function is not properly calculated:

equidistribution principle says:

$\omega(\phi_{j+1/2})(x_{j+1}-x_{j})-\omega(\phi_{j-1/2})(x_{j}-x_{j-1})=0$

so that my GS iteration is (using $n,n+1$ as time's index)

$\omega(\phi^n_{j+1/2})(x^n_{j+1}-x^{n+1}_{j})-\omega(\phi^n_{j-1/2})(x^{n+1}_{j}-x^{n+1}_{j-1})=0$

now if I want to use:

$\omega=\sqrt{1+\alpha|\phi_x|^2}$ , $\alpha>0$

and a central differencing for derivatives $\phi_x=\frac{1}{2}\left(\frac{\Delta^+_x\phi_j}{\Delta^+_xx_j}+\frac{\Delta^-_x\phi_j}{\Delta^-_xx_j}\right)$ , what's my linear system to solve?

A=\begin{bmatrix}1 & 0 &0 &0 & \dots & 0\\\dots \\0 &\omega_{j-1/2} & -\omega_{j+1/2}-\omega_{j-1/2} & \omega_{j+1/2} & \dots & 0 \\\dots\\0& 0 & \dots & 0 &\dots &1\end{bmatrix}
x=\begin{bmatrix}x^{n+1}_1 \\x^{n+1}_2 \\\dots \\x^{n+1}_j\\\dots \\x^{n+1}_M  \end{bmatrix}
b=\begin{bmatrix}-1 \\0 \\\dots \\0\\\dots \\1  \end{bmatrix}

Ax=b? I tried this with

$\omega_{j+1/2}=\omega\left(\frac{1}{2}\left(\frac{\phi(x_{j+1})-\phi(x_{j})}{x_{j+1}-x_{j}}+\frac{\phi(x_{j})-\phi(x_{j-1})}{x_{j}-x_{j-1}}\right)\right)$

$\omega_{j-1/2}=\omega\left(\frac{1}{2}\left(\frac{\phi(x_{j})-\phi(x_{j-1})}{x_{j}-x_{j-1}}+\frac{\phi(x_{j-1})-\phi(x_{j-2})}{x_{j-1}-x_{j-2}}\right)\right)$

maybe this are wrong? my reference is this article (2.2.1 par)


问:
Stability analysis for a hyperbolic PDE on staggered grid
I am trying to understand the stability of a finite difference equation on the staggered grid.

I could understand the Von Neumann stability analysis for the collocated grid for a simple acoustic equation  as explained in this tutorial and in this  Von Newman stability analysis for 2D acoustic wave equation explicit.

However the problem is to analyse the wave equation which is presented in coupled form on staggered grid such as :

http://geodynamics.usc.edu/~becker/teaching/557/reading/Virieux1987.pdf

Above equations are better shown as in appendix A  of following

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.460.9512&rep=rep1&type=pdf

Any suggestions for stability analysis of this kind pde problems?


问:
How to form the stiffness matrix for the Poisson equation using a spectral method
This is a follow-up question to How do I form the Chebyshev differentiation matrix in MATLAB? The goal of the following code is to solve the Poisson problem:

```<code>function [PerInfError]=Poisson2Dn0dc1sev12(N)
%Solve -nabla^2u=f(x,y) on Omega=[-1,1]^2, with u(Gamma)=x^3,  leading
%to I(U) =int(0.5(nabla.u)^2-fu)dxdy
%Exact solution u(x,y) = sin(pix)sin(piy)+x^3 +y.^3
%Current example has f(x,y) = 2pi^2sin(pix)sin(piy)+6x+6y
%Using GLL nodes over 1 spectral element

%------------GLL Nodes, weights and differentiation matrix-----------------
N1=N+1;
cheb=cos(pi*(0:N)/N)';
unif=linspace(-1,1,N1)';
if N<3
    x=cheb;
else
    x=cheb+sin(pi*unif)./(4*N);
end
P=zeros(N1,N1); %N1xN1 zero matrix
xold=2;

while max(abs(x-xold))>eps
    xold=x;
    P(:,1)=1;
    P(:,2)=x;
    %set first collumn entries to 1 (P(:,1) = 1);
    %set second collumn entries to x (P(:,2) = x);
    % and so on, the collumns of this matrix are the legendre polynomials
    for k=2:N
        P(:,k+1)=( (2*k-1)*x.*P(:,k)-(k-1)*P(:,k-1) )/k;
        %%Bonnets formula ;)
    end
    % Roots of (1-x^2)L'_N
    x=xold-( x.*P(:,N1)-P(:,N) )./( N1*P(:,N1));
end

%---chebyshev differentiation matrix-------------------------------------------------
x=-x;
X=repmat(x,1,N1);%set every collumn of X to x
Xdiff=X-X'+eye(N1); % X MINUS ITS TRANSPOSE YOU IDIOT
L=repmat(P(:,N1),1,N1);
L(1:(N1+1):N1*N1)=1;
D=(L./(Xdiff.*L'));
D(1:(N1+1):N1*N1)=0;
D(1)=-(N1*N)/4;
D(N1*N1)=(N1*N)/4;
%---GLL Quad Weights-------------------------------------------------------

w=2./(N*N1*P(:,N1).^2);

%---matrix of GLL x and y coords-------------------------------------------

Xcoord = repmat(x(2:N),1,N-1)';
Ycoord = repmat(x(2:N),1,N-1);

%---defining the known exact solution uEx----------------------------------

UEx=sin(pi.*Xcoord).*sin(pi.*Ycoord);

%----- RHS contribution --- zero dirichlet boundary conditions-------------
b=zeros((N-1)^2,1);
A= zeros((N-1)^2,(N-1)^2);
for k = 1:(N-1)^2;
    if mod(k,N-1)==0;
        n = N-1;
    else
        n = mod(k,N-1);
    end
    m = 1+floor((k-1)./(N-1));

    %----Full RHS for non zero dirichlet---------------------------------------

    b(k,1)= w(m+1).*w(n+1).*(2*(pi.^2).*sin(pi.*Xcoord(m,n))...
                                      .*sin(pi.*Ycoord(m,n)));

    %--- LHS contribution---split into the three contributions-----------------
    AY1=zeros(N-1,N-1);
    AX1=zeros(N-1,N-1);

    for ii = 1:N-1;
        AX1(ii,n)=sum(w(n+1).*w(:).*D(:,m+1).*D(:,ii+1)); %for A=A1=A2
        AY1(m,ii)=sum(w(m+1).*w(:).*D(:,n+1).*D(:,ii+1)); %for A=A1=A2
    end

    AXYSum=AX1+AY1;
    A(k,:)=reshape(AXYSum,1,(N-1)^2);
end
%---Solving the resulting linear system----equating the Error--------------

UEstV= linsolve(A,b);
UEst=reshape(UEstV,N-1,N-1);
Error=abs(UEst-UEx);
PerInfError=norm(UEx-UEst,inf)./norm(UEx,inf);

end
</code>```

Here is what I think the LHS should be.

$$\int_{\Omega}\nabla u \nabla v\approx w_l\sum_{i = 0}^{N}\sum_{m = 0}^{N}w_mD_{mi}D_{mk}u_{il} + w_k\sum_{j = 0}^{N}\sum_{n = 0}^{N}w_nD_{nj}D_{nl}u_{kj}$$

Where $D$ is the differentiation matrix mentioned in the original post. I got to this result by letting $u(x,y) = \sum_{i = 0}^{N}\sum_{j = 0}^{N}u_{ij}h_i(x)h_j(y)$ and $v(x,y) = h_k(x)h_l(y)$

Now the matrix, from this, should look something like

The reason for some blocks being diagonal and others being "Full" is that the first term in the above expression only makes sense when $i=k$, this is because it represents the derivative in the x direction and the second term when $j=l$ for the same reason but in the y direction.

I would like to understand how the author has implemented this. I am being quite vain here, assuming that he has done it this way.


问:
Gradients of non-uniformly sampled data in 3D space
I have measurements of magnetic field on a 3d grid. My measurements are distributed on four x-y planes similar to what is shown in the image below. The measurements roughly follow a Cartesian grid but have varying spacing and gaps. The measurements are more sparse in the z-direction but I expect the magnetic field to vary by an order of magnitude more in the x and y direction. Each vector component of the magnetic field ($B_x$, $B_y$, $B_z$) is measured at different points.

In the image below I show the measurement points for $B_x$. The points for $B_y$, $B_z$ are distributed similarly but not at the exact same points. The whole domain is not convex but I am only looking at the region inside the black rectangle.

I would like to calculate the curl of the magnetic field to determine the current density. E.g. $j_z = 1/\mu_{0} (\partial B_y / \partial x - \partial B_x / \partial y)$.

Since the data is non-uniform and $B_x$ and $B_y$ are not measured at same locations I can not use simple finite difference relations. I have tried three approaches

interpolate the $B_x$ and $B_y$ with scipy.interpolate.gridata(..., method='linear') to a common grid and then take 1st order differences between the point and divide by the spacing.

interpolate the $B_x$ and $B_y$ with scipy.interpolate.gridata(..., method='linear') to a common grid and then take 2nd order differences between the points with gradient and divide by the spacing.

Triangulate the $B_x$ and $B_y$ positions with  scipy.spatial.Delaunay and pass the triangulation and data to matplotlib.tri.LinearTriInterpolator which has a gradient method.

I plotted the magnetic field and $j_z$ for each of the three methods. Only the third (left to right) plot seems to be a curl of the field and gives me the magnitude that I expect from the experiment.

scipy.interpolate.griddata uses scipy.interpolate.LinearNDInterpolator which triangulates the positions and uses linear barycentric interpolation. The same is true for matplotlib.tri. Are these results expected? If all methods are using linear interpolation why is the third method more accurate than taking differences of the interpolated values?


问:
Is the numerical resolution of this huge sparse polynomial system tractable?
I'd like to find numerically a solution to a sparse system of 2000000 polynomial equations of degree 3 with 50000 variables and integer coefficients (or at least to decide whether or not a solution exists).

Question: Is there an algorithm and a supercomputer for which this computation is tractable today?


问:
Best way to to find fitting parameters for time series of decaying-growing oscillator type
I have discrete time series emerging from the numerical simulations.
It means that the time series can be slightly noisy.
The time series should "obey" to the following formula:
$$ \psi(t) = \sum_{i=1}^{N} a_i \exp(\alpha_i t) \sin (\omega_i t + \phi_i) $$
where $N$ is not known in advance (it can be 1, 2, or 20) depending on a given simulation.

UPDATE 2016-08-05: Unknown parameters here are $a_i$, $\alpha_i$, $\omega_i$, $\phi_i$, and $N$, although the only thing I need is to extract $\alpha_i$ and $\omega_i$ from these time series.

Previously, when $N=1$ was guaranteed, I used nonlinear least-squares algorithm. It turned out that this algorithm requires an extremely good initial guess to find a right fit.
I'm not sure how to provide a good initial guess when $N$ is not known in advance.

What is the most robust and best way to extract parameters I need from these time series?
Thank you!

UPDATE 2016-08-05: The following is the code that demonstrates an example time series. This example is purely artificial, though, as I don't have a real time series emerging from numerical simulations right now.

alpha(1) = 0.03;
alpha(2) = 0.05;
alpha(3) = -0.01;

a(1) = 1e-5;
a(2) = 1e-5;
a(3) = 1e-5;

omega(1) = 0.9;
omega(2) = 0.125;
omega(3) = 0.001;

phi(1) = pi/2;
phi(2) = pi/3;
phi(3) = pi;

t = linspace(0, 100, 1000);
d = 0 * t;

for i = 1:3
    d = d + a(i) * exp(alpha(i)*t) .* sin(omega(i)*t + phi(i));
end

plot(t, d, 'o-');
```


问:
Book Recommondation for the time-dependent Schrödinger equation
Can you recommend a good book that discusses several methods for the numerical solution of time-dependent and time-independent Schrödinger equation? I have searched the internet several times but could not find good references.

It would be very helpful if the book also contained matlab/python examples.

I am trying to solve the TDSE in strong laser field, and trying to implement the split operator method or Crank-Nicholson method.

OwnerDisplayName: hsinghal


问:
Stable computation of $\log\sum x_i$ from $\log x_i$, with many terms
Kahan's summation algorithm is a method to compute sums:

$$\sum x_i$$

with many terms, without significant error.

I want to do this with very large numbers, and instead of the numbers themselves, I am given their logarithms, $\log x_i$ ($x_i > 0$), and instead of the sum, I want to compute the logarithm of the sum,

$$\log \sum x_i$$

Is there a variation of Kahan's summation for this case?


问:
logsumexp with one very large term and many very small terms
I want to compute an expression of the form:

$$L = \ln\sum_i e^{x_i}$$

Suppose that there are many small terms, say $e^{x_i} \approx \epsilon$. If there are $N_\epsilon$ such terms, their contribution to the sum is $\approx \epsilon N_\epsilon$.

In the usual logsumexp trick, we let $a=\max x_i$ and compute

$$L = a + \ln \sum_i e^{x_i - a}$$

to avoid overflows. However, $\epsilon$ could be so small that $\epsilon / e^a$ becomes zero in floating-point (underflow). Even then, if $N_\epsilon$ is large enough, a significant contribution $\epsilon N_\epsilon$ to the sum will be lost.

How to deal with this case? A priori I may not know about the existence of these very small terms in the sequence. Applying the logsumexp trick in this case may lead to large errors. Is there a better algorithm?


问:
Numerical method to find stationary solution of an PDE in the phase space
I'm relatively new in numerical methods for PDEs, so I have a question about solving the following PDE numerically for $W(x,p)$ in the phase space:
$$-p\frac{dW}{dx}+x\frac{dW}{dp}+\frac{d}{dx}\left(xW \right)+2 \frac{d}{dp}\left(pW \right) + \frac{d^{2}W}{dx^{2}} +\frac{d^{2}W}{dp^{2}}-\frac{d^{2}W}{dxdp} = 0.$$
What is the simplest (yet reasonably accurate) numerical method to solve the equation above? I expect to see a Gaussian in $x,p$ and I wish to visualize the numerical solution in the window $(x,p,W(x,p))$.</p>


问:
Numerical methods for calculating the inverse CDF when closed form approximation not available
I need to calculate the inverse CDF for a probability distribution, however there is no closed form approximation available in the literature.

The distribution I am working with is the Normal Inverse Gaussian distribution, however I would like to make my question more general than this. Of course, mathematical analysis could be used in principle to derive an approximation for any given distribution but I am looking for a more general technique. To give you an idea, the PDF function is

$${\frac {\alpha \delta K_{1}\left(\alpha {\sqrt {\delta ^{2}+(x-\mu )^{2}}}\right)}{\pi {\sqrt {\delta ^{2}+(x-\mu )^{2}}}}}\;e^{\delta \gamma +\beta (x-\mu )}$$

where $K_1$ is the modified Bessel function of the third kind.

The implementations I have seen for this distribution simply

Use numerical integration to define the CDF; then

Use root finding to define the inverse CDF

There must be a better way to calculate the inverse CDF than this. It's not hard to imagine improvements, for example you could pre-evaluate the PDF at a list of points, use quadrature on this list to evaluate the CDF and then use interpolation to approximate the inverse CDF.

Are there any established methods or algorithms for calculating the inverse CDF in this scenario?


问:
Updating factorization of Laplacian (add/remove edges)
For a graph $G=(V,E)$, recall that the unweighted Laplacian is $L:=D^\top D$, where $D\in\{-1,0,1\}^{|E|\times|V|}$  is the graph "gradient" operator that subtracts adjacent vertex values onto edges.

I'm writing an algorithm that adds/removes edges one-at-a-time to a graph.  In each iteration, I need to apply $L^+$ to a few vectors, where $L^+$ is the Moore-Penrose pseudoinverse of $L$ for the current graph.  I could recompute $L$ and its pseudoinverse in each iteration, but this is computationally expensive!

Adding/removing an edge from $E$ makes a rank-1 change to $L$.  So, it seems like we could update a sparse factorization (Cholesky?  LDLT?) of $L$ and use this to accelerate application of $L^+$.  Are there easy-to-use algorithms for this?

NOTES:

I've struggled to apply standard Matlab calls.  Challenges include the fact that $L$ always has a null space (constant functions), $G$ may not be a connected graph, and removing an edge is a rank-1 but subtractive update.

Extra credit if you have a pointer to existing library of code!

This paper shows how to update the pseudoinverse of the Laplacian when adding an edge.  But it does not show what to do when removing an edge (maybe it's an easy extension, but I didn't see it!), and the pseudoinverse is a dense matrix that I would like to avoid storing.


问:
Seismic Wave modelling: Elastic Wave or Acoustic Wave?
I am modelling a seismic wave equation using FEM. In the few papers that I read, I understand the following: (Kindly correct me if you disagree)

A shear (secondary wave - no change of volume) is more important to model than the compression (primary wave - no rotation) as it brings more devastation than the latter.

The elastic wave equation is a better representative of a seismic wave than the acoustic wave equation because simulating it, we propagate both 'p' and 's' waves in the mesh instead of just 'p' which we get from simulating acoustic equation. (p - primary; s- secondary)

The elastic wave equation also known as Navier's equation of elastodynamics can be multiplied by curl operator to give the Shear wave equation, or it can be multiplied by divergence operator to give the acoustic wave equation. Both shear and acoustic wave equations are very similar except for the expressions for their speed. This expression is what differentiates both the equations and is responsible for the shape that they take. I say this because I found both expressions to be exactly same except that the acoustic wave also uses Lame's parameter $\lambda$ in addition to shear modulus $\mu$ and mass density $\rho$.

My main concern is to know which waves are propagating in my mesh if I am simulating the elastodynamic or elastic wave equation. Is it both 'p' and 's'?

Thank you!


问:
Journals that publish theoretical analyses of existing algorithms
In computational physics, the vast majority of papers has the following structure

They propose some new algorithm, or improvement to an existing one

They give numerical examples, ideally comparing relevant characteristics of their new method with existing methods

I am wondering which journals would publish a theoretical analysis of an existing algorithm, e.g. "Algorithm X was previously hypothesized to have stability condition Y (but this was not proven), and this is indeed true because (math proof)".

Such a paper would deviate from the structure above in that it neither proposes a new algorithm, nor provides numerical examples. It "merely" proves that an existing algorithm has certain characteristics.


问:
Rank sufficiency concept in Finite Element
From section 19.4.1, Rank Sufficiency:

The element stiffness matrix must not possess any zero-energy
  kinematic mode other than rigid body modes.

This can be mathematically express as follows....

I have no idea what does this section about Rank Sufficiency mean. What is rank sufficiency in physical term, and why is it so important for stiffness matrix? If rank deficiency occurs, then what will happen to the stiffness matrix and the solution of it?


问:
Differential eigenproblem with eigenvalue in boundary condition

Statement of the problem

I need to (numerically) solve an  eigenproblem of the type

$$-\omega^2\mathcal{D}_1\vec{x}=\mathcal{D}_2\vec{x}$$
on the interval $[-1,1]$, where $\mathcal{D}_1$ and $\mathcal{D}_2$ are differential operators with a generalised boundary condition of the form

$$\omega^2\mathcal{D}_3\vec{x}(±1)+a~\omega\mathcal{D}_4\vec{x}(±1)+b~\mathcal{D}_5\vec{x}(±1)=0~.$$

where $\mathcal{D}_3$, $\mathcal{D}_4$ and $\mathcal{D}_5$ are more differential operators.

The critical feature here being the fact that the eigenvalue appears in the boundary condition.

My personal take on it

SLEPc offers a wide range of method for dealing with generalised eigenvalue problems including the so-called polynomial eigenvalue problem. I used a spectral representation to cast my differential problem into a regular quadratic matrix eigenvalue problem of the form

$$(A\omega^2+B\omega+C)\vec{x}=0$$

I have used row replacement to include the boundary condition into my matrices. Since the first power of $\omega$ only appears in the boundary condition, the $B$ matrix is mostly empty.

For clarity, here is an idea of what the matrices look like

I expect that, if everything is done correctly, feeding this system into a dedicated solver for quadratic eigenvalue problems should spit out the expected eigenvalues (in some easy case, my particular system has analytical solutions). However, I haven't been able to recover these.

My question

I have already spent a decent amount of time trying to find bugs in my code and will continue to do so. But in the meantime, I figured I might ask if anyone had a similar problem and had a different approach. I am especially anxious to know if perhaps there are some results unknown to me that might render my whole approach too naive and hopeless.

Any help would be greatly appreciated.

Thanks lads !


问:
Generate gaussian random field in spherical polar coordinates
I would like to generate an isotropic gaussian random field  described by a power spectrum $P(k)$ on a 3D grid which represents spherical polar coordinates (i.e., the angular separation between pixels at each slice in the radial direction is equal and each pixel represents a solid angle).

I am currently generating a random field on a cartesian grid by

Filling a 3D grid in Fourier space $(k_x, k_y, k_z)$ with random numbers drawn from a normal distribution with standard deviation $\sigma = \frac{P(k)}{2}$ where $k = \sqrt{k_x^2 + k_y^2 + k_z^2}$

Imposing hermitian symmetry to ensure a real field

Fourier transforming into real space

My only thought is that I could do this process on a high resolution cartesian grid and average into spherical bins, though this seems convoluted.

Is there a more efficient way that I could generate this directly?


问:
Direction-splitting for SSP-RK schemes
What are the implications of applying a direction-splitting within each stage of an SSP-RK scheme?

For instance, given a standard advective transport type equation:
$$
\partial_{t}Q + \operatorname{div}(\mathbf{u}Q) = 0
$$
it's popular to make use of an explicit SSP-RK style time-stepping scheme, such as the well-known SSP-RK2 method:
$$
Q^{(1)} = Q^{(n)} - \Delta t \operatorname{div}(\mathbf{u}Q^{(n)})
\\[1ex]
Q^{(2)} = \tfrac{1}{2}Q^{(n)} + \tfrac{1}{2}Q^{(1)} - \tfrac{1}{2}\Delta t \operatorname{div}(\mathbf{u}Q^{(1)})
\\[1ex]
Q^{n+1} = Q^{(2)}
$$
But if the flux terms  can only be evaluated via a direction-splitting:
$$
\operatorname{div}(\mathbf{u}Q) = \partial_{x}(uQ) + \partial_{y}(vQ)
$$
what are the implications of advancing each RK stage according to a simple direction split:
$$
Q^{(1a)} = Q^{n}
\\[1ex]
Q^{(1b)} = Q^{(1a)} - \Delta t \, \partial_{x}(uQ^{(1a)})
\\[1ex]
Q^{(1c)} = Q^{(1b)} - \Delta t \, \partial_{y}(vQ^{(1b)})
\\[2ex]
Q^{(2a)} = \tfrac{1}{2}Q^{n} + \tfrac{1}{2}Q^{(1c)}
\\[1ex]
Q^{(2b)} = Q^{(2a)} - \tfrac{1}{2}\Delta t \, \partial_{y}(vQ^{(2a)})
\\[1ex]
Q^{(2c)} = Q^{(2b)} - \tfrac{1}{2}\Delta t \, \partial_{x}(uQ^{(2b)})
\\[2ex]
Q^{n+1} = Q^{(2c)}
$$
Is the order of the method reduced? Does it loose its SSP property? Is it possible to combine such an approach with SSP-RK3 to get a third-order accurate method?

Title: Direction-splitting for SSP-RK schemes


问:
Spherical Harmonics: band-limited representations of a vector field on a sphere
I have used pyShtools in the past to expand scalar functions to spherical harmonics and to synthesize band-limited representations of them. However, I am not too sure how to achieve this for a vector field on the surface of a sphere, where $u(\theta, \phi)$ and $v(\theta, \phi)$ are the co-latitudinal and longitudinal components, respectively.

I have also used spherepack in the past for decomposing vector fields into poloidal/toroidal components; but unfortunately I couldn't find what I am looking for.

I would very much appreciate any help in pointing me in the right direction.


问:
Preconditioning matrix with known spectrum
Assume I know all eigenvalues of a matrix $A$ fall into a certain set $\Omega \subset \mathbb{C}$. Is there any way I can exploit this knowledge to design a preconditioner for $A$?

Some further remarks:

The eigenvalues fill $\Omega$ uniformly, i.e. there are no isolated eigenvalues.

I know the eigenvalues but not the eigenfunctions.


问:
Mutual information of Multivariate Gaussian Distribution
Consider N random variables $\mathbf{x}=[x_1,x_2,...,x_N]^T$ with joint density $p(\mathbf{x})$, and let $\mathbf{R}\in \mathbb{R}^{N\times N}$ denote the covariance matrix. How to derive the mutual information for these multivariate Gaussian distributions?

Or do you know which literature provide the related derivation?
Thanks a lot.


问:
How to get the eigenvalues of Hamiltonian in an over complete basis
Let $|\psi_i\rangle$, $i=1...N+m$, be a set of overcomplete basis vector in a $N$-dim Hilbert space. The following are known: (Einstein's summation convention assumed)
$$\hat{H}|\psi_i\rangle=H_{ji}|\psi_j\rangle$$,
$$\lambda^\alpha_j|\psi_j\rangle=0,\alpha=1..m$$
$$\langle \psi_i|\psi_j\rangle=S_{ij}$$
We already know the numerical values of $H_{ij}, S_{ij}, \lambda^\alpha_j$. Then how to find the eigenvalues of $\hat{H}$ using the common Matlab or Mathematica functions?</p>

OwnerDisplayName: Lagrenge


问:
Formation sketching; trouble with sampling unit optimization
I'm looking for assistance in understanding a concise portion of a paper that involves a quadratic equation, so that I can properly implement it.

There are 2 versions of the paper. Each paper utilizes a slightly different quadratic form to achieve the same goal, but I've found success with neither. Here are the papers (along with the corresponding sections) in question:

http://graphics.cs.uh.edu/website/satelliteLinks/GI2011_formationSketching_preprint.pdf

(Section 3.6)

http://graphics.cs.uh.edu/website/satelliteLinks/2011_CGA-crowd-formation-generation-preprint.pdf

(Section 'Boundary sketches'. The quadratic can be found on page 18)

My math expertise is limited. I am primarily confused when I am asked to solve the quadratic for a single value, since I'm only aware of how to solve quadratic equations in a way that there are 2 possible solutions. How do I solve the quadratic in a way that I would find a single value for b in the context of the application?

EDIT: Here are the equations and a simple description:

Given an arbitrary shape and an arbitrary integer "sampling unit", the sampling unit is used to both place points along the boundary of the shape, as well as to employ a flood fill algorithm to fill the rest of the shape, generating a formation for a group of AI agents $n$. Solving either quadratic below will yield a more optimized/desired boundary point count $b$, which is then used to compute a more optimized sampling unit for a more realistic formation based on the number of agents $n$.

I'm looking to understand how to properly solve for $b$.

$b$ = Desired number of boundary agents

$n$ = Number of AI agents

Version 1:

$b{}'$ = Number of points on the boundary of the initial formation

$f{}'$ = Number of points inside the boundary of the initial formation

$b^{2}+(2\times (n+\frac{\sqrt{f{}'}}{b{}'}))\times b +n^{2} = 0$

Version 2:

$r$ = Ratio between boundary points and inside points of the initial formation

$b^{2}+r^{^2}\times b-r^{2}\times n=0$

NOTE: I'm aware this may not be the right place for this question. I'd appreciate references to any other more appropriate places.


问:
A linear polynomial that is the minimax (best) approximation to $f(x)$ on the interval $[-1,1]$
Suppose you have the following analytical information about a function $f(x)$ on $[-1,1]$:

\begin{align*}
f(x) &= \frac{1}{x+3}\\
f'(x) & = -\frac{1}{(x+3)^2}\\
f^{''}(x) &= \frac{2}{(x+3)^3}
\end{align*}

We are asked to find the linear polynomial that is the minimax (best) approximation to $f(x)$ on the interval $[-1,1]$

Thus we must find three points that have equal by alternating sign of the error. Let $q_1(x) = \alpha x + \beta$. Since $f(x)$ is monotonically increasing on $[-1,1]$ i.e. $f'(x) < 0$ we therefore I am told that we need to look for an interior point $c$ such that $$e'(c) = f'(c) - \alpha = 0 = -(c+3)^{-2}$$

I have been trying to understand why this is the case as in why do we need to find an interior point when $f'(x) < 0$? What would happen if $f(x)$ was monotoncally increasing instead?

Any suggestions or references are greatly appreciated. I am using the Numerical Analysis by Quarteroni, et al. but this particular section is short and does not go into much detail.


问:
Find alternative path closest to original
I have a set of points in a 2D space. I want to connect the outer points so I get the convex hull. The problem here is that there is a limit to the distance between two points. Let me clarify that with a picture:

The green lines are lines that are good. The distance between the two points of the convex hull doesn't exceed the maximum distance (the outer circle of a point). But the red lines aren't correct. They exceed the maximum distance.

What I need is a line, through multiple points, correctly drawn (e.g. all green lines), as an alternative for the red lines, while keeping the path a kind of convex hull.

I've been looking in concave hulls, and that might do the trick, but I haven't been able to find a good, preferably in Javascript, algorithm.