问:
High-order continuous Galerkin Navier-Stokes solver in deal.II
I looked at the list of publications related to deal.II, but it seems that when it comes to Navier-Stokes equations, everyone speaks about Discontinuous Galerkin solver in their papers. I know there's also a Continuous Galerkin method implemented in deal.II. Has it been used for high-order accurate solution of Navier-Stokes equations? Could someone comment on performance/accuracy differences between CG and DG for fluid flow problems in deal.II?</p>


问:
Sparse matrix factorization of a rank deficient matrix by decomposition into linearly independent components
I've got a little conjecture I need to prove for a theoretical result related to causal Bayes net search with latent variables under sparsity constraints.  If you're interested in the application please feel free to email me at adambrod@andrew.cmu.edu.

Here's the problem:

Suppose we have a $p\times q$ real matrix $A$ such that 1) $A$ is rank deficient so $\textrm{rk}(A)=r<\min(p,q)$, and 2) there are no zero elements in $A$.

I want to prove that there is no factorization of $A$ into two matrices with
fewer total non-zero elements than some factorization of $A$ into two
matrices $B$ and $C$ of dimensions $p\times r$ and $r\times q$.

The proof strategy I'm employing is proof by induction on rank.

The base case, $\textrm{rk}(A)=1$ is trivial ($B$ and $C$ must have at least $p+q$ non-zero
elements between the two of them in order for every element of $A$ to be
non-zero and that is the maximum number of elements in a pair of matrices
with dimensions $p\times r$ and $r\times q$).

The induction step is as follows:  Subtract a matrix $D$ of rank 1 from $A$
such that the resulting matrix $S$ is of rank $r-1$.  This corresponds to
removing all of the "influence" of some specific row vector $v$ from $A$,
such that $v$ is linearly independent of the remainder, $S$.  $D$ can be expressed
as the outer product of some column vector $w$ with $v$.  Suppose $wv$ is the
sparsest factorization of $D$ (the possibility of this is implied by our
base case).  Let $RT$ be the sparsest factorization of $S$ of dimensions $p\times
r-1$ and $r-1\times q$ respectively.  My conjecture is that the sparsest
factorization of $A$ is achieved by appending these two sparsest
factorizations, that is, the sparsest factorization of $A$ is the product of
$R$ with $w$ appended as the $r$th column and $T$ with $v$ appended as the $r$th row.
Because $D$ and $S$ are linearly independent this seems quite plausible, but a
rigorous demonstration eludes me.

Any solutions, suggestions, or references would be greatly appreciated, and of course, if you provide me with personal information (name/department affiliation), I would be happy to acknowledge your assistance in any publications I apply this result in.

Thanks!

ContentLicense: CC BY-SA 3.0


问:
Givens method for sparse matrix
I have a (very) large sparse matrix in CSC form and I'm supposed to factorize it. I've read that between Givens and Householder transformations, Givens is better for a sparse matrix. The problem is that I don't really know how I should do it (everything's fine for a full matrix, but now all I have are three vectors).

I don't want to use built-in Matlab functions.

Title: Givens method for sparse matrix


问:
How can simulation be fast when Electric Field is being measured with probe in CST Microwave Studio?
I am measuring the electric field in the closure box with aperture and using probe to achieve that. But I have a problem. When I use the probe for measurement, results are not accurate in transient solver. So, I changed the settings in the transient solver > Specials and I use the program with high number of pulses option ( For example, 100000 and more pulses ). But simulation time is really bad and it takes nearly four days to complete with my desktop pc ( and I know that my pc is not as fast as enough, but this is my situation). How can I drop the simulation time to measuring the E-Field with probe?

Thank you.


问:
Optimal partition - variable number of parts
Suppose I have a box $D \subset \Bbb{R}^2$ (compact set). Denote $\mathcal{P}= \{ (\Omega_1,...,\Omega_n) : \bigcup_{i=1}^n \Omega_i = D,\ \Omega_i \cap \Omega_j =\emptyset\}$ the family of partitions of $D$ into open sets (the sets can overlap on their boundaries...)

Consider a functional $F$ defined on the subsets of $D$ such that $\lim_{|\Omega|\to 0}F(\Omega)=+\infty$. Suppose that there exists a partition $(\Omega_i)$ such that $\sum F(\Omega_i) < F(D)$ (so that the whole box is not optimal).

Then I consider the problem of finding the partition $(\Omega_i)$ which solves
$$ \min_{(\Omega_i) \in \mathcal{P}}\sum_{i=1}^n F(\Omega_i).$$

The difference between this problem and the problems I've seen and worked until now is that here the number of members of the partition is not fixed.

The first condition tells us that $n$ is bounded above, since many parts means smaller volumes, which mean great $F(\Omega)$. The second condition says that $D$ is not the solution. In what I have in mind $F(\Omega)=G(\Omega)+\alpha \frac{Per(\Omega)}{|\Omega|}$ where $G \geq 0$ is something complicated.

And now the question...

Do you know some references of articles which treat this kind of problems of optimal partitioning with variable number of parts?

Do you know any available programs/toolboxes which are good for manipulating partitioning problems?


问:
Strategy for solving a non-trivial differential equation
I would like to numerically solve an equation of a type as shown below. Does anyone of you have an idea how to approach such a problem? Any links to literature or for further reading would be greatly appreciated.

Thank you!

I search a function $f(\vec{x})$ with $f \in \mathbb{R}$ and $\vec{x} \in \mathbb{R}^n$, where $n \in\{2,3\}$.

$f$ is positive:

$f > 0~ \forall ~\vec{x}$

and radial-symetric with $r=\left|\vec{x}\right|$:

$f(\vec{x}) = f(\left| \vec{x}\right|) = f(r)$

It is bound in infinity:

$\lim_{r\rightarrow \infty} f(r) = 1$

The differential equation looks like:

$\Delta f(\vec{x}) = \nabla \cdot \left[ a \vec{G}(\vec{x})f(\vec{x}) + b f(\vec{x}) \left( \left( \vec{G}(\vec{x})f(\vec{x}) \right) *  f(\vec{x})  \right) \right]$,

where $\Delta$ is the Laplace Operator, $\nabla$ the Nabla operator, both acting on $\vec{x}$. $\cdot$ denotes the dot-product, and $*$ a convolution, with $a,b \in \mathbb{R}$, and $\vec{G}(\vec{x}) \in \mathbb{R}^n$


问:
Adaptive Mesh vs Uniform Mesh for multiple source/boundary/initial Data
I'm going to ask some beginners' questions.

Adaptive mesh can save many DOFs than a uniform mesh. But it also needs to solve linear systems changing with mesh adaptive process. Is this not problem?

Now suppose we are going to solve the problem a lot of times with different multiple source/boundary/initial data. What techniques are useful for adaptive mesh to further reduce costs?

I once heard a lecture given by Hackbusch who proposed to just use a very fine uniform mesh and make use of some low rank property of the matrix to compute its inverse. The number of DOFs sounds crazy but the method seems to work well.

Would someone shares his/her experiences?


问:
Why not use this simpler variant of Stepwise Regression?
In stepwise regression, you step predictor by predictor, each time selecting the one with the greatest correlation with the measurement, subtracting greedily to leave a residual with no correlation to any of the previous predictors.

Why go stepwise? Why not just pick the predictors with best correlation to the measurement right now, however many you need, and use those? Why is a stepwise approach better? Does it yield smaller 1-norms?


问:
Stationary phase approximation for an integral with infinity saddle points?
I need a hand with the numerical evaluation, in Mathematica, for this integral:

$$f(t)=\int_{-\infty}^\infty Exp\{it(\omega_H-\omega_l-\omega_k) - \sum _{j\neq(l,k)} S_j [1-e^{-it\omega_j}]\}\, dt$$

This is the plot of the real and imaginary part:

<img src="https://i.stack.imgur.com/EVtAP.jpg" alt="Real and Imaginary plots of $f(t)$ from $-150 < t < 150$">

And this are the values for the parameters of $f(t)$ employed in the graph:

$$\omega_H=0.408194 ~~;~~ \omega_l= \omega_k =0.102036 ~~;~~\omega_j=0.169499~~;~~S_j=1.05569$$

Where the all parameters are reals; $t$ is the integration time (s), $\omega$ are oscillation frequencies (1/s) and $S_j$ is an adimensional factor.

I find that this type of non-convergent integrals are approximated by the steepest descendent method or stationary phase approximation.

Employing this approximation, the integral can be approximated to:

$$f(t) \approx \sqrt{\left(\frac{2 \pi}{\sum _{j\neq(l,k)} S_j \omega_j^2 e^{-it^*\omega_j}} \right)}  Exp\{it^*(\omega_H-\omega_l-\omega_k) - \sum _{j\neq(l,k)} S_j [1-e^{-it^*\omega_j}]\}$$

Where the saddle points value $t^*$, are determined numerically by the transcendental equation:

$$\omega_H-\omega_l-\omega_k=\sum _{j\neq(l,k)} S_j \omega_j e^{-it^*\omega_j}$$

Solving this transcendental equation in Mathematica, for the parameters above, I get this result:

$$C[1] \in Integers \, \&\& \, t ==(0.- 5.89974 i) (0.131673 + (0.+6.28319 i) \, C[1])$$

This means that the transcendental equation has number of infinity solutions, which also means that this equation has an infinity number of saddles points. Then the unity partition is used to get a convergent approximation to the integral $f(t)$. In this link make an appointment for this case (steepest descendent method)

So, can anyone give me a hand with formation of the cut off functions, using the partition unity, to get a convergent form of $f(t)$?

Thanks!!


问:
Segregated solving of a tightly coupled system of PDEs
To compute the evolution of a free surface between two incompressible, immiscible liquids, two tightly coupled equations have to be solved, the volume fraction advection and the Navier-Stokes equations:

$$\frac{\partial\alpha}{\partial t}=\nabla \cdot (U\alpha)$$
$$\frac{\partial \rho U}{\partial t} + \nabla \cdot (\rho U U)-\nabla\cdot (\mu \nabla U) = -\nabla p + \sigma\kappa\nabla\alpha$$
with $\kappa = -\nabla \cdot \frac{\nabla \alpha}{|\nabla \alpha|}$.

These equations are tightly coupled, and by far the most common approach is to couple them <s>explicitly</s> using a very small time step, i.e. solve for $\alpha$ once, then once for $U$ and then proceed to the next time step. Regardless of the time marching scheme used to solve $\alpha$ and $U$ very small time steps are required to maintain the coupling of the two equations.

Now this is probably a naive question, but why is it not helpful to wrap the solution in an outer loop, relax U and $\alpha$ and solve with larger time steps?


问:
Wanted: sequences of linear systems for recycling Krylov solver analysis
In the solution of sequences of linear systems
$$A_ix_i=b_i\quad\text{for}\quad i=1,2,\dots$$
with Krylov subspace methods, data can be recycled from already solved linear systems in order to speed up the computation of subsequent linear systems (complementary to preconditioning). For the analysis of recycling strategies I developed in my PhD thesis and which are implemented in KryPy, I'm looking for interesting cases of sequences of linear systems that satisfy the following conditions:

The (possibly preconditioned) matrices exhibit only a few eigenvalues close to zero that are separated from the remaining spectrum. Often, the presence of a few eigenvalues of small magnitude is responsible for an initial slow convergence of CG, MINRES or GMRES.

The matrices are sparse.

The problem can be boiled down to small scale (dim<1000).

The matrices and/or the right hand sides of the linear systems may change but the difference should be within limits.

The problem can be implemented in a few lines of code, preferably based on well-known Python packages such as NumPy, SciPy, FEniCS, PyAMG, ...

Symmetric/Hermitian matrices would be nice but are not a must.

