问:
Mapping to a computationally less expensive basis when employing Newton's method
I'm looking for advice, or references, for a change of basis to my dependent variables that leads to a less computationally expensive scheme when solving a system of coupled polynomial equations. Below is an explanation of the problem, which will be somewhat involved.

The general idea is this. I have a system of equations I'm solving for the dependent variables $a=(a_1,...,a_N)$, all of them are a function of a parameter $\mu$. That is, my system takes the form ( I will call this method 1)

$$F(a; \mu)=0.$$

For a fixed (prescribed) value of $\mu\in(0,1)$, I solve this system using a Newton-Raphson method (implemented in SUNDIALS), increasing $N$ until my solutions converge to some given precision. For large values of $\mu$, $N$ grows very large and at a certain point becomes very impractical (see my note at the bottom for more about this).

Now, these variables $a_i$ correspond to Fourier coefficients of a physical boundary displacement. In particular, they are Fourier coefficients to two parametric functions $(X(x),Y(x))$, with $X= x+\sum_{n=1}^N a_n \sin n x$, $Y=\sum_{n=1}^N a_n \cos n x$, and $x\in (-\pi,\pi)$.

For a simple example, I have plotted $(X,Y)$, $(x,X(x))$, and $(x,Y(x))$ in black in the figure.

Okay, so there exists another way to formulate this problem, using a different set of variables that leads to a very different set of equations that returns the same information about the boundary displacement. Using ginput in Matlab, I can extract the boundary displacements based on results in the literature that employ this other method (call it method 2), to find a set $(U(x),V(x))$, $(x,U(x))$, and $(x,V(x))$, where the independent variable $x$ is exactly the same as before. These are the red plots in the figure.

Now, we note that $(U,V)$ and $(X,Y)$ give approximately the same plot, but $U$ and $V$ are not nearly as localized in space as $X$ and $Y$. This implies that $U$ and $V$ need much lower spectral resolution than $X$ and $Y$.

This motivates me to look for a change of basis in my original problem, to one that needs lower resolution. In particular, this later approach makes me think that I'm working too hard to get at these solutions, and this becomes particularly apparent as I try to describe the system with higher values of $\mu$.

Now, I realize I should just figure out a mapping between the two methods, but this is very difficult due to the fact that both methods use a series of mappings to get to the final set of equations. So, this information is not available.

My question: Is there any general way, or are there any instances where people have found a way, to do a simple transformation of my variables $a$, to converge for smaller $N$. (I'm not sure if this is the same thing as asking about $preconditioning$ the jacobian of $F$, to allow for more rapid convergence of the Newton-Raphson scheme?). If so, how does one look for this transformation ?

I apologize for the vague question, and please feel free to ask any questions of clarification.

Note: I'm actually solving an ODE with a nontrivial set of coefficients on the highest order derivative, so I'm really solving a system $F(a, t_i)$ for each time step $t_i$. This makes it very important for me to optimize the time taken to solve $F(a)$, as I will do this many, many times.


问:
Bracket Algebra, Straightening Algorithm
My apologies if the question is simple. I need to write a code for straightening algorithm. Which includes defining bracket algebra. I tried to write it in CoCoA-5, but it wasn't possible because CoCoA-5 does not accept brackets like [1,4,6,8] as ring's indeterminates. I did some Googleing and got more confused.
Can anyone introduce a Computer Algebra System which is capable of defining bracket algebra, computing G-Basis, defining ideals and other abilities which are required in Lie algebra and commutative algebra?

Thanks

ContentLicense: CC BY-SA 3.0


问:
Large residual when integrating 2nd order ode close to singularity with SciPy ode / ODEPACK
I am trying to integrate a 2nd order ODE with a singularity at close to the initial condition.
Why do I get large residuals when I plug-in the result of my integration back into the ODE?

The equation is Newcomb's Euler-Lagrange equation from the field of plasma physics.

$$\frac{d}{dr}(f \frac{d\xi}{dr}) - g \xi = 0$$
or as a set of first order ODE's:
$$y_0 = \xi $$
$$y_1 = f \xi' $$
$$y_0' = \frac{y_1}{f} $$
$$y_1' = y_0 g $$

$f$ and $g$ are complicated expressions of magnetic fields and pressure gradients.
However, I even have trouble for simple cases, such as:
$$f=r$$
$$g=-1 + r + \frac{1}{r} $$

The ODE is always singular at $r=0$
From Fobenius expansion the solution close to $r=0$ for this case should be $\xi \propto r^1 \approx 0$ and $\xi' \propto r^0 \approx 1$.

In the code below I set up the problem and integrate with scipy.integrate.ode using the lsoda solver.
I plug the result back into the 2nd order ODE form:
$$f \xi'' + f' \xi' -g \xi $$
I have to use numpy.gradient to calculate $\xi''$.

I get residuals on the edge order of order of $10^2$. Why are they so large?
If I use constants for $f$ and $g$ I get residuals on the order of my step size.

```<code>import numpy as np
from scipy.integrate import ode

# Setup ODE system
def f(r):
   return r

def g(r):
    return -1 + r + 1./r

def der(r, y):
    y_der = np.zeros(2)
    y_der[0] = y[1]/f(r)
    y_der[1] = g(r)*y[0]
    return y_der

#Integrate
integrator = ode(der)
integrator.set_integrator('lsoda')
r_init = 1E-3
init = [r_init, 1.]
integrator.set_initial_value(init, t=r_init)
r = np.linspace(r_init, 1., 100)
results = np.zeros((2, r.size))
results[:,0] = init
for i, position in enumerate(r[1:]):
    integrator.integrate(position)
    results[:, i+1] = integrator.y

# Print residual
dr = r[1] - r[0]
print(f(r)*np.gradient(results[1]/f(r))/dr +
      np.gradient(f(r))/dr*results[1]/f(r) - g(r)*results[0])
</code>```

The output I get looks like the array below. There are large residuals at the edge.

```<code>array([  5.70419912e+02,  -1.59592700e+02,  -9.00242498e-01,
        -1.50959245e-01,  -4.59057150e-02,  ...,
        -1.13773852e-03,  -1.11566938e-03,  -1.07162899e-03,
        -1.79645207e+00])
</code>```


问:
Compute eigenvalues with Arpack
I am using Arpack to compute the eigenvalues of the problem $\lambda Mx = Ax$ with reverse shift method with complex shift. $A$ and $M$ are real, $M$ is symmetric. Then, I use znaupd e zneupd. I use superLU to solve the system. Everything seems to be ok. I wrote a function that takes the eigenvalues and the eigenvectors computed with Arpack and put them in the original problem: $\lambda Mx -Ax = \text{residual}$. The norm of the residual is supposed to be almost zero or very small. Anyway, I obtain an error of $10^{-3}$ that is quite high! Do you have any idea where the problem can come from? I hope it is quite clear.

Thanks,

Nora

Title: Compute eigenvalues with Arpack


问:
What is the Radius of Convergence for analytic functions?
The radius of convergence of any power series can be found by simply using the root test, ratio test etc.

I am confused as to how to find the radius of convergence for an analytic $f$ such as

$f(z)=\frac{4}{(z-1)(z+3)}$.

I can't imagine that I would have to find the power series representation of this, find the closed form, and then use one of the convergence tests. I am fairly certain that the radius of convergence would have to do with the singularities at $1$ and $-3$, however, I can't find a formula for the radius of convergence..


问:
approximation of nonlinear time-dependent system with history
I have two time-dependent coupled equations. One of which is several orders of magnitude more computationally demanding than the other. I am trying to use machine learning to reproduce the behavior of the more expensive equation.

equation 1

```<code>input: c(t), a(t)
output a(t+dt)
</code>```

equation 2

```<code>input: a(t)
output: c(t+dt)
</code>```

So essentially I want to reconstruct the response of equation 2. Keep in mind that internally there are variables in equation 2 which retain 'memory' of the previous states. So the response depends on the history of input.

Any advice on where to start or what methods have been developed for this type of system? OR if there is a more appropriate place to post this?

edit: some more details, this is a multiscale simulation

equation 1 is a simple finite difference equation

$a(x,t+dt) = 2a(x,t) - a(x,t-dt) + \frac{dt^2}{dx^2} \left[ a(x+dx,t)-2a(x,t) + a(x-dx,t) \right] + dt^2 c(x,t)$

for the second part at each x I have a time-dependent set U(t) to propogate to U(t+dt). This propagation depends on an input a(x,t) and produces c(x,t+dt) to be fed back into the first equation. The details of this part are a bit convoluted/involved, but the essential point is that I want to avoid explicitly storing or propagating U (very very very expensive e.g. 10,000+ cpu cores needed)

EDIT2:
A NARX network seems to be able to do almost what I want. However, I have a number of different 'examples' which I want the network to learn from. Maybe the only way I can do it is to stitch everything together into one big (input, output) set?

http://www.mathworks.com/help/nnet/ug/design-time-series-narx-feedback-neural-networks.html


问:
Estimating the second largest eigenvalue
I am currently dealing with the following problem. I am given a matrix $A$ of order $n \times n$ where $n \leq 20.$

The principal $(n-1) \times (n-1)$ matrix of $A$ is symmetric and contains only $\{0,1\}$ as its entries. The $n'$th row and column contains a scaled row sum of this principal sub-matrix. In particular $A$ is not symmetric but its eigenvalues are real.

I would like to determine whether the second largest eigenvalue $\lambda_2$ of $A$ is greater than $2.$

Its perfectly ok for me if the procedure says this is not so (while in fact it is) as long as I do not get false positives.

Right now I am simply computing the eigenvalues of $A$ and checking the condition but it takes a lot of time given that $A$ is not symmetric.

Hence my question is

Are there any better (faster) ways to determine if $\lambda_2 > 2?$

I suppose I could evaluate some determinants of the form $A+\epsilon+2I$ and compare signs, but this approach looks really scary to me.


问:
How to fix time intervals to store data in a stochastic simulation (continous time markov chain)
I am using FORTRAN to implement Gillespie's stochastic simulation algorithm. I would be running many simulations in parallel (both parallel instances with different seed and parallel functions); if I have to store every time interval then a lot of memory would be consumed. So I thought I can have some fixed number of time points where the data would be saved (Allocate space beforehand either dynamically or statically).

How do I decide the size of the time intervals? Should the intervals have a fixed size or it should be based on local dynamics(some regions would have more dynamics)? Is there any way to predict that based on the deterministic function without actually running a deterministic simulation?

I have this vague idea that the tau-leap formula can be used to predict time intervals where n number of a major reaction would occur. But this still doesn't help me predict the number of steps beforehand.


问:
Equilateral triangle based mesh generation by intersection
In work I am currently working on I need to mesh some structure with equilateral triangles to study it using a kind of discrete element method known as spring networks or Lattice model.

To mesh the structure, first I generate the equilateral triangular grid then I do the intersection between the structure and the grid.
My problem lies in the second step, the intersection between the structure to mesh and the grid, I found difficulties to solve it.

Any help or suggestion will be highly appreciated.


问:
HDF5 Many Writers, One External Reader
I have an application where I want to have 3 different machines writing to their own independent HDF5 files at a high rate, but have a way to read data from them externally and combine it. Is HDF5 a good way to do that or should I look into another format?

My main purpose is I have 3 independent packet capturing boxes that will be capturing at about 10Gbps each. They all filtered for specific UDP messages and drop all the rest of the traffic. Inside of these UDP messages there are 3 different types of data, which can be thought of as C structs. These UDP messages need to be written to disk (and possibly formatted differently) in a time series format so that I can search them later. My ultimate goal is to be able to load them into MATLAB quickly, which likely means HDF5 since that's the native format there. I was looking in to PyTables and TsTables, but I'm not quite sure if these do what I want. PyTables doesn't look like it's for time series, and TsTables looks fairly immature.

Title: HDF5 Many Writers, One External Reader


问:
Solving condensate density problem in MATLAB
I want to solve for $n_{0}$ for a fixed value of $n$, lets say $n=1$
 $$ n= n_{0}+ \dfrac{1}{2}\int_{-1/2}^{1/2}dq\left(\dfrac{e_{q}+Un_{0}}{\hbar\omega}-1\right) $$
where $e_{q}=2[1-cos(2\pi q)] $ and $\hbar\omega=(e_{q}^{2}+2Un_{0}e_{q})^{1/2} $

Then I want to plot between  $U$ and $n_{0}$. The following is the code in MATLAB.

```<code> eq = @(q) 2*(1-cos(2*pi*q));
 hq = @(q,U,n0) ((eq(q))^2+2*U*n0*(eq(q))).^(1/2);
 y = @(q,U,n0) (((eq(q))+(U*n0))/hq(q,U,n0))-1;
 a = -0.5;
 b = 0.5;
 v = @(U,n0) quad(@(q) y(q,U,n0),a,b);
 cv =@(U,n0) n0+(0.5*v)-1;
 n0 = 0.1;
 options = optimset('Display', 'iter');
 n = @(U) fsolve(vv,n0,0.1);
 U = 0:0.1:20;
 plot(U,n(U))
</code>```

edit 2 code corrected
edit 1.

Edit 3

But I still facing the following errors:

Error using  ^
Inputs must be a scalar and a square matrix.
To compute elementwise POWER, use POWER (.^) instead.

Error in @(q,n0)((eq(q))^2+2*Un0(eq(q))).^0.5

Error in @(q,n0)(((eq(q))+(U*n0))/hq(q,n0))-1

Error in @(q)y(q,n0)

Error in quad (line 72)
y = f(x, varargin{:});

Error in @(n0)quad(@(q)y(q,n0),a,b)

Error in @(n0)n0+(0.5*v(n0))-1

Error in simulv1 (line 12)
n = fsolve(cv(n0),n0,0.1,options);


问:
Calculation of phonon dispersion
I have a couple of conceptual questions regarding the calculation of phonon dispersions. I wish to calculate the phonon dispersion of Nickel with my lattice vectors given by [111], [1-10] and [11-2] (though my lattice vectors are non-orthogonal to minimize the number of atoms. My lattice vectors are 2.48901587 0 0, 1.244507935 2.155550974 0, 0 0 40). I also have ~15 A vacuum in my cell. I seem to be getting incorrect phonons for this structure but my results are correct if I choose a FCC unit cell with 4 atoms. My questions are summarized below:

Do I have to make any changes to my calculations because I have vacuum in my cell?

I see that calculating force constants (using small displacement method), and hence phonons, requires the creation of supercells. What should they ideally be? Would it be wrong if I just use 1x1x1 cell?

I have attached the results I'm obtaining. The first figure is for the Ni unit cell and the second figure is for the Ni cell with vacuum.

Title: Calculation of phonon dispersion


问:
Resources for viscous behavior in simple FEM
I am working on a simple explicit-integration lumped-mass elastic FEM code which implements CST+DKT triangles (plate+shell) and constant-strain tetrahedra (http://woodem.eu/doc/theory/membrane-element.html, http://woodem.eu/doc/theory/tet4-element.html). The code focuses on contact dynamics, so FEM is there only to model flexible boundaries. I would like to add some kind of viscous damping to the model, and I am looking for some resource which is not overly complicated.

I independently thought I could use the elastic stiffness matrix $\mathbf{K}$ (as in $f=\mathbf{K}u$), scaled by some viscosity factor $\eta'$ (that would be computed from material's $E$ and $\eta$), to compute viscous resisting force as $f_v=-\eta'\mathbf{K}\dot u$. It this formulation something known in literature? Or is it plain wrong?

Thanks for pointers.


问:
What is the difference between SSPRK3 and RK3 time discretization methods?
Here, SSPRK3 refers to third order strong stability preserving Runge-Kutta and RK3 refres to regular third order Runge-Kutta method. The meaning of the method is obvious from the name. However there is a restriction while one is choosing the time step. It is mentioned in "Strong stability preserving high-order time discretization methods" by Gottlieb, Shu , and Tadmor. Since my problem is nonliner systems of hyperbolic PDEs (partial differential equations), related part about my problem is in the section 4.1. The question is that what is the interaction between CFL condition and strong stability preserving coefficient ?</p>

ContentLicense: CC BY-SA 4.0


问:
Alternative to (costly) matrix multiplication
Consider the integral

$$2\pi T = -\frac{1}{2} \int_0^{2\pi} A\frac{\partial B}{\partial \xi} \ d\xi$$

where $$A= \sum_{-N}^N i \ sign(n) \ B_n e^{-in\xi}; \quad \quad B= \sum_{-N}^N \ B_n e^{-in\xi}.$$ These expansions imply

$$2\pi T = \sum |k|B_kB_{-k}$$.

Furthermore,
 $$B_k = \frac{i}{k}\left(-\dot{Y}_k + \sum_{l+j=k}\dot{Y}_lY_j\ j (sign(l)-sign(j))\right),$$

with dots indicating temporal derivatives.

Now, the canonical variables of the system end up being the complex valued variables $Y_{\pm 1},Y_{\pm 2},..,Y_{\pm N}$, and so the Euler-Lagrange equations involve terms like

$$ (1) \quad \quad \frac{\partial T}{\partial Y_m}; \quad \quad \frac{d}{dt}\left(\frac{\partial T}{\partial \dot{Y}_m}\right).$$

Therefore, to derive the governing equations, I must explicitly expand terms like these and write them as some function $F(Y_{\pm 1},Y_{\pm 2},..., \dot{Y}_{\pm 1},...)$

Now, the way I'm currently computing this term is to note

$$T = \sum_{k,l} Q_{kl} \dot{Y}_k\dot{Y}_l,$$

where $$Q_{kl} = \frac{1}{4}\sum_{m=1}^N P_{mk}P_{-ml}$$

and $$P_{mn} = (sign(m-n)-sign(n))Y_{m-n} -\frac{1}{2} Y_mY_{-n}.$$

To set up my governing equations via the Euler-Lagrange equations, I need to compute terms like the ones mentioned above in (1). They both lead to the need to calculate

$$\frac{\partial Q_{kl}}{\partial Y_n} = \sum_{m=1}^N \left(\frac{\partial P_{mk}}{\partial Y_n}P_{-ml}+P_{mk}\frac{\partial P_{-ml}}{\partial Y_n}\right).$$

I need to compute this for all n, so I end up needing to perform a series of 2N+1, $((2N+1),N)\times (N,2N+1))$ matrix multiplications. For large $N$, this computation is absolutely slaughtering my run times, and generally take up +80% of the run time for each time step (depending on the size of N). This is my problem.

