问:
Rotational kinematics problem @ $\theta = 0$
I'm simulating a magnetic dipole that is subjected to an evolving magnetic field.  In ISO convention ($\theta$ is the polar angle and $\phi$ is the azimuthal angle), my equation of motion that is giving me issues is

$$\ddot\phi = \frac{\textrm{garbage}}{\sin{\theta}}.$$

So this equation makes perfect sense, as $\theta$ goes to one of the two poles, $0$ or $\pi$, then the azimuthal angle has to quickly change (if $\theta$ went directly through the pole then $\phi$ would change almost instantly to be $\phi_0 \pm \pi$)  My question is not understanding what is going on; rather, how is the best way to approach this problem since $\ddot\phi$ blows up and the simulation stops at the poles?

So far I have tried moving my coordinate system so the poles were in a place where possibly the trajectory wouldn't pass through it then convert my coordinates back (this worked but occasionally I do get a trajectory close and it blows up).  I've tried decreasing step size, which makes it better, but still if you get too close it blows up.  I've tried looking on-line since I'm sure I'm not the first to come across this problem, but was unable to find anything.  Any help would be great!


问:
Rank filter on an nXm array using python
I would like to apply a rank filter on an nXm numpy array. Let's say I have this array:

```<code>a = np.random.randint(100,size=(5,5))

array([[ 9, 71, 55, 72, 43],
       [46, 56, 44, 17, 64],
       [53, 51, 17, 71, 36],
       [12, 28, 97, 46, 35],
       [29, 55, 27, 48,  5]])
</code>```

I would like to run a rank filter (kernel window of 3X3) that will give each central pixel in the kernel a value\rank (from 1-9 because it is a 3X3 window) according its position in the histogram. The largest number will get the value of 9 and the lowest number will get the value of 1. For example, the value 56 in position 1,1 will get the value 8 because according to its 8 neighbors, it is the second largest number. The value 44 in position 1,2 will get the value of 3 because it is the 3rd lowest number, there are only two values below it both are 17. (let's ignore for now the edges issues and the tied values issues).

I found that scipy.stats.rankdata can achieve my desired results but the problem is that it is not a convolved filter, meaning, that it works on the entire array and not according to a kernel method, also it gives an output of a 1Xm array, meaning that it does not keep the original nXm array shape.

So I think I need somehow to use this scipy.stats.rankdata filter as a moving kernel of a size of my choise.

I know that scipy has convolve methods such as ndimage.convolve, which I used in the past, and this method works as a running kernel window. But I can't figure out how to combine them in order to get a "proper" kernel ranking filter.

I also saw that scipy has scipy.ndimage.filters.rank_filter but I don't understand how the ranking is actually work there.

ContentLicense: CC BY-SA 3.0


问:
Are the No Free Lunch Theorems Useful for Anything?
I have been thinking about the No Free Lunch (NFL) theorems lately, and I have a question which probably every one who has ever thought of the NFL theorems has also had. I am asking this question here, because I have not found a good discussion of it anywhere else.

The NFL theorems are very interesting theoretical results which do not hold in most practical circumstances, because a key assumption of the NFL theorems is rather strong. This assumption is, roughly speaking, that the performance of an algorithm is averaged over all problem instances drawn from a uniform probability distribution. In realistic applications the problems an algorithm typically encounters are NOT drawn from a uniform distribution, and are instead drawn from what is likely a very interesting and complicated distribution specific to the general problem setting.

So, while the NFL theorems are quite interesting results, do they have any practical implications? Or are they merely theoretical results?


问:
Finding "hidden" subassemblies: Suggestions for algorithms
I'm faced with an interesting problem and need some help finding existing best techniques to solve it. The setup is that we're analyzing a large system trying to find what I call "hidden subassemblies." Using round numbers, there are 1,000 different basic parts that are used in 10,000 different combinations or "assemblies". (I call a unique full combination an "assembly.") One assembly might have 150 parts, another 700.

What we're trying to do is efficiently detect "hidden subassemblies." Namely, groups of parts that frequently show up together. I'm sure there's an existing body of research and practice, algorithms and statistical methods for exactly this problem...but I don't know what it is, can't invent it myself (sad but true), and don't know the terminology to look it up. Can anyone point me in the right direction? This seems like a problem that would be found in manufacturing, anything related to groups of people, and genetics.

In case I haven't been clear enough, what I'm trying to find are groups of parts that are commonly (or always) used together. Say that there's a #3 bolt and it's always used along with a #3 nut. That's what I mean by a "hidden subassembly." In practice, we're actually likely to find much larger hidden subassemblies, but they may not be 100% the same across all assemblies.

I'd be incredibly grateful for suggestions. While I'm no good at inventing mathematical solutions, I can usually bang them out (eventually) if pointed in the right direction.

Thanks very much!


问:
Use of multidimensional FFTW and normalisation factor
I am using the FFTW MPI in C.

I have a simple question.

Quoting from  fftw.org

The multi-dimensional transforms of FFTW, in general, compute simply the separable product of the given 1d transform along each dimension of the array. Since each of these transforms is unnormalized, computing the forward followed by the backward/inverse multi-dimensional transform will result in the original array scaled by the product of the normalization factors for each dimension (e.g. the product of the dimension sizes, for a multi-dimensional DFT).

Imagine we have $\tilde{A}$ in the Fourier space with the coresponding $A$ in the physical space. I want to compute $\tilde{B}$ that is $\tilde{A}^2$. The normalization factor is $F_{norm}$.

The idea is $A=IFFT(\tilde{A})$ (inverse fourier transform, A now is real)

So

$B=A^2$

$A=A*F_{norm}$, $B=B*F_{norm}$ (we have to normalize  before to transform back to Fourier space.)

$\tilde{A}=FFT(A)$

$\tilde{B}=FFT(B)$

The questions are:

1) Do I normalize A before or after the computation of B ?

2) Does anything change if B is used directly in physical space ?


问:
Solving linear elasticity equations for a cantilever problem
Has anyone successfully used the linear elasticity equations to model a cantilever problem?

The equation I am solving is simply:
  $$\nabla \cdot \boldsymbol{\sigma} = 0$$

subject to boundary conditions at the fixed end, free surfaces, as well as where the load is applied. I apply the load as a shear stress/traction boundary condition on the end face.

I am having no luck getting my results to match the analytical solution. I plotted the lateral displacement along the length of the beam and the shape of my curve is approximately linear, whereas the analytical solution is a cubic. When comparing max displacement in my numerical result vs the analytical result, I am off by orders of magnitude.

I've studied my code and numerics carefully and I do not see any mistakes, but I could be wrong. Is the cantilever problem able to be modeled with the equation above?

Editing this post based on the comments: In addition to the above equation, I am also using Hooke's law and the stress strain displacement relationship:

$$\boldsymbol{\sigma} = 2\mu \boldsymbol{\epsilon} + \lambda tr(\boldsymbol{\epsilon}) I$$
$$\boldsymbol{\epsilon} =0.5[\boldsymbol{\nabla u} + \boldsymbol{(\nabla u)}^T]$$

Editing the post to include the analytical solution and dimensions.
The analytical solution that I am using is (http://ruina.mae.cornell.edu/Courses/ME4735-2012/Rand4770Vibrations/BeamFormulas.pdf):
$$v = \frac{Fx^2}{6EI}(3L-x)$$
Where v is the lateral displacement or the y-displacement in the link. The beam I am currently modeling has L = 10m
with a square cross section of 1m x 1m. So the moment of inertia, $I=\frac{1}{12}$ for this problem.
I apply a force of $10^6$N at the free end, which is equivalent to a $10^6$Pa shear stress at the free end. I prescribe this as a traction boundary condition so that the boundary condition at the free surface is:
$$\sigma = \begin{bmatrix}
0 & -10^6 & 0 \\
-10^6 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}$$

or in terms of tractions:
$$
T = \sigma n = \begin{bmatrix}
0 \\
-10^6 \\
0
\end{bmatrix}
$$

Attached is the result I am getting for lateral displacement for $E=2*10^{11} Pa$.


问:
Huge accelerations in plasma simulation
I'm trying to make a numerical simulation of pulsar magnetosphere using FDTD on a log-spherical Yee lattice for fields and PIC for plasma particles. Field part is working like charm but issue arises when adding the plasma.

To simulate the dynamics of charged plasma particles I use Newton-Lorenz equation:

$\frac{d\gamma\vec{v}}{dt} = \frac{e}{m}\left( \vec{E} + \frac{\vec{v}}{c} \times \vec{B} \right)$

Adopting CGS system of units, characteristic values for $\vec{E}$ and $\vec{B}$ are of order of magnitule $10^{12}$, charge to mass ratio for electron is of order of magnitude $10^{17}$, so the value for acceleration is ~$10^{29}$.

In the code itself I normalized all variables in the following manner:

$
\vec{u} = \frac{\gamma\vec{v}}{c}\\
\vec{E}^{*} = \vec{E}/B_0\\
\vec{B}^{*} = \vec{B}/B_0\\
\tau = \frac{ct}{R}\\
$

Where $B_0=10^{12}$, $R=10^6$ (radius of the pulsar) and $c$ is the speed of light. As a result, my equation looks as follows:

$\frac{d\vec{u}}{dt} = \frac{eRB_0}{mc^2}\left(\vec{E}^{*}+\frac{\vec{u}\times\vec{B}^{*}}{\gamma}\right)$

Normalization factor $\frac{eRB_0}{mc^2}\sim 10^{15}$, acceleration is still huge.

I may choose different normalization for time as follows:

$t = t_0\tau\\
\frac{eB_0}{mc}t_0 = 1$

But with $\tau$ chosen in that way it takes literally forever for fields to get to stationary configuration.

What can I do?

OwnerDisplayName: Tajimura

Title: Huge accelerations in plasma simulation