Candidates for such sequences are time-dependent problems, nonlinear problems (e.g. with Newton's method), parameter-dependent problems, shifted linear systems, ...

One application that satisfies the above criteria and benefits from recycling Krylov solvers are nonlinear Schrödinger equations, see this paper on arxiv. I'm now looking for a simpler model problem in order to demonstrate the mathematical idea without being computationally cumbersome.


问:
Center of Mass for molecule in Gaussian09
Gaussian09 prints coordinates where the local coordinate system as center at

$C = (\Sigma n_i*X_i)/(\Sigma n_i)$

where $n_i$ is the Mass Numer of the atom. However, the Mass Number is not actually the mass of the nucleus. If we instead choose $m_i$, the true mass of each atom, as the weighting factor, we do not get the origin as the center of mass.

Is it ok to choose the origin that Gaussian09 gives us as the center of mass? Why or why not?


问:
Interface Formulation at Finite Volume Boundaries when using the Dual Mesh
When using the dual mesh (vertex-centered) for finite volume methods, you end up with a cell center at the boundaries between materials.  It is possible that the equations being solved in each material are different, or at the very least don't have continuous derivatives at the boundary (for example: poisson's equation for electric potential moving from gas to metal).

My question is, how do I handle these cells.  Do I split the cell in half so each partial cell is entirely contained in one material or do I treat the boundaries of the cell as the material boundaries (i.e. purely gas cells see the cell boundaries as metal and metal cells see the cell boundaries as gas)?  If I split the cell, do I move the cell center for each off of the new boundary?  Any book or reference handling this would be greatly appreciated, since all the references I have found for dual mesh FVM limit themselves to dirichlet conditions, and don't hanlde interfaces at all.


问:
Keplerian Orbit Estimator
While calculating the eccentric anomaly of a satellite using the Newton-Raphson method, what initial estimate of the root for E (eccentric anomaly) should be taken for least number of iterations(best approximation of root) to take place and why?</p>

Title: Keplerian Orbit Estimator


问:
finite volume for diffusion equation with anisotropic (tensor) coefficient
Consider the scalar  PDE for $u$ with Dirichlet boundary conditions:

$\mathrm{div}(\mathcal{K}\nabla u) = f\; \forall x\;  \in \Omega \subset R^2$,

$u = 0 \; \forall \; x\;\in \partial\Omega$

where $\mathcal{K} \in R^{2\times 2}$ is positive definite symmetric.

For a start I assume $\Omega$ is the unit square with a uniform rectangular mesh.

EDIT: I edited the question following Jed's comment to make it more specific.

If $\mathcal{K}$ is diagonal, it is simple enough. When it is not, I see that there it will require  a computation of the "secondary gradient" (gradient along the face). How is that usually done ? And what I should do when I am near the boundary in this case ?


问:
How to solve singular non symmetric poisson equation with Neumann boundary condtions?
I am trying to solve 2D Poisson equations with Neumann boundary conditions. When the mesh is uniform, Poisson equation is singular and symmetric, so the method listed in Null Space Projection for Singular Systems works well. But when the mesh is non uniform and sharp immersed boundary method is implemented to treat the curved boundary, discretized Poisson equation is singular and asymmetric based on finite difference method. In this case, the pervious method seems not working well. So can someone help me and tell me how to solve this kind of problem, or just provide me some references? Thanks a lot!


问:
CAS to perform "simple" tensor operations in index notation
Heres a couple of questions,

1) Does a CAS exist that can manipulate Cartesian based tensors and produce a simplified output that is hopefully in index notation?
2) If not, is there some indirect method to accomplish this?

I'm talking about doing rather simple manipulations such as computing,

$$
\nabla^{2}_{\mathbf{r}}\left(\frac{a}{r}\left[\mathbb{I} + \frac{\mathbf{r}\mathbf{r}}{r^{2}}\right] + \left(\frac{a}{r}\right)^{3}\left[\mathbb{I} - 3 \frac{\mathbf{r}\mathbf{r}}{r^{2}}\right]\right)
$$

in a 3D Cartesian basis where $\mathbf{r} = \langle x,y,z\rangle$, $\mathbb{I}$ is the identity 3x3 matrix, $r = |\mathbf{r}|$ etc...

What I've tried:
I've tried using Maxima to do this the "pedestrian" way (i.e. making vectors and explicitly computing matricies etc...).  This just results in a complex mess of a big 3x3 matrix with long expressions.  This "works" but the main problem is that simplification is not easy in this case (or maybe I'm not doing it right).

I'm aware of other "Tensor software" (e.g. http://en.wikipedia.org/wiki/Tensor_software) that does these sorts of things but they all seem to be made for more general or complex uses and its not clear to me how to control them for my use.

So I ask if anyone has any experience with this and could perhaps provide an example where simple manipulations like this in a CAS can be performed.

Thanks

OwnerDisplayName: kmarsh


问:
Beam finite element stiffness matrix from section constitutive matrix
I'm trying to construct the 12 x 12 beam element stiffness matrix from a section constitutive matrix (6 x 6 with shear stiffnesses, axial stiffness, bending stiffnesses and torsional stiffness on the diagonal). I can get a 6 x 6 beam element matrix as would be used in a multibody formulation using the method described here:

Malcolm, David J., and Daniel L. Laird. "Modeling of blades as equivalent beams for aeroelastic analysis." ASME 2003 Wind Energy Symposium. American Society of Mechanical Engineers, 2003.

However, I want to use the 6 x 6 matrix that comes out of that process in a conventional finite element model which needs a 12 x 12 matrix. Is there a simple method of converting between the two? I've found information for a simple stiffness matrix but mine is fully populated as it's for a composite beam with lots of material couplings.

For a two-noded beam element, the element stiffness matrix $K$ relates the generalized forces $f$ to the generalized displacements $u$, i.e., $f = Ku$.
In matrix form
$$
  \begin{bmatrix} f_{x1} \\ f_{y1} \\ f_{z1} \\ m_{x1} \\ m_{y1} \\ m_{z1} \\
        f_{x2} \\ f_{y2} \\ f_{z2} \\ m_{x2} \\ m_{y2} \\ m_{z2}
  \end{bmatrix}
  = \begin{bmatrix} K_{1,1} & K_{1,2} & \dots & K_{1,12} \\
                    K_{1,2} & K_{2,2} & \dots & K_{2,12} \\
                    \vdots  & \vdots  & \ddots & \vdots  \\
                    K_{1,12} & K_{2,12} & \dots & K_{12,12} \\
    \end{bmatrix}
    \begin{bmatrix} u_{x1} \\ u_{y1} \\ u_{z1} \\ \theta_{x1} \\ \theta_{y1} \\ \theta_{z1} \\
        u_{x2} \\ u_{y2} \\ u_{z2} \\ \theta_{x2} \\ \theta_{y2} \\ \theta_{z2}
  \end{bmatrix}
$$
where $f_x, f_y, f_z$ are forces, $m_x, m_y, m_z$ are moments, $u_x, u_y, u_z$ are displacements, and $\theta_x, \theta_y, \theta_z$ are rotations.

The section stiffness $S$ relates the generalized forces at a point in the beam to the generalized strains at that location, i.e., $f = S\varepsilon$. In matrix form,
$$
 \begin{bmatrix} f_{x} \\ f_{y} \\ f_{z} \\ m_{x} \\ m_{y} \\ m_{z}
  \end{bmatrix} = \begin{bmatrix} S_{1,1} & S_{1,2} & \dots & S_{1,6} \\
                    S_{2,1} & S_{2,2} & \dots & S_{2,6} \\
                    \vdots  & \vdots  & \ddots & \vdots  \\
                    S_{6,1} & S_{6,2} & \dots & S_{6,6} \\
    \end{bmatrix}
    \begin{bmatrix} \gamma_x \\ \gamma_y \\ \varepsilon_z \\ \kappa_x \\ \kappa_y \\ \kappa_z \end{bmatrix}
$$
where $\gamma_x,\gamma_y$ are shear strains, $\varepsilon_z$ is a normal strain, and $\kappa_x,\kappa_y,\kappa_z$ are curvatures.  These strains are defined as
$$
 \gamma_x = \frac{\partial u_x}{\partial z} - \theta_y \ , \
 \gamma_y = \frac{\partial u_y}{\partial z} + \theta_x \ , \
 \varepsilon_z = \frac{\partial u_z}{\partial z} \ , \
 \kappa_x = \frac{\partial \theta_x}{\partial z} \ , \
 \kappa_y = \frac{\partial \theta_y}{\partial z} \ , \
 \kappa_z = \frac{\partial \theta_z}{\partial z} \ .
$$

Integration of the above relations and then evaluation at the two nodes will probably lead to the beam stiffness matrix, but it is not clear how that should be done.


问:
SciDB vs OLAP cube
I am working on good retrieval mechanism for big data systems. I am of view that multidimensional structures could be of great use for aggregations and querying .  Multidimensional structures  as mentioned in theory/products could replace the need of OLAP CUBE.

Here is some data size info residing into data-warehouse and it is approx. around 1 trillion per month and some time even more and having rolling window of 1 month or greater, it depends per client. I used to write reports,sub-reports and calculate aggregates using sub-reports , and also used analytical and analytical aggregate function to accomplish the reports/query.

Can some one share the thoughts on it.

Regards,

Title: SciDB vs OLAP cube


问:
Inverse problem with a rank-1 update
I hope you can help me out with this. I have to find the solution x to an inverse system

$$
x=A^{-1}b
$$

This inverse problem is basically a least square problem with a rank-1 update.

$$
x=[uv^{T}+H^{H}H]^{-1}H^{H}b
$$

There are two objectives, one can be achieved with the least square problem and the other with the rank-1 update. I need to find the best tradeoff between the two.

The least square solution (without the rank-1 update) needs to be regularized, so the smallest singular values of the matrix $[H^{H}H]$ would need to be replaced by zeroes.

Normally, I should introduce a parameter (real positive scalar) in front of the rank-1 update $uv^{T}$ and hand-tune it until the solution leads to the best compromise between the two objectives. Nevertheless, I am thinking: could I just do a singular value decomposition of the matrix $[H^{H}H]$ and replace one of my useless elements (one that will be replaced by zero after regularization) by my rank-1 update?* Wouldn't that save me from having to include a parameter that has to be tuned? Once I have my new matrix A, I think the solution could be easily regularized by applying Selective Singular Value Decomposition...

The reason why I think this might be done is that there is only one non-zero singular value in $uv^{T}$ and, at the same time, some (the smallest) singular values of $[H^{H}H]$ will be rejected, so why not create a new matrix A such that one of its svd components is the (non-zero) $svd(uv^{T})$?

Thanks in advance

Title: Inverse problem with a rank-1 update


问:
Error analysis and the Model Problem
In numerical methods for ODE's, the model problem

```<code>y' = cy
</code>```

where c is complex is regarded as sufficient in performing error analysis for different methods in solving ODE's.

Is it sufficient, and why? In my opinion, it overlooks many potential functions that can have different error forms.

Title: Error analysis and the Model Problem


问:
What's the most efficient way to calculate the Wiger quasiprobability distribution?
I want to calculate the Wigner quasiprobability distribution function of a particular wavefunction. The definition suggests a few straightforward ways of calculating it, but I was wondering if there's something more efficient than the naive approaches of either directly evaluating the integral or taking the fft of the autocorrelation. </p>


问:
Mass-conservative reprojection (on a sphere)
I have a 2D distribution of mass on a sphere given as a matrix of masses in latitude-longitude grid cells. I need these masses projected to another grid
on the same sphere with different location of pole and different cell sizes, again as masses in grid cells.

The projection has to be:

fast

monotonous, i.e. not create substantial artificial maxima/minima of concentration (mass/cell area)

strictly mass-conservative (to numerical accuracy)

I could not find any ready made tool for that, though
I believe, there should be something available (C, fortran or whatever).
Any clue?

Thank you!


问:
Convolution of two radially symmetric functions for double logarithmic plot

Question: I have problems with sidelobes when computing the convolution of two radially symmetric functions. I wonder whether I should just switch from second order to fourth order and refine the grid aggressively, or whether there are more appropriate methods for computing this convolution.

I have radially symmetric convolution kernels, given in sufficiently well specified representations like sum of Gaussians, piecewise constant, or piecewise polynomial in $r^2$. I have routines to analytically evaluate the Fourier transform of these kernels either pointwise or integrated over a ring. See appendix if interested in these details.