The way I've been trying to overcome this is to focus my attention on using the best optimized BLAS libraries to compute the matrix multiplications, and exploiting symmetries of the matrices to try to reduce computation time. However, I think I've milked this line of reasoning for all it's worth without getting too deep into the weeds, so I would like to atleast think about alternatives.

I am very new to numerical methods, and might be completely off base, but it seems like there should be a way to manipulate the FFT here, to avoid doing these matrix multiplications as these Cauchy products are similar to convolutions (ie the pseudo spectral method).

However, I don't see how one can explicitly do this, especially as I need to differential these products with respect to the dependent variables of the system.

If anyone has any insight, it's greatly appreciated.


问:
Piecewise linear fitting of a curve in loglog plot
Iam basically trying to fit a curve similar to that shown in fig (a) with 7 straight line segments (ie. 6 knots) as shown in fig (b) and I need to get the optimum knot positions. This is a tripartite plot in which there are 4 axes (vertical , one inclined at 45 degrees with horizontal, other inclined at 135 degrees with the horizontal ) and the horizontal axis and all are in log scale. Constraints are 3 rd line should be parallel to 45 degree lines, 4 th line should be parallel to the horizontal, 5th line should be parallel to 135 degree lines as shown in the attachment. –  Aiswarya 11 mins ago   .


问:
Creating FEM mesh for image region — what is the most suitable shape function?
I wish to create a FEM mesh to solve an inverse elasticity problem, for an irregular domain. This domain is given by a medical image, so it is discretised and each square on the grid has one scalar value.

I assume the discretized and band-limited nature of the domain can inform my choice of FEM mesh, but have not seen this particular case discussed in my source books.

Consider the 1D case. I know from Lagrange that for a set of regular samples of length n I can construct a polynomial of order n−1. This is one way of expressing the bandwith and information limitations of the image.

My gut tells me that for this reason, I can mesh with squares at the same resolution as the image, using linear shape functions, and capture all the information available in the image. But I lack the math chops to derive this result. Is it correct?


问:
Visualization of solution for a MAC grid
So I implemented the projection method as in Chorin (1969) for the solution of the Navier-Stokes Equations with the Boussinesq Approximation, using a 2nd Order accurate FD scheme. The primitive variables are arranged using a MAC scheme. Basically, what I want is to be able to generate a file (or files) to see the solution on a rendering software (a.k.a, Tecplot, paraView, etc). How can I arrange the data? I mean, from a pseudo-code kind of view.</p>


问:
MCNPX for nuclear simulation
is it legal here to ask this problem?
I need someone help.
I run mcnpx. I try to calculate f4:n using nps:1000. But, it stop directly before culculate f4. There is comment: stack trace terminated abnormally. And there is warning: cross-section file xsdir and bertin does not exist.
Anyone can help what the problem in my work? I appreciated every suggestion. Thanks. Rose</p>

Title: MCNPX for nuclear simulation


