问:
How to express the polynomial in an graded lex (with the descending power of multiple variables)?
```<code>from sympy import *
from sympy.abc import A,B,C,D,E,F,x,y,U,V,W,X,Y,Z

curve = A*x**2+B*x*y+C*y**2+D*x+E*y+F
</code>```

Now that we know such codes above will form the polynomial curve like
$$A x^{2} + B x y + C y^{2} + D x + E y + F$$
However, the six terms of the polynomial are sorted only relevant to the upper cases. To be sure of this, we could just substitute the coefficients like:

```<code>curve.subs({A:Z, B:Y, C:X, D:W, E:V, F:U})
</code>```

And then the output polynomial comes to be totally out of the order of x, y like
$$U + V y + W x + X y^{2} + Y x y + Z x^{2}$$
Then, we know that class Poly can be used to pull x, y out, making the terms of the polynomial arranged only by the times of x and y, and arrange x first, then arrange y when the terms include x to the same power.

```<code>Pcurve = curve.as_poly(x,y)
</code>```

$$\operatorname{Poly}{\left( A x^{2} + B xy + D x + C y^{2} + E y + F, x, y, domain=\mathbb{Z}\left[A, B, C, D, E, F\right] \right)}$$

The way to keep sorting is easy, just redefine x and y via

```<code>x,y = symbols('\ x, \ y')
</code>```

where the spaces with the smaller ASCII code of 32, give the x and y priority. However, like the Poly result,  x and  y are not in the same priority. the  x is preferred.
$$ x^{2} Z +  x  y Y +  x W +  y^{2} X +  y V + U$$

Then, how to just collect the terms, while the times of x and y can be treated together, to make $Xy^2$ in front of $Wx$, while the terms are sorted by the preferred x, y, which are presented after the coefficients?

In fact, what I really did is

```<code>x,y = symbols('\ \ x, \ \ y')
x0, y0 = symbols('\ x_0  \ y_0')

# dpp = Symbol('(x-x0)**2 + (y-y0)**2')
dpp = sqrt((x-x0) ** 2 + (y-y0) ** 2)
dpl = (A * x + B * y + C)

elps = x**2 / a**2 + y ** 2 / b ** 2 - 1
Lagrange = Add(dpp ** 2, - e ** 2 * dpl ** 2, - lamda * elps, evaluate=False)
Lpoly = Poly(Lagrange, (x,y))
Lag_exp = Lagrange.expand()
Lag_coll = Lag_exp.collect([x,y])
</code>```

Where the last four variables are assigned by

1.
$$- e^{2} \left(  x A +   y B + C\right)^{2} - \lambda \left(\frac{  x^{2}}{a^{2}} + \frac{  y^{2}}{b^{2}} - 1\right) + \left(\left(  x -  x_{0}\right)^{2} + \left(  y -  y_{0}\right)^{2}\right)$$

2.
$$\operatorname{Poly}{\left( \frac{- A^{2} a^{2} e^{2} + a^{2} - \lambda}{a^{2}}   x^{2} -  2 A B e^{2}   x  y + \left(- 2  x_{0} - 2 A C e^{2}\right)   x + \frac{- B^{2} b^{2} e^{2} + b^{2} - \lambda}{b^{2}}   y^{2} + \left(- 2  y_{0} - 2 B C e^{2}\right)   y +  x_{0}^{2} +  y_{0}^{2} - C^{2} e^{2} + \lambda,   x,   y, domain=\mathbb{Z}\left(a, b, e,  x_{0},  y_{0}, A, B, C, \lambda\right) \right)}$$

3.
$$-   x^{2} A^{2} e^{2} +   x^{2} - \frac{  x^{2} \lambda}{a^{2}} - 2   x   y A B e^{2} - 2   x  x_{0} - 2   x A C e^{2} -   y^{2} B^{2} e^{2} +   y^{2} - \frac{  y^{2} \lambda}{b^{2}} - 2   y  y_{0} - 2   y B C e^{2} +  x_{0}^{2} +  y_{0}^{2} - C^{2} e^{2} + \lambda$$

4.
$$  x^{2} \left(- A^{2} e^{2} + 1 - \frac{\lambda}{a^{2}}\right) +   x \left(- 2   y A B e^{2} - 2  x_{0} - 2 A C e^{2}\right) +   y^{2} \left(- B^{2} e^{2} + 1 - \frac{\lambda}{b^{2}}\right) +   y \left(- 2  y_{0} - 2 B C e^{2}\right) +  x_{0}^{2} +  y_{0}^{2} - C^{2} e^{2} + \lambda$$

And my ideal order for the letters are: x = y > a = b = A = B > C > e > λ. The normal 'lex' or 'inv-lex' (Lagrange.expand().as_ordered_terms('rev-lex')) both do not fit me.

I had a simple and crude way to realise MY goal, but still not a general solution. Could it be enhanced?

```<code>from sympy import *
from sympy.abc import A,B,C,D,E,F,a,b,e,x,y,lamda
x0, y0 = symbols('x_0  y_0')

dpp = sqrt((x-x0) ** 2 + (y-y0) ** 2)
dpl = (A * x + B * y + C)
elps = x**2 / a**2 + y ** 2 / b ** 2 - 1

Lagrange = dpp ** 2 - e ** 2 * dpl ** 2 - lamda * elps
Lexpd = Lagrange.expand()

standard_curve = A*x**2 + B*x*y + C*y**2 + D*x + E*y + F

Ldict = {}
for i in [x**2, x*y, y**2, x, y]:
    coeff = Lexpd.coeff(i)
    if len(coeff.as_terms()[0]) > 1:
        Ldict[standard_curve.coeff(i)] = '\\left(' + latex(Lexpd.coeff(i))+ '\\right)'
    else:
        Ldict[standard_curve.coeff(i)] =  latex(Lexpd.coeff(i))
    Lexpd -= Lexpd.coeff(i) * i
Ldict[standard_curve.coeff(1)]= latex(Lexpd.expand())

print(latex(standard_curve, symbol_names=Ldict))
</code>```


问:
How can I optimize the solution of linear systems in my Finite Element Method code for solving nonlinear computational mechanics problems?
Currently, I use Eigen3 for linear algebra operations on sparse matrices and either UMFPACK or CHOLMOD from SuiteSparse to solve sparse linear systems. However, as my model grows larger and the need for iterative methods to solve nonlinear equations arises, I am facing a bottleneck in the code's performance, particularly in the solution of linear systems at each iteration.

I would like to improve computational efficiency and reduce execution time by implementing a parallel solver. However, I'm unsure if libraries like PETSc and Trilinos are suitable for use on laptops or if they are primarily designed for large clusters. I would appreciate recommendations for parallel solver libraries that are well-suited for solving sparse linear systems efficiently on a laptop, taking advantage of the available resources such as my i7-11390H processor and 16GB RAM, to achieve optimal performance.

ContentLicense: CC BY-SA 4.0


问:
Why are the structures in Gaussian software contracting while performing optimization calculation?
I was performing the optimization for the following Ru metal complexes using DFT method in Gaussian software. The structures contract and the calculation terminates with convergence failure error.
I am using LanL2Dz for Ruthenium and 3-21G for other atoms. 3-21G basis set as beginning step. I would like to use 6-31G+(d,p) further.
Could anyone please help how I can resolve this issue?

Thank you.


问:
How to derive a specific Gibbs sampling scheme?
I have a bio-informatics related question regarding gibbs sampling and how to derive it. ANY help would be highly appreciated as i have no clue at all how to start.


问:
How to calculate the change in water level height after each step in a water displacement simulation?
I'm trying to write a very simple simulation of volume/water displacement. Every timestep a sphere gets lowered a bit further into a cylindrical container filled with water.
While the physics behind this should be quite clear, I just can't wrap my head around how I would go about calculating the change in water height per timestep.

If the sphere is lowered by $y$ below the initial water level, then the water level will rise by $h$, resulting in the sphere being submerged by $y + h$, which would in turn mean it would displace more water, making the water level rise some more and so on and so on...
This of course is wrong somewhere, but I fail to come up with a better idea. So, where do I err?

Assume $r_{C}$to be the radius of the cylindrical container, $h_{0}$the
original water level without a sphere in it, $h_{n}$ the water level
after the nth step $r_{S}$the radius of the sphere and $y_{n}$ how
high the spherical cap of the submerged part (below the original water
level) of the sphere is after step $n$.

From the balance of volumes before and after insertion for step 1 we get $V_{o}=V_{1}-V_{S}$ with $V_{0}$ being the original water volume, $V_{1}$ the volume after inserting part of the sphere in the first step, and $V_{S}$ the volume of the inserted spherical cap.

$$V_{0}=r_{C}^{2}\pi h_{0}\\
V_{1}=r_{C}^{2} \pi h_{1}\\
V_{S}=\frac{\pi}{3}h_{1}^{2}(3h_{1}-h_{1})$$

this is fine on its own and can be solved by solving a cubic equation,
but nowhere in this does $y_{n}$ appear, and that's the part where
I get hung up. How to transform these equations into something dependent
on $y_{n}$, e. g. how to transform them into something dependent
on how much of the sphere is below the original water level?


问:
Has the arithmetic for exotic (unsigned float with positive exponent) number format been solved?
The data type is a doubly unsigned float. This is where the value and exponent are both strictly positive. The range of this number should include $0$ and $[1, ~2^{\text{exponent}})$, skipping all values between 0 and 1.

My use case for this would be implementing a certain incremental game (Cookie Clicker) for an older game console (NES) where space might be an issue, the numbers could get large, and more complex operations are likely unnecessary.

I have looked at Steve Wozniak's and another person's implementation of floats in 6502 assembly and while they may work for me, I can't shake that this specific idea has to have been someone's PhD project, right?


问:
What is causing the incorrect solutions in my program while using scipy.integrate?
I'm trying to program the following algorithm

and I used scipy.integrate to evaluate the integrals. This algorithm should converge towards u and v which should look like this

Here's the program I used.

```<code>import numpy as np
import scipy.integrate as integrate
import scipy.special as special
import matplotlib
import matplotlib.pyplot as plt
import math

N = 20

u0 = np.ones(N)*0.08
v0 = -np.ones(N)*0.6

bta = 2.1
eps = 0.1
sgm = 1.5

def K(t, s, lbda, T):
    L1= np.exp(lbda*(T-s))*np.exp(lbda*t)/(1-np.exp(lbda*T))
    if t>=s:
      L2 = np.exp(lbda*(t-s))
    else :
      L2=0.0
    return L1+L2

def f(s,u,v):
  return np.sin(s)**2-bta*(u/(1+v))
def g(s,u,v):
  return -np.cos(s)**2+sgm*(u/(1+v))+eps/(1+v)

def sigma(s):
    return 1+np.sin(np.pi*s)**2

def integrand_u(s, t, lbda, T,u,v,lin):
    return K(t, s, lbda, T) * (sigma(s) + lbda*np.interp(s,lin,u) + f(s,np.interp(s,lin,u),np.interp(s,lin,v)))

def integrand_v(s, t, lbda, T,u,v,lin):
    return K(t, s, lbda, T) * (sigma(s) + lbda*np.interp(s,lin,v) + g(s,np.interp(s,lin,u),np.interp(s,lin,v)))

def integral_u(t, lbda, T,u,v,lin):
    result = integrate.quad(integrand_u, 0, T, args=(t, lbda, T, u, v, lin))[0]
    return result

def integral_v(t, lbda, T,u,v,lin):
    result = integrate.quad(integrand_v, 0, T, args=(t, lbda, T, u, v, lin))[0]
    return result

# Define the parameters
lbda = 10
T = np.pi
t = np.linspace(0,T,N)
err=1
epss=1e-5
n=0
while err>epss:
  u = [integral_u(a, lbda, T, u0, v0,t) for a in t]
  v = [integral_v(a, lbda, T, u0, v0,t) for a in t]
  print("u = ",u)
  print("v = ",v)
  err=np.sqrt(np.linalg.norm(np.array(u)-np.array(u0),2))+np.sqrt(np.linalg.norm(np.array(v)-np.array(v0),2))
  u0 = u
  v0 = v
  print('err = ',err)

plt.plot(t,u)
plt.plot(t,v)
plt.show()

</code>```