To compute the convolution of two such kernels, I evaluate the Fourier transform pointwise on a sufficiently fine grid in Fourier space, multiply the two Fourier transforms pointwise with each other. I use a logarithmically spaced grid (a factor of 10 is sampled by 200 points), with an additional point at the origin. Then I use some interpolation scheme (currently just linear interpolation in $\rho^2$, but I wonder whether I should switch to a cubic spline interpolation) to get a piecewise polynomial in $\rho^2$ representation of the product in Fourier space. Then I evaluate the Fourier transform of this representation integrated over a ring to get a piecewise constant approximation of the convolution of the two kernels.

My problem is that a double logarithmic plot of the result can reveal spectral leakage in the form of side-lobes. What I don't like about this is that I can make the grid in Fourier space finer to reduce the amplitude of the side-lobes, but a double logarithmic plot can easily cover 12 orders of magnitude, and it seems excessive to me to compute the convolution to an accuracy of $10^{-12}$, just because of artifacts in a double logarithmic plot which are irrelevant to the accuracy of relevant results from the application of the resulting convolution kernel. On the other hand, the artifacts might be more dangerous than I realize, because my kernels are all expected to be positive, but the resulting convolution kernel might have negative values due to the artifacts from the spectral leakage.

As an example, I computed the convolution of a sum of a 2 nm wide Gaussian and a 10 um wide Gaussian with two different very small (0.35 nm and 1 nm) Gaussians. In the images below, the red line is this sum these two Gaussians. The result should be nearly identical to this, i.e. a sum of two Gaussians, where the smaller is 2.03 nm or 2.236 nm wide.

Image of sidelobes after convolution with a 0.35 nm wide Gaussian:

Image of sidelobes after convolution with a 1 nm wide Gaussian:

Appendix
Here is some more of the math (equations) involved: The Fourier transform of $\exp(-r^2)$ is $\pi\exp(-\pi^2\rho^2)$, the Fourier transform of $\operatorname{circ}(r)$ is $\frac{J_1(2\pi\rho)}{\rho}$, and the Fourier transform of $(1-r^2)^n\operatorname{circ}(r)$ isn't much more complicated either. The antiderivative of $\rho$ times these functions can also be computed analytically. Hence it should be clear that it is possible to analytically evaluate the Fourier transform of these kernels either pointwise or integrated over a ring.


问:
Solving "virus", a flash game, in the least amount of attempts
We am currently working on a project related to computational optimization, in which I have tasked myself with calculating the most efficient solution to one of the game's levels.

Link to the flash game: http://trackmill.com/virus-2/

The hexagon grid is randomized in color, and the player starts out in one of the corners of the grid. The player selects a color from a short set of predefined colors corresponding to those present in the grid. When the color of the starting hexagon changes by this action, every neighboring hexagon with the same color now also changes, when the player selects a new color. In this fashion, the player must make the entire grid the same color, within a limited number of moves.

For simplification, we will use squares instead of hexagons, and use only three colors. The algorithm should be able to calculate order in which we should select colors resulting in a win in the least amount of turns.

So far, we are trying to calculate the most efficient solution by way of brute force, testing every single imaginable combination, but how can we do this smarter? Which algorithms can we use / adapt?


问:
Why does my Finite Difference approximation not work?
I am trying to find out the magnitude of the acceleration of my object based on non-uniformly sampled 3D position data. I'm using the standard approximation of the 2nd order derivative on a non-uniform grid (see e.g. here, page 121):

$$
f'' \approx 2 \cdot\frac{h^+ f_0 - (h^- + h^+) f_1 + h^- f_2}{h^+ h^- (h^+ + h^-)}
$$

where $h^-$ is the first sample interval, $h^+$ is the second sample interval, and $f_n$ are the sampled points.

I'm trying to implement it in the following way:

```<code>class ComputeAcceleration
{
    Vector3D prevPosition;
    Vector3D prevPrevPosition;
    float prevTimeStep;

    float acceleration;

    // This function gets called every time step with the time step
    // since the last data input as argument.
    void Update(const float timeStep)
    {
        float doubleTimeStep = prevTimeStep + timeStep;
        Vector3D currentPosition = getCurrentPosition();
        Vector3D accelerationVector = prevTimeStep * currentPosition -
                                      doubleTimeStep * prevPosition +
                                      timeStep * prevPrevPosition;
        accelerationVector  *= 2 / (doubleTimeStep * prevTimeStep * timeStep);
        acceleration = sqrt(accelerationVector.X * accelerationVector.X +
                            accelerationVector.Y * accelerationVector.Y+
                            accelerationVector.Z * accelerationVector.Z);

        prevPrevPosition = prevPosition;
        prevPosition = currentPosition;
    }
}
</code>```

I can't find the error in this, but when I plot this in real-time, I get a curve that doesn't resemble acceleration at all, but rather speed, in some jerky way. I plotted it against the speed (which I know to be correct), and it looks roughly like this:

Am I doing something fundamentally wrong here?


问:
Implementing the transition matrix for page rank
I'm trying to implement PageRank. I'm reading the description here: http://nlp.stanford.edu/IR-book/html/htmledition/markov-chains-1.html

Everything is very clear to me, however I'm concerned about the construction of the matrix $P$. I find that constructing $P$ the naive way would be very expensive. For example: to implement step 1, one would need to check every row of $A$ and then check every element of that row to see if all elements are zero. For step 2 one would need to compute the number of ones for each row. I can imagine my code to have nasty slow loops. I was wondering if there are smart linear algebra techniques that could efficiently construct $P$. I will be using python numpy for my coding.

EDIT: one way I'm thinking now to solve this is by doing a summation element wise over the columns of $A$. By that I would have a column vector. Now I will go through each element of this vector to check which elements are zeros. Thus I can now know which rows has no 1s and I can multiply those rows with $1/N$.


问:
PageRank computation as the left principal eigenvector
I'm following the book here to compute PageRank: http://nlp.stanford.edu/IR-book/html/htmledition/markov-chains-1.html

I built the transition matrix P as indicated in the book. Then I computed the PageRank scores by raising the transition matrix P to higher powers, as said in the book. To make sure my computations are correct I followed this example: http://nlp.stanford.edu/IR-book/html/htmledition/the-pagerank-computation-1.html

The transition matrix is:

```<code>P = np.array([[1./6, 2./3, 1./6],
              [5./12, 1./6, 5./12],
              [1./6, 2./3, 1./6]])
</code>```

The results of the PageRank scores are the following (by computing P^40):

```<code>[0.28, 0.44, 0.28]
</code>```

as in the book.

But the book also says that one can find the PageRank scores by computing the principal left eigenvector of the transition matrix P. So I tried this with numpy:

```<code>v, V = np.linalg.eig(P.T)
left_vec = V[:, 0].T
print left_vec
</code>```

and I got:

```<code>[-0.46829291 -0.74926865 -0.46829291]
</code>```

so what's wrong?!


问:
Finite difference mixed derivatives on quadrilateral (non-orthogonal) grids
Could someone tell, please, has someone constructed finite difference mixed derivatives on quadrilateral (non-orthogonal) grids, already? Without using aconformal mapping to an orthogonal grid. I mean we have the mixed differential operator and try to represent it through a difference relation.

I searched with the help of Google Scholar, but found nothing.

I tried to constuct it myself, but have problems with the symmetry of mixed second difference derivative: $$ f_{xy} = f_{yx} $$


问:
boundary conditions with non-constant coefficients in cell centred finite volume method

Suppose am solving the heat conduction equation in 1d with
Dirichlet boundary conditions. The thermal conductivity $k$ is a non
constant function.
So $-(k(x)u'(x))' = f(x)$

The value of $k$ enters the discretization only through its values at
the faces of the control volume $k_{i+1/2}, \; k_{i-1/2}$, and not at the cell centres $k_i$.

But most books and notes available online suggest that we store the
value of $k$ at cell centres a priori and use it in the computation.

If $k(x)$ has jumps  across the faces, I understand one uses the harmonic average.  What about control volumes near the boundary ? What is the  standard way to approximate $k(x)$ there; harmonic averaging with neighbour is not possible. Can I simply use the cell centre value for these cells ?

I am unable to see how approximating $k(x)$ will eventually affect the order of
accuracy of the method in these  cases.

My question is about accuracy, not exactly the same as How should non-constant coefficients be treated with finite-volume first order upwind scheme?,

EDIT:
When $k(x)$ is smooth, by expanding  the terms $k_{i+1/2}, \; k_{i-1/2}$ in a Taylor's series around $k_i$ in the standard 3 point scheme, I see that using the cell centre  value directly is atleast $O(h^2)$ or second order accurate. But that would violate the conservation property of the scheme since $k_{i+1/2}^+ \neq k_{i+1/2}^-$. But we can  use linear interpolation  with $k_{i}, \; k_{i+1}$  to calculate $k_{i+1/2}$ and maintain $O(h^2)$ accuracy.

But in case of discontinuous $k$, I am still puzzled by order of accuracy of the  harmonic mean approximation for $k_{i+1/2}, \; k_{i-1/2}$, and how to  calculate $k_{i+1/2}, \; k_{i-1/2}$ at the boundary  cells.


问:
One cancels the other qsub semantics atop SGE?
Is there a way to achieve "one cancels the other" semantics (http://www.investopedia.com/terms/o/oco.asp) for qsub jobs submitted to Sun Grid Engine? That is, I submit two (or more) jobs and when one of them runs it cancels the others?

Two examples...

Say I can run for 4 hours on N ranks or 2 hours on 2N ranks. If both jobs can start soon, I'd rather get the answer 2 hours earlier by running the bigger job. If I'm going to wait more than 2 additional hours to run the 2N-rank job, however, it would make sense to run the smaller job first and cancel the bigger one.

Say I can run part of a long running problem for 8 or 16 hours. There's some startup cost amortization for running longer batch jobs so 16 hours is a mildly better use of SUs. But if the 8 hour job can start sooner, I'd rather make progress now than wait for a 16 hour window.

Possibly another way to accomplish these things would be for me to submit a single job with a range of acceptable node counts (hypothetically, "-pe 12way N,2N") or wall times ("-l h_rt 16-8:00:00"). I haven't been able to tease this out of the SGE man pages, however.


问:
A Question About Weak Forms in Fenics
Is it possible to use test and trial functions from two different function spaces (defined over two different meshes) in a single weak form? Under what conditions can I do this (eg., each term in the form has functions belonging to only one of the spaces)?

Essentially, I need this to model heat transfer in a block with a channel passing through it - the block is solid and the channel has water flowing through it. The solid domain and fluid domain are two different meshes. Basically I need to implement two equations like

$$
\int_L \mathbf{u}.\nabla T_l \,\tilde{T_l} + \alpha_l\int_L \nabla T_l. \nabla \tilde{T_l} + \frac{\alpha_l}{k_l} \int_{\Gamma_w} h(T_l - T_s) \tilde{T_l} = 0
$$

$$
\alpha_s\int_S \nabla T_s. \nabla \tilde{T_s} + \frac{\alpha_s}{k_s} \int_{\Gamma_w} h(T_s - T_l) \tilde{T_s} - \alpha_s \int_{\Gamma_1} Q \tilde{T_s} = 0
$$
where $T_l$ and $T_s$ are liquid and solid temperatures respectively, $\mathbf{u}$ is the liquid velocity (for which I'm solving a Stokes equation seperately), $L$ and $S$ are liquid and solid domains respectively, $\Gamma_w$ is the wall between the solid and liquid domains and $\Gamma_1$ is a surface of the solid domain where some heat $Q$ is supplied.
The water inlet temperature is fixed (a Dirichlet BC).

To have a steady state solution, I need to add these equations to get one. Can anyone tell me how I can go about implementing this in Fenics?

Title: A Question About Weak Forms in Fenics


问:
Algorithms for adding hydrogens on a molecule
I want to add hydrogens to some linear polymer molecules (polyethylenes).

I know some working methods like using PyMOL internal function h_add. This method works, but hydrogens are added at distances from the backbone carbon atoms that vary from 1.07 to 1.1 Ångstroms. I need to control the hydrogen distances from the backbone atoms and also the HCH and CCH angles to constant values that I choose every time (for example, H–C = 1.5 Å, ∠HCH = 108.5°, ∠CCH = 109.5°).

Is there a geometry algorithm which I can use to understand the process of hydrogen adding? I know that it can be implemented nicely to a Fortran code but is there an easier way to do this, e.g. a script in Python, or better yet, a bash script?

The molecule coordinate files I am using are in the PDB format.


问:
Solving the Helmholtz equation numerically
I am trying to solve the Helmholtz equation in some complex unbounded 2D domain:

$$\left( \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} + E \right) u(x, y) = 0$$

However, to start with, I tried doing that in the simplified domain $\Omega = [-\infty, \infty] \times [0, H]$, that is, just a tube, with Dirichlet boundary conditions on the walls: $\psi(x, 0) = 0, \psi(x, H) = 0$.

It is known that solutions of this problem are of form:

$$\psi_{n \rightarrow}(x, y) = \sin(\frac{\pi n}{H} y) e^{i \sqrt{E - (\pi n / H)^2} x}$$

$$\psi_{n \leftarrow}(x, y) = \sin(\frac{\pi n}{H} y) e^{-i \sqrt{E - (\pi n / H)^2} x}$$

, for $n = 1 \dots \infty$.

To solve numerically, I have to bound the domain somehow, so I restricted it to $[0, L] \times [0, H]$ for some big $L$.

We want the solution to be unique up to normalization, I am doing that by:

assuming $\left( \frac{\pi}{H} \right)^2 < E < \left( \frac{2 \pi}{H} \right)^2$, that is, getting solution of form $A \psi_{1 \rightarrow}(x, y) + B \psi_{1 \leftarrow}(x, y)$

getting rid of the $B$ coefficient and using the following ansatz in the asymptotic region $x \le 0$: $\psi(x, y) = A \sin(\frac{\pi}{H}y) \sqrt{E - \pi^2/H^2} x)$

So now we should get the solution $\psi_{1 \rightarrow}(x, y) = \sin(\frac{\pi}{H}y) \sqrt{E - \pi^2/H^2} x)$ up to the normalization after solving numerically.