问:
Stabilizing online average calculation
In Knuth, the following method for computing an average is presented:
\begin{align*}
M_{n} = M_{n-1} + (x_{n} - M_{n-1})/n
\end{align*}
(See here, if you don't have TAOCP.)

Assuming the samples all lie a bounded distance away from the average, then as $n \to \infty$, $(M_{n-1}-x_n)/n$ becomes very small, and eventually $M_{n-1} \oplus (x_{n}\ominus M_{n-1}) \oslash n$ becomes bitwise equal to $M_{n-1}$. The average computed by this formula then stops evolving in response to new data.

What are the known techniques to avoid this?

Update: I believe that this scheme has error $|\hat{M}_n -M_{n}| \le n\mu_{M} \underline{M}_{n}$, where
\begin{align*}
\underline{M}_{n} := \frac{1}{n} \sum_{i=0}^{n-1} |x_{i}|,
\end{align*}
and $\mu_{M}$ is the unit roundoff, but I couldn't prove it. This won't directly answer the question, but might put me on the track to get improvements. Any ideas?

Title: Stabilizing online average calculation


问:
How to account for the interface between two different phases in a discretized diffusion model?
I have tried to set up a model for the diffusion of a gas into a liquid. The two media are next to each other and the geometry is spherical because the system should simulate the diffusion out of a gas bubble.
I have made a diffusion model for the gas phase (the left compartment) and a corresponding diffusion model for the liquid phase (the right compartment).
The two diffusion models are as follows:
\begin{align}
  \frac{\partial c_\mathrm{left}}{\partial t}
  &= D_\mathrm{left} \cdot \frac{1}{r^2} \frac{\partial}{\partial r}
     \left(r^2 \frac{\partial c_{left}}{\partial r} \right)\\
  \frac{\partial c_\mathrm{right}}{\partial t}
  &= D_\mathrm{right} \cdot \frac{1}{r^2} \frac{\partial}{\partial r}
  \left(r^2 \frac{\partial c_{right}}{\partial r} \right)\\
\end{align}

Then I have discretized the volume. So the PDEs are turned into ODEs of the form:
\begin{align}
  \frac{\mathrm{d}c_{i,\mathrm{left}}}{\mathrm{d}t}
  &= D \left(\frac{c_{i+1,\mathrm{left}}
     - 2 \cdot c_i + c_{i-1,\mathrm{left}}}{\Delta r^2}
     + \frac{2}{r_{i,\mathrm{left}}}
     \cdot \frac{c_{i+1,\mathrm{left}}
     - c_{i-1,left}}{2 \cdot \Delta r} \right)\\
  \frac{\mathrm{d}c_{i,\mathrm{right}}}{\mathrm{d}t}
  &= D \left(\frac{c_{i+1,\mathrm{right}}
     - 2 \cdot c_i + c_{i-1,\mathrm{right}}}{\Delta r^2}
     + \frac{2}{r_{i,\mathrm{right}}}
     \cdot \frac{c_{i+1,\mathrm{right}}
     - c_{i-1,\mathrm{right}}}{2 \cdot \Delta r} \right)
\end{align}

When discretizing I got the following expressions:
\begin{align}
  \frac{\mathrm{d}^2c_i}{\mathrm{d}r^2}
  &= \frac{c_{i+1} - 2 \cdot c_i + c_{i-1}}{\Delta r^2}\\
  \frac{\mathrm{d}c_i}{\mathrm{d}r}
  &= \frac{c_{i+1} - c_{i-1}}{2 \cdot \Delta r}
\end{align}

But how does one incorporate the following two boundary conditions at the gas-liquid-interface, where $H$ is Henry's constant?
\begin{align}
  c_\mathrm{left,interface} &= H \cdot c_\mathrm{right,interface}\tag1\\
  J_\mathrm{left,interface} &= J_\mathrm{right,interface}\tag2
\end{align}

Meaning that there is equal flux across the interface between the two compartments.

The concentration in the last compartment of the gas phase (before the interface) is approximated by linear extrapolation.
This means that the concentration just before the interface is given by:
$$c_{3} = \frac{c2-c1}{r2-r1} \cdot r_3 + c_1 - r_1 \frac{c_2-c_1}{r_2 - r_1}$$

In this example it is assumed that only three compartments are present on each side of the interface.


问:
Discretization used for steady linear elasticity - followup to previous question
This is a follow-up of my previous scicomp question (https://scicomp.stackexchange.com/posts/28863/edit). I figured I'd start a new thread on this as the question is a bit different from my previous question.

I am solving the linear elasticity equations, using Hooke's Law and the strain-displacement relationship.

My discretization is as follows using a finite volume formulation:

$$\nabla \cdot \boldsymbol{\sigma} = \nabla \cdot (\mu \boldsymbol{\nabla u} + \mu(\boldsymbol{\nabla u})^T + \lambda\boldsymbol{I}tr(\boldsymbol{\nabla u})) = 0$$

In integral form, we have:

$$\iint_A [\mu \boldsymbol{\nabla u} + \mu(\boldsymbol{\nabla u})^T + \lambda\boldsymbol{I}tr(\boldsymbol{\nabla u})] \cdot dA = 0$$

I evaluate this surface integral using:
$$\iint_A [\mu \boldsymbol{\nabla u} + \mu(\boldsymbol{\nabla u})^T + \lambda\boldsymbol{I}tr(\boldsymbol{\nabla u})] \cdot dA = \sum\limits_{k=1}^N |A_k|n_k\cdot[\mu (\boldsymbol{\nabla u})_k + \mu(\boldsymbol{\nabla u})_k^T + \lambda\boldsymbol{I}tr(\boldsymbol{\nabla u})_k] $$

I set the above equation to $f(\boldsymbol{u})$ for simplicity. I can expand this term using a Taylor series expansion of the form:
$$f(\boldsymbol{u}) = f(0) + \frac{df}{d\boldsymbol{u}}\boldsymbol{u} = f(0) + \boldsymbol{Ju}$$

I compute the jacobian, $\boldsymbol{J}$, using a forward differencing scheme:

$$J = \frac{d\boldsymbol{f}}{d\boldsymbol{u}} = \frac{f(\boldsymbol{u}+h)-f(\boldsymbol{u})}{h}$$

So this gives us a linear system of the form:
$$\boldsymbol{Ju} = \boldsymbol{b} = -f(0)$$

Since this is a finite volume formulation, my mesh consists of cells. I compute the cell gradients using the discrete green gauss method:

$$\boldsymbol{\nabla} u = -\frac{1}{V}\sum\limits_{j=1}^N \hat{u}_jA_j\boldsymbol{n_j}$$

Where $\hat{u}=\frac{u_L+u_R}{2}$, i.e., the face displacement is obtained by averaging the neighboring cells' displacement.

I compute face gradients using a simple average and applying a correction term:
$$\nabla u_{face} = \frac{1}{2}(\nabla u_L + \nabla  u_R) - \underbrace{\hat{\boldsymbol{d_{lr}}} [\frac{1}{2}(\nabla u_L + \nabla u_R)\cdot\hat{\boldsymbol{d_{lr}}} ] + (u_L - u_R)\frac{\hat{\boldsymbol{d_{lr}}}}{|\boldsymbol{d_{lr}}|}}_{\text{Correction term}}$$

So this is the gist of my formulation.

I tried to model a cantilever problem, and it ended in a disaster (see previous question). I have tested out this formulation on the classical 1-D compression problem and was able to match analytical results.

Does anyone see anything wrong with this formulation?


问:
Non-reflecting boundary conditions for compressible Navier-Stokes equations
I have some questions about the implementation of non-reflecting OUTFLOW boundary condition for Navier Stokes equations.
Following

Poinsot, Lele "Boundary Conditions for Direct Simulations of Compressible Viscous Flows"

Pirozzoli, Colonius "Generalized characteristic relaxation boundary conditions for unsteady compressible flow simulations"

both authors suggest an equation like

\begin{equation}
\frac{d\mathbf{u}_b}{dt} + \mathbf{d(u_b)} + \mathcal{\mathbf{T(u_b)}} = \mathbf{S(u_b)}\quad\quad (1)
\end{equation}

where $\mathbf{u}_b = (\rho, \rho u, \rho v, \rho w, \rho e)$ are the conservative variable at the bound b; $\mathbf{d(u_b)}$ is a certain characteristic treatment of x-flux; $\mathcal{\mathbf{T(u_b)}}$ are the transverse therms; $\mathbf{S(u_b)}$ are the source therm.

Now consider a finite difference method and suppose that the bound is located at $n+1/2$ (as usual in this problems); considering furthermore $n+1, n+2,\dots, n+GN$ ghost points.
Equation (1) can be applied just in $n+1/2$, so my questions are:

How can I discretise the derivatives in $n+1/2$ with some FD scheme in order to solve (1) at this location.

How can I treat the ghost nodes?


问:
Stable Method of orthogonal projection onto a subspace with the help of Moore-Penrose inverse,
Projection of a vector $v$ onto the column space of a matrix $A$ is given by $AA^\dagger v$. From the definition of Moore-Penrose Inverse we know that $AA^\dagger v = (A^T)^\dagger A^T v $.

Below is the code for implementing the projection of a random vector onto the space of a random matrix. This is related to my other question Backward stable algorithm to get orthogonal projection onto the column space of a matrix. I didn't get an answer to that question.

I would like to know why is there a huge difference between the two methods of calculating the projection.

```<code>  % testingprojfrostackexchange
  clear;
  M = 1400;
  N = 1300;
  r = 1;
  A = rand(M,N);
  u = rand(M,r);

  projLN = pinv(A')*(A'*u);%This is projection through Least Norm
  projLS = A*(pinv(A)*u);%This is projection through Least Square

  [Q R] = qr(A);
  Q = Q(:,1:N);
  z1 = Q*(Q'*u);%This is the actual projection

  display('(projection through QQT) - projLS');
  norm(z1-projLN)/norm(projLN)

  display('(projection through QQT) - projLN');
  norm(z1-projLS)/norm(projLS)
</code>```

Output

```<code>>> stackexchange
   (projection through QQT) - projLS

   ans =

   2.1569e-13

   (projection through QQT) - projLN

   ans =

   8.3546e-15
</code>```

The actual projection which is given by $Q(:,1:N)Q(:,1:N)^T$, where $Q$ is from the Housholder decomposition of $A$. We get the projections from $AA^\dagger$ and $(A^T)^\dagger A^T$ and find that $AA^\dagger$ is much better when $A$ is a tall matrix.


问:
overlapping additive schwarz
I am solving 1D laplace problem discretized with finite differences (3-point stencil).
I would like to use additive Schwarz method in classical form:

$U_{k+1}=U_{k}+M^{−1} r_k,$

where $r_k=F−A U_k$ and $  M^{−1} =  \sum_{i=1}^N R_i^T(R_i A R_i^T)^{−1} R_i$.

Restriction operators $R_i$ are defined as binary matrices, with ones marking existing parts of the domain.

Method converges nicely for non-overlapping case (block Jacobi). However, I am not able to obtain convergence in case of any overlap.

If I use restricted additive Schwarz, method is convergent again.
What is reason for such a behaviour?

Title: overlapping additive schwarz


问:
Kronecker product of matrices
I have to use Kronecker product of a matrix with a unit matrix is there any routine in ScaLapack or Lapack which can do so efficiently.</p>

Title: Kronecker product of matrices


问:
Radially symmetric system of PDEs in deal.II
I am trying to solve the radially symmetric polar form of the PDE with homogeneous Neumann BC in deal.II on a unit circle:
$$ u_t = \Delta u - \nabla \cdot (u \nabla h) $$
$$ h_t = \Delta h $$

I am not sure how to handle the origin & how the bilinear form change in this case.

Also, how to use the grid and triangulation in deal.II.

Could someone please explain or suggest similar literature?


问:
What is the inverse Laplace transform algorithm that is most accurate given the fewest frequencies considered?
Based on your empirical knowledge.

This paper suggests a nonlinearly accelerated Fourier series approach, such as the one proposed here, but I have one constraint: we should be able to express the method as a linear combination of the function $F$ at different frequencies $s$, that is

$f(t) \approx \sum\limits_{i=1}^N w_i^t \, F(s_i^t)$

In this equation I am looking for the most accurate algorithm that requires the lowest $N$.


问:
Procedure to identify characteristic properties of unknown functions in a DAE model
I have a system of 1st order odes given by

$$
\dot{x_1}(t) = \alpha_1 f_1(x_1,t) + \beta_1 u(t) \\
\dot{x_2}(t) = \alpha_2 f_2(x_2,t) + \beta_2 u(t)
$$

They are constrained by an algebraic equation

$$ x_1(t) + x_2(t) =  k $$

where $\left( \alpha_1,\alpha_2, \beta_1,\beta_2 , k \right) \in \mathbb{R}$  are known constants (i.e. parameters). $f_1(t)$ and $f_2(t)$ are both unknown.

Starting from a rich set of input-output noise-free data available from simulating a complex proxy system, what would be the best procedure to identify (even a subset of repeatable/characteristic properties) the unknown possibly time-varying functions $f_1(x_1,t)$ and $f_2(x_2,t)?$ I am almost certain that $f_1(x_1,t)$ and $f_2(x_2,t)$ are both linear.

I am looking for a suitable approach that shall work well to arbitrary excitations in all future simulations. (NOT merely a curve-fitting procedure to match a specific input-output dataset)


问:
Levenberg-Marquardt for root-finding: just square the function?
This question might be so obvious and trivial that I'm having a hard time googling it.

I have a multivariate root finding problem that I'm trying to solve in C# and the library that I'm trying to use is alglib.  That package doesn't have any dedicated root solvers in it, but has Levenberg-Marquardt.  To do root finding, do I just take the squared vector result and feed it to the LM solver?  Then simply check the final result against the zero vector to ensure that it found the global minimum I'm looking for?


问:
Alternatives to Newton-Raphson for nonlinear elasticity via finite element
As far as I have seen, solving problems of nonlinear elasticity using the finite element method proceeds by linearizing, either around the initial configuration (total Lagrangian approach) or around the most recently solved configuration (updated Lagrangian approach).  This produces a system of linear equations, and represents an example of the Newton-Raphson method.

My question is about possible generalizations to other methods of solving systems of nonlinear equations.  I realize that for engineering problems, the linearization may not cause any problems, because the goal is often to study deformations starting from a known, unstressed configuration, which can then be gradually loaded.  In other words, there is a good small parameter, so linearization is fine.  However, I'm interested in the case that external forces (loadings) depend on the deformation, and the dependence can exhibit strong nonlinearities.  In this case, the global solution may not be anywhere near the undeformed configuration.

I think that attempts to address this problem are known as "globalization strategies", and I'm interested in the application of these to nonlinear elasticity.  I know about methods like arc-length control, but I'm interested in other methods such as full multigrid.

I don't know enough about the lingo or the literature to even search the literature properly - so far my naive attempts have not yielded anything.  I also don't know good examples of problems with strong nonlinearities that might be solved with finite element - maybe Navier-Stokes?  I don't even see how to properly formulate the problem as a system of nonlinear equations - do I just numerically integrate the principal of virtual work?


问:
Can I expect reasonable results when expanding in scaled Legendre polynomials?
Imagine I want to compute the eigenvalues of an operator $\hat O$ defined on $ L^2(\mathbb{R})$, however using a properly scaled N-dimensional polynomial basis of $ L^2([-a,a]) $ which fulfills Dirichlet boundary conditions: $P_n(\pm a)=0$ (e.g. scaled Legendre polynomials), i.e. $ O_{i,j}=\langle P_i|\hat O|P_j\rangle$.

Can I expect convergence to the exact result for $N, a\to\infty$?


问:
Conservation in finite element codes
Typical finite volume methods are conservative, because fluxes (of e.g., mass or energy) are always between neighboring cells. Is the same generally true for finite element codes? Do I correctly assume that implicit time stepping would then be 'globally' conservative, except for flows through the domain boundaries?

Background information: The reason I'm asking is that colleagues who use the plasma module of Comsol do not observe charge conservation in their simulations. The equations are of the following form:
$$\partial_t n_1 + \nabla \cdot \vec{f} = S \, n_1\\
\partial_t n_2 = S \, n_1,$$
where $\vec{f}$ and $S$ non-linearly depend on other variables. The quantity $\int n_2 - n_1 dV$ is not conserved, but the conservation error shrinks for smaller time steps. I believe that this module internally works with $log(n)$ instead of $n$, where $n$ is the density, to avoid certain numerical problems.

Title: Conservation in finite element codes


问:
Monotone, monotonicity preserving, LED, TVD, NVD, bounded, stable and stability preserving discretization schemes
When it comes to discretization schemes for finite volume method, the following terms can be found in literature:

monotone schemes

monotonicity preserving schemes

local extremum diminishing schemes

TVD schemes

NVD schemes

bounded schemes

stable schemes

stability preserving schemes

A single work doesn't use all of them, rendering it more difficult to understand how these terms compare.

Could you please explain the relationship between the terms? For example: are some of them synonyms, subclasses of others, equivalent provided that some other condition is met?


问:
Meaning of a symbol in a research paper
The 2014 paper "Iteration-Free Computation of Gauss--Legendre Quadrature Nodes and Weights" by I. Bogaert (https://doi.org/10.1137/140954969) contains the following expression in Appendix A:

What is the meaning of the "LL" here? It is not defined in the body of the paper and it's such a commonly occurring character sequence that Googling for it has proven futile.

Title: Meaning of a symbol in a research paper


问:
Genetic Algorithm: Need some clarification on selection and what to do when crossover doesn't happen
I'm writing a genetic algorithm to minimize a function. I have two questions, one in regards to selection and the other with regards to crossover and what to do when it doesn't happen.

Here's an outline of what I'm doing:

```<code>while (number of new population < current population)
    # Evaluate all fitnesses and give them a rank. Choose individual based on rank (wheel roulette) to get first parent.
    # Do it again to get second parent, ensuring parent1 =/= parent2

    # Elitism (do only once): choose the fittest individual and immediately copy to new generation

    Multi-point crossover: 50% chance
    if (crossover happened)
        do single point mutation on child (0.75%)
    else
        pick random individual to be copied into new population.
end
</code>```

And all of this is under another while loop which tracks fitness progression and number of iterations, which I didn't include. So, my questions:

As you can see, two parents are chosen randomly in each
iteration until the new population is filled up. So, the two same
parents may mate more than once and surely several fit parents will
mate many more times than once. Is this in any way bad?

In the obitko tutorial, it says if crossover doesn't
happen, then child is exact copy of parents. I don't even understand
what that means, so, as you can see, I just picked a random parent
(uniformly; no fitness considered) and copied to new population.
This seems weird to me. Whether I actually do this or not, my results really don't change that much. What's the proper way to handle the case
when crossover doesn't happen?


问:
Left eigenvectors using ARPACK
I'm trying to find both the dominant $k$ left and right eigenvectors, that is,

$$V_L\mathcal{A} = \Lambda V_L\\
\mathcal{A}V_R = V_R\Lambda\\
V_LV_R = I_{k\times k}$$

 $V_L$ being the $k\times N$ matrix of left eigenvectors,

 $V_R$ being the $N\times k$ matrix of right eigenvectors, and

 $\Lambda$ being the diagonal $k\times k$ matrix of eigenvalues.

I'm using ARPACK (via scipy.sparse.linalg). scipy.sparse.eigs gives me  the right eigensystem $V_R$ and $\Lambda$; I've tried two methods to get $V_L$:

Apply scipy.sparse.linalg.eigs to $\mathcal{A}^T$ to get left eigenvectors and eigenvalues (call them $V_{L,0}$ and $\Lambda_L$)

If $V_{L,0}V_R$ is singular or $\Lambda_L$ is too different from $\Lambda$, reduce $k$; if that doesn't work, restart

Set $V_L \leftarrow (V_{L,0}V_R)^{-1}V_{L,0}$; or

 In principle this is reasonably fast (2x the running time of just getting the right eigenvectors), but in practice, in the cases where it needs to be restarted, it may need to re-run many times.

 For each of the $k$ eigenvalues $\lambda_i$, apply scipy.sparse.linalg.lsmr to the overdetermined system
$$\lambda_i^{-1}v_L^{(i)}\mathcal{A} = v_L,\\
v_L^{(i)}v_R^{(j)} = \delta_{i,j}$$
 This is much more reliable than the other method, but it has a much worse best-case running time and it also has a quite limited accuracy.

I imagine it should in principle be possible to get left and right eigenvectors via IRAM; is this available in some way on ARPACK? (As mentioned I've been using scipy.sparse.linalg, but I'm amenable to accessing ARPACK more directly as a last resort). Is there some way to repackage this approximate eigendecomposition as something that ARPACK can handle (sort of like the $AA^\dagger$ trick for SVD)? Or is there some other way to get $V_L$, after having gotten $V_R$ and $\Lambda$, that's reasonably efficient, reliable, and accurate?

On the matrices I'm trying to diagonalize:

 Complex, non-Hermitian (also not symmetric, unitary, normal, etc.)

 Represented via scipy.sparse.linalg.LinearOperator, not technically a sparse matrix

 Equally efficient to perform left- and right- vector multiplication

 Magnitude of eigenvalues should theoretically decay fast (between polynomially and exponentially); initial eigenvalues tend to decay even faster

 For an example see arXiv:1510.00689; operator is defined in Fig. 2 and eigenvalue problem is described in Fig. A.4.b.

Title: Left eigenvectors using ARPACK


问:
Magnetostatic modelling Radia: Increasing distance between magnets a produces positive force until a certain point, beyond which force goes haywire
I have two sets of magnets. One set consists of two electromagnets (Shown below: Blue) and the other set consists of two NdFeB N40 permanent magnets. Both sets of magnets lie on the same plane. I want to plot the force of attraction/ repulsion between the two sets of magnets as the distance between them varies from 10cm to 100cm (Y-axis below is in mm - and shows the initial condition of 10cm (100mm) separation).

I have modeled this setup using Radia Free software found here. You must have Mathematica (I am currently using version 11.2) to use Radia.

My problem: When Radia solves the model and produces an output of Force vs. Distance I expect ALL the force measurements to be EITHER positive OR negative - meaning either attraction or repulsion between the two sets of magnets. What actually happens though is after about 50cm separation Radia produces a mixture of positive and negative force results. The First plot looks accurate (yellow and blue lines) but on closer inspection some of the values are negative and this becomes glaringly obvious when the Force axis (Y) is made logarithmic.

And the logarithmic version:

Included in the graph (orange) is the solution I obtained (Force vs. Distance) from a 2D magnetostatics software called FEMM. This is purely a reference (which I believe is what the yellow and blue lines should look like especially in the logarithmic version) showing a Force>0 for all distances from 10cm to 100cm.

I don't know where I am going wrong and would really appreciate some help. I have read on the Radia documentation that correct subdivision is crucial - especially to iron components. I have tried to subdivide my iron core using radObjDivMag[] Radia Functions but I can't get it to work (keep in mind I used the same function on a bar magnet (cuboid) and the function worked perfectly, I have also tried to use the function argument {"cyl",....} but it will still not subdivide my component, if you know where I am going wrong please comment).

I would like to attach my Mathematica notebook (heavily commented however if you have any questions please ask me) which I used to create and solve the model and the Excel spreadsheet used to plot my data but don't seem to be able to - Please message me and I can send it to you.

Has anyone used Radia before? Should I be using different software?


问:
Can post-processing use topoSet and createPatch without screwing up the results?
I would like to use patchAverage to obtain the average pressure on an object that consists of multiple patches (due to different boundary conditions), however patchAverage only averages over single patches. As long as the constituent patches aren't shared with other objects, I could simply take the area-weighted average of each patch's average.

However, due to a dynamic geometry where the individual patches change area in time, this would become quite messy. So instead, my idea is to use topoSet to define a faceSet and turn that into a new object patch via createPatch (replacing the old ones, which are no longer needed after simulating).

Would createPatch correctly map the simulated fields to the new patch or would I end up with values being mapped to the wrong vertices? Or is there another way to achieve this? (I also considered using sample, which can obtain the values for multiple surfaces at once, but then I'd have to manually average, and I'm not sure this way I'd weigh the patch faces correctly.)


问:
fitting exponential versus exponential w/ power
I have two models which I would like to investigate for my data.  One form is:
\begin{equation}
\label{one}
    f(r) = A e^{-B r}
\end{equation}
and the second is:
\begin{equation}
\label{two}
    g(r) = \frac{A}{r^{C}} e^{-B r}
\end{equation}

Does anyone have any advice on how to tell the difference between these two?  I typically have about ~ seven points that I can fit in all my situations.  The problem is these seven points don't span a particularly large range of $r$ and so it's hard to see the difference between $g$ and $f$.  The least-squares fits in the case of $g$ are more unstable, but I assume it's because it's very similar to $f$ and the parameters have more freedom...

Any advise or wisdom would be welcome, and I can clarify if there are any questions.


问:
How does the MADS algorithm work in practice
Mesh Adaptive Direct Search (MASH) is an algorithm for black box optimization

I want to understand an implement this method to solve some 2D multivariate blackbox function $f(x,y)$, but am having trouble understanding the algorithm.

What I understand so far:

Start with some initial trial value $x_0, y_0$

Create the 'Mesh'

Search a finite number of points on the mesh

If the search fails, perform a 'Poll' step $P_k = x_k + \Delta_k d$

Where I am confused

How do we 'create' the Mesh? How does the mesh change with each iteration?

How do we perform a search step? How do we choose which points on the mesh to search?

What does it mean to perform a poll step? How should $\Delta_k$ change? What is meant by $d$?

Thanks in advance for an help


问:
Open source implementation of Multiscale Combinatorial Grouping
I would like to use Multiscale Combinatorial Grouping for my PhD research. However, I am restricted to use open-source implementations and this one runs on Matlab.

Does anyone know of an equivalent contour detection method that works on large images that is implemented open-source, preferably in Python?

Here some more background info. Equivalent methods would be:

Contour detection and hierarchical image segmentation -> cannot handle images of more than 1000x1000 pixels and runs on Matlab

SLIC Superpixels -> creates too many redundant segments

I am interested in detecting contours of visible objects in drone imagery with a resolution of 2-5cm per pixel and extents of up to 1kmx1km for the purpose of automating cadastral mapping. Here an example of what the contour detection output should look like:


问:
Stokes flow around rigid body
I'm trying to simulate Stokes flow in 2D around an arbitrary polygon (representing a rigid body).  I'd like to get both the effect of the body on the flow velocity and the forces on the body by the fluid.  Basically, like how panel methods are used to examine inviscid flow around airfoils, I'm trying to examine "airfoils" in a Stokes flow.

If I'm reading it right, in The Method of Regularized Stokeslets, section 2.2, it suggests doing this by discretizing the surface of the body and solving for the "external forcing" given the velocities of the fluid at the surface of the body.

That is, the velocity of a Stokes flow at $\vec{x}$ given point forces $\vec{f_k}$ acting at $\vec{x_k}$ is:

$$u(x) = \sum_{k=1}^N \frac{-\vec{f_k}}{4 \pi \mu} \ln|r| + \vec{f_k} \cdot (\vec{x} - \vec{x_k}) \frac{\vec{x} - \vec{x_k}}{4 \pi \mu |r|^2}$$

With $r = |\vec{x} - \vec{x_k}|$.  The paper suggests inverting that to find the forces given the velocities.  However, I don't entirely understand what is meant by "external forcing" in the case of flow around a body. Is it the same as the fluid's force on the body (equal but opposite)? If so, it doesn't seem to match the pressure forces, given that the paper says the pressure is given by

$$p(x) = \sum_{k=1}^N \frac{\vec{f_k} \cdot (\vec{x} - \vec{x_k})}{2 \pi (\vec{x} - \vec{x_k})^2}$$

Am I misunderstanding something?

Title: Stokes flow around rigid body


问:
Why would someone use empirical sum instead of numerical integration methods?
In the context of a scientific computing application, using data coming from (powerful) embedded systems, acquiring raw data (but from calibrated acquisition electronics), I have been asked to integrate these data (in RT), with an empirical sum, instead of making use of usual numerical methods. (Important note : even if we are dealing with embedded systems, computational power is not an issue).

When I was asking why not using these usual numerical methods (e.g. Simpson's, midpoint rule, trapezoidal rule, etc.), I could not get any answer. Then, I am wondering why someone would want to perform it this way.

The acquisition application on the embedded system is in C++03, compiled with GCC 4.4.7, the average being performed looking somewhat like this (this code is executed each time the acquisition is performed) :

```<code>float sum = 0.0;
for (int i = 0; i < DATA_BLOCK_LENGTH; i++) {
    sum += data[i]
}
sum /= DATA_BLOCK_LENGTH;
</code>```

This snippet gives you an idea of the code performing the average from a trace of the signal received ; the signal usually looking like this most of the time (one or two spikes along the whole signal trace).

The previous constant, DATA_BLOCK_LENGTH, is usually set to the maximum (2048, which is the size of electronics buffer). The measures are equally-spaced samples.


问:
Finite differenced eigenvalue prob. of inhomogeneous boundary conditions?
I am basically asking about eigenvalue problems of differential equations using some finite difference method (FDM). Usually the system is subject to some boundary conditions (BC), e.g., Dirichlet or Neumann. Often it is homogeneous BC like $u(0)=0,u'(R)=0$, which seems to be fine to deal with. One way is to introduce and then eliminate ghost points with the BC.

But if it is inhomogeneous BC (e.g., for a Dirichlet $u(0)=a,u(R)=b$), the nonzero constants have a sting in the tail. Effectively they're put into an inhomogeneous column as shown below, which is also clearly introduced in this note. For instance, the eigenvalue probelm $Lu=\lambda u$ of the operator $L=-\frac{\partial^2}{{\partial x}^2}$ becomes an inhomogeneous eigenvalue problem $Mu+R=\lambda u$ with
$$M=\frac{1}{h^2}\begin{bmatrix}
2 &-1 & & & \\
-1 &2 &-1 & & \\
&\ddots  &\ddots &\ddots & \\
& & -1 &2 &-1\\
 & & & -1 &2
\end{bmatrix},
R=\frac{1}{h^2}\begin{bmatrix}
a \\
0 \\
\vdots \\
0\\
 b
\end{bmatrix}.$$

But after this, practically how to proceed to solve the original problem? Any other formulation? Or we simply don't need to think about such problem for some reason?

This may not be a problem for usual DE solving (I naively assume one-sided stencil often suffices.). But for eigenvalue case, I can't find any way out. The inhomogeneous eigenvalue problem doesn't have any standard or general solution as far as I've searched.


问:
Padua-type pointset for functions singular on line $x=y$
The Padua points $\mathrm{Pad}_{n} \subset [-1,1]^{2}$ are a unisolvent pointset with optimal growth of Lebesgue constant, described in detail here. With some work they can be used to generate a multivariate Chebyshev series representation of the function they attempt to interpolate.

I would like to use the Padua points to generate a multivariate Chebyshev series of a function singular along the line $x=y$, but this is impossible because two of the Padua points are consecutive vertices of the square.

Is there a way to generate a multivariate Chebyshev polynomial interpolant of a function $f\colon [-1,1]^2\to \mathbb{R}$ which is singular on the line $x=y$, perhaps using a modification of the Padua points?

Update: The goal is to write a piece of code which projects as many 2D functions defined on rectangular regions as possible onto Chebyshev polynomials and allow for interpolation and recovery of the coefficients for perhaps other uses. This is similar to the goals of chebfun in 2D.


问:
Numerical analysis: Chebyshev coefficient representation error
If $x_k$ are the Chebyshev nodes, that is for $n \in \mathbb{N}^*$, we have $x_k = \cos(\pi \frac{2k + 1}{n})$. Now suppose you have approximated values of $x_k$ and $x_{k'}$ for $k,k' \leq n$. In the worst case, what is the number of significant values we may loose when we calculate $x_{k'} - x_k$?

Now, taking into account the fact that each real is represented as $x \approx m \times b^e$, we then have a formula that  for $x'$ the approximated value of $x$ and $y'$ approximated value of $y$, $\Delta (x' + y') \leq \epsilon(\Delta| x| + \Delta |y| + |x| + |y|)$ where $\epsilon = b^{1-N}$, where $N$ is there number of significant numbers after the comma.

So for a fixed $N \in \mathbb{N}^*$ we have $\Delta (x_{k'}' - x_k') \leq \epsilon (\Delta |x_{k'}| + \Delta |x_k| +  |x_{k'}| + |x_k| ) \leq 4 \epsilon $

Thus in worst case we have a loss of precision of $\frac{4}{b^{N-1}}$. Thus we lose one significant term?


问:
Global truncation error behavior at fixed time step
I am trying to solve the following diffusion equation problem:

$\frac{\partial f}{\partial t}=\frac{\partial (D\frac{\partial f}{\partial x})}{\partial x}+S$

$D=1+x^{2}+\sin(x)$

$f(x,0)=1  ,  f(0,t)=f(1,t)=1$

$S $ is choosen adequately in order to have $f_{ana} = 1+\sin(\pi x)\sin(\pi t)$ as the analytical solution of the equation to be compared with the numerical one.

The equation was discretized using a finite difference/ Central method in space (2nd order in space) and Euler implicit in time (1st order in time).

The global error behavior should follow asymptotically $\dfrac{A}{dx^{2}}+\dfrac{B}{dt}$ with $A$ and $B$ that do not depend on $dt$ and $dx$.

Now after implementing the resolution algorithm on Python, and constructing the error plot variation in terms of $dx$  (2nd norm and $\infty $ norm ) for a fixed $dt$, I was expecting a behavior following $C+\dfrac{A}{dx^{2}}$, with a "smooth" decrease in the convergence rate as C would become more predominant for smaller $dx$ till reaching the constant. What I got instead (and what I don't understand) is a small increase in the convergence rate, before recovering to the constant value ("reversed bell" behavior).

Is it normal to have an "acceleration" of the convergence rate in this case ?
Is it a problem related to rounding effect of the very small space error ?

Here is the ready-to-use Python code :

```<code>import numpy as np
import matplotlib.pyplot as plt
import scipy as sci
from scipy.sparse.linalg.dsolve import spsolve
import scipy.sparse as sp
import sympy as sym

#Symbolic calculation
x, t= sym.symbols('x t')
Dxx_1=1+x**2+sym.sin(x)

fi_1=1+sym.sin(sym.pi*x)*sym.sin(sym.pi*t) #Analytical solution
Soi_1=-sym.diff(Dxx_1*sym.diff(fi_1,x),x)+sym.diff(fi_1,t)#Source

Dxx_2=sym.lambdify((x),Dxx_1,"numpy")

#Conversion from symbolic to function
Soi=sym.lambdify((x,t),Soi_1,"numpy")
fi=sym.lambdify((x,t),fi_1,"numpy")

#Lists that will contain the trucation errors for plot
Err00=[]
Err2=[]

#list of the number of points in the domain
N=np.array([10,20,40,60,80,100,120,140,160,180,200,320,640])

#Loop on Nx
for Nx in N:

    #Space domain definition
    lx=1.0
    ox=0.0
    dx=(lx-ox)/(Nx)

    x1=np.linspace(ox+dx,lx-dx,Nx-1)#Points of the domain where f is unknown
    x2=np.linspace(ox,lx,Nx+1)#All points of the domain

    #Diffusion coefficient
    Dxx=Dxx_2(x2)

    #Initialisation
    f=fi(x1,0)

    #Time grid
    t1=0.5
    dt=0.0001
    Nt=int(t1/dt)

    #Implicit matrix definition - Sparse method
    diag0= -(2*Dxx[1:-1]/(dx**2))
    diag1= ((Dxx[2:-1]-Dxx[:-3])/(4*dx**2))+((Dxx[1:-2])/(dx**2))
    diag_1= (-((Dxx[3:]-Dxx[1:-2])/(4*dx**2))+((Dxx[2:-1])/(dx**2)))
    data = [-dt*np.append(diag_1,[0]),1-dt*diag0,-dt*np.append([0],diag1)]     # Diagonal terms
    offsets = np.array([-1, 0, 1])                     # Their positions
    J = sp.dia_matrix((data, offsets), shape=(Nx-1,Nx-1))

    #Iteration and resolution
    for k in range(1,Nt+1,1):
        #Boundary condition vector
        S1=np.zeros(Nx-1)
        S1[0]=fi(0.0,k*dt)*(-(Dxx[2]-Dxx[0])/(4*dx**2)+(Dxx[1])/(dx**2))
        S1[-1]=fi(1.0,k*dt)*((Dxx[-1]-Dxx[-3])/(4*dx**2)+(Dxx[-2])/(dx**2))
        #2nd member vector
        S=Soi(x1,k*dt)
        f=spsolve(J,f+dt*S+dt*S1)
        sol_ana=fi(x1,k*dt)
        print('iteration '+str(k)+' - t = '+str(k*dt)+' - Implicit - Nx ='+str(Nx))

    print('Nx = '+str(Nx))
    Err2.append(100*np.linalg.norm(f-sol_ana)/np.linalg.norm(sol_ana))
    print('Err2 = ' + str(Err2[-1]))
    Err00.append(100*(abs(f-sol_ana).max()/abs(sol_ana).max()))
    print('Err00 = '+ str(Err00[-1]))

plt.figure()
plt.plot(N,Err00,'o',N,Err2,'o',N,N[0]**2*Err00[0]*N**float(-2),N,N[0]*Err00[0]*N**float(-1))
plt.title(' Err = f(Nx) - 1D - DF - Dirichlet  - Impl - dt='+str(dt)+' - t= '+str(t1))

plt.xscale('log')
plt.yscale('log')
plt.legend(['$Err_{\infty}$','$Err_{2}$','$\dfrac{1}{Nx^{2}}$','$\dfrac{1}{Nx}$'])
plt.xlabel('Nx')
plt.ylabel('Err')
plt.annotate('$f_{ana}='+sym.latex(fi_1)+'$', xy=(0.4, 0.14), xycoords='axes fraction')
plt.annotate('$D_{xx}='+sym.latex(Dxx_1)+'$', xy=(0.4, 0.1), xycoords='axes fraction')
</code>```


问:
Neumann boundary condition FD implementation for instationnary diffusion equation
I am trying to solve this diffusion equation :
$\dfrac{\partial D\dfrac{\partial f}{\partial x}}{\partial x}+S = \dfrac{\partial f}{\partial t}$ ($D$ is not constant and varies according to $x$) with the following BC: $f(x,0)=1 , f(0,t)=0, \dfrac{\partial f}{\partial x}(1,t)=0$

I am using a central finite difference scheme (2nd order) for space and Euler explicit for time (1st order).
The discretized $[0..1]$ domain contains the $x_{i}$ points, $i \in [0..N]$
 ($x_{0}=0 , x_{N}=1$).

The implementation of Dirichlet condition at $x=0$ is simple.
For the Neumann boundary condition, I saw a lot of references that treated simple cases, adequate for constant diffusion coefficients and stationnary cases, like the ghost point, that I can't consider in my case(Diffusion coefficient can't be evaluated out of the domain, ghost point envolves evaluation the equation on a point that is not inside the interior domain... ).

So I am trying something at my own and I want confirmation or correction if I omitted an important aspect in my work:

The discretized equation evaluated at $x_{N-1}$ involves the value of $f_{N}$. If we had Diricihlet condition at $x_{N}$, the problem would have been solved directly. But in our case, $f_{N}$ is unknown.

What I did is that I treated the Neumann BC  with a backward scheme (at the 2nd order to preserve the general order of the whole scheme) in order to have an expression of $f^{k}_{N}$:

$\dfrac{\partial f}{\partial x}(x_{N},k\Delta t)=\dfrac{-3f^{k}_{N}+4f^{k}_{N-1}-f^{k}_{N-2}}{2\Delta x}=0$

$f^{k}_{N}=4/3 f^{k}_{N-1}-1/3 f^{k}_{N-2}$

then I injected the new expression of $f^{k}_{N}$ in the discretized equation at $x_{N-1}$. Meaning that I did not add an equation to the problem as I would do for a stationnary problem but I injected the expression of $f^{k}_{N}$ inside the last equation.

So the unknown vector $f^{k}$ in a case of matrix writing would only contain $f^{k}_{i}, i\in [1..N-1]$ (like in a full Dirichlet situation)

Is it correct ?

For information, the matrix form of the discretized problem is :
$f^{k+1}=Af^{k}+S+B$, with $B$ the vector containing the B.C , and $S$ the 2nd member vector


问:
Finite Element Stabilization for Drift-Diffusion/Advection-Diffusion Equations
I've tried my best to look through the relevant suggested similar questions when posting this, and hopefully this contains enough new material to not be considered a duplicate.

I'm currently trying to draw up a semiconductor drift-diffusion simulation framework, and the work that I'm supposed to be working off of, even though it's called 'finite element' in the section heading, is control-volume finite element, or vertex-centered finite volume (depending on which field you come from, one may be more common notation than the other).

As best as I understand it right now, the finite volume method preserve fluxes by construction, and is therefore preferred for implementation of conservation laws. It also allows the Scharfetter-Gummel Finite Box method to be implemented for current densities across the finite volume elements. The SG-FB method is a generalization of the original SG formulation, where the finite differences between carrier densities (the $\nabla n$ term below) are calculated using a special exponential smoothing scheme. This is done because $n$ can often change discontinuously between two points in space, and the SG method is a stable difference scheme.

However, I would like to implement this in finite elements, not volume, almost entirely because of the libraries available for FEM (deal.II, libMesh, FEniCS, ...). The relevant current density equations for electrons are below (these, combined with an identical equation for holes, and a nonlinear Poisson equation, make up the system)

$$
\frac{\partial n}{\partial t} = \frac{1}{q} \nabla \cdot \mathbf{J_n} + G_n - R_n
$$

$$
\mathbf{J_n} = qD_n\nabla n - q\mu_ nn\nabla \phi
$$

For the questions:

Are there any flux/mass-conserving schemes in FEM that have properties similar/identical to the flux/mass-conserving properties of FVM? I keep seeing in a few places that discontinuous-Galerkin conserves flux/mass, but I can't quite get a consensus on it.

What would be appropriate stabilization scheme(s) to apply to this system, given the rapid fluctuations in $n$ (and $p$, identically, for the hole continuity equation)? I keep reading about streamline-upwind Petrov-Galerkin, Galerkin least Squares, and shock-capturing terms, but this seems to be a generally open question.

Given that I keep seeing shock capturing terms in my research, is it reasonable to claim that rapid, discontinuous changes in carrier density are identical to 'static' shock waves/fronts?

I'm still quite new to these equations, having dealt with them only for the past few months. I'd be happy to provide any more information, if required.


问:
fourth order Poisson iterative solver --in Matlab
I want to calculate the stream function $\psi$ starting from a velocity field $(u,v)$ (such that $u=-\frac{\partial\psi}{\partial y}$ and $v=\frac{\partial\psi}{\partial x}$). I thus calculate the vorticity $\omega$ (such that $\omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$) and then solve the Poisson equation $\nabla^2 \psi=\omega$.

To easily specify the boundary conditions, I choose a successive over-relaxation (SOR) approach.

The code seemingly runs smoothly. However, when I check if the solution is correct by differentiating $\psi$, its derivatives do not match the initial velocities $(u,v)$

Is this due to some intrinsic limitation of the SOR?

```<code>close all
clearvars
clc

%coordinates
x=1:10:600; x = repmat(x',1,50);
y=1:10:500; y = repmat(y,60,1);
dx = 10*ones(60,50); %x spacings
dy = 5*ones(60,50); %y spacings

%define velocities
u = y.^3;
v = x.^3;
figure;pcolor(x,y, u);title('u')
figure;pcolor(x,y, v);title('v')

%vorticity
[dudy,~] = gradient(u, dy(1,1),dx(1,1));
[~,dvdx] = gradient(v, dy(1,1),dx(1,1));
figure;pcolor(x,y, dudy);title('dudy')
figure;pcolor(x,y, dvdx);title('dvdx')

omega = dvdx - dudy;

figure;pcolor(x,y, omega);title('omega');colormap(jet(10));colorbar

%iterative solver (SOR)
psi=rand(62,50); %initial guess

beta=dx./dy; %spacing ratio

%boundary conditions in y
omega(:,1:2)=0;
omega(:,end-1:end)=0;

%periodic boundary conditions in x -- append/prepend values on the other side
omega = [omega(end-1:end,:); omega; omega(1:2,:)];
beta = [beta(end-1:end,:); beta; beta(1:2,:)];

[M,N]=size(psi);

w=1.4; %relaxation factor

figure
for iter=1:1000
  for ii=3:M-2
    for jj=3:N-2

%%fourth order formula
psi(ii,jj) = (1-w)*psi(ii,jj) + w*...
                    ( 16*( psi(ii+1,jj) + psi(ii-1,jj) ) - psi(ii+2,jj) - psi(ii-2,jj) +...
                    beta(ii,jj)^2*( -psi(ii,jj+2) + 16*(psi(ii,jj+1) + psi(ii,jj-1)) -psi(ii,jj-2) ) +...
                    12*dx(ii,jj)^2*omega(ii,jj) )...
                   /(30*(1+(beta(ii,jj)^2)));

%%second order formula
%          psi(ii,jj) = (1-w)*psi(ii,jj) + ...
%              w*( (psi(ii+1,jj)+beta(ii,jj)^2*psi(ii,jj+1)+psi(ii-1,jj)+beta(ii,jj)^2*psi(ii,jj-1)) -...
%              2*dx(ii,jj)^2*omega(ii,jj) )...
%              /(2*(1+(beta(ii,jj)^2)));

    end
  end

      if mod(iter,10)==0
        pcolor(x,y, psi(2:end-1,:))
        colorbar
        title(num2str(iter))
        drawnow
      end

end
psi = psi(2:end-1,:); %resize back

figure;pcolor(x,y, psi);title('psi')

%%%test
[dpsidy,dpsidx] = gradient(u, dy(1,1),dx(1,1));

figure;hold on
subplot(1,2,1)
pcolor(x,y, u); colorbar; title('u orig')
subplot(1,2,2)
pcolor(x,y, -dpsidy); colorbar; title('u reconstructed')

figure;hold on
subplot(1,2,1)
pcolor(x,y, v); colorbar; title('v orig')
subplot(1,2,2)
pcolor(x,y, dpsidx); colorbar; title('v reconstructed')
</code>```


问:
Jacobian Elements for Coupled Drift-Diffusion System using Vertex-Centered Finite Volume
I'm trying to solve the fully coupled drift-diffusion system using Newton's Method. Although I eventually plan to potentially use a Jacobian-Free Newton-Krylov approach, this is still something that I want to implement and look into.

A bit of notation for the equations presented below:

$d_{ij}$ and $l_{ij}$ are control area/volume constants, not physical entities. The $i$ and $j$ indices are for the central point being considered ($i$), and all neighbouring points that are connected ($j$'s). $A_{\Omega_i}$ is the control area/volume associated to each vertex/node, and is also constant.

All quantities that are explicitly required during calculation are at node points. In the equations $B(\Delta_{ij}) \equiv B(\phi_i - \phi_j)$, the Bernoulli function, used as a smoothing method for the differences. This difference should technically exist between node points, but that's for another time.

Here are the equations: $\phi$ is electrostatic potential, $n$ and $p$ are carrier concentrations. These are the independent variables to solve for. Everything else is either a constant, or a function of one of these variables.

$$\sum_{j\ne i} \frac{d_{ij}}{l_{ij}}(\phi_i - \phi_j) = -\frac{q}{\epsilon}(p_i - n_i - N_{D_i}^+ - N_{A_i}^-)A_{\Omega_i}$$

$$\frac{\partial{n_i}}{\partial t}A_{\Omega_i} = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{n_iB(\Delta_{ij}) - n_j B(-\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i}$$

$$\frac{\partial{p_i}}{\partial t}A_{\Omega_i} = D_p \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{p_iB(-\Delta_{ij}) - p_j B(\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i}$$

In residual form, and assuming equilibrium conditions so all time derivatives go to 0, these equations are:

$$\Phi = \sum_{j\ne i} \frac{d_{ij}}{l_{ij}}(\phi_i - \phi_j) +\frac{q}{\epsilon}(p_i - n_i - N_{D_i}^+ - N_{A_i}^-)A_{\Omega_i} = 0$$

$$N = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{n_iB(\Delta_{ij}) - n_j B(-\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i} = 0$$

$$P = D_p \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{p_iB(-\Delta_{ij}) - p_j B(\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i} = 0$$

Now, my question is whether I have these Jacobian elements correct. Here is the Jacobian element formula that I've tried to come up with. Only a few equations have been derived, since the remainining elements are just trivial variations of these 3 formulations. This post is mostly about getting more experienced eyes to look at these equations, and potentially see something that I haven't quite seen yet. Are these formulations for the Jacobian elements correct, or is there something wrong with them that I'm missing?

$$\frac{\partial\Phi}{\partial\phi_i} = \sum_{j\ne i} \frac{d_{ij}}{l_{ij}} +\frac{q}{\epsilon}(\frac{\partial p_i}{\partial \phi_i} - \frac{\partial n_i}{\partial \phi_i})A_{\Omega_i}$$

$$\frac{\partial{N}}{\partial \phi_i} =  D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \left\{\frac{\partial{n_i}}{\partial \phi_i}B(\Delta_{ij}) + n_i \frac{\partial B(\Delta_{ij})}{\partial \phi_i})\right\}\right]
+ \frac{\partial(G_i - R_i)}{\partial \phi_i}A_{\Omega_i}$$

$$\frac{\partial N}{\partial n_i} = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} B(\Delta_{ij})\right] + \frac{\partial (G_i - R_i)}{\partial n_i}A_{\Omega_i}$$


问:
Solving numerically a linearized system of elliptic (?) Navier-Stokes equation (Shallow Water Derived)
For my PhD Thesis, my advisor asked me to build a solver inspired from the article "Optimal Control Theory Applied to an Objective Analysis of a Tidal Current Mapping by HF Radar, J-L Devenon, 1989". In this article, in the numerical fluid-solving part, the author linearize the shallow water 2-D equations by removing the advection terms and using a linear friction coefficient. Then he uses a Fourier transform on the equation to finally end up with the system:

$$
\begin{aligned}
-i\sigma\vec{U}+\vec{f}\times\vec{U}+g\vec{\nabla}h&=-r\vec{U}\\\\
i\sigma h+\vec{\nabla}\cdot\left(H\vec{U}\right)&=0
\end{aligned}
$$

Here $i$ is the basis for the imaginary number, $\sigma$ the frequency of the given Fourier component, $\vec{f}$ the Coriolis force vector, $g$ the gravitational acceleration, $h$ the local sea-level elevation compared to "mean local level", $r$ the linear friction coefficient, $H$ the local depth and $\vec{U}$ the local velocity.

Then, this system is simplified in a singled vector equation of the form :

$$
F(\vec{V})=\lambda\vec{V}+\vec{f}\times\vec{V}+\vec{\nabla}\left(\vec{\nabla}\cdot\left(\gamma\vec{V}\right)\right)=\vec{0}
$$

Where:

$$
\lambda = rH-i\sigma,\quad \gamma=\frac{-igH}{\sigma}
$$

And $\vec{V}$ is the transport, given as the product of the speed (supposed constant on the entire water-column) and the local height of the water-column( supposing the fluctuations $h$ are small compared to the overall depth).

He then solves the Dirichlet problems with specified boundary conditions. For this, he uses an iterative relaxation method (from Smith, G.D, Numerical solution of partial differential equations: Finite differences methods, 1978). And use only centered discretisation schemes.

I then wrote my own code for solving it. My first attempts where on a rectangular bassin with fixed boundary conditions. However It did not seem to work with any condition that differed from 0. I also tried to figure out analytical solutions but the only (very obvious) solution I found was the null solutions everywhere. Here is an example of the code I used:

```<code>############################################
# First test of implementing Devenon's     #
# Equation system and solving them in a clo#
# -sed basin                               #
############################################

#*****Packages import**********************#
import numpy as np

#****Functions*****************************#

def  relaxation_scheme_Jean_Luc_overrelaxed(U,V,dy,dx,nx,ny,gamma,lamda,omega):
 """Function use to solve 2D linearized equation
 of Navier-Stokes for a rigid lid. From Devenon 1989
 . We're trying an overrelaxation technique with omega
  as the relaxation parameter.

 INPUT:
 -----
 U,V: Zonal and meridional speed, size nx*ny, also the output,complex valued
 nx: Number of grid cell in x direction
 ny: Number of grid cell in y direction
 dx: Size of a grid cell in x direction
 dy: Size of a grid cell in y direction
 b: Forcing source function, of size nx x ny
 gamma: Local complex phase speed
 lamda : Local complex wave vector
 omega: Relaxation factor, must belong to [1:2]

 OUTPUT:
 -------
 U,V: Already specified earlier """

 # First we compute the non-homogeneous terms

 #U_b=f*V[1:-1,1:-1]+ (gamma[:-2,:-2]*V[:-2,:-2] +
 #       gamma[2:,2:]*V[2:,2:]-gamma[2:,:-2]*V[2:,:-2]-
 #       gamma[:-2,2:]*V[:-2,2:])/(4*dy*dx)

 #V_b= -f*U[1:-1,1:-1]+ (gamma[:-2,:-2]*U[:-2,:-2] +
 #         gamma[2:,2:]*U[2:,2:]- gamma[2:,:-2]*U[2:,:-2]-
 #         gamma[:-2,2:]*U[:-2,2:])/(4*dx*dy)

 U_b=V[1:-1,1:-1]*0.
 V_b=V[1:-1,1:-1]*0.
 # Now we can use relaxation technique with an external term
 U[1:-1,1:-1] = 1/lamda[1:-1,1:-1] * (omega*gamma[2:,1:-1]*U[2:,1:-1]+
                2*(1-omega)*gamma[1:-1,1:-1]*U[1:-1,1:-1]+(omega-1)*gamma[:-2,1:-1]*U[:-2,1:-1]) \
                /dx**2  +1/lamda[1:-1,1:-1]*U_b

 V[1:-1,1:-1] = 1/lamda[1:-1,1:-1] * (omega*gamma[1:-1,2:]*V[1:-1,2:]+
                 2*(1-omega)*(gamma[1:-1,1:-1]*V[1:-1,1:-1])+(omega-1)*gamma[1:-1,:-2]*V[1:-1,:-2]) \
                 /dy**2 +1/lamda[1:-1,1:-1]*V_b

 return U,V

#*****Parameters and variables*************#
nx = 40
ny = 40 # Number of grid points in each direction
nt  = 100 # Number of iteration of the relaxation method
xmin = 0.
xmax = 12000.
ymin = 0.
ymax = 12000. # Dimension (in meters) of the grid

dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)

# Initialization
# I understood that each speed would be complex-valued
# However, I wondered if U and V could only be the real and imaginary
# Part of the same variable in this case
U,V  = np.zeros((nx,ny),dtype=complex),np.zeros((nx,ny),dtype=complex)
x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(xmin, xmax, ny)
X, Y = np.meshgrid(x, y)

H=10 # Setting constant depth for testing
sigma=2.31484e-5 # Approximately semidiurnal frequency,  for tide (in hertz)
g=9.81 # Gravity acceleration
K=2.5e-2 #Friction coefficient
f=2e-4
lamda=np.ones((nx,ny),dtype=complex)*(K-sigma*1j)# (Friction and wavenumber)
gamma=np.ones((nx,ny),dtype=complex)*-1j*g*H/sigma #(Gravity and wavenumber)

#Initial conditions
U[0,:]=1.
U[-1,:]=1.
U[:,0]=1.
U[:,-1]=1.

V[0,:]=1.
V[-1,:]=1.
V[:,0]=1.
V[:,-1]=1.

#********Solving****************************#

for i in range(nt):
    U,V=relaxation_scheme_Jean_Luc_overrelaxed(U,V,dy,dx,nx,ny,gamma,lamda,2.)
</code>```

For any conditions, I used except for the constant ones I found whether a very stable behaviour (the system stayed exactly in its initial position) or a very unstable one (I have instabilities growing from the sides that propagate, very quickly giving me aberrant values). This behaviour seems to be driven by the parameters lambda, gamma as well as the size of the cells. However no combination I have tried yet helped me finding a solution. This leads me to believe that I did a very bad job at coding this.

I tried searching the site ( and the web in general for answers) and found many links for nice EDP and even Navier-Stokes solvers. Though I did not get anything that I could link to my problem, may be due to a poor comprehension of it. This being said, I have another constraint that my advisor would like this to become an educational and simple problem. I, therefore, have to understand the internal functioning of my solver.

Then I have a few questions for you :

1) The first is obvious: Why doesn't it work? Did I do a bad job at coding? What should I correct to have a correct working algorithm?

2) The second is less obvious: Can I write my vector V in a complex manner, using one component of the transport as the real part and the second as the imaginary part? This seems obvious to my advisor but it does not look obvious to me at all.

3) The paper refers to it as an "elliptic" problem, close to a Helmholtz one. However, I just see two ODE's with a coupling term. Can I talk about an elliptic problem in this case?


问:
Acceleration of matrix geometric series
Suppose we want to find $x$ such that:

$$x=b+Ax$$

where $A$ is a large sparse square matrix with eigenvalues in the unit circle.

There are two representations of the solution:

1)
$$x=(I-A)^{-1}b,$$

2)
$$x=\sum_{k=0}^\infty{A^k b}.$$

Since we cannot practically store powers of the matrix $A$ (due to memory limits: indeed in practice, we just have a function that computes $Ax$ given $x$), using the representation 2) means evaluating the recurrence:

$$x_t=b+A x_{t-1}.$$

However, since $A$ may have eigenvalues close to $1$ in magnitude, this will converge slowly.

Is there a multivariate series acceleration method which would handle such recursions? (E.g. a multivariate version of the Aitken "delta-squared" process/Steffensen's method.)

Alternatively, is there a technique for solving 1) directly that is guaranteed to have good performance? Do our assumptions on $A$ imply anything about the performance of GMRES / CGS / BICGSTAB(L) / QMR etc.?

Note: This question: "Is there a faster method to compute the geometric series of a matrix?" is similar but seems to be interested in finite sums of small dense matrices, rather than infinite sums of large sparse matrices.

Title: Acceleration of matrix geometric series


问:
Looking for a version of DSDP that is less prone to integer overflows than the original
I am working on a problem that involves semidefinite programming (constrained optimization of fairly large positive definite matrices). The software is written in C++ and calls DSDP 5.8 to solve the SDP problem. It works as advertised on small/medium test problems (up to 200 x 200 matrices, ~$10^5$ constraints), but dies due to indexing errors on larger tasks. The reason is, of course, that the authors of DSDP used int variables for indices (and int[]-s for indexing arrays) instead of e.g. size_t or long and with the larger problems the 4 byte int-s eventually overflow.

Before I start modifying a C source that has not been updated since 2006, I thought I might ask the community here: does anyone have, by any chance, a modified version of DSDP that uses size_t or long for indices (and size_t[] for the index arrays)? These types are usually 8 bytes long on modern 64-bit platforms.

My other option would be to swap out DSDP and use another SDP solver instead. Please understand that I am reluctant to do this because I would have to adapt my software to a new API. Having said that, feel free to recommend open-source SDP solvers for first-order cone programming that can deal with large problems. I am aware of CVXOPT but that is written in Python; and MOSEK is not FOSS as far as I can tell.


问:
2D reaction-diffusion system simulation
I am a complete beginner in numerical simulation and I am pretty lost about how to tackle this problem.

I have been trying for some time to find the steady state (or simulate), the following system

$$ \dot S = D_S \Delta S + \rho_S \frac{ F^n } { F^n + 1} - \delta_S S -\delta_C P^2 S $$

$$ \dot P  = D_P \Delta P + \rho_P - \delta_P P + (v - 2 \delta_C) P ^ 2 S $$

$$ \dot F = \Delta + \rho_F \frac{ 1 } { (P^2S)^n + 1 } - \delta_F F $$

on a mesh that looks like this (with no boundary conditions given):

I realize that this systems contains the temporal derivatives, so I'm not explicitly going for the steady state, but when I took out the TransientTerm it complained that it was Exactly singular.

```<code>S = CellVariable(mesh = mesh, hasOld=True)
P = CellVariable(mesh = mesh, hasOld=True)
F = CellVariable(mesh = mesh, hasOld=True)

sdot = TransientTerm(var=S) == DiffusionTerm(coeff=Ds, var=S) + \
rs * pow(F,2) / (pow(F,2) + 1) - ds * S - dc * pow(P, 2) * S

pdot = TransientTerm(var=P) == DiffusionTerm(coeff=Dp, var=P) + \
rp  - dp * P  + (v - 2. * dc) * S * pow(P, 2)

fdot = TransientTerm(var=F) == DiffusionTerm(coeff=1, var=F) + \
rf / (pow(S * pow(P, 2), 2)  + 1.0 ) - df * F

dt = 0.005
coupled.sweep(dt = dt)

for i in range(0, 100000):

    v1.plot('plots/SHH_'+str(dt * i)+'.png')
    v2.plot('plots/Ptc_' + str(dt * i)+'.png')
    v3.plot('plots/FGF_'+str(dt * i)+'.png')

    S.updateOld()
    P.updateOld()
    F.updateOld()

dt = 0.005

for i in range(0, 100000):

    v1.plot('plots/SHH_'+str(dt * i)+'.png')
    v2.plot('plots/Ptc_' + str(dt * i)+'.png')
    v3.plot('plots/FGF_'+str(dt * i)+'.png')

    S.updateOld()
    P.updateOld()
    F.updateOld()
    coupled.sweep(dt = dt)
</code>```

The simulation initially appears to be proceeding as normal, but eventually, regardless of the size of my time-step, the $S$ values and $P$ values drop below zero and subsequently $S$ blows up do its large production rate $\rho_s = 1600$, causing an overflow error.

I should mention that when I simulated this problem in 1D using pdepe, the simulation was find and converged to a steady state.

How do I debug this, and is there other problems with my simulation? Should I try a different package such as Firedrake? If so, how do I solve the convergence problems there?

Thanks very much in advance.

Title: 2D reaction-diffusion system simulation


问:
Are planar 2nd order NURBS curve segments exactly the segments of conic sections?
I think it is well known that 2nd order NURBS curve can describe any segment of a conic section. Such a segment has 7 degrees of freedom, take an ellipse arc for example: the ellipse itself has 5 degrees of freedom (center, half-axes and rotation), and the two points on the ellipse (2 angles) give the remaining 2 degrees of freedom. A 2nd order NURBS curve segment has less than 9 degrees of freedom: 3 control points and 3 weights, but the scaling of the weights is irrelevant. Can one discuss away the remaining additional degree of freedom?

(I think I can show that a planar 2nd order B-SPLINE segment is a segment of a parabola or of a line. The argument that both have 6 degrees of freedoms applies, but the algebra is also easy to work out.)


问:
Data transfer in the context of adaptive mesh in finite element method (FEM)
In the context of adaptive mesh in FEM, after a new mesh is created, the data on the old mesh are to be transferred to the new mesh. For the data on the integration points (IPs), it seems the usual way is to (1) transfer the data on IPs of the old mesh to the nodes of the old mesh, (2) use the shape function to transfer the data from the nodes of the old mesh to the nodes of the new mesh, (3) then use the shape function to transfer the data on the nodes of the new mesh to the IPs of the new mesh.

There are different techniques for the first step (see (1) above), and the most popular one seems to be the patch method or superconvergent patch recovery. In this method the IPs of the old mesh are used as the sampling points.

My question would be why not just simply use the patch method to transfer the data on the IPs of the old mesh to the IPs of the new mesh directly, as the IPs of the old mesh are used as sampling points anyway.


问:
Method of manufactured solutions - choice of type of boundary conditions
I am attempting to use the method of manufactured solutions (MMS) for code verification for linear elasticity. However, this is more of a general question regarding the general use of MMS.

In MMS, the manufacturer makes up some solution that satisfies the PDE, and then we obtain values for BCs and sources that we feed into the numerical code and solve. Does the manufacturer have the freedom in choosing whether a boundary is a Dirichlet or Neumann boundary condition, and then evaluating the manufactured solution at the boundary in accordance with their choice of the boundary condition type?

For example, in 1-D steady heat-conduction, I manufacture some solution $T(x)$. Do I have the freedom of using either Dirichlet BC at $x=0$ or Neumann BC? For a Dirichlet BC, I simply evaluate $T(x)$ at $x=0$ and that's my BC. On the other hand, I could instead use a Neumann BC at $x=0$ so I would then evaluate $\frac{dT}{dx}$ at $x=0$. Is the manufacturer in full of control of what type of BC is to be prescribed?


问:
Limit to volume change in a discretized mathematical model?
I have set up a mathematical model describing the diffusion of ozone out of a gas bubble. The bubble is surrounded by a thin gas film. So actually, the model describes the diffusion of ozone through this gas film.
The mathematical model is created by discretizing the volume from $V_0$ to $V_1$. Where $V_0$ represents the outer volume of the bubble, and $V_1$ represents the outer volume of the gas film (surrounding the bubble). The discretized scheme consists of $N$ equally spaced volume elements from $V_0$ to $V_1$ (finite difference method):
\begin{equation}
    \Delta V = \frac{V_1 - V_0}{N}
\end{equation}
The volume of the bubble changes as a function of the amount of ozone leaving the bubble, which in turn changes as a function of time.
The volume of the bubble and the amount of ozone inside the bubble is linked by the ideal gas law:
\begin{equation}
    V_0 = \frac{n_\text{total}(t) \cdot R \cdot T}{P}
\end{equation}
$n_\text{total}$: the total amount of gas inside the bubble, $T$: the temperature, $R$: the gas constant, and $P$: the pressure.

The bubble does not only contain ozone. It also contains inert gases, so that:
\begin{equation}
    n_\text{total} = n_\text{ozone}(t) + n_\text{inert}
\end{equation}
$n_\text{inert}$ will remain constant and only $n_\text{ozone}(t)$ will change over time.

There should be a limit for how much the volume can change before numerical errors will start to occur. Beyond this limit, the discretization scheme should break down and cause errors. How do I express this limit?

Is the limit given by:
\begin{equation}
    Ratio = \frac{V_{0,initial}}{V_{1,initial}}
\end{equation}

So that the volume change must not exceed the ratio of the two initial volumes of the bubble and the gas film?


问:
Iteratively solving a sparse, ill-conditioned system
I have a sparse (density = 0.2%), ill-conditioned system that I am trying to solve, with no luck.

Background

I have a sequence of sampled data, where two of every 8 samples have been zeroed due to a bug. The sequence is two sinusoids in the presence of white noise.

What I'd like to do is determine the Signal to Noise Ratio of the original sequence. I understand that there is information that is forever lost, and so I won't get an exact answer, but that's acceptable, as long as I can get something reasonably close using the existing data, without any filtering (i.e. no interpolation).

Strategy

Having a sequence with missing samples as I do is like multiplying the original sequence with a square wave (75% duty cycle, in my case). In the frequency domain, this is equivalent to a periodic convolution of regularly-spaced impulses with the DFT of my original sequence. Essentially, its like doing circular shifts on the input sequence, multiplying particular points by certain weights according to the square wave DFT, and summing them together. This can be represented as a linear system.

So I've generated the sparse matrix representing the convolution, and am now trying to solve Ax=b where A is the sparse matrix, x is the DFT of the original sequence, and b is the DFT of the corrupt sequence.

Unfortunately, the sparse matrix is ill conditioned. Using matlab, neither mldivide nor svds provided useful answers. The thing is, I know what the input sequence (and its DFT) look like. It's a DC component, two bin-centered sinusoids, and white noise. So I have an excellent initial guess, I just can't figure out a way to iterate on that. Also worth noting is that although the DFT is complex-valued, and these operations are all on complex numbers, I only care about the magnitudes.

If I solve mag(A)x=mag(b) I get a decent result, but its slightly too optimistic, so I'm hoping an iterative solution on the complex values will be more accurate.

My criteria for an acceptable solution is one that minimizes the mean square error between the magnitudes of  Ax and b.

Any help would be greatly appreciated! I have little experience in these sorts of linear systems, but am doing my best to learn the appropriate methods.

This is the code I'm using to test my method: https://pastebin.com/F3bArgnR


问:
Decrease in slope during convergence analysis
I am using the method of manufactured solutions to perform the order of accuracy testing. I am using a cube for the testing. The cube is size 1m on all sides.

I used 5 refinements:

$dx = dy = dz = 0.5$ ($8$ cells)

$dx = dy = dz = 0.5/2$ ($8^2$ cells)

$dx = dy = dz = 0.5/2/2$ ($8^3$ cells)

$dx = dy = dz = 0.5/2/2/2$ ($8^4$ cells)

$dx = dy = dz = 0.5/2/2/2/2$ ($8^5$ cells)

I first performed the analysis for a pure Dirichlet BC problem and observed 2nd order accuracy, as expected.

Next, I changed one of the cube surfaces to be a Neumann condition, the other 5 surfaces are still Dirichlet conditions. The results for this are shown below.

So the slope here goes from about 2 to 1. I am still expecting 2nd order accuracy with the Neumann condition.

Assuming that my math is correct and the scheme that I am using is 2nd order accurate, does the decrease in slope suggest that I have a problem with my implementation?

Could it also signify the onset of some other sort of error (perhaps algebraic?) taking over the discretization error? I would think not given that we did not observe this in the pure Dirichlet case.

I am puzzled as to why the first segment of this 2nd figure indicates 2nd order accuracy and then gradually decreases. Any insight would be greatly appreciated.

ContentLicense: CC BY-SA 4.0


问:
computing dual matrix trace norm and tensor gradient in python
I'm trying to write the following function in python:
$$
f_\mu(\mathcal X) = f_0(\mathcal X) + \sum_{i = 1}^n \max_{||\mathcal Y_{i(i)}|| \leq1} \alpha_i\langle \mathcal X_{(i)},\mathcal Y_{i(i)} \rangle - \frac{\mu_i}{2}||\mathcal Y_{i(i)}||_F^2 \\ \mathcal X, \mathcal Y \ \text{are tensors (MxNx3), } \\ \mathcal X_{(i)} \ \text{is the tensor unfolded along mode i (a matrix), } \\||.||_F \ \text{ is the Frobenius Norm}, \\ \alpha, \mu \in {\mathbb R}, \\ n = 3
$$

What i have so far:

```<code>import tensorly as tl
import numpy as np
LA = np.linalg

def f_m(a,X,Y,m):
   result = np.zeros((3,3))
   zz = []
   for i in range(3):
      for j in range(3):
        yy = tl.unfold(Y[j],j)
        xx = a[j]*tl.unfold(X,j)
        zz = np.multiply(yy,xx)
        zz = LA.norm(zz, 'nuc')
        result[i,j] = zz-0.5*m[j] * LA.norm(tl.unfold(Y[j], j), 'fro') ** 2
   F1 = max(result[0,:])
   F2 = max(result[1,:])
   F3 = max(result[2,:])
   F = max(F1,F2,F3)
return F

def grad_f_m(m,a,X):
   grad = a*truncate_op(a/m * X, 1)
return grad

def truncate_op(X,tau):
   u , sig, v = LA.svd(X, full_matrices=False)
   t = np.zeros(len(sig))+tau
   np.minimum(sig,t,sig)
   trunk = np.dot(u[:,:], np.dot(np.diag(sig[:]), v[:,:]))
return trunk
</code>```

First question:

What does $ f_\mu(\mathcal X)$ return? A scalar? A matrix? A tensor? I'm thinking its a scalar.

Second question:
$$
g(X) := ||X||_* = \max_{||Y||\leq1} \langle X,Y \rangle
\\||.||_* \ \text{ is the Trace Norm}
$$
If i compute the trace norm like this:

numpy.linalg.norm(zz,'nuc')
(which computes the nuclear norm according to the documentation)

Is that equivalent to $$ \max_{||\mathcal Y_{(i)}||\leq1}\alpha_i\langle \mathcal X_{(i)},\mathcal Y_{(i)} \rangle $$

Third question:

$$
\nabla f_\mu(\mathcal X) = \nabla f_0(\mathcal X) + \sum_{i = 1}^n \alpha_i T_1(\frac{\alpha_i}{\mu_i}\mathcal X_{(i)})
$$
$$
T_{\tau} = U\sigma_{\overline \tau}V
$$

$$
\sigma_{\overline \tau} = diag(\min(\sigma,\tau))
$$

Where

$$
U,\sigma, V = SVD(\mathcal X_{(i)})
$$

Since $\ T_{\tau}$ will result in a different sized matrix for each $ \mathcal X_{(i)}$ how can I sum them across $ n $ ?


问:
Book Suggestion for Approximating Integrals using Random Partitions
Suppose I want to approximate the integral $\int_0^1 x^2\,dx$ using Riemann Sums or Darboux sums over random partitions of the interval $[0,1]$, Like in the image below:

Here, A "random" partition of $[0,1]$ containing $16$ partition points,

$$P=\{0, 0.048499, 0.097571, 0.25324, 0.28515, 0.4087, 0.45946, 0.53824,  0.62423     0.67352, 0.72305, 0.84254, 0.85201, 0.87255, 0.91006, 1\}$$

was created randomly with corresponding lower Darboux sum (showed in red) = $0.292536$ and upper Darboux sum (showed in blue) = $0.377436$, And we've got an approximation of  $\int_0^1 x^2\,dx = \frac{1}{3} = 0.333333\ldots$

Which Books or other resources do you suggest to read that treats rigourously problems like this, Both theoretical (Just the mathematical theory used to address such problems) and practical (Approximating integrals on a computer in such manner)?

There are some challenges I was thinking about that need addressing:

First, What is a random partition? It needs to be rigorously defined...

Second, What about partitions that their points are distributed according to other probability distributions such as the Gaussian distribution or others, And how can we know which probability distribution will give the best result (probabilistically) to approximate an integral?

Third, How can we answer questions such as, "What is the probability that the approximated value of the integral is within $0.5$ from the exact value of the integral?", i.e. How can we calculate probabilities like $$P\left\{\int_0^1 x^2 \,dx - 0.5 < \text{LOWER SUM}<\int_0^1 x^2\,dx < \text{UPPER SUM} < \int_0^1 x^2\,dx + 0.5\right\} \text{?}$$

Thanks...


问:
How to make supercell of partial occupancies for DFT?
In the case of partial occupancies (for example, this CIF), a supercell should be made for DFT calculations. However, numerous supercells can be built from a lattice (in the above example, 200 different supercells).

Which supercell should be used for the DFT calculations?

Many DFT works (if not most) are based on partial occupancies (e.g., dopants). How is a typical supercell chosen for the DFT calculation?


问:
Fixed point iteration reduction factor
In a book for solving a nonlinear differential equation with $N+1$ points, $u_{xx} = e^{u}, u(-1)=u(1) = 0$, in $[-1,1]$ with homogeneous Dirichlet boundary conditions, the fixed point iteration is used. The MATLAB code is as follows:

```<code>u = zeros(N-1,1);
change = 1; it = 0;
while change > 1e-15
       unew = D2\exp(u);
       change1 = change;
       change = norm(unew-u,inf);
       change/change1
       u = unew; it = it+1;
end
</code>```

where $D^2$ is a differential matrix. In each step the error is reduced by a factor of (change/change1) 0.2943. In the book it is asked how this factor is obtained? and what is its relationship with the eigenvalues of the $D^2$ matrix?

I evaluated the eigenvalues of the matrix $D^2$, but could not find a relation. How is the reduction factor in fixed point method obtained?

Thanks in advance.

Title: Fixed point iteration reduction factor


问:
Algorithm for group forming: as individual or in a preformed group
I have 20-80 users and 5-10 events with varying ranges of minimum and maximum number of free seats (2-4, 3-5, 2-6...). For example, with a range of 3-5 it is acceptable to only assign three users to the event. The users can choose weighted preferences (1,2,3) for 1-3 events. The problem is to assign as many as possible users to the events in an optimal way.

I'm using Hungarian algorithm and the JavaScript implementation is done using munkres-js. If there is no solution, events and finally users can be dropped until there is a solution. This produces satisfactory non-optimal solutions.

My current problem is how to allow the users to form predefined groups of varying sizes (2-8) before the optimization? The group cannot select an event if the maximum number of seats is smaller than the size of the group. This can be handled before the algorithm execution. The group can choose three weighted preferences in the same way as individual users. For example, out of 40 users there are three predefined groups of 4, 5, and 6 people and the remaining 25 users will be assigned individually.

My initial idea was to assign all users to predefined groups and then assign the groups to the events but I could not figure how to do it with varying groups sizes and event seats.

Non-optimal approximations and constraints that don't change the scope too much are also fine. There is a time constraint of ~30 seconds.

Any ideas?

// EDIT

Proposed constraints:

Limit group size to two people. This is feasible because common scenario is that people want to join an event with one friend.

// EDIT2

List of requirements:

Goal

The goal is to assign as many users as possible to the events.

User preferences must be honored.

Events

Events have minimum and maximum number of seats (for example: 3-5).

If an event doesn't receive the minimum number of users, the event must be cancelled.

Individual users

Users can choose up to three (1-3) weighted preferences (1,2,3) for as many different events (1-3).

Users can only be assigned to events they have chosen.

If there are more users than free seats, the excessive users must be removed fairly.

If users cannot be matched to their preferences, they must be removed. This should be the last measure.

Groups

Users can form groups of varying sizes (2-8).

Groups can choose three weighted preferences like individual users.

The group cannot select an event if the maximum number of seats is smaller than the size of the group. This can be handled before the algorithm execution.

Other

Algorithm must be able to execute in under 30 seconds.

Non-optimal approximations are fine.

New constraints can be introduced if they don't change the scope too much.


问:
Space covering optimization
I have the following problem:

In the space $E=\{1, 2, \dots, N_x\} \times \{1, 2, \dots, N_y\}$, I want to define $N_R$ rectangles $R_k=\{x_k^0, \dots, x_k^1\}\times\{y_k^0, \dots, y_k^1\}$ which cover part of my space $E$.

The total size of the union of all rectangles has a upper bound $C$ (constraint i). And the rectangle size is of course lower-bounded to $s_x$ on the 1st dimension (constraint ii) and $s_y$ on the 2nd one (constraint iii).

Each point of the space $E$ is associated with a score value by the function $f: E \rightarrow \mathbb R_+$.

I want to maximize this objective function:

$$\max_{N_R}\max_{R_1, \dots, R_{N_R}\\ s.t. |\bigcup_{k=1}^N R_k|\le C \> (i) \\ |x_k^1 - x_k^0| \ge s_x \quad (ii) \\ |y_k^1 - y_k^0| \ge s_y \quad  (iii)} \sum_{z\in E}\sum_{k=1}^{N_R}1_{z\in R_k}f(z)$$

To express it in a non-mathematical way, I want to find a set of rectangles such that the sum of the scores of the points they contain is maximal.

In my case, $E$ is quite big and I must do that for several different score functions f, so that it is impossible to brute-force all possibilities.

I do not know how to wisely optimize this function. Do you have some pieces of advice? Or could you lead me to existent algorithms to do that?

I have thought of (and implemented) suboptimal processes such that fixing a rectangle size (10 x 10 for instance) and solving the problem in this simpler setting.

Title: Space covering optimization


问:
Finite difference methods for coupled 2nd order nonlinear pdes
I have a system of coupled nonlinear PDEs that I cannot figure out how to solve in a smart way using FDM, so I was hoping someone here might have a clue.
The equations go as:

\begin{align*}
\frac{1}{n_\alpha}\frac{\partial}{\partial t} \left[n_\alpha-A_\alpha\left(\nabla^2\phi+\nabla^2 p_\alpha\right) \right] &=\text{RHS}_{n_\alpha} \\
\frac{\partial}{\partial t} \left[\ln(p_\alpha)- B_\alpha\left(\nabla^2\phi+\nabla^2 p_\alpha\right) \right] &=\text{RHS}_{p_\alpha}\\
\frac{\partial}{\partial t} \left[\sum_\alpha C_\alpha\left(\nabla^2\phi+\nabla^2 p_\alpha\right)\right] &= \text{RHS}_\text{vort}\\
\end{align*}
Where $A_\alpha, B_\alpha, C_\alpha$ are constants and $\alpha$ indicates the "species". I.e. it could range from $1$ to $n$ and so there would in the end be $2n+1$ equations that need to solved.

I suppose one would start with time stepping the equations, so one "only" needs to solve what seems to be coupled with nonlinear Helmholtz-type equations. But then what?


问:
Efficient initial identification of solid or liquid domains for a block structured Cartesian grid generation system
INTRO

Within the last 5 days I was able to generate a block structured Cartesian grid
generation system with a combination of Fortran,C++ and Python.

I am running intersection tests of the imported geometry with the bounding box
it is surrounded by and if the intersection test succeeds then I am subdividing the parent bounding box into 4 children boundary box.

As it can be understood for this 2D case, I am using Quad Tree decomposition.

And now I am planning to code a FVM / FEM solver for it.

PROBLEM

However, I am unable to think of an algorithm that would help me identify cells which are SOLID (i.e inside the body) and FLUID ( i.e outside the body) at a fast speed.

I can run a double ray intersection test by casting two opposite rays out of all the cells and see if both of the rays intersect with the boundary of the geometry or not to identify the solid domain. However such a test will be computationally very expensive to run after we have decomposed the quad tree
to a certain mesh level.

The figure given below has 40,577 quad elements, so it is not very high.
But future meshes will have more elements. And hence this particular problem can become a performance bottleneck.

FIG 1 : NACA 2412 Aerofoil meshed

FIG 2 : Zoomed view of trailing edge demonstrating mesh distribution

QUESTIONS

i) Can I perform any kind of initial tests to reduce my computational load ?

ii) Is there any other algorithm that would be more suited for this problem ?


问:
Fast Algorithms for the Simplicial Decomposition of a Convex Polytope in N-Dimensions
I'm in the process of constructing an algorithm which computes the Voronoi diagram of a set of points, but I now need a method to decompose each Voronoi cell into simplices. The information we have is:

The vertices of the polytope

The faces of the polytope and their vertices

The location of a central point

The radial lines passing through the faces

The fact that the polytope forms a convex hull

The main issues are:

The faces of the polytope can be non-simplicial

We want the total algorithm to scale well with increasing dimension, but each Voronoi polytope has an average number of vertices which increases with dimension. Unfortunately, I can't find useful information on exactly how.

EDIT: While I have no formula for the increase in the number of polytope vertices with dimension, testing shows it to be $\mathcal O(5^d)$.

If all faces were simplicial I'd draw simplices from the central point to each face. I can triangulate the non-simplicial faces and apply the above algorithm, but I still have the same question, just in a lower dimension, as each face has the above properties, except that we lack knowledge of its faces.

I'm using the Delaunay triangulation via the Bowyer-Watson algorithm at the moment, as this was the method I used for computing the Voronoi diagram. There are clearly much faster approaches, but before I spend too long implementing one of them, I felt it wise to check whether there's a well-known method which I'm missing.


问:
Getting started with FEM: Ill-conditioned matrix when evaluating flux terms in conservation law?
I have a system of conservation laws of the form

$$ \frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathbf{F}\!\left(\mathbf{q}\right) = 0 $$

I want to use finite elements to solve this system. As an initial choice, I am using Legendre polynomials $\phi$ for my function basis. Let $j$ be a polynomial index. If I substitute the decomposition of my unknown function $\mathbf{q}$ and the flux function $\mathbf{F}$,

$$ \mathbf{q}\!\left(t, x\right) = \sum\limits_j \hat{\mathbf{q}}_j\!\left(t\right) \phi_j\!\left(x\right) $$

$$ \mathbf{F}\!\left(\mathbf{q}\right) = \sum\limits_j \hat{\mathbf{F}}\!\left(t\right) \phi_j\!\left(x\right) $$

and finally write out the weak form (given some test volume $\Omega$ with boundary $\partial \Omega$)

$$ \sum\limits_j \left(\frac{\partial \hat{\mathbf{q}}_j}{\partial t} \int\limits_\Omega \phi_i \phi_j + \hat{\mathbf{h}}_j \oint\limits_{\partial \Omega} \phi_i \phi_j - \hat{\mathbf{F}}_j \cdot \int\limits_\Omega \phi_j \nabla \phi_i\right) = 0 $$

Here, I replaced $\mathbf{F} \cdot \hat{\mathbf{n}} = \mathbf{h}$.

My problem occurs when I actually try to evaluate the surface second and third terms inside the parentheses. Because my problem is 1-dimensional, the second term amounts to just evaluating $\phi_i \phi_j$ at 1 and -1 and taking the difference ($\hat{\mathbf{n}} = \pm \hat{\mathbf{x}}$ is the outward-facing normal). In analogy with the stiffness matrix, we could write this as a matrix $\mathbb{B} = \{b_{ij}\}$. Any given element $b_{ij}$ will be either 0 or 2 because the Legendre polynomials take the value of 1 or -1 at the edges of the domain. This matrix is sparse but singular. The third term exhibits a similar problem if we try to write it as a matrix $\mathbb{C} = \{c_{ij}\}$.

The only way I know to solve for the $\hat{\mathbf{F}}_j$ and $\hat{\mathbf{h}}_j$ coefficients is by solving a matrix equation, but all my attempts to solve them (using PETSc) converge slowly and may not even be correct. I would like to use a preconditioner, but I get errors any time I try (PCSETUP_FAILED due to SUBPC_ERROR), even though I have specified the null spaces of these matrices.

Have I approached this problem incorrectly altogether? Or is there a common way of overcoming this challenge, in PETSc or otherwise?


问:
Efficiently solve linear system with matrix quadratic form
Take the system $$A^TCAx=b$$
where $$A\in\mathbb{R}^{n\times m},\;C\in\mathbb{R}^{n\times n},\;x,b\in\mathbb{R}^{m}, \;m\leq n$$ and $$A^TA=I$$ and $$Cy=d$$ can be solved efficiently in general (specifically $C$ is SPD and it/it's inverse already has Cholesky decomposition). Clearly, if $m\ll n$, then it probably best to just generate and factor the new matrix, but if $m$ and $n$ are both large and $A,C$ are sparse/structured, generating this (typically dense) matrix may not be possible. While I could just use iterative techniques (CG, etc.) at this point, I was curious if there is any to use the knowledge of the inverse of $C$ to develop a more direct (and potentially more efficient) method? </p>

<p>For more context, $A,C$ are Kronecker products in my particular problem. </p>


问:
Numerical solution to Time-dependent Schrodinger equation with time dependent hamiltonian
Currently I am facing the problem to solve numerically the following equation for a double well harmonic potential:

$iℏ\frac{\partial}{\partial t}ψ(x,t)=−\frac{ℏ}{2m}\frac{\partial ^2}{∂x^2}ψ(x,t)+V(x,t)ψ(x,t)$

where

$V(x,t)=\frac{mω^2_0}{2}[a(t)x^4−b(t)x^2]$

Where both $a(t)$ and $b(t)$ are of the form $\frac{1}{2*\sqrt{const_1}}(const_1-const_2 \cos(\omega t))$. I started to use the split-step Fourier method algorithm and here is my code:

import numpy as np
import matplotlib.pyplot as plt

tfinal=10000#total simulation time
dt=0.001#time step
M=np.round(tfinal/dt)
time_steps=M.astype(int)# number of time steps

b=10.# final x
a=-10.# initial x
N=512# num of pieces

dx=(b-a)/N` # x-space step size
dk=2*np.pi/(b-a)
n=np.linspace(-N/2,N/2 -1,N) # this is gonna be my n's in x_n=a+dx*n

x=n*dx
k=n*dk
t=0

alpha=0.0005
beta=0.0001
acoef=(alpha-beta*np.cos(t))/(2*alpha**0.5)
bcoef=(alpha-beta*np.cos(t))**0.5 / (2*alpha**0.5)
V=0.5*(acoef*x**4 -bcoef*x**2)

u0=np.exp(-(x+4.72)**2 /2)*(1/np.pi)**0.25
u=u0

for i in xrange(1,time_steps):
    u=np.exp(-1j*dt*V/2)*u
    c=np.fft.fftshift(np.fft.fft(u))
    c=np.exp(-1j*dt*k**2 /2)*c
    u=np.fft.ifft(np.fft.fftshift(c))
    u=np.exp(-1j*dt*V/2)*u
    t=t+dt

plt.plot(x,u*np.conjugate(u),'r--')
plt.xlim(-10,10)
plt.ylim(0,0.6)
plt.show
```

However the problem that arises it that it takes a very long time to execute in python. I am not sure what am I doing wrong since I am new in python or if I am implementing correctly the algorithm. I would appreciate any insight or suggestion.


问:
$ A * B $ computation when B is a symmetric matrix in armadillo
Is there any way to multiply a symmetric matrix by a dense one in armadillo(and use the fact that we have a symmetric matrix)? I know about DSYMM Routine in BLAS,but the matrices I'm dealing with are mat type(they are not arrays),So I wonder if there is some thing in armadillo using DSYMM Routine from BLAS .
I couldn't find any thing in armadillo documentation.


问:
Domain for convex perspective function
The perspective of a function $f : \mathbb{R}^n \to \mathbb{R}$ is the function $g: \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}$ where $g$ is defined as $$g(x,t) = tf(x/t)$$ with $$\mathbf{\text{dom}}\, g = \{(x, t) : x/t \in \mathbf{\text{dom}}\, f, \,t > 0\}.$$

If I'm modeling a convex constraint of the form $g(x,t) \leq 0$ in a convex solver software, how do I ensure the domain constraint $t>0$ is satisfied? I tried $t \geq 0$, but for some problems, my $t^* = 0$.

Title: Domain for convex perspective function