And while the error does seem to converge, I still get the same wrong solutions at whatever step I stop the loop.

Can anyone point out where I was wrong ? Thanks.


问:
Plot 2D piecewise function in MATLAB
I'm trying to use this method Plot 2D piecewise constant in matlab in a finite elements mesh for my data. However, I have problems with the values of n1, n2 and n3. I guess this is why it doesn't show a figure for what I have written.

What values would be appropriate for n1, n2, n3?

Title: Plot 2D piecewise function in MATLAB


问:
Existence theorem for Nonlinear BVP
I am recently trying to have a existence result for a nonlinear partial differential equation with boundary value.
What are the existence theorem for a nonlinear boundary value problem?</p>

Title: Existence theorem for Nonlinear BVP


问:
How can you plot a value in PYTEARCAT?
I've been working with Tensor calculus using the Python module PYTEARCAT. (See: https://arxiv.org/abs/2106.15016 and https://github.com/pytearcat/pytearcat.) It calculates quantities such as the Ricci curvature and Einstein tensor from a metric tensor, and this part works just fine for me. However, I can't seem to find a way to plot any quantities or components of any of the calculated tensors with respect to some coordinates (Such as plotting the 0,0 component of an Einstein tensor with respect to x and t). I'm wondering if there's some way to use Sympy (or some other module) to plot these quantities, as doing so would be very helpful. On a related note, I have also tried printing individual components of a tensor using indices (i.e. Ein[0][0]), but it just says that 'Tensor objects are non-iterable', which might be a problem as well since these are the objects I want to plot. So, if anyone has experience with PYTEARCAT--or is just good with code and doesn't mind looking at the links above--how would you go about plotting tensor components calculated by this module?</p>

Title: How can you plot a value in PYTEARCAT?


问:
Simpson rule and trapezoidal rule
Pls I am not that good in mathematics,  but I am solving an equation using trapezoidal rule and simpson rule , but according to simpson rule, can someone assist me on how to easily know the Yodds and Yeven .

Thank you.

Title: Simpson rule and trapezoidal rule


问:
fdtd ansys Lumerical vs fresnel solution
I am comparing reflectance intensity of silicon structure with a hole. I am getting oscillating solution for 0 roughness factor in equation based model. but spectrum shape of fdtd and equation based model solution starts matching if i increase the roughness factor. does this mean there is some default roughness factor in fdtd Lumerical solution?
thanks

enter image description here


问:
When is hyper-reduction needed in reduced order models?
Consider a high-dimensional ODE, which comes from semi-discretization

$$
\frac{d\mathbf{u}}{dt} = f(\mathbf{u}), \qquad \mathbf{u}\in\mathbb{R}^N \tag{1}
$$

If we want to build Reduced Order Models (ROM), we typically say

$$
\mathbf{u} = \mathbf{u}_{ref}(\mathbf{x}) + \mathbf{\Phi}{\mathbf{\hat{u}}}, \qquad {\mathbf{\hat{u}}}\in\mathbb{R}^n,   \mathbf{u}_{ref}\in\mathbb{R}^N,\mathbf{\Phi} \in \mathbb{R}^{N \times n}, n \ll N \tag{2}
$$

Using (2) in (1) we get

$$
\mathbf{\Phi}\frac{d\mathbf{\hat u}}{dt} = f( \mathbf{u}_{ref}+\mathbf{\Phi}\mathbf{\hat u}) \tag{3}
$$

Premultiplying by $\mathbf{\Phi}^T$ (Galerkin projection) and using $\mathbf{\Phi}^T\mathbf{\Phi}=1$ we get,

$$
\frac{d\mathbf{\hat u}}{dt} = \mathbf{\Phi}^Tf( \mathbf{u}_{ref}+\mathbf{\Phi}\mathbf{\hat u}) \tag{4}
$$

Since the RHS contains $\mathbf{\Phi}\mathbf{\hat u}$ and $\mathbf{\Phi}^T f$ the computation scales with $N$ and so is undesirable in a ROM (we want computations to scale with $n$). To reduce this cost, by making an approximate computation of the RHS which scales with $n$ instead of $N$, algorithms such as DEIM have been developed.

I understand the development so far.

However, consider the following concrete example, from Brunton and Kutz (Data-Driven Science and Engineering: Machine Learning, Dynamical Systems and Control) pg 425 , where $u=u(x,t)$ is a scalar variable

$$
\frac{du}{dt} = u^3 \tag{5}
$$

If we expand $u$ in terms of, say, two orthogonal POD modes, we get

$$
u(x,t)=a_1(t)\psi_1(x) + a_2(t)\psi_2(x) \tag{6}
$$

If we use (6) in (5) we get

$$
\psi_1\frac{da_1}{dt} + \psi_2\frac{da_2}{dt} = a_1\psi_1^3 + 3a_1^2a_2\psi_1^2\psi_2 + 3a_1a_2^2\psi_1\psi_2^2 + a_2^3\psi_2^3
$$

Multiplying by $\psi_1$ and integrating in space, and using orthogonality of $\psi_1$ and $\psi_2$ we get

$$
\frac{da_1}{dt} = a_1^3(\psi_1,\psi_1^3) + 3a_1^2a_2(\psi_1,\psi_1^2\psi_2) + 3a_1a_2^2(\psi_1,\psi_1\psi_2^2) + a_2^3(\psi_1,\psi_2^3) \tag{7}
$$

where the bracket $(\cdot,\cdot)$ denotes an inner-product in space. In the equation above, all the inner-products on the RHS can be pre-computed in the offline phase. Hence the cost of running the ROM (online phase) will not scale with $N$.

So, my question is why is hyper-reduction needed? Is it meant for non-linearities which cannot be explicitly written in polynomial form as above? Or is it simply meant to leverage the fact that we have a FOM (full order model) code which is capable of computing the RHS $f$ ? Is it because that the number of terms increases quickly in the RHS (multinomial expansion), rendering it impossible to compute even in the offline phase?


问:
Oscillation in non-linear porous flow solved by finite difference
I am trying to solve numerically the flow of a gas through a porous, spherically symmetric body. The non-dimensional equations read:
$$
\frac{\partial\rho}{\partial t}+\tau\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\rho v\right)=g_P\\
v=-\frac{\partial \rho}{\partial r}\\
$$
where $\rho$ is the gas density, $t$ time, $r$ radius (between $0$ and outer radius $R_p$), $v$ gas (radial) velocity, $g_P$ gas production rate, and $\tau$ a non-dimensional parameter accounting for permeability among other things (actually all quantities are non-dimensional). Essentially, the first equation is the conservation of mass and the second is Darcy flow. There is also a perfect gas law hidden in the second equation. Anyway, physics is not the point here. The boundary conditions:
$$
\frac{\partial\rho}{\partial r}(r=0)=0\\
\rho(r=R_p)=\rho_\infty
$$
In order to solve this non-linear system, I use a Picard iteration where the second equation is solved for $v$ using finite difference, the first equation is then solved for $\rho$ also by finite difference using the calculated $v$ and an explicit time-stepping (checking for CFL condition). The process is then iterated until $\rho$ has converged, and time is then incremented.

It seems pretty straightforward, yet when it comes to running the code, oscillation rapidly appear, propagating from the outer boundary, where a pressure (i.e. density in the non-dimensional form) gradient builds up, as you can see on the figure below:
<img src="https://i.stack.imgur.com/Q5kgT.png" alt="$rho$ profile at $t=1$ for $\tau=0.1$, $\rho_\infty=0.1$, a time step of 1e-5 and a radial resolution of 100 grid points" />
I had the same issue with an implicit time-stepping, so I thought I'd stick to the simple version. The amplitude of the oscillations is not sensitive to the timestep, and shows little sensitivity to the radial resolution. I am using first order central FD for both equations. I thought it might be an accuracy issue and went for second order central FD for the first equation, which to some extent helped, but not enough. Also, I am not sure about how to go about implementing a second order scheme given my set of boundary conditions; what I did is switch from a central scheme to a forward/backward one when getting to the boundaries, but is doesn't feel very neat.

I can elaborate on the model if needed, I just didn't think it was required, please let me know if I should. What I need is a diagnostic really, because so far I have been trying different ways without much success! The code is provided below.

Thank you for the help!

```<code>"""
Non-dimensionalization:
r[non-dim]   = r[m]/R_out
rho[non-dim] = rho[kg/m³]/(deltaP/R/Tref)
p[non-dim]   = p[Pa]/deltaP
v[non-dim]   = v[m/s]
gP[non-dim]  = gP[kg/m³/s]/g0
t[non-dim]   = t[s]/tau_g

Timescales:
tau_g = R/M*T*deltaP/g0
tau_p = viscosity/permeability*R_out^2/deltaP

The physical parameter is tau_g/tau_p,
which quantifies the gas production
time-scale compared to the pressure
gradient-driven outgassing time-scale.

"""

import numpy as np
import matplotlib.pyplot as plt
import time
import copy

# Parameter
tau_g_tau_p = 0.1
P_out       = 0.1 # aka rho_infinite

# All quantities are non-dimensional
N        = 100               # resolution
r        = np.linspace(0,1,N) # radii
dr       = 1/(N-1)            # radius step
rho      = np.ones(N)*P_out   # density
p        = np.ones(N)*P_out   # pressure
gP       = 1.                 # gas production rate
t_max    = 1                  # final time
dt       = 1e-5               # time-step
t        = 0.                 # initial time

while t<t_max:

    rho_prev   = rho
    res = 1
    tol = 1e-5
    while res > tol:

        ##### Picard iteration #####

        ##### Solve Darcy Equation #####

        # calculate p on the staggered grid
        p_stag = (p[:-1]+p[1:])/2 # pressure on the staggered grid

        # compute Darcy flux CFD from staggered p
        v       = np.zeros(N)                  # initialize Darcy flux vector
        v[0]    = 0                            # center BC
        v[1:-1] = -(p_stag[1:]-p_stag[:-1])/dr # central finite-difference
        v[-1]   = -(P_out-p_stag[-1])/(dr/2)       # outer BC

        # CFL check
        if max(v)*dt>dr:
            raise Warning('Courant criterion not respected')

        ##### Solve Mass Conservation Equation ######
        rho_new = rho
        rho_new[1:-1] += gP*dt - tau_g_tau_p*dt/r[1:-1]**2*(r[2:]**2*rho_prev[2:]*v[2:]-r[:-2]**2*rho_prev[:-2]*v[:-2]) # central finite-difference

        # boundary conditions
        rho_new[0]  = rho_new[1] # center BC
        rho_new[-1] = P_out      # outer BC

        ##### Perfect gas law #####
        p = rho_new # non-dimensonal rho and p are the same

        # calculate residuals
        res = np.linalg.norm(rho-rho_new)
        rho = rho_new

    t += dt # time increment

# Plot final state
plt.figure()
plt.plot(r,p,label=r'<span class="math-container">$dt$</span>='+str(dt)+',N='+str(N))
plt.xlabel('r [non-dim]')
plt.ylabel('rho [non-dim]')
</code>```


问:
Numerical Simulation of a Quadratic MIP with a highly rational term
I am interested in solving the following minimization problem:
$$
\begin{array}{cl}
\displaystyle\min_{x,y}&\displaystyle\frac{1}{K}\sum_{i=1}^{K}\left(\frac{x_{i}}{y_{i}}-\frac{X}{Y}\right)^{2} \\[0.4cm]
\mathrm{s.t.}&\sum_{i}x_{i}=X \\[0.2cm]
&\sum_{i}y_{i}=Y \\[0.2cm]
&x_{i}\leqslant y_{i},\quad i=1,2,\ldots,K \\[0.2cm]
&x_{i},y_{i}\in\mathbb{Z}^{+},\quad i=1,2,\ldots, K
\end{array}
\tag{1}
$$
with $X$ and $Y$ being fixed.

I attempted to simulate $(1)$ using $\textsf{cvxpy}$. However, it did not work as the quadratic term $\frac{x_{i}^{2}}{y_{i}^{2}}$ is neither convex, nor concave, and $\textsf{cvxpy}$ has no framework for such cases.

Question: How can we reformulate this problem so that it can be simulated by $\textsf{cvx/cvxpy}$, or what QMIP solver is able to handle an objective function of this form being highly rational and quadratic?