I am substituting the ansatz, getting boundary conditions:

$\psi(0, y) = \sin(\frac{\pi}{H}y)$, (no $x$ dependency due to the freedom in normalization)

$\frac{\partial \psi}{\partial x}(0, y) = \sin(\frac{\pi}{H}y) i \sqrt{E - \pi^2/H^2}$

Next, to solve using finite differences, I used the grid (I will assume cell size is equal to $1$ for simplicity), so $\psi_{n, m}$ now approximates $\psi(n, m)$, and the central differences scheme, rewriting the equation as:

$$\psi_{n, m + 1} + \psi_{n + 1, m} + \psi_{n, m - 1} + \psi_{n - 1, m} - 4 \psi_{n, m} + E \psi_{n, m} = 0$$

I know the values of two first columns of the grid from my boundary conditions:

$$\psi_{0, m} = \sin\left( \frac{\pi}{H}y \right)$$

$$\psi_{1, m} = \sin\left( \frac{\pi}{H}y \right) e^{i \sqrt{E - \pi^2/H^2}}$$

Now that we have two consecutive columns $n - 2$, $n - 1$, we should be able to obtain all the columns up to  $L$, since:

$$\psi_{n, m} = 4 \psi_{n - 1, m} - E \psi_{n - 1, m} - \psi_{n - 2, m} - \psi_{n - 1, m - 1} - \psi_{n - 1, m + 1}$$

However, for some reason, this computation diverges and grows exponentially as $n$ goes from $0$ to $L$. Am I doing something fundamentally wrong (wrong boundary conditions, or I shouldn't use central differences), or that should work?

P.S. I tried the same for 1D Helmholtz in the domain $[-\infty, \infty]$ and $\psi$ diverges as well in this case, so it does seem as I am missing something fundamental.

P.P.S. Turns out that I was just using large step size in the 1D case. However, setting step sizes smaller and smaller in the 2D case does not help :(


问:
About Convex Geometry
A consistency notion in constraint programming:

Let $P = (X, D, C)$ be a CSP.

Given a set of variables $Y \subseteq X$ with $|Y| = k -1$, a locally consistent instantiation $I$ on $Y$ is $k$-consistent iff for any $k$th variable $x_{i_{k}} \in X \setminus Y$ there exists a value $v_{i_{k}} \in D(x_{i_{k}})$ such that $I \cap \{ (x_{i_k}, v_{i_k}) \}$ is locally consistent.
    The CSP $P$ is $k$-consistent iff for any set $Y$ of $k-1$ variables, any locally consistent instantiation on $Y$ is $k$-consistent.

For example, 3-consistency ensures that any instantiation any pair of variables can be extended to an instantiation involving any third variable without violating any constraint. It is equivalent to path consistency. Similarly, 2-consistency is also known as arc consistency.

My question is how convex geometry can ensure global consistency, because it seems that there is a relationship between consistency and convexity.
Because Row Convex Constraint with path(3)-consistency guarantees global consistency.

Can any one help me finding a track to research or a way to prove this relationship?

Title: About Convex Geometry


问:
Why would this a pair of E5-2670 cpus be faster than 4 E5-2640 cpus for large matrix diagonalisation problems?
Not really much more to say;

Of several computers available for use, these are two of the larger ones; one has 2x E5-2670, and the other has 4x E5-4640. The problems we're looking at essentially boil down to large matrix diagonalisation problems (i.e. i have been dealing with 125K rank matrices, and in the past have looked at larger ones. Both sparse and non sparse, using PARPACK).

All intuition would tell me that the Quad CPU 4640 should be much faster, but it's not, and I can't work out why at all. Both have hyperthreading enabled, both have 250GB ram.

Could it be that the gains from the parallelisation of PARPACK beyond 16 cores just start to diminish, and are not worth adding the additional processors due to the increase in inter cpu communcation required?

The codes are written in fortran if that makes much difference.


问:
Integration of nonlinear PIDE via spectral methods
At the mean-field level, the dynamics of a polariton condensate can be described by a type of nonlinear Schrodinger equation (Gross-Pitaevskii-type), for a classical (complex-number) wavefunction $\psi$. Its form in momentum space reads (in 2+1 dimensions):

\begin{multline}
i \frac{d}{dt}\psi(k) =\left[\epsilon(k)
-i\frac{\gamma(k)}{2}\right] \psi(k)
+F_{p}(k)\,\, e^{-i\omega_{p}t}  \\
+  \sum_{q_1,q_2} g_{k,q_1,q_2}\, \psi^*(q_1+q_2-k)
\, \psi(q_1)\, \psi(q_2).
\end{multline}

The function $\epsilon(k)$ is the dispersion of the particles (polaritons). The polaritons are a non-equilibrium system, due to their finite lifetime (damping rate $\gamma$). Therefore, they need continuous pumping with amplitude $F_p$ at energy $\omega_p$.
Finally, there exists a momentum-dependent nonlinear interaction $g_{k,q_1,q_2}$ that depends of the so-called Hopfield coefficients $X$ (simple functions of momentum) as:

\begin{equation}
g_{k,q_1,q_2}=g\, X^*(k)\, X^*(q_1+q_2-k)\, X(q_1)\, X(q_2)
\end{equation}

Equivalently, the real-space correspondent of this equation is:

\begin{multline}
i\frac{d}{dt}\tilde{\psi}(r)=\left[\epsilon(i\partial_{r})-i\frac{\gamma(i\partial_{r})}{2}\right]\tilde{\psi}(r)+F_{p}e^{i(k_{p}r-\omega_{p}t)}\\
+\iiint dr_{123}\,\tilde{g}(r_1-r,r-r_2,r-r_3)\tilde{\psi}^*\!(r_1)\tilde{\psi}(r_2)\tilde{\psi}(r_3)
\end{multline}

with

\begin{multline}
\tilde{g}\left(r_{1}-r,r-r_{2},r-r_{3}\right)=g\int dq_{3}X(q_{3})e^{iq_{3}\left(r-r_{3}\right)}\int dq_{1}X^*(q_{1})e^{iq_{1}\left(r_{1}-r\right)}\\
\times \int dq_{2}\, X^*(q_{2}+q_{3}-q_{1})\, X(q_{2})\, e^{iq_{2}\left(r-r_{2}\right)}
\end{multline}

where i have used that

\begin{equation}
\psi(k)=\int dr e^{-ikr}\tilde{\psi}(r)
\end{equation}

My question is how to best solve this equation numerically. Boundary conditions are not important, one can use PBC if that makes it easier.

In fact, I was thinking that it would maybe be easier to solve the equation in k-space, wouldn't that be equivalent to using a spectral method? I am not sure how to do it in practice though, and especially how to take care of the nonlinear interaction term containing the 2 sums over $q_1$ and $q_2$.


问:
Computing a description of a set of lattice points
I have an infinite subset $S$ of $\mathbb{N}^k$ which I know (or at least strongly suspect) is a finite union of regions, each of which is defined by finitely many linear equalities and inequalities (with integer coefficients). I would like to find an explicit description of $S$, i.e., sets of (in)equalities that define it. Is there a good algorithm for this?

A little more information. The representation I have for $S$ is its indicator function, i.e., I have a function $f \colon \mathbb{N}^k \to \{0,1\}$ such that $f(x)=1$ if and only if $x \in S$. In my case, $k \approx 10$ and $f$ can be computed reasonably efficiently. Obviously, in this generality one cannot verify in a finite computation that a hypothetical description of $S$ is correct. However, for my $S$ I expect the answer to be reasonable, and I have other means to prove a purported description is correct.


问:
advection diffusion equation
Would you know what is the condition for stability for the advection-diffusion equation where we treat the diffusion part using Crank-Nicholson and the advection part using FCTS (forward in time centered in space)? I am applying von Neumann analysis but I am not sure about the final condition for stability. Do you know where I could find the proof for the stability of that scheme? Thank you</p>

Title: advection diffusion equation


问:
How does one handle the source term in the Shallow Water Equations when using the discontinuous galerkin method?
I use the discontinuous galerkin method to solve the steady flow 1D shallow water equations with a bump at the bottom. This flow is frictionless.

I use the runge-kutta method to approximate the time derivative of these equations and gauss-lobatto quadrature to approximate the integral equation in the weak form of the source term.

The solution that I get always blows up and doesn't correspond with the exact solution, I think that I may be approximating the source term incorrectly but I don't know how to fix this problem.

Any suggestions?


问:
learning Gurobi for MIP calculation
I am a new learner of Unix, python and Gurobi for MIP calculation. I basically finish reading the quick start guide in Gurobi website but still feel that I need more training to use this combination to automatically write MIP input files. Do anyone have some suggestion on any available trainings (especially in using Unix or Python to automatically generate the input file for Ising model with large grid sizes). Thank you :)</p>

Title: learning Gurobi for MIP calculation


问:
Polynomial Fitting with Least Squares using Numpy and Scipy
I am trying to fit data to a polynomial using Python - Numpy. The points, with lines sketched above them are as in the picture.

I am trying to fit those points to a polynomial of 4. or 5. degree. However, all I can get is nothing more than a line. The coefficients of other than linear terms are too small. I tried to do that both with Numpy and Scipy. In Numpy, I have used the polyfit function and followed the example given in this link and in Scipy used the optimization.curve_fit function and followed the page in the comment below. (I cannot post more than two links!).

My question is: why do I get a linear fit for such data? I provide high enough data (about 70 points). And how can I improve my fit?

Any suggestions are highly appreciated!


问:
Estimating the parameters of the DACE stochastic model (EGO optimization algorithm)
Good day.

I am trying to implement the EGO optimization algorithm. The algorithm itself is rather long to describe here in full. It is presented here, with another example of usage here.

But before I implement the algorithm, I have decided to test the DACE model on a simple 1-d function.

I am receiving strange results; I believe I have located the problem, but am uncertain, whether my implementation is wrong, or if what I found is a weak point in the model.

To initiate the algorithm, we need first to determine a set of parameters $\theta = \{\theta_h\}$ and $p = \{p_h\}$ from some sample data we measure.

The way we do so is as follows:

Assume we have $X = \{X^1, ..., X^n\}$ a set of points we have evaluated the function at, and $Y = \{y_1, ..., y_n\}$ the values at those points. Denote $X^i = (x^i_1, ... x^i_k)$.

Then define
$$ d_{ij} = \sum_{h=1}^k \theta_h |x_h^i - x_h^j|^{p_h}$$
$$R_{ij} = e^{-d_{ij}}$$

We estimate the values $\theta, p$ by maximizing the likelihood function
$$\frac1{(2\pi)^{n/2}(\sigma^2)^{n/2}|R|^{1/2}} e^{-\frac{(...)R^{-1}(...)}{\sigma^2}}$$

(The exact expression in the exponent does not matter for what I am about to say, but is present in the files linked to in the beginning. Also, $\sigma^2 = \sigma^2(\theta,p)$).

I carry the maximization of the likelihood by the Nelder-Meads algorithm (fminsearch in MATLAB). But here is my problem: in some cases, when the initial data set is too large (I am minimizing $-sin(x)$ on $[0,\pi]$, and 'too large' is ~10) what happens, is that the values the algorithm goes for are such, that the matrix $R$ becomes nearly singular ($cond(R)$ ~ $1e^{16}$). This can be explained, since we have a determinant of the matrix in negative power in the likelihood function. Thus, for a singular matrix the determinant becomes zero, and the likelihood function itself diverges to infinity. But because my calculations also feature $R^{-1}$, having a nearly-singular matrix is problematic for the algorithm (and that is indeed the case, the returned estimated function behaves nothing like the original).

Note, that this does not happen always, and for some initial data of ~4 points, the aprroximation to the function behaves very well, so I have reason to believe that the mistake is not simply a coding mistake.

Thus I ask again: is there another way to estimate the parameters? Am I doing something wrong? Is this a weak point in the algorithm?


问:
Algorithm to explore a 3D scalar potential
I have to write an algorithm/code to reconstruct an unknown 3D scalar potential $ V(\mathbf{r})$ up to a specified threshold $E_c$ (assume that the potential is monotonic)

The problem has the following rules:

You start from the minimum of the potential and define a grid

you can only do "small" steps, i.e. you can only move to the neighboring points with respect to where you are (random sampling is not possible)

at the end of the procedure, all the points $\mathbf{r}_i$ of the grid having $V(\mathbf{r}_i)<E_c$ must be known.

Is the a algorithm or a keyword that I can look for in google for this kind of problem?

Thanks


问:
How should I choose the knot sequence when using B-splines as a basis for solving a PDE?
I'm looking to solve the Schrödinger equation with a basis made of a tensor product of basis splines. A number of papers describe calculations made with a program designed this way, but they never describe how they got the knot sequence they use to define the basis.

How does one go about choosing a knot sequence for a set of B-splines to be used as a basis for the solution of a PDE?


问:
Elliptical problem with Robin BC
Working in finite differences, I am using a transformation on the temperature variable

$\Theta = \int_{T0}^T \kappa(T)dT$

to linearize the steady-state heat equation into a Poisson equation

$-\nabla\cdot(\kappa\nabla T)=-S \rightarrow \nabla^2 \Theta= -S$

The Dirichlet and Neumann boundary conditions transform in a way such that they have the same form.  However, it seems that a Dirichlet condition for a heat sink pulls out too much heat from the system.  So I'm trying a convective boundary condition, but this Robin boundary condition causes a problem

$ -\kappa \nabla T = h(T-T_f) \rightarrow -\nabla \Theta=h[T(\Theta) - T_f] $

where $h$ is a heat transfer coefficient, $T(\Theta)$ is the inverse transform back to the temperature variable, and $T_f$ the ambient temperature outside the boundary.  The inverse transform is nonlinear making the boundary condition nonlinear as well.

I'm just trying to find some ideas on how I could possibly implement this sort of boundary condition in this case.  Can I Just discretize the derivative in the boundary condition as in a Neumann except that instead of a ghost point I have a point outside the domain where the temperature is known ($T_f$) and then solve iteratively?

Title: Elliptical problem with Robin BC


问:
Tecplot Span-Wise Loading Visualization
all.

I have data from a CFD simulation for a wing, and I need to obtain a visualization of the span-wise load distribution. I'm attempting to do this using Tecplot. I've tried using 'slices' and 'vectors', but got stuck. Any ideas on how to do this? Any input from engineers would be greatly appreciated.

Title: Tecplot Span-Wise Loading Visualization


问:
Floating-point arithmetic in scientific computations rules of thumb
I am looking for a nice reference (a review, tutorial, or maybe a book) that has tips and their explanations about general issues of floating-point arithmetic in scientific computations.

Some that I've seen are Numerical Recipes and Goldberg, D. E. (1991). What every computer scientist should know about floating-point arithmetic. ACM Computing Surveys, 23(1), 5.

Numerical Recipes is a very large book with specific algorithms. I am looking for general rules of thumb, of general applicability.

The article by Goldberg has a bit too much emphasis on multiple hardware implementations. I really only care about C/C++ programming on modern PCs, which has a pretty standard way of dealing with floats.


问:
NURBS surface fitting for a closed region on mesh
I'm developing a tool that allows users to select a closed boundary (a polygon) on the triangle mesh and then from this boundary, generate a NURBS surface fitting the original mesh surface.

My idea is as following:

Generate a list of NURBS curves from the closed boundary (polygon)

Combine such curves together using NURBS loft surface technique (Reference: The Nurbs book - Les A. Piegl, Wayne Tiller)

My problem arise at the first step, where I am trying to interpolate a list of NURBS curves, which require an input of sample points. To get such sample points, I am planing to generate a point grid on region inside the closed boundary, and then, for each row of this point grid, I will take sample points from this. But I don't know how to create a point grid on the mesh. Could you help me out?

This is just my idea, and if you know more about NURBS fitting for a region of mesh surface, please suggest me some material about it.

Thanks in advance.


问:
Matlab: why I am getting the block diagonal form of this antisymmetric matrix always as zero?
```<code>    function [R, RI , S ] = Matrix_block (A,N)
[U,D]= schur (A);
E=ordeig(D) ;
[R, S]= ordschur (U,D, abs (E) <1000* eps ) ;
RI=R';
 m= max(max(R* RI-eye (2*N) ) ) ;
if m >100* eps;
 sprintf( 'R not unitary , N: %d , Maximum : %d \n ' ,N,m);
 end
C = R*S*RI-A;
m= max(max( abs (C) ) ) ;
if m >10* eps ;
sprintf ('Rtransformeriet , N: %d , Maximum : %d \n', N,m)
end
</code>```

I am using the above code to get the block diagonal form of a certain antisymmetric matrix with non-zero entries. I am getting the block diagonal form to be always a null matrix. Is this a precision error ?

How do I correct this ? Here is the matrix if needed :- https://drive.google.com/file/d/0BxqrCZfy_mcXLW5PbkVUeHVSODA/edit?usp=sharing
Thanks a lot in advance :-)


