问:
scipy.integrate.ode ignores boundary conditions
I am trying to solve the 1-dimensional diffusion problem numerically using method of lines:

$$ \frac{\partial c}{\partial t} =D \frac{\partial^2 c}{\partial z^2},$$

where the right hand side is discretised and the left hand side is handled with a standard ode solver.

I have two boundary conditions. One boundary condition dynamically changes with every timestep (dcdt[0]), while the other is always fixed at a constant value (dcdt[-1]). There is also one initial condition (concentration profile) at time zero, which is the line in dark blue in the picture shown.

Unfortunately, what happens is that scipy.integrate.ode ignores both boundary conditions!

```<code>def diffusion(t,c, h):
    dcdt = zeros(Nz)
    dcdt[0] = funcA(h)
    dcdt[-1] = B # a constant
    for i in range(1,Nz-1):
        dcdt[i] = D * (c[i+1] - 2*c[i] + c[i-1]) / (dz*dz)

    return dcdt
</code>```

I use scipy.integrate.ode to integrate the discretisation equation. I have defined the diffusion function to set the boundary condition every time the ode(diffusion) is called.

```<code>ode = integrate.ode(diffusion).set_integrator('vode')
ode.set_initial_value(c_initial, 0).set_f_params(h0)
k = 0

while ode.successful() and k < Nt:
    dcdz = (y[1] - y[0]) / dz
    h[k] = h[k-1] + Q*(h[k-1]**2) / func(h[k-1]) * dt * dcdz
    c_out[:,k] = ode.y
    t[k] = ode.t

    ode.set_f_params(h[k]) # Change the b.c.
    ode.integrate(ode.t + dt)
    k += 1
</code>```

I have checked the following things:

By inserting print h into the diffusion function I have checked that set_f_params has successfully passed h into the solver.

There are no numerical instabilities. I have also checked other solvers dopri5, stiff/non-stiff options, etc.

I have switched the order of functions in the while loop around. It seems that the problem is internal, i.e. ode.integrate(ode.t + dt) solves the equation and ignores the boundary conditions I want to impose before returning the output to me. I'm puzzled why it would do that.

One thing I have thought about is that I could use set_initial_value to set the boundary condition. But as far as I know, I have to provide a full solution vector $c(z,t)$ at a certain time, and ode.integrate is already ignoring my boundary condition at time zero.


问:
Selecting a Finite Element discretization for subsurface flows
Are there guidelines to follow in selecting the right FE discretization. Specifically, for subsurface flow models like Darcy's equation (both single and multi-phase), Richard's equation, etc.

I know one of the most important issues in the fluid mechanics community is the issue of local mass conservation, which the standard classical Galerkin method does not satisfy. That said, there are so many discretizations available out there (e.g., DG, RT0, P2/(P1-P0), least squares FEM, etc). What factors should be taken into account when deciding on a FE formulation to use for my problem? (assuming that they are compatible with the governing equations)

I had these in mind:

1) User implementation - probably the most obvious one. How easy it is to implement in your specific application? Why spend weeks or months attempting to implement a certain formulation when one can achieve the same solution using another formulation that's easier to use?

2) Numerical accuracy - what is the mathematical rate of convergence for each method? This may be important for problems in which there is no analytical solution to compare to.

3) Problem size - that is, how many degrees of freedom each formulation ends up having? If one does not have access to super computers, this may lead to memory issues when running on a laptop or PC

4) Solvers/preconditioners - certain formulations (namely the classical mixed ones) result in saddle-point structures and render standard iterative solvers/preconditioners useless.

5) Parallel performance - How do they scale across multiple processing cores? I think this would be a factor to consider particularly during the assembly, because certain formulations like DG are completely local whereas nodal degrees of freedom may be shared by multiple processors, hence communication overhead. Or would this be more of a user/software implementation issue?

That's all the main ones I could think of for now. Are any of the above valid points to consider? And/or are there other important factors that I missed?

ContentLicense: CC BY-SA 3.0


问:
How to keep velocities in check in molecular dynamics simulation?
I am trying to make a very simple molecular dynamics simulator with Reflective Boundary Conditions. I am assigning the initial positions in a cube randomly while making sure they are not too close to each other and keeping the initial velocities zero. I am using the Lennard-Jones potential and Velocity-Verlet algorithm. But as soon as I start the simulation, Velocities of some of the particles reach very high value and total energy also goes to very high value(not consistent with the initial total energy). Can anyone point out what are the possible reasons for such a behaviour?

Here is my code (in Fortran).

```<code>program mol_dyn_ref
  implicit none
  double precision,allocatable,dimension(:) :: posx,posy,posz,velx,vely,velz,ax,ay,az,tempx,tempy,tempz
  integer,allocatable,dimension(:) :: seed
  double precision :: x,y,z,m,rad,eps,kin,pot,k,p,avgvx,avgvy,avgvz
  double precision :: xinit,yinit,zinit,xlen,ylen,zlen,mindist
  integer :: i,j,st,seedsize,n,iter,t
  double precision :: tot_t,dt,dist,r,Fx,Fy,Fz

  interface
     subroutine acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
       implicit none
       integer :: i,j,n
       double precision,dimension(n) :: posx,posy,posz,ax,ay,az,tempx,tempy,tempz
       double precision :: r,dist,Fx,Fy,Fz,eps,m
     end subroutine acc

     subroutine pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
       implicit none
       double precision,dimension(n) :: posx,posy,posz,velx,vely,velz,ax,ay,az
       integer :: n
       double precision :: tot_t,dt,xlen,ylen,zlen,rad
     end subroutine pos_upd

     subroutine vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)
       implicit none
       integer n
       double precision,dimension(n) :: velx,vely,velz,ax,ay,az,tempx,tempy,tempz
       double precision :: dt
     end subroutine vel_upd
  end interface

  n=50         !Number of particles
  eps=1.0d0
  m=1.0d0
  xlen=20.0d0
  ylen=20.0d0
  zlen=20.0d0
  xinit=0.0d0
  yinit=0.0d0
  zinit=0.0d0
  mindist=1.5d0
  rad=mindist/2.0d0
  tot_t=10.0d0
  dt=0.0001d0

  iter=int(tot_t/dt)

  allocate(posx(n),posy(n),posz(n),velx(n),vely(n),velz(n),ax(n),ay(n),az(n),tempx(n),tempy(n),tempz(n))

  open(100,file="pos.dat",status="replace")
  open(200,file="vel.dat",status="replace")
  open(300,file="acc.dat",status="replace")
  open(400,file="energy.dat",status="replace")

  call random_seed(size=seedsize)
  allocate(seed(seedsize))
  do i=1,seedsize
     call system_clock(st)
     seed(i)=st
  enddo
  call random_seed(put=seed)

  !Assigning initial position to first particle
10 call random_number(x)
  posx(1)=xinit+x*xlen
  call random_number(y)
  posy(1)=yinit+y*ylen
  call random_number(z)
  posz(1)=zinit+z*zlen

  if(posx(1)<rad .OR. posx(1)>xlen-rad) goto 10
  if(posy(1)<rad .OR. posy(1)>ylen-rad) goto 10
  if(posz(1)<rad .OR. posz(1)>zlen-rad) goto 10

  !Assigning initial position
  do i=2,n
20   call random_number(x)
     posx(i)=xinit+x*xlen
     call random_number(y)
     posy(i)=yinit+y*ylen
     call random_number(z)
     posz(i)=zinit+z*zlen
     if(posx(i)<rad .OR. posx(i)>xlen-rad) goto 20
     if(posy(i)<rad .OR. posy(i)>ylen-rad) goto 20
     if(posz(i)<rad .OR. posz(i)>zlen-rad) goto 20
     do j=1,i-1
        if (dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))<mindist) goto 20
     enddo
  enddo
  print*, "Position initialisation finished"

  !Assigning initial velocities
  do i=1,n
     velx(i)=0.d0
     vely(i)=0.d0
     velz(i)=0.d0
  enddo
  print*, "Velocity initialisation finished"

  !Calculating initial acceleration
  do i=1,n
     ax(i)=0
     ay(i)=0
     az(i)=0
     do j=1,n
        if (j==i) cycle
        r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
        ax(i)=ax(i)+ Fx(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        ay(i)=ay(i)+ Fy(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        az(i)=az(i)+ Fz(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
     enddo
  enddo
  print*, "Acceleration initialisation finished."

  !Molecular Dynamics Simulation
  do t=1,iter
     print*, t
     do i=1,n
        write(100,*) posx(i),posy(i),posz(i)
        write(200,*) velx(i),vely(i),velz(i)
        write(300,*) ax(i),ay(i),az(i)
     enddo

     k=kin(velx,vely,velz,m,n)
     p=pot(posx,posy,posz,eps,n)
     write(400,*) t,k,p,k+p

     call pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
     call acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
     call vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)
  enddo
  call system("gnuplot --persist plot.gp")
endprogram mol_dyn_ref

!Updating Acceleration
subroutine acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
  implicit none
  integer :: i,j,n
  double precision,dimension(n) :: posx,posy,posz,ax,ay,az,tempx,tempy,tempz
  double precision :: r,dist,Fx,Fy,Fz,eps,m
  do i=1,n
     tempx(i)=ax(i)
     tempy(i)=ay(i)
     tempz(i)=az(i)
     ax(i)=0.d0
     ay(i)=0.d0
     az(i)=0.d0
     do j=1,n
        if (j==i) cycle
        r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
        ax(i)=ax(i)+ Fx(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        ay(i)=ay(i)+ Fy(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        az(i)=az(i)+ Fz(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
     enddo
  enddo
end subroutine acc

!Updating Position
subroutine pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
  implicit none

  integer n,i
  double precision,dimension(n) :: posx,posy,posz,velx,vely,velz,ax,ay,az
  double precision :: dt,xlen,ylen,zlen,rad
  do i=1,n
     posx(i)=posx(i) + velx(i)*dt + (ax(i)*(dt**2))/2
     posy(i)=posy(i) + vely(i)*dt + (ay(i)*(dt**2))/2
     posz(i)=posz(i) + velz(i)*dt + (az(i)*(dt**2))/2

     if(posx(i)<rad .OR. posx(i)>xlen-rad) velx(i)=-velx(i)
     if(posy(i)<rad .OR. posy(i)>ylen-rad) vely(i)=-vely(i)
     if(posz(i)<rad .OR. posz(i)>zlen-rad) velz(i)=-velz(i)
  enddo
end subroutine pos_upd

!Updating Velocity
subroutine vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)

  implicit none
  integer n,i
  double precision,dimension(n) :: velx,vely,velz,ax,ay,az,tempx,tempy,tempz
  double precision :: dt
  do i=1,n
     velx(i)=velx(i) + 0.5d0*(ax(i)+tempx(i))*dt
     vely(i)=vely(i) + 0.5d0*(ay(i)+tempy(i))*dt
     velz(i)=velz(i) + 0.5d0*(az(i)+tempz(i))*dt
  enddo
end subroutine vel_upd

function pot(posx,posy,posz,eps,n)
  implicit none
  double precision pot,r,dist,eps
  integer i,j,n
  double precision, dimension(n) :: posx,posy,posz
  pot=0.d0
  do i=1,n
     do j=1,n
        if(i==j) cycle
        r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
        pot= pot + (4.d0*eps*((1.d0/r)**12 - (1.d0/r)**6))                 !r is relative distance. x,y,z are components of r.
     enddo
  enddo
end function pot

function kin(velx,vely,velz,m,n)
  implicit none
  double precision :: kin,m
  integer :: i,n
  double precision, dimension(n) :: velx,vely,velz
  kin=0.d0
  do i=1,n
     kin = kin + ((velx(i)**2.d0 + vely(i)**2.d0 + velz(i)**2.d0)/(2.d0*m))
  enddo
end function kin

function dist(x1,y1,z1,x2,y2,z2)
  implicit none
  double precision :: dist,x1,y1,z1,x2,y2,z2
  dist = sqrt((x1-x2)**2.d0 + (y1-y2)**2.d0 + (z1-z2)**2)
end function dist

function Fx(r,x,y,z,eps)
  implicit none
  double precision :: Fx,r,x,y,z,eps
  Fx = 4.d0*eps*((12.d0/r**14.d0) - (6.d0/r**8.d0))*x
end function Fx

function Fy(r,x,y,z,eps)
  implicit none
  double precision :: Fy,r,x,y,z,eps
  Fy = 4.d0*eps*((12.d0/r**14.d0) - (6.d0/r**8.d0))*y
end function Fy

function Fz(r,x,y,z,eps)
  implicit none
  double precision :: Fz,r,x,y,z,eps
  Fz = 4.d0*eps*((12.d0/r**14) - (6.d0/r**8))*z
end function Fz
</code>```


问:
Arpack and Matlab give different values for eigenvalues
I am solving a generalized eigenvalues problem with inversed complex shift:

$$(M-\sigma J)^{-1}J \boldsymbol{x} = \boldsymbol{x} \nu \enspace .$$

