问:
Condition number of finite element stiffness matrix
In the FEM, there are certain applications in which the condition number $(\kappa)$ of the overall global stiffness matrix $(A)$ is computed or reduced by some preconditioner $(P)$. Are there any applications where the condition number of the tangent stiffness matrix for a single finite element (or some subregion) in the mesh would be useful?


问:
How to increase the stability of a DAE solver?
I am trying to solve a set of linear PDEs of the form
$$F\left(\vec{y},\frac{\partial \vec{y}}{\partial x},\frac{\partial^2 \vec{y}}{\partial x^2},\frac{\partial \vec{y}}{\partial t}\right)=0.$$
To solve it I am using the method of lines where I discretize in $x$ to get a set of differential-algebraic equations (DAEs) which I solve using the IDA solver shown in this tutorial. It seems to be working okay, however, it is highly numerically unstable. I find that if I set atol=rtol=1e-1then the code works okay. If I reduce the error tolerance anymore the solution blows up to infinity. Do you know how I can increase stability? I am very new to solving DAEs so any advice is appreciated.

Edit:

The equations are given, here:
$$\frac{\partial u}{\partial t} = -\left[i\omega b + ik_\perp v - \sin\alpha\frac{\partial v}{\partial x}\right],$$
$$\frac{\partial b}{\partial t} = -\frac{i}{\omega}\left[\cos^2\alpha\frac{\partial^2 u}{\partial x^2}+\frac{\omega^2u}{v_A^2(x,t)}\right],$$
$$\cos^2\alpha\frac{\partial^2 v}{\partial x^2}+\frac{\omega^2v}{v_A^2(x,t)}=i\omega\left[ik_\perp b - \sin\alpha\frac{\partial b}{\partial x}\right],$$
where
$$\vec{y}(x,t)=\begin{pmatrix}
u(x,t) \\
b(x,t) \\
v(x,t)
\end{pmatrix}.$$

ContentLicense: CC BY-SA 4.0


问:
Abnormalities when using SOR to solve the Poisson Equation
I am trying to solve the Poisson equation for Heterostructures using SOR. The equation to solve looks lik

I have discretized the Poisson equation using finite difference and my code is written in FORTRAN 90. I have the following questions.

Is it correct to apply SOR to above equation (I have discretized the 2nd order differential using FD method).

My solutions are depending on the initial guess. I get one potential profile when I take initial guess as zero and a different profile when I have a different initial guess based on charge neutrality. I am not sure why I am getting this.

My solutions are also depending on the relaxation parameter. For a relaxation of 1.45, I get "expected" results but for a relaxation of 1.90, my solver converges very rapidly (I get zero residual basically) and I get wrong results. Is this behavior expected and if so how do I make sure that I always get correct results?


问:
The error in SOR algorithm suddenly falls to zero when it reaches 1e-7 range
I am solving the Poisson equation for heterojunction using Fortran90. I use the SOR algorithm to arrive at the potential profile.

I see the weird behavior where the error (the difference between the $n$-th and $(n-1)$-th iteration of the solution) abruptly falls to zero when the error reaches 1e-7 range. Is this behavior expected?

Also, the solution is depending weakly on the choice of the mesh, e.g. potential profiles are off by about 0.02 eV when I choose a coarser mesh. Is this behavior to be expected?

Let me know if I need to provide any additional details.


问:
Fortune algorithm for voronoi diagram
Although there are many algorithms to construct Voronoi diagram, some of them are faster than others. Based on my knowledge Fortune algorithm is fastest for construct Voronoi diagram either in two dimensional or three dimensional. I search many forums but can't find any Fortune algorithm(Fortune's algorithm) code or a hint for Fortran.

So my question is: Is there any free Fortran code for Fortune algorithm in two dimensions or three dimensions?

Also, any kind of hint that makes the code writing easy for me is acceptable, for example flowchart. I don't know really it right to ask free source code here, so if I do wrong forgive me.

Title: Fortune algorithm for voronoi diagram


问:
1-D Conduction Steady state heat transfer using FD
I have tried to solve 1-D Conduction Steady state heat transfer problems in Matlab (see below).

Here is the 1-D model:

```<code>T''[x] == 0, T[0] == 100, T'[100] == 0
</code>```

Analytical Solution:

```<code>T[x] = 100 for x: [0, 100]
</code>```

I am so suprised that so many iterations (2,500,000) are required for a correct solution (here:T=100 for every node)...

Is it true that we need so many iterations for such a very simple ODEs using FD? Or?

Update:

@ConvexHull: the numerical solution is equal to the analytical solution if k = 2,500,000 is used. (Just run it )

The 1-D Conduction Steady state heat transfer in Matlab code looks like:

```<code>clc
clear all

L=100;
N=1000;
dx=L/(N-1);

T = zeros(N,1);
Tb = 100; % initial temperature

k = 1500000;
T(1,1)=Tb;

for j=1:k

    for i=2:N-1

        T(i,1)=(T(i+1,1)+T(i-1,1))/2;

    end

    T(N,1) = T(N-1,1);

end

 xx= 0:dx:L;
 plot(xx,T)
 ylim([0 100])
</code>```


问:
Set of integrators do not consistently solve an equation in Python
I must solve the following second order differential equation:

$\delta \phi^{''}_{\mathbf{k}}+(3-\epsilon)\delta \phi^{'}_{\mathbf{k}}+\left(\frac{k^2}{a^2 H^2}+\frac{V_{,\phi\phi}}{H^2}-6\epsilon +4\epsilon\eta -2\epsilon^2\right)\delta \phi_{\mathbf{k}}=0;$

which is that of scalar field perturbations at linear order in cosmic inflation. $\epsilon$, $\eta$, $\phi$, $V_{,\phi\phi}$, $H$ and $a$ denote background quantities which can be solved separately (i.e., their corresponding field equations do not contain field perturbation terms). $k$ represents the Fourier mode to which the field perturbation is ascribed and can be viewed as constant.

This equation can also be written as follows:

$\frac{1}{a^3 \phi^{'2}H}\left[a^3 \phi^{'2} H \left(\frac{\delta \phi_{\mathbf{k}}}{\phi^{'}}\right)^{'}\right]^{'} =-\frac{k^2}{a^2H^2} \frac{\delta \phi_{\mathbf{k}}}{\phi^{'}}.$

When I proceed to solve the first equation above, I use $\delta \phi_{\mathbf{k}}$ to solve the second one at the same time. However, it turns out that I do not get to solve the equations consistently (i.e., $a^3 \phi^{'2} H \left(\frac{\delta \phi_{\mathbf{k}}}{\phi^{'}}\right)^{'}$ obtained from the first equation does not agree with the one got from the last one, by far).

To solve the equations, I am employing various SciPy integrator solvers (RK45, RK23, LSODA, BDF, etc...), and none of them sorts out the discrepancy successfully.

Is there any suggestion about how I could get a consistent (more accurate) solution from those equations? Let me know if further details are needed.


问:
Producing Voronoi diagram in three dimensional
A Voronoi diagram is a kind of tesselation that divided the medium into polygons in 2D and polyhedrons in 3D.

Although there are many algorithms to construct a Voronoi diagram, some of them are faster than others. Based on my knowledge, Fortune's algorithm is fastest to construct the Voronoi diagram either in two dimensional or three dimensional.

Fortune's algorithm is based on the sweep-line strategy and balanced binary search tree.

Now, I want to implement Fortune's algorithm into three dimensional. I want to know which basically variation is made to fit the algorithm to three dimensional?


问:
Understanding STRUCTURED_POINTS dataset in VTK
Based on my readings this is what I have gathered:

Dimensions: Number of data points in the x, y, z directions

Spacing: Space between points in each dimension

Origin: Location of origin

Cell data: Number of total points

I must have a misunderstanding somewhere, as the file I am working with has dimensions 13 13 13, yet my cell data is 1728. Would anyone be able to assist in pointing out my misunderstanding?

Also, following all the above declarations I see the VTK file outputs a lot of numbers. What do these numbers represent exactly?

Thank you.


问:
Quadratic optimization with nonlinear vector term
I wish to minimize the quantity $$W=1/2x^TAx-x^Tg(y)$$ with respect to $x$ and $y$, which are vectors of unknowns. $A$ is a sparse square symmetric positive definite matrix and $g(y)$ is a vector with non-linear dependence on $y$. Everything is real-valued.
Finding $x$ given $y$ is straight forward:
$$Ax-g=0$$
$$x=A^{-1}g$$
To find the optimal value of $y$, I am trying to use gradient descent (I have analytic expressions for the derivatives of $g$):

$$\nabla W(y)=-\nabla g^Tx - g^T \nabla x$$

Is the correct approach to substitute in an expression for $x$, and try the usual gradient descent tools (line searches, etc) to update $y$?
$$\nabla W(y)=-\nabla g^TA^{-1}g - g^T A^{-1}\nabla g=-2\nabla g^TA^{-1}g$$

Because I can calculate the derivatives of $g$ quite easily, I would like to try to use this information for Newton's method, or something similar. I can calculate the components of $\nabla^2g$, which I suppose is an order-3 tensor. However, I am unsure how to apply this object to a matrix or vector. I know there must be some contraction that doesn't require me to form the tensor explicitly. Is this block matrix expression correct?
$$(\nabla (\nabla g)^T)x=
    \left(
    \begin{matrix}
    \nabla^2_1g^T\ x & \nabla^2_2g^T\ x \dots \nabla^2_kg^T\ x
   \end{matrix}
\right)$$
where column k is $\nabla^2_kg^T\ x=(\frac{\partial }{\partial y_k}\nabla g^T)x$.
The dimensions seem to be compatible, though I can't tell if the simplicity is the product of my wishful thinking.


问:
What are the practical differences between a Galerkin finite element method vs an orthogonal collocation method?
What are the practical differences between a Galerkin finite element method vs an orthogonal collocation method?  The best reference I've found thus far as been a paper titled The Performance of the Collocation and Galerkin Methods with Hermite Bi-Cubics by Dyksen, Houstis, Lynch, and Rice.  To me, it seems like:

When comparable elements used, same asymptotic convergence rate

Assembly for collocation slightly easier since it can occur one row at a time (rows correspond to collocation points) rather than in Galerkin methods where an element stiffness or mass matrix is assembled and then summed into the global matrix

Bandwidth slightly lower on collocation, but a nonsymmetric system is generated even when the original problem is self-adjoint

Galerkin methods generate a symmetric system when the original operator is self-adjoint

Does this sound correct?  Are there other practical differences such as preconditioner suitability, problems with ad-hoc methods like mass lumping, flexibility enforcing boundary conditions, or other intangibles that may favor one method over another given a certain situation?


问:
Method of Lines (MOL) for coupled PDEs
I am trying to solve the set of equations posted here- Solving the heat diffusion equation with source term
One of the methods I have learnt to solve second order PDEs is the 'method of lines', where we discretize all independent variables except one. So for this example, I can discretize the spatial domain as follows-
$$\frac{\partial^2 T}{\partial x^2} = \frac{T_{i+1}-2T_i+T_{i-1}}{(\Delta x)^2} \hspace{11pt} \forall \hspace{3pt}i = 2,3,4...N-1$$
And I can write the time derivative as-
$$\frac{\partial T}{\partial t} = \frac{dT_i}{dt}$$

Ordinarily without the source term, this would have given me a system of ODEs which I could have solved using the ode45/ode15s command in MATLAB. With the addition of a non-linear source term, I now get a coupled system of PDEs. I am not sure how to incorporate the $\frac{\partial m}{\partial t}$ term in the system of ODEs formed by the variables $T_i$. Thanks in advance

Title: Method of Lines (MOL) for coupled PDEs


问:
N-dimensional optimization of moments of distribution function
Let's say I have a bag of marbles. Each marble has several attributes (color, diameter, surface roughness, weight).

I know that there is a statistical relationship between various marble attributes and weight. I am trying to optimize for the weight distribution (either maximize/minimize the mean or the variance) by only making cuts in the other parameters (e.g., pick colors and range of diameters to achieve the desired result). The cuts can be single-ended (e.g., diameter > x or diameter < y; need not be between x and y, but that's a bonus).

I don't have analytic closed-form equations. Instead, I have a large number of discrete samples. I can calculate the weight distribution numerically for the entire population, or for a sub-population based on pre-specified cuts.

Mathematically, the weight distribution $w$ can be thought of as a joint PDF of the form:

$$f_w = f(w,x_0,x_1,...x_n)$$

Where $w$ is the weight of the marble, and $x_0 ... x_n$ are other attributes. Find $x_{0,A}, x_{0,B} ... x_{1,A},x_{1,B}$ such that the following moment sum is minimized (for example, to minimize the mean weight):

$$\min\left( \sum \limits_{w} w \sum \limits_{x_0 = x_{0,A}}^{x_{0,B}} ... \sum \limits_{x_n = x_{n,A}}^{x_{n,B}} f(w,x_0,x_1,...x_n) \right)$$

The brute-force method involves picking values of $x_{0,A}, x_{0,B} ... x_{n,A},x_{n,B}$ and computing the moment sum of interest until the optimum value is found. This matrix-calculation involves overlapping subproblems. In two or three dimensions (2 or 3 parameters), it is straightforward to avoid the overlap (e.g., by using prefix sums/weighted sums). This technique gets unwieldy for N-dimensions.

Is there a well-known algorithm or group of algorithms that can achieve this?

EDIT: August 26, 2020

The function f is constructed from a fixed number of prescribed samples that are available before the calculation begins. A number of bins and/or bin sizes are selected for each of the parameters ($w$, $x_0$ ... $x_n$) available. An ($n+1$)-dimensional matrix is created (based on the selected bin/size) and the count for f[w][x0]...[xn] is incremented as each sample is processed to build the joint PDF (a normalization constant is retained separately).

For single-ended cuts, one can compute cumulative sums (even weighted cumulative sums, per above) starting from either low or high indices to the maximum index for each dimension. For example, starting at low indices and working our way up in a cumulative fashion, we are inherently making cuts at each iteration where we only retain values of x_i less than some value (where i is the dimension or i$^{th}$ parameter which runs from 0 to n).

For two parameters (two dimensions) f becomes a matrix. The fastest way to implement the computation is to use a dynamic programming approach, where a separate matrix x is kept. For each element of x, we compute the cumulative sum corresponding to the same indices of f.
So in two dimensions, the computation becomes x[i][j] = x[i-1][j] + x[i][j-1] - x[i-1][j-1] + f[i][j].  As the iteration proceeds, on can keep track of the maximum or minimum value of the sum over f for a given cut (indices i and j). Reversing the traversal order of the indices allows single-ended cuts in the other direction ("four corners" are possible with 2 parameters, with cuts below or above certain values for each of the parameters).

I've derived a formula that extends the above to four dimensions, which now involves 16 terms (not a trivial task). I'm looking for a generalization to N dimensions, or an alternative approach.


问:
Binary data clustering by Matrix factorization
I have read an article talking about binary clustering using Matrix factorization(see attached), but i would like to understand some optimization concepts:

Is it reasonable to use a Frobenius norm in such optimization?

What does it mean the centroid constraint: C^T*1=0, and why it is equivalent to K-means especially in this case (When rho is large)?

Is there other constraints that could improve clustering optimization?

Is there other optimization technique in place of DPLM(Discrete Proximal Linearized Minimization)?


问:
Openmp critical section inside if conditional statement
I have a function like this

```<code>void foo(x, *y) {
    some calculations

    if (x > *y) {
        *y = x/2;
    }

    more calculations
}
</code>```

My parallel section is

```<code>int x[] = some array;
double y=0;
#pragma omp parallel for
for ( int i  = 0; i<=n; i++) {
    do something
    for ( int j = 0; j<=n; j++) {
           foo(x[i], &y);
}
</code>```

If I rewrite the if statement as

```<code>*y = x > *y ? x/2 : *y;
</code>```

I can do

```<code>#pragma omp critical
*y = x > *y ? x/2 : *y;
</code>```

But I don’t want to change the code. How do I create a critical section for the conditional
y-update block?


问:
ode15s command for coupled PDEs
I am trying to solve the coupled PDE system given here-Solve System of PDEs, using the method of lines and the ode15s command. Referring to variables $u1$ and $u2$ as $u$ and $v$ respectively, I have discretized the x domain as-
$$\frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1}-2u_i+u_{i-1}}{(\Delta x)^2} \hspace{10pt} \forall i = 2,3,4...N-1$$ assuming there are N grid points. Discretizing the boundary conditions, I get-
$$u_1 = u_2 $$ $$v_1 = 0$$ $$v_{N-1} = v_N$$ $$u_N = 1$$ and the initial conditions would be- $$u = 1 \hspace{10pt} \forall \hspace{2pt} i=1,2,3...N$$ $$v = 0\hspace{10pt} \forall \hspace{2pt} i=1,2,3...N$$
I pack my solution vector in the form $[u \hspace{4pt}v]^T$ such that it looks like- $[u_2\hspace{4pt} u_3 \hspace{4pt}.... \hspace{4pt}u_{N-1},\hspace{4pt} v_2 \hspace{4pt} v_3 ....\hspace{4pt} v_{N-1}]$. So in total, I am solving for $2(N-2)$ variables and the values at the boundaries I can get manually. Here is the code I have written so far-

```<code>clc
clear all
tspan1 = 0:0.005:2;
xmesh1 = 0:0.005:1;

N = size(xmesh1);
N = N(1,2);

N2 = size(tspan1);
N2 = N2(1,2);

y01 = ones(N-2,1);
y02 = zeros(N-2,1);
y0 = [y01; y02];

[t,y] = ode15s(@pde2, tspan1, y0);
u1 = y(:,1:N-2);
u2 = y(:,N-1:2*(N-2));
v1 = u1(:,1);
v2 = ones(N2,1);
w1 = zeros(N2,1);
w2 = u2(:,N-2);
u1 = [v1 u1 v2];
u2 = [w1 u2 w2];

surf(xmesh1, tspan1, u2)
shading interp

function dydt = pde2(t, y)
    xmesh1 = 0:0.005:1;
    N = size(xmesh1);
    N = N(1,2);
    dx = 0.005;
    dydt = zeros(2*(N-2),1);
    dydt(1) = 0.024*(y(2)-y(1))/(dx*dx)-(exp(5.73*y(1))-exp(-11.47*y(N-1)));
    dydt(N-1) = 0.170*(y(N)-2*y(N-1)+0)/(dx*dx)+(exp(5.73*y(1))-exp(-11.47*y(N-1)));
    for i=2:N-2-1
        dydt(i) = 0.024*(y(i+1)-2*y(i)+y(i-1))/(dx*dx)-(exp(5.73*y(i))-exp(-11.47*y(i+N-2)));
        dydt(i+N-2) = 0.170*(y(N-1+i)-2*y(N-2+i)+y(N-3+i))/(dx*dx)+(exp(5.73*y(i))-exp(-11.47*y(i+N-2)));
    end
    dydt(N-2) = 0.024*(1-2*y(N-2)+y(N-3))/(dx*dx)-(exp(5.73*y(N-2))-exp(-11.47*y(2*(N-2))));
    dydt(2*(N-2)) = 0.170*(-y(2*(N-2))+y(2*(N-2)-1))/(dx*dx)+(exp(5.73*y(N-2))-exp(-11.47*y(2*(N-2))));
end
</code>```

It's a bit inefficient, but I think it should get the job done. The problem is that I get a very different surface plot from the one given in the above link for $u1$ and $u2$ (or $u$ and $v$). For comparison, the plots of $u2$ have been attached.

This is the original plot for $u2$. This is what I get-

The shape looks a bit similar but the values are quite far off. I checked the equations in the pde2 function and they looked OK to me. The plot and values for $u1$ are also very different. Could someone point out the error in this? Thanks in advance.

EDIT: Does it have something to do with stiffness of the equations? It has been mentioned that the solution changes rapidly at small $t$. However I thought that using a stiff solver would circumvent this problem. I am able to solve the 1D heat equation without source term just fine with this approach but this particular problem eludes me.

Title: ode15s command for coupled PDEs


问:
Balanced constraint
In several optimization problems, we found what we called the balanced constraint(c^T • x = z )
For exp: C^T • 1 = 0 (C is a binary Matrix)

Can I have an intuitive explanation about the concept and how it is balanced?

Title: Balanced constraint


问:
Crank-Nicolson solution of parabolic PDE with Newumann boundary conditions
I am trying to solve the non-linear parabolic PDE in $c(t,r)$
$$c_t=\frac{1}{r}(rDc_r-\alpha r^2 c)_r$$
with initial  condition $c(0,r)=f(r)$
and boundary conditions $c_r(t,r_1)=\alpha r_1c_1/D$ and $c_r(t,r_n)=\alpha r_nc_n/D$; in other words, there is no flux across the boundaries.

This the Lamm equation, a 1-dimensional problem in cylindrical coordinates. In the simplest case, $D$ and $\alpha$ are constants.

I have set up a uniformly-spaced $r$ grid with $n$ points, so $r_j=r_1+(j-1)\Delta r$, with $1<=j<=J$. I'm using the Crank-Nicolson method with forward difference for $c_t$, and central differences for $c_{rr}$ and $c_r$. At the boundaries, I have added "ghost points" $c_0$ and $c_{J+1}$.

The difference equation resulting is, for $j=1,\ldots, J$
$$\frac{c_{n+1,j}-c_{n,j}}{\Delta t}=
\frac{D}{2}\left[\frac{c_{n,j+1}-2c_{n,j}+c_{n,j-1}}{(\Delta r)^2}+\frac{c_{n+1,j+1}-2c_{n+1,j}+c_{n+1,j-1}} {(\Delta r)^2}\right]
 +\left(\frac{D}{ 2r_j}- \frac{\alpha r_j}{ 2}\right) \left[\frac{c_{n,j+1}-c_{n,j-1}}{ 2\Delta r}
+\frac{c_{n+1,j+1}-c_{n+1,j-1}}{ 2\Delta r}\right]
 -\alpha\left[c_{n,j}+c_{n+1,j}\right]$$
where now $c_{n,j}$ is the value of $c$ at the radial point $j$ and the $n$-th time step. Values for the "ghost points" are obtained from the above equation applied at $j=1$ and $j=J$ and the boundary conditions. They are then substituted into the difference equation.

This formulation leads to the matrix equation ${\bf Fc}^{n+1}={\bf Gc}^n$, where ${\bf F}$ and ${\bf G}$ are tridiagonal, and ${\bf c}^k$ is a column vector of $c_{k,j}$ at time step $k$.

Solving for ${\bf c}^{n+1}$ gives a result that does not conserve mass (i e, $c$). Since the boundaries have 0 flux, $\int_{r_1}^{r_J} r dr$ should remain constant, but it does not.

Can anyone advise me if there is something I have done incorrectly, or where to look for an error in the computations? I did not include the elements of ${\bf F}$ and ${\bf G}$, but I could make those available.


问:
2-norm of solution update suddenly becomes zero after a few iterations
I am trying to solve the Poisson equation in 2D for heterostructure devices. I have linearized the equation and discretized it using FDM. I am using BiCGStab to iteratively solve for the solution as follows:

Step 1: Solve $A^{(i-1)} x^{(i)} = b^{(i-1)}$ for $i = 1,\dots,N$ where convergence is reached

Step 2: Use $x^{i}$ to update the central coefficients of $A^{i-1}$ to get $A^{i}$ and similarly  update $b^{i-1}$ to get $b^{i}$

Step 3: If $||x^{i}-x^{i-1}||_2$, the 2-norm of the solution update, is greater than a tolerance, then go back to Step 1  to solve the new system of equations using BiCGStab. Else, exit the loop.

I am facing the following problem with this procedure:
The 2-norm of the solution update is suddenly becoming zero after a few iterations in some cases. This behaviour is leading to "false convergences", in the sense that the solutions obtained are not physical.
Similar behaviour was observed when I used SOR instead of BiCGStab. At this point, I am starting to suspect if it is wrong to use linear solvers on the Poisson equation which is a nonlinear equation (although linearized). If you could please comment on this, that would be very helpful.

Any help with this problem is greatly appreciated. Please let me know if you need any further information.


问:
Which algorithm is capable of solving a combinatorial optimization problem like this?
The problem I have at hand is a regression problem where each of $p$ inputs, $x_1, x_2, x_3, \cdots, x_p$, needs to undergo a variable transformation using one of $q$ basis functions from a set of available functions $F = \{ f_1, f_2, f_3, \cdots, f_q \}$. The cost function is the Mean Square Error (MSE) defined as $MSE = \frac{1}{n} \sum_{k=1}^n (y_k - \hat{y}_k)^2$, where $n$ is the number of samples and for each $k$, $y_k$ is an experimental value and $\hat{y}_k$ is the predicted value computed by some method (e.g., as from a Gaussian Process Regression) using the transformations ($f_j(x_i)$) as input variables. This method must be also selected from a set $M = \{ M_1, M_2, M_3, \cdots, M_m \}$ of widely different available methods. The objective is to decide which function transforms each one of the $p$ inputs and which method to use such that $MSE$ is minimum. There are no constraints whatsoever.

Any advise would be much appreciated.


问:
How would one discretize the dynamical equations for Kerr spacetime?
Using the Hamiltonian for a test particle in Kerr spacetime, we arrive at the following equations for generalized position and momenta (in natural units, $G = c = M = 1$):

\begin{align}
\dot{r} &= \frac{\Delta}{\Sigma}p_r\\
\dot{\theta} &= \frac{1}{\Sigma}p_\theta\\
\dot{\phi} &= \frac{1}{\Sigma\Delta}\left(2arE + (\Sigma - 2r)L\csc^2\theta\right)\\
\dot{p}_r &= \frac{1}{\Sigma\Delta}\left[(r - 1)((r^2 + a^2)\mu - \kappa + r\Delta\mu + 2r(r^2 + a^2)E^2 - 2aEL\right]-\frac{2p_r^2(r-1)}{\Sigma}\\
\dot{p}_\theta &= \frac{\sin\theta\cos\theta}{\Sigma}\left[L^2\csc^4\theta - a^2(E^2 + \mu)\right]\\
\dot{p}_\phi &= 0
\end{align}
where, $\kappa = p_\theta^2 + L^2\csc^2\theta + a^2 - (E^2\sin^2\theta + \mu)$, $\Sigma = r^2 + a^2\cos^2\theta$, $\Delta = r^2 - 2r + a^2$; and $\mu$, $E = -p_t$ and $L = p_\phi$ are the test particle mass, energy and angular momentum, respectively. Also, $\mu = -1$ for massive particles, while for massless particles, $\mu = 0$, and the derivatives are with respect to an affine parameter, say $\lambda$.

I want to implement a symplectic leapfrog/verlet solver to solve this system. But I don't understand, how this set of equations should be discretized. Since $p_i$ are the momenta and not velocities, does the usual Taylor series approach make sense here? If not, could someone help me understand, how a verlet/leapfrog scheme can be used here?


问:
More efficient way to calculate magnetic field using Biot-Savart
I am writing a program in python that is supposed to calculate the magnetic field along a conducting coil that is made up of a bunch of points, and the magnetic field is generated by other conducting coils which were made by rotating the points that make up the original conducting coil. I have tried doing this using Biot-Savart law (code below), but each coil is made up of ~300 points, and there are 10 coils in total, so my method is called ~810,000 times, but I think creating and subtracting 2 vectors and then crossing them with each iteration makes it take a very long time to run, so I was wondering if someone knows of a more efficient way to do Biot-Savart, or if I should try a different method. Any insight would be very much appreciated, thank you!

My Biot-Savart method:

```<code>def biotsav(l,lprev,pointpos,cur,mu0=1e-7):

    #l is the point along the wire generating the magnetic field
    #lprev is the previous point along the wire that generated a magnetic field
    #pointpos is the position where the magnetic field will be measured at
    #cur is the current running through the wire

    #if the point causing the magnetic field is the point where the magnetic field is being measured at, return 0 otherwise it would return infinity
    if l[0] == pointpos[0] and l[1] == pointpos[1] and l[2] == pointpos[2]:
        return np.array([0,0,0])

    dl = np.subtract(l,lprev) #vector going from previous point along wire to current point
    r = np.subtract(pointpos,l) #vector going from point causing magnetic field to point where magnetic field is being measured
    rmag = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2) #magnitude of r
    return mu0*cur*np.cross(dl,r)/(rmag**3) #biot savart equation
</code>```

The loop where the method is called:

```<code>#array that will hold magnetic field values along each point on the coil
Bfield = []

#Caclculates the magnetic field at a point on the TF caused by all points on all coils
for i in range(numpts): #this for loop runs through all points along original coil

    Btemp = [0,0,0] #placeholder array for updating magnetic field measured at a point along original coil

    for j in range(len(allxs)): #this for loop goes through all coils created through symmetry

        for k in range(len(allxs[j])): #this for loop goes through all points along coils created by symmetry
            currentpt = [allxs[j][k],allys[j][k],allzs[j][k]]
            prevpt = [allxs[j][k-1],allys[j][k-1],allzs[j][k-1]]
            coilpt = [xlist[i],ylist[i],zlist[i]]
            Btemp = np.add(Btemp,biotsav(currentpt,prevpt,coilpt,cur)) #adds magnetic field contributions from all points in coils made by symmetry

    Bfield.append(np.array([Btemp[0],Btemp[1],Btemp[2]])) #updates array with net magnetic field measured at point along original coil
</code>```

**Note: allxs, allys, and allzs are arrays that contain the x, y, and z coordinates of the 9 coils created by symmetry. Also, all coils have equal current running through them

For reference, the red coil is the original coil and the one I want to measure the magnetic field caused by the other 9 white coils along the original coil.


问:
Which class does this NP problem belong to?
Suppose $n$ inputs ($x_1, x_2, x_3, \cdots, x_n$) can take on any of $m$ values, say $\{ k_1, k_2, k_3, \cdots, k_m \}$, and that there is a cost function $y = f(x_1, x_2, x_3, \cdots, x_n)$. For instance, for $n = 3$ ($x_1, x_2, x_3$) and $m = 2$ ($\{ k_1, k_2 \}$) we have a total of $8$ possible triplets $S = \{ (k_1k_1k_1), (k_1k_1k_2), (k_1k_2k_1), (k_1k_2k_2), (k_2k_1k_1), (k_2k_1k_2), (k_2k_2k_1), (k_2k_2k_2) \}$. The objective is to find the triplet that minimizes $y$ without enumerating all possibilities.

I tried to relate this NP problem to various problems in combinatorial optimization (e.g., TSP, 0/1 Knapsack, Assignment Problem, etc.), but found none that could barely resemble it.

Which class does this NP problem belong to?

EDIT: The problem I have at hand is a regression problem where each of $n$ inputs, $x_1, x_2, x_3, \cdots, x_n$, needs to undergo a variable transformation using one of $q$ basis functions from a set of available functions $F = \{ f_1, f_2, f_3, \cdots, f_q \}$. The cost function is the Mean Square Error (MSE) defined as $MSE = \frac{1}{p} \sum_{k=1}^p (y_k - \hat{y}_k)^2$, where $p$ is the number of samples and for each $k$, $y_k$ is an experimental value and $\hat{y}_k$ is the predicted value computed by some method (e.g., as from a Gaussian Process Regression) using the transformations ($f_j(x_i)$) as input variables. This method must be also selected from a set $M = \{ M_1, M_2, M_3, \cdots, M_m \}$ of widely different available methods. The objective is to decide which function transforms each one of the $n$ inputs and which method to use such that $MSE$ is minimum. There are no constraints whatsoever.

Any advise would be much appreciated.


问:
Problems with simulation of a spatial filter 4f setup (Python)
I have a question about my code which computes numerically the output field of a 4f setup with a pinhole in the middle which works as a spatial filter.
My setup consists of two lenses with 50mm focal length (distance 2f), and a pinhole between both lenses. The input field has a Gaussian field distribution, where the Gaussian beam has a waist of 4mm and in addition I add to the input field some noise. My goal is to look how good I can filter out the noise for different pinhole diameters.

I wrote a program that simulates a spatial filter process in a 4f setup, and in general the spatial distribution of the output field has the shape that one would expect (Gaussian field goes in and Gaussian field comes out). The only problem that I have is that the field amplitudes of the output field are really high compared to the amplitudes of the input field (max values of fields 4 vs. 2000).

In my program I'm using the Fraunhofer approximation and solve the Fraunhofer diffraction integral (for a field starting in the front focal plane of a lens) numerically by using the FFT/IFFT functions from python. My code is inspired by the book: "Numerical simulation of optical wave propagation with examples in MATLAB" from Jason D. Schmidt (chapter 4). I guess there is something wrong with the scaling factors for the fft/ifft, but I have no idea where the error in my code is, since I'm using the exact same scaling factors like in the schmidt book.

(fft2 -> dl^2,   ifft2 -> (N*dl_f)**2,   where dl_f = 1/(N*dl))

The functions of my code are:

import numpy as np
from matplotlib import pyplot as plt

def lens_FT(U_in, lambda_0, L, N, focal):

    '''Computes the field dirstribution of a field starting in the front focal
       plane of a lens,and is then focused into the back focal plane of the
       lens using Fraunhofer diffraction.
       We assume a square shaped computational domain in this function.

    Arguments
    ---------
        U_in : 2d-array
            Initial_field dirstibution
        lambda_0 : float
            Vacuum wavelength of the initial field in mm
        focal : float
            Focal length of the lens
        L : float
            length of computational domain at one side
        N : int
            Number of grid points for one side of the computational domain

    Returns
    -------
        U_out : 2d array
            Field distribution in the Fourier plane of the initial field after
            interaction with lens
        kx : 1d array
            Coordinates in the output plane in x-directiom
        ky : 1d array
            Coordinates in the output plane in y-directiom
    '''
    pass

    # vacuum wavenumber
    k_0 = 2*np.pi/lambda_0

    # grid spacing of the initial field
    dl = L/N

    # define grid for output plane
    kx = np.fft.fftfreq(N, d=dl)
    ky = np.fft.fftfreq(N,d=dl)

    # sort arrays
    kx = np.asarray(np.append( kx[int(N/2)::], kx[0:int(N/2)] ))*lambda_0*focal
    ky = np.asarray(np.append( ky[int(N/2)::], ky[0:int(N/2)] ))*lambda_0*focal

    # define mesh for output field
    kxx, kyy = np.meshgrid(kx, ky)

    # shift U_in in order to perform the fft correct
    U_in = np.fft.fftshift(U_in)

    # compute output field with Fraunhofer integral
    U_out = (np.exp(1j * k_0 * 2*focal) / (1j * lambda_0 * focal) *
             np.fft.fftshift(np.fft.fft2(U_in)) * (dl)**2)

    return U_out, kx ,ky

def lens_inv_FT(U_in, lambda_0, L, N, focal):

    '''Computes the inverse Fourier tranform of a field which interacts with a
       thin lens using Fraunhofer diffraction.
       We assume a quadratic computational domain in this function.

    Arguments
    ---------
        U_in : 2d-array
            Initial_field dirstibution
        lambda_0 : float
            Vacuum wavelength of the initial field in mm
        focal : float
            Focal length of the lens
        L : float
            length of computational domain at one side
        N : int
            Number of grid points for one side of the computational domain

    Returns
    -------
        U_out : 2d array
            Field distribution in the Fourier plane of the initial field after
            interaction with lens
        x : 1d array
            Coordinates in the output plane in x-directiom
        y : 1d array
            Coordinates in the output plane in y-directiom
    '''
    pass

    # vacuum wavenumber
    k_0 = 2*np.pi/lambda_0

    # grid spacing of the initial field in spatial domain
    dl = L/N
    # grid spacing in spatial frequency domain
    dl_f = 1/(N*dl)

    # define mesh for output field
    x = np.linspace(-L/2, L/2, N, endpoint=False)
    y = np.linspace(-L/2, L/2, N, endpoint=False)
    xx, yy = np.meshgrid(x, y)

    U_in = np.fft.ifftshift(U_in)
    U_out = (np.exp( 1j * k_0 *2*focal) / (1j * lambda_0 * focal) *
             np.fft.ifftshift(np.fft.ifft2(U_in))* (N*dl_f)**2)

    return U_out, x ,y

def pinhole_filter_fct(cx, cy, r, N ):

    '''Creates filter function with pinhole shape

    Arguments
    ---------
        cx : int
            matrix x index where pinhole is centered
        cy : int
            matrix y index where pinhole is centered
        r : int
            radius of the pinhole  in mm
        N : int
            number of grid points of one side of the computational domain

    Returns
    -------
        filter_fct : 2d-array
            filter function with circular pinhole shape (matrix with 0 and 1 elements)

    '''
    # convert radius length into pixel
    pixel_pinhole_rad = round(radius*N/L)

    # Define matrix for filter fct
    x = np.arange(-N/2, N/2)
    y = np.arange(-N/2, N/2)
    filter_fct = np.zeros((N, N))

    # creating pinhole in
    mask = (x[np.newaxis,:]-cx)**2 + (y[:,np.newaxis]-cy)**2 < pixel_pinhole_rad**2
    filter_fct[mask] = 1

    return filter_fct

```
<p>With these functions, I then calculate the spatial filtering process in the 4f setup as follows:</p>
<pre class="lang-py prettyprint-override"><code>   # define parameter for simulation and physical parameter
   # number of grid points per side
   N = 602
   # length of one side of the computational domain in mm
   L = 15
   # Gauss width in mm
   w = 4
   # radius of pinhole in mm
   radius = 0.1

   # vacuum wavelength of initial field in mm
   lambda_0 = 0.00081
   # focal length of the lens in mm
   focal = 50

   # define initial field (gaussian beam field distribution)
   # define grids
   x = np.linspace(-L/2, L/2, N)
   y = np.linspace(-L/2, L/2, N)
   xx, yy = np.meshgrid(x, y)

   # build 2D gaussian field distribution array
   gaussian_beam_in = 3*np.exp(-(xx**2 + yy**2)/(2*w**2))

   # add noise to input field
   noise = np.random.uniform(0,1,(N,N))
   gaussian_beam_in = gaussian_beam_in + noise

   # Compute field in Fourier domain
   gaussian_beam_out, kx, ky = lens_FT(gaussian_beam_in, lambda_0, L, N, focal)

   # compute filter function
   pinhole_filter = pinhole_filter_fct(0, 0, radius, N)
   # apply filter function by elementwise matrix multiplication with field in Fourier plane
   gaussian_beam_out_filtered = np.multiply(gaussian_beam_out, pinhole_filter)

   # compute field in output plane at 4f
   gaussian_beam_out_4f, x, y = lens_inv_FT(gaussian_beam_out_filtered, lambda_0, L, N, focal)

</code>```

These are the plots of my results for a 0.1mm pinhole diameter:

Initial field

Filtered field in Fourier plane

Output field

(I can only give you these stack overflow links as I am not allowed to post images yet)

As you can see the filtering process works in my simulation, but the amplitudes of the output field are really high compared to the input field. Has someone an idea what I'm doing wrong? I'm looking for a solution of this problem for a view days now but without success. As I said before, I guess there is probably something wrong with my scaling factors for the fft.

I'm really grateful to anyone who wants to help me, because right now I have no idea what I'm doing wrong in my code.


问:
How to calculate the interior value of triangular element in edge (vector) finite element?
I was using an edge (vector) finite element to solve electromagnetic diffusion (two-dimensional cases). The element that I used was a triangular element. I have got the result of the finite element in each edges element of the triangular element. But I don't know how to extract the value interior the triangular element.

How can I get the value interior the triangular element (in $x$ and $y$ directions)??


问:
computing time scale and steady state concentration in microfluidic channels
I have been performing convection-diffusion transport studies on microfluidic channels like the following

The inlet concentration is specific and I obtain the time-dependent concentration profiles of the
solute at the discretization nodes in the domain. Using this data, I would like to know how to compute the time scale required for the system to reach steady-state concentration at the channel outlet.

image source

EDIT1: The concentration profiles are obtained after performing numerical simulations. The equation that is solved is the 1D convection-diffusion equation

EDIT2: Yes, I consider an incompressible flow. As suggested below, I could compute Peclet number(Pe) to find out whether convective transport dominates or advective transport dominates. But I am not sure whether Pe number has to be determined for each branch. Pe for mass transport is given by,
$$Pe  = \frac{u L}{D}$$

i.e is L the effective length of the whole channel?.  $D$ is the diffusion coefficient. What happens when convection-dispersion is studied? Can $D$ in the expression for Pe be replaced with the dispersion coefficient which is a function of diameter?

Is the diameter halved at every bifurcation point?
I'm sorry if the image that I've used for illustration gives that impression. In the actual system, the diameters aren't halved at a junction point (which could be a bifurcation or trifurcation point). The ratio of diameters of the branches aren't constant always. The diameters range from 2-10 micrometers and lengths range from 50-150 micrometers. The diffusion coefficient is ~1E-10 $m^2/s$

EDIT3:
From the explanations provided in the comments below and this post I understand the characteristic length scale is the smallest length scale.

For the characteristic length within the Peclet number you may use
your discretization width, or possibly the width of your channel.

I would like to understand if the choice of the characteristic length(discretization length or the diameter channel) depends on the direction of diffusive transport or is it always the smallest length irrespective of whether diffusion occurs in the radial/axial direction.


问:
How to calculate the average velocity over a surface through time using paraview?

My problem is a little more complex, but here is a simplified similar model.

Let us consider a fire source inside a box and the hot gases (burnt or unburnt) from the combustion are going out of the exhaust (white surface) at a certain velocity. We have the velocity field of the entire domain and can be visualized using ParaView. We also have density and temperature fields among others.

Now, using ParaView I would like to calculate the average velocity of the gases coming out of this surface (white) over time (unsteady regime).

Note: The individual boundaries are not loaded in ParaView. That is outlet as a separate surface, the domain as a separate surface. The entire volume is treated as a single entity with specified boundary conditions.

Any help to calculate this is appreciated.


问:
I'm trying to add a colour to the LogicalVolume in Geant4. I get the error that fVisAttributes and visAttributes are undefined. How can I fix this?
I'm trying to add a colour to the LogicalVolume which has been renamed to TrackerLogical in Geant4. I get the error that fVisAttributes and visAttributes are undefined although I have defined them in the #include section as .hh files.

The primary framework I have used is that of exampleB1 and I am trying to add color using exampleB5 's approach by applying G4VisAttributes. I have referred to this as a guide: https://github.com/Geant4/geant4/blob/master/examples/basic/B5/src/B5DetectorConstruction.cc

I have included the necessary bits in the #include section as shown in example B5, and have included a line for fVisAttributes in the section B1DetectorConstruction::B1DetectorConstruction() but I get the error messages when I run -jN after adding the following section of code before the last line ( return physWorld) :

```<code> visAttributes = new G4VisAttributes(G4Colour(0.9,0.9,0.9));
 visAttributes->SetVisibility(false);
 trackerLogical->SetVisAttributes(visAttributes);
 fVisAttributes.push_back(visAttributes);
</code>```

This is the error message I get:

```<code>home/neha/G4WORk/B1/src/B1DetectorConstruction.cc: In constructor ‘B1DetectorConstruction::B1DetectorConstruction()’:
/home/neha/G4WORk/B1/src/B1DetectorConstruction.cc:53:3: error: class ‘B1DetectorConstruction’ does not have any field named ‘fVisAttributes’
   53 |   fVisAttributes()
      |   ^~~~~~~~~~~~~~
/home/neha/G4WORk/B1/src/B1DetectorConstruction.cc: In member function ‘virtual G4VPhysicalVolume* B1DetectorConstruction::Construct()’:
/home/neha/G4WORk/B1/src/B1DetectorConstruction.cc:202:2: error: ‘visAttributes’ was not declared in this scope; did you mean ‘G4VisAttributes’?
  202 |  visAttributes = new G4VisAttributes(G4Colour(0.9,0.9,0.9));
      |  ^~~~~~~~~~~~~
      |  G4VisAttributes
/home/neha/G4WORk/B1/src/B1DetectorConstruction.cc:205:2: error: ‘fVisAttributes’ was not declared in this scope; did you mean ‘G4VisAttributes’?
  205 |  fVisAttributes.push_back(visAttributes);
      |  ^~~~~~~~~~~~~~
      |  G4VisAttributes
make[2]: *** [CMakeFiles/exampleB1.dir/build.make:89: CMakeFiles/exampleB1.dir/src/B1DetectorConstruction.cc.o] Error 1
make[1]: *** [CMakeFiles/Makefile2:104: CMakeFiles/exampleB1.dir/all] Error 2
make: *** [Makefile:130: all] Error 2
</code>```

What is causing the issue and how can I resolve this? I am still a beginner with Geant4 and may need some fundamental explanations (in simple/layman's terms preferably). Thank you!


问:
Nondifferentiable coordinate transforms
Suppose that we have coordinates $u=u(x,y)$ and $v=v(x,y)$ in $\mathbb{R}^2$ so that $v$ is not differentiable when $u(x,y)=u_0$ where $u_0$ is a constant. Can we solve a differential equation, such as the Schrödinger equation, using these coordinates? Can we compute integrals involving the Jacobian determinant in these coordinates?

Title: Nondifferentiable coordinate transforms


问:
numerical solution to pde on an ellipse
Looking for advice on discretization (preferably finite difference) schemes for pdes on curves in general, but in this case it is an ellipse (so given by $(a\cos(r), b\sin(r)$). The problem is the advection equation

$$\frac{ \partial u_s}{\partial t}+ \frac{ \partial u_s}{\partial s}=0,$$
with $s$ corresponding to arc length and the initial condition $u(s,0)=\cos^2(2\pi s/L)$. The solution should be  $u(s,t)=\cos^2(2\pi(s-t)/L)$. I don't want to do this using the closest point method, I want to keep the curve parameterized and just discretize on the curve but I'm not sure how to proceed.

Title: numerical solution to pde on an ellipse


问:
Computing Singular Value Decomposition of small ($4\times 4$) matrices
I need to compute the Singular Value Decomposition (SVD) of many $4 \times 4$ matrices. I'm looking for SVD algorithms specialized for small matrices. I've read that the std::vector sorting algorithm in the C++ STL is often hardcoded for short vectors for optimal performance. I want to hard code the SVD for the case $4 \times 4$ matrices.

What information is there about optimally solving small cases of SVD?


问:
How to take convolution of two arrays in Python by using NumPy?
Generally, we know that if we have this relation between Fourier transforms of three functions in frequency domain as:

$$\mathfrak{F}\{\mathsf{P}(t)\} = \mathfrak{F}\{\mathsf{Z}(t)\}\mathfrak{F}\{\mathsf{Q}(t)\}$$

We should have this relationship in the time domain as:

$$\mathsf{P}(t) = \int \mathsf{Z}(\tau) \mathsf{Q}(t-\tau) d\tau$$

Now I have this in Python:

import numpy as np

time = np.genfromtxt('time-data.txt',delimiter=',').T

fft = np.genfromtxt('fft-data.txt',delimiter=',',dtype=np.complex_).T

Z_t = time[0]
Q_t = time[1]

Z_fft = fft[0]
Q_fft = fft[1]
P_fft = fft[2]

#Shows that indeed Z_fft is the Fourier transform of Z_t
assert (np.fft.irfft(Z_fft,n=len(Q_t)) == Z_t).all()

#Shows that indeed Q_fft is the Fourier transform of Q_t
assert (np.fft.irfft(Q_fft,n=len(Q_t)) == Q_t).all()

#Shows that indeed Z_fft = P_fft / Q_fft
assert (np.divide(P_fft,Q_fft) == Z_fft).all()

P_t = np.fft.irfft(P_fft,n=len(Q_t))

P_t_from_convolve = np.convolve(Z_t,Q_t,mode='same')

#Expect to see that P_t = P_t_convolve but it's not the case!
assert (P_t == P_t_from_convolve).all()
```

Is there any reason that $\int \mathsf{Z}(\tau) \mathsf{Q}(t-\tau) d\tau \neq \mathfrak{F^{-1}}\{\mathfrak{F} \{ \mathsf{P}(t) \} \}$?

These are the data files here:

time-data.txt

fft-data.txt


问:
Norm estimates if adjoints can't be computed
Assume that you have two linear maps $A$ and $V$. For a given $x$ (of appropriate dimension) you can compute $Ax$ numerically, and for any $y$ (of appropriate dimension) you can calculate $V^Ty$ numerically. But you do not have access to both $A^Ty$ and $Vx$ on your computer. Hence, you can compute things like $$\langle Ax,y\rangle\quad\text{and}\quad \langle x,V^Ty\rangle$$
and you may even be able to estimate
$$\sup_{x,y}|\langle Ax,y\rangle -  \langle x,V^Ty\rangle|.$$

My question is:

Is there a way to reliably compute $\|A-V\|$? (Any norm would be fine, spectral norm desirable. An upper bound would also be ok.)

I'd like to add, that this situation is somehow always the case in floating point arithmetic, because each time you compute A*x or A'*y, different rounding errors occur, hence the results of (A*x)'*y and (x')*(A'*y) are different, too. (Btw, does anybody know an estimate for this quantity?) In double precision, this is probably not much of a problem, but in single or half precision, the errors seem to become non-neglectable.


问:
How to solve for underlying function from discrete data set containing integral of that function
New to Computational Science, I hope I'm on the right exchange network for this question.

I have a time series data set that contains the sum of a source data set representing an exponential decay curve at discrete points in time. In other words, the time series data set is the integral over discrete time bins. An example, written in R:

```<code>expFunc <- data.table(year=seq(1,20,1),
                      value=2*exp((-0.1*seq(1,20,1))+4)-2,
                      variable="Desired Exponential Function")
fiveYear <- data.table(year=c(1,5,10,15,20),
                       value=c(expFunc[1,value],
                               sum(expFunc[2:5,value]),
                               sum(expFunc[6:10,value]),
                               sum(expFunc[11:15,value]),
                               sum(expFunc[16:20,value])),
                       variable="Data That I Have")
</code>```

demonstrated in this plot

The data that I have follows the time steps in my example, where the first value is always the value of the desired function, and the remaining values represent the sum of the desired function over sequential time bins. Essentially, the data that I have I think can be represented as

$$
\int_{t_{1}}^{t_{2}} a e^{-bt+c}dt
$$

over different time intervals $t_{1} \rightarrow t_{2}$. I was imagining that I could solve for variables $a$, $b$, and $c$ using three equations, three unknowns but I'm either stuck in the algebra or am missing something.

Does anyone have an idea of how could obtain the underlying exponential function, either as

a discrete data set at a one year time resolution, or

an explicit function?

Thanks in advance for your ideas.


问:
Reference request: Textbook similar in structure to "Computational Materials Science: An Introduction" by June Lee for Computational Biology
In the Computational Materials Science: An Introduction by June Lee, he discusses molecular dynamics and density functional theory with examples from LAMMPS and QuantumEspresso, and explains LAMMPS and QuantumEspresso in detail.

I was wondering if there is a similar text in Computational Biology for molecular dynamics but in GROMACS/NAMD, I can use tutorials but I would like to find a textbook with more examples and explanations.


问:
Why does COMSOL treat a well as a rectangular cube?
I'm trying to learn COMSOL for a graduate research project and am struggling through a flow and transport in 3D tutorial. As such I'm walking through my boundary conditions and trying to figure out if I can find something that is wrong in my model.

I noticed that a cylinder, which I've assigned to be a well is treated as an elongate rectangular cube with all 12 edges, and 6 sides as selectable entities. I'm not sure if this is even the route of my problem, but my question is why treat the well as a rectangle instead of a cube?

Some background if it's helpful: my model consists of a simple 100mX100mX10m granite substance with the goal of modelling flow and transport using the well as a mass flow injection well, injecting at a rate of 1X10^6 kg/s. The well is 0.1m in diameter and 10m long.


问:
Hi I am trying to model a 2D Lug angle using Gmsh 4.6. How can I combine transfinite quad and regular full quad meshes in the following geo file?
I need transfinite mesh a small section of the bolt hole to insert a crack. However, The transfinite mesh and regular full quad mesh seem being incompatible and throwing errors. How can I combine transfinite quad and regular full quad meshes in the following geo file?
I used the following meshing options:

Packing of parallelogram (experimental)

Simple ful quad recombination

Recombine all triangular meshes (checked)

subdivision algorithm (none)
I require a transfinite mesh in Region 1 so that I can do J integrals.

```<code>// Gmsh project created on Wed Sep 02 17:32:13 2020
ri=19.05;
ro=3*ri;
w=2*ro;
b=200+ri;
lc1=20;
lc2=0.5;
lc3=2;

Point(1) = {0, 0, 0, lc1};
Point(2) = {b, 0, 0, lc1};
Point(3) = {b, w, 0, lc2};
Point(4) = {0, w, 0, lc1};
Point(5) = {b, ro, 0, lc1};

Line(1) = {1, 2};
//+
Line(2) = {1, 4};
//+

cx=b; // center of hole
cy=ro;
l1=cx-ri;
r1=cx+ri;

//+
Point(6) = {l1, cy, 0, lc3};
//+
Point(7) = {r1, cy, 0, lc3};
//+
ct=0.98480775301220805936674302458952;
st=0.17364817766693034885171662676931;
Point(8) = {cx-ri*st, cy+ri*ct, 0, lc2};
//+
Point(9) = {cx+ri*st, cy+ri*ct, 0, lc2};
//+
Circle(6) = {7, 5, 9};
Circle(7) = {9, 5, 8};
Circle(8) = {8, 5, 6};

Point(10) = {cx+ri*st, cy-ri*ct, 0, lc3};
Circle(5) = {6, 5, 10};
Circle(9) = {10, 5, 7};

/* Point(11) = {cx-ri*st, w, 0, lc2}; */
Point(11) = {cx-ro*st, w, 0, lc2};
Point(12) = {cx+ro*st, cy+ro*ct, 0, lc2};

Line(3) = {4, 11};
Circle(4) = {2, 5, 12};

Line(10) = {3, 11};
Circle(11) = {12, 5 ,3};

Line(12) = {8, 11};
Line(13) = {9, 12};

//+
Curve Loop(1) = {-12, -7, 13, 11, 10};
Plane Surface(1) = {1};

Transfinite Curve{7} = 15;
Transfinite Curve{10} = 8;
Transfinite Curve{11} = 8;

Transfinite Curve {12,13} = 60;

Transfinite Surface{1}={8, 9, 12, 11};
Recombine Surface{1};

Curve Loop(2) = {1, 4, -13, -6, -9, -5, -8, 12, -3, -2};
Plane Surface(2) = {2};

</code>```


问:
Is there any function to calculate condition number of sparse matrix in Eigen libraray?
The function JacobiSVD and BDCSVD can calcuate condtion number of a dense matrix via singular values.
However I need to know condition number of a sparese matrix due to slow computation speed using iterative linear solver. I can't find such function. If the sparse matrix is converted to one dense one in advance, the memory allocation occur because of large size of the matrix(about 10,000 x 10, 000). Many people may encounter this problem. Is there one way there to calculate condition number for a large sparese matrix in Eigen?</p>


问:
Split of complex parts in weak form
I am working on a numerical model to simulate the acoustic and elastic wave propagation in frequency domain via the Finite Element Method. Basically, the problem is to solve the Helmholtz equation in the fluid domain, the Navier-Cauchy equation in the solid domain and a couple of interface conditions.

More precisely, these interface conditions are the continuity of the normal component of the velocity and the continuity of traction, both on the interface $\Gamma$:

$$\cases{
\partial_{\mathbf{n}}p = \rho_F \omega^2 \mathbf{n} \cdot \mathbf{u}
\\
-p \mathbf{n} = \sigma(\mathbf{u}) \cdot \mathbf{n}
}$$

Here, $p$ is the deviation from ambient pressure, $\sigma$ the Cauchy stress tensor, $\rho_F$ fluid density, $\omega$ the angular frequency, $\mathbf{u}$ the displacement vector in $\Omega_S$ and $\mathbf{n}$ the unit normal vector of $\Gamma$. Their weak form are respectevely:

$$\int_{\Gamma} \rho_F \omega^2 (\mathbf{n} \cdot \mathbf{u}) q ds = 0
\quad (1)
$$
and
$$\int_{\Gamma} (p \mathbf{n}) \cdot \mathbf{v} ds = 0
\quad (2)
$$

where $q$ and $\mathbf{v}$ are test functions. Now I want to split $p$, $\mathbf{u}$, $q$ and $\mathbf{v}$ into its real and imaginary parts. I simply did $p = p^R + ip^I$, $\mathbf{u} = \mathbf{u}^R + i\mathbf{u}^I$, $q = q^R + iq^I$ and $\mathbf{v} = \mathbf{v}^R + i\mathbf{v}^I$. By apllying the distributive property, I rewrote $(1)$ and $(2)$, respectively, as

$$\int_{\Gamma} \rho_F \omega^2 [\mathbf{n} \cdot \mathbf{u}^R q^R - \mathbf{n} \cdot \mathbf{u}^I q^I + i(\mathbf{n} \cdot \mathbf{u}^R q^I + \mathbf{n} \cdot \mathbf{u}^I q^R)] ds = 0$$

and

$$\int_{\Gamma} [(p^R \mathbf{n}) \cdot \mathbf{v}^R - (p^I \mathbf{n}) \cdot \mathbf{v}^I + i((p^R \mathbf{n}) \cdot \mathbf{v}^I - (p^I \mathbf{n}) \cdot \mathbf{v}^R)] ds = 0$$

Well, since I have little experience on this topic, I am afraid it is not as simple as that. Can someone confirm it or correct me if I am wrong?

Title: Split of complex parts in weak form


问:
How can I practice multivariable root-finding?
Recently, I've been reading up on various root-finding / optimization algorithms such as the Levenberg-Marquardt method, Gauss-Newton, Conjugate Gradient, trust-region and trust-region-dogleg.

I've then been practicing using Matlab's fsolve function for some easy problems such as solving g (g(1), g(2)) = (0,0), where g's component functions are quadratic in two variables.  I can guess the roots easily, since I'm making up these practice problems.  Afterwards, I learn how to read all of fsolve's outputs and flags, and its optimization metrics and convergence analysis.

However, this is far too simple for my needs, and I want to solve a rootfinding problem in many more inputs and outputs.

Do you have a suggestion on how I can gain better practice to eventually solve my real problem in several input and output variables?

I need an intermediate form of practice but am struggling to come up with what to work on.

Any particularly good algorithm to look into would be appreciated  too.  It seems that Matlab's fsolve uses both the trust-region-dogleg and Levengberg-Marquardt algorithms.

Thanks,


问:
How to describe function convergence and function tolerance for numerical root-finding?
I'm currently doing some practice problems on root-finding and am writing up some notes / comments on my code.  In my solver code, if my function value is below the tolerance that I've set, should I write,

"the function converged to a root, since the function evaluation is below the function tolerance"

or

"the function converged to a root, since the function tolerance is achieved"?

Is there a better way to word it?


问:
FEM solution becoming wider as number of nodes increase
My FEM scheme uses a 4-node quadrilateral element with bilinear shape functions. The simple problem I'm solving is. $\nabla ^2 f = 5$ But as I increase the number of nodes, the plot of the solution becomes flatter. I'm stumped as to why this might happen.Two plots below show how when I change the size of the mesh from a 7x7 to a approx. 30x30, the plot becomes much more spread out. Also the boundary condition is that $f = 0$ on the boundary

FIXED, had some typos in my code


问:
Numerical Issues with DDE from SEIRU$\delta$ model
I'm new in this community. I moved this question from Math community. I'm reading the following article

Article Here

and my target is to replicate the results for a project.

SEIRU Model: I obtained very nice results though I not used article data but my data. Here is the evolution image I obtained playing around with the parameters.

SEIRU$\delta$ Model: Started with the same parameters values as the SEIRU Model. Cannot obtain a nice picture. The model has the same behaviour as the standard one, but I think there are some initial conditions errors. For small $\delta$ it is completely out of scale. To obtain the same numbers as the SEIRU I must use $\delta=81.5$ which is a very big number. Moreover with that $\delta$ I obtain a reproduction number $\mathcal{R}_0=2.6626e+05$ which again is impressive and not correct.

<img src="https://i.stack.imgur.com/0RV9i.jpg" alt="SEIRU$\delta$ Results" />

Never done numerical simulations with DDE. Maybe I did not manage very well the DDE system. I wrote my code in MATLAB, now I share it with you. Here is the links.

SEIRU Code Link

SEIRU$\delta$ Code Link

In the code there isn't the real data, only the fitted parameters the models need. Moreover the pictures you will see launching the code correspond to the picture here. I would say the target is to obtain the same results using the same fitted parameters. In this case they are $\chi_1$ and $\chi_2$, $\chi_3=1$ is fixed. Hope you can understand the code. Some easiest ways to manage DDE? Hope some of you could help me. Thanks in advance.


问:
Inverse problem with uncertain forward operator
Suppose I want to solve a linear inverse problem. In this example we take a convolution with the kernel:
$$\frac{1}{(y^2+z^2)^{3/2}}$$
We only take a fixed $z$ for the computation and convolve with respect to $y$. The associated matrix $A$ is dependent on $z$, denoted as $A_z$.

The output vector $b$ is (spatially) located at the same $y$ values, i.e. measurements done along $y$. Furthermore, we impose a nonnegativity constraint on the sources $x$.

We want to estimate the sources $x$ via minimizing the cost function
$$||A_zx-b||^2+\lambda\mathcal{R}(x) \,\,\,\,\,\,s.t.\,\,\,x\geq0$$
With $\mathcal{R}(x)$ being the regularizer.

My question: If I do not know $z$ exactly, am I somehow able to both correctly estimate $x$ and $z$ simultaneously? I could impose that the solution $x$ lies on some ball via the regularizer, such that $z$ could (maybe) be estimated correctly. But is this necessary, or can I somehow use less information than knowing the norm of $x$ beforehand?

I also think this may be dependent on the spatial frequencies present in the vector $b$. That is, if $b$ has high spatial frequencies with strong amplitude present, then we must (I think) have a minimum value for $z$, because too large $z$ can not result in $b$ having these high frequencies amplitudes.


问:
Multiplying by E[xy'] where only some statistics of xy' are known
(cross-posted on crossvalidated)

For random variable $(x,y)$ in $\mathbb{R}^{d}\times \mathbb{R}^{d}$ and vector $v \in \mathbb{R}^d$, I need to perform the following matrix vector multiplication.

$$T(v)=E[xy']v$$

The issue is that the expected matrix $E[xy']$ is too large to represent in memory ($d\approx $1 million), so I can only afford to store $O(d)$ worth of statistics and perform $O(d^{1.5})$ worth of operations. Three such statistics are $E[x]$, $E[y]$ and $E[x\odot y]$, where $\odot$ refers to element-wise multiplication.

If I only had the first two, one could argue that the following modification of $T$ is appropriate, representing an unbiased guess subject to these constraints

$$T(v)\approx E[x]E[y']v$$

What's an appropriate way to incorporate $E[x\odot y]$?


问:
Solution of non-linear Poisson equation does not match reference
I'm trying to solve the non-linear Poisson equation as a first step to solve the drift-diffusion equations for semiconductors. For reference, I'm using a preprint from the Weierstrass Institut (which I will refer to as [W]) that can be found here. All the model parameters and equations can be found there. I'll make references to specific parts of this paper.

I started implementing the algorithm in 2.4 in [W]. The discretization is done using FVM and the solution is then calculated using the Newton method. I ran into two issues:

The local electroneutrality in my case is the negative of what they get in [W]. What could be the cause of that? Do I have to switch the sign?

The non-linearity (_assemble_nonlinear) does not contribute correctly to the solution. Comparing my solution to the one in [W], Figure 2, it is clear that something went wrong. I'm suspecting, that the problem lies in the non-linearity since without it, the solution is a homogenous solution to the Poisson equation.

Here is my code with annotations refering to [W]:

import numpy as np
import sys
from typing import Callable
from scipy.constants import e as q, epsilon_0 as eps0
from scipy.constants import physical_constants
from scipy.linalg import solve_banded
from scipy.sparse import diags
import matplotlib.pyplot as plt
import matplotlib as mpl

np.set_printoptions(suppress=True,linewidth=np.nan,threshold=sys.maxsize)

kB = physical_constants["Boltzmann constant"][0]
kBev = physical_constants["Boltzmann constant in eV/K"][0]

# Setup for plotting
plt.style.use("ggplot")
plt.rcParams["text.usetex"] = True
plt.rcParams["font.family"] = "serif"
cmap = mpl.cm.get_cmap("summer")
plt.rcParams["axes.prop_cycle"] = mpl.cycler(
    color=[cmap(v) for v in np.linspace(0, 0.7, 10)]
)

class PoissonSolver:
    def __init__(self,
                 grid : np.ndarray,
                 C : Callable,
                 T : float,
                 Ev : float,
                 Ec : float,
                 Nv : float,
                 Nc : float,
                 epsr : float
        ):
        """
        grid : np.ndarray
            Grid for the calculation
        C : Callable
            Doping profile
        T : float
            Temperature in K
        Ev : float
            Valence band energy in eV
        Ec : float
            Conduction band energy in eV
        Nv : float
        Nc : float
        epsr : float
            Relative permettivity
        """
        self.N = len(grid)  # Number of grid points
        self.xks = grid
        self.xkk1s = np.pad(1 / 2 * (self.xks[1:] + self.xks[:-1]), (1, 1), constant_values=(self.xks[0], self.xks[-1]))
        self.hs = self.xks[1:] - self.xks[:-1]
        self.omegas = self.xkk1s[1:] - self.xkk1s[:-1]
        self.C = C(grid) # Doping profile
        self.T = T  # Temperature in K
        self.U_T =  kB * self.T / q  # Thermal voltage
        self.Ev = q * Ev  # Valence Band energy in J
        self.Ec = q * Ec  # Conduction Band energy in J
        self.Nv = Nv
        self.Nc = Nc
        self.epss = epsr * eps0  # Permettivity
        self.psis = []

    def _calc_psi0(self):
        # Intrisic carrier density squared from Eq. (9)
        Nisq = self.Nc * self.Nv * np.exp(-(self.Ec - self.Ev) / (kB * T))
        quad = np.zeros_like(self.C)
        nvals = np.where(self.C <= 0)
        pvals = np.where(self.C > 0)
        # Numerically stable quadratic formula used in Eq. (8)
        quad[nvals] = -2 * np.sqrt(Nisq) / (self.C[nvals] - np.sqrt(self.C[nvals] ** 2 + 4 * Nisq))
        quad[pvals] = (self.C[pvals] + np.sqrt(self.C[pvals] ** 2 + 4 * Nisq)) / (2 * np.sqrt(Nisq))
        # Actual Eq. (8)
        psi0 = - (self.Ec - self.Ev) / (2 * q) - 1 / 2 * self.U_T * np.log(self.Nc / self.Nv) + self.U_T * np.log(quad)

        # Plot potential
        #fig = plt.figure()
        #ax = fig.add_subplot(111)
        #ax.plot(self.xks, psi0)
        #ax.set_xlabel("Position [$\mu m$]")
        #ax.set_ylabel("Initial Potential [V]")
        #plt.show()

        return psi0.reshape((self.N, 1))

    def _assemble_matrix(self):
        M = np.zeros((3, self.N))  # Only store diagonals
        K = np.zeros((self.N, self.N))  # Store entire matrix
        for i in range(self.N):
            # M is Jacobian of Eq. (12)
            # K is the matrix for the linear part of Eq. (12)
            if i == 0:
                M[1, i] = 1
                M[2, i] = - self.epss / self.hs[i]
            elif i == self.N - 1:
                M[0, i] = - self.epss / self.hs[i - 1]
                M[1, i] = 1
            else:
                # Derivatives of the mobile carrier concentration
                En = q * self.psis[-1][i][0]
                dpdpsi = -Nv / self.U_T * np.exp((self.Ev - En) / (kB * T))
                dndpsi =  Nc / self.U_T * np.exp((En - self.Ec) / (kB * T))

                M[1, i] = + self.epss * (1 / self.hs[i] + 1 / self.hs[i - 1]) + q * self.omegas[i] * (dpdpsi - dndpsi)
                if i != 1:
                    M[0, i] = - self.epss / self.hs[i]
                if i != self.N - 2:
                    M[2, i] = - self.epss / self.hs[i - 1]

                K[i, i - 1] = - self.epss / self.hs[i - 1]
                K[i, i] = + self.epss * (1 / self.hs[i] + 1 / self.hs[i - 1])
                K[i, i + 1] = - self.epss / self.hs[i]

        return M, K

    def _assemble_nonlinear(self):
        # Non-linear part of Eq.(12)
        S = np.zeros((self.N, 1))
        for i in range(1, self.N - 1):
            # Mobile carrier concentration
            En = q * self.psis[-1][i][0]
            p = Nv * np.exp((self.Ev - En) / (kB * T))
            n = Nc * np.exp((En - self.Ec) / (kB * T))

            S[i] = -q * (self.C[i] + p - n) * self.omegas[i]
        plt.show()

        return S

    def solve(self, relax=1.0):
        self.psis = [self._calc_psi0()]
        for n in range(9):
            # Newton iterations
            dFdpsi, K = self._assemble_matrix()
            S = self._assemble_nonlinear()
            F = np.dot(K, self.psis[-1]) + S
            dpsi = solve_banded((1, 1), dFdpsi, -relax * F)
            psi = self.psis[-1] + dpsi
            self.psis.append(psi)

        self.psis = np.array(self.psis)
        for psi in self.psis:
            plt.plot(self.xks, psi)
        plt.show()

# Same parameters as in Fig. 2
T = 300
Ev = 0
Ec = 1.424
Nv = 9e18
Nc = 4.7e17
epsr = 12.9
NA = 1e17
ND = 1e16
grid = np.linspace(0, 10e-6, 51)

def C(x):
    return -ND * np.heaviside(-(x - 1e-6), 1) + NA * np.heaviside((x - 9e-6), 1)

ps = PoissonSolver(grid, C=C, T=T, Ev=Ev, Ec=Ec, Nv=Nv, Nc=Nc, epsr=epsr)
#fig = plt.figure()
#ax = fig.add_subplot(111)
#ax.plot(ps.xks, ps.C)
#ax.set_xlabel("Position [$\mu m$]")
#ax.set_ylabel("Doping concentration [idk]")
#plt.show()
ps.solve(relax=1.0)
```

Thank you in advance for your help!


问:
Which library to use in C++ for fast DCT's?
I need to apply a 1D type-II DCT (Discrete Cosine Transform) to each column of a matrix while normalizing the result in C++. Effectively, I want to multiply my matrix with an orthogonal type-II DCT-matrix (while avoiding the performance cost of matrix multiplication, of course).

Since this seems like an important computational kernel in many applications, isn't there a highly optimized library for this purpose out there? After searching for a bit I can find a couple of implementations, but I don't know if they're good. I was expecting a function somewhere in a well-known library such as Boost or Eigen. After all, in Python everything I need to do is this (using SciPy): dct(x, norm='ortho')


问:
Remez algorithm convergence
I have implemented the Remez algorithm in Python where all calculations were done with the Python mpmath library. I have noticed that sometimes the $|E_{max}|$ and $|E_{min}|$ do not monotonically increase anymore but recurrently iterate over the same values. Mainly, my stopping condition $|E_{max}| - |E_{min}| < u$ where $u$ is round-off error of the working precision never occurs. Increasing the working precision of mpmath did not solve this behavior.

For all steps of the Remez algorithm, I use the default methods for obtaining solutions for each step: mpmath's LU decomposition for solving a linear system of equations, the secant method for root finding and golden section search for extrema.
Is using different methods for the different steps, possibly more precise, will improve this behavior, or is there something else that I miss?
Is my stopping criterion is practical or I must relax it for general use?

Title: Remez algorithm convergence


问:
Finite elements algorith for a fluid in a tube with an elastic obstacle (Fluid-Solid coupling)
I want to solve the model of a tube with an elastic obstacle, something like a simple model of an vessel with a valve. The fluid is given by an evolutionary incompressible Navier--Stokes equations, and the obstacle for elasticity equation.

In the step 46 of dealii appears a simple way to manage this kind of problems

If I understand correctly, the algorithm proposed is the following:

Solve the fluid problem assuming that on the elastic obstacle the velocity is equal to zero.

Use the velocity fluid solution on the boundary of the solid as a boundary condition in the elasticity problem, and then solve the elastic problem.

Update the fluid-solid mesh.

I'm understanding correctly? In that page says that is only an academic algorithm, with the goal of explaining the implementation, and is not a realistic FSI model.

I have two questions:

This algorithm looks easy to program because it uses existing solvers. Do you know if it is an algorithm used for FSI? (a one-way algorithm perhaps).

and my main question:

2. Which is the Finite Elements algorithm for this kind of coupled problems? I searched a lot on internet, because surely someone already proposed an algorithm, but I could not find an algorithm.


问:
Collaborative coding platform in the style of Google Colab
(I hope this question is fitting for this SE, I couldn't find very relevant tags)

I am looking for an online platform in the style of Google Colab with the following requirements:

possibility to run notebooks like in JupyterLab

necessity of at least python and julia

and following wishes:

privacy-focused (even if I am not coding for sensitive projects I don't like the idea that my code is hosted at Google)

possibility to have a team of 5 (at least) students and coworkers (preferably 10) and the ability to share and version the code

possibility to pay for more computational power

I am considering free and paid options.


问:
Save font selection for plots in Paraview
I adjusted plotting fonts in a Paraview scene, saved as .pvsm file, and loaded it in a different Paraview version. The fonts changed. Can this be avoided? Also ticks (I set it to maximum of 3 ticks) settings are gone. Is .pvsm not the correct way to save a complex scene?


问:
Numerically estimate the Jacobian from a scalar time-series
I'm trying to numerically estimate the Jacobian from a time-series. Following the paper, https://journals.aps.org/pra/abstract/10.1103/PhysRevA.43.2787.

Consider that I have a scalar time series $x = (x(1), x(2), \dots, x(N - 1), x(N))$. In total there are $N$ data-points. So the steps that I'm following are as follows.

Compute the time-delay (call it $\tau$).

Compute the embedding dimension (call it $d_{emb}$).

Compute the delay vectors (call it $\bar y_{i}(k)$). $\bar y_{i}(k) = [x(i), x(i + \tau), x(i + 2\tau), x(i + 3\tau), \dots, x(i + (d_{emb} - 1) \tau]$.

Choose some random initial point in the time series and call the delay vector corresponding to that point, $\bar y(k)$.

Look at the neighbours of $\bar y(k)$ and store them in ascending order of their Euclidean distance in $\bar y^{(r)}(k)$. Say we keep $N_{b}$ neighbours.

Store the distance between the delay vector and its neighbours in $\bar z^{(r)}(k) = \bar y^{(r)}(k) - \bar y(k)$.

The basic idea is that the underlying dynamics is governed by a map $\bar F$ which takes $\bar y(k) \rightarrow \bar y(k + 1)$.

Compute the next evolution of $\bar y(k)$. Therefore, $\bar y(k+1) = \bar F(\bar y(k))$. Similarly do it for the $r$th neighbour of $\bar y(k)$. One gets, $\bar y(r; k+1) = \bar F(\bar y^{(r)}(k))$.

Compute the distance between these vectors, $\bar z(r; k+1) = \bar y(r; k+1) - \bar y(k+1)$.

Now, $\bar z(r; k+1) = \bar y(r; k+1) - \bar y(k + 1) = \bar F(y^{(r)}(k)) - \bar F(y(k)) = \bar F(\bar y(k) + \bar z^{(r)}(k)) - \bar F(\bar y(k))= M_{\alpha \beta} z^{(r)}_{\beta}(k) + Q_{\alpha \beta \gamma} z^{(r)}_{\beta}(k) z^{(r)}_{\gamma}(k) + \cdots$.

We assume that $\bar z^{(r)}(k)$ are small so that one can do a Taylor expansion about it. $M$ is the Jacobian matrix of dimension, $((d_{emb} - 1)\tau, (d_{emb} - 1)\tau)$. $Q$ is the Hessian matrix.

Idea now is to compute the unknown coefficients of $M$ using least-square fit keeping $H$ (for accuracy and dropping it finally).

The authors suggest that if $\tau$ be the step-size for the next iteration, so $\bar y(k) \rightarrow \bar y(k + 1)$, then $\bar z(r; k + 1) = M_{\alpha \beta} z^{(r)}(k)$. And then the only components unknown are the bottom row of $M$ Rest are $1$'s along the superdiag, and rest are all $0$.

What I don't understand it is that $\bar z^{(r)}(k) \rightarrow \bar F(\bar z^{(r)}(k)) = \bar F(\bar y^{(r)}(k) - \bar y(k)) \ne \bar F(\bar y^{(r)}(k)) - \bar F(\bar y (k)) = z(r; k+1)$. So how do they write -- $\bar z(r; k+1) = M_{\alpha \beta} z^{(r)}_{\beta}(k)$?


问:
Error curve does not oscillate in between reference points using Remez
Using the Remez algorithm, implemented using multi-precision library, in certain functions that I want to approximate, the error curve does not oscillate in between reference points, and so no roots can be found (subsequently no extrema).

For example, this happens when I try to approximate $ln$ (natural logarithm) over the interval of $I=[a=0.5, b=1]$, with a rational polynomial $P/Q$ of degrees $n=deg(P)=6$ and $m=deg(Q)=5$.

For the first set of reference points, Chebyshev nodes are computed using the following:
$$
x_i^{(0)} = 0.5 \cdot (a + b) + 0.5 \cdot (b - a) * \cos(\frac{(2i-1)\cdot \pi}{2n})
$$
$$
i = 1,\ldots,N=n+m+2
$$

The figure below shows the first iteration of the error function with the first set of nodes marked in circles.

Next, the extrema are found:

The extrema, together with the reference points of interval bound, i.e., $x_a$ and $x_b$ compose the next set of reference points ${x^{(1)}}$

The following figure shows the error function in the second step, with the reference points marked with circles.

Here, the 7th and the 8th reference points are both below the x-axis and so no root will be found in between them. Regardless, this error curve has 10 extrema. Together with the two bounding reference points, there will be at most 12 new reference points which is less than 13, the minimum required for this approximation. So it seems that the error curve, in this step, does not oscillate enough. How this can be fixed so that more oscillation occurs?


问:
Ill-conditioned stiffness matrix
I am writting a Fem code in c++ for a 2d plane stress model. My question is regarding the assembly stiffness matrix.I noticed that some elements of the matrix are not exactly zero but insted a number close to zero for example  -1.7677e-15. Can floating point arithmetic erros result to an ill-conditioned stiffness matrix? If yes how can I avoid it?</p>

Title: Ill-conditioned stiffness matrix


问:
Open-source iterative solvers robust to noise?
I need to solve $Ax=b$ in about 1 million dimensions. Furthermore, A is only accessible through matrix vector products, and these are noisy/inexact.

Is there any solver with Python interface I can try on this problem?

Basically I want to check if there's anything in the linalg community that might work to train a neural network better than current method. The current method is SGD which is equivalent to linearizing the problem and applying first-order Richardson iteration.


问:
Scaling tensor approximation by symmetric tensor decomposition with SciPy's L-BFGS-B
I am trying to approximate a symmetric tensor of which the values are in the range of [1e-7,1e-4], by a symmetric tensor decomposition of lower rank. For this I am using the L-BFGS-B method in SciPy's optimization package.

The relative errors (based on the Frobenius norm) I am obtaining are quite big (greater than 1). After some research I have come to the conclusion that there is need for scaling, as my tensor is 'poorly scaled'. When I multiply all values of X with 1e7, I do obtain a smaller relative error (in the range of [1e-4,1e-3]), but only when the modes of X have a small dimension in comparison with the chosen rank of the symmetric tensor decomposition.

Since I am not that familiar with the concept of scaling in numerical optimization, is there a better way to solve this scaling problem, so that I can obtain a small relative error even when the dimension of the modes of X is large in comparison with the chosen rank?

I am also doing some research about the existence of a low rank approximation, as this may not even exist. While this may be the case, reproducing the first experiment of Sherman & Kolda (2020), does not give me the same order of magnitude of the relative errors. This means that there should exist some improvements to be made to my implementation, of which I think only the scaling aspect can be improved.


问:
2nd-order TVD criteria for flux-limiter
Consider a nonlinear hyperbolic conservation equation:
$$
\partial_{t}u = -\partial_{x}f(u)
$$
The spatial derivative of $f(u)$ may approximated after a spatial discretization by $x_{j}=j\Delta x$
$$
\partial_{x}f(u_{j}) \approx \frac{1}{\Delta x}\left(F_{j+\frac{1}{2}}-F_{j-\frac{1}{2}}\right)
$$
where $F$ represents a numerical flux function that may depend on neighboring points and holds the same form from one cell to the next.  In order to avoid Gibbs phenomena, the idea is to combine a low-order flux and high-order (2nd or greater) flux and combine them with a limiter function $\phi$
$$
F = F_{L} - \phi(r)(F_{H}-F_{L})
$$
where
$$
r_{i} = \frac{u_{i+1}-u_{i}}{u_{i}-u_{i-1}}
$$
The goal is preserve monotonicity by keeping the total variation (TV) of the solution non-increasing.
I understand how the application of the TVD condition leads to the condition on the limiter:
$$
\begin{cases}
0 \leq \phi(r_{i}) \leq 2r & r<1 \\
0 \leq \phi(r_{i}) \leq 2 & 1 \leq r
\end{cases}
$$

But I do not understand at all the conditions required to keep the scheme 2nd-order and TVD.  Could someone please show me how that region (see image) is derived?

Title: 2nd-order TVD criteria for flux-limiter


问:
Most scalable distributed consensus mechanism based on message complexity?
One of the most challenges in distributed consensus mechanisms is both time complexity and message complexity.

For example, PBFT message complexity is O(n^2) that means that it is only scalable to tens of nodes. Thus, classical BFT (Byzantine Fault Tolerant) consensus mechanisms are not scalable to the large networks.

Nevertheless, I am looking for the most scalable BFT consensus  among all.

In other words, is there a distributed BFT consensus with less than O(n^2) message complexity?

P.S. If you think that there a better and more relevant Stack Exchange site to bring up this question, please let me know.


问:
Convergence of Conjugate Gradient Algorithm
I am trying to solve a linear elasticity model using finite element discretization in a rectangle domain [0,1]x[0,1]. For the solution of the the linear system $Ku=F$ I am using the CG algorithm. However, I have noticed some weird behaviour of the algorithm. For small grids (e.g 2x2 or 3x3 or 2x1 ...) the algorithm fails and doesn't converge while if the grid becomes larger(10x10 or 25x25 ...) the algorithm converges and the results seem reasonable.My first thought was If the application of homogeneous BC leads to an ill-conditioned matrix ,since I zero the columns and rows accosiated with BC and put 1 in diagonals, but this make no sense if the algorithm converges for some grids right? Any idea why is this happening?


问:
Implementation method selection for sparse constrained linear least squares or quadratic programming
I need to slove one optimization problem of quadratic programming. The number of optimization variables is about 16,000. The constraints include equality constraints and inequality constraints.

I have no such practical experiences before.
There are are three choices for the probelm after reading some materials:

```<code> active set method,
 interior point method
 augmented Lagrangian method
</code>```

I need implement the optimization algorithm on my own.

Active set method is not suitable for such problem size.

Interior point method is fast but difficult to implement from the link:

rank-deficient NNLS

A first order method (Augmented Lagrangian, ADMM, split Bregman, etc.) These are possible to implement yourself without needing to use a packaged library.

So augmented Lagrangian method will be my choice.
How about my analysis?
Most materials I find is about augmented Lagrangian method with equality constraints. Can you recommand any links or materials on augmented Lagrangian method with inequality constraints?


问:
Haw to apply central difference to viscous flux in energy equation?
In many modern papers Navier-Stokes equations are solved with finite-difference or finite-volume methods using WENO reconstruction for non-viscous fluxes and central differences for viscous ones. It seems pretty straightforward to apply central difference to "simple" second derivative like $\frac{\partial^2 u}{\partial x^2}$, but I can't find detailed examples for more complex ones like $\frac{\partial}{\partial x} \left( u\frac{\partial u}{\partial x}\right)$.

Here are some details:

Suppose we solve 1D transient compressible Navier-Stokes equations with constant viscosity and heat conductivity coefficients:
$$
\frac{\partial U}{\partial t} +
\frac{\partial F}{\partial x} = 0; \qquad
U = \left( \begin{array} {c} \rho \\ \rho u \\ e  \end{array} \right), \;\;
F = \left( \begin{array} {c} \rho u \\ \rho u^2 + \sigma \\ (e + \sigma)u - k \frac{\partial T}{\partial x}\end{array} \right)
$$
where $\sigma$ is viscous stress:
$$
\sigma = p - \frac{4}{3}\mu \frac{\partial u}{\partial x}
$$

If we implement explicit Maccormack method according to his 1975 paper, we compute the following predictor and corrector steps:
$$
U^{(p)}_i = U^n_i - \frac{\Delta t}{\Delta x} \left( F_i - F_{i-1}\right) \\
U^{n+1}_i = \frac{1}{2} \left[ U^n_i + U^{(p)}_i - \frac{\Delta t}{\Delta x} \left( F^{(p)}_{i+1} - F^{(p)}_{i}\right) \right]
$$
and $\frac{\partial u}{\partial x}$ terms inside $F$ and $F^{(p)}$ vectors are simply substituted with forward and backward differences:
$$
\frac{u_{i+1} - u_{i}}{\Delta x},\;\; \frac{u^{(p)}_{i} - u^{(p)}_{i-1}}{\Delta x}
$$
for predictor and corrector steps, correspondingly. Choosing difference direction opposite to used in predictor and corrector steps effectively generates 2nd-order central difference analogue.
Pretty straightforward and seems to work well in my tests.

Suppose now we switch to scheme with WENO flux reconstruction and 2nd or 3rd-order Runge-Kutta time stepping. Then terms like $\rho u$, $\rho u^2 + p$ and $(e+p)u$  in vector $F$ may be reconstructed with WENO. Terms like $\frac{\partial^2 u}{\partial x^2}$ or $\frac{\partial^2 T}{\partial x^2}$ extracted from $\frac{\partial F}{\partial x}$ may be approximated with central difference:
$$
\frac{u_{i+1} - 2u_{i} + u_{i-1}}{\Delta x^2}.
$$
and treated as source terms in Runge-Kutta steps.

But what to do with $\frac{\partial}{\partial x} \left( u\frac{\partial u}{\partial x}\right)$ term? Should we use something like
$$
\frac{\frac{u_{i+1} + u_{i}}{2}\frac{u_{i+1} - u_{i}}{\Delta x} - \frac{u_{i} + u_{i-1}}{2}\frac{u_{i} - u_{i-1}}{\Delta x}   }{\Delta x}?
$$

And what to do if we want to use higher-order central difference?


问:
Negative binomial expansion of general symbolic polynomial
Using Sympy, I would like to compute the negative binomial expansion of a general symbolic polynomial, e.g., $(x_1 + x_2 + x_3 + 4 x_4)^{-1}$. I understand that I can go by recursively partitioning the expression inside the parenthesis into two parts and computing their expansion like $x_1$ and $x_2 + x_3 + 4 x_4$, then, $x_2$ and $x_3 + 4 x_4$, and so on.

Is there a more efficient way to do this? Thanks.


问:
Time complexity of derivation, gradient,differential, jacobian matrix
what is the time complexity of gradient $\nabla_{f}$ using the $\mathcal O$-notation?
what is the time complexity of jacobian matrix using the $\mathcal O$-notation?
who knows some references to introduce the time complexity of gradient?
I make reply to reviewers for computational complexity of gradient and matrix operations. I need some references to cite.


问:
How to use FEniCS to calculate the electric field of an isolated charged sphere
Initially I thought that this is the kind of question which ought to have already been answered in the form of an example online, but so far I haven't found one. I will admit that I am very new to both FEniCS and the use of FEM software for solving PDEs in general, so my goal here in solving this problem is to a gain a better understanding of the boundary conditions (B.Cs.) used in FEM software (specifically FEniCS). Based on that I suppose my question can be broken down into the following parts:

For isolated electrostatic systems (basically an object/objects with a non-zero charge) do we need to somehow include the B.C. that $v(r) \rightarrow 0$ as $r \rightarrow \infty$ ? If not, what B.Cs. do we choose to make sure that the FEM solution approximates the analytical solution?

If we do need to pass the information that the potential goes to zero as $r$ goes to infinity how do we do so on a finite domain? What type of B.Cs. are we supposed to use and how do we implement them?

As far as answering this question I would appreciate anything from useful links to books/websites to explanations with equations. Ideally, I would like it if someone could show me how to implement this particular problem in FEniCS but I don't know how practical that is. Thank you for your time.

Edit:
Is it possible to use something like the Boundary Element Method (BEM) to solve this type of problem more exactly? I found a software package that is written in Python that connects to FEniCS called Bempp (Link to Bempp website). A basic description of the BEM comes Wikipedia (Link to BEM Wikipedia page).


问:
Bipartite Euclidean Matching simple to implement approximate algorithm
I am looking for a simple to implement algorithm for the bipartite euclidean matching problem (or an implementation of any practical algorithm). I am aware of Agarwal's paper, but I would like to implement something simpler. The simplest idea I could come up with was starting with some initial matching (possibly random) and then minimizing the energy through simulated annealing. In my case the vertices of each class are ~ $2^{14} - 2^{20}$ so brute-forcing the solution is not realistic. A benefit of the stochastic method is that I also do not necessarily need to compute the possibly $2^{20} \times 2^{20}$ distance matrix, since I can compute the distances for each swap on demand.
Any suggestions/alternatives are welcome.


问:
Dolfin convert : How to interpolate data at vertices of (3D) cells?
I hope that one of you guys can help me because i have been stuck here for a week.
I am trying to read a gmsh file (.msh) using dolfin convert to XML and then download it with dolfin.

The thing is when i assemble the stiffness matrix and the mass matrix and try to generate an exemple (MonteCarlo) using these matrices , i have figured out that they are not the good ones .

I think that dolfin when downloaded the XML file which has already convert it from MSH , he changes the ordering of the numbering of the elements in ascending way.
I am sure that he does this, because i have tried my code on a box generated by BoxMesh of fenics and my program worked and if i check mesh.cells() and my XML file i can see that the ordering had been changed, So that’s why i am not having my good matrices

And to be honest, i don’t know how to fix this thing? Does dolfin destroys the ordering of any mesh and build another ordering (in an ascending way) ? how do you guys get the good stiffness matrix and mass one when you download a mesh of your type?

I have read the same problem founded by another person but still cannot understand how to fix that.

Thank you very much in advance

this is my code

```<code>from dolfin import *

import os

import numpy as np

from numpy import linalg

import scipy

from scipy.sparse import csr_matrix

from scipy.sparse.linalg import spsolve

import scipy.stats as stats

from scipy.optimize import minimize

alpha=1.22

beta_x=2

beta_y=0.01

beta_z=0.01

mu=-4.25

str_os = ‘dolfin-convert ./mesh/calcul_little_specimen.msh
./calcul_little_specimen.xml’

os.system(str_os)

mesh = Mesh(“calcul_little_specimen.xml”)

V = FunctionSpace(mesh, “Lagrange”, 1)

u = TrialFunction(V)

v = TestFunction(V)

parameters[‘linear_algebra_backend’] = ‘Eigen’

S=as_matrix([[beta_x**2,0,0],[0,0,0],[0,0,0]])

kx=inner(dot(S,grad(u)),grad(v))*dx

Kx=assemble(kx)

row,col,val = as_backend_type(Kx).data()

Kx_sp=scipy.sparse.csc_matrix((val,col,row))

S=as_matrix([[0,0,0],[0,beta_y**2,0],[0,0,0]])

ky=inner(dot(S,grad(u)),grad(v))*dx

Ky=assemble(ky)

row,col,val = as_backend_type(Ky).data()

Ky_sp=scipy.sparse.csc_matrix((val,col,row))

S=as_matrix([[0,0,0],[0,0,0],[0,0,beta_z**2]])

kz=inner(dot(S,grad(u)),grad(v))*dx

Kz=assemble(kz)

row,col,val = as_backend_type(Kz).data()

Kz_sp=scipy.sparse.csc_matrix((val,col,row))

m=uvdx

M = assemble(m)

row,col,val = as_backend_type(M).data()

M_sp=scipy.sparse.csc_matrix((val,col,row))

from sksparse.cholmod import cholesky

factor=cholesky(M_sp,ordering_method=“natural”)

L=factor.L()

W=np.random.randn(M_sp.get_shape()[0],1)

y=L.dot(W)

#generation of an example

e=spsolve(M_sp+Kx_sp+Ky_sp+Kz_sp,alpha*y)

nu=np.ones((M_sp.get_shape()[0],1))*mu

e=e.reshape(-1,1)

e=e+nu

file_mean_kalman = File(‘results/realisation.pvd’)

func = Function(V)

func.vector().set_local(e)

func.rename(‘Kal_mean’,‘Kal_mean’)

file_mean_kalman << func
</code>```


问:
Solution predictors for accelerating convergence in nonlinear FEM
I am looking for the details of commonly-used predictors for accelerating the convergence of iterations using Newton-Raphson scheme for nonlinear problems in FEM. I am looking specifically for static problems.

Note: Such predictors based on velocity and/or acceleration are used for transient problems but they are not applicable for pure static problems.

Appreciate any inputs.