问:
Wavelets frame for $L^2[0,\infty)$
I need a wavelet frame for $L^2[0,\infty)$. Moreover, the wavelet should be twice differentiable and with continuous second order derivatives. Hopefully, the wavelet should have compact support (either in the time-domain or in the frequency-domain). Is there such a wavelet?

Previously I had asked similar questions here and here. However I think this question fits well here.

Thanks.

Title: Wavelets frame for $L^2[0,\infty)$


问:
Algorithm for optimizing graph interconnectivity
I have a partiuclar kind of graph problem and (not having a background in graph algorithms) I would like to know how this kind of problem is called in the literature and what algorithms exist for solving it.

I have an undirected graph consisting of vertices $V$ and edges $E$. Each edge has a potential $p_e \in [0..1]$

Independently of the edge potentials there is also edge weights $w_e \in \{0,1\}$ that turn on/off an edge and connect or disconnect two vertices. The number of 1-weight-edges of a vertex let's call $n$, which is the "connectedness" (or "subgraph-degree") of a vertex.

In our problem edge potentials are given and edge weights are what we look for.

Let's then call $s_e=p_e*w_e$ the "strength" of an edge. Also, for any two vertices, if there's an indirect connection the strength of the indirect connection is the product of all strengths $s_e$ along the path. The strength of the shortest connection between two vertices let's call $t_e$.

So, the problem is: for a given global $n$ (same for all vertices) how can I find weights $w_e$ for each edge resulting in the "best" interconnectivity for the graph, the one that maximizes the sum of $t_e$?

In other words: Having a given (usually small) $n$, I would like to choose $n$ out of many existing edges for each vertex such that the resulting graph is "nicely" connected. That means that the algorithm would favor connecting two vertices with a somewhat lower $p_e$ over connecting vertices with a higher $p_e$ but who already have a pretty good indirect connection.

Update: According to the comment below I am looking for the $n$-regular edge-induced subgraph that maximizes the sum of all $t_e$.

I have a hunch that this might be NP hard but any approximate solutions would be appreciated too.


问:
Understanding the use of Minimum Translation Vector (MTV) for an elastic collision response
So, I'm trying to simulate polygon rigid body physics, both linear and angular.

To calculate the MTV is easy in SAT (Minimum Translation Vector), and to use it to adjust position is also easy:

```<code>PolygonA.center -= MTV * .5
PolygonB.center += MTV * .5
</code>```

But if you want to adjust velocity and angular velocity with it how would you do that?

This is what I've thinked about, I haven't tested it yet though in my program:

```<code>AccelerationA = -MTV * (2 * PolygonB.mass / (PolygonA.mass + PolygonB.mass))
AccelerationB = MTV * (2 * PolygonA.mass / (PolygonA.mass + PolygonB.mass))

PolygonA.velocity += AccelerationA
PolygonB.velocity += AccelerationB

PolygonA.angularVelocity = Cross(Normalize(AccelerationA), PolygonA.velocity) / PolygonA.radius
PolygonB.angularVelocity = Cross(Normalize(AccelerationB), PolygonB.velocity) / PolygonB.radius
</code>```

Where radius is the longest distance from the center to the fartherst vertex, the rotation axis length of that polygon.

Now is my idea correct?, and if not, how do I fix it?

OwnerDisplayName: Murplyx


问:
Basis for spatially localized function
I'm solving an ODE, where currently the dependent variables are the (time dependent) spatial Fourier coefficients. It turns out that the phenomena I'm interested in describing is spatially localized, and it takes an impractical number of Fourier coefficients to fully describe it. Therefore, I'm interested in looking at different basis to use, instead of the canonical trig functions $e^{ikx}$ (for $k\in \mathbb{Z}$).

However, for this model I need these basis functions to have certain properties. In particular, I need these functions to be differentiable, as well as periodic. Also, I need the functions to have relatively simple Hilbert Transforms, since this is exploited in solving for the governing equations (and is one reason why the trigonometric basis was desirable). Finally, the basis should converge more rapidly than the Fourier series for spatially compact scenarios. Hence, I would like to reduce the number of dependent variables that I'm solving for in my ODE.

Does anyone know of any functions fitting the bill? Any tips/references are appreciated.

Nick

Title: Basis for spatially localized function


问:
Lattice Boltzmann Method
I have done Molecular Dynamics Simulation and now want to venture into Lattice Boltzmann Method. What would be the best reference book/lecture notes/videos for a beginner?</p>

Title: Lattice Boltzmann Method


问:
Question about ellipsoid method
I have some technical question concerning the ellipsoid method Referring to the paper :

http://paswkshop.comm.utoronto.ca/~weiyu/01658226.pdf

It is mentioned in p.1317 at the last line in the left column that "When the subgradient of the dual objective function has a bounded norm [26], the ellipsoid update is guaranteed to converge"