My matrices are obtained from a finite element discretization with Getfem++. The jacobian $J$ depends on some parameters. I computed the eigenvalues both with Matlab and directly with Arpack in the same C++ code where I compute the matrices $M$ and $J$.
When I use Matlab to compute the eigenvalues I load the matrices $M$ and $J$ previously saved in Market Matrix format with Getfem++ and I use the function eigs.
Now, in some cases the difference between Arpack eigenvalues and Matlab eigenvalues is really big. In particular, when the values of my parameters are such that I have an eigenvalue with almost zero real part.
With Arpack the real part is of the order of $10^{-10}$ with Matlab is of the order of $10^{-4}$. In addition, in these cases Matlab gives me some warning about the matrix is ill-conditioned.
Have you ever experienced such a problem?

P.S. : I also tried to change default options of eigs but I always obtain the same result!


问:
Rank deficient Jacobian in discretized periodic solutions to autonomous ODE
I'm trying to numerically find periodic solutions to different systems of autonomous nonlinear ordinary differential equations. I decided to use a finite difference scheme and solve the resulting system of equations iteratively, which works reasonably well using a sufficiently good initial guess.

Now most of the systems considered depend on an additional parameter $\lambda$, with a solution in the limit $\lambda\to 0$ easily found, so I would like to use some sort of numerical continuation method to generate a family of solutions for different values of $\lambda$. All literature I found on the subject, e.g. http://dx.doi.org/10.1137/1.9780898719154, assumes the Jacobian to have full rank, which it does not in my case: There is a zero eigenvalue corresponding to a shift in the parametrization of the solution to the original system of ODEs.

The ODEs typically look like
$$
 \ddot{x} = - \sqrt{\dot{x}^2+\dot{y}^2}\ f(x, y, \lambda)\ \dot{y},\\
 \ddot{y} = \sqrt{\dot{x}^2+\dot{y}^2}\ f(x, y,\lambda)\ \dot{x},
$$
with the constraints $x(0) = x(1)$, $y(0) = y(1)$.

Is there either

a good way to somehow fix the parametrization of the solution to remove this degree of freedom or

a numerical continuation method that works with a rank deficient Jacobian?


问:
Large meshing with tetgen
So I have a point cloud that I am creating a 3D flat rectangular surface from.  I'm then turning it into a hollow box and connecting the corners by just dropping this surface mesh down.  I need it to be consistently a meter thick throughout so it is important to drop the surface mesh down so it matches the 1 meter delta z from the top.  So I end up with something that looks like this:

The problem is that tetgen even when I'm allowing it to doesn't want to seem to throw away points.  So I would like a very refined mesh at top and a coarse one at the bottom.  It wants to keep all of the bottom points, so if I drop the mesh one time 1 meter down and give it a high volume constraint I get something like this:

As you can see it is refined at the top but it needs to go back to being refined at the bottom to match the surface mesh.  I don't want tetgen to keep these facets.  I'm wanting it to get coarser as it goes down.  Does anyone know how to fix this problem?  Or maybe some other software that is better suited to my needs?  Thanks for your help and for taking a look.  I'm most proficient in python, so I have used meshpy in the past as well.

Title: Large meshing with tetgen


问:
Conservation at grid interface in adaptive mesh refinement
I am using adaptive mesh refinement to solve one dimensional inviscid Burgers equation. However I am facing some difficulty to handle grid interfaces which are not uniform (coarse-fine grid interface). There has to be conservation of variable at the grid interfaces which I am not able to incorporate in my code.

If anyone has any idea about it or has done the similar work based preferably on finite difference schemes, it could be a great help to me.

Reference: http://epubs.siam.org/doi/abs/10.1137/0724063


问:
Precession of Mercury Python simulation
I was trying to simulate the precession of Mercury based on the perturbed solution:
$$\frac{1}{r}=\frac{m}{B^{2}}(1+e\cos\phi+3\frac{m^{2}}{B^{2}}(1+e\phi \sin\phi +e^{2}(\frac{1}{2}-\frac{1}{6}\cos2\phi))),$$ where $e$ is the eccentricity of mercury, B is a constant which is the angular momentum per mass.
here is my code in Python: http://sprunge.us/TFhU