问:
Help about Fluid-Fluid coupling techniques
I need a little help and advice with a project I want to do: the idea is to "couple" (I don't know whether I can call it like this) a conservative Navier-Stokes Solver (Fractional-step, 2nd order FDM) with itself. Roughtly speaking, I would like to split the domain in halves so to be able to solve one fluid on one side and another fluid on the other side. The interface is supposed to be fixed at all times, and for the moment I don't need to have different discretizations for each domain. The problem is that one of the fluids is around 1000 times lighter that the other fluid, so (I think) the transmission conditions for the stresses are not going to be as straightforward. Anyhoo, I would like someone to point me in the right track (state-of-the-art literature, jargon, "big names",etc) and shoot some ideas if you like. Ciao & grazie</p>


问:
Am I using the incorrect implementation of the fast Chebyshev transform?
I was told that the fast Chebyshev transform has superior spectral convergence, but I am unable to verify its rumored convergence. I was given plots of its spectral convergence, where the signal's 40th Chebyshev coefficient shrinks to $10^{-16}$ from its largest. The signal was a Gaussian sampled evenly distributed 64 data points. My implementation of the fast Chebyshev transform only yields the 40th Chebyshev coefficient to be $10^{-2}$ smaller than the largest. Unfortunately, I have not been given permission yet to post the plots of the supposed convergence. I wrote code in Python to reproduce the results for myself, to convenience myself of the power of the fast Chebyshev transform. Below is my Python code for implementing the fast Chebyshev transform and sorting the log absolute of the Chebyshev coefficients. Am I using the incorrect implementation of the fast Chebyshev transform?

```<code>def fct(xp,fp): # fast Chebyshev transform
    # size of array
    N = xp.size
    # define Chebyshev grid and scale/translate to fit inside window
    x = xp[0]+(xp[n-1]-xp[0])*(1-np.cos(np.pi*(np.arange(N)+1/2)/N))/2
    # interpolate signal from Chebyshev grid and compute Chebyshev coefficients
    a = ((-1)**(np.arange(N)))*dct(np.interp(x,xp,fp),type=2)
    a[0]=a[0]/2
    return a
def sortnormlogspec(a,ee=1e-100): # normalize data to log absolute
    # log10 abs norm with -100 floor
    b=np.log10(ee+np.abs(a))
    # sort from largest to smallest coefficient
    b=sorted(b,reverse=True)
    return b
</code>```

Below are my plots of the results.


问:
Is there an easy way to read a PetscBag into a python dict?
I'm using a PetscBag to store the input parameters of my program. At some point, I'm going to need to use python to plot these parameters against some output parameters, and dicts along with matplotlib seem to be the way to go with this.

Basically, I want to use PetscBagView() or something like it to write a petsc bag to a file, and then use some python thing (built-in part of the standard library, external package especially for this purpose, whatever works) to read it into a dict. The idea is to make it easy to write matplotlib code that plots various input parameters against certain output parameters.

If dicts are the wrong thing to use for plotting the inputs of several simulations against their outputs, what should I use?

Is there an easy way to read a (binary or ascii) serialized PetscBag into python as a dict ?


问:
Using HMM for speech synthesis
I am trying to make a speech synthesizer using a Hidden Markov Model. I read the wikipedia page and several whitepapers and presentations (such as this one) but I am still not sure how this algorithm can be used to generate waveforms. Can someone explain to me how the algorithm can be applied to this problem?</p>

Title: Using HMM for speech synthesis


问:
numerical analysis of a partial integro-differential equation
I have to numerically solve a nonlinear partial integro-differential equation. This is my equation,

$$\frac{\partial y(x,t)}{\partial t}=\int_{-1/2}^{1/2} \frac{\pi\cos u}{\sin\pi u-\sin\pi x} \frac{\partial y(u,t)}{\partial u}\mathrm du-\sin\,y(x,t)$$

The kernel has a singularity at $x=u$ and vanishes at the boundaries ($u=-1/2,1/2$). Answer $y$ is a function of $x,t$.
At the moment I am not worried about the initial conditions.

Notice, the derivative inside the integral is a first derivation in spacial component and the equation has a non-linearity of $\sin(y(x,t))$. No boundary condition is needed to solve this equation (It's all been taken care of by the kerenl).

What would be the possible numerical recipes to solve this equation?


问:
How to choose (and many) threshold in Eigenface?
I edited my question trying to make it as short and precise.

I am developing a prototype of a facial recognition system for my Graduation Project. I use Eigenface and my main source is the document Turk and Pentland. It is available here: http://www.face-rec.org/algorithms/PCA/jcn.pdf.

My doubts focus on step 4 and 5.

I can not correctly interpret the number of thresholds: If two types of thresholds, or only one (Notice that the text speaks of two types but uses the same symbol). And again, my question is whether this (or these) threshold(s) is unique and global for all person or if each person has their own default.

I understand the steps to be calculated until an matrix O() of classes with weights or weighted. So this matrix O() is of dimension M'x P. Since M' equal to the amount of eigenfaces chosen and P the number of people.

What follows and confuses me. He speaks of two distances: the distance of a class against another, and also from a distance of one face to another. I call it D1 and D2 respectively. NOTE: In the training set there are M images in total, with F = M / P the number of images per person.

I understand that threshold(s) should be chosen empirically. But there must be a way to approximate. I was initially designing a matrix of distances D1() of dimension PxP. Where the row vector D(i) has the distances from the vector average class O(i) to each O(j), j = 1..P. Ie a "all vs all."

Until I came here, and what follows depends on whether I should actually choose a single global threshold for all. Or if I should be chosen for each individual value. Also not if they are 2 types: one for distance classes, and one for distance faces.

I have a theory as could proceed but not so supported by the concepts of Turk:

Stage Pre-Test:

Gender two matrices of distances D1 and D2:
In D1 would be stored distances between classes, and in D2 distances between faces. This basis of the matrices W and A respectively.

Then, as indeed in the training set are P people, taking the F vectors columns D1 for each person and estimate a threshold T1 was in range [Min, Max]. Thus I will have a T1(i), i = 1..P

Separately have a T2 based on the range [Min, Max] out of all the matrix D2. This define is a face or not.

Step Test:

Buid a test set of image with a 1 image for each known person
Itest = {Itest(1) ... Itest(P)}

For every image Itest(i) test:

Calculate the space face Atest = Itest - Imean

Calculate the weight vector Otest = UT * Atest

Calculating distances:

dist1(j) = distance(Otest, O (j)), j = 1..P

Af = project(Otest, U)

dist2 = distance(Atest, Af)

Evaluate recognition:

MinDist = Min(dist1)

For each j = 1..P

If dist2 > T2 then "not is face" else:

If MinDist <= T1(j) then "Subject identified as j" else "subject unidentified"

Then I take account of TFA and TFR and repeat the test process with different threshold values until I find the best approach gives to each person.

Already defined thresholds can put the system into operation unknown images. The algorithm is similar to the test.

I know I get out of "script" of the official documentation but at least this reasoning is the most logical place my head. I wondered if I could give guidance.

EDIT:

i No more to say that has not already been said and that may help clarify things.

Could anyone tell me if I'm okay tackled with my "theory"? I'm moving into my project, and if this is not the right way would appreciate some guidance and does not work and you wrong.


问:
interpolate the circle going through three points using splines
It is known from Geometry that any 3 points determine a circle.  This is the problem of Appolonius

Programming solutions try to approximate using Bezier curves.  Can you not draw conics as Bezier curves?  I thought you could get an exact answer using enough control points.

https://stackoverflow.com/questions/28229695/draw-2-bezier-curves-through-3-points-to-mimic-a-circle-arc

https://stackoverflow.com/questions/4376749/drawing-a-circle-given-three-points [closed]

https://stackoverflow.com/questions/20314306/find-arc-circle-equation-given-three-points-in-space-3d

I would like to be able to input the points $A,B,C$ and return the list of points on the arc $\overline{AB}$ on the circle.

It is possible to compute this formula exactly, such as using the determinant method (or others), but here I am interested in a spline based solution.  This will be useful when I have many many measurements, known to fit a circle, but no one circle exactly fits all of them.

$$ \left| \begin{array}{llll} x^2 + y^2 & x & y & 1 \\
x_1^2 + y_1^2 & x_1 & y_1 & 1 \\
x_2^2 + y_2^2 & x_2 & y_2 & 1 \\
x_3^2 + y_3^2 & x_3 & y_3 & 1 \\\end{array}\right| = 0$$


问:
What library can you use to detect personalities from comments?
I have a bunch of comments from my website that I want to use to build a profile of someone's personality based on the topics they chose to write in and the their specific comments.

I know that this entails natural language processing but I have no idea where to start. I am  a python enthusiast so prefer the library being python.

I'm open to commercial libraries as well but prefer to learn how to do this successfully.


问:
Isolated solving of coupled system of PDE
I would like to solve 3 differential equations for 3 unknowns. So I wrote a MATLAB code, which solves (using the '\' operator) these equation using a linear system of equations (in which the 3 equations are coupled). So the linear system  has a size of (ny*nx*3,ny*nx*3) [nx/ny are the # of nodes in x/y - direction]. I used a finite-difference method.

Solved for $F$:

$\nabla u - \nabla a \nabla F =  \nabla a B_1$

where $u$ is the unknown of a different equation, and $a$, $\alpha$, $t$ and $B$ are parameters.

The parameter $a$ produces very different results of $F$, as shown in the figure below [dotted-black line-$a$ parameter, which is here normalized to the range of $F$; solid red line - solution of $F$]. These results seem to be correct, since they reflect the solution of a benchmark.

Now I would like to isolate these equations and pass the solution of one equation to the the other one (as a parameter on the right-hand-side of the equation). So I store the solution of $u$ from the previous ('coupled') approach and pass it as a parameter to the equation, which is written above. Again I use MATLAB, finite differences and the '\' operator. But this time the linear system of equation only has a size of (ny * nx,ny * nx), since I solve only for one unknown ($F$). I expect the solution of this isolated system to be very similar to the one of the coupled system, since I solve the same equation and have the same values for $u$, $a$, $B$, spatial step size ... !

Unfortunately, the results are very different from the previous ('coupled') approach, if $a$ is very low (see the solid blue line on the lower figure). Especially where $\nabla a$ becomes very low, there is a great discrepancy.
Why is such an approach numerical instable? Is it because of floating point inaccuracy? How could I make the problem work?


问:
Spin-spin correlation functions in the Ising Model Monte Carlo
I'm using the Metropolis algorithm for 2D up to 5D for the Ising Model and I want to compute the spin-spin correlation function.

$$c(r)=<s_is_r>−〈s_i〉〈s_r 〉$$ but I'm not sure how to estimate

$$<s_is_r> = \frac{1}{N}\sum s_is_r$$

What can I do for $s_is_r$ numerically?


问:
Can variational formulations be solved using series solutions?
What I specifically mean is, given some functional $F\left[\mathbf{x}\right]$ which is stationary with respect to $\dot{\mathbf{x}}=f(\mathbf{x})$ and some boundary or initial conditions, can one choose:
$$
\mathbf{x}(t)=\mathbf{x}_0+\mathbf{x}_1t+\mathbf{x}_2t^2+\dots
$$
And substitute this into the functional (or it's first variation) and some how solve for the coefficients?

How would the variations be expressed in this case, would they also be:
$$
\delta\mathbf{x}(t)=\delta\mathbf{x}_0+\delta\mathbf{x}_1t+\delta\mathbf{x}_2t^2+\dots
$$
And then I would collect with respect to them if I was substituting into the first variation?


问:
Potential Reduction and Primal Path following methods
In both the potential reduction and primal path following interior point methods for linear programming, a barrier function is constructed which contains the terms $-\sum \log x_j$ where $x_j$ are the variables. This is to keep the variables from being 0 and thus keep the current solution feasible.

However, in potential reduction, the new direction to move the current solution to is found by solving a linear program which contains as the single inequality constraint $||\mathbf{X^{-1}d}|| \leq \beta$ where $\beta < 1$ and $\mathbf{d}$ is the direction. This constraint is so that the new solution $\mathbf{x}+\mathbf{d}$ remains $> 0$ and thus feasible. $\mathbf{X}$ is a diagonal matrix with the elements of the variable vector $\mathbf{x}$ forming the diagonal.

However in the primal path method (related approach), one can just solve the Lagrangian dual of the barrier function problem and the new solution is guaranteed to be feasible - i.e. the property $\mathbf{||X^{-1}d}|| \leq \beta$ is guaranteed to be satisfied.

Is there some good intuition as to why primal path can get away with just solving the Lagrangian dual (and thus not have to deal with inequalities) while potential reduction cannot? The best I can think of is that the barrier function in the potential reduction algorithm contains the term $\log \mathbf{s'x}$ where $\mathbf{s'x}$ is the duality gap ($\mathbf{s}$ is the slack vector for slacks in the dual version of the problem). This term thus strongly pulls $\mathbf{x+d}$ into a potentially infeasible region. The barrier function for the primal path only has the term $\mathbf{c'x}$ (where $\mathbf{c'x}$ is the objective function of the original problem) in addition to the log terms in the barrier function so it does not exert as strong a pull (as it will not for instance go to $-\infty$ when $\mathbf{c'x}$ goes to 0 unlike $\log \mathbf{s'x}$ which will). However I don't know how to formalize this.


问:
Iteratively finding both left and right eigenvectors for non-symmetric complex matrix
I have a complex, non-Hermitian matrix $\mathbf{A}$, for which I need to find a few eigenvalues and eigenvectors in the generalised eigenvalue problem:

$$\mathbf{A}\cdot \mathbf{x} = \lambda \mathbf{B} \cdot \mathbf{x},$$

with real symmetric positive definite matrix $\mathbf{B}$. I am currently using scipy.sparse.linalg.eigs, which is a Python wrapper around ARPACK routines.

I have successfully found both the eigenvalues and right eigenvectors, however I also need the corresponding left eigenvectors satisfying:

$$ \mathbf{y}\cdot\mathbf{A} = \lambda \mathbf{y}\cdot \mathbf{B}.$$

(According to my notation there is no complex-conjugation required on these left eigenvectors, but obviously undoing a complex-conjugate is easy enough).

As far as I can see, neither scipy nor ARPACK provides a ready-made solution to this problem, not could they be easily modified to give both sets of vectors. The possible solutions I have come across are:

Transpose my matrix $\mathbf{A}$ and make another call to eigs. This is obviously inefficient, and I'm not sure if the returned eigenvalues are guaranteed to have the same values and be in the same order, so matching the left and right eigenvectors may not be reliable.

Use inverse iteration, extrapolating from Numerical Recipes 2nd ed Section 11.7. Solve

$$(\mathbf{A}^T - \lambda\mathbf{B})\cdot \mathbf{y} = \mathbf{b}$$

iteratively for $\mathbf{x}$, starting from some random vector $\mathbf{b}$. The problem with this approach is that it may not find a full set of eigenvectors in the case of degenerate eigenvalues, which I expect to see in some cases.

Are there any better solutions available, or ways that my proposed solutions could be fixed to overcome their (perceived) limitations?


问:
1 D Diffusion equation FDM with different layers
I'm trying to solve this particular equation

$\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \big[D_{i}(x)\frac{\partial u}{\partial x} \big] + S(x,t)$

where the $i$ index denotes different layers. Let's say that $ i \in \{1,2,3\}$.

The interface conditions are:

$-D_1(x)\frac{\partial u}{\partial x}\bigg|_{x=x_1^-} = -D_2(x)\frac{\partial u}{\partial x}\bigg|_{x=x_1^+}$

$-D_2(x)\frac{\partial u}{\partial x}\bigg|_{x=x_2^-} = -D_3(x)\frac{\partial u}{\partial x}\bigg|_{x=x_2^+}$

and at the boundaries:

$\frac{\partial u}{\partial x}\bigg|_{x=0}=0$

$\frac{\partial u}{\partial x}\bigg|_{x=L}=0$

with initial conditions:
$u(x,t_0) = u_0$

where $[0;L]$ denotes the spatial domain over which I want to solve my problem. Every single layer has its own length, let's define with $L_1,L_2$ and $L_3$ and $L = \sum_i L_i$.

I'm using the MOL approach in order to obtain a continous time-discrete space formulation.

Moreover let's assume that the source term $S(x,t)$ is not zero within the layers 1 and 3, but it's null within the second layer.

My current workflow has been:

1 - Divide every signel layer into $N_i$ points with the FDM method

2 - Treat the interfaces points as "common" points; therefore the total number of discretized points will be $N_1+N_2+N_3-2$

3 - At the interfaces points ($N_1$ and $N_1+N_2-1$) enforce the interface BCs. In $N_1$ I am treating that point as part of the layer 1, while in $N_1+N_2-1$ I'm treating that point as part of the layer 2 (is that correct?)

4 - At the end points I'm enforcing the null-flux conditions.

I am able to solve the equation, but I'm experiencing an increase of the overall mass in the system. As far as I know this is a common problem when dealing with conservation of mass problems and FDM; theoretically I could adopt the FVM method in order to avoid this behaviour, but I was courious to understand if I was doing right in this way or wrong. Understand if it was up to the numerical approach or up to some kind of mine bug.

Moreover: the term $S(x,t)$ I'm assuming to be non-zero in $[1;N_1]$ and $[N_1+N_2;N_1+N_2+N_3-2]$

Thanks


问:
state of art non smooth convex optimization
Basically, I am trying to implement non-smooth convex optimization in c++. I am wondering what is the state of art condition of non-smooth convex optimization. For example, what's the best method currently in practice? What is the best solver available? And what is the most recent theoretical progress/ or the most recent review of this topic? Thank you :D</p>


问:
Computing accuracy of my finite difference scheme for uniform grid on a non-uniform grid
I have a ADI finite difference scheme for the 2D Navier-Stokes equations that uses a second order accurate (central) approximation for the advective terms. I am ignoring the diffusive terms for now. I am trying to check the order of accuracy of the scheme by using a test case with analytical solution. For implementation reasons, the grid is non-uniform in the upper half (see footnotes for the reasons for it).

My question is: can I still get second order of accuracy? It is known that for non uniform grid the central scheme reduce to first order (see slides 17-18 here). Note that I use equation 32 without the two rightmost terms, clearly. However  the mentioned slides say that if grid is slowly expanding, i.e., $\Delta x_i = r \Delta x_{i−1}$ with the grid expansion ratio $r$ not too far from unity, the leading order error can then still be essentially second order. Any experience with that? How small has to be $r$? I guess it depends from the problem and the value of the second derivative, since the first truncated term is proportional to it.

PS: The grid is non-uniform in the upper half because it is the only way to have a circular periodic boundary condition without rewriting all the code, i.e., I stretch the grid and I use two lines of cells as the halo for periodicity). The upper half have basically 4 line of cells missing, so I need either to:

stretch the remaining smoothly

have a discontinuity on the grid size at 2 single interfaces.

I guess option (1) is the better, since I make the error smaller locally, right?


问:
Can Gauss-Seidel/SOR (preconditioned?) be applied to a non-diagonally dominant matrix?
After applying finite difference method to a Laplace/Poisson problem always arises a diagonal dominant system of equations that can be solved with Gauss-Seidel or SOR methods. If the original PDE does not have a second order derivative (Laplace operator), but  a first order operator (divergence or curl), the arising system of equations is not diagonal dominant, thus we need to solve it with another solver.

There any way to change it (pivoting, preconditioning, ...) so that it can be solved with Gauss-Seidel or SOR methods?


问:
Determining Youngs Modulus of defined material in FEA (ABAQUS)
I am quite new to FEA but need to determine if I am performing my simulations correctly.

I have a cube of material with defined elastic properties (Youngs Modulus and poison ratio). I perform a simple stretching experiment ($F=-kx$) to determine the effective spring constant $k$ of the material. I can convert $k$ into Youngs modulus $E=k\frac{A}{L}$ where $A$ is area over which force is applied and $L$ is the original length before stretching.

I can stretch it in two ways with the bottom face as an encastre boundary condition

i) Apply a pressure to the top face and monitor displacement

ii) Prescribe a displacement to the top face and monitor the reaction forces