Now, the question is what if the subgradient doesn't have a bounded norm ? Can convergence to the optimal multiplier value still take place ? In my problem, there is no upper bound on the multiplier value like what the authors did in the attached paper in equation (23).

Title: Question about ellipsoid method


问:
Python - calculation time derivative and laplacien by finite differences
I would like to determine a temporal derivative and a Laplacian by the finite differences method to solve the heat equation in a 1-dimensional case. My aim is to get the sources term that is why I need to get all left terms of the equation.

When I plot the sources term as a function of time, the result is surprising... very far from expected. Is this a problem with the np.reshape?  Thanks for help ;)

```<code>#profil1D_theta.shape[0] : time
#profil1D_theta.shape[1] : space

#profil1D_theta_reshape.shape
#Out[986]: (328, 411)

T_tempo = 1./115. #time period

#Time derivative

dt_theta_list = []

for i in range(profil1D_theta.shape[0]-1):
    for j in range(profil1D_theta.shape[1]-1):
        dt_theta = (profil1D_theta[i+1,j] - profil1D_theta[i,j])/T_tempo
        dt_theta_list.append(dt_theta)

dt_theta = np.asarray(dt_theta_list)
dt_theta_reshape = dt_theta.reshape(328,411)

#Laplacian

T_spatial = 0.000021661 #spatial period

dx2test_theta_list =[]
for i in range(0,profil1D_theta.shape[0]-1):
    for j in range(0,profil1D_theta.shape[1]-1):
        dxtest_theta = (profil1D_theta[i,j+1] - 2*profil1D_theta[i,j] + profil1D_theta[i,j-1])/T_spatial**2
        dx2test_list.append(dxtest_theta)

dx2test_theta = np.asarray(dx2test_theta_list)

dx2_theta_reshape = dx2test_theta.reshape(328,411)

rho = 7850.
C = 500.
k = 42.
D = k/(rho*C)
e = 3.*0.001
b = 6.*0.001
Boltz = 5.679373*10**(-8)
Emissivite_echant = 0.97 #
T0 = 273.15 + 26.
h = 5.

tau_th = rho*C*e/(2*h+8*Boltz*Emissivite_echant*T0**3)

#Sources

profil1D_theta_reshape = profil1D_theta[0:328,0:411]

w1_chaleur_list=[]

for i in range(0,profil1D_theta_reshape.shape[1]):
    for j in range(0,profil1D_theta_reshape.shape[0]):
        w_chaleur1D = rho*C*(dt_theta_reshape[j,i]+1./(tau_th)*profil1D_theta_reshape[j,i] - D*dx_theta_reshape[j,i])
        w1_chaleur_list.append(w_chaleur1D)

d,e = dt_theta_reshape.shape
w1_chaleur = np.asarray(w1_chaleur_list)
w1_chaleur_reshape = np.reshape(w1_chaleur,(e,d))

clim_min = 0
clim_max = 1*10**6

plt.imshow(w1_chaleur_reshape);plt.clim([clim_min,clim_max]);plt.colorbar();
plt.ylabel(r'Axe longitudinal')
plt.xlabel(r'Temps (images)')
plt.show()
</code>```


问:
Fast way to compute barycentric lagrange interpolation
Is there any fast way to compute the barycentric Lagrange interpolation using matlab?
something more faster than using repmat instead of for loops</p>


问:
Generating Hessian of the Lagrangian with dynamic number of nonlinear constraints in Fmincon
I'm using interior point algorithm for solving a nonlinear optimization problem and want to provide Hessian of the Lagrangian as part of fmincon to speed up the process (running couple of thousand variations for different sets of parameters)

The workflow that I'm designing has the (nonlinear) objective function, form of linear and nonlinear constraints each in separate function files. Since i'm running different variations of the problems, sometimes there may be zero nonlinear equality & zero nonlinear inequality constraints and in other scenarios, nonzero number of nonlinear equality/inequality constraints. In short the number of nonlinear constraints are dynamic and stored in the 2 variables "NumOfNonLinEqConstr" and "NumOfNonLinInEqConstr"

I have taken a look at Generating Hessian using Symbolic toolbox and few other web-pages but cannot see an example where the Hessian of the Lagrangian is constructed for dynamic number of constraints. In the referred matlab webpage example, like in one variation I tried replacing 10 with NumOfNonLinInEqConstr bu it doesn't work as matlabFunction does not work on cell data type. Can anybody provide a working example for constructing hessian of the lagrangian with dynamic number of nonlinear constraints

Relevant code lines from the matlab link provided below: My question is how could I avoid hard-coding 10 in the below code (and also manually add the product of hessian of the constraint and lambda.ineqnonlin 10 times)?

```<code>hessc = cell(1, 10);
for i = 1:10
    hessc{i} = jacobian(gradc(:,i),x);
end

for i = 1:10
    ii = num2str(i);
    thename = ['hessc',ii];
    filename = [currdir,thename,'.m'];
    matlabFunction(hessc{i},'file',filename,'vars',{x});
end

function H = hessfinal(X,lambda)
%
% Call the function hessenergy to start
H = hessenergy(X);

% Add the Lagrange multipliers * the constraint Hessians
H = H + hessc1(X) * lambda.ineqnonlin(1);
H = H + hessc2(X) * lambda.ineqnonlin(2);
H = H + hessc3(X) * lambda.ineqnonlin(3);
H = H + hessc4(X) * lambda.ineqnonlin(4);
H = H + hessc5(X) * lambda.ineqnonlin(5);
H = H + hessc6(X) * lambda.ineqnonlin(6);
H = H + hessc7(X) * lambda.ineqnonlin(7);
H = H + hessc8(X) * lambda.ineqnonlin(8);
H = H + hessc9(X) * lambda.ineqnonlin(9);
H = H + hessc10(X) * lambda.ineqnonlin(10);

end
</code>```


问:
Generating start point in a systematic manner for fmincon
I'm trying to generate start points for my optimization problem in Matlab. At this point Im not worried about feasibility but only a fast way to generate the points from which I could test the performance of different algorithms available within fmincon. Though am not worried about feasibility I don't want to use random number generator and would prefer to use a more "systematic" approach as below.

As a simple example, assume I have 2 decision variables x1 and x2 and I can provide the min/max ranges for each of these 2 variables.

Like,

```<code>xlow=[2,9];% 2 represents the min possible value of decision variable x1 and 9 represent the minimum possible value of x2;
xhigh=[5,16];% 5 represents the max possible value of decision variable x1 and so on;
</code>```

One start point that I can use is

```<code>Random_Point1=(x1ow+xhigh)/2;% centre of the rectangle formed by max/min of decision variables;
</code>```

After running the optimization algorithm with above start point, let's assume I want to generate more start points. Now I can generate 4 additional start points based on the mid-point Random_Point1 and the 4 original vertices (another way to think about this is by dividing the original rectangle in to 4 equal size smaller rectangles and finding the centre of each of the rectangle)

```<code>Random_Point21 = (Random_Point1 + [xlow(1),xlow(2)])/2;
Random_Point22 = (Random_Point1 + [xlow(1),xhigh(2)])/2;
Random_Point23 = (Random_Point1 + [xhigh(1),xlow(2))/2;
Random_Point24 = (Random_Point1 + [xhigh(1),xhigh(2)])/2;
</code>```

Now I again run the fmincon with above 4 start points.

After looking at the results, I decide to generate more start points. This time I can generate 16 start points by using the previously generated above start points. The way to think about this is dividing the original rectangle in to 16 equal pieces along the respective axis and finding the centre of each of the 16 rectangles

I would continue the process of generating start points and running fmincon till my solution is deemed satisfactory (I have some criteria on when I want to stop optimization process)

My question is, what approach should I use to generate these start points in the most efficient manner (quick). Whatever I wrote (similar to above hard-coding) is definitely not the right approach. I guess I need to use recursion

Stating the problem in general terms:-

a) Assume I have k decision variables and for each decision variable I have the min/max ranges

b) I want to start with generating the mid-point (easy) and then in first iteration generate 2^k points, in 2nd iteration generate 2^(2k) points, in 3rd iteration generate 2^(3K) points and so on.

Note, in my application the number of decision variable is not expected to be large (may be 3 or maximum 4)..the complexity I have is more in functional form of objective rather than # of decision variables. I'm stating this so that there are no worries about generating unpractical number of start points (like I will apply a condition that if in the pth iteration if 2^((p-1)k) was greater than 10K abort the generation process.

Note: I do have a working model albeit of another approach where for let's say 2 decision variables with given min/max I want 400 start points. I divide each of the individual decision variables to Square root of 400 and for all the resulting combinations I will have 400 points. I would prefer to use alternate approach and one of them is described above.


问:
Issues with Open MPI
Not sure if this question is more appropriate for StackOverload, but I'm guessing a bunch of people here have experience with MPI2, so thought I would try here first. I have some large dense complex symmetric matrices that I would like to factor and wanted to try out Elemental for its MPI support. I compiled Open MPI 1.8.1 and a simple "hello world" example, but for some reason when I try to run: mpirun -n 2 -host localhost,machine-2 hello_world it just exits without outputting anything or even starting the hello_world program. If I run it on just localhost (i.e. mpirun -n 2 -host localhost hello_world), it works just fine. Do you use Open MPI or a different flavor of MPI2? I've searched everywhere online and can't find a working solution or even just a cause for this problem. (I have ssh setup between all the machines such that I do not need to use a password.) The machines are running RedHat 5, and I do not have root access to them. I have read that it might be a firewall problem (what ports does OpenMPI use outside of ssh?). Do I need to have anything placed in /etc/hosts? Thanks so much in advance!

Title: Issues with Open MPI


问:
The solution found by the optimization function when exitflag is negative
I'm using Matlab's fmincon to solve a constrained nonlinear optimization for varying sets of parameters in the objective function. I have 2 parameters and have 400 unique combination of these 2 parameters. Within a particular combination of parameters I use a range of start points (100 start points). Trying to do rough simulation of multi-start functionality using optimization toolbox.

For a given combination of parameters in objective function, majority of the start points would solve the optimization problem and exit with a positive exitflag like 1 and in some cases 2. But for very small number of cases, the program would finish with an exitflag of -3 which corresponds to objective function at current iteration went below options.ObjectiveLimit and maximum constraint violation was less than options.TolCon.. In such cases, the output message mentions about problem might be unbounded. One example message is below

Problem appears unbounded.

fmincon stopped because the objective function value is less than
  the default value of the objective function limit and constraints
  are satisfied to within the default value of the constraint tolerance.

Stopping criteria details:

fmincon stopped because the objective function value, -5.374294e+40, is less than
  options.ObjectiveLimit = -1.000000e+20, and the relative maximum constraint
  violation, 0.000000e+00, is less than options.TolCon = 1.000000e-06.

Optimization Metric                                         Options
  objective value =  -5.37e+40               ObjectiveLimit =  -1e+20 (default)
  relative max(constraint violation) =   0.00e+00    TolCon =   1e-06 (default)

Here are my questions: -

a) If the problem was able to find local minima for certain start points (as indicated by exitflag 1) within a particular combination of parameters, why does few other start points mention about problem being "unbounded" (exitflag -3).

b) When the exitflag is -3, the objective function value returned by fmincon does NOT match the explicit calculation of function value at the x point returned by fmincon. Matlab link on Current Point and Function value mentions "The function value is the value of the objective function at the current point." which does not state restriction on exitflags. On the other hand, Output Arguments states this for x "The solution found by the optimization function. If exitflag > 0, then x is a solution; otherwise, x is the value of the optimization routine when it terminated prematurely." possibly implying that if exitflag is negative then The function value is NOT the value of the objective function at the current point. I'm not sure whether it is inconsistent documentation or a mis-interpretation on my part?

