问:
Prediction model on GPU
I am doing a small project at school. I have done my code implementations in CUDA and did some performance measurements with real values, i.e. running the program with different number of threads, changing the size of the problem, changing both, etc.

The results of (time, speedup) look really nice. But then I was asked to do some theoretical analysis of the algorithms, specifically do a prediction of the time it will take my implementation to run as the number of threads increases. The problem was that I received some documents where they use "Timing diagrams" to try to model this, but this approach was for problems solved using MPI and I'm being asked to use exactly the same stuff to create a prediction model for CUDA. Is it possible to do this?

How can I do a prediction of how a given CUDA program (kernel) will run as I increase the number of threads used for computation? Say that my GPU has 240 cuda cores, now for sake of simplicity I am taking that if I launch a kernel with one thread that kernel will use only 1 core, if I launch 4 threads it'll use 4 cores up to a limit of 240 in that case when I launch 500 threads there won't be no more than 240 cores in use. It is highly likely that this is somehow wrong but It's something to start with.

Could you please give me some ideas to sketch this prediction model? I am really cracking my head because of this, I cannot came to a simple prediction model.

UPDATE

The very first problem was how to simulate a sequential running on the GPU in CUDA. As mentioned before what I did was launch a kernel with only one block and one thread, the result was what I expected as it turned out that the kernel ran slower. Then increasing the number of threads launched reduced the running time so I was getting benefit from the parallel architecture of the GPU. Now, the key question is launching K threads in CUDA, can I be sure that the number of cores used in the GPU equals the number of K threads? Somehow this is true because of what I said before.

Now moving forward I have to make a prediction model, based on timing diagrams. So first the sketch of the timing diagram is

What my diagram shows (or at least tries to show) is that the critical path for the application to run depends on

Sending data from CPU to GPU

CUDA initialization and kernel launch configuration

The time it takes the last warp to execute

Sending data from GPU to CPU

The time of (1) and (4) can be straight forward computed: I know the GPU is connected to the CPU via PCI-Express 2nd gen and knowing the size of the problem I can know the actual time of transfers. Good :)

The time of (2) is hardware dependent, but for what I've been reading it takes around 60 to 65 ms.

Finally (3) the hardest part to predict is the actual running time of the computation (kernel runtime for a given number of threads). The way I modeled this was Tprocessing = Tcomputation + Tmemory

Tcomputation : number of cycles one thread does to compute 1 pixel of the set. The key question here is How to count the number of cycles? I naively counted the number of sums and multiplications and then multiplied it by the cycles each of them takes according to the architecture of the GPU. I am completely uncertain about this.

Tmemory : number of cycles one thread takes to go to global memory to read/write. From various references (including Udacity course on GPU) is no average 500 cycles.

Then we know the frequency at which the GPU clock works and summing all I would have an upper bound of the time it will take my program to run.

So basically everything is based on the idea that if I run one thread only one CUDA core is used, if I run ten threads then ten CUDA cores are used. Then I can predict what will is likely to happen if we have a GPU with 2048 CUDA cores as the GTX Titan.

I test my model with the Mandelbrot Set and well my timing was much bigger than the actual one, so something is missing.

Title: Prediction model on GPU


问:
Benchmarks or generic configurations for optimal flow control
I am about to test my algorithms for solving optimal control problems of type:

Find an input $u$, such that for a time interval $(0,T]$ the cost functional
$$J(v,u) = \mathcal M(v(T)) + \int_0^T\mathcal K(v(t),u(t)) dt, $$
becomes minimal, where $v$ is the velocity of flow modelled by the incompressible Navier-Stokes equations. The functionals $\mathcal M$ and $\mathcal K$ measure the state trajectory, e.g. in terms of how well a given target flow is approximated.

For 2D, typical test problems I have seen comprise the backward facing step, driven cavity (as in the pics below) or the cylinder wake. However, I haven't found discussions on

What is a good configuration, or even a benchmark, for testing numerical solution approaches to optimal flow control?

I want to do my first tests on a very basic level. That's why I think the setup should be

reasonably sized (~$10^5$ DOFs),

laminar,

and for distributed control.

Here, distributed means that the control appears in the right hand side of the Navier-Stokes equations and not in the boundary conditions.

ContentLicense: CC BY-SA 3.0


问:
How to apply periodic boundary conditions to a Raviart-Thomas finite element space in Fenics?
I'm trying to use periodic boundary conditions within a Raviart-Thomas finite element space in Fenics (dolfin 1.2.0) in a Ubuntu 12.04 (amd64) machine (with python 2.7). If other FE space is used, everything works fine, but if a Raviart-Thomas is considered, python crashes (due to libdolfin errors). In what follows, it is attached a test code to reproduce the error:

```<code>from dolfin import *

# Periodic boundary condition class
class PeriodicWall(SubDomain):
    def inside(self, x, on_boundary):
        return bool(on_boundary and near(x[0], 0.5))
    def map(self, x, y):
        y[0] = x[0]
        y[1] = x[1] + 1.0

# Mesh and periodic boundary conditions
mesh = RectangleMesh(-0.5,-0.5,0.5,0.5,5,5,"left")
periodicwall = PeriodicWall()

# Raviart-Thomas space with periodic boundaries
V = FunctionSpace(mesh, "RT", 1, constrained_domain=periodicwall)
</code>```

Any kind of help would be helpful: are there any errors in the test code? or is it a known bug in Dolfin?


问:
What is the best OpenFOAM RAS turbulence model for a motorbike/vespa problem?
I am learning OpenFOAM as a hobby and using my Vespa racing as the topic to apply it to. The objective is to produce modifications that improve the top speed (as well as getting some values such as lift). I am basing my simulation on the motorbike tutorial.

The motorbike tutorial uses the kOmegaSST turbulence model but someone had suggested to me in the past that SpalartAllmaras might be a better choice. I have access to a super computer so saving CPU cycles is not much of an issue for me. Below is a picture of my model (which I am cleaning up). Is kOmegaSST a good enough choice or is there another worth focusing on?