The solution to this minimization problem all boils down to $x_{i}=\lceil y_{i}\frac{X}{Y}\rceil$ or something close to that. For what its worth, I attempted to reformulate this problem by adding an envelop through $f_{i}:=\frac{x_{i}}{y_{i}}$. This leads to the following minimization problem:
$$
\begin{array}{cl}
\displaystyle\min_{f}&\displaystyle\frac{1}{K}\sum_{i=1}^{K}\left(f_{i}-\frac{X}{Y}\right)^{2} \\[0.4cm]
\mathrm{s.t.}&y_{i}f_{i}-x_{i}=0 \\[0.2cm]
&\sum_{i}x_{i}=X \\[0.2cm]
&\sum_{i}y_{i}=Y \\[0.2cm]
&f_{i}\leqslant 1,\quad i=1,2,\ldots,K \\[0.2cm]
&x_{i},y_{i}\in\mathbb{Z}^{+},\quad i=1,2,\ldots, K
\end{array}
\tag{2}
$$
However this formulation in $\textsf{cvxpy}$ was to no avail as there is no guided way to handle cases with bilinear constraints.


问:
How can the current-voltage relationship of a series of tunneling junctions be most easily computed?
I'm working on a research project with my professor where we're trying to figure out how to determine the necessary semiconductor composition for a solar cell to have a particular bandgap in order to maximize the power efficiency.  Specifically, we're attempting to simulate the effect of using alternating layers of AlGaAs and GaInAs, with varying percentages of Al and In.

The problem is that the Python library we're using (Solcore) only seems able to calculate the current for a given semiconductor structure if given the voltage or vice versa.  But what we're interested in is the general $IV$ relationship for a given structure (assume we know the doping concentrations, thicknesses of each side of the junction, etc.). How can this relationship, or at least a reasonable approximation of it, most easily be determined?

From speaking with the professor, as well as searching online, it seems the usual approach for this type of thing is to use the WKB approximation.  But he said the problem with doing it that way is that you need to know what the potential function is, and that it will be pretty complicated in this case because the extreme thinness of the barriers means there will be a lot wavefunction overlap.

Here's a plot of the wavefunctions  for one of the structures we're trying to simulate so you can see what I mean:

This is the Python code used to create it:

```<code>from solcore import si, material
from solcore.structure import Layer, Structure
import solcore.quantum_mechanics as QM

# First we create the materials we need
top = material("GaAs")(T=293, strained = False)
barrier = material("GaAsP")(T=293, P = 0.32, strained = False)
well = material("GaInAs")(T=293, In = 0.11, strained = False)
bottom = material("AlGaAs")(T=293, Al = 0.5, Strained = False)

# As well as some of the layers
top_layer = Layer(width=si("10nm"), material=top)
barrier_layer = Layer(width=si("2.5nm"), material=barrier)
well_layer = Layer(width=si("8.5nm"), material=well)
bottom_barrier_layer = Layer(width=si("50nm"), material=bottom)
bottom_layer = Layer(width=si("50nm"), material=top)

# Layered quantum wells
layers = [top_layer, barrier_layer] + 10 * [well_layer, barrier_layer] + [bottom_barrier_layer, bottom_layer]

test_structure_1 = Structure(layers, substrate=top)
output_1 = QM.schrodinger(test_structure_1, quasiconfined=0, graphtype='potentials', num_eigenvalues=20, show=True)
</code>```

Is there a reasonably efficient way to at least approximate the $IV$ relationship for this sort of structure, either by perhaps approximating the potential function with a relatively simple function or using a method that doesn't require using the potential function?


问:
How to show that the solution of the following quadratic programming is non-negative
I have the following quadratic problem:

$max$ $a^Tx+0.5x^TAx$

$s.t:    1^Tx=1$

in which $a=[a_1, a_2,...,a_n]$ is a non-negative vector, and $1^T=[1,1, ..., 1]$. The hessian matrix $A$ has the following structure:
\begin{bmatrix} -n\alpha\beta/a_1&\beta_{12}&\beta_{13}&...&\beta_{1n}\\\beta_{21}&-n\alpha\beta/a_2&\beta_{23}&...&\beta_{2n}\\\beta_{31}&\beta_{32}&-n\alpha\beta/a_3&...&\beta_{3n}\\...&...&...&...&...\\\beta_{n1}&\beta_{n2}&\beta_{n3}&...&-n\alpha\beta/a_n\end{bmatrix}
in which $\alpha=max(a_1,a_2,...,a_n)$ and $\beta=max(\beta_{ij}), i,j=1,...,n$.
Recall that $A$ is symmetric, and $\beta_{ij}=\beta_{ji}$ and that $\beta_{ij}\ge0$.

It follows that $A$ is negative definite; $A$ is an $M$ matrix and $A$ is strictly diagonally dominant matrix. Hence it is guaranteed that $A^{-1}\le0$. From $KKT$ conditions, the optimal solution is as follows:

$$x = (\frac{1^TA^{-1}a+1}{1^TA^{-1}1})A^{-1}1-A^{-1}a\\=\delta x_1+(1-\delta)x_2.$$

where $$x_{1}=\frac{A^{-1}a}{1^TA^{-1}a}$$ and $$x_{2}=\frac{A^{-1}1}{1^TA^{-1}1}$$

in which $\delta =-1^TA^{-1}a\ge0$. When $\delta\le1$ it follows that $x\ge0$. For $\delta>1$, I failed to prove that $x\ge0$. However, as I checked based on hundred thousands of random problem instances, it seems that $x\ge0$ for any value of $\delta$.
Please help me to prove or disprove that $x\ge0$ for any value of $\delta$.
I also could prove the followings that may be useful.

$\delta = \frac{a^T(x_1-x_2)}{-(x_1-x_2)^TA(x_1-x_2)}$

$\delta\le\sqrt{(1^TA^{-1}1)(a^TA^{-1}a)}$


问:
Verification of a Function Definition in Python
I want to write a function $f$ and it is defined as $f = - \nabla \cdot(|\nabla u|^{p-2} \nabla u) $ and I exact solution $u(x) = \tilde{u}(r) = 1 - \frac{p-1}{p-2} \left( s^{p/{p-1}} - (1-s)^{p/{p-1}} + 1 - \left(\frac{p}{p-1}\right)s \right)$ if $r<R$ and 0 if $r \geq R$. The following is the code block where I defined $f$ and $u$.

```<code>def grad_u(r, s, p):
    du_dr = np.where(r < 0.75,
                     ((p-1)/(p-2)) * ((p/(p-1))*s**(p/(p-1)-1)*(1/0.75) + (p/(p-1))*(1-s)**(p/(p-1)-1)*(1/0.75) - (p/(p-1))*(1/0.75)),
                     0)
    return du_dr

    def f_fun(input_xy, p):
        x = input_xy[:,0]
        y = input_xy[:,1]
        r = np.sqrt(np.maximum(x**2 + y**2, 1e-8))  # calculate the magnitude of the vector (x, y)
        s = r/0.75
        du_dr = grad_u(r, s, p)
        F_value = np.abs(du_dr)**(p-2) * du_dr
        # calculate the derivative of rF_value with respect to r and divide by r to get the divergence
        F = - np.divmod(np.abs(du_dr)**(p-2) * du_dr, np.maximum(r, 1e-8))[1]  # calculate f
        F[r>=0.75] = F[np.abs(r-0.75).argmin()]  # make f continuous
        return F
</code>```

and the exact solution I defined as

```<code>    R = 0.75
    p = 4 # Insert your value of p here
    s = r/R
    U_exact_all = np.zeros_like(r)
    cond1 = r<R
    cond2 = r>=R
    U_exact_all[cond1] = (1 - ((p-1)/(p-2)) * (s[cond1]**(p/(p-1)) - (1-s[cond1])**(p/(p-1)) + 1 - (p/(p-1)) * s[cond1]))
    #U_exact_all[cond1] = (1 - ((p-1)/(p-2)) * (s**(p/(p-1)) - (1-s)**(p/(p-1)) + 1 - (p/(p-1)) * s))
    U_exact_all[cond2] = 0

    U_exact_boundary = U_exact_all[:num_boundary]
    U_exact_uniform = U_exact_all[num_boundary:]
</code>```

While I am implementing it, the loss function is high. Could anyone please let me know the mistake in these two blocks?


问:
How can I incorporate angular momentum in this code?
I'm currently working on the 3-body problem, and I was writing a code to plot the trajectories of all 3 bodies while also manipulating the angular momentum of the system. I found a code online and tweaked it a little bit to match what I want. I was successful in plotting their trajectories; however, I couldn't measure their angular momentum.

Angular momentum is defined as the cross product of position $\vec{r}$ and momentum $\vec{p}$ (I have set the masses to be $1$, so momentum is the same as velocity in my case). I have tried to concatenate the arrays which represent the $x$ and $y$ components of the position and velocity of each body and used np.cross to calculate the cross product, but it didn't work. I'm not sure if my reasoning is wrong or I have just written the code wrong. Does anybody have a suggestion on how I could do this?

Here is the code:

```<code>import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

#Inital Conditions
y0 = [-1        ,#x1
  0             ,#y1
  1             ,#x2
  0             ,#y2
  0             ,#x3
  0             ,#y3
  0.306893      ,#vx1
  0.125507      ,#vy1
  0.306893      ,#vx2
  0.125507      ,#vy2
  -2*0.306893   ,#vx3
  -2*0.125507]   #vy3

# Time steps and the total number of steps
N = 10000
T = 0.001

#Definition of the Function
def ThreeBody(t,y):
 f = np.zeros(12)

 # The velocities of the three bodies
 f[0] = y[6]
 f[1] = y[7]
 f[2] = y[8]
 f[3] = y[9]
 f[4] = y[10]
 f[5] = y[11]

 # The x and y positions of each object respectively
 f[6] = -(y[0]-y[2])/(((y[0]-y[2])**2+(y[1]-y[3])**2)**(3/2)) \
       -(y[0]-y[4])/(((y[0]-y[4])**2+(y[1]-y[5])**2)**(3/2))

 f[7] = -(y[1]-y[3])/(((y[0]-y[2])**2+(y[1]-y[3])**2)**(3/2)) \
       -(y[1]-y[5])/(((y[0]-y[4])**2+(y[1]-y[5])**2)**(3/2))

 f[8] = -(y[2]-y[0])/(((y[2]-y[0])**2+(y[3]-y[1])**2)**(3/2)) \
       -(y[2]-y[4])/(((y[2]-y[4])**2+(y[3]-y[5])**2)**(3/2))

 f[9] = -(y[3]-y[1])/(((y[2]-y[0])**2+(y[3]-y[1])**2)**(3/2)) \
       -(y[3]-y[5])/(((y[2]-y[4])**2+(y[3]-y[5])**2)**(3/2))

 f[10]= -(y[4]-y[0])/(((y[4]-y[0])**2+(y[5]-y[1])**2)**(3/2)) \
       -(y[4]-y[2])/(((y[4]-y[2])**2+(y[5]-y[3])**2)**(3/2))

 f[11]= -(y[5]-y[1])/(((y[4]-y[0])**2+(y[5]-y[1])**2)**(3/2)) \
       -(y[5]-y[3])/(((y[4]-y[2])**2+(y[5]-y[3])**2)**(3/2))

 return f

 #Solving for the positions of all bodies
 t = np.linspace(0,N*T,N)
 solution = solve_ivp(ThreeBody,[0,800],y0,t_eval=t,rtol=1e-12)

 #Evolution in Position with respect to Time
 plt.plot(solution.y[0],solution.y[1],'-g') #Positions body 1
 plt.plot(solution.y[2],solution.y[3],'-r') #Positions body 2
 plt.plot(solution.y[4],solution.y[5],'-b') #Positions body 3
 plt.ylabel("Position(y)")
 plt.xlabel("Position(x)")
 plt.show()

 plt.plot(t,solution.y[0])
 plt.ylabel("Position (x)")
 plt.xlabel("Time")
 plt.show()
</code>```


问:
Convergence rate of Monte-Carlo variance estimate
What is the convergence rate for Monte-Carlo variance estimate for a random variable $X \in {L^q}(\Omega ,R),2 < q < 4$?</p>

ContentLicense: CC BY-SA 3.0