Lastly, will draw attention to the statement Exit flags are not infallible guides to the quality of a solution. Many other factors, such as tolerance settings, can affect whether a solution is satisfactory to you. You are responsible for deciding whether a solver returns a satisfactory answer. Sometimes a negative exit flag does not correspond to a "bad" solution. Similarly, sometimes a positive exit flag does not correspond to a "good" solution.. I can understand +ve exitflags not being necessarily good in all scenarios, but why would a negative exitflag not necessarily mean a "bad" solution?


问:
algorithm Matlab function 'stream2' and 'stream3' use to compute streamline
Can anyone explain how the MATLAB functions: stream2(x,y,z,u,v,startx,starty) and stream3(x,y,z,u,v,startx,starty) compute stream lines from vector data u and v? In the stream2.m file, it calls the stream3c function. I assume it is a .c file but unable to see it. Thanks in advance.</p>


Title: Finite differences and Neumann boundary conditions

Body:
I am dealing with a highly nonlinear system of two PDEs. I already have a code to solve the system in case of Dirichlet boundary conditions. The explicit system is:

$$
\begin{eqnarray*}
\partial_{t}u & = & D_{11}\partial_{x}^{2}u+\partial_{x}D_{11}\cdot\partial_{x}u+D_{12}\partial_{x}^{2}v+\partial_{x}D_{12}\cdot\partial_{x}v\\
\partial_{t}v & = & D_{21}\partial_{x}^{2}u+\partial_{x}D_{21}\cdot\partial_{x}u+D_{22}\partial_{x}^{2}v+\partial_{x}D_{22}\cdot\partial_{x}v
\end{eqnarray*}$$

with boundary conditions:

$$u(0,t)=0 \ (\partial_x v)(0,t) = 0$$

And the discretized form of the scheme are:

$$\begin{eqnarray*}
-D_{21}(x_{j})u_{n+1,j-1}-D_{22}(x_{j})v_{n+1,j-1}\\
{}[\mu+2D_{11}(x_{j})]u_{n+1,j}+2D_{12}(x_{j})v_{n+1,j}\\
-D_{11}(x_{j})u_{n+1,j+1}-D_{12}(x_{j})v_{n+1,j+1} & = & \frac{1}{4}[4\mu-D_{11}(x_{j+1})+D_{11}(x_{j-1})]u_{n,j}\\
 &  & -[D_{12}(x_{j+1})-D_{12}(x_{j})]v_{n,j}\\
 &  & +\frac{1}{4}[D_{11}(x_{j+1})-D_{11}(x_{j-1})]u_{n,j+1}\\
 &  & +\frac{1}{4}[D_{12}(x_{j+1})-D_{12}(x_{j-1})]v_{n,j+1}
\end{eqnarray*}$$

$$
\begin{eqnarray*}
-D_{21}(x_{j})u_{n+1,j-1}-D_{22}(x_{j})v_{n+1,j-1}\\
2D_{21}(x_{j})u_{n+1,j}+[\mu+2D_{22}(x_{j})]v_{n+1,j}\\
-D_{21}(x_{j})u_{n+1,j+1}-D_{22}(x_{j})v_{n+1,j+1} & = & \frac{1}{4}[4\mu-D_{21}(x_{j+1})+D_{21}(x_{j-1})]u_{n,j}\\
 &  & -[D_{22}(x_{j+1})-D_{22}(x_{j})]v_{n,j}\\
 &  & +\frac{1}{4}[D_{21}(x_{j+1})-D_{21}(x_{j-1})]u_{n,j+1}\\
 &  & +\frac{1}{4}[D_{22}(x_{j+1})-D_{22}(x_{j-1})]v_{n,j+1}
\end{eqnarray*}
$$

That is, I am solving a linear system of the form:

$$A\boldsymbol{x}=\boldsymbol{b}$$

where $A$ has $2J\times 2J$ entries, according to the structure:

$$\left(\begin{array}{cccccccccc}
\mbox{1st equation}\rightarrow & | & c_{u}^{j} & c_{v}^{j} & c_{u}^{j+1} & c_{v}^{j+1} & 0 &  & ... & 0\\
\mbox{2nd equation}\rightarrow & | & d_{u}^{j} & d_{v}^{j} & d_{u}^{j+1} & d_{v}^{j+1} & 0 &  & ... & 0\\
\vdots & | & . & . & . & . & . & . & . & .\\
\mbox{1st}\rightarrow & | & 0 & 0 & c_{u}^{j-1} & c_{v}^{j-1} & c_{u}^{j} & c_{v}^{j} & c_{u}^{j+1} & c_{v}^{j+1}\\
\mbox{2nd }\rightarrow & | & 0 & 0 & d_{u}^{j-1} & d_{v}^{j-1} & d_{u}^{j} & d_{v}^{j} & d_{u}^{j+1} & d_{v}^{j+1}\vdots
\end{array}\right)$$

where the $c$'s and $d$'s are the the coefficients corresponding to the scheme above.

Now, I am trying to implement the Neumann boundary condition (the one on $v$), using the ghost point method as discussed, for instance, here for a simpler equation.

Thus, I would like to change my code in order to add some extra terms to the array, but I can't work out what to put in them. Also, in the respective components of the matrix $A$.
Because when I write the condition:

$$\frac{u_1-u_{-1}}{(\Delta x)^2}=0$$
(and same for $v$)

I need to eliminate $u_{-1}$ and $v_{-1}$, using also the third and fourth equations I wrote in this post. If I inverted for $u_{-1}$ and $v_{-1}$ (assuming I know how to), I do not see how I could write down a code which works for any kind of scheme I have, without putting the explicit inversion in the matrix.


问:
Parameter ESTimation (PEST) algorithm? Steps and purposes?
I'm reading up on the Parameter ESTimation algorithm, (Link to PEST homepage 1), and having trouble understanding the broad strokes, as I keep getting lost in null spaces, different regularized ion techniques, etc.

Can someone who is familiar with this algorithm, often used in the calibration of environment models, provide a broad outline of its steps and the purpose of each step?

It is based on the Gauss-Marquadt-Levenberg algorithm and uses many different types of "regularization".

Edit: I understand the part of the algorithm relate to the method of Steepest Descent, but not the Gauss-Newton algorithm part or how it combines these two methods into one.

Much thanks.


问:
Algorithm for merging mesh with cad file
I am writing the pre processing program for a porosity based cfd project we just started. Basically I have a 3 dimensional mesh made of cubes and I need to import a STL file over it, and calculate how much space of each cell is being used by the object. The STL file is a CAD format that uses triangles to represent a surface, but I will assume that it represents a closed surface (a solid).

I made a drawing of the 2D problem:

On the STL object is over the mesh. Each cell of the mesh receives a value ranging from 0 to 1 based on how much of the object is over it. Cells 1,2,3,4,5 have a value of 0 while cell 6, have a value of ~.15 and, cell 7 have a value of 1.

So my question is, where do I start? Is there any "surface casting" algorithm? How do I deal with the triangle-cubes conversion?


问:
Efficient assembly of finite element matrix(coupled equations case)
I noticed this post, where spalloc and sparse are recommended for efficient assembly in Matlab. I personally use sparse assembling for simple cases.

However, when it comes to the case of coupled PDE, say, 3-PDE coupled, then the scalar unknown becomes a $3\times 3$ tensor. In this case, I can't figure out a way to exploit sparse(), and have to use for-loops to assemble.

When the final assembled sparse matrix is as large as $30\text{k}\times 30\text{k}$(on PC), the assembling process becomes really slow (~10min), while the final matrix solving step is still fast(~less than 3 seconds).

Is there any suggestion? Any generic solution is appreciated (not necessarily Matlab)

ContentLicense: CC BY-SA 4.0


问:
What disables Gradual Overflow?
After studying the process of Gradual Underflow, I'm left a little curious as to why machines don't implement Gradual Overflow; where numbers exceeding the overflow level would be stored as un-normalized floats.

Just as the precision of the mantissa is decreased for representation of a larger (in magnitude) exponent (and thus a smaller number in scientific notation) for subnormal numbers, so too could precision be traded for numbers above the overflow.

The only conceivable reason I can place for this is based on number format preservation during gradual underflow; the mantissa digits are 'right-shifted' (preceding zeros are adding after the decimal point) so although the representation is no longer normalized, it still consists of a single digit before the decimal point and the same storage spaces of mantissa and exponent.

The best conception of a Gradual Overflow implementation I can muster would involve either expanding the exponent storage space (at the cost of mantissa storage space and thus precision), or 'right-shifting' the decimal point down the mantissa (again, a loss of mantissa precision). This would seemingly break the mantissa format of a single digit before the decimal point (or otherwise change the mantissa / exponential allocated storage spaces).

I suppose I can understand that in a binary representation (with no delimiters between representation components (ie sign, mantissa, etc); distinction only made via the constant lengths of the components), breaking the storage pattern as Gradual Overflow would require is unsupportable. This assumes that there isn't a better implementation than the one I personally conceived that can preserve the storage pattern.

Is this line of thinking correct, or is there a more pertinent reason as for the lack of Gradual Overflow?

Thanks!

(I can't really tag this appropriately, since there is a rather small pool of existing tags and I don't have the reputation to create any)

Title: What disables Gradual Overflow?


问:
Principal Components Analysis Not Behaving as I Expect it to
I have a bunch of points in $\mathbb{R}^3$ that I would like to translate and rotate so that their center is at the origin and the variance along the $x$ and $y$ axes are maximal (greedy, and in that order). To accomplish this I am trying to use python's principal components analysis algorithm. It is not behaving as I expect it to, most likely due to some misunderstanding about what PCA actually does on my part.

The Problem:
When I center and then rotate the data, the variance along the third component is greater than along the second. This means that, once centered and rotated, there is more variance in the data along the $z$ axis than there is along the $y$. In other words, the rotation is not the correct one.

What I am Doing:
Python's PCA routine returns an object (say myPCA) with several attributes. myPCA.Y is the data array, but centered, scaled, and rotated (in that order). I do not want the data to be scaled. I simply want a translation and a rotation.

import numpy as np
from matplotlib.mlab import PCA

# manufactured data producing the problem
data_raw  = np.array([
                     [80.0, 50.0, 30.0],
                     [50.0, 90.0, 60.0],
                     [70.0, 20.0, 40.0],
                     [60.0, 30.0, 45.0],
                     [45.0, 60.0, 20.0]
                     ])

# obtain the PCA
myPCA = PCA(data_raw)

# center the raw data
centered = np.array([point - myPCA.mu for point in data_raw])
# rotate the centered data
centered_and_rotated = np.array([np.dot(myPCA.Wt, point) for point in centered])
# the variance along axis 0 should now be greater than along 1, so on
variances = np.array([np.var(centered_and_rotated[:,i]) for i in range(3)])
# they are not:
print(variances[1]>variances[2]) #False; I want this to be True