In the code I have to multiply a large number(SF) so that the precession is slightly visible(or you can make a bigger SF value then the precession would be much easier to see but it would distort the orbits).
I guess it is okay if I use this as a illustration for the mercury precession, however, I want to be able to obtain the value of the precession rate which is known as 43 arcseconds per century, does anyone know how to do that in the code or anyone has python code that does a better job? (I also came through a model based on these two equations: $$r_{n+1} = 2r_{n}-r_{n-1}+\Delta \tau^{2}(-\frac{GM}{{r_{n}}^{2}}+\frac{B^{2}}{{r_{n}}^{3}}-\frac{3GMB^{2}}{{r_{n}}^{4}}),\phi_{n+1}=\phi_{n}+\Delta\tau\frac{B}{[\frac{1}{2}(r_{n+1}+r_{n})]^{2}},$$ where B is same as above, but I couldn't get it to work.)

(Update: I used Vpython and Leapfrog method to do it and the simulation worked out all right. here is the code and modify it as however you want to,
  http://sprunge.us/ehgS (in linux the sphere doesn't leave trail probably because the vpython in linux is very out of date).
  However, I still need to find a way to output the precession rate 43" out of the data.)

Title: Precession of Mercury Python simulation


问:
Eigenvalue problem of the symmetric real operator which corresponds to the symmetric positive definite matrix
I have a real symmetric function $C(x,y)$ defined on $x,y\in[0,\infty)$, i.e. $C(x,y)=C(y,x)$.

I want to solve the eigenvalues problem, i.e. find eigen values and eigen functions:
$$\lambda \psi(x)=\int_0^\infty C(x,y)\psi(y)dy$$

This operator corresponds to the covariance matrix such as $C_{i,j}=Cov[x_i,x_j]$, i.e. it's positive definite, real and symmetric. In the functional case, I'm not sure what is the analog of positive definite. I'm trying to extend the covariance matrix into the infinite dimensional (continuous) case from the matrix version.

I'm looking for references on eigenvalues problem solutions, software packages or papers.


问:
Execute commands when starting ipdb
I usually debug a python script by putting the following line into the source code:

```<code>import ipdb; ipdb.set_trace()
</code>```

Then when I run the script, ipdb starts. Very often I need to plot some arrays in an interactive graph. I use the following commands to make it possible inside the ipdb:

```<code>import matplotlib.pyplot as plt
plt.ion()
</code>```

My question is whether it is possible to run these two commands automatically when ipdb starts.

P. S. I put this question here, in Computational Science, as I believe that ipdb is mostly used by people who do scientific computing.

Title: Execute commands when starting ipdb


问:
MD simulation - Water Vapour
I wanted to know the procedure of starting up a water vapor calculation in LAMMPS. I want to calculate properties of water (i.e. in vapor state) at a higher temperature than 100 Celsius. Currently I am setting up the simulation same as liquid water but with the density of $0.001\, \text{g/cm}^3$.  Initially my box has 1000 times less molecules than that of liquid water and it has 800 atoms in $200\times200\times200\ \mathring{\text{A}}^3$ box.

I am running this simulation in LAMMPS 2014. I have tried all ensembles (NVT, NPT and NVE) with various time intervals. The simulation starts running for only NVE but for a very small time-step of $0.001\, \text{fs}$. If I increase the time interval, the simulation stops with an error "Out of range atoms - cannot compute PPPM". The error also comes if I run a NPT simulation.

I would appreciate if anybody can help me with this problem. I can also attach my input scripts and data file if needed.

Title: MD simulation - Water Vapour


问:
What is the computational cost of using complex numbers in contrast to real numbers in matrix operations, e.g. $LU$ or $LDL^T$ factorizations?
I am curious about how much one loses in terms of computational cost when complex numbers are used instead of real numbers?

I guess the number of floating-point operations and memory doubles, but I am not sure if the wall time doubles necessarily. I am especially asking this for matrix factorizations, $LU$, or $LDL^T$ factorizations of complex symmetric (non-Hermitian) vs real-symmetric matrices. I am not sure if the answer changes based on using sparse or dense storage, but I'd appreciate if you comment on both situations. Obviously symmetric factorization cost should not change, but numerical factorizations should have a performance penalty.

ContentLicense: CC BY-SA 4.0


问:
"Damping factor" for a set of non-linear ODEs
I have a set of four non-linear ODEs representing a negative feedback. I have done parameter variation by random sampling to study the sensitivity of steady state and other dynamic properties to parameter fluctuations.

From the analytical expressions of steady state I have calculated jacobian and its eigenvalues for each parameter set and based on that I could observe the fraction of oscillating cases (complex eigenvalues). I wish to calculate something like damping factor for these cases. I am not sure what to use as a metric. I think of three choices:

Least eigenvalue (absolute value of the real part). This would be rate limiting.

Biggest eigenvalue (absolute value of the real part). This might dictate the initial damping (which is what I am basically interested in).

Determinant of the jacobian (?)

Please suggest me what is the right metric or let me know if there are better alternatives.


问:
How can I efficiently solve $Ax$=$b$ given $A$ is symmetric and contains very small (even negative) eigenvalues using EIGEN
Currently I am using the EIGEN C++ library to try to solve $x$ from the equation $Ax$ = $b$. One problem I encountered is that the matrix $A$ is a correlation matrix with size  > 5000 and can sometimes be rank deficient and have small, even negative eigenvalues. As a result of that, the matrix tends to have huge condition numbers 100~1000+.

Special properties of the matrix $A$ include: symmetric (because it is a correlation matrix) and the correlation decreases as the two entries get further away (e.g. $R_{1,2}$ > $R_{1,6}$).

The current way I used to solve the equation is to use the truncated SVD method to obtain a pseudoinverse of $A$. I also tried to exploit the fact that the matrix is symmetric by using the self-adjoint eigen solver instead of SVD which speed up the process considerably. However, it still remains rather slow and from what I read, people always said one should NEVER EVER try to solve the inverse whenever possible.

As a result of that, I wonder if there is anyway for one to speed up this process even further. I have read that one can use the QR method to obtain the pseudoinverse. My understanding is that for $A$ = $QR$, I can replace $R$ with $R_k$ the when the diagonal at $k$ is the last element bigger than a threshold and any subsequent row will be replaced by 0.

So, my question is: is this QR method correct? Or is there any method that can better solve this kind of equation?


问:
Simple MCMC Algorithm in Matlab
I would be really glad to get some specific advise on how to implement a simple MCMC algorithm (in Matlab, if possible). I'm not yet too familiar with optimization methods. My problem goes as follows:

The function $s(c_1,c_2)$ is a simulation which depends on the two constants $c_1$ and $c_2$, and $f(s(c_1,c_2))$ is the "goodness function" which determines the error of the simulation. This means I want $f$ small in the end. Furthermore, I have a range for both $c_1$ and $c_2$; all values in these ranges are assumed equally probable at the beginning.

Is there anything more one would need to know? My question is now if anyone can provide a simple algorithm or pseudo-algorithm for this specific problem. Thanks in advance.

Title: Simple MCMC Algorithm in Matlab


问:
Need a smart way to numerically take residues in a multidimensional integral
I'm trying to do an integral of the form $\int_C f(u,v) $, where $C$ is a set of contours in $u$ and $v$. In particular, each variable's contour starts at $-\infty+i \epsilon$, goes around a branch cut at $-2g$ (where $g$ is a parameter of the problem), then goes below the real line and out to $-\infty-i \epsilon$ before circling around to $\infty-i \epsilon$, coming back around a branch cut at $2g$ and ending up at $\infty+i \epsilon$. If this description is a bit unclear, there's a drawing of the contour in Figure 5 of this paper.

$f(u,v)$ has poles when $u=v-2i$ and when $v=u-2i$. While the paper suggests integrating above the real axis and replacing the integration below the real axis with analytic residues around these poles, that method isn't feasible for the case I'm working on. Instead, I'm trying to integrate along this contour numerically.

The problem I'm having is that when I try to approximate the (infinite) contour with a finite one, I inevitably end up with a contour that crosses the poles, since the poles in each variable depend on the other one.

I suspect that there's a smarter way to approach all of this, but I can't find much literature on taking residues numerically at all, let alone anything applicable to this particular situation. Any recommendations would be appreciated!

Edit: One thing that I should clarify, I need to have fixed contours in $u$ and $v$, neither can depend on the other. If I had a different $v$ contour for each value of $u$ then taking the $v$ residues analytically would be faster, and I'm trying to get a speedup over doing that.


问:
Stability analysis for explicit time discretization in the Finite Element Method
I have been looking for stability analysis of general reaction-diffusion problems, of the form

$\frac{\partial u}{\partial t}=\nabla\cdot D\nabla u-k\,u$ ,

to be solved using the standard Finite Element in 2D with an explicit Euler time stepping. I haven't found any estimate or proof for the stability criteria for $\Delta t$ for this case, most of the literature seems to point to the solution of the Heat equation by Finite Differences in 1D.

I would like to know if someone can point me how to find the stability criterion in this case, especially for quadrilateral elements.


问:
Implicit time integrator for Chebyshev collocation method for linear hyperbolic system
I want to solve linear hyperbolic system using Chebyshev collocation method. As this method puts severe constraint on the time step for the explicit time integration, I decided to switch to implicit time integration (backward Euler). I use time integrator vode with options method='bdf' and order=1 from Python package scipy.integrate. I do not supply Jacobian to the time integrator, only the vector of right-hand side. In the beginning it was great: the computation with resolution $N=100$ gave me solution 108 times faster with implicit time integration, then with explicit (Runge-Kutta 4th order).

Unfortunately, It turned out that if I increase the grid resolution, then computations become unstable—in initial time period solution is stable, but then it starts to "tremble" and develops instability.

This figure shows the numerically stable solution obtained with resolution $N=200$:

This figure shows the numerically unstable solution obtained with resolution $N=400$:

Could you please recommend to me implicit time integrators that are good to use with Chebyshev collocation method?


问:
Nonlinear 2D modeling of Neural Electromagnetic field in Matlab
I am trying to replicate the MATLAB simulation presented in this paper.  More specifically, I have to code the solution to this equation

$$\frac{a}{2R_i}\frac{\partial^2 V_m}{\partial x^2} - \left[g_\text{Na}m^3 h(V_m - E_\text{Na}) + g_\text{K} n^4(V_m -  E_\text{K}) + g_\text{S}(V_m - E_\text{S})\right] = C_m \frac{\partial V_m}{\partial t} + \frac{a}{2R_i}\frac{\partial E_x}{\partial x}(x,t) \enspace .$$

Here $m$, $n$, and $h$ are gating variables that are ordinary differential equations and are functions of time, and I have solved for them already. I have also written the difference equation for this PDE.  I would like to set $\frac{\partial E}{\partial x}(x,t)$ equal to 1 (denoted as $g$ in the script) for now (and tackle that next).  I am using a finite differences method (I think)

```<code>clear all;
T=1; % Max time (s)
L=1; % Max length of neuron
dt=0.04; %time step
dx=0.01; %space step
t=[0:dt:T]; %time vector
x=[-L:dx:L]; %Space vector

M=length(x);
g=1; %The dE/dx bit
a=5; %Neuron radius
R_i=1; %Intracellular resistance
A=a./2*R_i;
I=0.1; %External Current Applied
ENa=55.17; % mv Na reversal potential
EK=-72.14; % mv K reversal potential
El=-49.42; % mv Leakage reversal potential
g_Na=1.2; % mS/cm^2 Na conductance
g_K=0.36; % mS/cm^2 K conductance
g_l=0.003; % mS/cm^2 Leakage conductance
b=4;
Cm=0.01; %Membrane capacitance
C=Cm;
%Solve for m,n,h

%Solve for m, n, h
V=-60; % Initial Membrane voltage
m=am(V)/(am(V)+bm(V)); % Initial m-value
n=an(V)/(an(V)+bn(V)); % Initial n-value
h=ah(V)/(ah(V)+bh(V)); % Initial h-value
y0=[V;n;m;h];
tspan = [0,max(t)];
%Matlab's ode45 function
[time,V] = ode45(@HH,tspan,y0);
V_hh=V(:,1);
n=V(:,2);
m=V(:,3);
h=V(:,4);
N=length(time);
plot(time,m,time,h,time,n);
V_initialgrid=zeros(M,N);
figure();
plot(time,V_hh)
V_initialgrid(1,:)=V_hh';

xlabel('time');
ylabel('length');
for j=2:N-1 %Stepping through time
    B=g_Na*m(j).^3*h(j)+g_K*n(j).^4+g_l;
    for i=2; %stepping thru space
        V_initialgrid(i,j)=(A* V_initialgrid(i-1,j)+B* V_initialgrid(i,j)+C* V_initialgrid(i-1,j)+g);

        for i=3:M-1 %stepping thru space
            V_initialgrid(i+1,j)=(-dt*(g-(A.*(( V_initialgrid(i,j-1)-2.* V_initialgrid(i,j)+ V_initialgrid(i,j+1))./(dx.^2))-B.*-60))./C)+ V_initialgrid(i,j);
        end
    end
end
    imagesc(log(abs(V_initialgrid)))

    xlabel('Time');
    ylabel('L');
</code>```

The boundary conditions are as follows:

$$\frac{\partial V}{\partial x}=\frac{\partial m}{\partial x}=\frac{\partial h}{\partial x}=\frac{\partial n}{\partial x}=0\quad \text{for } x=\pm L$$

and

$$\frac{\partial V}{\partial t}=\frac{\partial m}{\partial t}=\frac{\partial h}{\partial t}=\frac{\partial n}{\partial t}=0 \enspace .$$

I know my error lies within the indexing of the solution of the difference equation.


问:
Sparse Linear Algebra vs Dense Linear Algebra
I am interested in a reference in the literature that discusses the performance of Dense Linear Algebra (blas routines) and dense linear algebra (sparse blas routines).

I am interested in knowing for what combinations of size and density sparse routines outperform dense routines.  I am mainly interested in a shared memory computer.  Does anyone know of a reference?

Thanks,


问:
Relaxation - spreadsheet solution to recursive algorithm
I am setting up a simple model that uses recursion to iterate between known constraints - relaxation, in other words. I need a spreadsheet that allows:

```<code>A0 50

A1 =(A0+A2)/2

A2 =(A1+A3)/2

A3 =(A2+A4)/2

A4 100
</code>```

It should be simple, but OpenOffice and Apple Numbers specifically forbid this. If you know of configuration setting to override this, that would be great too.


问:
Assigning new values based on original guesstimates and ranking / ordering?
Lets say we have two things as input, $N$ scalars (measurements) that we know are erroneous to some degree (i.e. the correct values are somewhat similar).

In addition, we also have a roughly more reliable order $O$ among the measurements. Using the ordering:

$O = <u_1, u_2, u_3 ... u_n>$ we can deduct that with some known probability, $p > 0.5$ the ordering is correct, $[i \geq j \rightarrow u_i \geq u_j]$ is true. Now, how would you then adjust the scalar values accordingly?

Currently:

Creating a new set of values, cloning the original scalars.

Making the minimal needed changes in order to satisfy the order. i.e. slightly increasing or decreasing the original values until we get a corrctly ordered $\{v_1, v_2, v_3 ... v_n\}$ such that $\sum\limits_{i=1}^n (u_i - v_i)^2$ is minimal.

Doing a weighted average between the newly generated values and the original ones, using $p$ as the weight for the new value.

While $p$ is unknown, we can guesstimate it from data that is manually inspected.


问:
matplotlib contourplot for $\log z$ in the Complex Plane $\mathbb{C}$
I tried using Python's matplotlib on the logarithm and here is what I got, a kind of starburst pattern.  Since the angle jumps between $\theta = 0$ and $\theta = 2\pi$, contour assumes there is a contour line between those those two points.  Here in in fact $\bbox[#EEEDAA,1pt]{0=2\pi}$

Perhaps instead of matplotlib.pyplot.contour the contour lines can be connected with meshgrid and numpy.where statements.

Contour Plot of $\log z$ on $z \in [-1,1] + i[-1,1]$

```<code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize']=4,4

plt.contour(
    np.arange(-1,1,0.01),
    np.arange(-1,1,0.01),
    np.angle(t[...,None] + 1j*t[None,...]) ,
    levels = 2*np.pi*np.arange(-1,1,0.025)
)
</code>```

Contour Lines Accumulate at the Negative Imaginary Axis $-i \, \mathbb{R}_{\geq 0}$

```<code>t = np.arange(-1,1,0.025)
plt.contour(t,t, np.angle(t[...,None] + 1j*t[None,...]) , levels = 2*np.pi*np.arange(-1,1,0.025))

L = 0.25
plt.xlim([-L,L])
plt.ylim([-L,L])
</code>```

Code Inspection

```<code>import inspect
print inspect.getsource(plt.contour)

@_autogen_docstring(Axes.contour)
def contour(*args, **kwargs):
    ax = gca()
    # allow callers to override the hold state by passing hold=True|False
    washold = ax.ishold()
    hold = kwargs.pop('hold', None)
    if hold is not None:
        ax.hold(hold)
    try:
        ret = ax.contour(*args, **kwargs)
        draw_if_interactive()
    finally:
        ax.hold(washold)
    if ret._A is not None: sci(ret)
    return ret
</code>```

```<code>print inspect.getsource(plt.Axes.contour)

def contour(self, *args, **kwargs):
    if not self._hold:
        self.cla()
    kwargs['filled'] = False
    return mcontour.QuadContourSet(self, *args, **kwargs)
</code>```

Finally one can call on the function containing the actual source code:

```<code>import matplotlib
print inspect.getsource(matplotlib.contour.QuadContourSet)
</code>```


问:
Solving a system of 4 coupled PDEs representing variable diffusivity
I have four partial differential equations representing mass conservation of two compressible fluid phases (marked by subscripts $p1$ and $p2$) in two different continuum media (marked by subscripts $c1$ and $c2$) and they are coupled together by the term $T_{c1c2}$ as shown below.

Continuum 1:
$$\begin{align}
\text{Phase 1: } \frac{\partial }{\partial t}\left(\frac{\phi_{c1}S_{p1c1}}{B_{p1c1}} \right) = \nabla\left(\lambda_{p1c1}\nabla{P_{p1c1}} \right) - \frac{T_{c1c2}}{V} - \frac{q_{p1c1}}{V}, \\
\text{Phase-2:} \frac{\partial }{\partial t}\left(\frac{\phi_{c1}S_{p2c1}}{B_{p2c1}} \right) = \nabla\left(\lambda_{p2c1}\nabla{P_{p2c1}} \right) - \frac{T_{c1c2}}{V} - \frac{q_{p2c1}}{V},
\end{align}
$$
Continuum 2:
$$\begin{align}
\text{Phase 1: }\frac{\partial }{\partial t}\left(\frac{\phi_{c2}S_{p1c2}}{B_{p1c2}} \right) = \nabla\left(\lambda_{p1c2}\nabla{P_{p1c2}} \right) - \frac{T_{c1c2}}{V} - \frac{q_{p1c2}}{V}, \\
\text{Phase 2: } \frac{\partial }{\partial t}\left(\frac{\phi_{c2}S_{p2c2}}{B_{p2c2}} \right) = \nabla\left(\lambda_{p2c2}\nabla{P_{p2c2}} \right) - \frac{T_{c1c2}}{V} - \frac{q_{p2c2}}{V},
\end{align}$$

Additionally, there are 4 more linear equations which are:
$$S_{p1c1} + S_{p2c1} = 1$$
$$S_{p1c2} + S_{p2c2} = 1$$
$$P_{p2c1} - P_{p1c1} = P_{cap1}$$
$$P_{p2c2} - P_{p1c2} = P_{cap2}$$

Variables: In these 4 linear equations, $P_{cap1}$ and $P_{cap2}$ are known.
Except for $\phi_{c1}$, $\phi_{c2}$ and $V$, which are constants, all other variables vary with time $t$ and distance $(x,y)$. Also,
$$\frac{1}{B_{p1c1}} = 1+c_{p1}\left(P_{p1c1} - P^{STC}\right) \text{and} \frac{1}{B_{p2c1}} = 1+c_{p2}\left(P_{p2c1} - P^{STC}\right)$$ where $c_{p1}$, $c_{p2}$ and $P^{STC}$ are constants.

Initial Condition: Dirichlet constant pressure

Boundary Condition: Neumann no-flow boundary condition

Objective: To solve for saturations and pressures as a function of space and time, i.e.

$S_{p1c1}$, $S_{p2c1}$, $S_{p1c2}$, $S_{p2c2}$ and
$P_{p1c1}$, $P_{p2c1}$, $P_{p1c2}$, $P_{p2c2}$,

This system of equations could be discretized into linear system of equations and then solved for each time step, however, that is a complex task and would take a lot of time. I would, however, like to use some form of "tool" which can solve this coupled system of equations in a mesh grid format. I saw an answer in other post suggested to use FiPy, but the feedback from the OP was that it's way too slow. Please suggest what could be a good tool to solve this problem.


问:
Method of lines: explanation for high precision low resolution flux and vice versa
Could someone please explain what is meant with e.g high precision low resolution flux. Does it mean, that I use some higher order method (still first order) but with a more complex stencil and consequently lower resolution ? While low precision and high resolution means essentially a forward or backward difference formula ? As a consequence, it is on me to choose the appropriate combinations of difference methods and flux limiter (e.g by solving an analytical form of my PDE prior to the more complex one ?)</p>


问:
1D fit on a surface in matlab
I have a problem fitting a $N\times N$ ($32\times 32$ in this case) matrix in Matlab.

Here is the deal

I have this image, and I want to create 1 dimensional fit between the points and then be able to interpolate it in order to make it smoother.

Does anyone have any idea if that can be done?

Title: 1D fit on a surface in matlab


问:
Return derivative vector from odeint scipy function
I have the following function that I want to implement in scipy.integrate.odeint

```<code>def numericSolver(rv_,t,alpha=1.,gamma=1.):
    BVal = B_f(alpha,gamma,t) #this function returns a vector
    BdotVal = Bdot_f(alpha,gamma,t) #this function returns a vector
    matDynamics = np.array([[0.,1.,float(BdotVal[2]),0.,0.,0.],
            [0.,0.,0.,float(BVal[2]),0,-float(BVal[1])],
            [0.,0.,0.,1,0.,0.],
            [0.,-float(BVal[2]),0.,0.,0.,float(BVal[0])],
            [0.,0.,0.,0.,0.,1.],
            [0.,float(BVal[1]),0.,-float(BVal[0]),0.,0.]])
    derivatives = matDynamics.dot(rv_)
    return derivatives
</code>```

In order to implement this function I call odeint as

```<code>solution = odeint(numericSolver,IC,t=t,args=(alpha,gamma))
</code>```

where IC is the initial condition vector, t is the time vector of interest, and alpha and gamma are inputs for numericSolver.

For a more clear explanation, derivatives is the solution, $\dot{\vec{X}}$, to the equation

$
\dot{\vec{X}} = M \vec{X},
$

where $M$ is matDynamics and rv_ is $\vec{X}$.

Is there a way that I can get odeint to return the values of derivatives at each value in t? Or is it best to calculate the derivatives values manually using the values returned in solution?


问:
Computing entries of a matrix in parallel using GPU
I am trying to solve $Ax=b$ a system of linear equation. The matrix $A=(a_{ij})$ is constructed as $$a_{ij}=\int_0^1\int_0^1 K(x,y)h_i(x)h_j(y) dxdy.$$ The evaluation of integral takes a lot of time, and overall computation is very slow. Can I evaluate these integrals in parallel using GPU cores? Please suggest how can I compute these integrals simultaneously??

thanks


问:
C++: How to find the roots of polynomial modulus N
I need to know what library, given a vector of coefficients and modulus N return the roots of polynomial modulus N. The library should support big integer. I'm coding in C++. </p>


问:
Monotonic convergence of Newton's method for boundary value problems
I’m interested in solving nonlinear elliptic boundary value problems of the type
$$
-a\Delta u + f(u) = 0,
$$
$$
u|_\Gamma = u_0
$$
by Newton’s method when its convergence is global and monotonic.
Could you advice some references concerning this problem, containing proofs of global convergence?

Newton's method takes the form
$$
-a\Delta u + f(\widetilde u) + f'(\widetilde u)(u - \widetilde u) = 0
$$
where $\widetilde u$ is the previous approximation for the solution.


问:
ANSYS Fluent: Defining External Force on a Cell Zone
I am using ANSYS Fluent 15 to simulate a case of fluid flow. In my case, I have three cell zones in my simulation, let's call them top, middle and bottom.

My requirement is to put a vertical constant force in the "middle" cell zone. I can not use gravity because it is then acted in all other cell zones as well. My purpose is to define this force exclusively in the middle cell zone only.
Apart from it, I tried to use the "Y-momentum" source term in cell zone conditions panel but it a did not seem to work either.

Any suggestion or help would be great. Thanks


问:
How can i get this number set to be recooperated from a product sum?
W==a(x)=2*(W* x)/41;

I'm attempting to make a new type of compression algorithm. A rubber band ball type.

I'm trying to make a software that extracts, from that formula, the secret volumes of number sets within. By taking the the opposite formula, b(x) i should get the answer. But unfortunately, i dont. I have to search 0-255(+1). Through the first generation the numbers would be like as follows
(W starts at 1)
    a(82*1)==2
    b(2/1)==82

I need the numbers reappear for me. How can I do this without cheating?

Q: What is the best if statement to get my numbers to reappear?

My code as is right now.


问:
Graph optimization for parallel processing
Consider the following example structure of overlapping images marked A,B,C,D. The possible overlaps are marked by gray color:

The structure can be represented by a weighted undirected graph (images being the nodes, edges representing overlaps and the weights representing overlap areas or time needed to process them):

The images are to be rendered on surface while overlaps require extra processing (e.g. avergaing the pixels for transparency).

For example, we can decide to merge A,B first. The situation becomes:

After each merge, the graph has to be updated.
In this case, only images C,D are to be merged and the algorithm finishes with single composite image.

You can observe that if the merging algorithm runs in parallel, it can now process C and D simultaneously.

However, if we start merging A,C first, then B and D are blocked and the algorithm can run only sequentially (considering that individual image can be used by single thread / cannot be shared).

Given the graph, I would like to schedule parallel computation for $N$ computing units (e.g. CPU cores).

The schedule for $N=2$ can look like this

```<code>iteration #1
    core 1: A-B
    core 2: waiting
iteration #2
    core 1: (A-B)C
    core 2: (A-B)D
</code>```

I think the scheduler should aim for two objectives:

Minimize the number of iterations (i.e. utilize maximum number of cores at any given time)

Balance load to all cores at any given iteration (i.e. assign them approximately same areas of overlap so that they will finish approximately at the same time; unnecessary waiting would spoil the benefits of parallel computation).

The first objective is important while the second is "nice to have" (the time needed to process an overlap is proportional to its area).

I am stuck with two major problems:

What algorithm to use? Backtracking with storing current best solution?

Would it be easier to schedule tasks for $N=\infty$ (maximize parallelism) rather that $N$ being corresponding to actual number of processing units available?

Only the overlap areas between pairs of images are known, but since the image are merged in the process, how to compute new overlaps (e.g. C(A-B))?

The greedy version of the algorithm would look like this:

```<code>Grab as much overlapping image pairs as possible such that no two pairs themselves are overlapping.
Sort the pairs in N bins such that each bin contains roughly same total area.
Merge images in parallel.
Update graph and repeat until single node remains.
</code>```

I know this question is maybe too broad. I believe this problem has been certainly solved many times in parallel computation projects so many people from parallel computing areas should have some immediate ideas and hints.

Once solved, I plan to write an elaborate answer to benefit you as well. But first I need something to start with. Thanks.


问:
Petrov-Galerkin enrichment method for Darcy equation
I was reading about Petrov-Galerkin Enrichment Method for Darcy equations. Here are a couple papers that discuss this in detail:

A Petrov–Galerkin enriched method: A mass conservative finite element method for the Darcy equation.

A Symmetric Nodal Conservative Finite Element Method for the Darcy Equation.

However, I am not sure I understand what they are doing completely and would like some help understanding how they are constructing these "multiscale" or enrichment functions.

From what I understand, it seems they are "enriching" the $P_1$ velocity by an $RT_0$, but what exactly does this mean? Specifically, what is meant when they "correct edge residuals"? Because I would like to implement this in something like FEniCS (if possible) but in order to do that I need to understand this method completely.


问:
Efficient principal pivots
It was suggested I should try posting this question here from Mathematics

Background

I'm working on a numerical linear algebra package in C#.  I'm trying to implement a variety of "principal pivoting" methods to solve optimization problems (specifically linear complementarity problems, but I think similar schemes are used in linear and quadratic programming).  These methods require the ability to multiply a column vector against the "principal pivotal transform" of a matrix.

The principal pivot transform (PPT) of a matrix $M$ is defined as:

$$PM^'P^T = \left (\begin{array}{cc} M_{\alpha\alpha}^{-1} &\ -M_{\alpha\alpha}^{-1} M_{\alpha\beta} \\ M_{\beta\alpha}M_{\alpha\alpha}^{-1} &\ M_{\beta\beta} - M_{\beta\alpha} M_{\alpha\alpha}^{-1}M_{\alpha\beta} \end{array} \right )$$

Where $\alpha$ is the set of pivoted row/column indices, and $\beta$ is the set of unpivoted row/column indices.  $P$ is a permutation matrix that keeps the pivoted and unpivoted rows/columns together to make the block matrix form above cleaner.  $M_{\alpha\beta}$ is a submatrix of $M$ with the rows in $\alpha$ and the columns in $\beta$, and the other submatrices of $M$ are defined similarly.  Also, $M_{\alpha\alpha}^{-1} = (M_{\alpha\alpha})^{-1}$ and not $(M^{-1})_{\alpha\alpha}$

...

Each iteration of a principal pivot method swaps a single index from $\alpha$ to $\beta$ or vice-verse.  This single pivoting operation can be expressed as a low rank update of the matrix.  Say we're pivoting index $r$ and using $M^'$ to represent the updated matrix:

$$\begin{matrix}
M_{rr}^' & = & & \frac{1}{M_{rr}} & \\ \\ \\
M_{ir}^' & = & & \frac{M_{ir}}{M_{rr}}, & (i \neq r) \\ \\ \\
M_{ri}^' & = & & \frac{-M_{ri}}{M_{rr}}, & (i \neq r) \\ \\ \\
M_{ij}^' & = & M_{ij} - & \frac{M_{ir} * M_{ri}}{M_{rr}}, & (i \neq r, j \neq r)
\end{matrix}
$$

This is what I'm currently doing.  Unfortunately, this is basically calculating $M_{\alpha\alpha}^{-1}$ using something like Gaussian elimination.  I have a few test cases where the results are suffering from numerical problems even on on small problems, even using double precision floating point.  To help combat this, if I remove an index from $\alpha$ (that is, if I unpivot an index), I recalculate the matrix from scratch using one pivot at a time with the above scheme.  But this makes the Principal Pivot algorithm I have quite slow.

Question

I know there are reasonable ways to calculate the product $M_{\alpha\alpha}^{-1} z$ for some column vector $z$ without calculating $M_{\alpha\alpha}^{-1}$ directly (such as LU decomposition or QR decomposition).  I could calculate this decomposition fresh every iteration from the current pivot set ($\alpha$), and probably get a pretty robust method.  But then each iteration of the pivoting algorithm would take $O(n^3)$; that's too slow.  I could do a low rank update of the LU or QR decomposition, but since most of the algorithms pivot a single index at a time, I lose most of the benefit of these decomposition schemes after a few iterations.  I could perform the full decomposition only every few iterations, but then the question becomes how often should I do so?  And it still would be quite slow.

Is there something more clever that can be done?  I think I basically need to find the decomposition of principal submatrices.  Is there any literature on the subject?  Any other ways to approach the problem?

Title: Efficient principal pivots


问:
C60 orbital calculation
I am currently trying to reproduce the results published by Hands et al.1 using MATLAB. They calculated the bases of the C60 wave functions of HOMO, LUMO and LUMO+1.

I did the following:

I implemented these functions into MATLAB syntax, e.g.

Tp1uxL5 = @(x,y,z) 35.*z.*(x.^4-6*x.^2.*y.^2+y.^4)-5*x.*(1-14*z.^2+21.*z.^4);
Tp1uyL5 = @(x,y,z) 5*y.*(28.*z.^4-28*x.*z.*(x.^2-y.^2)-(1-7*z.^2).^2);
Tp1uzL5 = @(x,y,z) 2*z.*(15-70*z.^2+63*z.^4)+7*x.*(5*(x.^2-y.^2).^2-4*x.^2);

for the pentagon prone $T^p_{1ui} (L=5)$ (i=x,y,z) functions, multiplied by an exponential function $\exp(-k R)$ utilizing $k = \frac{Z_\text{eff}}{2} a_0$ with the effective charge $Z_\text{eff}=3.14$, the Bohr radius $a_0 = 52 \text{ pm}$, and $R = \sqrt{x^2+y^2+z^2}$ as proposed by the authors.

I normalized the functions by integrating in three dimensions using MATLAB's integral3 function, i.e.

A = integral3( FKT , -B,B , -B,B , -B,B )

with the boundary $B$ which i varied between $0.5 \text{ nm} \le B \le 10 \text{ nm}$ with no change in result. I also tried just generating the 3D "wavefunction matrix" in those boundaries and taking the sum over all values.

I equated each base $T$ times exponential times normalization factor $\frac{1}{\sqrt{A}}$, squared for the propability amplitude, added the corresponding functions for the probability density $\Psi^2$ of HOMO, LUMO and LUMO+1 as well as pentagon, hexagon and double-bond prone configuration.

I generated plots:

A) "Constant height" STM images at $z=z_0$ by equating $\Psi^2 (x,y,z_0)$

B) "Constant current" STM images equating $\Psi^2 (x,y,z) = const.$ and visualizing the $z$ map.