问:
Orthogonal vs general curvilinear coordinates
Solutions to PDEs over irregular domains can be computed using the finite difference method by the introduction of so called body fitted coordinate systems where the coordinate lines are aligned to the domain boundary. The benefit of the idea is that the finite difference method can be applied directly to the problem in a new coordinate system (with certain modifications regarding the computation of the derivatives).

As these modifications seem ridiculously complicated for practical consideration (http://hdl.handle.net/2060/19800017591 p. 18) I am looking for ways to make the problem easier by having a less generic coordinate representation.

My question is if using orthogonal coordinate systems has any advantages over a generic 2D curvilinear coordinate system (one to one mapping)? Orthogonal coordinate systems can be generated for a specific domain based on the Scwarz-Christoffel mapping (conformal mapping), what I would be interested in is if this approach would have any practical benefits, i.e. if orthogonality can somehow be exploited for this problem?


问:
Find maximum distance between elements given constraints on some
I have a list of numbered elements 1 to N that fit into positions on a number line starting with 1. I also have constraints for these elements:

The element 1 is in position 1, and element N must be at a position >= the position of element N-1. (i.e. element 2 could be at position 1, element 3 at position 7, and element 4 at position 8 (but not position 5))

Some elements must be within a certain distance from each other on the line.

Some elements must be at least a certain distance from other on the line.

My objective is to return an integer that represents the maximum span between element 1 and element N. If no lineup is possible, return -1, and if the elements can be any distance apart, return -2.

I am given:

The number of elements

A withinArray[][] where withinArray[x][y] = the distance elements x and y must be within on the line. Any zero values represent no constraints.

An atLeastArray[][] where atLeastArray[x][y] = the distance elements x and y must be apart on the line. Any zero values represent no constraints.

An example input would be: 4 elements, withinArray1[3] = 10, withinArray[2][4] = 20, and atLeastArray[2][3] = 3. (all other array values are zero).

The return value for this input would be 27. (element 1 at position 1, element 2 at position 8, element 3 at position 11, and element 4 at position 28)

The problem was first posted here by someone else. I'd like to figure out an elegant solution to it programmatically. Though I've been working on it for a whole day, I still have no luck coming up with a good solution. I feel that I need to use dynamic programming techniques, but have a hard time finding a good substructure.

I am not able to locate the same example on the web. Any pointer to such materials online would be appreciated. It's even better if you are an expert for such kind of question and can outline the solution in detail. Executable code is a plus.

EDIT: I know how to formulate this as a linear programming problem and I know how to solve it using LP solver. I just prefer an end-to-end solution with input and output as specified. General answers might not help me much.


问:
How to accurately decompose positive semidefinite matrix and use the lower triangular part in linear equations
I have $n$  arbitrary $p\times 1$ vectors $x_i$, and $p\times k$ matrices $A_i$, and $n$ $p \times p$ positive semidefinite matrices $S_i$, where some (often most) of the $S_i$'s are same (for example only two different $S$ matrices, one positive definite which applies to $i=1,\ldots n-1$ and semidefinite $S$ for $i=n$). I want to make following linear transformations for all $x_i$ and $A_i$:

$x_i^* = L_i^{-1}x_i$ and
$A_i^* = L_i^{-1}A_i$, where

$L_i$ is a lower triangular matrix so that $L_i L'_i=S_i$ (or even better $L_i D_i L'_i=S_i$ where $D_i$ is diagonal matrix.)

$n$ can be range from few to thousands or even more, $p$ is usually less than 10 and $k$ is usually less than 100. $S$ can contain rows and columns with zeros, and it is originally formed as $BB'$ where $B$ is a lower tridiagonal matrix with nonnegative diagonal elements.These $B$'s are not available though.

I'm wondering what would be the optimal way of accomplishing this in terms of speed and accuracy, with more weight on accuracy?

Currently I'm using my self written LDL decomposition function for different $S$, invert $L$ and compute the transformations, as I have beem dealing with cases with small number of different $S$ and large $n$. If I have understood correctly, it would be wiser to just solve the linear equations without explicitly inverting the matrix in terms of accuracy, but for speed it seems to be better option?

I'm using Fortran, and I would rather use some LAPACK function for the factorization instead of my own LDL decomposition which I cannot guarantee to be stable etc., but I'm not sure which would be the good way.

edit: This may be better suited for stackoverflow, so this can be removed by whoever has the power to do so.


问:
Conservative field mapping between two topologically disconnected surface meshes
Some background: the Front-Tracking method uses a triangular surface mesh to describe the boundary between two immiscible fluids. To deal with the breakup and coalescence of the fluid interface, direct topological operations are executed on the surface mesh. Since in this case, there exists a history of topological operations between the old and the new mesh, it is possible to conservatively map a field stored on the surface mesh before and after topological operations. This is used e.g. to solve transport equations for a surfactant (a field that describes surface-active agents whose concentration on the fluid interface changes the surface tension).

What I would like to know if there is an established algorithm that can logically connect two completely topologically disconnected surface meshes that could serve as a mapping basis?

E.g. something like the nearest element search algorithm: for each triangle of surface mesh A, find the nearest triangle of surface mesh B.


问:
Enforcing continuity conditions in pseudospectral domain decomposition methods for time dependent PDEs
I have a partial differential equation of the form
$$ \frac{d}{dt}f(x,t) = \Theta(x) f(x,t) \qquad \Theta(x) \sim \left[\frac{d^2}{dx^2} + k^2(x)\right] $$
subject to $f(x,t=0) = f_0(x)$, and $f(x=0,t) = f(x=1,t) = 0$. The function $k(x)$ is piecewise constant over $x\in[0,1]$. Call the $x$ intervals where $k$ is constant $I_j$. I would like to locally represent the solution $f$ as a sum of (say) Chebyshev polynomials within each interval $I_j$ and enforce continuity conditions at the boundaries of the intervals (continuity of $f$ and $\frac{\partial f}{\partial x}$).

Now, since the RHS operator is time independent, I could just solve this using an exponential: $f=e^{\Theta t}f_0$. My question is how I would compute the matrix elements of $\Theta$ to take into account the subdomain continuity conditions.

Alternatively, I could use a time marching method and make this a true domain decomposition method, but I still don't know how I would enforce the continuity conditions. I don't want to use a penalty method to enforce the constraints only approximately; I want them enforced as exactly as possible, and I really don't care about parallelism (so perhaps calling this DD is misleading). The literature on this probably exists, but I simply have no idea what I'm looking for.

I only give a rough description of $\Theta$, but it comes from the Helmholtz equation so it's hyperbolic, and a $2\times2$ block operator, with $f$ actually being a 2-vector of a function with its dual. This isn't terribly relevant to the question, since I just want to know how to enforce internal consistency.


问:
PCA performed on a configuration with scaled axes
Suppose a configuration $X\in\mathbb{R}^{n\times 2}$ is output of PCA on some high-dimensional data $Y\in\mathbb{R}^{n\times h}$. Note that this PCA is performed by $$X=Y\cdot U,$$ where columns of $U$ are eigenvectors of $Y^TY$ corresponding to its 2 dominant eigenvalues. Since $X$ is a result of PCA, then $X^TX=diag(a, b)$, for some $a, b\in\mathbb{R}$. If one now performs PCA on $X$, then the matrix $V$ containing eigenvectors of $X^TX$ would be columns of an identity matrix, and the result would be intact.

Suppose that a configuration $X$ is normalized (axes of unit length) by $\hat{X}=X\cdot diag^{-1/2}(a, b)$. If one then performs PCA on $\hat{X}$, then $\hat{X}^T\hat{X}$ would be an identity matrix. In other words, by normalizing columns of $X$ one cannot affect the PCA result except for the aspect ratio change (stretching of axes). If fact, any arbitrary stretching of axes of $X$ would have no effect since also then $\hat{X}^T\hat{X}$ would be a diagonal matrix. Am I missing something?


问:
Hamiltonian Matrix Size in Schrodinger Equation
I'm attempting to solve the particle-in-a-box problem using Scipy (with the help of http://www.physics.buffalo.edu/phy410-505/2011/topic4/app2/index.html). At first, I used a 16x16 matrix to model the Hamiltonian, like the link, and my results corresponded to theirs. However, when I used a larger matrix (50x50), I found many extraneous eigenvalues in my results due to the larger matrix size.

Why does a 16x16 matrix produce exactly the correct eigenstates while a larger one produces extraneous ones? When using a larger matrix (which I thought would increase the accuracy of the simulation due to fewer basis elements being omitted), how can I tell which elements correspond to actual eigenstates and which are extraneous?

My code is below:

```<code>from scipy import linalg, mat, matrix

def Sfun(m,n):
    if (m+n)%2==0:
        v1 = 2/(m+n+5)
        v2 = 4/(m+n+3)
        v3 = 2/(m+n+1)
        return v1 - v2 + v3
    else:
        return 0

def Hfun(m,n):
    if (m+n)%2==0:
        return -8*(1-m-n-2*m*n)/((m+n+3)*(m+n+1)*(m+n-1))
    else:
        return 0

Slist = []
Hlist = []

for m in range(0,16):
    tlist = []
    for n in range(0,16):
        tlist.append(Sfun(m,n))
    Slist.append(tlist)

for m in range(0,16):
    tlist = []
    for n in range(0,16):
        tlist.append(Hfun(m,n))
    Hlist.append(tlist)

Smat = matrix(Slist)
Hmat = matrix(Hlist)

vals,vecs = linalg.eig(Hmat, Smat)

for i in range(0,16):
    print('Vector: ', end="")
    print(vecs[i], end="")
    print("                     Value: ", end="")
    print(vals[i])
</code>```


问:
computation complexity of OLS in estimating a VAR model
Could someone tell me the computation cost of using Ordinary least squares in estimating a Vector autoregression model?
I am thinking the cost is O(n) where n is the number of the training instance.
But I am not sure about this. Thank you.</p>


问:
System of non-linear ODEs and estimating unspecified initial conditions on Maple 12
I have the following 1st order equations and need to solve them using Maple 12.
There are unspecified initial conditions and can only be estimated through the Newton raphson method. My problem is how do I implement it so that the equations are numerically solved using rk4-order. Tried out with BVPsolve but doesn't work.
this is in line with works by Makinde.
"On MHD boundary-layer flow and mass transfer past a vertical plate in a porous medium with
constant heat flux."</p>

```<code>> k1 := diff(X[1](t), t) = X[2](t);
> k2 := diff(X[2](t), t) = M*(X[1](t)-1)-(2*(eta+b))*X[2](t);
> k3 := diff(X[3](t), t) = X[4](t);
> k4 := diff(X[4](t), t) = (2*Sc*Du*(eta+b)*X[6](t)-Du*lambda*X[5](t)-2*Pr*(eta+b)*X[4](t)-   Pr*Ec*X[2](t)^2-Pr*Ec*M*(X[1](t)-1)^2)/(1-Du*Sr);
> k5 := diff(X[5](t), t) = X[6](t);
> k6 := diff(X[6](t), t) = (lambda*X[3](t)+2*Pr*Sr*(eta+b)*X[4](t)+Pr*Sr*Ec*X[3](t)^2+Pr*Sr*Ec*M*(X[1](t)-1)^2-2*Sc*(eta+b)*X[6](t))/(1-Du*Sr);
> ICS := X[1](0) = 0, X[2](0) = S[1], X[3](0) = 1, X[4](0) = S[2], X[5](0) = 1, X[6](0) = S[3];
</code>```

OwnerDisplayName: Saronie


问:
exact area resampling
I do image processing, and right now I need to resample some images taken from slightly different perspectives so I can match up features. The pixel intensities have scientific significance, so I want to do something that will preserve local flux and not "miss" anything. What I want to do is exact area resampling, as discussed here:

http://www.cl.cam.ac.uk/techreports/UCAM-CL-TR-261.pdf

Except the coordinate change I'm doing is not just a translation; at the very least it involves translation, rotation, and perspective change.

I can't find anything on file exchange or by googling. Any ideas?

Title: exact area resampling


问:
Solving PDE or eigenvalue problems without FEM
Do you know any methods/solvers for PDE or eigenvalue problems like

$\begin{cases} \Delta u= 0\ (\text{ or }\lambda u) & \text{ in }\Omega \\ u =0 & \text{ on }\partial \Omega \end{cases}$ (Dirichlet)