(http://www.openfoam.com/features/RAS.php)!


问:
Sine series using exponential based FFT
I have such a problem - I would need to expand a discrete function in a sine fourier series but I would like to use exponential based library for FFT (I will use CUDA to compute it).
What have I to do with the signal I'm sending to the procedure so this is treated as sine fourier? I think I should add zeros but I am not sure where they have to be.
Many thanks</p>

Title: Sine series using exponential based FFT


问:
Stochastic Galerkin projection approach for using generalized polynomial chaos expansion (GPCE) in solving PDE
I want to know if there is any way to define the test and trial function in the way that I want instead of using the default functions. So if I want define the polynomial and basis and coefficient, how can I do it in Fenics? Here is the concept in general:

I am trying to solve stochastic partial differential equation which is like
$$\Pi(x,t,w;u)=f(x,t;w)$$   $$x \in \Re^d  \qquad t \in [0,T] \qquad w \in Ω $$
where is the solution is $U:=U(x,t;w)$ and $f(x,t;w)$ is the source term.
We can discretize the stochastic space using a finite number of independent random variables.
$$\xi(w)=\{\xi(w_1),\ldots,\xi(w_n)\}$$
This allows us to rewrite the governing equation as a dimensional
differential equation in the strong form:
$$\Pi(x,t,\xi;u)=f(x,t;w)\qquad (\xi,x)\in \Gamma \times D$$
Thus the solution which is regarded as a stochastic process, can
be expanded by the generalized polynomial as:
$$u(x,t;w)=\sum_{i=0}^P u_i(x,t) \phi_i (\xi)$$
where $P$ is the number of GPCE terms.

Then we substitute $U$ into the weak form.
Choosing the test function the same as the generalized polynomial
in order to let the residual to be orthogonal to the subspace, thus
minimize the error, we obtain:
$$  \left\langle \Pi\left(x,t,\xi;\sum_{i=0}^P u_i(x,t) \phi_i (\xi)\right),\phi_k (\xi)\right\rangle = \left\langle  f(x,t;\xi),\phi_k (\xi)  \right\rangle   $$   $K=0,1,.. P$

Results in a set of coupled deterministic equations. One can
employ e.g. classical finite element methods to seek the
deterministic coefficients.

In order to make the problem more clear there is well defined problem related to 2D Poisson Equation: here on page 7.

Since we only have one random variable, the number of random dimensions is 1. Polynomial chaos expansion value for my problem with polynomial of order of 3 will be a PCE=3*3 by matrix.

Problem

I need to modify the V = VectorFunctionSpace(mesh, "Lagrange", 1)so in this problem in order that V will be the Kronecker tensor product of V and PCE which will lead to a 10*3 ,10*3, 10*3 and based on that matrix I can modify my equation to calculate the statistics of the solution. Is there any way that I can calculate the Kronecker tensor product of V and another matrix and replace it with my V?

More specifically  I am trying to solving this equation with uncertain location load with this method.

Here is the my deterministic code,Fenics: Result of Steady state dynamic linear elastic doesn't match with actual values  how I can make this deterministic code to stochastic with modifying test and trial function in Fenics.


问:
Computation of potential flow using dirichlet conditions
I know I already posted this question and I thank you for your answers, but unfortunately I didn't find what I was looking for among them. Anyway now I'll rewrite the question more clearly, because before perhaps I've not been specific enough.

I've seen that there is a way to use the finite differences method, on a cartesian orthogonal grid, to perform calculations on potential flow about an obstacle without using the Neumann conditions, but only Dirichlet conditions. It is done using the stream function, i tested it and it seem to work.
This method works this way:
You take a rectangular box, then you put two arbitrary values of the stream function (but the same in all the side) in the upper and lower side; then you linearly vary it along the lateral sides (so there is continuity with the potential on the upper and lower sides). In this way you define the undisturbed flow far from the obstacle. Then you place the symmetric obstacle right in the middle of the box, so you know for sure there is an horizontal streamline that divides the box in two and follows the obstacle boundary, so its potential is exactly the average value between the upper and lower side. So you have defined the Dirichlet conditions in all the box boundaries and in the obstacle boundary, and the problem is well defined and can be easily solved by using the finite difference scheme: you solve the Laplace equation for the stream function.

But this method I used has a big drawback: it can be used only with obstacles that are symmetric along an horizontal axis (the direction of the undisturbed stream). If the obstacle is asymmetric, the streamline that follows the obstacle is no longer horizontal outside the obstacle, and I don't know where it hits the box wall; then I don't know what is the stream function value that I have to put in the obstacle boundary.

Do you now any workaround or method to use only Dirichlet conditions, but on an arbitrary shaped obstacle (a wing, for example) ?

Before, someone advised me to use panel methods, coordinate transformations, etc. It would be good, but I'm looking for something simpler, if possible. It's just for demonstration of the finite differences method, it doesn't have to be particularly efficient of fast; I'd prefer just the simplest solution you can think. Something I can possibly add to my program editing it the less possible. If instead you know a simple way to implement the Neumann conditions, it would also be good.

Thank you in advance, your help would be greatly appreciated.


问:
Strong stability preserving RK scheme
For the ODE

$$
\dot{x} = f(x)
$$

we have the 2-stage, second order SSP RK scheme (Shu, Osher, Gottlieb)

$$
x^{(0)} = x^n
$$

$$
x^{(1)} = x^{(0)} + \Delta t f(x^{(0)})
$$

$$
x^{(2)} = \frac{1}{2} x^n + \frac{1}{2}[ x^{(1)} + \Delta t f(x^{(1)})]
$$

$$
x^{n+1} = x^{(2)}
$$

What is the scheme in case of the following ODE

$$
\dot{x} = f(x,t)
$$

In particular I am looking for papers which deal with this case.

Title: Strong stability preserving RK scheme


问:
Data to test optimization algorithms
I plan to implement some optimization algorithms in Python, say to follow the sequence here or here.

What would be a systematic approach to testing the algorithms I write, using either well-known test datasets or creating my own?

Title: Data to test optimization algorithms


问:
Non-convergance when calculating temperature/heat flows through a section of rock
I am attempting to calculate temperature of section of rock in the earth as a function of vertical position in the rock and time.  Along with it I am calculating the heat flow through the rock as a function of vertical position in the rock and time.  Note that this is a one-dimensional examination of the problem- imagining that the rock is uniform in every direction but up-down.

Please be patient with me: I currently do not have the vocabulary nor the background to talk at the level of familiarity with FEM and differential equations I see in other questions on this site.  I am a software developer with some physics background, which has allowed me to get this far.

Description of methodology: I am currently splitting up the section of rock into horizontal slices (depth steps), and then calculating the temperature and heat flow of each depth step at specific time-steps.  I am assuming that heat only flows upward from below the section of rock.

As I understand it, the temperature difference between the top and bottom of a depth step would be expressed as:
$$
\Delta T = \frac{Q\Delta Z}{\lambda}
$$
Where $Q$ is the heat flow at the middle of the depth step in $mW/m^2$, $\Delta Z$ is the thickness of a depth step at this time, and $\lambda$ is the thermal conductivity of the rock (constant).

Additionally, the change in the heat flow between the bottom and top of a depth step could be expressed as:
$$
\Delta Q = - C_p \rho \Delta Z \frac{\Delta T}{\Delta t}
$$

where $C_p$ is the heat capacity of the rock (constant), $\rho$ is the density (constant), and $\frac{\Delta T}{\Delta t}$ would be the change in temperature of the depth step over the course of this time-step.

As a starting point for the calculations, I know the temperature at the top of the rock layer as a whole (but not the temperature at the bottom of the rock layer). I also know the heat flow value just below the rock layer (but not the heat flow value at the top).

The key piece of complexity to this is that the rock layer is increasing in thickness through time (aka deposition), so that I am adding to the $\Delta Z$ of the top-most depth step and then adding more depth steps to the top of the section of rock as time progresses.  As a consequence, individual depth steps get more and more buried, and thus hotter through time, which thus causes the heat flow to decrease as you move upwards in the rock layer (as per equation 2 above).

My method of calculating the temperatures and heat flows for each time-step is as follows:

Increase the thickness of the rock layer, adding depth steps to the top as necessary.

Calculate the temperatures at the top and bottom of each depth step going down the layer (since I know the surface temperature, and thus the starting point), using the heat flow values from the previous time-step.

Calculate the heat flow values at the bottom and top of each depth step going up the layer (since I only know the heat flow at the very bottom of the rock layer), using the temperatures of each depth step calculated in step #2, along with the temperature in the previous time-step.

Repeat steps 2-3 until the solution converges. Note that in step #2 I can now use the heat flows calculated in step #3.

The Problem:  This solution does not always converge.  As the rock layer gets thicker and thicker through time, each time-step takes longer and longer to converge, until it reaches a tipping point and starts to diverge.  What can I do to force this to converge?

I am suspicious that the reason for divergence is because I'm forced to calculate temperature down the depth steps and heat flows up the depth steps.  This means that the bottom-most depth step's temperature is sensitive to ALL the changes in temperature above it as per the first equation above.  But this then causes the heat flow at the bottom to change from the previous iteration by a fair amount, which then propagates all the way to the top of the rock layer via the second equation above, which then affects the temperature at the top strongly, etc.

Warning: The above problem is actually a major simplification of a larger problem I'm trying to solve- the constants in the equation above are actually temperature dependent and depth dependent in a complicated, non-continuous way.  However, the problem I am facing still exists in the simplification as described above.


问:
Block Backward Differentiation Formula (BBDF), on order 4 formula
I am trying to implement a program the numerical method to solve ODE called Block BDF as explained in this article: https://waset.org/journals/waset/v38/v38-49.pdf

As it is variable step-size, I need to compute the LTE (local truncation), that is $LTE=y_{n+2}^{(5)}-y_{n+2}^{(4)}$. There is a formula in the article for $y_{n+2}^{(5)}$ but I cannot find one for $y_{n+2}^{(4)}$; where $y_{(n+2)}^{(k)}$ means the $k-$order approximation to $y_{(n+2)}$.

Can someone please help me finding the formula for fourth order approximation? (As I'm only trying to implement, information only about the formula can be helpful, but explanation is encouraged).

Thank you very much!

PS. The more general method is BDF: http://en.wikipedia.org/wiki/Backward_differentiation_formula but this method works in block of 2 (works for $y(n+1)$ and $y(n+2)$ at the same step).


问:
Which turbulence model: k-epsilon or k-omega?
I'm modelling a tidal harbour in Ansys Fluent and want to include a turbulence model in my simulation. I want to predict the water movements around a pontoon which is situated in the harbour. Which turbulence model should I use, a $k$-$\varepsilon$ model or a $k$-$\omega$ model?

The harbour is modelled as a rectangular ($120\times 80\times 6$ m LxWxH) box and the water flows in through an entrance channel ($40 \times 10 \times 6$ m). To simulation flow around a pontoon, a rectangular box ($40\times 3\times 5$ m) is placed at the long side of the harbour and water can flow around an underneath the pontoon. Velocities are 0.3 m/s and lower.


问:
Operator Splitting methods for DAEs
After doing some research, I've found that most of the literature on operator splitting methods (e.g. Strang Splitting, Fractional Step, etc.) are specifically designed for a standard problem type of the form:

$$u'=(A+B)u.$$

In my application, I have a slightly different form:

$$Mu'=(A+B)u$$

where $M$ is not the identity operator.  Can the same methodologies be applied when an operator is applied to the $u'$ term?  Are there any other numerical methods to handle this case specific case there $M$ is not the identity operator?  In particular, I'm looking at a case where $M$ is singular.

Update

To put my particular problem into greater context, my operator is

$$M=\left[\begin{array}{cc} 0 & 0 \\ I & I\end{array}\right]$$

and

$$u=\left[\begin{array}{c} u_1 \\ u_2 \end{array}\right]$$

Title: Operator Splitting methods for DAEs


问:
Hatree-Fock, reasons for convergence/ non-convergence
I'm new here so please forgive me if I lack proper stack exchange etiquette.
So, I was wondering if anyone here could provide insight on a problem that I am running into with with a Hartree-Fock calculation. I wrote a code from scratch that attempts to calculate the HF energy and basis for a 1D potential that is supposed to model the nucleon-nucleon two body interaction:

$$
V(x)= 60 e^{-25x^2} - 15 e^{-1.56x^2}
$$

If I run the code for so called "spinless" fermions, it converges and gives results similar to those given in literature on this potential. However, if I add any spin degrees of freedom to the problem, the SCF procedure does not converge. Instead it arrives at a group of two or more sets of eigenvectors which when used to calculate the new fock-matrix will eventually reproduce one another.     I want to believe that this is because I am trapped in a local minimum in the HF functional.   So my main question is: Have any of you ever encountered such a problem with a basic HF code? Considering that in some cases I get good results and in other cases no results, I feel that the problem lies in the physics of the potential, rather than in my code.   Is this a reasonable statement in your eyes?


问:
How do I implement thin plate splines with barriers?
I want to implement thin spline interpolation of scattered elevation data $ \{z_i(x_i,y_i)\}_{i=1..n} $ in C++.

This seems fairly simple using Radial Basis Functions:

$$ z(x,y) = p(x,y) + \sum_i l_i\phi(r)$$
where p(x,y) is a polynomial and $\ \phi= r^2\ln(r) $ .

Requiring the interpolated surface to go exactly trough the n datapoints + ortogonality of polynomials leads to a matrix equation. This is then solved to find the coefficients $l_i$.

However I want the implementation to take into account "barriers":

As you can see in this figure: points on one side of the "barrier" should not be influenced by points on the other side of the "barrier".

As an example consider this simple form of interpolation:
the elevation in a grid node is found as an average over the elevation of points within a circle of radius R. Now consider a grid node very close to the "barrier". This would be influenced by points on both side of the "barrier" and as a result the interpolated elevation would be somewhere in between. The resulting surface would not show a discontinuity ("cliff") instead it would show a "hill". In this case things could be remedied by only using points from the same side of the "barrier" as the grid node when doing the averaging.

NB! The "barriers" do not necessarily split the data into two sets like in the image.

How do I do this?


问:
Simulink PDE (MATLAB)
I need to build a system in Simulink that solves a PDE, but I can't find any literature or books where it is described how to do it (especially any stuff according to modeling PDE in Simulink).

Please give some advice where I can find such literature.

Thanks beforehand!

Title: Simulink PDE (MATLAB)


问:
Error in treating non-uniformly spaced data as uniform
In discussing smoothing filters,
Numerical Recipes p. 772 says

... irregularly sampled data, where the values $f_i$ are not uniformly spaced ...

  one can simply pretend that the data points are equally spaced ...
  a rough criterion:

  If the change in $f$ across the full width of the $N$ point window is
  less than $\sqrt{N/2}$ times the measurement noise on a single point,
  then the cheap method can be used.

Where does this come from ?
It looks odd:
for $f(x) = \text{sin} ( 2 \pi f x )$ on [$-N \dots N]$,

change in $f = 2 < \sqrt{N/2}\ \sigma$ would require bigger $N$ for smaller $\sigma$ . Or have I misunderstood ?

Also, does this work for smoothing data on non-uniformly spaced 2d squares, 3d cubes ... ?

Added: https://gist.github.com/denis-bz/5957279 is a short Python program
that tries to elucidate this, for quantization noise
(sample at integers +- random-uniform( -1/2, 1/2 )).


问:
discrete versions of Lp norm
The discrete analogue of the $L_p$ norm for the mesh function $V$ is $$\|V\|_{l^p(\bar{\Omega}^N)}=\left(\sum_{i=0}^NV_i^p\bar{h}_i\right)^{1/p}$$
where $\bar{\Omega}^N$ is an arbitrary mesh, $\bar{h}_i=(h_{i+1}+h_i)/2$ and $h_i=x_i-x_{i-1}$. The error $e$ of a certain approximation is such that $e_i=-e^{-x_i/\epsilon}$ for $1\leq i\leq N-1$ and 0 for $i=0$ or $i=N$

For any $p$, $1\leq p\leq \infty$ the $l^p(\bar{\Omega}^N)$ norm of the error on a general non-uniform mesh is
$$\|e\|_{l^p(\bar{\Omega}^N)}=\left(\sum_{i=1}^{N-1}e^{-px_i/\epsilon}\bar{h}_i\right)^{1/p}$$
My first question is why is this expression $O(N^{-1/P})$.

On a uniform mesh, why is $$\left\vert {e_i}^p\right\vert \leq N\sum_{j=1}^{N-1}\left\vert e_j\right\vert^ph$$ How does this imply $$\|e\|_{\bar{\Omega}^N}\leq N^{1/P}\|e\|_{l^p(\bar{\Omega}^N)}$$

Title: discrete versions of Lp norm


问:
Estimating eigenvalues from time-dependent non-linear operator
I have a very sparse non-linear system $N(u) = 0$ that can be solved as a time-dependent ODE, $\frac{du}{dt} = N(u)$, and explicitly integrated until $\frac{du}{dt} = N(u) = 0$, e.g. by forward euler, $u^{n+1} = u^n+\Delta t N(u^n)$.

In the process of this numerical integration, is there any way to recover the eigenvalues of the system by examining something about the growth/response in $u$ and it's dependence on the timestep? Gershgorin's theorem provides a bound on the spectrum, but in practice this bound's range is too large to be useful. I don't need exact eigenvalues, but relatively tight approximations would be great. I suspect there's something about this in the literature of linear time-invariant (LTI) theory, but I don't know where to start looking.


问:
Generating pseudo-random orthonormal bases for random projection
I am performing series of random projections i.e. projecting the input matrix onto randomly generated orthonormal bases (of much lower dimensionality). The projection is just a matrix multiplication like $A \cdot B = P$, where $B$ is the basis composed of several mutually orthogonal column vectors of unit length.

I have come to conclusion that there is no need to generate new basis before each projection because simple permutation of its rows should be sufficient to keep the randomness of the whole process. If this is true, doing so would be very profitable in terms of computational complexity, since row-permutation also keeps orthonormality of the basis.

The problem is that I am a programmer and no matematician, so I'm not sure if indeed projections over such bases are random and thus independent of each other, or they are somehow similar and thus limited. Is my idea right or wrong, and why?

Edited:

The matrix $A$ consists of $n$ row vectors of binary values ($0$ or $1$), where the vector dimension is close to $n$ which is between $10^5$ and $10^7$.

The process of repeating random projections aims to find well-clusterable projection, or at least avoid worst cases. Currently the clusterability is defined as variance.

In order to minimize complexity of row permutation I'm just rotating row indices, so the original basis is not altered at all.

The final goal is to reduce dimensionality of the input vectors and then cluster them in this low-dimensional space with some distance or density based algorithm.

$$
B = \begin{pmatrix}
1 & 0 \\
0 & 1 \\
0 & 0 \\
\end{pmatrix}
$$

The above counterexample clearly shows potential risk of generating a basis which is insensitive to certain permutation - exchange rows 1 and 2, and the result is the same as changing the order of basis vectors, which does not change the projection at all.

Question 1: What is the probability of generating a random basis with such an inherent flaw.

I see that row permutation somehow reflects the basis, but I don't know which bases do not reflect well. I do not expect exact numbers, since it would require precise definition of words "have", "flaw" and "well". Any help in expressing the problem more profesionally is greatly welcome.

Question 2: Is it the only sort of flaw that can reduce randomness of the entire process?


问:
Curvature of level-sets of CG1 function
Let us have simplicial mesh and continuous function $u$ which is piece-wise linear and non-constant on every cell. Then normal vector to level-sets of $u$ is given $$\mathbf{n}=\frac{\nabla u}{|\nabla u|}$$ and is piece-wise constant function. What is then possible definition of curvature $\kappa$ which would be for $u\in{\cal C}^2$ $$\kappa = \mathrm{div}\,\mathbf{n}.$$ This definition is not suitable as it yields $\kappa$ being zero in every cell and a multiple of Dirac $\delta$ on every facet.

In conclusion: what is a suitable finite-element space for $\kappa$ and weak formulation for $\kappa = \mathrm{div}\,\mathbf{n}$ with given piece-wise constant $\mathbf{n}$?

EDIT: I got another idea: first project $\mathbf{n}$ on some $H(\mathrm{div})$-conforming space; then it is no problem calculating $\mathrm{div}\,\mathbf{n}$ directly. I did some numerical experiments with $u$ being projection of
$$ u^\mathrm{exact} = (1+x)(1+y) $$
to CG1 (continuous, piece-wise linear) on domain $\Omega=(0,1)\times(0,1)$. Exact curvature for $u^\mathrm{exact}$ is
$$ \kappa^\mathrm{exact} = \frac{-2(1+x)(1+y)}{[(1+x)^2+(1+y)^2]^{3/2}} .$$
Taking Raviart-Thomas degree 1 space (piece-wise linears, continuous on facet midpoints) or CG1 space as space for $\mathbf{n}$ method seems to converge in norms but with annoying wiggles on $\partial\Omega$. Here is calculated $\kappa$

and error $\kappa-\kappa^\mathrm{exact}$

Can you see the problem?

Title: Curvature of level-sets of CG1 function


问:
Reduction for NP-hardness
Consider the following optimization problem:
\begin{align}
\text{Min}_{i\neq j\neq s\neq t} |x_i x_j-x_sx_t|\\
s.t: Ax=b\\
x\geq 0;
\end{align}
This problem can be seen as an instance of non convex optimization problems. How can the NP-hardness be shown using a suitable reduction?</p>

Title: Reduction for NP-hardness


问:
speckle interferometry (or other interferometry) software
This might not be the best place to ask, and if its not could someone please point me to where would be better, but here we go anyway.

Does anyone out there know of open source software which can be used for speckle interferometry with normal digital cameras? While speckle seams to be the standard, im happy to use other systems as well, but don't have easy access to lasers etc. It is to measure stress/strain curves of 15x15mm orthotropic samples but I don't need the strain fields within the samples so something less precise than speckle may well be appropriate.

Thanks for any ideas


问:
Network Flow - based optimization in renowned paper appears to have a wrong step
This is probably going to be an easy question for most people and embarass me. I've been reading the Zhang et al. 08 paper on detection-based tracking through network flows: http://iris.usc.edu/outlines/papers/2008/zhang-li-nevatia-cvpr08.pdf. In this paper, equation (11), the authors seem to go from: $-\log(1-\beta_i)\cdot f_i - \log\beta_i \cdot (1-f_i)$ to $C_i \cdot f_i$ where $C_i = \log\frac{\beta_i}{1-\beta_i}$. I've been going over this for about an hour now, shifting logarithms and exponents around but the derivation doesn't seem possible.

What am I missing here?


问:
Order of convergence of Scrodinger eq. with CN scheme
I'm trying to solve numerically the 1-dim time dependent Schrodinger equation using the Crank Nicolson scheme and the Thomas algorithm to solve the tridiagonal matrix.
The physical system consists of a free particle moving in a region of costant (zero) potential, far away from the boundary of my domain.
At this time I am interested in checking the accuracy of the method experimentally, so I am studying the L2 error for different mesh spacing.
The initial wave function and the analytical solution of the Scrhodinger eq. are the same that one can find in Cohen- Tannoudji (complement G1).

The Matlab code is here: http://tny.cz/1bd1bf94.
and the code for the Thomas algorithm (already tested) here: http://tny.cz/83c04868.

The problem is that I don't get the result I should (in particular, if I study the order of convergence experimetally I get inconsistent results and not the second order accuracy in space and time that I expect).
Does anyone has any ideas?
Thank you


问:
Energy anomolies in many body simulation
I am trying to simulate the gravitational interaction between many bodies. I am using a direct PP force calculation and a 4th order symplectic integrator with a variable step size.

The energy of the simulation is calculated by this function:

```<code>real energy(body* system,int n){
  int i,k,j;
  real energy_var=0;
  for(i=0;i<n;i++){
    energy_var += 0.5*system[i].mass*pow(mag(system[i].vel),2);
    for(j=i+1;j<n;j++){
      if(i==j){continue;}
      real position_vect[3]={0,0,0};

      for(k=0;k<3;k++){
    position_vect[k] = system[i].pos[k] - system[j].pos[k];}

      real distance = mag(position_vect);

      energy_var += -Gravitational_constant*system[i].mass*system[j].mass/distance;
    }
  }
  return energy_var;}
</code>```

In a test run with 100 bodies the energy( calculated every 100000's ) appears to be constant apart from one anomaly.

```<code>Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -4.24807e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
</code>```

If it was a numerical or coding error i.e. the timestep being too high then I would not expect the energy to return to the same constant value. Is there a known phenomenon that might cause such an error?


问:
Lapack++ for QR algorithm
I have recently started using Lapack++ which I found convenient for my programming purpose, in general.

Now, I need to solve a matrix using QR algorithm. I've searched the user manual
and I found a LaGenQRFactComplex Class for the purpose. But it need a complex matrix type as input parameter (as per my understanding).

Is there any other way for me to use the same for solving a real valued matrix??

Title: Lapack++ for QR algorithm


问:
Using MINPACK for curve fitting: implementation?
I need to implement a non-linear fitting algorithm in Fortran and chose to use MINPACK's flavor of the Levenberg-Marquardt algorithm as a basis for the least-squares stuff. However, I seem to misunderstand something very basic about how these routines are to be used. My understanding is that I have to call the routine lmdif1 and define a function fcn to be called by the MINPACK routines that calculates the square of the difference between the data to be fitted and the model function; lmdif1 would then return the best solution for the free fitting parameters in an array x.

I have written a small test program (given below) to try this out. The program produces a set of "data" by calculating values from the model function (exp(ax) in this case) with a given parameter a (=-0.5 in this case) and changing them by a small random "error". This dataset is then fitted with the same model function, with a being the free parameter. But it doesn't really work, the fits are usually very bad, unless I almost give the "correct" solution as a starting guess.

Considering that Levenberg-Marquardt is said to be quite robust and seeing how well it works e.g. in Gnuplot, I wonder what I'm doing wrong here. Any ideas? The code below should work when linking to the MINPACK routines.

Tom

```<code>! test lmdif
program testanf
implicit none
integer, parameter :: m=10,n=1,lwa=m*n+5*n+m
integer :: info,iwa(n)
double precision :: x(n),fvec(m),wa(lwa),tol=1d-3
external :: fcn
x(1)=-0.5d0
call lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa)
call translate_info(info,n,tol)
write(*,*) 'a=',x(1)
end program testanf

subroutine fcn(m,n,x,fvec,iflag)
implicit none
! m: number of data points
! n: number of fitting parameters
integer, intent(in) :: m,n,iflag
integer :: i
real, allocatable, save :: t(:),fdata(:)
real :: tmin,tmax,dt,rand
double precision, intent(inout) :: x(n)
double precision, intent(out) :: fvec(m)
character(len=1) :: c
logical, save :: init=.true.
if (init) then
!  create "measured data"
   tmin=0.
   tmax=10.
   dt=(tmax-tmin)/m
   call random_seed(put=(/1,2,3,4,5,6,7,8/))
   allocate(fdata(m),t(m))
   do i=1,m
      t(i)=i*dt
      fdata(i)=fmodel(-0.5d0,t(i))
   end do
   tmax=0.1*maxval(fdata)
   do i=1,m
      call random_number(rand)
      fdata(i)=fdata(i)+(rand-0.5)*tmax ! add measurement errors
      write(10,*) t(i),fdata(i)
   end do
   init=.false.
end if
! minimize the residuals
do i=1,m
   fvec(i)=(fdata(i)-fmodel(x(1),t(i)))**2
end do
if (iflag < 0) deallocate(t,fdata)
return
contains
! model function
  real function fmodel(a,t)
  implicit none
  real, intent(in) :: t
  double precision, intent(in) :: a
  fmodel=exp(a*t)
  return
  end function fmodel
end subroutine fcn
</code>```


问:
What is average pooling
I was reading an article where they have a set of pages. For every page they extract a feature vector of dimention 300. Then they fuse these vectors into one vector by using a method called average pooling: Quoting : " to encode sets of pages into a compact document representation. The main idea consists of performing average pooling of page representations to obtain a single representation of an entire document"

P is the set of pages P = {p1,p2,p3,...,pn} and pi is the vector

F(p) = 1/|P| * Somme(pi) Such that p E P

Title: What is average pooling


问:
PETSc KSP Tutorial ex7 Output
I am learning PETSc following the tutorials of the distribution (version 3.4.1). My question is about the KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); (line 216) of ex7.c (KSP tutorial).
It teaches how to use PCBJACOBI, setting different solvers for each PC block. When I use

mpirun -n 4 ./ex7 -m 20 -ksp_type cg

the output is kind of a mess, like if there wasn't order on the output (and difficult to understand). I just figured out this when I did the serial test

mpirun -n 1 ./ex7 -m 20 -ksp_type cg

Someone can tell me what is going on, since when using PETSC_VIEWER_STDOUT_WORLD the process [0] should print the outputs of the others processes in some (random) order, but without mixing the outputs?

Title: PETSc KSP Tutorial ex7 Output


问:
Malliavin Derivative with Mathematica is it possible?
From mathematica.stackexchange:

Is it possible to define a  Malliavin calculus with Mathematica 9?

Consider a random variables on the Wiener-space $\Omega=\mathcal{C}([0,1])$ of the form $$F=F(\omega)=\displaystyle\int_{0}^{T}h_{t}dW_{t}$$ for some deterministic $h(.)$ in $L^{2}[0,T].$
Here, $\omega$ is a path in the Wiener-space, $\omega \in \Omega$.

We define the Malliavin derivative by

$$
DF=h,\text{ and }D_{t}F=h_{t}.
$$
Let say MalliavinDer is this derivative. We have the following rules

MalliavinDer[F_ + G_] := MalliavinDer[F]+MalliavinDer[G];

MalliavinDer[F_G_]:=MalliavinDer[F]G+ MalliavinDer[G]F

MalliavinDer[f(F)]:=f'[F]MalliavinDer[F].

On the other hand, we have for example:

MalliavinDer[\int_0^T h(t)dWt] := h(t)

MalliavinDer[W_s] := Piecewise[{{1, s <t}, {0, s >t}}]

MalliavinDer[f(W_s] := f'[W_s])Piecewise[{{1, s <t}, {0, s >t}}]


问:
Understanding the meaning of Computational Order of Convergence
I am a postgraduate student with interest in numerical methods for solving nonlinear systems of equations. I have read some papers that discussed about 'computational order of convergence' for some iterative methods for example Increasing the order of convergence of iterative schemes for solving nonlinear systems, but still I can't explain it or implement it practically using MATLAB.

Please, can somebody help me with a good reference to understand how to compute the computational order of convergence for my proposed iterative method?


问:
Patch priority and its effect on Criminsi's Exemplar Based Inpainting
I am trying to implement exemplar based in-painting, as proposed by Criminsi et. al. in its original format (before moving into further optimisations).

I have a few theoretical doubts which I would like to clarify.

I am not sure how the patch priority and the determination of fill-order works to propagate linear structures into the user-selected target region.

From what I understand, Criminsi suggests the following:

Determine the fill front dΩ.

Compute the patch priorities [C(p)D(p)] on patches lying on the fill
front.

Propagate texture into the patches - in order of their patch
priorities.

Update confidence values.

Re-compute fill front dΩ from the remaining unfilled section.

[Repeat these 5 steps till all of the target region is filled].

Now, if this understanding of the algorithm is correct (and, please do correct me if the interpretation is wrong), I have the following doubt:

If I'm always filling patches from the fill front only, and if the
patches are fixed in their spatial location (w.r.t other patches on
the fill-front), then I don't see how the order of filling-in
matters. The patches are filled in by direct sampling of the source
region, which stays constant for all iterations. Hence, I do not
understand how the order of filling will result in any pixel level
difference on the propagated textures. From what I gather, the best match from the source region should remain unaltered, irrespective of the order. What am I missing here?

Based on this video (and assuming it is correct), the patch
filling does not always occur on the fill-front only. How does this
work?  Does this mean we are to compute the priorities for all patches in the target region at once and go on filling them based on their priorities? Then, how does the fill-front determination matter and why does Criminsi determine dΩ at all? I understand Criminsi's approach was a step ahead of the traditional onion-peel technique. But how, exactly?

I really wish to clarify these questions before I delve into 2-3 days of intense coding.

Please help.

Thanks in advance.


问:
GMRES implementations
I am looking for simple GMRES implementations. Unfortunately, I have a few requirements:

should be C or Fortran, Fortran preferred

needs a license that permits inclusion in an MIT/BSD licensed package

needs to deal with complex data

So far, the "templates" code on netlib seemed to fit well. However, it doesn't deal with complex data (although I think changing a few declarations from real to complex could fix that), and also it obviously implements operations better left to BLAS (I think e.g. this code's Givens rotations that'd probably better be drot() calls). Another alternative, PETSc, seems overkill.

I am looking forward to your suggestions. Thank you!

Title: GMRES implementations


问:
Which non-linear conjugate gradient possess finite termination property
There are many variants of non-linear conjugate gradient method available ( Flatcher-Reeves, Polak-Rebiere, Dai-Yuan). In case of minimization of quadratic function when first search direction is chosen to be negative of gradient then all of these algorithms are  equivalent and possess finite termination property.

However, in case where first search direction is arbitrary (not equal to steepest descent direction), all these algorithms are not equivalent anymore. Moreover, their convergence properties varies as well. I am interested to know what happens to finite termination propety in this case? In particular  which non linear CG still possess finite termination property under this condition.


问:
How to estimate if a velocity field is statistically homogeneous?
For a 3D velocity field $\mathbf{u}$ obtained by direct numerical simulation.

Assuming that the field is defined on a periodic domain $\mathcal{P}$ of periodicity $L_x$ in the $x$ direction,  $L_y$ in the $y$ direction and  $L_z$ in the $z$ direction.

Assuming that the simulation is computed from time $t_1$ to $t_2$.

How practically estimate the degree of statistical homogeneity of the velocity field?


问:
Environment for implementing/testing Computer Graphics algorithms
I need to code up a computer graphics algorithm for Surface Registration.

Briefly surface registration is the process of finding "optimal" one-one correspondence
between surfaces, where the meaning of optimal depends on the specific problem.

I do not have much experience with Computer Graphics software, and so I was looking around
the web for a good software environment to code my algorithm up in.

One promising option I hit upon was Blender, which to my delight also has a Python interface to it.

My questions are

Is Blender + Python a good software environment for implementing and testing computer graphics algorithms?

If not, are there better alternatives to it (hopefully using a Pythonic interface)?

Are there any particular points from your experience
which I should keep in mind while using Blender (or whatever package you recommend) while implementing/testing my algorithm?


问:
Issues with CVX package for optimization
I am trying to use the cvx package for optimization. However, I am having some issues with it. I have a variable X which is a matrix but I cannot add $X^{-1}$ in the objective function. What should I do? Is it even possible to use it or it violates the definition of semi definite programming and hence cannot be included</p>


问:
Partitioning Data for Multiple Regression Lines
We're all familiar with traditional least-squares method for constructing a straight line through a set of data points. The question is: suppose I show you a scatter plot which clearly is suggestive of more than one line (suppose two). I need to split the data points into two groups and build a regression line for each.

How do I split the dataset into two so groups so as to make optimal split? What if I allow to leave out some data points entirely from either dataset? Is there a formal name for this problem in research? Are there any known solutions?

Thanks


问:
References on the topic of DEM and XDEM
DEM: discrete element method.
XDEM: extended discrete element method.

For my current project of furnace simulation with granular materials, I am interested in the methods mentioned above. I have not studied these method in the past and now I am looking for general introductory materials. Either monograph or scientific papers would be fine.

Title: References on the topic of DEM and XDEM


问:
Spectra of the energy flux VS transfer spectra?
In direct numerical simulation of turbulence what is the advantage to look at the spectra of the energy flux rather than the transfer spectra?

The transfer spectra corresponds to the triple-velocity correlation coming from the non-linear interaction of the Navier-Stokes equations.
The transfer spectra is generally noted:
\begin{equation}
T(\mathbf{k},t)
\end{equation}
The energy flux through wave number $\mathbf{k}$ is defined from the transfer spectra by
\begin{equation}
\Pi(\mathbf{k},t)=\int^{\infty}_{\mathbf{k}}T(\mathbf{p},t)d\mathbf{p}=-\int^{\mathbf{k}}_{0}T(\mathbf{p},t)d\mathbf{p}
\end{equation}


问:
Good approximate solutions for a MILP problem
The company I work for has been developing an application for real-time control of sewer networks. Every 5 minutes, a MILP problem is built or updated, then solved using Gurobi.

For mid-sized cities, this works well. But for larger cities, the MILP problems get bigger and  of course more time is required to find an acceptable solution.

My question : can anyone suggest me methods to find a good approximate solution instead of THE optimal solution? The method must be able to find a solution fast, so calculation time is an issue as important as the quality of the solution.

Actually, we can solve a MILP problem to its optimal value $F$, or we can find an approximate solution by doing this :

We replace half of the binary variables with continuous variables
varying between 0 and 1;

We find a first optimal solution to this modified problem. Let's define $F_1$ as the optimal value of the objective function.

In the original problem, we set the binary variables to their optimal values from the previous problem;

The approximate solution is found by solving this last problem. Let's define $F_2$ as the optimal value of the objective function.

The quality of the solution depends on the range between $F_1$ and $F_2$ since  $F_1 \leq F \leq F_2$. This method has been giving relatively good results, but I'm pretty sure there may be other methods to find an approximate solution. I googled "approximate solution MILP", but I'm not sure how reliable my findings are, and how fast a sub-optimal solution can be found.

Thanks.


问:
Preconditioner for large size hermitian eigenvalue problems
Basically I try to compute several smallest eigenvalues of some sparse 50k*50k eigenvalue problems using matlab.

$$Ax = \lambda Bx$$

With matlab eigs, it's not as fast as I expected. So I tried some preconditioned routine like lobpcg, however it said "The problem size is too small, relative to the block size. Try using eig or eigs instead".

Anyone can give some suggestion on how can I circumvent this with lobpcg or using some other routine?


问:
Optimizing display of millions of triangles per second
I am working on a web-based molecular visualization program for a research project (e.g., one that can display XYZ and PDB files).  The current options are outdated (for instance, Jmol uses Java).  My advisor uses MDL Chime but almost no browser supports that plugin nowadays.

My idea has been to use modern web tools so the visualization program is future-proof for a while.  I'm using WebGL through Three.js and it seems to be working pretty well.

However, I'm already running up against some performance bottlenecks.  I need to be able to display systems of 500,000 to 2,000,000 atoms in realtime (60 fps).  These systems should be rotatable/zoomable/etc.

The entire system is fed into WebGL as one geometry using a BufferGeometry.

500,000 atoms * 24 triangles per sphere * 60 = 720 million triangles per second

My MBP-retina (GT 650m) chokes on this, so I have no doubt most people's computer would do the same.

However, the MDL Chime program I mentioned before renders this scene fine on my advisor's 10 year old laptop without using the GPU.  What kind of trick could it be using to achieve this performance?  The spheres don't look low-quality; there's no way it could be rendering that many triangles per second, right?


问:
Solving nonlinear wave equation in a dispersive infinite waveguide
I would like to solve a three-dimensional nonlinear wave equation in an infinite cylindrical waveguide numerically. Since the waveguide is dispersive, shocks are less likely to form. Both the convective and the materila nonlinearity (from the state equation) are considered. Is there any open source CFD codes/softwares (like CLAWPACK) to assist me with this task? Also, please help with identifying appropriate schemes.</p>


问:
Time-stepping for coupled nonlinear PDEs
What are good references for time-stepping of the coupled incompressible Navier-Stokes-heat equation (Boussinesq flow),
$$
\begin{cases}
\rho\left(\dot{\mathbf{u}} + \mathbf{u}\cdot\nabla \mathbf{u}\right) = -\nabla p + \mu\Delta \mathbf{u} + \rho(\theta) \mathbf{g},\\
\nabla\cdot \mathbf{u} = 0,
\end{cases}\\
\begin{cases}
\rho c_p \left(\dot{\theta} + \mathbf{u}\cdot\nabla\theta\right) = \nabla(\kappa\nabla\theta) + f
\end{cases}
$$
While the solution of either of the equations is well elaborated upon, and (the pitfalls of) the spatial discretizations are understood by and large, articles on time-stepping are harder to find. It is not clear to me, for example, when it makes sense to treat the nonlinear terms implicitly.

This question may or may not be treated in a more general sense, e.g.,
$$
\dot{u} = \mathcal{L}_1(u, v) + \mathcal{N}_1(u, v),\\
\dot{v} = \mathcal{L}_2(u, v) + \mathcal{N}_2(u, v).
$$
(with linear terms $\mathcal{L}_*$  and nonlinear terms $\mathcal{N}_*$).


问:
Stochastic Collocation for time evolving ODE
For an Stochastic Differential Equation, e.g.,

$$
\frac{du}{dt} = \alpha*\sin(u*t)
$$
where $\alpha$ is normally distributed with nonzero mean, I am trying to use a stochastic collocation approach to find simple statistics such as the mean and the variance.

From what I understand, $\alpha$ is sampled a few times, and the deterministic solution is solved from $t_0$ to $t_F$. Each of these samples is used in a collocation scheme to form the interpolating polynomial. Since in this case, there are multiple time steps between $t_0$ and $t_F$, don't I need to form (possibly) different interpolating polynomials for each of the time points chosen between $t_0$ and $t_F$? Would this not be very slow/inefficient?

Thank you,

Raghu


问:
Find a function's minimizing set of parameters
I have a function defined by three parameters (or variables): $f(a,b,c)$ This function is not explicitly defined but is actually a piece of code which returns the result of fitting a curve $g(a,b,c)$ to a set of $(x,y)$ points in a 2D space:

O = {$(x_1,y_1,w_1), (x_2,y_2,w_2), ...$}

where $w$ are weights associated with each point. The smaller $f$, the better the fit of $g$ with those points. Each one of these points has a weight associated to it.

Each of these three parameters can only assume values from a given discrete set:

A = {$a_1, a_2, a_3, ..., N_a$}

B = {$b_1, b_2, b_3, ..., N_b$}

C = {$c_1, c_2, c_3, ..., N_c$}

Every set $(a_i,b_j,c_k)$ implies a value for $f(a_i,b_j,c_k)$ and as I said, I need to find the $(a_i,b_j,c_k)$ set that returns the minimum value for $f$.

What I do to evaluate $f$ is:

Take a set of $(a_i,b_j,c_k)$ values and construct a $g_{ijk} (a_i,b_j,c_k)$ curve in 2D space.

For each $(x,y)$ point in $O$, find the point in $g_{ijk}$ that is closest to it (euclidean distance)

Multiply that distance by the weight ($0<w<1$) associated to that $(x,y)$ point. This gives me a weighted distance for each point.

After processing all points in the $O$ set, sum all weighted distances.

The sum obtained above is what I call $f(a_i,b_j,c_k)$, ie: a measure of how good $g_{ijk} (a_i,b_j,c_k)$ fitted $O$.

So far I've managed to do this through a brute force approach where I test every possible combination of $a,b,c$ (ie: $N_a*N_b*N_c$ combinations) one after another and at the end I simply isolate the minimum.

Is there a better/more efficient algorithm to perform this minimization, given the restrictions on my function and parameters?

Add

After a bit of research it looks like this falls into the category of combinatorial optimization problems for which methods such as Cross-Entropy seem to be the standard approach.

OwnerDisplayName: Gabriel


问:
Quadratic optimization without any cross terms
I need to solve a quadratic program in Java which minimizes a sum of squares. The problem is that there are constraints to be satisfied, and so a quadratic optimizer must be used. I've been thinking of using JOptimizer for this.

However, the problem is that the $x^T A x$ term in the quadratic program. There are approximately 1000 variables, so the $A$ matrix will be dense with a million terms. What a waste of memory! Is there a way around having to make this big matrix in any other Java library that solves quadratic programs?


问:
What are the governing equations solved in coupled atmosphere-ocean models?
In (hydrostatic) atmospheric general circulation models, for example the so-calle Primitive Equations, consisting of the horizontal momentum equation, the hydrostatic balance, the continuity equation, and the energy balance

$$
\begin{eqnarray}
\nonumber
\frac{\partial}{\partial t}\vec{v} & = &  \vec{v}\times(f+\xi)\vec{e}_z -  \omega\frac{\partial}{\partial p}\vec{v}
                                           -\nabla \left(\Phi+\frac{\vec{v}^2}{2}\right) + \vec{R} \\
\frac{\partial }{\partial p}\Phi & = &  -\frac{1}{\rho} =  -\frac{RT}{p}\\
\nonumber
0 & = & \nabla\cdot\vec{v} +  \frac{\partial}{\partial p} \omega \\
\nonumber
c_p\frac{d}{d t} T & = & \frac{\omega}{\rho} + Q_{rad} + Q_{lat} + \eta + \epsilon .
\end{eqnarray}
$$

are solved numerically. How are these equations modified and complemented, when the atmospheric model is to be coupled to a full ocean general circulation model, what additional effects and equations have to be included?

As the salinity plays an important role in driving the oceanic circulations among other things, for example a balance equation for the salinity, appropriately coupled to the other governing equations will certainly be needed.


问:
Recursive-Multiply Analysis
I am not sure if I can post these kinds of questions on this site, but I was wondering if you could help me with a question that's been bothering me.

We have started learning about analysis of recursive algorithms and I got the gist of it. However there are some questions, like the one I'm going to post, that confuse me a little.
Can you please help me?

Thank you for your time.

Question:

Consider the problem of multiplying two big integers, i.e. integers represented by a large number of bits that cannot be handled directly by the ALU of a single CPU. This type of multiplication has applications in data security where big integers are used in encryption schemes. The elementary-school algorithm for multiplying two n-bit integers has a complexity of . To improve this complexity, let x and y be the two n-bit integers, and use the following algorithm

```<code>Recursive-Multiply(x,y)
  Write  x = x1 * 2^(n/2)+x0  //x1 and x0 are high order and low order n/2 bits
       y = y1 * 2^(n/2)+y0//y1 and y0  are high order and low order n/2 bits
  Compute x1+x0  and y1+y0
  p = Recursive-Multiply (x1+x0,y1+y0)
  x1y1 = Recursive-Multiply (x1,y1)
  x0y0 = Recursive-Multiply (x0,y0)
  Return  x1y1*2^n + (p-x1y1-x0y0)*2^(n/2)+x0y0
</code>```

(a)Explain how the above algorithm works and provides the correct answer.

(b)Write a recurrence relation for the number of basic operations for the above algorithm.

(c)Solve the recurrence relation and show that its complexity is O(n^lg3)

Title: Recursive-Multiply Analysis


问:
Assignment optimization problem
I have to find a solution to a problem and I don't know where to start.

Lets imagine I have a list of n containers and m processes. Each process resides already within a container and have 4 attributes, ex: cpu, memory, IO and network utilization.

The sum of each attribute type of all processes within a host cannot exceed the physical limitation of resources of that host, ex: 5 processes consuming 10% of cpu within a host can be together within a host, but 15 can't.

I would like to find the best (or a good one) solution that assigns those processes to the hosts with the following restrictions:

Use the minimum containers/hosts possibles.

Minimize processes exchanges, from one host to another, if a process is already in one host it would be desirable for it to stay there.

The processes attributes can changes dynamically but there might be a trend based on time or prior behavior.

Maximize a host CPU utilization

Maximize memory utilization

Minimize IO utilization within a host so processes don't conflict to each other.

Minimize network utilization.

I have seen a proposed solution of a simpler version of this problem (just using CPU attribute), using Markov chains evaluating the past values of the CPU utilization.

I first thought applying a multi-objectives Knapsack approach.

Any help please.

Thanks.

Title: Assignment optimization problem


问:
Numerically evaluate 1D inhomogeneous wave equation solution
I am trying to solve the following 1D inhomogeneous wave equation.
Forgive me if I a miss any rigorous mathematical concept.

$$ \frac{\partial^2 u}{\partial x^2} - \frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = \delta(x)f(t) $$

where $c$ is the wave velocity, $f(t)$ any force function and $\delta (x)$ is the dirac delta function. The $\delta(x)$ is used to simulate a point source at $x=0$. I used the Fourier Transform $\mathcal{F}$ first to get to:

$$  \frac{\partial^2 U_{\omega}(\omega,x)}{\partial x^2} + \frac{\omega^2}{c^2} U_{\omega}(\omega,x) =  \delta(x)F(\omega)$$

$$  \frac{\partial^2 U_{\omega}(\omega,x)}{\partial x^2} + k^2 U_{\omega}(\omega,x) =  \delta(x)F(\omega)$$ (A)

with $ k^2 = \frac{\omega^2}{c^2} $ and the subscript $ U_{\omega}(\omega,x) = \mathcal{F_t}(u(t,x)) $ meaning Fourier Transformed from $t$.

Since this is the Helmholtz equation I used the already known Green's function solution:

$$ u(t,x) = \mathcal{F}^{-1}\{ G(x) \ast \left[ F_{\omega} \delta(x) \right] \} $$

$$ u(t,x) = \mathcal{F}^{-1}\{F_{\omega} \frac{j e^{jkx}}{2k}\} $$ (B)

Where $j = \sqrt{-1} $ imaginary unit and $ G(x) = \frac{j e^{jkx}}{2k} $ is the Green's function for the inhomogeneous Helmoltz equation.

Forgive-me for the simple question, but what would be the best way to evaluate B numerically, for any discrete source f(t) in a discrete grid?

I was using FFT and IFFT but I was getting weird results, your guidance will help a lot.

Edit:

Another condition I want the energy vanishes at infinity $ x^+_-$ .


问:
Crank-Nicolson for 2nd- and 4th-order finite differences
I modeled the heat equation,
$$
u_t = au_{xx}
$$
using the common 2nd-order Crank-Nicolson scheme,
$$
\frac{u^{n+1}_i-u^{n}_i}{dt} = \frac{a}{2\,dx}\left(u_{i-1}^{n+1}+u_{i+1}^{n+1}-2u_i^{n+1} + u_{i-1}^n+u_{i+1}^n-2u_i^n\right)
$$
and the standard 4th-order Crank-Nicolson scheme,
$$
\frac{u^{n+1}_i-u^{n}_i}{dt} = \frac{a}{24\,dx}\left(-u_{i-2}^{n+1}+16u_{i-1}^{n+1}-30u_{i}^{n+1}+16u_{i+1}^{n+1}-u_{i+2}^{n+1}\right) \\
   \quad+\frac{a}{24\,dx}\left(-u_{i-2}^{n}+16u_{i-1}^{n}-30u_{i}^{n}+16u_{i+1}^{n}-u_{i+2}^{n}\right)
$$
Outside of figuring out how to solve the pentadiagonal matrix (without resorting to libraries like LAPACK), this was fairly straight-forward.

But my results were strange: the 2nd-order scheme was more accurate to the exact solution (in 3 different test cases!) than the 4th-order scheme. Is there any reason we should expect a 2nd-order accurate scheme to perform better than the 4th-order scheme?

If it matters, $a=1$, $x\in(0,1)$ and $u(0,t)=u(1,t)=0$ for all examples:

$u(x,0)=6\sin(\pi x)$ leading to the analytic solution $u(x,t)=6\sin(\pi x)\exp(-\pi^2t)$

$u(x,0)=12\sin(9\pi x) - 7\sin(4\pi x)$ leading to $u(x,t)=12\sin(9\pi x)\exp(-9\pi^2t) - 7\sin(4\pi x)\exp(-9\pi^2t)$

$u(x,0)=10$ leading to $u(x,t)=\sum_{k=1}^{100}\frac{40}{(2k-1)\pi}\exp\left(-(2k-1)^2\pi^2t\right)\sin\left((2k-1)\pi x\right)$

With this last one, though, both did pretty poorly at the boundaries, but the 2nd-order method was closer to the exact in the range $x\in(0.1,0.9)$.

EDIT

You can see that it's an order of magnitude difference between the two schemes.

two sine term

single sine term

constant 10

EDIT 2

I ran this with different $N$'s and still found the 2nd order version to be more accurate than the 4th order.


问:
analytic or numeric integral of diverging function
I'm trying to carry out the following integral numerically

$$\int_{r_\mathrm{in}}^{r_\mathrm{out}} \Sigma\left(r'\right) \frac{r'}{r} \left( \frac{1}{r-r'}\, E(L) + \frac{1}{r+r'}\, K(L) \right) \mathrm{d}r',$$

where

$$L = \frac{4 r r'}{(r+r')^2},$$

and $K$ and $E$ are the complete elliptical integrals of the first and second kind. $r$, $r_\mathrm{in}$, and $r_\mathrm{out}$ are all $>0$, $r_\mathrm{in}\leq r \leq r_\mathrm{out}$ and $\Sigma(r)$ is some positive, well behaved function (no singularities, or discontinuities, ...).

In the end, I need to evaluate the result numerically for a given $\Sigma(r')$, but my usual python integrators scipy.integrate.quad or scipy.integrate.romberg fail to converge or give unreliable results due to the two singularities:

$1/(r-r')$ in the first summand

$K(L)\rightarrow \infty$ as $r'\rightarrow r$ in the second summand ($L \rightarrow 1$)

Can anyone help me to evaluate this numerically or to rewrite it such that the usual methods don't fail?

OwnerDisplayName: John Smith


问:
Negative viscosity stabilized by fourth order terms
I am trying to solve a "Navier-Stokes"-type problem where the viscosity is negative. Of course this renders the equation unstable and thus I add a fourth order term, so the entire equation becomes:

$$
\frac{dv}{dt} + (v\cdot\nabla)v = - \nabla p+\nu_1 \nabla^2 v - \nu_2 \nabla^4 v
$$
with $\nu_1<0, \nu_2>0$.

Now the fourth order term should stabilize the solution with negative viscosity - I can show this by stability analysis - but I am not able to solve it numerically. I am using the finite element code Fenics.

I am able to solve the system with $\nu_1>0, \nu_2 = 0$ using a simple Picard technique, see e.g. http://fenicsproject.org/documentation/dolfin/dev/python/demo/pde/navier-stokes/python/documentation.html

However, I am not able to make the system stable with negative viscosity. I believe this is a problem with my Finite element approach, but I have no idea how to fix it.

Edit in response to comments: My implmentation is done via fenics. I have tried with both a Lagrange and a CG space. I've done only very simple implementations. I've tried doing a the fourth order derivative on the previous time step's solution. I've also tried doing partial integration twice.

By stable I mean that the field diverges. The velocity field quickly becomes infinite - just as it would if $\nu_2$ is set to zero. It basically behaves just as a Navier Stokes system with negative viscosity and not stabilized by the fourth order term.


问:
Will Sumatra interfere with git?
As far as I understand, Sumatra does some sort of versioning of simulation scripts and output files. If I start to use Sumatra inside a git working directory, can Sumatra end up interfering with git or vice versa?</p>

Title: Will Sumatra interfere with git?


问:
Practical implementation of spatial binning for rectangular range queries
I have a bunch of polygons and a coarse uniform grid. I want to implement two different range queries, for a rectangle aligned with the uniform grid:

Does the rectangle intersect with any polygon at all?

Give me all the polygons which intersect the rectangle.

A simple data-structure would be to have for each cell a list of pointers to the polygons that overlap the cell. This data-structure is appropriate for the first range query. It seems less appropriate for the second range query, because the same polygon may overlap multiple cells, and we don't want to report the same polygon multiple times. One solution would be to use a std::set<polygon_t*> container to collect the pointers to polygons that are delivered by a query. However, let's assume I mistrust std::set in the sense that it will do too many memory allocations even if I reuse the same container for multiple queries. (This mistrust may be unfounded, so I should definitively try this first.)

I wonder whether the problem can be solved more efficiently, if I replace the polygons by their "grid aligned" bounding boxes. During a range query, the knowledge of the lower left point of the bounding box should be sufficient for determining whether the polygon corresponding to the box has already been reported. But even in this case, the complexity of a query returning $m$ polygons won't be $O(m)$, because some work is done for every cell intersecting the bounding box of the polygon. I wonder whether creating three additional lists with the polygons whose bounding box starts in this cell in "x-direction",  "y-direction" and "x- and y-direction" allowing to avoid this additional work would be a good idea. Instead of creating new lists, I could also try to just partition the existing list appropriately.

In general, the coarse uniform grid based data-structures discussed above offer essentially constant lookup time, at the cost of potentially high memory consumption. What is the state of the art regarding such data-structures? Is there a way to avoid the potential memory bloat caused by elements with huge bounding boxes which overlap many grid cells? Also, I don't like the fact that only the simple data-structure can handle polygons directly, while the more efficient data-structures only handle the bounding boxes of the polygons.


问:
How can i find the coordination number with voro++?
Could anybody please help with voro++ cause i am new to this software?My problem is how can i find the coordination  number of an atom ,cause i checked it for a bcc lattice and it gave me 24 for all atoms which is correct cause the wigner-seitz polyhedron of bcc is a truncated octahedron with indeed 24 vertices. But what if i don't know the initial lattice and i want the coordination number for example for bcc is 8?Does anynody knows the option to do this???

Thanks in advance


问:
Finding roots of systems of equations with a Jacobian that is singular everywhere
Given $a\in\mathbb{R}^{mn\times n}$, find a $C\in\mathbb{R}^{n}$, $x\in\mathbb{R}^{m\times n}$ such that
$$
0 = f_{k}(\boldsymbol{C}, \boldsymbol{x}):=\sum_{i=1}^{m} C_{i} \left(\prod_{j=1}^{n} a_{kj}^{x_{i,j}} - 1\right)
$$
for all $k\in\{1,\dots,mn\}$.

(The rows of $a_{kj}$ are vectors of solutions to a system of differential equations, and the goal of finding the roots of the above system is to find invariant quantities of those diff. eq.)

To show that the system has a singular Jacobian, I have a test setup where $a_{k1} = (1-2t_{k})^{-0.5}$, $a_{k2} = (1-2t_{k})^{-1}$ and the rest zero (only $i=1$, but it is still singular if you include the sum). It is in this case easy to check that the Jacobian is singular irregardless of the values of $\boldsymbol{C}$ and $\boldsymbol{x}$. The solution is any $x$ with $2x_{11} - x_{12}=0$ and $C$ anything, and I'm guessing that the Jacobian is singular because this does imply infinite solutions.

I just need to be able to find a couple of them (that are non-trivial such as $C_{i} =0$, and I can go from there. Most root-finding methods I found that can deal with singular Jacobians only mention singularities at the solution, not everywhere. Are there any methods that can tackle this problem? I should also mention that the system should eventually get quite large (order of 1000 equations).


问:
Choosing good basis functions to approximate a Lipschitz function
Let $D = \left\{0, t_1, t_2, \ldots, t_n\right\} \times [0,1]$
and
$$
f:  D\to [0,1],
$$ be a function of time and a one-dimensional space.

There is no analytical formula for $f$, but $f(t_i, \cdot)$ is Lipschitz.

The function $f$ can be computed with arbitrary precision
at a finite number of points $\{x_{t_i, j}\}_{j=0}^m$, $(t_i, j) \in D$.
Alternatively, the Fourier coefficients
$$
\hat f(t_i,k) = \int_0^1 f(t_i,x) e^{2\pi i kx} dx
$$
can be computed for each time slice $f(t_i,\cdot)$, up to a finite order $K$.

I know that the time slices $f(t_i,\cdot)$ look like this

What will be a good choice of basis functions to approximate as accurately as possible $f$?

Is the approximation better with a finite sample $\{f(x_{t_i, j})\}_{j=0}^m$,
or with a truncated Fourier series?

Can the particular basis functions used give bound on the approximation error?

I think a truncated Fourier series is the way to good
as the error is of order $O(1/K)$.

But I'm not too familiar with gibbs phenomenom, so I don't know if the graph of $f$ discard Fourier series at once.

Also, I know Chebyshev polynomial are for non periodic function.

Would they perform better?

Numerically, the easiest thing is to use a polynomial of some order, and use a collocation methods on the finite sample $\{f(x_{t_i, j})\}_{j=0}^m$.

Is this viable?

I think the main problem of this approach is that there is no bound on the approximation error.

My function $f$ is the solution of a global optimization problem,
which is Lipschitz itself
$$
|V\left(f_1\right) - V\left( f_2\right) | \le L \|f_1 - f_2\|_{L^2(D)}.
$$
I know $f$ lives in a compact subset of $L^2(D)$,
so by constructing a fine enough grid, $f$ can be approximate arbitrarily well.


问:
Prove $T_n(x)$ of Chebyshev Polynomial given the recurrence relation
Using the recursion formula for Chebyshev polynomials, show that $T_n(x)$ can be written as
$$T_n(x)=2^{n-1}(x-x_1)(x-x_2)...(x-x_n)$$
where $x_i$ are the $n$ roots of $T_n$

The recurrence relation: $T_0(x)=1$,$T_1(x)=x$, and $T_{n+1}(x)=2xT_n(x)-T_{n-1}(x)$

I have an intuition that I need to use induction,but what is my inductive hypothesis?

Here is my work:

Assume $T_n(x)=2^{n-1}(x-x_1)(x-x_2)...(x-x_n)$

Prove $T_{n+1}(x)=2^{n}(x-x_1)(x-x_2)...(x-x_n)(x-x_{n+1})$

Use the definition of Chebyshev polynomial, $T_n(x) = cos(n\theta),x =cos(\theta)$

$T_0(x)=1$ and $T_1(x)=x$

And then use the recursion formula $T_{n+1}(x)=2xT_n(x)-T_{n-1}(x)$, I get

$T_{n+1}(x) = 2x[2^{n-1}(x-x_1)(x-x_2)...(x-x_n)]-[2^{n-2}(x-x_1)(x-x_2)...(x-x_{n-1})]$

Then

$T_{n+1}(x) = 2x(2^{-1})(x-x_n)[2^n(x-x_1)(x-x_2)...(x-x_{n-1})]-2^{-2}[2^n(x-x_1)(x-x_2)...(x-x_{n-1})]$

Rearrange I get

$T_{n+1}(x) = [2x(2^{-1})(x-x_n)-2^{-2}][2^n(x-x_1)(x-x_2)...(x-x_{n-1})]$

What should I do next(assume what I did above is correct) so that I can get $T_{n+1}(x)=2^{n}(x-x_1)(x-x_2)...(x-x_n)(x-x_{n+1})$ ??


问:
Image processing problems that can be solved using parallel computing?
Many of the problems in Image processing are embaressingly parallel. But I want to know is it worth solving Image processing using parallel computing? Are there real practical problems in image processing that can benefit form parallel computing? If so what are those? I am not taking about mliticore CPUs, instead I am talking about parallelism offered by GPUs, and clusters of GPUs.    </p>


问:
Finite Difference for Hamilton Jacobi Belman
I have hjb equation where $V=V(x,t)$ and $u=u(x,t)$

$V_t + \sup(u) [A(x,u)V_x + B(x,u)V_{xx}]=0$ for $x$ in $[0,1]$ and $t$ in $[0,1]$

I have been able to successfuly resolve it numerically having theese boundary conditions:

$V(0,t)=0$ all t in [0,1]

$V(1,t)=1$ all t in [0,1]

$V(x,1)=x$ all x in [0,1]

I used finite differences to compute $V_t$, $V_x$ and $V_{xx}$ and iterative process at each temporal step to find the optimal control $u^*$ recomputing the solution $V$ until differences among the steps are negligible.

Now i would like to solve it on a infinite space boundary, so letting $x$ go from $-\inf$ to $+\inf$.

I have the condition:
$V(x,1)=max(x,0)$ for all $x$

I thought to use:
$V(x,t) \rightarrow 0$ for $x \rightarrow -\inf$ for all $t$ in $[0,1]$
$V(x,t)/x \rightarrow 1$ for $x \rightarrow \inf$ for all $t$ in $[0,1]$

how do you think is the best way to "discretize" these latter two conditions, if possible?

Many thanks in advance, and i am sorry for the mistakes or imprecisions i can have made.


问:
Increase convergence of non-linear equations resulting from ODEs
I am trying to solve a set of couple ODEs:

$V_l(r) - r W_l(r) -  f1(r) W_l'  = 0\tag 1$

$r^2 h''_l(r) + f2 r h_l'(r) + f3 h_l(r) - f4 U_l(r) = 0 \tag 2$

$\kappa (U_l + h_l) + V_{l+1} + W_{l+1}  = 0\tag 3$

$\kappa (U_{l-2} + h_{l-2}) + f5V_l' + f6W_l' + f7 U'_{l-1} +f8U'_{l-2} = 0 \tag 4$

for $V_l,W_l,h_l,U_l$ and $\kappa$. These functions and the $fn$ are only functions of $r$. Note: These are not the actual ODEs, just a representation of them. There are actually not 4 equations, but $4*l$, where $l$ is the one that appears in the spherical harmonics $Y^l_m$, and $l$ must start at $l=2$.  Also I know what the $fn$ functions are, I do not need to solve for them.

The method I have used to solve this is a spectral approximation. I take know and unknown functions in the ODEs and expand them in Chebyshev polynomials, for example

$f1(r) = \Sigma_{i=0}^{\infty}f1_i T_i[y] - \frac{1}{2} f1_0 $

and

$U_2 = \Sigma_{i=0}^{\infty}U_{2_i} T_i[y] - \frac{1}{2} U_{2_0} $

Applying this to equation (2) I get something like:

$[\kappa * (U_{2_N} + h_{2_N}) + V_{3_N} + W_{3_N} ]*T_N + \ldots - 1/2[\kappa * (U_{2_0} + h_{2_0}) + V_{3_0} + W_{3_0} ]*T_0  $

where $N$ is the maximum order of the Chebyshev polynomial. So if I extract the coefficients of the Chebyshev polynomials I get $4*l*N$, non-linear equations. I say nonlinear because there will be terms like $\kappa * U_{2_N}$, etc. From now on I will refer to the vector of $U_{2_N}, h_{2_N}, \ldots$ as $X$.

In order to solve this non-linear set of equations I am first using the condition that the matrix of equations, $a(\kappa,X)$ needs to satisfy $det(a) = 0$. To makes use of this I am using a singular value decomposition (SVD) to find a $\kappa$ that makes $a$ as singular as possible. To get the values for $X$ I make a new matrix $a'$ in which I remove a row from $a$ and replace it with a "normalization" condition

$\Sigma U_{2_i} T_i = 1$

and then do LinearSolve[a',B]. All the elements of B are zero, expect the one that corresponds to the row that I delete, that element of B is 1.

The main issue I am having is that my value for $\kappa$ does not converge very well. That is, when I increase $l$ and increase $N$ I would expect to get one value of $\kappa$, however at best I can get about 4 digits of convergence, but no matter what I try to do the fifth digit is just random

I've tried messing with the physical model, and trying to increase the accuracy of the calculations, but nothing seems to help. The matrix is poorly conditioned: 1^8 to 1^10. And so I'm sure this is causing some issues.

My attempts at a solution
1st Attempt:

I've tried  to implement a Newton-Raphson method to try and fix the convergence. What I do is start with the matrix $a$, and treat $\kappa$ as a variable like one of the variables in $X$. Then I create the Jacobian

$J = \frac{\partial A_i }{\partial x_j}$
where $A_i$ represents an element of a vector that is made from the rows of my matrix $a$ with the additional equation that comes from the "normalization" condition. I include this so that I have $4*l*N +1 $ equations since now I'm trying so solve for $\kappa$ and $X$ all at once. And $x_j$ has all the elements of $X$ and $\kappa$.

Then I use the values of $\kappa$ and $X$ I got from doing the SVD and row replacement stuff above to use as the initial guess for the Newton's method, and then do

$z = LinearSolve[J, -A]$, and I take $z$ and add it to guess $X$, and repeat. However there is no convergence, in fact the the vector $X$ goes to zero like I'm stuck at the root, basically the $z$ I get is just $-X$. I tried giving LinearSolve the Krylov option, and there the $z$ that I get just flies off to infinity.

I'm not sure where the problem is. I'm thinking it has to do with the fact that my Jacobian might be close to singular ($det (J) \approx.00001$). I read some other stackexchange forums that said to use a QR decomposition instead of LinearSolve, and I tried that and got no convergence there either.

2nd Attempt:
I tried using a QR algorithm for finding eigenvalues on matrix $a(\kappa,X)$. The idea is that I will look for the value of $\kappa$ which give me an eigenvalue as close to zero as possible. Again, this method does not do any better than the SVD stuff I was already using.

Long, long, story short, I really need some help, figuring out just how to make this result converge better. I would want something like 10 digits, but seriously 6 or 7 would be a miracle. What's more I don't care how slow the code gets, if it works it works.

1st Edit
The functions are defined on the domain ($0,R$), where physically zero is the center of a sphere and $R$ is the radius.

The boundary conditions are:
$W_l(r = R) = 0$
$h_{in}(r=R) = h_{out}(r=R)$

and
$h_{l_{out}(r)} = constant*(R/r)^l$

and
$h_l \rightarrow (r/R)^{l+1}$ as $r \rightarrow 0$

$U_l \rightarrow (r/R)^{l+1}$ as $r \rightarrow 0$

$V_l \rightarrow (r/R)^l$ as $r \rightarrow 0$

$W_l \rightarrow (r/R)^l$ as $r \rightarrow 0$


问:
Constraint solving over modular domains
I have a set of constraints over modular domains e.g. $\exists a \in A_i : x \equiv a \pmod{n_i}$ for all $i=0,\ldots,k$

The question is, does such an $x$ exist? I've been pointed to method of successive substitution which ended up with me implementing the following python code.

import fractions
def exist(ruleList, acc = [0,1]): # currentCandidate , currentDomain
  if len(ruleList) == 0:
    return True
  else:
    rule = ruleList[0]
    n = rule[1]/fractions.gcd(rule[1], acc[1])
    newDomain = acc[1] * n
    for candidate in (acc[0] + acc[1] * i for i in xrange(n)):
      if candidate % rule[1] in rule[0]:
        if exist(ruleList[1:], [candidate, newDomain]):
          return True
    return False

rules = [([0,1,2],5),   # A_1,n_1
         ([4,5,6],9),   # A_2,n_2
         ([0,3,6],10)]  # A_3,n_3

print (exists(rules),)  # True

rules = [([0,2,3],5),   # A_1,n_1
         ([1,2,8],9),   # A_2,n_2
         ([4,6,9],10)]  # A_3,n_3

print (exist(rules),)   # False
```

Is their a better way of implementing this or better algorithmic approach?

Title: Constraint solving over modular domains


问:
Sound Waves Simulation in 3D Environment
I want to do a simulation of sound waves including wave propagation, absorption, and reflection in 3D space.

I did some research and I found this question in stackoverflow but it talks about electromagnetic waves not sound waves.

I know I can reimplement the FDTD method for sound waves but how about the sources and does it act like electromagnetic waves? Are there any resources to start with?

Thanks in advance.


问:
N-body simulation for particles
I would like use N-body simulation for measure a torque on an object in a liquid (helium). It's a 2D study. On the site : http://www.browndeertechnology.com/docs/BDT_OpenCL_Tutorial_NBody-rev3.html I can understand the method. But it's only for free object. How can I do for compute particles when they reached fixed wall ? and for compute object that have only one degree of freedom ? Could you help me to start the equations for the simulation ?

Consider red lines with no mass, black disk with no mass. All walls are not computed. Object I would like to compute torque is composed of one part with vacuum and another part with helium (exactly density of helium if there was no walls). Gravity of Earth is perpendicular to the screen, and it is considered perfectly the same (parallel). So like there is only one layer of particles to compute it's possible to forget gravity from Earth I think.

I would like to compute helium :

Density (the object change the density in helium due to the lack of attraction), it's very important

Movement of particles

Torque

I don't compute (superfluid):

friction (helium)

temperature (helium)

walls

The worst is the density, how take it in equation ?

Maybe this sort of program exist ? Could you help me ?

Title: N-body simulation for particles


问:
Size reduction of matrices in dispersion curve calculation
I have an energy dispersion curve obtained from the eigenvalues of
$$E(k) = \text{eig}(T e^{ik} + T^H e^{-ik} + H_0),$$

where $H_0$ and $T$ are $N\times N$ square matrices, $T^H$ is the Hermitian transpose $T$, and $k$ is the wave-vector. So, there are just 2 matrices $T$ and $H_0$ that determine the problem.

I want to reduce the size of matrices $T$ and $H_0$ in a way that the approximation results in the same eigenvalues for a limited energy range (i.e. $E = [1:3]$). Is it possible to construct new matrices $T_1$ and $H_1$ such that $E_1(k) = \text{eig}(T_1e^{ik} + T_1^H e^{-ik} + H_1)$ be close to $E(k)$ for a limited energy range? How?


问:
How to print out a network in peersim?
Peersim is a peer-to-peer network simulator.  I'm simply trying to get the program to print out the network (an edge list is fine).  Presumably this is possible, since there is an in-built writeEdgeList function.

I suspect there is an easy way to do this by modifying this line in the example: include.control io degree  (Section 5 in tutorial 2) which instead prints out a table that lists the number of nodes with certain degrees.  But I don't know what to change "degree" to.

Question: How do I print out an edge list of a network in Peersim?

Title: How to print out a network in peersim?


问:
Texture analysis methods modern survey paper
I want to study the methods of analyzing textured images. So i searched google scholor but only found very old papers

statistical and structural approaches to texture 1979 haralick

Image Texture Analysis Techniques - A Survey (1980)

Texture analysis — a survey  July 1980

Texture Analysis Methods – A Review 1998

Does anyone know of a more recent survey paper ?