In both cases I apply low enough pressures/displacements so that I don't accidentally go beyond the elastic regime (<0.05% strain in total).

Whenever I do this, I get a Young's modulus that is wrong by about 10%, and sometimes even higher.Should FEA not give an effectively perfect result?

A few notes:

-One thing I am aware of is that I use the original $A$ and not the $A$ at each step of the simulation - so there will be a small error associated with that (perhaps not as small as I expect?)

-My results do not depend on mesh density

-My problem is quasi-static and the solutions do not change if my step set to solve as "Static" or "Dynamic".


问:
Boundary equations for constant right hand side in Poisson equation (FD)
I am setting up to solve a 2D Poisson equation $\nabla^2u = T$ by the finite difference method.

I am using the multigrid method. On the coarsest mesh, where the grid size is 2x2, I want to set up boundary conditions so that the right hand side, $T$, is constant. I do not want a Dirichlet boundary condition, where $u$ is constant, nor a Neumann boundary where the derivative of $u$ is constant.
I do have Dirichlet boundaries. See edit, below.

Here is the base mesh:

Now I set up boundary "ghost" nodes as so, which will hold the constant values of $T$ and influence the 4 "actual" nodes.

So now, using the discretized Poisson equation for $u_{0,0}$,

$\frac{1}{h^2}(-4u_{0,0} + u_{0,1} + u_{1,0} + u_{-1,0} + u_{0,-1}) = T_{0,0}$

I think this is straightforward. The question I have is - what equation(s) to use for the boundary nodes?

To further clarify with a physical correlation, I am attempting something similar to convection, where $T$ represents a temperature analogue, and $u$ an analogue of vorticity, which will then be used to calculate velocities, which then calculates temperature changes.

However, I want these boundary nodes to be set at a constant temperature, in essence "heated from below" where $T_{0,2}$ and $T_{1,2}$ are set to $1.0$ and all other boundary nodes to $0.0$.

Obviously when calculating temperature change for the next time step, it is trivial simply not to update the $T$ values on the boundary and leave them constant. But I'm not sure what to do with the boundaries as part of solving the Poisson equation exactly.

Is it necessary to make the boundary nodes influenced by each other, i.e.
$\frac{1}{h^2}(-4u_{0,-1} + u_{1,-1} + u_{0,0}) = T_{0,-1} = 0.0$

The left hand side of that equation is probably totally wrong, I could forward difference it I suppose in the y-direction, but then do I need to add a corner boundary node, $u_{-1,-1}$ for the x-direction?

Or does a boundary node, being outside of the "fluid" (base mesh) so to speak, not need to be influenced by other nodes also outside the "fluid"? Could I simply write the equation as a 1D forward difference:
$\frac{1}{h^2}(u_{0,-1} - 2u_{0,0} + u_{0,1}) = T_{0,-1} = 0.0$

In summary, although there are plenty of resources describing how to set up Dirichlet and Neumann boundaries, I couldn't find much on how to set up a constant right-hand side boundary. I'm sure it's trivial and I'm missing something obvious, but I appreciate the help, especially since I have no formal training in any scientific computation.