$\begin{cases} \Delta u= 0\ (\text{ or }\lambda u) & \text{ in }\Omega \\ \frac{\partial u}{\partial n} =0 & \text{ on }\partial \Omega \end{cases}$
(Neumann)

$\begin{cases} \Delta u= 0\ (\text{ or }\lambda u) & \text{ in }\Omega \\ \frac{\partial u}{\partial n}+\beta u =0 & \text{ on }\partial \Omega \end{cases}$
(Robin)

in 2D and 3D that do not involve the Finite Element Method?

In 2D there is at least an option (for Dirichlet and Neumann): mpspack which uses a basis for the space of harmonic functions and works only with coefficients in this basis. In 3D I do not know anything that doesn't use Finite Element Method.

I am interested in a different method than FEM because in shape optimization domains change at every iteration and so does the mesh, the discretization, etc.


问:
Solving diffusion PDE using finite differences
I need some hints on how to solve this diffusion equation ($\alpha, k_1,k_2$ and $k_3$ are constants):

$$ {\partial P \over \partial y} + k_1 {\partial P \over \partial t} + \alpha P  = {1 \over k_2} {1 \over x}{\partial  \over \partial x} \left(x k_3 {\partial P \over \partial x}\right) $$

With boundary conditions:

$P(x=1,y,t)=0$

$k_3 {\partial P \over \partial x}|_{x=0}=0$

$P(x,0,t)=p(x)$

I've been told that DuFort-Frankel is a proper method for this example, so I tried (not sure if it is correct):

$${P_{i,j+1}^n-P_{i,j-1}^n \over 2 \Delta y}+k_1{P_{i,j}^{n+1}-P_{i,j}^{n-1} \over 2 \Delta t}+\alpha P_{i,j}^n = {k_3 \over k_2x}\left({P_{i+1,j}^n-P_{i-1,j}^n \over 2 \Delta x}+x{P_{i+1,j}^n-(P_{i,j}^{n+1}+P_{i,j}^{n-1})+P_{i-1,j}^n \over \Delta x^2}\right)$$

But I don't know how tie this with boundary conditions so I can solve it using recursive functions... It is supposed to be pretty easy, am I missing something? DuFort-Frankel is not necessary, if You know how to solve it using Taylor, Leapfrog, Richardson or any other method, I would be very grateful for any hints...

OwnerDisplayName: apm


问:
Numerics for heat equation
I want to simulate on my computer the solution of the heat equation in 3 space dimensions with Cauchy initial data, that is
$$\partial_t u=Tr[A(x)\cdot \Delta u], u(0,x)=u_0(x) $$

where $u_0\in C(\mathbb{R}^3,\mathbb{R})$.

Even if $A$ is constant I'm not sure what's numerically the best way to do this but I'm sure this must be in standard toolboxes of Matlab, Scilab etc. however I have quite some trouble with the documentation.

Could someone give me pointers in the right direction?
More general, what is a good book or lecture notes that deal with numerics for parabolic PDEs (with multidimensions in space)?

OwnerDisplayName: john

Title: Numerics for heat equation


问:
Parallel algorithm to use in place of PORTA?
We currently use PORTA software to find the list of facet-defining inequalities (FDI) for polytopes that we work with.  For certain polytopes, PORTA works fine.  But because it is a serial algorithm (uses Fourier-Motzkin) we sometimes encounter cases where the analysis simply takes too long and we have to force the program to stop.  We are wondering if there exists a parallel algorithm that would allow us to complete these computations faster? Ideally we would prefer to have exact solutions like PORTA gives; but in lieu of such software programs we would be open to trying a method that may not yield exact results.  (And if it matters:  our input coordinates are almost always discrete, either 0 or 1.)</p>

OwnerDisplayName: cez


问:
An asynchronous version of the Covariance Matrix Adaption (CMA)?
Does anyone know of a variant of the CMA algorithm that is suitable for an asynchronous parallel implementation?

The conventional version of the algorithm allows one to evaluate the fitness function at different points in parallel, but one has to wait until all computations are finished before starting the next iteration. If the computational cost of the fitness function is heterogeneous, the parallel efficiency quickly degrades. Ideally an asynchronous version should be such that parameters are not disadvantaged when the fitness function takes longer to compute.