C) 3D plots of the electron density (comparable to Fig. 4,6,8,11 in 1) by MATLABs isosurface(X,Y,Z,DOS,DOSvalue).

Unfortunately, I fail to reproduce the pictures of the report. Below you see a calculated constant current image; constant height and the 3D probability density fronts show the same structure, so the visualization part of the code I guess is fine.

I checked and checked the functions again, I am starting to get mad over this. I am not reproducing the results of the paper. What am I doing wrong?

Fig. 1. Simulated constant current images. First row: LUMO, second row: HOMO, third row: LUMO+1. First column: pentagone prone, second column: hexagone prone, third column double-bond prone.

View image here.

1 "Calculation of images of oriented C60 molecules using molecular orbital theory" Phys. Rev. B 81, 205440 (2010), http://core.ac.uk/download/pdf/97962.pdf

Title: C60 orbital calculation


问:
Pastix installation
I'm trying to install Pastix library, but I'm still having problems with missing libraries and I'm really helpless now.

I started with copying LINUX-GNU.in as config.in, because I'm working on Lubuntu 15.04. I've changed paths to Scotch and Blas.

But now, when I try to make examples_s, it appears ok, but no compiled files appear in the bin folder.

When I try to make any other example, e.g. make examples_z, eventually I'm getting the error:

```<code>/usr/bin/ld: cannot find -lblas
</code>```

I know that Blas library couldn't be found, but I don't see the mistake I made. Blas libraries are located in my /usr/lib folder like this:

```<code>├── libblas
|     ├── libblas.so.3 -> libblas.so.3.0
|     └── libblas.so.3.0
├── libblas.so.3 -> /etc/alternatives/libblas.so.3
├── libblas.so.3gf -> /etc/alternatives/libblas.so.3gf
├── libgslcblas.so.0 -> libgslcblas.so.0.0.0
└── libgslcblas.so.0.0.0
</code>```

In config.in I've set BLAS_HOME variable like this:

```<code>BLAS_HOME=/usr/lib/libblas
</code>```

So, do you know, what I'm doing wrong? I have no experience with Pastix and so I'm very confused right now.

Solution

I've finally found out, that Pastix needs specifically libblas.a library.

So, the solution is:

Find the package which contains libblas.a - apt-file search libblas.a

Choose one of the packages - e.g. libblas-dev

Install the chosen package - sudo apt-get install libblas-dev

config.in

```<code>HOSTARCH    = i686_pc_linux
VERSIONBIT  = _32bit
EXEEXT      =
OBJEXT      = .o
LIBEXT      = .a
CCPROG      = gcc -Wall
CFPROG      = gfortran
CF90PROG    = gfortran -ffree-form
MCFPROG     = mpif90
CF90CCPOPT  = -ffree-form -x f95-cpp-input
# Compilation options for optimization (make expor)
CCFOPT      = -O3
# Compilation options for debug (make | make debug)
CCFDEB      = -g3
CXXOPT      = -O3
NVCCOPT     = -O3

LKFOPT      =
MKPROG      = make
MPCCPROG    = mpicc -Wall
MPCXXPROG   = mpic++ -Wall
CPP         = cpp
ARFLAGS     = ruv
ARPROG      = ar
EXTRALIB    = -lgfortran -lm -lrt
CTAGSPROG   = ctags

VERSIONMPI  = _mpi
VERSIONSMP  = _smp
VERSIONSCH  = _static
VERSIONINT  = _int
VERSIONPRC  = _simple
VERSIONFLT  = _real
VERSIONORD  = _scotch

###################################################################
#                  SETTING INSTALL DIRECTORIES                    #
###################################################################
# ROOT          = /home/osboxes/pastix
# INCLUDEDIR    = ${ROOT}/include
    # LIBDIR        = ${ROOT}/lib
# BINDIR        = ${ROOT}/bin
    # PYTHON_PREFIX = ${ROOT}

###################################################################
#                  SHARED LIBRARY GENERATION                      #
###################################################################
#SHARED=1
#SOEXT=.so
#SHARED_FLAGS =  -shared -Wl,-soname,__SO_NAME__
#CCFDEB       := ${CCFDEB} -fPIC
    #CCFOPT       := ${CCFOPT} -fPIC
#CFPROG       := ${CFPROG} -fPIC

###################################################################
#                          INTEGER TYPE                           #
###################################################################
# Uncomment the following lines for integer type support (Only 1)

#VERSIONINT  = _long
#CCTYPES     = -DFORCE_LONG -DINTSIZELONG
#---------------------------
#VERSIONINT  = _int32
#CCTYPES     = -DINTSIZE32
#---------------------------
#VERSIONINT  = _int64
#CCTYPES     = -DINTSSIZE64

###################################################################
#                           FLOAT TYPE                            #
###################################################################
CCTYPESFLT  =
# Uncomment the following lines for double precision support
VERSIONPRC  = _double
CCTYPESFLT := $(CCTYPESFLT) -DPREC_DOUBLE

# Uncomment the following lines for float=complex support
#VERSIONFLT  = _complex
#CCTYPESFLT := $(CCTYPESFLT) -DTYPE_COMPLEX

###################################################################
#                          MPI/THREADS                            #
###################################################################

# Uncomment the following lines for sequential (NOMPI) version
#VERSIONMPI  = _nompi
#CCTYPES    := $(CCTYPES) -DFORCE_NOMPI
    #MPCCPROG    = $(CCPROG)
#MCFPROG     = $(CFPROG)

# Uncomment the following lines for non-threaded (NOSMP) version
#VERSIONSMP  = _nosmp
#CCTYPES    := $(CCTYPES) -DFORCE_NOSMP

# Uncomment the following line to enable a progression thread,
#  then use IPARM_THREAD_COMM_MODE
#CCPASTIX   := $(CCPASTIX) -DPASTIX_THREAD_COMM

# Uncomment the following line if your MPI doesn't support MPI_THREAD_MULTIPLE level,
#  then use IPARM_THREAD_COMM_MODE
#CCPASTIX   := $(CCPASTIX) -DPASTIX_FUNNELED

# Uncomment the following line if your MPI doesn't support MPI_Datatype correctly
#CCPASTIX   := $(CCPASTIX) -DNO_MPI_TYPE

# Uncomment the following line if you want to use semaphore barrier
# instead of MPI barrier (with IPARM_AUTOSPLIT_COMM)
#CCPASTIX    := $(CCPASTIX) -DWITH_SEM_BARRIER

# Uncomment the following lines to enable StarPU.
#CCPASTIX   := $(CCPASTIX) `pkg-config libstarpu --cflags` -DWITH_STARPU
    #EXTRALIB   := $(EXTRALIB) `pkg-config libstarpu --libs`
# Uncomment the correct 2 lines
#CCPASTIX   := $(CCPASTIX) -DCUDA_SM_VERSION=11
    #NVCCOPT    := $(NVCCOPT) -maxrregcount 32 -arch sm_11
#CCPASTIX   := $(CCPASTIX) -DCUDA_SM_VERSION=13
    #NVCCOPT    := $(NVCCOPT) -maxrregcount 32 -arch sm_13
CCPASTIX   := $(CCPASTIX) -DCUDA_SM_VERSION=20
    NVCCOPT    := $(NVCCOPT) -arch sm_20

# Uncomment the following line to enable StarPU profiling
# ( IPARM_VERBOSE > API_VERBOSE_NO ).
#CCPASTIX   := $(CCPASTIX) -DSTARPU_PROFILING

# Uncomment the following line to disable CUDA (StarPU)
#CCPASTIX   := $(CCPASTIX) -DFORCE_NO_CUDA

###################################################################
#                          Options                                #
###################################################################

# Show memory usage statistics
#CCPASTIX   := $(CCPASTIX) -DMEMORY_USAGE

# Show memory usage statistics in solver
#CCPASTIX   := $(CCPASTIX) -DSTATS_SOPALIN

# Uncomment following line for dynamic thread scheduling support
#CCPASTIX   := $(CCPASTIX) -DPASTIX_DYNSCHED

# Uncomment the following lines for Out-of-core
#CCPASTIX   := $(CCPASTIX) -DOOC -DOOC_NOCOEFINIT -DOOC_DETECT_DEADLOCKS

###################################################################
#                      GRAPH PARTITIONING                         #
###################################################################

# Uncomment the following lines for using metis ordering
#VERSIONORD  = _metis
#METIS_HOME  = ${HOME}/metis-4.0
    #CCPASTIX   := $(CCPASTIX) -DMETIS -I$(METIS_HOME)/Lib
    #EXTRALIB   := $(EXTRALIB) -L$(METIS_HOME) -lmetis

# Scotch always needed to compile
#SCOTCH_HOME ?= ${HOME}/scotch_5.1/
    SCOTCH_INC ?= /usr/include/scotch-int64  #$(SCOTCH_HOME)/include
SCOTCH_LIB ?= /usr/lib  #$(SCOTCH_HOME)/lib
    # Uncomment on of this blocks
    #scotch
    #CCPASTIX   := $(CCPASTIX) -I$(SCOTCH_INC) -DWITH_SCOTCH
    #EXTRALIB   := $(EXTRALIB) -L$(SCOTCH_LIB) -lscotch -lscotcherrexit
    #ptscotch
    CCPASTIX   := $(CCPASTIX) -I$(SCOTCH_INC) -DDISTRIBUTED -DWITH_SCOTCH
    #if scotch >= 6.0
    EXTRALIB   := $(EXTRALIB) -L$(SCOTCH_LIB) -lptscotch -lscotch -lptscotcherrexit
    #else
    #EXTRALIB   := $(EXTRALIB) -L$(SCOTCH_LIB) -lptscotch -lptscotcherrexit

###################################################################
#                Portable Hardware Locality                       #
###################################################################
# By default PaStiX uses hwloc to bind threads,
# comment this lines if you don't want it (not recommended)
HWLOC_HOME ?= /opt/hwloc/
HWLOC_INC  ?= $(HWLOC_HOME)/include
    HWLOC_LIB  ?= $(HWLOC_HOME)/lib
CCPASTIX   := $(CCPASTIX) -I$(HWLOC_INC) -DWITH_HWLOC
EXTRALIB   := $(EXTRALIB) -L$(HWLOC_LIB) -lhwloc

###################################################################
#                             MARCEL                              #
###################################################################

# Uncomment following lines for marcel thread support
#VERSIONSMP := $(VERSIONSMP)_marcel
    #CCPASTIX   := $(CCPASTIX) `pm2-config --cflags` -I${PM2_ROOT}/marcel/include/pthread
    #EXTRALIB   := $(EXTRALIB) `pm2-config --libs`
# ---- Thread Posix ------
EXTRALIB   := $(EXTRALIB) -lpthread

# Uncomment following line for bubblesched framework support (need marcel support)
#VERSIONSCH  = _dyn
#CCPASTIX   := $(CCPASTIX) -DPASTIX_BUBBLESCHED

###################################################################
#                              BLAS                               #
###################################################################

# Choose Blas library (Only 1)
# Do not forget to set BLAS_HOME if it is not in your environnement
BLAS_HOME=/usr/lib/libblas
#----  Blas    ----
BLASLIB  = -lblas
#---- Gotoblas ----
#BLASLIB  = -L${BLAS_HOME} -lgoto
    #----  MKL     ----
    #Uncomment the correct line
    #BLASLIB  = -L$(BLAS_HOME) -lmkl_intel_lp64 -lmkl_sequential -lmkl_core
#BLASLIB  = -L$(BLAS_HOME) -lmkl_intel -lmkl_sequential -lmkl_core
    #----  Acml    ----
    #BLASLIB  = -L$(BLAS_HOME) -lacml

###################################################################
#                             MURGE                               #
###################################################################
# Uncomment if you need MURGE interface to be thread safe
# CCPASTIX   := $(CCPASTIX) -DMURGE_THREADSAFE
    # Uncomment this to have more timings inside MURGE
    # CCPASTIX   := $(CCPASTIX) -DMURGE_TIME

###################################################################
#                          DO NOT TOUCH                           #
###################################################################

FOPT      := $(CCFOPT)
    FDEB      := $(CCFDEB)
CCHEAD    := $(CCPROG) $(CCTYPES) $(CCFOPT)
    CCFOPT    := $(CCFOPT) $(CCTYPES) $(CCPASTIX)
CCFDEB    := $(CCFDEB) $(CCTYPES) $(CCPASTIX)
    NVCCOPT   := $(NVCCOPT) $(CCTYPES) $(CCPASTIX)

###################################################################
#                        MURGE COMPATIBILITY                      #
###################################################################

MAKE     = $(MKPROG)
    CC       = $(MPCCPROG)
CFLAGS   = $(CCFOPT) $(CCTYPESFLT)
FC       = $(MCFPROG)
    FFLAGS   = $(CCFOPT)
LDFLAGS  = $(EXTRALIB) $(BLASLIB)
CTAGS    = $(CTAGSPROG)
</code>```

Title: Pastix installation


问:
1D k-epsilon turbulence model in a turbidity current
I am trying to implement a 1D k-$\varepsilon$ turbulent model for a turbidity current, hence the conservation equation for $c$. I'm solving for the variables $u,c,k$ and $\varepsilon$. The remaining terms are parameters.
$$
\newcommand{\cmu}{c_{\mu}}
\newcommand{\cmuc}{c_{\mu c}}
\newcommand{\dervz}[1]{\frac{\partial #1}{\partial z}}
\newcommand{\eddyvis}{\cmu \frac{k^2}{\varepsilon}}
\newcommand{\eddyvisC}{\cmuc \frac{\hat{k}^2}{\hat{\varepsilon}}}
$$

$$
\frac{\partial u }{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \cmu \frac{k^2}{\varepsilon} \right) \frac{\partial u }{\partial z} \right]
+ \frac{\rho_s - \rho_w}{\rho_w} c g_x
$$

$$
 \frac{\partial c}{\partial t} =
 v_s \dervz{c} + \frac{\partial}{\partial z} \left[ \left( \nu + \cmuc \frac{k^2}{\varepsilon} \right) \dervz{c} \right]
$$

$$
 \frac{\partial k}{\partial t} =
  \frac{\partial}{\partial z}  \left[ \left( \nu + \frac{\cmu}{\sigma_k} \frac{k^2}{\varepsilon} \right) \dervz{k} \right] + \eddyvis \left(\dervz{U}\right)^2 - \varepsilon
  + \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} \eddyvisC \dervz{c}
$$