EDIT: In reviewing my question, I realized that I do have an implicit Dirichlet boundary of $u=0$ anywhere outside of the boundary nodes, if I choose a 4x4 mesh, and the equation for the outer nodes is, e.g. $\frac{1}{h^2}(-4u_{0,-1} + u_{-1,-1} + u_{1,-1} + u_{0,0}) = T_{0,-1}$, because $u_{0,-2}$ would be 0 in this case (sorry about the negative indexes, hopefully it's not too confusing).

Which would mean I just stick to central differences. But I still want to define my source term on the outer nodes to get the "heated from below" effect. Should I just redefine $T$ as needed at the beginning of each time step?


问:
Using Centroid decomposition instead of SVD
This paper says centroid decomposition (CD) is an approximation to singular value decomposition (SVD). First I do not understand CD yet, since code is available I just want to try it out how it works but I cannot figure out how to use it, it looks preety much different than SVD. Given a CD how do I get left/right singular vectors and singular values which I can use in my code?  </p>


问:
Finding most efficient route (distance/number of nodes) that uses nodes at least X amount away from another node
I'm trying to find the "looped" route with the lowest value of D/n, where D=Distance, and n=Number_of_points, for a cloud of 3D points. However there are a few conditions.

Conditions:

Each point must have one or more points in the route that are greater than $X$ distance away.

Not every point needs to be passed through.

I will try and explain what I mean:
Imagine a system of delivery routes, whereby you get paid only if the delivery point is greater than $X$ meters away from the pickup point.

So, if you travel point $A \rightarrow B$ where the distance is $X$, you get paid, but you have been paid once in $X$. If you have say 10 equidistant points on a circle with circumference $X$. The distance in a line between them is $X/10$. So you are stopping once every $x/10$ to pick up a new package. Now, once you travel half way around the circle to point 5, you can deliver the package from the starting point and get paid. This means once the route has been completed halfway, you are getting paid every $X/10$ instead. Much more efficient.

Now here is the kicker, my data is not in a neat circle, but rather a cloud of 3D points.

Which algorithm or set of algorithms could be used to determine the most efficient route between the number of points?

Sorry if I horribly explained my problem, from what I understand the problem might be NP-complete but if an approximate solution can be calculated, that would work.


问:
Calculate Integral Using Gauss Jacobi Quadrature or otherwise
I need to integrate the following integral:

\begin{align}
I = \int^z\frac{1-\zeta^2}{(1+\zeta^2)(\zeta-\zeta_l)(1-\zeta_l\zeta)}\prod_{k=2}^{n-1}\left ( \frac{\zeta-z_k}{1-\zeta z_k} \right )^{-\beta_k} d\zeta
\end{align}

along any complex path within the upper half of the unit circle in the argand plane, where $z_k$, $\beta_k$ and $\zeta_l$ have prescribed values.

I am wondering if anyone can suggest a way to robustly calculate this integral with their own code, or comment on if what I am doing below is sensible:

The integral arises in Kirchoff flow problems and is part of the Schwarz Christoffel Parameter problem. I have modified the following code from Toby Driscoll in an attempt to compute my integral. Note that this assumes that qdat which is a matrix of the Gauss Jacobi Quadrature weights is available.

As explained in the comments below, the idea is to refine the interval on which quadrature is applied such that any singularities $z_k$ (other than the end point of integration, are farther away than at least half the integration interval. This is coined "adaptive or compound" Gauss Jacobi Quadrature. Because of the form of the integrand (essentially it is of the type $\int f(x) (1+x)^\alpha(1-x)^\beta dx$) it is amenable to Gauss Jacobi Quadrature.

```<code>function I = modified_hpquad(z1,z2,sing1,z,beta,qdat)
%HPQUAD (not intended for calling directly by the user)
%   Numerical quadrature for the half-plane map.

%   Copyright 1998 by Toby Driscoll.
%   $Id: hpquad.m 298 2009-09-15 14:36:37Z driscoll $

%   HPQUAD(z1,z2,sing1,z,beta,qdat)
%   z1,z2 are vectors of left and right endpoints.  sing1 is a vector of
%   integer indices which label the singularities in z1.  So if sing1(5)
%   = 3, then z1(5) = z(3).  A zero means no singularity.  z is the
%   vector of finite singularities; beta is the vector of associated
%   turning angles.  qdat is quadrature data from SCQDATA.
%
%   Make sure z and beta are column vectors.
%
%   The integral is subdivided, if necessary, so that no singularity
%   lies closer to the left endpoint than 1/2 the length of the
%   integration (sub)interval.

%Some parameters for modified SCT:
sl = 0;

nqpts = size(qdat,1);
% Note: Here n is the total number of *finite* singularities; i.e., the
% number of terms in the product appearing in the integrand.
n = length(z);

bigz = z(:,ones(1,nqpts));
bigbeta = beta(:,ones(1,nqpts));
if isempty(sing1)
  sing1 = zeros(length(z1),1);
end

%I = zeros(size(z1));
I = 0;
%nontriv = find(z1(:)~=z2(:))';
nontriv = find(z1~=z2);

if nontriv == 1

    za = z1;
    zb = z2;
    sng = sing1;

    if sng > n
        sng = 0; %No singularity
    end

    % Allowable integration step, based on nearest singularity.
    dist = min(1,2*min(abs(z([1:sng-1,sng+1:n])-za))/abs(zb-za)); % z([1:sng-1,sng+1:n]) jumps over sing
    zr = za + dist*(zb-za);
    ind = rem(sng+n,n+1)+1;

    % Adjust Gauss-Jacobi nodes and weights to interval.
    nd = ((zr-za)*qdat(:,ind) + zr + za).'/2; % G-J nodes
    wt = ((zr-za)/2) * qdat(:,ind+n+1);     % G-J weights
    %terms = nd(ones(n,1),:) - bigz;
    terms = (nd(ones(n,1),:) - bigz)./(1-bigz.*nd(ones(n,1),:)); %Modified SC
    if any(terms(:)==0)
        % Endpoints are practically coincident.
        I = 0;
    else
        % Use Gauss-Jacobi on first subinterval, if necessary.
        if sng > 0
        terms(sng,:) = terms(sng,:)./abs(terms(sng,:));
        wt = wt*(abs(zr-za)/2)^beta(sng);
        end
        %I(k) = exp(sum(log(terms).*bigbeta,1))*wt;
        I = exp(log((1-nd.^2)./((1+nd.^2).*(nd-sl).*(1-sl.*nd)))+sum(log(terms).*bigbeta,1))*wt;
        while dist < 1
            % Do regular Gaussian quad on other subintervals.
            zl = zr;
            dist = min(1,2*min(abs(z-zl))/abs(zl-zb));
            zr = zl + dist*(zb-zl);
            nd = ((zr-zl)*qdat(:,n+1) + zr + zl).'/2;
            wt = ((zr-zl)/2) * qdat(:,2*n+2);
            %terms = nd(ones(n,1),:) - bigz;
            terms = (nd(ones(n,1),:) - bigz)./(1-bigz.*nd(ones(n,1),:));
            I = I + exp(log((1-nd.^2)./((1+nd.^2).*(nd-sl).*(1-sl.*nd)))+sum(log(terms).*bigbeta,1))*wt;
        end
    end
end
</code>```

The code above is a modified version of a code which does exactly what I want for the integral

\begin{align}
I = \int^z\prod_{k=1}^{n-1}\left ( \zeta-z_k \right )^{-\beta_k} d\zeta
\end{align}

and I have copied that code down for reference here:

```<code>function I = hpquad(z1,z2,varargin)
%HPQUAD (not intended for calling directly by the user)
%   Numerical quadrature for the half-plane map.

%   Copyright 1998 by Toby Driscoll.
%   $Id: hpquad.m 298 2009-09-15 14:36:37Z driscoll $

%   HPQUAD(z1,z2,sing1,z,beta,qdat)
%   z1,z2 are vectors of left and right endpoints.  sing1 is a vector of
%   integer indices which label the singularities in z1.  So if sing1(5)
%   = 3, then z1(5) = z(3).  A zero means no singularity.  z is the
%   vector of finite singularities; beta is the vector of associated
%   turning angles.  qdat is quadrature data from SCQDATA.
%
%   Make sure z and beta are column vectors.
%
%   HPQUAD integrates from a possible singularity at the left end to a
%   regular point at the right.  If both endpoints are singularities,
%   you must break the integral into two pieces and make two calls, or
%   call HPQUAD(z1,z2,sing1,sing2,z,beta,qdat) and accept an automatic
%   choice.
%
%   The integral is subdivided, if necessary, so that no singularity
%   lies closer to the left endpoint than 1/2 the length of the
%   integration (sub)interval.

if nargin==7
  % Break into two pieces with recursive call.
  [sing1,sing2,z,beta,qdat] = deal(varargin{:});
  mid = (z1+z2)/2;
  mid = mid + 1i*abs(mid);
  I1 = hpquad(z1,mid,sing1,z,beta,qdat);
  I2 = hpquad(z2,mid,sing2,z,beta,qdat);
  I = I1-I2;
  return
else
  [sing1,z,beta,qdat] = deal(varargin{:});
end

nqpts = size(qdat,1);
% Note: Here n is the total number of *finite* singularities; i.e., the
% number of terms in the product appearing in the integrand.
n = length(z);
bigz = z(:,ones(1,nqpts));
bigbeta = beta(:,ones(1,nqpts));
if isempty(sing1)
  sing1 = zeros(length(z1),1);
end

I = zeros(size(z1));
nontriv = find(z1(:)~=z2(:))';

for k = nontriv
  za = z1(k);
  zb = z2(k);
  sng = sing1(k);

  % Allowable integration step, based on nearest singularity.
  dist = min(1,2*min(abs(z([1:sng-1,sng+1:n])-za))/abs(zb-za));
%%  if isempty(dist), dist=1; end
  zr = za + dist*(zb-za);
  ind = rem(sng+n,n+1)+1;
  % Adjust Gauss-Jacobi nodes and weights to interval.
  nd = ((zr-za)*qdat(:,ind) + zr + za).'/2; % G-J nodes
  wt = ((zr-za)/2) * qdat(:,ind+n+1);   % G-J weights
  terms = nd(ones(n,1),:) - bigz;
  if any(terms(:)==0)
    % Endpoints are practically coincident.
    I(k) = 0;
  else
    % Use Gauss-Jacobi on first subinterval, if necessary.
    if sng > 0
      terms(sng,:) = terms(sng,:)./abs(terms(sng,:));
      wt = wt*(abs(zr-za)/2)^beta(sng);
    end
    I(k) = exp(sum(log(terms).*bigbeta,1))*wt;
    while dist < 1
      % Do regular Gaussian quad on other subintervals.
      zl = zr;
      dist = min(1,2*min(abs(z-zl))/abs(zl-zb));
      zr = zl + dist*(zb-zl);
      nd = ((zr-zl)*qdat(:,n+1) + zr + zl).'/2;
      wt = ((zr-zl)/2) * qdat(:,2*n+2);
      terms = nd(ones(n,1),:) - bigz;
      I(k) = I(k) + exp(sum(log(terms).*bigbeta,1)) * wt;
    end
  end
end
</code>```


问:
Fast Forward Laplace transform
There are examples for fast numerical inversion of the Laplace transforms. For example here:
http://www.mathworks.com/matlabcentral/fileexchange/32824-numerical-inversion-of-laplace-transforms-in-matlab

My question is: are there implementations of Fast Laplace transform?
For example if one takes and array of values $[a_1 \ldots a_n] $ and transforms then to  $[b_1 \ldots b_n] $.

Title: Fast Forward Laplace transform


问:
Conditions for always positive gradient of heat field in evolutionary thermo-elastic system
I am investigating stability and convergence of series of approximations for coupled thermoelasticity problem yielded by one-step recurrent time-integration scheme.

I've managed to show that the one-step time integration scheme ( OSTIS ) is stable and convergent when the gradient of thermal field is positive (energy of heat of a body at every step is not less than energy on the previous time step of OSTIS). But now I wonder if this is sufficient result or this is unnatural condition?

Generally speaking can someone tell the example of thermo-elastic
process that occurs in nature (or in machinery) where the gradient of
temperature of a body is always positive?

P.S. As I understand this occurs when the system/body is constantly heated by external source, so that more heat comes is than comes out in every time step.

UPDATE:

I need to know when this is true:
\begin{align}
 & \tfrac{1}{2}\Delta {{t}^{-1}}s({{z}^{j+1}}-{{z}^{j}},{{u}^{j+1}}+{{u}^{j}})\ge 0 \\
\end{align}
Here I user the following notation:
\begin{align}
  & \,\theta =\theta (x,t) \\
\end{align} is the temperature field over the body in every time;
\begin{align}
 & \frac{\partial }{\partial t}{{\left. [{{\theta }^{\Delta t}}(t)] \right|}_{{{t}_{j}}}}={{\left. {{[{{\theta }^{\Delta t}}(t)]}^{\prime }} \right|}_{{{t}_{j}}}}={{z}^{j}} \\
\end{align} and

\begin{align}
 & s(\theta ,\xi )=\int_{\Omega }{\theta \xi dx} \\
\end{align}

So basically, my question concerns term

\begin{align} \Delta {{t}^{-1}}[{{z}^{j+1}}-{{z}^{j}}]\ge 0 \end{align}

Are there physical processes where this occurs?

P.S. Meaning of indexes on the top of variables can be understood from here:

\begin{align}
  & {{\theta }^{\Delta t}}(t)={{\theta }^{j+1/2}}+\Delta t[{{\omega }_{j}}(t)-\tfrac{1}{2}]{{z}^{j+1/2}}+\tfrac{1}{2}\Delta {{t}^{2}}{{\omega }_{j}}(t)[{{\omega }_{j}}(t)-1]{{{\dot{z}}}^{j+1/2}}, \\
 & {\theta }'(t)={{z}^{j+1/2}}+\Delta t[{{\omega }_{j}}(t)-\tfrac{1}{2}]{{{\dot{z}}}^{j+1/2}} \\
 & {\theta }''(t)={{{\dot{z}}}^{j+1/2}}, \\
 & \forall t\in [{{t}_{j}},{{t}_{j+1}}]. \\
\end{align}


问:
Jacobian-Free Newton-Krylov vs explicitly forming jacobian in DG
For a given discontinuous galerkin (DG) implementation for Navier-Stokes, targeting 10,000 to 1,000,000 4th order cells in 3D, I'm using PETSc's suite of linear/non-linear solvers on the back-end. It seems to me that there are two broad methods of implicitly solving the discretization. One can explicitly form the Jacobian in a sparse format and use whatever tools are available to solve it, or use a Jacobian-Free Newton-Krylov method that approximates the jacobian-vector products by evaluating the right-hand-side by finite differences. Advantages and disadvantages I see of each:

Explicit Jacobian

Con - Requires a large amount of memory (precisely how much? can I do this with <100GB total memory?)

Con - Requires writing code to evaluate the Jacobian and store it

Pro - Allows the use of black-box preconditioners already present in PETSc

Jacobian-Free Newton-Krylov

Pro - Minimal memory needs

Con - Unknown how often to re-linearize or what finite difference steps to use, how often to restart, etc.

Con - must write own preconditioners.

I'd like some input from those more familiar with these methods as to the advantages and disadvantages of them. Am I right to assume the explicit jacobian would be much faster?


问:
How to solve a nonlinear diffusion equation?
Consider a thin film with a perpendicular applied magnetic field $H_a$ in $z$-axis. The nonlocal relation between $H_a$, the self-field $H_\text{self}$ (generated by the eddy current $J$) and the local magnetic flux density $B_z$ is given by :

$$B_z = \mu(H_\text{self} + H_a)\quad (1)\\
H_\text{self} = \frac{1}{4\pi} \int \frac{J(r')(r-r’)}{(r-r')^2} dr$$
that is Biot-Savart’s law in 2D with $r=(x,y)$.

One defines the local magnetization $m$ as : $J = \nabla\times m$.

As the Biot-Savart's law becomes highly singular when $r = r'$, one uses the Fourier transform. In fact, by inverting (1), one obtains

$$m = 2\mathcal{F}^{-1}[ k^{-1} \mathcal{F}(\mu B_z – H_a)]\quad (2)$$

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ are the Fourier and inverse Fourier transforms and $k$ is the wave vector.

By increasing the applied magnetic field $H_a$, one can write the time evolution of the local magnetization as

$$\frac{dm}{dt} = 2\mathcal{F}^{-1}\left[ k^{-1} \mathcal{F} \left(\mu  \frac{dB_z}{dt} – \frac{dH_a}{dt}\right)\right]\\
\frac{dB_z}{dt} = \nabla\times E \enspace ,$$
the Farady's law.

I think one can solve this diffusion equation using the weak form PDE interface in Comsol but I don’t know how to do it. However, I have no idea how to solve this in Matlab.

Nonlinear diffusion equation


问:
Help understanding the so-called "spectral method"
This is a follow-up question to an answer I read here. $M$ is some hermitian matrix and $V$ an vector.

Since the matrix is hermitian, you could use it as a hamiltonian to propagate it in imaginary time. That is, solve the following system of differential equations:

$$ i\frac{d \vec{V}}{dt}=M\vec{V}$$

The general solution to this is:

$$ V(t)=V_0e^{iMt}$$

Then you take your $\vec{V(t)} \cdot \vec{V(0)}$, fourier transform it, and the height and placement of the peaks will tell you the components along various eigenvectors and their associated eigenvalues. This is sometimes called "the spectral method" in ultrafast atomic physics.

I want to understand this method.

My thoughts on an intuitive level so far: If I take the second equation as given, $V(t)$ is a vector-valued process for $t$. The scalar product $\left<V_{0}e^{iMt},V_{0}\right>$ encodes the matrix $e^{iMt}$ in a scalar-valued object. The fourier transform for a scalar $e^{imt}$ is the delta distribution $\delta_{m}$, up to possibly rescaling by $2\pi$. That is, the fourier transform decodes information of test functions $m\mapsto e^{imt}$. So the foruier transform of $\left<V_{0}e^{iMt},V_{0}\right>$ might somehow decode information about the hermitian matrix $M$.

Is this intuition right so far? If so, any help on how to make this rigorous and how to understand what is going on? A small explicit example would be very helpful.

Any source such as a textbook, lecture notes etc. is most welcome.


问:
Resources for large-scale MILP optimization
With the advent of "big data" applications, different algorithms have to be used to efficiently solve optimization problems, even in the convex case (e.g. the recent success of stochastic gradient descent).

Now I was wondering what happens in cases where you have hard MILP instances with millions ($10^6$-$10^9$) of variables and a similar number of integer constraints (say $10^5$-$10^8$). Also assume that these problems have to be solved really quickly, in a few seconds or less, but a very good approximative answer suffices. In such cases, a standard branch-and-bound algorithm is bound to fail (setting a time limit on standard MILP solvers will usually produce a bad solution or no solution at all).

Are there any resources for dealing with such kind of problems (in the MILP context)? I am looking both for textbooks or for journal articles which describe methods to tackle such problems and especially cases where the large scale of the problem actually increases the power of these methods.

Examples of articles I am looking for, are the following ones by Robin Vujanic:

Vujanic, Robin, et al. "Large scale mixed-integer optimization: A solution method with supply chain applications." Control and Automation (MED), 2014 22nd Mediterranean Conference of. IEEE, 2014.

Vujanic, Robin, et al. "Vanishing duality gap in large scale mixed-integer optimization: a solution method with power system applications." submitted to Journal of Mathematical Programming (2014).


问:
How to curve fit an unknown function?
I have data which can be described by $y=f(x,z)$ where $z$ varies from 170 ~ 154. Now values given by $ks$ are known sample values that equals value given in the table header, $uks$ are unknown samples.
$$
\begin{array}
& z & x_{ks_{1}=0} & x_{ks_{2}=50} & x_{ks_{3} = 60} & \cdots & x_{uks_{1}}& \cdots & x_{uks_{n}}\newline
170 & 3.5 &  5 & 6 &\cdots & 5.3 & \cdots & 12 \newline
168 & 3.7 &  5.1 & 6.6 &\cdots & 5.6 & \cdots & 13.2 \newline
\vdots & \vdots &  \vdots & \vdots &\ddots & \vdots & \ddots & \vdots \newline
154 & 6.8 &  9.52 & 10 &\cdots & 11.5 & \cdots & 26.5 \newline
\end{array}
$$

Plotting values of $y$ vs. $x$ for various $z$, I get this plot:

and on a semi-log plot:

Note the values on $y-axis$ are the values given by the subscripts of $x$ in the table headers.

**My objective is to find the value of $y$ for the corresponding $x_{uks_{1}}, x_{uks_{2}}, \cdots, x_{uks_{n}}$ for different $z$.**There might be alternate methods, such methods are most welcome. And my approach is as follow:

I try to fit a linear function of the form:

$$y = \frac{A_{i}}{ln[x]^{3}} + B_{i},\quad \text{i for each z}$$

this gives me plots as follows:

Since $A_{i}$ and $B_{i}$ is different for $z$, I can  fit a function. Therefore $A_{i} = f(z)$ and $B_{i} = f(z)$. For step 1, I used the following function:

$$A = a_1 z^3 + a_2 z^2 + a_3 z + a_4$$
$$B = b_1 z^2 + b_2 z + b_3$$

Using the correlation:

$$y = \frac{a_1 z^3 + a_2 z^2 + a_3 z + a_4}{ln[x]^{3}} + b_1 z^2 + b_2 z + b_3$$

And substituting different values of $x_{uks}$ I get corresponding $y$, but when I plot them vs. $z$, I get plots of the form:

Now, my questions are as follows:

How am I getting $y$ to be a linear function of $z$?

Is this approach a valid one?

Title: How to curve fit an unknown function?


问:
Anybody using kriging for surface interpolation from huge point datasets?
Is anybody here using kriging for surface interpolation from huge datasets (> 1 million points)?
How large datasets are your datasets and what do you model?
Also what software are you using?</p>


问:
2D wave equation with Mur boundary condition - setting up the RHS and solving (time-steps)
I am trying to solve a 2D wave equation implicitly using FD with central approximations with the following boundary conditions

$$\begin{align}
&u=2\sin\left(\frac{2\pi}{5}t\right)\quad \text{at }x=0\\
&\frac{\partial u}{\partial y}=0 \quad \text{at }y=0\\
&\frac{\partial u}{\partial y}=0 \quad \text{at }y=2\\
&\frac{\partial u}{\partial t}=-\frac{\partial u}{\partial x} \quad \text{at }x=5
\end{align}$$

graphically  represented as

The boundary condition at x=0 generates the wave.
The boundary condition at x=5 refers to Mur boundary condition, i.e. the wave doesn't bounce back, but simply continues to move outside the domain.
The general equation is given as

The Mur boundary condition can be expressed mathematically as

The boundary conditions were evaluated analytically and these are the 7 equations following the entire domain:

Problem- I have managed to write some code, including setting up the [A] matrix (the LHS of the equations), but I have problems with formulating the [b] matrix (RHS) and solving it through the time-steps. I would appreciate if someone show me how to write it for the first equation - then I will be able to do it dor the remaining 6.

```<code>% length, time, height
L = 5; % [m]
h = 2; % [m]

d = 0.25; % space spacing
dt = 0.05; % time increment
M = 100;
T_max = M*dt; % [s]

nx = floor(L/d); % number of <x> samples
ny = floor(h/d); % number of <y> samples
k = floor(T_max/dt) + 1; % number of time samples

% Constants
alpha = dt/d;
r = dt^2/d^2;

% Number of grid points:
N=nx*ny;

% Initialise matrices
U = zeros(N,1);
b = zeros(N,1);
x = zeros(N,1);
A = zeros(N,N);

% Set numbering for the [A] matrix
num = 1;
for i=1:nx
    for j=1:ny
        number(i,j) = num;
        num = num + 1;
    end
end

% [A] matrix - boundary conditions
    % Case 1
    for i=2:nx-1
        for j=2:ny-1
            ii=number(i,j);
            A(ii, number(i, j)) = 1+4*r;
            A(ii, number(i+1, j)) = -r;
            A(ii, number(i-1,j)) = -r;
            A(ii, number(i, j+1)) = -r;
            A(ii, number(i, j-1)) = -r;
        end
    end

    % Case 2
    for i=1
        for j=2:ny-1
            ii=number(i,j);
            A(ii, number(1, j)) = 1+4*r;
            A(ii, number(2, j)) = -r;
            A(ii, number(1, j+1)) = -r;
            A(ii, number(1, j-1)) = -r;
        end
    end

    % Case 3
    for i=nx
        for j=1:ny
            ii=number(i,j);
            A(ii, number(nx, j)) = 1+alpha;
            A(ii, number(nx-1,j)) = 1-alpha;
        end
    end

    % Case 4
    for i=2:nx-1
        for j=1
            ii=number(i,j);
            A(ii, number(i, j)) = 1+4*r;
            A(ii, number(i-1,j)) = -r;
            A(ii, number(i+1,j)) = -r;
            A(ii, number(i, j+1)) = -2*r;
        end
    end

    % Case 5
    for i=1
        for j=1
            ii=number(i,j);
            A(ii, number(i,j)) = 1+4*r;
            A(ii, number(i+1,j)) = -r;
            A(ii, number(i, j+1)) = -2*r;
        end
    end

    % Case 6
    for i=2:nx-1
        for j=ny
            ii=number(i,j);
            A(ii, number(i,j)) = 1+4*r;
            A(ii, number(i-1,j)) = -r;
            A(ii, number(i+1, j)) = -r;
            A(ii, number(i, j-1)) = -2*r;
        end
    end

    % Case 7
    for i=1
        for j=ny
            ii=number(i,j);
            A(ii, number(i,j)) = 1+4*r;
            A(ii, number(i+1, j)) = -r;
            A(ii, number(i, j-1)) = -2*r;
        end
    end

% Inverse of [A] matrix
A_inv = inv(A);

% [b] matrix (RHS)
time = 0;
for t=1:M
    % Computing [b] matrix
    for i=1:nx
        for j=1:ny
            ii = number(i,j);
            % U(1,1) = 2*sin((2*pi/5)*(t*dt));

            % Case 1
            for i=2:nx-1
                for j=2:ny-1

                end
            end
            % or

            for i=1:nx
                for j=1:ny

                end
            end

        end
    end

    % The solution vector
    x = A_inv*b;

    % 2D matrix from the solution vector
</code>```


问:
Efficient Implementation of the Log Sum Exp Log trick for Preventing Underflow erros
So I'm constructing a Gibbs sampler for a moderately sized graphical model, and I started encountering underflow errors (the computed relative frequencies were all very small). To get around this, I've started using the LOG-SUM-EXP-LOG trick found here:

https://math.stackexchange.com/questions/648514/preventing-underflow-log-sum-exp-trick

Now however I have to compute a large number of logs and exponentials. Does anyone know how to do this efficiently in C/C++? I'm considering making a look-up table, but that's a lot of memory...


问:
Minimum effort merging of two sets
I have the following problem.

I have two sequences of elements $A = [a_1,a_2,\cdots,a_n]$ and $B = [b_1,b_2,\cdots,b_m]$. I can build a matrix $D[n \times m]$ where $d_{ij} = d(a_i,b_j)$

My greedy merging algorithm works as follows:

For each iteration it takes the shortest distance pair and merges them.

Then all the distances are recomputed (considering the merged pair)

The process is repeated until each element of the two sequences is merged.

Every time I merge a pair, I sum their distance. The final value of the algorithm is the total sum.

Now, I have to run it on a very large set of sequences.
I would like to have a fast, parallel, algorithm that provides an approximation or an upper bound of the sum.

I looked for eigenvector techniques, but I could not find any suitable algorithm.

Do you have any advice?

Title: Minimum effort merging of two sets


问:
Neural network performs worse when using more input variables
This question is based more on the theory of neural networks than my particular implementation. Therefore I will leave out my code unless requested.

I'm working on a project in C# which can create and train feed forward neural networks. It works well and can train networks to perform some simple functions such as squaring the input or other basic math functions.

My next task was to use this to train a network to predict an outcome of positive or negative based on 39 continuous input variables. I have a dataset comprised of around 7000 training examples.

I did a ROC curve analysis on each input variable and found that some are very good predictors and some are almost useless.

I sorted the variables by AUC, the best was around 78%. I created a network with only this variable as it's input and one hidden layer with one neuron. After training, the network performed better than random but not amazingly.

I next added the second best variable by AUC (75%) as another input, keeping everything else the same. After training the network's fitness score was about 1.5 times better than with one input.

Next I added the third best variable and trained. The networks fitness score was only around 85% of that when it had 2 inputs. I have tried adding more neurons to the hidden layer, running the training for more iterations and many different runs of initial random weights and biases and it seems that no matter what I do the network will perform worse with 3 inputs than 2.

I find this confusing as neural networks are theoretically able to replicate any function. So even if the 3rd variable I introduced was a rubbish predictor then the training should just be able to set the weight of that input to zero and it will perform the same as it did with two inputs.

Questions:

Is my assumption correct that adding inputs to a neural network (and given enough training iterations) it should never perform worse than with less inputs?

Do I require more neurons in the hidden layer when adding more inputs for assumption 1 to hold? Would a good guideline be one neuron in the hidden layer for every input?

Is it a fair assumption that, given everything I have tried, my training code has not been written well enough to train networks to their full capability? (Or in other words, it should be easy for the networks to just ignore the third input, so the problem is more likely to be that my training code is not good enough rather than the networks are not able to easily ignore the third input)?

Is it a good method to create networks with successively more inputs to see which variables improve the fitness the most, and then not use those that don't improve the network? Or should I just create a network with all 39 input variables, train it for a long time and then perform some kind of analysis on the resulting weights to see which inputs are being weighted close to zero and remove them? (as they are essentially being ignored by the network).

Extra details: Each neuron in the hidden layer is connected to every input and the output neuron is connected to every neuron in the hidden layer. All neurons use a sigmoid function with weights and a bias. If the output of the output neuron is >0.8 then this is taken as the network predicting positive, otherwise negative.


问:
Gauss Integration over Zero Order Element
I'm working with the Boundary Element Method and want to integrate an expression over a triangular region. I would like to use Gauss Integration to do this, but I'm having trouble since the triangular element is zero order(1 degree lower than a constant stress triangle). Unlike the higher order elements, the zero order triangle has only one node at its centroid. Normally, Gauss Integration is carried out over a triangular element by defining a natural coordinate system and corresponding shape functions/Jacobin. In the case of the zero order element, however, their is one shape function = 1. Therefore the shape functions can not be expressed in terms of the natural coordinates, and by extension, in terms of the Gauss Points. </p>

<p>Does anyone know how to integrate an expression over a zero order triangular element using Gauss Integration? If not does anyone know a different numerical integration scheme to integrate over a triangular region by just plugging in points and weights. I don't want to write out both integrals and then apply a numerical method to each one. I would also like to avoid defining the limits of integration. </p>


问:
Aliasing question for a spectrally derived Lagrangian
Consider two functions

$$X=\sum_{n=-N}^N i\ sign(n) y_n e^{-in\xi}; \quad \quad Y = \sum_{n=-N}^N y_n e^{-in\xi}.$$

where $y_n$ are the dependent variables of the system, and are functions of time. sign(n) is $1,0,-1$ for n positive, 0 and negative respectively.

Based on the physics of my system, I find that the Lagrangian contains terms like

$$\frac{1}{2\pi}\int_0^{2\pi} \frac{\partial X}{\partial \xi} Y^2 \ d\xi = \sum_{j+k+l=0} |j|y_jy_ky_l$$

with $(j,k,l)\in (-N,N)$.

I then apply the Euler-Lagrange equations, which take the form

$$\frac{d}{dt} \frac{\partial L}{\partial \dot{y_k}} - \frac{\partial L}{\partial y_k}=0 \quad \quad k=\pm 1,\pm 2,... ,\pm N,$$

for $L$ my Larangian to get my governing system of couples ODEs in the dependent variables $y_k$.

These are then solved numerically.

My question is this - does one need to worry about aliasing errors in the sample terms shown above as one normally does for pseudospectral models (even though here my adventure occurs entirely in spectral space)? e.g. should I be cleaving away the highest 1/3 of my modes?


问:
Is there a library which allows for converting 4D implicit surfaces to simplexes?
I have a 4D surface in the form of the spin quadric. I would like to display this surface onto the screen. One of the method which I would like to use to present surface onto the screen is projection. I would like to perform 4D->3D and then 3D->2D projections.

To be able to perform a projection I need to have a set of simplexes which belongs to the surface. Unfortunately, I have only implicit form of the surface and I can't find any library which allows me to convert 4D implicit surface to simplexes.

Do you know any kind of library which allows for converting 4D implicit surface to simplexes?

Thanks in advance.


问:
Tips on improving stability in numerical scheme for non-linear PDE
I am solving a non-linear second order system of PDEs in two variables. The equations are too complicated to write out here, but an essential feature is that there is a propagating wave which then bounces on a boundary.

The problem I have is that the numerics breaks down at the boundary, when the wave reaches this point. I have tried by "trial and error" to just change the way I compute derivatives at this point with different finite difference stenils, but with no luck (pseudospectral methods do not work btw), and I have no idea on how to systematically try to improve the stability and I have no intuition of what can work and what will not work.

Does anyone have any tips on what to try when one encounters such problems? How can I systematically move forward to try to make my numerical scheme stable at this boundary point? Basically, if someone just has a list of ideas to blindly try, that would also be great.

Edit: I should add that the "boundary" can also be interpreted as an origin in polar coordinates, so close to the boundary we have a cylindrical wave scattering at the origin.

Edit2: By request I have added a (reduced version) of the equations I am trying to solve

$\dot{c}=K'$

$\dot{K}=\frac{K^2f}{5r^2}+25(1-f)+\frac{5r}{4}(fc/r)'$

$\dot{f}=-\frac{f^2K'}{5r}$

The "polar coordinate origin" is at r=0, and the wave that propagates is in $c$ and $K$. I have removed many terms to obtain the above equations, and only kept the ones that I think are important. It seems that the instability is due to the $1-f$ term in $\dot{K}$ (if I remove it I get stable evolutions). The instability manifests itself by having $f$ and $K$ diverging with opposite signs at the points close to the origin. Around $r\sim 0$ we expect $K\sim r^4$ and $c\sim r^3$ and $f-1\sim r^2$.

The boundary conditions I put are $K=e^{-t^2}$ at some $r=r_0$, which induces a wave that should bounce at $r=0$. The initial conditions are $f=1$, $c=K=0$.

Some more details about my numerics: I have K and c on different grids (such that a point in $c$ lies between two points in $K$). This removed other instability issues I had. I use second/third order finite difference methods for computing derivatives and for interpolating between the grids. $c$ has a point at the origin, $K$ does not, and $f$ is on the same grid as $K$. At the origin I use the fact that the functions should be even when computing derivatives (by extending the radial coordinate to negative values), but this does not change much.


问:
Quick scheme for separable first-order ODE
I'm trying to integrate an incredibly simple ODE:
$$
y'(x) = -f(y),\quad y(0) = y_0 \ ,
$$
from $x=0$ to $x=1$. This is a decay type of equation, $f$ is the (variable) decay rate and $y$ is the abundance.  Both $y$ and $f(y)$ are positive everywhere.

I'm faced with a couple of constraints.  Sometimes $f(y) \gg 1$, which makes this a stiff problem on the timescales I'm interested in.  I can evaluate $f(y)$ exactly, but it can be costly. So I'd like to minimize the number of times I have to call it.  Luckily, the result doesn't have to be super accurate.  This is a tiny routine that will be called often in a much larger code.  The current method is obviously horrible: we just take a single forward Euler step:
$$
y(1) \approx y(0) + f(y_0) \ .
$$
So literally anything suggested will be an improvement :)

Obviously, I could use an RK scheme, probably with adaptive step sizing, or better yet an implicit scheme.  However, given the simplicity of the equation I thought there may be something more elegant.  For instance, we could rewrite the problem as:

$$
\int_{y_0}^{y_1} \frac{dy}{f(y)} = -1
$$

and attempt to solve for $y_1$, but I'm not sure of the best way to approach this numerically.

Any suggestions appreciated!


问:
How does an unpivoted QR fail to reveal rank?
An unpivoted QR factorization produces a triangular factor $R$. A rank-revealing QR factorization is typically done with column pivoting. My question is, how does an unpivoted QR factorization fail to reveal rank? If a matrix is singular, then a diagonal element of $R$ must be zero or tiny. And this element is almost always (as far as I can tell) the bottom left element.</p>


问:
Implementing a Hill-Type Muscle Model
I'm interested in implementing the muscle model used in Geijtenbeek and Wang et al's work.

Both papers link to the paper by Geyer and Herr, which describes this model:

However, the paper on this model assumes you already know how to model a Hill-Type muscle, and are simply augmenting it and placing many copies of it together in a biologically accurate way.

So in trying to learn how to implement the model itself, I found Gribble's Tutorial which was informative, but a very simplified model of what is used in practice.

I also found Prochazka et al's paper, which seems to be the best detailed reference on implementation, however it describes it in the context of a visual language for describing dynamic systems (Simulink). So I have been learning simulink to do this, but once I reproduce the models it will require some kind of Code Generation tool to convert into c code.

I can do this, but I was wondering if there was a better, more standard way to go about implementing a Hill-Type muscle model because the models presented there seemed overly complicated.

I wouldn't mind a description of how to implement one here as an answer, but a simple reading reference will also suffice.

Title: Implementing a Hill-Type Muscle Model


问:
Monte Carlo sampling of particle system for velocity dependent potential
When sampling configurations of point particles in statistical mechanics by means of (Markov chain) Monte Carlo, we may work with a potential that only depends on the positions of the particles, and to generate new states we assign probabilities to configurations according to the Boltzmann distribution, which (to my understanding) arises from the (temperature dependent) distribution of velocities of the particles. </p>

<p>What can we do when we want to work with a (generalized) potential that also depends on the velocity of the particles? I would assume that in that case we need to consider the velocity of each particle, along with its position, to be part of the state, that statistics alone doesn't suffice anymore. Or does it? How do we assign relative probabilities to different configurations? Can we still use the Boltzmann distribution? </p>


问:
Numerical Implementation of "integrates to some values" type constraint in convex solvers?
I am maximizing a linear functional subject to an integrates to one constraint. More explicitly, my problem is
$$\begin{align}
&\max_{x \in \mathbb{R}^n}\quad c \cdot x\\
&\text{subject to} \int_{a}^{b} \prod_{i=1}^{n} e^{x_i \cdot y^i} dy <=d
\end{align}$$
where $a, b, d \in \mathbb{R}$, $c \in \mathbb{R}^n$, and $n \in \mathbb{N} $ are constants.

I have plugged this into the convex solver SCS, relaxing my integral constraint to a crude Riemann-sum type estimate. Because the sum of exponential constraints is so large in this expression, it takes a very long time to solve. This is my first experience working numerically with a constraint like this, and I would like to increase the speed and/or accuracy of this solution. Is there any more efficient, or computationally tractable, method to implement such an integration constraint like the one above?

I suppose I am hoping for some "convergence in epigraph" type result for approximating this constraint. Is this something that is realistic to hope for?


问:
Optimization of nonlocal stencil-like operator on subset of regular grid
I am trying to optimize the execution time for this particular piece of fortran code.

Details:
i_gc is a (ngpts, 3) array of containing (i,j,k) indices for each grid point. This is a subset of the full NX x NY x NZ grid. The subset is connected in the sense that density of sampling varies, but the points aren't scattered randomly throughout the domain.

Essentially I have a non-local operator on the grid which is a tensor product of 3 one dimensional operators. If I had the data on the full 3d grid it would be pretty easy to write this code without all the if statements in the inner loop. Also, I would like to avoid storing a set {i},{j},{k} indices for each grid point, but maybe this is the way to go?

I am wondering if there is a way to do much better than this? I do have some freedom to reorder the points if needed. Maybe I can exploit that somehow? I suppose I could try something like a Hilbert curve, but I'm not sure how well that could apply to my case where I have an oddly shaped subset of a regular grid.

edit:
mathematically the operation is like this

$Y_{ijk} = \sum^{N^x}_{q=1} M^x_{qi} X_{qjk} +
\sum^{N^y}_{q=1} M^y_{qj} X_{iqk} +
\sum^{N^z}_{q=1} M^z_{qk} X_{ijq} $

where the sums are only over points in the subset (i_gc)

```<code>subroutine apply_lf4(ngpts, i_gc, mx,my,mz,x,y)
  implicit none
  integer, intent(in) :: ngpts
  integer, intent(in), dimension(:,:) :: i_gc
  real(8), intent(in), dimension(:,:) :: mx,my,mz
  real(8), intent(in), dimension(ngpts) :: x
  real(8), intent(out), dimension(ngpts) :: y

  integer :: ig,jg,i1,i2,i3,j1,j2,j3

  y=0.0d0

  do ig=1, ngpts
     i1=i_gc(ig,1); i2=i_gc(ig,2); i3=i_gc(ig,3)
    inner: do jg=1, ngpts
        j1=i_gc(jg,1)
        if (j1==i1) then
           j2=i_gc(jg,2)
           if (j2==i2) then
              j3 = i_gc(jg,3)
              y(ig) = y(ig) + mz(j3,i3)*x(jg)
              cycle inner
           end if
           j3=i_gc(jg,3)
           if (j3==i3) then
              y(ig) = y(ig) + my(j2,i2)*x(jg)
              cycle inner
           end if
        else
           j2 = i_gc(jg,2)
           if (j2 == i2) then
              j3 = i_gc(jg,3)
              if (j3 == i3) then
                 y(ig) = y(ig) + mx(j1,i1)*x(jg)
                 cycle inner
              end if
           end if
        end if
     end do inner
  end do

end subroutine apply_lf4
</code>```


问:
system of coupled nonlinear ODEs with complex coefficients
I am interested in numerically solving the following system of coupled ODEs

$$\left(i-\frac{1}{2\Omega}f_{m,n}\right) \frac{d a_{m,n}(t)}{dt} =E_{m,n}^{\text{kin}}(t) + \left(\omega_{0}+V_{m,n}-\frac{i}{2}\gamma+\frac{i}{2}f_{m,n}\right)a_{m,n}(t)
+\left(g-\frac{i}{2}r\right)|a_{m,n}(t)|^2 a_{m,n}(t)$$

Here $m$ and $n$ are integers, and represent indices on a 2D lattice, as in fact we are dealing with a system of coupled ODEs, one for each lattice point. The lattice is NxN, with N or the order of 30.

The unknown is the complex function $a_{m,n}(t)$. $f_{m,n}$ is also complex, while all the other parameters are real.

The kinetic and potential energy terms are

$$  E_{m,n}^{\text{kin}}(t) = -J [ e^{-i\phi_{m,n}^x}a_{m+1,n}(t)+e^{i\phi_{m,n}^x}a_{m-1,n}(t)
+e^{-i\phi_{m,n}^y}a_{m,n+1}(t)+e^{i\phi_{m,n}^y}a_{m,n-1}(t) ]$$
and
$$V_{m,n} = \frac{1}{2} \kappa b^2 \left[ (m-m_0)^2 + (n-n_0)^2 \right]$$
while
$(\phi_{m,n}^{x},\,\phi_{m,n}^{y})=(-\pi\alpha nb,\,\pi\alpha mb)$, $\alpha = \frac{p}{q}$ and $p,q \in \mathbb{Z}$.

Could someone suggest a suitable numerical scheme for tackling this problem?


问:
Modeling simple laser induced population transfer via adiabatic passage in python
I'm trying to model adiabatic passage between two levels in a three-level atom interacting with two laser fields using Scipy and Numpy..

I'm not sure if my model is wrong due to my incorrectly formulating the equations in python or the Scipy ODE solver picking up errors.

To get the coefficients of each of the three states I use

$$
i\frac{d}{dt}\left( \begin{array}{c}
c_1 \\
c_2 \\
c_3  \end{array} \right) = \left( \begin{array}{ccc}
0 & \Omega_1(t) & 0 \\
\Omega_1(t) & 0 & \Omega_2(t) \\
0 & \Omega_2(t) & 0 \end{array} \right)\left( \begin{array}{c}
c_1 \\
c_2 \\
c_3  \end{array}\right)
$$

Where $\Omega_{1,2}(t)$ give the (slow) time dependence of the amplitudes of the electric fields generated by the two laser beams. The initial conditions are $c_1(-\infty) = 1, c_2(-\infty) = c_3(-\infty) = 0 $.

I assume the Gaussian forms of $\Omega_{1,2}(t)$:

$$
\Omega_{1}(t) = exp(-(t-\tau)^2/T^2)
$$

and

$$
\Omega_{2}(t) = exp(-(t+\tau)^2/T^2)
$$

where $\tau$ is the delay time such that a positive value of $\tau$ means that field 2 is switched on before field 1. $T$ sets the duration of the pulses.

```<code>import numpy as np
import scipy as sp
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from math import exp
from collections import defaultdict
import pylab as p

tau = 2
T = 1

def omega1(t,tau=tau,T=T):
    return exp(-(t - tau)**2/T**2)

def omega2(t,tau=tau,T=T):
    return exp(-(t + tau)**2/T**2)

def g(y, t,tau=tau,T=T):
    c1_i = y[0]
    c2_i = y[1]
    c3_i = y[2]

    g1 = omega1(t, tau, T)*c2_i
    g2 = omega1(t, tau, T)*c1_i + omega2(t)*c3_i
    g3 = omega2(t, tau, T)*c2_i
    return [g1, g2, g3]

# initial conditions
c1_0 = 1
c2_0 = 0
c3_0 = 0

y0 = [c1_0, c2_0, c3_0]

step_size = 1000
t_max = 10.
t_min = 0.
delta_t = (t_max - t_min)/step_size

t = np.linspace(t_max, t_min, step_size) # time grid
soln = odeint(g, y0, t)
c1 = soln[:, 0]
c2 = soln[:, 1]
c3 = soln[:, 2]

# plot results
plt.figure()
plt.plot(t, np.abs(c1), label='c1')
plt.plot(t, np.abs(c2), label='c2')
plt.plot(t, np.abs(c3), label='c3')

plt.xlabel('time')
plt.ylabel('energy level population')
plt.title('initial attempt')
plt.legend(loc=0)
</code>```

The result seems wrong since it shows the coefficients having values above one. Furthermore, I can't seem to find values of $\tau$ and $T$ that yield a non-zero final $c_3$ value.

One issue that might be contributing is that currently my answer doesn't include the multiplication by $i$. I'm not sure how to include that in my python equations.

Thanks for the help!

OwnerDisplayName: Thunder Shiviah


问:
Elemental vs DPLASMA
I want to use one of these two libraries into my C++ project to basically invert a dense matrix (with Cholesky). Of course, I am interested in a distributed environment.

Both libraries seem nice so choosing between them is hard. Does anybody have experience on them? Can anybody provide a comparison for them?

So far I have:

Better support when you have a problem from Elemental, rather than
DPLASMA.

Relevant question: Recommendations for a usable, fast C++ matrix library?

Title: Elemental vs DPLASMA


问:
Is this problem statement good for a GPU?
I am used to using GPU hardware for large scale matrix operations and vectorizing mathematical operations on a continuous space which has been discretized for numerical computation, but this is a completely different problem and I would like to know if it is suitable for a GPU.

Suppose I have a string X of length N and a string Y of length M, such that N < M. I can use, for example, the Python statement if( X in Y ) to find if X is a substring of Y. Now, suppose that I have a string X and 10 million strings {Y_i}. Is it possible (worth it) to put this type of problem on a GPU? Rather than having this be 10 million operations, we could do it in one operation on a GPU, right?

The reason I am asking is that, in reality, I have maybe 10 Billion strings {X_i} and need to, for each string, compare it to every string in {Y_i} and am wondering if I can make it so that finding X in any of {Y_i} will be a GPU kernel for when these sets get incredibly large.


问:
radially averaged power spectrum of a binary image does not look like power law
It has been known that Fourier power spectrum somehow obeys power law therefore the slope of the spectrum can be used to calculate the fractal dimension of an image. Many people have used it for texture analysis.

I wanted to calculate the Fourier fractal dimension of a binary vessel image and have produced radially averaged power spectrum, using the matlab code here: http://www.mathworks.com/matlabcentral/fileexchange/6167-calculation-of-fractal-dimension-of-fractal-surfaces-using-fft

the spectrum I have plotted is like this,

However, I was told that this curve does not look like power law. I want to ask whether the comment is true and why, and is it because I'm using binary image rather than gray-scale?

I was also told that power law plot should be like a straight line so I further extended my experiments and added Gaussian noise to my binary image, like this.

the plot now much more like a power law so the problem is that I did it with binary image, but adding the Gaussian noise destroy the image and the results is not valid.

I want to keep the result I had before but is there anyway to explain my plot to be valid?

Thanks


问:
Numerical integration when solving PDE: Simpsons rule and high frequency noise
I am solving a PDE, and one of the intermediate steps is to numerically integrate a function over a compact interval. The function is represented on a linearly spaced grid. I am using Simpsons rule (so I use the value of the integrand at $x_i$,$x_{i-1}$ and $x_{i-2}$ as well as the value of the integral at $x_{i-2}$ to obtain the integral at $x_i$) , but sometimes I get a zigzag noise after integrating. I know why (in the Simpsons rule the even numered points $x_{2i}$ and the odd $x_{2i+1}$ are sort of obtained separately, and depend a lot on the error of the first point $x_0$ and $x_1$. But the error there is greater since one must use a lower order routine since simpsons rule requires three points).

How can I fix this? What are the state of the art integration routines on linearly spaced grids?


问:
Integration in MATLAB: Complicated integrand
Can I solve the integral given below using Matlab?

$$\frac{C(J)}{C_0} =
\frac{2e^J}{\pi}\int\limits_{0}^{\infty} \frac{e^{\frac{\gamma}{2}\left[1 - \sqrt{\rho}\cos(\theta/2)\right]}}{a_1^2 + a_2^2} [a_1\cos(ZJ-w) + a_2\sin(ZJ - w)]\, dZ$$
where
\begin{align}
&\theta = \arctan(v/u)\\
&u = 1 + \frac{4}{\gamma}\left[1 + \frac{ba + a(1 + Z^2)}{(1 + b)^2 + Z^2}\right]\\
&v = \frac{4Z}{\gamma}\left[1 + \frac{ab}{(1 + b)^2 + Z^2}\right]\\
&\rho = \sqrt{u^2 + v^2}\\
&b = \frac{af}{1 - f}\\
&a_1 = 1 + \sqrt{\rho}\cos(\theta/2) - Z\sqrt{\rho}\sin(\theta/2)\\
&a_2 = Z\left[1 + \sqrt{\rho}\cos(\theta/2)\right] + \sqrt{\rho}\sin(\theta/2)\\
&w = \frac{\gamma}{2}\sqrt{\rho}\sin(\theta/2)
\end{align}

Relation between $Z$ and $J$ is as follows:
$$Z = \frac{af}{1 - f} (J-y) \enspace .$$
$y$ varies in $[0, 1]$; $C_0=1$; $a= 0.065$; $f=0.7$.
I am new to Matlab. Any help will be highly appreciated. I ran the following code (gamma is not y and gamma=m=60, J=1/0.7, O=theta, C=C(J)/C0), but matlab is busy from last 12-13 hours!!

```<code>syms z
m=60;
a = 0.065;
f = 0.7;
b = a*f/(1-f);
u = 1+((4/m)*(1+(b*a+a*(1+z.^2))/((1+b).^2)+z.^2));
v = (4*z/m)*(1+(a*b/(1+b).^2+z.^2));
p = sqrt(u.^2+v.^2);
O = atan(v/u);
a1 = 1+sqrt(p)*cos(O/2)-z*sqrt(p)*sin(O/2);
a2 = z*(1+sqrt(p)*cos(O/2))+sqrt(p)*sin(O/2);
w = (m/2)*sqrt(p)*sin(O/2);
fun = (exp(m/2*(1-sqrt(p)*cos(O/2)))/(a1.^2+a2.^2))*(a1*cos((z/0.7)-w)+a2*sin((z/0.7)-w));
I = int(fun, z, 0, inf);
C = 2*exp(1/0.7)*I/pi;
display(C)
plot(z,C)
xlabel('Z')
ylabel('concentration')
if C>1
    stop
end
</code>```