问:
Lax equivalence theorem for integro-differential equation
Can the Lax equivalence theorem (http://en.wikipedia.org/wiki/Lax_equivalence_theorem) be applied to the discretization of integro-differential equations, or does a similar theorem exist for them?</p>


问:
A cost function minimisation genetic algorithm which includes a measure of search completeness?
I am looking at ways of finding the global best-fitting model to a set of data. The model is highly non-linear and has anywhere between 5 and 8 parameters.

I have used particle swarm optimisation with some success, as well as Metropolis-Hastings Markov Chain Monte Carlo.

At the end of any given run, there is still doubt as to whether the parameter space has been completely  - or at least sufficiently - searched.

Are there any algorithms for model fitting which includes some measure of how thoroughly the parameter space has been searched?

For example, considering PSO, is some way of estimating how "sensitive" any given particle is to its present neighbourhood? Let us consider a PSO particle, passing a local cost minimum in the N-dimensional manifold at a certain velocity and a certain distance from the minimum. Can we quantify the probability that this particle would be attracted to and "investigate" the minimum?


问:
Propogated Error in Mesh Interpolation
I am working with a code that solves diffusion/reaction equations on a 2D unstructured mesh.  Due to the stiffness of some of the processes, I start with time steps near 1e-13, and end with a final time in the millisecond range with time steps near 1e-8.

The issue appears in that the problem is also, what I would like to call geometrically stiff.  During the first ~250 ns of simulation I require a very fine mesh to resolve a geometric characteristic of the forcing function.  After this time, however, the forcing of the system is approximately zero, and I no longer need the spatial resolution.

The fine mesh takes about 2 weeks to run through that first 250 ns.  The coarse mesh takes 2-4 weeks to run the entire simulation out to 2.5 ms (an order of magnitude speed up).  Running an adaptive mesh is not possible without a serious rewrite of the code (currently beyond my skillset as well), however there are utilities for me to restart a simulation.  What I am considering doing is running the fine mesh out until it is no longer needed, taking the restart dump file, interpolating it onto the coarse mesh, then restarting to let it run out to the final time we want.

I can estimate the error in my interpolation fairly easily depending on my interpolation method, however I am unsure of how to estimate the error that is propogated once I begin time-stepping again.  I am particularly worried about the smoothing on variables such at temperature will have, since many of my process rates are exponentially dependent on it.  My question is, how do I estimate how some fixed error propogates through my stiff system without running multiple cases (computationally infeasible at the present time).

Title: Propogated Error in Mesh Interpolation


问:
Modified heat eq. in 3D
please how can I write gradient numerically?
I have an eq.
$\dfrac{\partial X}{\partial t}=\nabla^2 X + u\nabla X$
I wonder how write this in 3D.
Many thanks for any idea... I can find how to solve the "pure" heat eq. but so I don't know.. </p>

Title: Modified heat eq. in 3D


问:
FFT - function only in sine series? Can be done with MKL / Lapack?
please can I ask, how one can make from function sine series (Fourier transform) with MKL? I can do "normal" exponential FFT with MKL (Lapack of course), how can I say that I want only sine series?
Many thanks</p>


问:
How to compute the sum of a power series in a more robust way?
In order to compute the sum of a power series, we can use for loop, while loop or the analytic formula. I am wondering what is difference between those algorithms and how to improve the robustness of this question?</p>


问:
Adaptive mesh data structure for Fast Marching Method to overcome RAM limit
On an uniform mesh of positions in space $\ (x_i,y_j,z_k)$:

$$\ x_i = x_0 + i\Delta x,\quad i=0,\ldots,n_x$$
$$\ y_j = y_0 + j\Delta y,\quad j=0,\ldots,n_y$$
$$\ z_k = z_0 + k\Delta z,\quad k=0,\ldots,n_z$$

I have a velocity function: $V(x_i,y_j,z_k)$.

Given a start position $(x_s,y_s,z_s)$ I can calculate the time it takes to reach all points in space: $T(x_i,y_j,z_k)$ by the Fast Marching Method (FMM).

However, I would like to do this on a mesh larger than can fit in RAM,
e.g. $n_x=n_y=n_z=2000$.

FMM needs access to ALL velocity data. It is NOT possible to run it on only a part of the mesh and then later combine results somehow.

The data structure storing the velocity data can, however, be sparse. As long as the FMM can ask for the value at any $(x_i,y_j,z_k)$ it does not care how it is stored internally.

I imagine setting blocks of values at a time:

```<code>struct SparseArray3D {
public:
    SparseArray3D(int nx, int ny, int nz);
    void setBlock(int i0, int j0, int k0, float block[64][64][64]); // i0 <= i < i0 + 64, j0 <= j < j0 + 64, k0 <= k < k0 + 64
    float& operator()(int i, int j, int k); // needed by FMM
};
</code>```

Do you have any tips on which data structure I could use?

(As you probably already have guessed) I am programming in C/C++.


问:
What are Implications of Commutative Diagrams?
This question may be too broad. But I really want to know some concrete explanations.  I often find various commutations appear here and there, which concerns the application order of two operators.  For examples, between exterior calculus of the differential forms and discretization,  optimization-then-discretization or discretization-then-optimization, and inner-outer iterations.  Do you have some understanding of concrete examples where the commutation plays an important role? </p>


问:
Discontinuous Galerkin for flow through porous media
I am struggling with DG methods for 2 phase flow through porous media. I managed to get the global pressure, total flux equations to work with an unconditionally stable mixed FE DG formulation as described here:
http://math.aalto.fi/reports/a531.ps

but the saturation equation is giving me problems. Basically I tried LDG and IP and both gave pretty weird unphysical results.

Does anyone have any good references for DG with flow through porous media? I found at least one good compositional flow paper, but not much about blackoil models.

Thanks in advance.

PS: I probably should add that I am writing the code with FEniCS and CBCPDESys. So far it has saved me years in development time.


问:
Create 2D plot from two data sets in ParaView
I'm currently running a CFD simulation for an internal flow and need to compare my computed results to experimental data. The experimental data is saved as a comma separated file with columns for x,y,z coordinates and pressure. The computed data is saved in Ensight format on an unstructured 3D mesh.

As an example, I have the following experimental data:

-0.08573,0.00000,0.09154,83600.0

-0.08573,0.00000,0.08154,83610.0

etc.

I would like to plot the pressure along the y-axis at x=-0.08573 as a function of z, and overlay crosshairs to show the experimental data. I can create the two plots separately, using the "Plot Over Line" filter for the CFD data, and the "Plot Data" filter for the experimental data.

I have tried using the "Group Datasets" filter to combine the two sources, but when I try to plot I run into problems.

If I use "Plot Over Line", the experimental data doesn't show in the plot. I think this is due to rounding in the coordinates for the line vs. the data.

If I use "Plot Data", I can get both to show up, but only if I set the X axis data to "Use Array Index From Y Axis Data". When I do this, the experimental data aren't plotted over the correct part of the CFD data. Switching to "Use Data Array" and Points(2) for the z-coordinates, I get an error "No X column is set (index 0)."

This seems like it should be a simple thing to accomplish, and a fairly common task, but I've had no luck finding a tutorial or explanation online. I know that I can just export the data and use gnuplot or xmgrace, but I would rather do all the processing in one program if possible.


问:
Cusp Library performance worse than PETSC (GMRES 200 iterations) Why?
I wanted to compare the speeds of the GMRES implementations in the CUSP and the PETSc libraries.

The matrix (A) used for testing was a 3d Laplacian matrix obtained by using the 7 point stencil on a 20x20x20 grid => the matrix is of size 8000 x 8000. This is a banded matrix with 7 diagonals.

The RHS (b) was obtained by multiplying A by a vector of 1's. I wanted to solve Ax=b
by GMRES

For both PETSC and CUSP  I timed only the time taken for 200 iterations of the GMRES solver
(ie time taken for reading in the matrix is not counted)

On Petsc the timings were done with PetscGetTime() and on Cusp with Cuda events.

%--------------------------

Timings

Petsc (1 cpu) took 0.0558379 sec  to do 200 iterations of GMRES.

Cusp          took 0.1183 sec  to do the same.

%----------------------------

As you can see Petsc is 2.14 times faster than the Cusp library even on a single processor.
The speed gap gets larger (obviously) on using more than 1 processor.

This experiment was done on a  GTX 570 card on an Ubuntu 10.10 running CUDA 4.0 having cusp v 0.3.1 and thrust 1.6.0

I have read this paper by one of the authors of CUSP. However, they have
not compared the performance of CUSP with Petsc in it.

The paper also says that the data-structure used for storing the sparse matrix matters.
Accordingly they have provided several formats like DIA, ELL, HYB, COO etc.
but even after I tried them all, the Cusp GMRES performance does not change and still takes 0.1183 seconds for the 200 iterations.

Here is the (quite-short) GPU code on pastebin I used while using CUSP. I am not sure if I am using the Cusp
library in an optimal manner.

Although it is possible PETSc is still genuinely better than CUSP I want to know if significant speedup is possible with CUSP if used properly.

If more information is needed to answer this question please let me know.
Thank you.


问:
Algorithm to distribute symbols uniformly in a square grid?
Given a square grid of size $n\times n$ and $m$ symbols (say for example, alphabets A, B, C...), having $N(i)$ number of $i$th symbol; $\sum_{i=1}^{m}N(i) = n\times n$. Is there any computationally tractable algorithm that can be used to distribute the symbols uniformly inside the grid?

Uniform distribution here could be "shortest-distance-pairs have similar values". An example with 3 symbols A,B,C with equal counts of 3 in a $3\times3$ grid, I would expect a distribution of something like $ \left[ \begin{array}{ccc}
A & C & B \\
C & B & A \\
B & A & C \end{array} \right]$. There could be many possible solutions, but the algorithm is expected to give one of them. And one more thing I am looking for is, the algorithm place symbol(s) with count value of 1 at (or near) the center of the grid, rather than at the corners.


问:
Robust Counterpart of an uncertain LP
Consider the following robust optimization problem:

min c'x

s.t.: $Ax\geq b    \;\;\forall (A,b)\in \mathcal{U}$.

Why can the robust counterpart of the problem be written in this form?
$min_x{\{ max_{(A,b)} Ax\geq b\}}$?

How does it work principally? Does it mean that I have to, for a fixed $x$ value, maximize the constraints $Ax\geq b$ for all realization of the uncertain data, and then minimize the value $x$?

Please make this point clear to me with a small example.

Title: Robust Counterpart of an uncertain LP


问:
Dissipation and symplectic manifolds
I'm working on an API for simulation of port-Hamiltonian systems. As far as I understand it, a Hamiltonian system is symplectic if it is power conserving, and so including resistive elements would thus destroy the symplectic nature. My question is, if I add an extra state variable, the dissipated heat, along with its adjoint, does this restore its symplecticity?

Further, if I didn't include the extra state variable, would I still benefit from using symplectic integrators? And do systems lying on a contact manifold still benefit from symplectic integrators?

Title: Dissipation and symplectic manifolds


问:
What libraries are available for solving problems at the college physics level?
Are there any libraries available for solving problems at the level of college physics?

In particular, I'm interested in libraries for general purpose programming languages that can calculate results both numerically and symbolically.

Libraries which could be used for implementing something like this are SymPy (Python), GiNaC (C++), or SymbolicC++ (C++).

I'm asking this question because I've started cooking up a C# library for symbolic computation and have used it to solve a few projectile motion problems, but would like to review other work out there.

Here's an example of a program which solves such a problem.

Here's what the program outputs when run:

Any pointers or suggestions welcome. Thanks!


问:
Memory allocation error with GSL ODE solver applied to system of 4 ODEs
I am trying to solve a (large) system of ODEs with GSL solvers. When I use driver method I get an error message of could not allocate space for gsl_interp_accel, when I define control, error and stepper manually, I get bad_alloc exception, which is caused, as far as I understand it, by the same thing that causes could not allocate space for gsl_interp_accel in the other case - the lack of memory.

I have consulted with other bad_alloc queries, such as this one, but I haven't found anything useful for my particular case there. Also, I have tried other ODE solvers, but they also end up with memory errors. I have also checked my program with valgrind to make sure there are no memory errors/leaks anywhere else apart from the solver.

Any solver has "limits of integration", and in my case program works fine for about 10% of the upper limit (which is large comparing to lower limit - I am pretty sure this is the source of errors I get - but I do need to integrate between these particular limits) and then terminates with one of the exceptions I've quoted above. I have tried various (fixed/adaptive) step size, but never reached more than 10% of what I want.

The piece of code that gives an exception is:

```<code>gsl_ode_struct inputstruct;  // Struct that contains parameters for ODEs
gsl_odeiv2_system sys = {func, NULL, 2*G.getsize(), &inputstruct};
const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk8pd;
gsl_odeiv2_step * stepper = gsl_odeiv2_step_alloc (T, size_of_the_system);
gsl_odeiv2_control * control = gsl_odeiv2_control_standard_new (1.e-6, 0.0, 1., 0.);
gsl_odeiv2_evolve * error = gsl_odeiv2_evolve_alloc (size_of_the_system);
double hparam = 1.e-6; // Initial step size
double t = xspan[0]; // Initial time
while(t < final_time){
    // Here the excpection comes
    int status = gsl_odeiv2_evolve_apply (error, control, stepper, &sys, &t, final_time, &hparam, array);
    if(status != GSL_SUCCESS)
        break;
    // Do some stuff that includes reading all the intermediate results to a container as I need them later.
    }
    gsl_odeiv2_evolve_free (error);
    gsl_odeiv2_control_free (control);
    gsl_odeiv2_step_free (stepper);
</code>```

So if I change final_time to final_time/10 code executes, but the result then does not make any sense. Even when nothing is done after solver, exception is still thrown, could not allocate space for gsl_interp_accel, though.

I have tried to split the loop on several (many) loops with erasing memory in between, but this didn't help much.

In case this is important, I use Ubuntu 12.10, compiled with GNU compiler and Intel C++ Composer. Also tested on Mac (don't know which version of OS) with the same result.

The question is: is there any way to "cheat" on the solver and make the program work properly?

P.S.: ODEint solver, that has smarter way of getting intermediate results, also throws an exception.


问:
Can I use RANS to see the effect of mixed convection?
my question is: can I use a RANS simulation to see the effect of mixed/natural convection?

Actually I have also a second question: I would like to do this in Comsol multiphysics, but it seems that it has no model that takes into account turbulence (like k epsilon, Reinolds stress, LES, ecc..).
Ps: this could be because I use a pretty old version of the software, but anyway I haven't found documentation about this topic so far.

Thank you for your answers!


问:
Pixel-To-Angle Transformation in Camera Image
I'm trying to localize points I see in a camera image in terms of azimuth and elevation and match points between shots.
Individual shots should differ only in rotation around the camera's center (roughly).
It's not working (as I'll explain below), and I'm looking for things I might be doing wrong.

I've tried different things, but this is what I'm doing right now:
Let's say I have tilted my camera (around the $x$-axis) by some angle $\alpha$ and rotated it (around the $y$-axis) by $\beta$.
Assuming my lens gives me a perspective projection with a focal length of $d$, I can represent pixel $(i,j)$ as a vector
$$ \vec{v}=d\cdot\vec{z}' + i\cdot\vec{x}' + j\cdot\vec{y}', $$
where $\vec{x}', \vec{y}', \vec{z}'$ are the unit vectors of the standard basis rotated by $\alpha$ and $\beta$ around x-axis and y-axis, resp.
$\vec{v}$ then represents the point in 3D space at which the photoreceptor cell in the CMOS was when it perceived the point in real space.

Say $\vec{v}=(v_1,v_2,v_3)$.
Then I should be able to calculate azimuth $\phi$ and elevation $\psi$ wrt. an unrotated take as:
$$\phi = \mathrm{arctan}(\frac{v_1}{v_3})\quad\text{and}\quad\psi = \mathrm{arcsin}(\frac{v_2}{|\vec{v}|}).$$

This doesn't work, though.
That is, without tilting the camera it more or less does.
But say I have a shot where both $\alpha$ and $\beta$ are not zero.
Then, the points which appear in that shot and in others for which $\alpha$ or $\beta$ are zero don't come out with the same azimuth and elevation.
It seems that they generally cluster too much around the center of that shot.

One of the problems is that there are a lot of things that could get in my way here, like any combination of:

a distorted camera image

wrong focal length

wrong rotation angles

rotation around a point too far from the focal point.

I'm looking for more potential sources of error and (in)validation for the approach outlined above.

Above is a visualization of the outcome of the procedure so far: the white horizontal and vertical axes are azimuth and elevation, resp.
The distorted rectangles are boundaries of individual shots.
It's hard to explain what exactly the points are, but in essence points with the same color are identical in the real world and should, therefore, appear at the same absolute azimuth and elevation.
These points form a distorted shape, not a rectangle.
You can see that at each level of elevation, the points from each camera image appear more and more contracted around the camera image's center.


问:
features recognition & reconstruction of 3d mesh delaunay matlab
I managed to display the coordinates of x,y and z into a 3D mesh by using delaunay function. The coordinates are in .obj format actually and i have read it into matrix form.

Now, i would like to extract certain features / regions from the 3D mesh, reconstruct a new mesh and get the coordinates of it (I would like to extract specific ROI or volume of interest (VOI) form the mesh).

For example:

A 3D facial image mapped into triangular mesh, and i want to extract the nose region. the same method cam be used as template for other data.

Is it possible? Thank you in advance.


问:
How to do FEM in sector elements?
Suppose we have an element in the polar coordinate as $0\leq r\leq 1$, $0\leq\theta\leq \alpha$. What are the correct basis functions in this element?  For this kind of mesh, there
are a lot of references on finite difference discretization. But I can not find a fem solution.</p>

Title: How to do FEM in sector elements?


问:
Probabilistic algorithms for matrix approximation
Considering regular matrix approximation inequality

|| $A - QQ^TA $|| < e

where we try to approximate matrix $A$ by a lower rank orthonormal matrix $Q$. I've read an article on probabilistic algorithms for such an purpose. You simply pick random (Gaussian etc.) vectors $\omega$, then generate vectors from the range of $A$ by $y =A\omega$. Finally, construct $Q$ by these vectors, $y$, after Gram-Schmidt etc. Details can be found in "Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions".

My question is the following:

Let's assume we have several matrices $A_i$ with "same structure" and each have the same rank $k$. I may give details on the structure if you interest. And we find $Q_i$ for each corresponding matrix $A_i$. Again assume we may define a distribution over $Q_i$ depending on the distribution we use during approximating matrices $Q_i$ (the distribution of $\omega$).

Then, we get a new matrix $A$ with again same structure and rank $k$. But this time we only have a submatrix $A(:,J)$ i.e. only several columns. I want to approximate $Q$ for this new matrix using only given columns.

To summarize we have two phase:

1. We somehow learn a model for $Q$ using given full matrices $A_i$

2. Given a subsample $A(:,J)$ of a novel matrix $A$ (that is not present in the first phase), we derive its corresponding $Q$

I hope, you may simply recommend me some useful references. I do not even know where to look.

Thank you for your time.


问:
How does GAMG in OpenFOAM really work?
I use OpenFOAM for CFD simulations. A very popular preconditioner is GAMG which needs a low number of iterations per a time step in SIMPLE or PISO solvers that are used to simulate the fluid flow.

I would like to understand how GAMG prepares the internal structures.

Does it create the internal data structures used in GAMG only once in a single iteration on the basis of the geometry and then the entries in these data structures are only updated?

OR

Does it create the interal data structures multiple times in a single iteration?
I read the documentation here but it does not have the answer.

Thanks for your answer.

OwnerDisplayName: luksmir

Title: How does GAMG in OpenFOAM really work?


问:
Finding quadrature weights for a given set of points? How to select points such that all weights are positive?
Currently, I fit a Finite Element solution of a PDE on a spectral basis. The matrices ($R^{25000\times 2000}$) of the corresponding system of linear equations are highly ill-conditioned ($\kappa \approx 10^{15}$). However, I am able to find a solution with Tikhonov-regularization, which has also a small residuum on other points in the domain.

Now, playing around with a toy system, I came up with the idea to find the quadrature weights on a given subset of all points, such that instead of over-fitting with 25000 points, I have to find 2000 points out of my 100000 FE points such that they have positive quadrature weights.

What strategy do you suggest for selecting those 2000 points from the 100000?

Currently I do the following:

Pick 2000

calculate the weights via qr

select all points with negative weights and replace them with new points

On my 1D toy system with chebyshev polynoms, I need to draw about 100000 new points for a spectral basis of 25 modes till all my weights are positive. So this approach is not going to work well.

Any suggestions or literature links are highly welcome.
Thanks.

Here is my toy system:

```<code># encoding: utf8
import numpy as np
import matplotlib.pyplot as pl

def fun(x):
    x = np.asarray(x)
    y = np.zeros(x.shape)
    y[np.abs(x) > 0] = np.exp((-(1. / x[np.abs(x) > 0]) ** 2))
    return y

def gauss_extrema(n):
    i = np.arange(n)
    return np.cos(((2 * (i + 1)) - 1) / (2.*n) * np.pi)

def main():
    Nx = 25
    nx = Nx
    x = 2 * np.random.rand(nx) - 1
#    x = gauss_extrema(nx)
    m = np.arange(Nx)
    mm, xx = np.meshgrid(m, x)
    L = np.cos(mm * np.arccos(xx))

    sp = np.zeros(nx)
    sp[0] = np.pi
    w = np.empty(nx)
    w.fill(-1)
    eps = 0.

    i = 1
    while np.any(w < eps):
        u, s, vh = np.linalg.svd(L.T, full_matrices=False)
        if i % 10000 == 0:
            print i, s[0] / s[-1], np.sum(w < eps)
        beta = np.dot(u.T, sp)
        w = np.dot(vh.T, 1. / s * beta)
        x[w < eps] = 2 * np.random.rand(np.sum(w < eps)) - 1
        mm, xx = np.meshgrid(m, x)
        L = np.cos(mm * np.arccos(xx))
        i += 1

    print w
    print x
    b = fun(x)                # setup right side b
    L = np.multiply(L, w)     # setup weights on all points on L
    u = np.linalg.solve(L, b) # solve L u = b for coefficient vector u

    nx = 100
    x = np.linspace(-1, 1, nx)
    mm, xx = np.meshgrid(m, x)

    L = np.cos(mm * np.arccos(xx))
    L = np.multiply(L, w)
    b = fun(x)

    fig = pl.figure()
    ax = fig.add_subplot(121)
    ax.plot(x, b)
    ax.plot(x, np.dot(L, u))
    ax = fig.add_subplot(122)
    ax.plot(x, np.dot(L, u) - b)

    pl.show()

if __name__ == '__main__':
    main()
</code>```


问:
Sign or cardinality constraint when solving for sparse signal
I'm currently learning about using linear and semidefinite programming to find sparse solutions to problems. In particular, finding sparse solutions where the sampling functions are sinusoidal (Dirichlet kernel).

I've implemented a linear programming solution with a discrete approximation using l1magic MATLAB library, and an SDP solution based on the MATLAB example in the paper "Towards a mathematical theory of super-resolution" someone suggested in a previous question.

The CVX code is:

```<code>cvx_solver sdpt3
cvx_begin sdp
    variable X(n+1,n+1) hermitian;
    X >= 0;
    X(n+1,n+1) == 1;
    trace(X) == 2;
    for j = 1:n-1,
        sum(diag(X,j)) == X(n+1-j,n+1);
    end
    maximize(real(X(1:n,n+1)'*y))
cvx_end
</code>```

In the above paper, if the minimum distance between spikes is greater than or equal to $2 \lambda_c$, where $\lambda_c$ is the wavelength of the cutoff frequency, this condition is sufficient to guarantee the solution is unique and can be obtained to infinite precision (page 5).

I've noticed that if all the spikes are the same sign then shorter distances still give me the correct result. If they are different signs then the above condition seems to apply.

Question:

Is there a way to introduce either a

cardinality constraint, where the number of spikes is known, or

sign constraint, where the signs of the spikes are known

that can give a solution where the minimum distance condition is not met?


问:
How to create synthetic data from known weights
I'm doing some machine learning where I have lots of data and through optimization I'm trying to learn the weights for the model.

I'd like to check that my learning actually works correctly. For that I can create a set of known weights. My question is from these weights, how can I create synthetic data? On this data I can perform the training and check if the weights in the end are the same as the weights I created.


问:
Stationary solution converge but time dependent doesn't
I've coupled a COMSOL model for fluid dynamics with a very simple pde that model the transport of humidity in air.
When I solve it for the stationary case, the solution converge easily, but when I try to solve the time dependent problem, the solution can't reach convergence.
This seems very strange to me, the fluid is laminar and incompressible, the mesh is well designed, the initial values are coherent, all the conditions are smooth.
Is it even possible that COMSOL isn't able to find a solution for the time dependent problem, while he found very easily the time independent one?

Thank you


问:
frozen coefficient vs. constant coefficient
This is a follow up to the question about the method of frozen coefficients. I think it deserves to be a separate question. The frozen coefficient problems are obtained by fixing the coefficients of the variable coefficient problem. For example, consider pde $$u_t=a(x)u_{xx}, x \in \mathbb{R}, u(0,x)=u_0(x)$$, then at every grid point $x_i$ we can have a frozen problem $$u_t=a(x_i)u_{xx}$$. However, even though it is a frozen problem, in a sense that if I write a finite difference scheme, then at every grid point the diffusion coefficient is constant, but this is still quite different from considering the problem $$u_t=au_{xx}, x \in \mathbb{R}, u(0,x)=u_0(x)$$ where $a$ is just a constant number. So, for the latter case, to compute the amplification factor of a finite-difference scheme, I need to write the equation in terms of Fourier transform only once at some grid point $x_i$, but if I freeze the coefficients, I have to consider that at every single node and recompute the amplification factor at each of them. So, if this is the way, a problem with frozen coefficients and a problem with constant coefficient are not equivalent, the constant coefficient problem is rather a subproblem of the frozen coefficient problem, is that correct?

And another related question: if I complicate things more and let $a=a(t,x)$, then by the definition I can't use the Fourier analysis, however, if I freeze the coefficients, it is perfectly applicable? Do I have to take care in a special way of the time-dependent coefficients? I am looking for some intuition behind the method of frozen coefficient. It seems to make everything much simpler but I don't want to abuse that and know my limitations.


问:
Up-/downdating methods for a series of normal equations
In an application I have to solve a series of positive definite linear systems of the form $A^TA x = A^Tb$ (i.e. normal equations). The next system is obtained from the previous one by adding and/or removing a small number of columns to $A$ (usually one or two). The vector $b$ stays the same throughout. It seems clear that some kind of up-/downdating strategy should be applied, however, my questions are:

What kind of up-/downdating shall I choose (QR, Cholesky, Matrix-Inversion-Lemma,others?,...)?

On what conditions and requirements shall I make my choice (high precision needed,
  ill-conditioned matrices, well-conditioned ones, speed,...)?

P.S.: You may assume that $A$ and $A^T$ can be applied fast.


问:
PETSC Memory Leak
I tried to use PETSC direct solver (SUPER LU) to solve the system of equations.
After a long time effort, the code is working properly and gives some good results. The only issue that came up recently is a bout memory leakage.
I found that after iterations RES memory is increasing which leads to the server crash.
I tried bunch of cures which I found on the web, but none of them worked.
I am using MatDestroy and VecDestroy at the end of PETSC part in the code to free all the allocated memories. This is not helpful.
I also tried iterative solvers like Richardson and other direct solver like spooles and the same problem occurred.
I would appreciate if you could help me in this regard.</p>

Title: PETSC Memory Leak


问:
CUDA Mandelbrot Set effective bandwidth and optimization
I was reading through this article (very good article and excellent blog BTW) to do some measurements in my (very simple) implementation of the Mandelbrot Set.

I'm using a Quadro 2000D card which has a theoretical bandwidth of approximately 20 GB/s ( = 1304 MHz * 128 bit-bus).

Mandelbrot image has 700x700 pixels (uchar = 1 byte per pixel). Launching the kernel with a grid(44,44) and block(16,16) the time elapsed to compute it is 1.00026 ms.

Now if my calculations are good then I have the following

$$(490000~bytes)(10^{-9}~GB)~~\big/~~(1.00026)(10^{-3~}s)=0.4898~GB/s$$

How should I interpret this result? I am only using approx 2.5% of the bandwidth, which is bad but since the algorithm is doing only a write per pixel there's no way to improve this, is there? Is there something I am not taking into account?

UPDATE

I did three types of kernels to see what was going on and reduced the size of the image to 250x250.

First the very naive, 1 thread doing all the work: 125 ms + 3.1% Global efficiency store + 2% occupancy. This is BAD.

Second, one thread per row: 0.20 ms + 97.7% Global store efficiency + 16% occupancy. Much better.

Third, one thread per element: 0.25 ms with 50% Global store efficiency + 55% occupancy. Worst than the previous.

So, from here what can I do to improve? Is there way to do it?
The only thing it occurs to me in the third type is to write inside the kernel to shared memory, then synchronize and after that all threads write to global memory so to have coalesced accesses.


问:
Reformulation of a Frobenius norm maximization problem
This question is extended from this question

The original problem is to maximize the Frobenius norm of a matrix product

$max_{X\in \mathbb{D}}||B^TXA  ||_F$, where $\mathbb{D}$ is the doubly stochastic matrix set

One challenge is that the size of matrix $X$ is big and I would like to reformulate it into a low memory requirement formulation.

I have tried to consider the $l_1$ version of the problem.

But yet the performance is not too good. (Mentioned in here)

Any good suggestion?


问:
Newton Iteration method convergence
I wrote a Python code which solves a second degree nonlinear differential equation using the Newton iteration method. The code converges to a 2-cycle within 50 or so iterations.  The cycle only appears in noticable magnitude in the last half dozen grid points (of ~120).

I use Central Differencing for solving a 2 boundary conditions second degree nonlinear equation. The number of gridpoints is 100 and there 3 parameters, $\alpha_i$:

$u'' + \alpha_1 \cdot u' + {\alpha_2 \cdot u \cdot (u')^{2/3}} + {\alpha_3 \cdot (u')^{2/3}} = 0$

The boundary conditions are:

$u(x=0) = 0$

$du/dx(x=N) = m$

The $\alpha$ parameters and m are all calculated from experimental data, and the $\alpha$'s are functions of x.

I use the boundary conditions on the Jacobian matrix at the first and the last line.
I use also these conditions on the calculation of F vector
if the iteration is:
$X = X - J^{-1} \cdot F$

I noticed that the last value of F array change between + and - constant values.  I have already tuned the initial iterate to be near the expected final solution (from experimental data as well).  What options do I have to avoid a 2-cycle in the iterates of newton's method when it appears near to the solution I wish to find?

Moreover, If the F last value (BC2) or one before the end is 0, then I have a convergence, with the expected u values at the last point or last two gridpoints respectively, but still there is a decrement after the #106 point going near to zero.

Title: Newton Iteration method convergence


问:
On solution of a class of discrete-time Lyapunov equation for systems with multiplicaitve noise
Let's consider the following equation
$$X=F_{1}XF_{1}^{T}+...+F_{p}XF_{p}^{T}+C$$
where $p$ is an positive integer and $C$ is a known positive semidefinite matrix. If we augment $F=[F_{1}...F_{p}]$ and $Y=diag (X...X)$, then the equation becomes
$$FYF^{T}−[I...0]Y[I...0]^{T}+C=0$$
seems like a generalized Lyapunov equation. However, there is a constraint on $Y$ for its diagonal form. How to compute $X$?

I met with this problem for dealing with stability analysis for dynamic systems with multiple multiplicative noise.


问:
How much better a bounded BFGS is compared to augmented Lagrangian method with BFGS?
I mean, in handling boxed constraints?

In terms of stability, and more importantly, the numerical performance?

I have already written some well-optimized and well-tested C/CUDA/C++ codes for several unconstrained optimzation methods.

And I use augmented largrangian to handle constraints, however, it looks like bounded BFGS is quite different than the standard BFGS, I cannot do it over my existing BFGS routine,it require almost a complete re-written, I dont really want to do that unless there is significant gain there.

Have anyone, with practical experiences, can tell me if the rewritten is worth doing or I can stick to the let Lagrangian to handle it?

Btw, I usually use BFGS to solve MLE type optimzation problems, where the evaluation of the objective function is costly.


问:
Explanation of LidarBoost Algorithm?
I am trying to understand the LidarBoost algorithm as explained in this paper (PDF warning). My apologies if this isn't the appropriate exchange; I've posted in DSP, but I didn't get a response there, and the FAQ said that some image processing algorithms are acceptable here.

I don't understand how they take the original depth-images $Y_k$ and transform them into the up-sampled images $D_k$. I get how optical flow is used to align the $Y_k$ into a chosen reference frame, but is the transform between $Y_k$ and $D_k$ just a standard image up-sampling with nearest neighbors for the "interpolation" step? If that's the case, can someone explain how they get the term $W_k$ in their data term of the energy function? They say this about it, but I don't understand how to construct $W_k$ for each $k$:

$W_k \in \mathbb{R}^{\beta n x \beta m}$ is a banded matrix that
  encodes the positions of $D_k$ which one samples from during
  resampling on the high-resolution target grid

I'm trying (so far frustratingly ineffectively) to get enough of a sense of the algorithm to try to implement it.

Title: Explanation of LidarBoost Algorithm?


问:
Eigenmode constraints on a complex wave field
I'm studying the usual wave equation:

$$ \frac{\partial^2 \psi}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 \psi}{\partial t^2} = 0$$

with complex field $\psi$, on a n-dimensional cube region of length $L$, which are $n-1$ spatial dimensions of length $L$, and one time dimension of duration $L/c$. Now i proceed to express the field in a Fourier expansion:

$$ \psi(x,y,..t) = \sum_{k, \omega} e^{I ( k_x x + k_y y + ... + \omega t)} \psi_{k,\omega}$$

If i do a cutoff in the expansion so that we only take $2M+1$ cells on each direction in k-space, that is, $- \frac{ \pi M }{L}, - \frac{ \pi (M-1) }{L}, .. 0 .. \frac{\pi M}{L}$ are the allowed values of the $k$ indices (and for $\omega$ too but with a $1/c$ factor somewhere) I am left with $(2M+1)^n$ complex fields to compute, because the field is not real, and the fourier expansion does not satisfy $\psi_{k_x, k_y, ...,\omega} = \psi_{-k_x, k_y, ..., \omega}^{*}$

When i apply the Helmoltz condition on the bulk of the cube, i am left with $(2M+1)^n$ (complex) conditions:

$$(-k_x^2 - k_y^2 - ... + \frac{\omega^2}{c^2} ) \psi_{k_x, k_y,...,\omega} = 0$$

Now, it seems that the Helmoltz equation imposes the same number of contraints as there are unknowns. My question: how can boundary conditions be imposed on this problem without producing an over-determined system


问:
SVD computation with "initial guess"
Suppose I have some matrix A whose SVD I know. Now, I am given B which is A plus some small perturbation. Is there any algorithm which can use the SVD of A as an initial guess to gain speedup in the computation of the SVD of B?

Title: SVD computation with "initial guess"


问:
differences between all the FFTs?
Wikipedia lists a lot of FFT algorithms:

Cooley–Tukey FFT algorithm, Prime-factor FFT algorithm, Bruun's FFT algorithm, Rader's FFT algorithm, and Bluestein's FFT algorithm

What are the pros and cons of each algorithms?
And which one do you prefer?

Title: differences between all the FFTs?


问:
Looking for a java library or algorithm for efficiently implementing the second Chebyshev function
Does anyone know of a java library or algorithm for efficiently implementing the second Chebyshev function?

To be clear, I'm refering to this expression:

$$\vartheta(x) = \sum_{p \le x}\log p$$

$$\psi(x) = \sum_{m=1}^{\infty}\vartheta\left(\sqrt[m]{x}\right)$$

where $\vartheta(x)$ is the first Chebyshev function and $\psi(x)$ is the second Chebyshev function?

Thanks very much,

-Larry


问:
Fenics: Result of Steady state dynamic linear elastic doesn't match with actual values
I solved the steady state dynamic linear elastic model in a solid. My equation is a function of frequency and the strong form is:
$$\operatorname{div}(\operatorname{stress}(\vec x, w)) + w^2  \rho u( \vec x,w)=0$$
where $\rho$ is the density.

The BC is
$$ \sigma(\vec x,w) \cdot n(x)=T(x,w),   \qquad   u(x,w)=U_0$$

The physical problem is a plate with dimension of $1 \times 1 \times 0.1$ with a harmonic load on the all the top and also the bottom of the plate is fixed.
I solved this as a test with the fixed frequency $(w=200)$.

The weak form of the this equation is:
$$\int \sigma(x,w)\cdot \operatorname{grad}(v)\, dv - \rho w^2 \int  v(x)\cdot u(x,w) \, dv - \int v\cdot T(x,w)\, ds(2) = 0 $$

where $v$ is test function, $ds(2)$ is the top of the plate.

I compared the results with the model with converged mesh from Abaqus, I  checked almost everything in the code the order of the displacement is mostly correct but the signs of the displacements are wrong and also the error is large between the solution of the Abaqus and this code.

I checked the variation form and I can't find any error in it. So the only thing that I suspected to be wrong is the Neumann boundary condition but I can't find any bug in it.

I would be very thankful if someone can help me to figure out what is wrong in my code that I get different answers from the code. Or I am not sure if this is a bug and my code is right.


问:
Simple turbulence model appropriate for buoyancy-driven cavity like problem
Which turbulence model is suitable for resolving incompressible buoyancy-driven flow of a fluid within an cylindrical ampoule?

I prefer turbulence model which is sufficiently simple so that fully coupled (UFL) variational form of Navier-Stokes-Fourier is changed in way that each added term

is treated by Newton solver, or

is precomputed prior to assembling tensors at every Newton iteration

so that no additional fixed-point iteration scheme is needed.

You can answer immediately or read description of problem below if useful.

This is coupled problem of flow of melt in the ampoule, flow of air in the furnace (outside of the ampoule) and heat conduction in rest parts of the system (including crystal on bottom of the ampoule).

For purposes of this question is relevant flow of the melt in the ampoule which is similiar to buoyancy-driven cavity (especially in top part) but

problem is cylindrically symmetric, symmetry axis is on the left

boundary conditions for temperature are more complicated than in classical benchmark; furnace wall (right boundary) has non-monotonic temperature profile (see image below - note that interior of the ampoule occupies cca $0.3<x<0.4$, $x$ being vertical coordinate) - increasing with altitude on most of the wall but decreasing near top; top and bottom boundaries of fluid have temperature continuous with outter wall (it is coupled to the rest of the system; bottom boundary is not planar)

almost all material coefficients are temperature dependent

Flow in the ampoule resembles buoyancy-driven cavity in top part where temperature profile induces unstable stratification of fluid. This cavity-like flow is bounded by stably-stratified fluid from below. This was achieved with two orders of magnitude higher viscosity than target one and is non-stationary solution.

It was tested that stationary solution with implicit Euler in time and Newton solver accounting for non-linearity can be computed with three orders higher viscosity than target one. For only two orders higher viscosity stationary solution seems to not beiing stable but time-marching with Crank-Nicolson scheme with small timesteps is possible. I'm currently going to switch to implicit Euler timesteping, lower timestep more, try adding SUPG/PSPG stabilization and check how small viscosity is manageable. But I guess that three orders lower viscosity than that enabling stationary solution will induce turbulence.

For the target viscosity Grashof anfd Rayleigh numbers are
$$ \mathrm{Gr = 3.6E8},$$
$$ \mathrm{Ra = 1.1E9},$$
taking whole length of ampoule as chracteristic length and $850-775 = 75\;\mathrm{K}$
as temperature difference scale; $775\;\mathrm{K}$ is temperature at bottom of the ampoule and $850\;\mathrm{K}$ is peak temperature of temperature profile. Note that this numbers could few orders of magnitude lower taking into account that:

one could take only depth of unstably-stratified region as a characteristic lenght

one could take smaller temperature difference corresponding to unstably-stratified region

peak furnace temperature $850\;\mathrm{K}$ is not reached inside the ampoule


问:
Simulation of asymmetric structures (occupancy = 0.5) unstable
I am trying to simulate a metal-organic framework in LAMMPS using the UFF potential.  It's working quite well for some structures where all molecules have an occupancy of 1.

However, when I have a solvent molecule that has an occupancy of 0.5, what I do is randomly select one of the two positions for each unit cell.  The structure collapses upon simulation.  I'm not sure how to fix this.  My MOF is about 6x3x3 unit cells big; would there be any beneficial reason to try a huge system?

I've also tried simulating at low temperatures, NPT, NVT; it falls apart everytime.  You're not really supposed to use partial charges with UFF, but I tried it with and without.

It's just amazing to me that one structure is so incredibly stable and then a slightly different structure that has 0.5 occupancies falls apart.  Does anyone have any experience with this?


问:
Accurate force-field for MD
I am still a fairly new graduate student, and I am having trouble simulating a metal-organic framework.  I am currently using UFF (I implemented it for LAMMPS), and it works well for some MOFs, but does not give empirically accurate results for others.  What I mean by this is that I have a CIF file for the experimentally obtained structure, and when simulate it, the structure changes.

UFF is fairly old (1992?); there's got to be something better that will allow this to match the experimental structure.  Any ideas?  Currently, I'm obtaining partial charges for the solvent molecules, copper ions, and chlorine from various papers, but I am considering trying to use Gaussian even though I have never done that before.  UFF is generally used without partial charges, but I think the strong electrostatics in my MOF necessitates their use.

The MOF contains carbon, nitrogen, hydrogen, chlorine, oxygen, and copper.

(I'm still working with the possibility that I may have screwed something up somewhere in the simulation process, but I would still think there's something newer/more accurate than UFF).

Title: Accurate force-field for MD


问:
Combining trend estimation and constrained Marquart fit
This title certainly needs some clarification: I need to compute parameters $a_i$ for a helper function $f(\vec{a};k)$ (for grid interpolation) which is fitted to a number of values $y_k$ which are subject to evolution dictated by a system of ODE's $d y_k/d t = F_k(t,\vec{y})$. The Levenberg-Marquardt fit of the values $y_k$ to the function $f(\vec{a};k)$ is repeated in each time step $t_j$. The precise values found for the parameters $a_{i,j}$, which are varied in this fit, does not matter, but I found that in the presence of adaptive step-size control for the ODE this simple algorithm can lead to instabilities as it influences the value of the rhs-functions $F_k(t,\vec{y})$ at higher order. For the present purpose I can assume that the evolution of the parameters is smooth plus some noise $a_{i,j} = A_i(t_j)+\delta_{i,j}$ and the model for the interpolating functions $f(\vec{a};k)$ does not matter. My idea to fix the stability issues is therefore to make the parameters $a_i$ subject to some smoothening procedure.

The concept for this procedure in turn is to estimate the new values $a_{i,J}$ based on the previous ones $a_{i,j}$, $j<J$ and the new $\chi^2$ fit by adding to the function to be minimized in the fit a contribution proportional to $\sum_l \sigma_l (A_l(t_J)-a_{l,J})^2$. Here $\sigma_l$ are positive constant coefficients and $A_l(t_J)$ are the values suggested by the trend of the previous ones (only loosely related to $A_l(t_J)$ above). The obtained $a_{l,J}$ would automatically be consistent with constraints in the fit. (Note that this could be different if I first fit and later smooth) The question is which method to use to predict (to extrapolate) the values $A_l(t_J)$. I wanted to avoid using a simple polynomial or spline fit, because I definitely need to avoid overshooting. On the other hand a second $\chi^2$ fit of a model function seems to overdo it. Also due to the fact that the previous values, $a_{i,j}$, $j<J$, where themselves found based on the combined fit and I am not sure that this does not introduce any bias. What do you think? I also wondered if I could miss use a Sawatzki-Golay smoothing filter (assuming as an approximation that $t_{i+1}-t_i=const$ to avoid the overhead of a full $\chi^2$ fit. But, can it be used for extrapolation?