# Now look at the PCA output, Y. This is centered, scaled, and rotated.
# The variances decrease in magnitude, as I want them to:
variances2 = np.array([np.var(myPCA.Y[:,i]) for i in range(3)])
# This looks good, but the coordinates have been scaled.
# Let's try to get from the raw coordinates to the PCA output Y
# mu is the vector of means of the raw data, and sigma is the vector of
# standard deviations of the raw data along each coordinate direction
guess = np.array([np.dot(myPCA.Wt, (xxx-myPCA.mu)/myPCA.sigma) for xxx in data_raw])
print(guess==myPCA.Y) # all true
```

The last two lines in the above show that we may take a point $\mathbf{x}$ from its representation in the raw data input into its representation $\mathbf{x}'$ in terms of the PCA axes via
$$
\mathbf{x}' = \mathrm{R}\cdot\left((\mathbf{x}-\boldsymbol{\mu}) / \boldsymbol{\sigma} \right)
$$

where $\mathrm{R}$ is myPCA.Wt, the weight matrix, $\boldsymbol{\mu}$ is the vector of means of the original data along each coordinate axis, $\boldsymbol{\sigma}$ is the vector of standard deviations of the original data along each coordinate axis, and the division is element-wise. In order to write this in standard mathematical notation, let's replace this division by multiplication:
$$
\mathbf{x}' = \mathrm{R}\cdot\left(\mathrm{D}\cdot(\mathbf{x}-\boldsymbol{\mu}) \right)
$$
where $\mathrm{D}$ is a diagonal matrix whose diagonal entries are $1/\sigma_i$.

This notation makes clear the problem: to undo the scaling, I need to act on the RHS above with $\mathrm{R}\mathrm{D}^{-1}\mathrm{R}^{-1}$. This will return me to the problem situation, in which the variance is greater along the $z$ axis than the $y$.

Is there a way to use PCA to get what I want, or do I need to use another method?


问:
Printing a plot in Octave-3.8.1 breaks LaTeX labels
Since upgrading to Octave 3.8 I have been unable to export plots with LaTeX labels. When the figure is displayed the labels render fine but when I try to export the figure using the print command the labels come out as a raw string.

For example:

```<code>x = 0:0.01:(2*pi);
y = sin(x.^2 + cos(x)).^3;
plot(x,y)
xlabel('x')
ylabel('sin^3(x^2 + cos(x))')
print('test1.png')
</code>```

produces,

and saves,

Does anyone know how to fix this or hack around it?


问:
Compute hypergeometric function ratio: $\frac{_{2}F_{1}(a+1,b;c;x)}{_{2}F_{1}(a,b;c;x)}$, for large positive $a,b,c$?
I need a numerically stable way to compute the following ratio:

$$\frac{_{2}F_{1}(a+1,b;c;x)}{_{2}F_{1}(a,b;c;x)}$$

All the parameters are real positive numbers, with $0 < a,b,c$ and $0 < x < 1$.

Right now I am using GSL's implementation of the hypergeometric function, but I keep getting under/over-flows.

Is there a simplification that I can use? Or an aproximation scheme that remains accurate for a wide range of the parameters (say $0 < a,b,c < 1000$)?


问:
Histogram separation energy equation based on Bhattacharyya
I am working in level set method, specially Lankton method paper. I try to implement Histogram Separation (HS) Energy problem (Part III.C). It based on Bhattacharyya to control the evolution of contour.

To understand it, the first consider the global method in which we are given an input image and a contour. The contour divides the image into an inside and outside region. The Bhattacharyya distance is calculated by

$$B_{global}=\sqrt{P_{u}*P_{v}}$$

where ${P_{u}}$ and ${P_{v}}$ are pdf of inside and outside regions.

Now we return the Lankton paper which is a local level set.  He divides the image into small region by Ball function B(x,y). Then, the contour will separate these regions into inside and outside region. Each small regions have ${P_{u}}$ and ${P_{v}}$. And we can calculate Bhattacharyya distance.

$$F_{HS}=B_{local}=\sqrt{P_{u,x}*P_{v,x}}$$

where $x$ is center point of $B(x,y)$ in image.

I understand that that term. But I cannot understand second term in the formula (integration term, in $\color{red}{\textrm{red}}$ below). Can I represent it by convolution operator?  Here is my source code. It works but the result is not correct.

$$
\begin{aligned}
\frac{d\phi}{dt}(x) =\ &\delta\phi(x) \int_{\Omega_y}\frac{\mathcal{B}(x,y)\delta\phi(y)}{2} \times \left[F_{HS}\left(\frac{1}{A_v} - \frac{1}{A_u}\right)\right. \\
&+\left.\color{red}{\int_z K(z - I(y)) \times \left(\frac{1}{A_u}\sqrt\frac{P_{v,x}(z)}{P_{u,x}(z)} - \frac{1}{A_v}\sqrt\frac{P_{u,x}(z)}{P_{v,x}(z)}\right) dz} \right] dy \\
&+ \lambda\delta\phi(x)\ \mathrm{div}\left(\frac{\nabla\phi(x)}{\left|\nabla\phi(x)\right|}\right)
\end{aligned}
$$

where $K$ is a Gaussian kernel and $A_v$ and $A_u$ are areas inside and outside of these regions.

```<code>for its = 1:max_its   % Note: no automatic convergence test

%-- get the curve's narrow band
idx = find(phi <= 1.2 & phi >= -1.2)';
[y x] = ind2sub(size(phi),idx);

%-- get windows for localized statistics
xneg = x-rad; xpos = x+rad;      %get subscripts for local regions
yneg = y-rad; ypos = y+rad;
xneg(xneg<1)=1; yneg(yneg<1)=1;  %check bounds
xpos(xpos>dimx)=dimx; ypos(ypos>dimy)=dimy;

%-- re-initialize u,v,Ain,Aout
Ain=zeros(size(idx)); Aout=zeros(size(idx));
B=zeros(size(idx));integral=zeros(size(idx));
%-- compute local stats
for i = 1:numel(idx)  % for every point in the narrow band
  img = I(yneg(i):ypos(i),xneg(i):xpos(i)); %sub image
  P = phi(yneg(i):ypos(i),xneg(i):xpos(i)); %sub phi

  upts = find(P<=0);            %local interior
  Ain(i) = length(upts)+eps;

  vpts = find(P>0);             %local exterior
  Aout(i) = length(vpts)+eps;

  %% Bha distance
  p = imhist(img(upts))/ Ain(i) + eps; % leave histograms unsmoothed
  q = imhist(img(vpts)) / Aout(i) + eps;
  B(i) = sum(sqrt(p.* q));
  term2= sqrt(p./q)/Aout(i) - sqrt(q./p)/Ain(i); %Problem in here===I don't know how to code the integral term
  integral(i) =sum(term2(:));%integeration can do by this code?
end

F =-B./2.*(1./Ain - 1./Aout) - integral./2;
</code>```

OwnerDisplayName: johnmark


问:
Accurate computation of $\frac{\mathrm{B}_{x,y}(\alpha + 1,\beta)}{\mathrm{B}_{x,y}(\alpha,\beta)}$ for large paramers?
I need to calculate these ratios:

$$\frac{\mathrm{B}_{x,y}(\alpha + 1,\beta)}{\mathrm{B}_{x,y}(\alpha,\beta)} \tag{1}$$

where $\alpha,\beta > 0$ and $0\le x\le y \le 1$. Here $\mathrm{B}_{x,y}(\alpha,\beta)$ is the incomplete Beta function, defined as:

$$\mathrm{B}_{x,y}(\alpha,\beta) = \int_x^y t^{\alpha-1}(1-t)^{\beta-1}\mathrm{d}t$$

For large values of the parameters $\alpha,\beta$ , the numerator and the denominator can be very $\mathrm{B}_{x,y}(\alpha,\beta)$ can get very large or very small, resulting in over/underflow.

Is there a way that I can compute the ratio (1) avoiding these over/underflows?


问:
linear solution of curve fitting on multiple linear functions differing by a multiplier
I am facing the following problem. I know nonlinear least squares can provide a solution but I am wondering if a linear way to solve this data fitting problem may exists.

This is my input dataset:
I've got three different dataset composed of scattered points

I know a linear equation in the form
$$
y(x) = bx^2 + cx
$$

can be used to explain any of my dataset. I know how to fit this function to a given dataset (e.g the blue one) linearly, so  could the above function separately for each dataset, but I am looking for something different.

In the specific case I have an additional constrain: I also know that the three functions describing the three dataset share the $b$ and $c$ parameters, while they differ by a multiplier, such that:

$$
\begin{cases}
y(x) = (bx^2 + cx) \; \text{explains the blue data}\\
y(x) = a'(bx^2 + cx) \; \text{for red data}\\
y(x) = a''(bx^2 + cx) \; \text{for green data}
\end{cases}
$$

I am looking for a (if it exists) linear way to solve for $a', a'', b$ and $c$ minimizing the overall sum of squares for the three functions.
Any idea on how this problem can be approached?
Also approximate solutions are welcome...


问:
Does it matter if I use principal component analysis on the transpose instead of the original matrix?
My data set is a 60x10 matrix. I performed principal component analysis of this matrix with matlab using the princomp(AdjustedData) after I adjusting my original data set by subtracting the mean of each column. Because I was hoping to obtain a data set described by the top five principal components, I selected the top five vectors in the COEFF output as my feature vectors. I firstly multiply the row feature vectors (vector on each row, 5x10) by my row adjusted data matrix (transposed adjusted data matrix with different dimensions on each row, 10x60) to obtain the so called final data.

```<code>           FinalData (5x60)=Row feature vector(5x10) * row AdjustedData(10x60)
</code>```

Then, I multiply the transposed row feature vector(10x5) by the FinalData(5x60) to get a matrix as the same dimensions as my transposed original data set.

TransOriginalAdjustedData(10x60)=TransposedRowFeatureVector(10x5)*FinalData(5x60)
Finally, I transpose the output matrix and add the mean of each column back to retrieve the final original matrix (PCA data, a data set I thought described by only 5 principal components) that I want to compare with my original data set.

```<code>    PCA data=OriginalAdjustedData(60x10)+ original column mean
</code>```

I sum up the square difference between the data points of my original data set and the PCA data set, and take a square root of the sum. The result was about 0.022.

However, if I transpose my original data matrix, 10x60 instead of 60x10, and performed the same procedures, the result was only 0.015. I know the dimension of the matrix is changed when performing PCA if I transpose the original matrix, but I cannot figure out what is causing the large difference between the two. The variables in my data set were not defined so I thought either way should work about the same.

My procedures are based on this PCA tutorial: http://www.cs.otago.ac.nz/cosc453/student_tutorials/principal_components.pdf

Could anyone help me with this? Thank you very much.


问:
Treatment of Neumann (Traction) boundary conditions using projection methods
I am looking to solve the incompressible Navier-Stokes equations in 3D, using an inflow boundary condition specifying a velocity:

$\mathbf{u} = \mathbf{g}_0 \,\, \forall \,\, \mathbf{x} \in \Gamma_u$

and I want to use a Neumann type boundary condition where I specify pressure at the outlet

$-p\mathbf{I} + \tau = p_0\mathbf{I} \,\, \forall \,\, \mathbf{x} \in \Gamma_h$

I am looking to use a projection method to solve for the combined velocity and pressure fields, however I am not certain as to how to incorporate the Neumann term above into such a scheme. Most of the literature talks about deriving a velocity based boundary condition for the discrete poisson equation - but that is all for problems where the Neumann boundary $\Gamma_h = \phi$.

I was wondering whether anyone here had any insights on this. I can give more details of the mathematics that I am looking at if needed. Any references, or web-pages, or some brief pointers on this will be very helpful.


问:
FEniCs: help in implementing the boundary condition for 1D problem
I have just started learning FEniCS and have used: http://www.scientificpython.net/pyblog/fenics-linear-two-point-bvp to write a script for solving:

```<code>           u'' + u = 1

           u(0) = 1, u'(1) = 0
</code>```

with exact solution,

```<code>           u(x) = exp(-x)/ exp(-1) + x + ( -1+exp(-1) )/ exp(-1)
</code>```

Clearly the weak formulation of the above problem is:

```<code>           -(u', v') + (u,v) = (g,v) ; with g = 1
</code>```

Here is the edited code:

```<code>    from dolfin import *

    # definig mesh
    mesh = IntervalMesh(20, 0, 1)

    # definig Function space on this mesh using Lagrange polynoimals of degree 2.
    V = FunctionSpace(mesh, "CG", 2)

    # definign boundary values
    #u0 = Constant(0)
    u0 = Expression("x[0]")

    # this functions checks whether the input x is on the boundary or not.
    def DirichletBoundary(x, on_boundary):
         tol = 1e-14
         return on_boundary and abs(x < tol)

     # Enforcing u = u0 at x = 0
     bc = DirichletBC(V, u0, DirichletBoundary)

     # Setting up the variational problem
     u = TrialFunction(V)
     v = TestFunction(V)
     f = Constant(1)
     g = Constant(1)
     a = -inner(grad(u), grad(v))*dx + inner(u,v)*dx
     L = f*v*dx

     # solving the variational problem.
     u = Function(V)
     solve( a == L, u, bc)

     # plotting solution
     plot(u, interactive = True)
</code>```

Clearly the solution plot does not incorporate the boundary condition u(0) = 1.

Would any body please help?