$$
 \frac{\partial \varepsilon}{\partial t} =
  \frac{\partial}{\partial z}  \left[ \left( \nu + \frac{\cmu}{\sigma_{\varepsilon}} \frac{k^2}{\varepsilon} \right) \dervz{\varepsilon} \right]
  + c_{\varepsilon 1} \frac{\varepsilon}{k}
  \left( \eddyvis \left(\dervz{u}\right)^2  + c_{\varepsilon 3}
  \frac{\rho - \rho_w}{\rho_w} g_z \frac{1}{\sigma} \eddyvisC \dervz{c} \right) - C_{\varepsilon 2} \frac{ \varepsilon^2}{k}
$$

So far, I am using finite difference. Upwind scheme for the first order derivatives. Because the paramter $v_s$ is positive, concentration travels downwards at a speed $v_s$, I'm implementing the upwind scheme as:
$$
\frac{c^{i+1}_n - c^{i}_n}{\Delta z}
$$
for the node $i$. I also use the same upwind scheme for all the first-order derivative terms, i.e. also the terms inside the diffusion terms, which are nonlinear. For the linear diffusion terms, I apply second order central difference. For the time integration, I use a Crank-Nicholson method treating all terms implicitly.

The boundary conditions are complex. I'm using a standard wall function, using the log law at a node far from the wall. Using the log law, I also impose a gradient boundary condition because the log law adds the shear velocity as a unknown, hence the additional equation/condition. For
$k$ and $\varepsilon$ I use conditions that depend on the shear velocity. That's for one side of the domain, at the other side, sufficiently far, I use zero value conditions.

I'm not seeing much success with this implementation, should I use different discretization schemes for the diffusion and the advective terms? Should I use a different method altogether such as Finite Volume? What confuses me are the diffusion terms. Which discretization should I use there?

With regards to the time discretization. Should I treat some terms explicitely? The source terms are nonlinear as well.

Thanks in advance.

EDIT: Response to the first comment

Sorry I was not clear there. The model blows up for certain initial values of the concentration (too much acceleration), although that might have to do with the wall function implementation I think. Another issue that I had  were the initial conditions for $k$ and $\varepsilon$. I'll explain. The initial conditions for the velocity and concentration do not span the entire domain. They have positive values only for one tenth of the lower part of the domain, where I implement the wall function boundary condition, eventually, they will diffuse upwards. However, they will not do so if the values for $k$ and $\varepsilon$ are zero above where $u$ and $c$ are. This is because they are in the diffusion term and I think because of the upwind scheme, which assumes the wave travels downwards.


问:
Monte Carlo simulation
I am wondering if I am thinking correctly about the following problem :

Define the box of the dimensions $(a,a,H)$ in the $X$,$Y$, and $Z$ directions, respectively.

Insert $n$ particles into the box using PBC in $X$ and $Y$ directions

Move the particles using a quite analogous condition as in the previous pgms (of course, PBC in $X$ and $Y$ directions should be taken into account)

To test the pgm use $a=10$, $H=6$, $n$ of the order of $100-500$

I have generated random numbers for the $x$ and $y$ axis from the interval $[0,a]$ and for $z$ $[0,h]$.

Then, I have applied the PBC to the $x$ and $y$ axis and shifted them with the following rule: xt=x(j)+shift*2*rand() - shift and also applied PBC to the shifted position.

Xt is the new position of the point. For the y axis the operation is similar.
shift is the parameter in the interval $[0.1,1]$ but in this case, it is $0.5$

Am I thinking correctly about this problem, or there is some other way that I should take into account?

Title: Monte Carlo simulation


问:
Differences in answers between Python and Fortran
I am translating a piece of Fortran code into Python and am testing my code with a certain test case.

All the results differ with 0.04% compared to the Fortran results. This is a very small difference and I am not sure if I can almost see this as machine zero or if there are mistakes in my Python implementation.

What is odd is that all the values differ with exactly 0.04...%  from the Fortran results.


问:
Fortran solver for the Sparse LSE problem
I was wondering if there is a Fortran library that contains a solver for the Sparse LSE(linear equality-constrained least squares) problem
$$
min_{x}\|Cx-d\|^2  \text{ subject to }   Ax=b
$$
where $A$ and $C$ are large and sparse.

For dense linear systems, LAPACK has the subroutine DGGLSE, but I am looking for a solver for sparse systems.


问:
Finite difference scheme for solving nonlinear least-squares problem
I am dealing with following problem:

$$
\min_{u,\gamma}\Bigg\{ \frac{1}{1000}  \iint_{S_2} {\gamma (x,y)^2 dxdy} + \iint_{S_2} {[u(x,y)  - u_0 (x,y)]^2 dxdy} + \iint_{S_2} {[\Delta u(x,y) - \gamma (x,y) u(x,y)]^2 dxdy}\Bigg\}
$$

where the unknown functions $u(x, y)$ and $\gamma(x, y)$ are defined on the unit square $S_2$ and the function $u_0(x, y)$ is defined on $S_2$ by $u_0(x, y) = \sin(6\pi x) \sin(2\pi y)$.

This problem should be discretized using 5-point finite differences.

I have no experience discretizing such problems. Could somebody help me out, or point me to right literature to look for some examples?

Note: this problem comes from the paper Recursive trust-region methods for multiscale nonlinear optimization (page 10)


问:
fixed point iteration to find out second order non-linear diff equations
I am working on some model analysis, getting two diff equations and
after I convert them into matrix form, I have equations looks like

$$ [A][X]=C\times\big(\exp([B][X])-1\big), $$

where $C$ is a constant
and Both $[A]$ and $[B]$ are found out from two diff equations and boundary conditions.

Now my question is how am I able to find the value of $[X]$, I am thinking using fixed point method, but this is in matrix form, could any one give me any clue or hint?


问:
Are there any standardized file formats for point group character tables?
Character tables are an important tool for symmetry analysis in many computational chemistry software packages. Are there any standardized file formats for point group character tables?

This may seem to be an odd question as I could (and I do) just retrieve this information from other sources. But for code interoperability/sharability purpose, code documentation and just for laziness (as it is a very common task), it would be nice to have a standard file format for point group character tables.


问:
Calculating theoretical order of accuracy of least squares fit advection scheme
I'm familiar with finding the order of accuracy using von Neumann analysis for finite difference schemes formulated using Taylor series expansions.

But is there a similar technique for finding the order of accuracy for schemes that use a least squares fit for approximating gradients?

Currently, I am solving a 1D advection equation in advective form.  I'm using Forward Euler timestepping and an upwind quadratic least squares fit with a three-point stencil.  My experiments suggest an order of accuracy of ~0.7 which suggests that I've either got bugs in my implementation or the scheme is inherently slow to converge.


问:
Numerically solving a system of partial integro-differential equations in Matlab
Given the following system of partial integro-differential equations
$$ \frac{dX(t)}{dt}=\Lambda-\mu X(t)-\beta X(t)Z(t),\\
\frac{\partial Y(t,\omega)}{\partial t}+\frac{\partial Y(t,\omega)}{\partial \omega}=-\delta(\omega) Y(t,\omega),\\
\frac{dZ(t)}{dt}=-\gamma Z(t)+(N-Z(t))\int^{\infty}_{0}\kappa (\omega)Y(t,\omega)d\omega $$
with the boundary condition
$$Y(t,0)=\beta X(t)Z(t)$$ and
initial conditions
$$X(0)=X_{0},Y(0,\omega)= \sigma(\omega),Z(0)=Z_{0} $$
where all parameters $\Lambda, \mu, \beta, \delta(\omega), \gamma, N, \kappa(\omega), X_{0},Z_{0}$ are positive.

How does one numerically compute a solution for this on Matlab?

Note : I was asked to repost this question in scicomp.se -

Numerically solving a system of partial integro-differential equations in Matlab


问:
Time advance in Adaptive Mesh Refinement method
I am working on solving complex system of 2D PDEs governing the behaviour of plasma in a gas lamp during discharge. Recent tests have shown that because of steep gradients in temperature field and overall explosive nature of a process I won't go anywhere further without adding gas dynamics to the model. As far as I know from my mentor, there will be supersonic shocks traveling across the 5mm lamp. As far as I understand, to even think about somewhat accurate solution I first have to implement some kind of mesh refinement algorithm.

I have read Berger and Colella's LAMR for Shock Hydrodynamics, Colella's short course of block-structured AMR  at Berkley, Structured AMR by Vakili and Martin, Colella, Anghel and Alexander's AMR for Multiscale Nonequilibrium physics. So I understand that there have to be performed somewhat untrivial operations during the interpolation from coarse grid onto finer ones, synchronisations between halfsteps of finer grids, averaging of finer grids onto coarse ones and some kind of flux correction to preserve conservation laws.

That brings to the question: What exactly happens during the time step in Adaptive Mesh Refinement algorithm? I would appreciate full explanation, but here are some specific questions that are likely to arise either way:

Consider parabolic problem on two-level grid, $L_0$ & $L_1$, $L_0$ being the coarse one.

Is it correct that to get least amount of errors, I have to first advance L0 problem by time step $\tau$ as if there was no mesh refinement at all, then interpolate resulting $L_0$ data to get boundary conditions for $L_1$ problem at $t+\frac{\tau}{2}$ and advance $L_1$ for $\frac{\tau}{2}$, then interpolate L0 data for $L_1$ problem at $t + \tau$ and advance $L_1$ fol half-step again, then average $L_1$ data onto L0 grid and correct $L_0$ points adjacent to L1 boundary for flux conservation using $L_1$ data? Do I understand correctly that those flux fixups do not propagate further into coarse grid?

1.a. During correction of $L_0$ data using $L_1$, do I recalculate all $L_1$-bordering L0 points at the same time, or am I allowed to use corrected data on $L_0$ $X$-point alongside with data averaged from $L_1$ grid during my correcting of $Y$-point adjacent to $X$? (To illustrate: I am talking about something similar to difference between synchronous and asynchronous Hopfield networks).

1.b. If it is even possible, how do I synchronize grids for elliptic task? Radiative transfer is a large part of my system, too.

How do I generate boundary conditions on fine-fine grid boundary?

2.a. In some of the papers I've read I noticed the phrase similar to "If possible, boundary conditions are taken from the same-level grid or finer". How can this be, if at the moment of me updating the first of $L_k$ grids I do not yet have any data from the other $L_k$ or $L_{k+1}$ grids?

2.b. If the very concept of fine-fine boundary is not possible or leads to significant errors, then do I have to merge those contacting fine grids and solve the problem on the non-rectangular grid?

2.c. If I am to solve for the non-rectangular grid, what SLAE solvers besides iterative are good for sparse matrices? I was using Thomas algorithm modified for block-tridiagonal matrices, but as far as I understand, matrix for generic non-rectangular grid won't be tridiagonal.


问:
Is it generally unstable to use in multidimensional simulations finite difference schemes with higher orders than 2?
I'm part of a team trying to generalize a 1D Advection-Diffussion-Reaction code we inherited by extending it to 2D by using dimensional splitting, i.e. solving advection and diffusion for x and y separately with 1D schemes.

Diffusion was (relatively) painless but there are issues with advection.

Originally the code had a flux limited, third order upwind-biased scheme for 1D advection which worked like a charm in the 1D case. However, when the same scheme is extended to 2D, crazy instabilities develop. What supposed to be a radial advection-diffusion profile wrinkles in the transverse direction and shoots a "jet-like" instability:

[

However when I use a first order upwind scheme I do get a nice, symmetric radial profile.  The problem with the first order upwind scheme solution is the excessive false diffusion which the higher order scheme does suppress. But the higher order scheme is even more problematic because of the aforementioned instability.

I was wondering if higher order schemes in general are a bad idea in a multidimensional case? Are multidimensional problems more succeptible to "oscillations" caused by high order schemes? What about flux limiters, should 1D flux limiters prevent this oscillations from happening (in my situation, it seems the flux limiter for the third order case doesn't help)? Should I stick to lower order schemes even if they bring false diffusion? This third order scheme I was using is also flux limited, but the flux limit doesn't stop the instabilities...


问:
Method for calculating the stopping distance using only integers
I am trying to find how much an engine would turn until it stops accelerating with an acceleration of $a = A \sin^2(x)$.

Acceleration

Integrating the acceleration I get velocity $v = A(x/2 - \sin(2x)/4) $

Integrating the velocity, I get the distance $d = A(x^2/4 - \cos^2(x)/4) $

Velocity and distance

Now, I need to do this calculation on a microcontroller, without using floating point and with a sine look-up table which is basically a $2048\sin(x\pi/512)$ function.

How can I find the distance correctly using the sine look up table? Or what other way could I use?


问:
Reducing oscillations a 3D Alternating direction explicit scheme for the diffusion equation?
Hi I have made a 3D alternating direction explicit scheme for solving the diffusion equation, which will eventually replace a FTCS scheme in model of bubble dynamics in tissue. I have been testing it and have realized that despite the unconditional stability I get oscillations in the system if the Courant number is  greater than ~4. Is there a way round this other than reducing the Courant number?

Here are some more details: The scheme I am using alternating direction explicit ADE is an explicit method which does two sweeps (one upwind and one downwind). In 1D this looks like:

Upwind:

```<code>C1_new[i]-C_old[i] = (dt/h^2)(C_old[i+1]-C_old[i]-C1_new[i]+C1_new[i-1])
</code>```

Downwind:

```<code>C2_new[i]-C_old[i] = (dt/h^2)(C_old[i-1]-C_old[i]-C2_new[i]+C2_new[i+1])
</code>```

Final Average:

```<code>C_new[i] = C1_new[i]+C2_new[i]/2
</code>```

Full details can be found at http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1552926

However, when I implement a 3D version and play around with time and space steps I see oscillations occurring from (dt/h^2) values of about 4.5 upwards.

This is a plot of the pt (5,5,5) from my $10\times10\times10$ grid for various Courant numbers. As you can see, there are oscillations for the higher numbers.

]


问:
Non-overlaping Domain decomposition - assemble of Laplacian
I am dealing with following 2-dimensional problem in the unit square domain $S_2$

$$- \Delta u (x,y) = f \  \text{in} \ S_2,  \hspace{1.5cm} u(x,y) = 0 \ \text{on} \  \partial S_2$$

where $f$ is such the analytical solution to problem is

$$u(x,y) = \cos [  2 \pi x  (y - x) ]^2.$$

Problem should be discretized by 5-point finite difference scheme, which will result in Ax = b system.
I need to decompose problem into 4 non-overlaping domains. Should I divide  my domain square into 4 parts and then assemble? or assemble and then divide? Main question here is, how the laplacians will look like.


问:
What is the source of the error in the Sherman-Morrison formula application?
The Sherman-Morrison formula
$$ (A+uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u} $$
results in small errors in relation to the standard matrix inverse operation after each application, as shown here. From what I understand, the Sherman-Morrison identity is exact, so the error must come from numerical problems. But where exactly? Also, are there any theoretical bounds for this error? Some sources argue that the error is related to the matrix condition number, but wouldn't that affect the standard inversion too?</p>


问:
C++ Library: What is the Most Efficient Library that Factorizes a Polynomial?
I need to know what library is the fastest one that can factorize a polynomial of large degree whose coefficients are big integers. </p>

<p>I have tried NTL library to factorize a monic polynomial but it is not efficient. </p>


问:
How to compute this double integral?
Let

$$T=1, K=100, S_0=100, \sigma=0.05, r=0.15. $$

Define $\nu:=\frac{2r}{\sigma^2}-1$
and $$H(y,z)=\frac{z e^{\pi^2 /4y}}{\pi \sqrt{\pi y}}\int_0^{\infty} e^{-z \cosh(u) -u^2/(4y)} \sinh(u) \sin(\frac{\pi u}{2y})du$$
then next define
$$ f(y,z)=\left( \frac{z}{S_0}\right)^{\nu/2}\frac{1}{4KT}\exp\left(-\frac{2(S0+z)}{KT\sigma^2}-\frac{\nu^2\sigma^2y}{8} \right) \times H(\frac{\sigma^2y}{8},\frac{4\sqrt{S_0 z}}{KT\sigma^2}) $$
and $$ g(y,z)=z\frac{e^{-ry}-e^{-rT}}{rT} $$
I would like to compute the following double integral
$$ I=\int_0^{\infty}\left[\int_0^T g(y,z)f(y,z)dy\right]dz. $$

When I tried this with Maple, please see the picture below,

it ran out of memory.

I would like to ask whether there is a way to compute the above integral? What should I do to compute the integral?

Title: How to compute this double integral?


问:
Symplectic Partitioned Runge Kutta method in Mathematica
I tried to solve Hamiltonian system ($Q$ is a vector of all generalized coordinates, $P$ - of generalized momentum)
$$
\frac{\mathrm{d} Q}{\mathrm{d} t}=\frac{\partial H}{\partial P} \\
\frac{\mathrm{d} P}{\mathrm{d} t}=-\frac{\partial H}{\partial Q}
$$
where
$$
H=\frac{1}{2} \sum_{i=1}^{n}\overrightarrow{p}_{i}^{2} + \frac{k}{2} \sum_{i=2}^{n}(\left \| \overrightarrow{q}_{i}-\overrightarrow{q}_{i-1} \right \|-rest)^2 + \frac{\kappa}{2} \sum_{i=2}^{n-1}\arccos^2\frac{(\overrightarrow{q}_{i}-\overrightarrow{q}_{i-1})\cdot (\overrightarrow{q}_{i+1}-\overrightarrow{q}_{i})}{\left \| \overrightarrow{q}_{i}-\overrightarrow{q}_{i-1} \right \|\left \| \overrightarrow{q}_{i+1}-\overrightarrow{q}_{i} \right \|}
$$
with "SymplecticPartitionedRungeKutta" method numerically using Mathematica (description could be found here http://www.sciencedirect.com/science/article/pii/S0377042701004927)

The code is below

```<code>\[CurlyKappa] = 20;
k = 20;
rest = Sqrt[5];
n = 3;
dim = 3;
Q = Table[Table[Subscript[q, j][i][t], {j, 1, dim}], {i, 1, n}];
P = Table[Table[Subscript[p, j][i][t], {j, 1, dim}], {i, 1, n}];
icsQ = {{-1, -1, 0}, {0, 1, 0}, {1, -1, 0}};
icsP = Table[Table[0, {j, 1, dim}], {i, 1, n}];

H = 1/2 (Sum[P[[i]].P[[i]], {i, 1, n}] +
     Sum[(Norm[Q[[i]] - Q[[i - 1]]] - rest)^2, {i, 2, n}]*k +
     Sum[(VectorAngle[Q[[i]] - Q[[i - 1]],
         Q[[i + 1]] - Q[[i]]])^2, {i, 2, n - 1}]*\[CurlyKappa]);
Q' = Flatten[D[#, t] & /@ Q];
P' = Flatten[D[#, t] & /@ P];

Conjugate'[y_[x_][t_]] := 1;
Abs'[x_] := Sign[x];

var = Flatten[Q~Join~P];
ics = Flatten[icsQ~Join~icsP];
eqn = Table[Q'[[i]] == D[H, var[[n*dim + i]]], {i, 1, n*dim}]~Join~
   Table[P'[[i]] == -D[H, var[[i]]], {i, 1, n*dim}]~Join~
   Table[(var[[i]] /. t -> 0) == ics[[i]], {i, Length[var]}];

method = {"SymplecticPartitionedRungeKutta",
   "PositionVariables" -> var[[1 ;; dim*n]]};
sol = NDSolve[eqn, var, {t, 0, 10}, Method -> method,
  WorkingPrecision -> 10, MaxStepSize -> 0.001,
  MaxSteps -> \[Infinity]]
</code>```

However, error pops out:

NDSolve::sprksep: "The Hamiltonian of the differential system in the
  method does not appear to be in separable form. Try using the method
  ImplicitRungeKutta with coefficients
  ImplicitRungeKuttaGaussCoefficients."

But it's easy to see that the Hamiltonian of the system definately has separate form H(p,q)=T(p)+V(q).

Can anyone help?

P.S. In order to prevent questions. You could notice that I defined two "strange" functions:

```<code> Conjugate'[y_[x_][t_]] := 1;
 Abs'[x_] := Sign[x];
</code>```

Idea: Mathematica doesn't take the derivative from non-analytical functions. And when you take derivative in form D[H,q[t]], terms of type D[Conjugate[q[t]],q[t]] appear. If I don't fix it, I will get an error:

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0


问:
Maximizing slow multi-parameter function
What's the best method to maximize the value of slow multivariate function? Here is what I know about the function:

Number of parameters is ~ 10.

It takes considerable amount of time to compute the value of the function for single variable set. Let's say 30 minutes.

Function seems to be concave against every individual parameter.

I can run computation in parallel.

All my code is in Java and I'm interested in coding crude version of the optimization algorithm myself rather than using some industry solver.


问:
"Tunneling" optimization algorithm in Matlab
I was wondering if someone has implemented the "Tunneling algorithm" for the global minimization of a single variable function in Matlab.

I am hoping to implement it on $[B, A - \lambda I]$. Where $A$ is and $N\times N$ matrix, B is an $N\times M$ matrix, $I$ is the identity and $\lambda$ is the unknown variable in the complex domain.

This problem is related to "the distance to uncontrollability of a dynamic system".


Title: How to adjust scikit-learn FactorAnalysis() settings to get scores output similar to one from factanal() in R?

Body:
I'm new to factor analysis and I need to compute factor scores using Python. I have R code to compute scores and I want to adjust my Python code (decomposition.FactorAnalysis() settings) to get similar output.

R code:

```<code>data <-read.csv("test.csv", header=F)
col_names <- names(data)
drops <- c()

for (name in col_names){
  st_dev <- sd(data[,name], na.rm = T)
  if (st_dev == 0){
    drops <- c(drops, name)
  }
}

da_nal <- data[,!(names(data) %in% drops)]
factanal(na.omit(da_nal), factors = 1, scores = 'regression')$scores
</code>```

Its output:

```<code>    Factor1
1   4.89102190
2   3.65004187
3   0.14628700
4  -0.20255897
5  -0.01565570
6  -0.16438863
7   0.40835986
8  -0.25823984
9  -0.20813064
10  0.09390067
11 -0.28891296
12 -0.28882753
13 -0.26624358
14 -0.25202275
15 -0.25181326
16 -0.15653679
17 -0.28702281
18 -0.28865654
19 -0.23251509
20 -0.28066125
21 -0.18714387
22 -0.24969113
23 -0.28302552
24 -0.28712610
25 -0.29196529
26 -0.28659988
27 -0.29502523
28 -0.15802910
29 -0.27440118
30 -0.29083667
31 -0.29548220
32 -0.29461059
33 -0.23594859
34 -0.29654336
35 -0.29759659
36 -0.29085001
37 -0.29539071
38 -0.29234303
39 -0.29702103
40 -0.27595130
41 -0.27184361
</code>```

Python code:

```<code>from sklearn import decomposition, preprocessing
import numpy as np

data = np.genfromtxt('rangir_test.csv', delimiter=',')
data = data[~np.isnan(data).any(axis=1)]
data_normal = preprocessing.scale(data)
fa = decomposition.FactorAnalysis(n_components = 1)
fa.fit(data_normal)
for score in fa.score_samples(data_normal):
    print score
</code>```

Its output is quite different from R's, but want I them to be similar.

Python output:

```<code>-69.8587183816
-116.353511148
-24.1529840248
-36.5366398005
-7.87165586175
-24.9012815104
-23.9148486368
-10.047780535
-4.03376369723
-7.07428842783
-7.44222705099
-6.25705487929
-13.2313513762
-13.3253819521
-9.23993173528
-7.141616656
-5.57915693405
-6.82400483045
-15.0906961724
-3.37447211233
-5.41032267015
-5.75224753811
-19.7230390792
-6.75268922909
-4.04911793705
-10.6062761691
-3.17417070498
-9.95916350005
-3.25893428094
-3.88566777358
-3.30908856716
-3.58141292341
-3.90778368669
-4.01462493538
-11.6683969455
-5.30068548445
-24.3400870389
-7.66035331181
-13.8321672858
-8.93461397086
-17.4068326999
</code>```

Here is the data itself:

```<code>49.6,34917,24325.4,305,101350,98678,254.8,276.9,47.5,1,3,5.6,3.59,11.9,0,97.5,97.6,8,10,100,0,0,96.93,610.1,100,1718.22,6.7,28,5
275.8,14667,11114.4,775,75002,74677,30,109,9.1,1,0,6.5,3.01,8.2,1,97.5,97.6,8,8,100,0,0,100,1558,100,2063.17,5.5,64,5
2.3,9372.5,8035.4,4.6,8111,8200,8.01,130,1.2,0,5,0,3.33,6.09,1,97.9,97.9,8,8,67.3,342.3,0,99.96,18.3,53,1457.27,4.8,8,4
7.10,13198.0,13266.4,1.1,708,695,6.1,80,0.4,0,4,0,3.1,8.2,1,97.8,97.9,8,8,45,82.7,0,99.68,4.5,80,1718.22,13.8,0,3
1.97,2466.7,2900.6,19.7,5358,5335,10.1,23,0.5,0,2,0,3.14,8.2,0,97.3,97.2,9,9,74.5,98.2,0,99.64,79.8,54,1367.89,6.4,12,4
2.40,2999.4,2218.2,0.80,2045,2100,8.9,10,1.5,1,3,0,2.82,8.6,0,97.4,97.2,8,8,47.2,323.8,0,99.996,13.6,24,1249.67,2.7,12,3
0.59,4120.8,5314.5,0.54,14680,13688,14.9,117,1.1,0,3,0,2.94,3.4,0,97.6,97.7,8,8,11.8,872.6,0,100,9.3,52,1251.67,14,14,2
0.72,2067.7,2364,3,367,298,7.2,60,2.5,0,12,0,2.97,10.5,0,97.5,97.6,8,8,74.7,186.8,0,99.13,12,57,1800.45,2.7,4,2
1.14,2751.9,3066.8,3.5,1429,1498,7.7,9,1.6,0,3,0,2.86,7.7,0,97.6,97.8,8,9,76.7,240.1,0,99.93,13.6,60,1259.97,15,8,3
1.29,4802.6,5026.1,2.7,7859,7789,6.5,45,1.9,0,3,0,2.5,8.2,0,98,98,8,8,34,297.5,0,99.95,10,30,1306.44,8.5,0,4
0.40,639.0,660.3,1.3,23,25,1.5,9,0.1,0,1,0,2.5,8.2,0,97.7,97.8,8,8,94.2,0,0,100,4.3,50,1565.44,19.2,0,4
0.26,430.7,608.1,2,33,28,2.5,7,0.4,0,6,0,2.5,8.2,0,97.4,97.4,8,8,76.5,0,0,98.31,8,54,1490.08,0,0,4
4.99,2141.2,2357.6,3.60,339,320,8.1,7,0.2,0,8,0,2.5,5.9,0,97.3,97.4,8,8,58.1,206.3,0,99.58,13.2,95,1122.92,14.2,8,2
0.36,1453.7,1362.2,3.50,796,785,3.7,9,0.1,0,9,0,2.5,13.6,0,98,98.1,8,8,91.4,214.6,0,99.74,7.5,53,1751.98,11.5,0,2
0.36,1657.5,2421.1,2.8,722,690,8.1,8,0.4,0,1,0,2.5,8.2,0,97.2,97.3,11,12,37.4,404.2,0,99.98,10.9,35,1772.33,10.2,8,3
1.14,5635.2,5649.6,3,2681,2530,5.4,20,0.3,0,1,0,3.1,8.2,0,97.7,97.8,8,11,50.1,384.7,0,99.02,11.6,27,1306.08,16,0,2
0.6,1055.9,1487.9,1.3,69,65,2.5,6,0.4,0,8,0,2.5,8.2,0,97.9,97.7,8,11,63,137.9,0,99.98,5.1,48,1595.06,0,0,4
0.08,795.3,1174.7,1.40,85,76,2.2,7,0.2,0,0,0,2.5,8.2,0,97.4,97.5,8,8,39.3,149.3,0,98.27,5.1,52,1903.9,8.1,0,2
0.90,2514.0,2644.4,2.6,1173,1104,5.5,43,0.8,0,10,0,2.5,13.6,0,97.5,97.5,8,10,58.7,170.5,0,80.29,10,34,1292.72,4,0,2
0.27,870.4,949.7,1.8,252,240,2.2,31,0.2,0,1,0,2.5,8.2,0,97.5,97.6,8,8,64.5,0,0,100,6.6,29,1483.18,9.1,0,3
0.41,1295.1,2052.3,2.60,2248,2135,6.0,12,0.8,0,4,0,2.7,8.2,0,97.7,97.7,8,8,71.1,261.3,0,91.86,4.6,21,1221.71,9.4,0,4
1.10,3544.2,4268.9,2.1,735,730,6.6,10,1.7,0,14,0,2.5,8.2,0,97.7,97.8,8,8,52,317.2,0,99.62,9.8,46,1271.63,14.2,0,3
0.22,899.3,888.2,1.80,220,218,3.6,7,0.5,0,1,0,2.5,8.2,0,97.2,97.5,8,8,22.5,0,0,70.79,10.6,32,1508.02,0,0,4
0.24,1712.8,1735.5,1.30,41,35,5.4,7,0.5,0,1,0,3.28,8.2,0,97.8,97.8,9,10,16.6,720.2,0,99.98,4.3,53,1324.46,0,4,2
0.2,558.4,631.9,1.7,65,64,2.5,7,0.2,0,5,0,2.5,8.2,0,97.7,97.5,8,8,60.7,0,0,99.38,6.1,52,1535.08,0,0,2
0.21,599.9,1029,1.1,69,70,3.7,85.7,0.1,0,12,0,2.5,8.2,0,97.4,97.5,8,8,48.6,221.2,0,100,5.4,40,1381.44,25.6,0,2
0.10,131.3,190.6,1.6,28,25,2.9,7,0.3,0,3,0,2.5,8.2,0,97.7,97.8,8,8,58.9,189.4,0,99.93,6.9,42,1525.58,17.4,0,3
0.44,3881.4,5067.3,0.9,2732,2500,11.2,10,1.5,0,5,0,2.67,8.2,0,97.4,97.3,8,11,14.5,1326.2,0,99.06,3.7,31,1120.54,10.3,10,2
0.18,1024.8,1651.3,1.01,358,345,4.6,35,0.3,0,2,0,2.5,8.2,0,97.8,97.9,8,10,15.9,790.2,0,100,4.3,48,1531.04,10.5,0,3
0.46,682.9,784.2,1.8,103,109,2.2,8,0.4,0,4,0,2.5,8.2,0,97.8,97.9,8,8,82.7,166.3,0,99.96,6.4,44,1373.6,13.5,0,2
0.12,370.4,420.0,1.10,28,25,3.4,10,0.1,0,6,0,2.57,8.2,0,97.6,97.8,8,11,51.6,120,0,99.85,8.1,40,1297.94,0,0,3
0.03,552.4,555.1,0.8,54,49,3.5,10,0.4,0,0,0,2.5,8.2,0,97.4,97.6,8,10,33.6,594.5,0,100,3.2,41,1184.34,6.6,0,3
0.21,1256.5,2434.8,0.9,1265,1138,6.3,20,1.3,0,2,0,2.6,8.2,0,98,97.9,8,9,20.1,881,0,99.1,3.9,31,1265.93,7.8,0,3
0.09,320.6,745.7,1.10,37,25,2.7,8,0.3,0,9,0,2.5,8.2,0,98,97.8,8,8,49.2,376.4,0,99.95,4.3,39,1285.11,0,0,3
0.08,452.7,570.9,1,18,20,4.7,9,0.6,0,2,0,2.45,8.2,0,97.1,97.1,8,8,19.9,1103.8,0,99.996,2.9,22,1562.61,21.9,0,3
0.13,967.9,947.2,1,74,65,4.0,25,1.4,0,6,0,2.5,8.2,0,98,98,9,11,30.1,503.1,0,99.999,3.4,55,1269.33,0,0,2
0.07,495.0,570.3,1.2,27,30,4.3,7,0.5,0,12,0,3.62,8.2,0,98.2,98.2,15,13,29.8,430.5,0,99.7,4.9,40,1461.79,14.6,0,2
0.17,681.9,537.4,1.1,113,120,2.9,12,0.4,0,8,0,2.5,8.2,0,98.2,98.3,8,8,24,74.3,0,100,5,43,1290.16,0,0,3
0.05,639.7,898.2,0.40,9,12,3.0,7,0.1,0,1,0,2.5,8.2,0,97.6,97.8,15,11,11.9,1221.1,0,99.996,1.7,40,1372,7,0,4
0.65,2067.8,2084.2,2.50,414,398,7.3,6,0.7,0,4,0,2.16,8.2,0,97.8,97.9,12,12,60.1,146.3,0,99.96,10.4,44,1059.68,7.4,0,2
0.12,804.4,1416.4,3.30,579,602,4.2,7,1.8,0,1,0,2.5,8.2,0,98.1,98.3,8,10,8.9,2492.3,0,95.4,2.2,34,1345.76,7,0,2
</code>```

P.S. Don't have enough rep to create factor-analysis tag.


问:
Broadening spectral data by using FFT's
I obtain numerical discrete data of the form

$$
S_{raw}(\omega) = \sum_{j}w_{j} \delta(\omega-\omega_{j})
$$

to compare the result with experimental data the delta peaks need to be broadened according to

$$
S_{broad}(\omega) = \int d\omega^{\prime} K(\omega,\omega^{\prime}) S_{raw}(\omega^{\prime})
$$

where $$ K(\omega, \omega^{\prime})$$ is the broadening kernel i.e.

$$
K(\omega,\omega^{\prime}) = \frac{\theta(\omega\omega^{\prime})}{\sqrt{\pi}b|\omega|}\exp\left[-\left(\frac{1}{b}\ln\left|\frac{\omega}{\omega^{\prime}}\right| - \frac{b}{4}\right)^2\right].
$$

If $$K(\omega,\omega^{\prime}) = K(\omega-\omega^{\prime})$$ the integral equation above can be interpreted as a convolution and can be efficiently calculated by means of the convolution theorem and FFTs. But the example kernel $K$ does not have this property.

I ask me now if there is still a way to compute the broadened function by exploiting the convolution theorem and fast FFT algorithms ?

Title: Broadening spectral data by using FFT's


问:
Software to simulate behavioral economic models?
I am looking for a way to simulate behavioral economic models to illustrate the effects of changing parameters. I am familiar with statistical analysis (especially R), but I only have few programming skills beyond R. I think that if there is a way to achieve my goal (I'm not really sure if it is possible at all!) that Octave should be a good choice.

Usually such a model defines some kind of game with actions, strategies and payoffs for every player. Then these actions get evaluated and normally an equilibrium between the players actions/strategies is determined.

One of my desired models is Rabin's idea of Fairness Equilibria. Pages 1286ff can give a short idea of how Rabin defines the game and determines the equilibrium. LINK

My goal is to have a function that depicts that model and enables me to illustrate the effects of different assumptions and situations within the frame of the specific model.

So my question is: Is it possible at all? Is Octave suitable or should I choose another language/environment?

Thanks in advance for any sort of advice!

EDIT:

I try to sketch the model:

Starting point is a two player normal form game (Prisoner's Dilemma) with the strategy sets $S_1, S_2$ for every player from which they choose their actions. The payoff function for each player is defined as $\pi_i = S_1 \times S_2$. So player $i$ chooses action $a_i \in S_i$. Now player $j$ has a believe about $i$'s strategy $b_i \in S_i$ and player $i$ has a belief about player $j$'s belief about $i$'s strategy $c_i \in S_i$.
So the feasible payoffs for $i$ depending on his beliefs are in the set $\Pi(b_j) = \{(\pi_i(a_i,b_j),\pi_j(b_j,a_i)) | a_i \in S_i\}$

That's the framework.

The routine now has to determine

$\pi_j^h(b_j)$ as the highest payoff for $j$ depending on $i$'s belief that is Pareto efficient in $\Pi(b_j)$

$\pi_j^l(b_j)$ as the lowest payoff for $j$ depending on $i$'s belief that is Pareto efficient in $\Pi(b_j)$

$\pi_j^{min}(b_j)$ as the worst possible payoff for $j$ depending on $i$'s belief

the reference point for $j$: $\pi_j^e(b_j) = \frac{\pi_j^h(b_j) + \pi_j^l(b_j)}{2}$

the kindness of $i$ towards $j$: $\kappa_i(a_i,b_j)=\frac{\pi_j(b_j,a_i) - \pi_j^e(b_j)}{\pi_j^h(b_j) - \pi_j^{min}(b_j)}$

believed kindness of $j$ towards $i$: $\bar\kappa_j(b_j,c_i)=\frac{\pi_i(c_i,b_j) - \pi_i^e(c_i)}{\pi_i^h(c_i) - \pi_i^{min}(c_i)}$

Then the Expected Utility is defined as: $U_i(a_i,b_j,c_i) = \pi(a_i,b_j) + \bar\kappa_j(b_j,c_i) \cdot [1+\kappa_i(a_i,b_j)]$

Now, to be a Fairness Equilibrium, a pair of strategies $(a_1^*,a_2^*)$ has to fulfill following conditions:

$a_i^* \in \text{argmax}_{a_i \in S_i} U_i(a_i,b_j,c_i)$

$c_i = b_j = a_i^*$

So my goal is to get a routine that determines these Fairness Equilibria when a certain game is given.


问:
Solving an LP greedily
I have the following LP:
$$
\begin{array}{ll}
\text{Minimize} & \sum_{j=1}^n x_j \\
\text{Subject to} & \sum_{j=1}^n a_{ij} x_j \geq b_i,~~~i\in\{1,\ldots,M\} \\
                  & 0 \leq x_j \leq 1,~~~j\in\{1,\ldots,n\}
\end{array}
$$
where all $a_{ij} \geq 0, \forall i,j$ and $b_i \geq 0,\forall i$, and these are known.

I am looking for a quick direct solution for this LP if it exists. We know that $n$ is much smaller than $M$. This seems similar to the fractional multi-Knapsack problem. I am wondering if there exists a greedy solution that I can use.

Any help is much appreciated.

Thank you

Title: Solving an LP greedily


问:
Fast library for parallel n-diagonal matrix inversion
I've tried looking for efficient parallel algorithms to solve banded systems. SPIKE is a good example but it however requires that the band size is large.

Is there any efficient library implementing a parallel algorithm for n-diagonal systems (where n is less than 10)

I haven't found any implementations of Thomas' algorithm which seems promising.


问:
Inverted value is not consistent with expectation
We have a group of observations

$$y = f(x_1, x_2, x_3) \enspace .$$

We have also a forward model $y = f(x_1, x_2)$. The forward model does not include $x_3$ because $x_3$ might include dozens of parameters which we do not understand clearly.

First, we inverted parameter $x_1$ with the forward model $y = f(x_1, 0)$ embedded in our inversion strategy. The best model can fit 80% of the observed data sets. Then we inverted both parameters $x_1$ and $x_2$ with the forward model $y = f(x_1, x_2)$. However, the misfit did not reduce and the inverted $x_2$ is always quite close to zero.

Based on previous laboratory and numerical experiments, $x_2$ should have significant effects on the observation and its value is obvious larger than zero. So we fix $x_2$ as the experimental observed value and inverted for the optimized model. The misfit from the best inverted model is larger than that with $x_2$ fixed at zero. How can we understand these inversion results? Do they mean that $x_2$ is balanced by $x_3$ or that there was overfitting in inversions with $x_1$ only?