问:
Using scipy.quad to calculate difficult integral
When evaluating the integral below in python using scipy.quad I get the following warning:

UserWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze
  the integrand in order to determine the difficulties.  If the position of a
  local difficulty can be determined (singularity, discontinuity) one will
  probably gain from splitting up the interval and calling the integrator
  on the subranges.  Perhaps a special-purpose integrator should be used.

Furthermore the real integration bounds should be zero to infinity (see below) but when I change the bounds to this my result looks very different from the finite values used in my example below. How should I go about calculating this integral?

I have the following complex integral to compute numerically:

\begin{align}
\int_0^{\infty} Re \left( \frac{t_1(s)-it_2(s)}{\sqrt{\zeta^2(s)-1}} \right) \exp(i\omega s/V) \, \mathrm{d}s
\end{align}

here $i=\sqrt{-1}$, $t_1$ and $t_2$ are real parametric functions of $s$, $V=\frac{\pi}{2(\pi+2)}$ is a constant, and $\zeta(s)$ is a complex function of $s$. The functions $t_1$, $t_2$, and $\zeta$ are not available in closed form, rather I have computed them numerically. A numpy array containing the sampled function values at various vales of s can be found to download at http://speedy.sh/prG9e/stackexchange.npy. The integrand is problematic because $\zeta(0) = 1$.

I have plotted the integrands for various values of omega:

The singularity in the real integrand is apparent. My work so far:

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

T2 = np.load('stackexchange.npy') #Data: [s, x1, x2, t1, t2]

def pyfunc(z):
    return np.sqrt(z**2-1)

def integrand(s, omega):
    x1 = np.interp(s, T2[:,0], T2[:,1] )
    x2 = np.interp(s, T2[:,0], T2[:,2] )
    t1 = np.interp(s, T2[:,0], T2[:,3] )
    t2 = np.interp(s, T2[:,0], T2[:,4] )
    zeta = x2+x1*1j
    sigma = np.pi/(np.pi+2)
    V = 1/(2*sigma)
    return (-t2*np.real(1j/pyfunc(zeta))+t1*np.real(1/pyfunc(zeta)))*np.exp(1j*omega*s/V)

def integral(omega):
    def real_func(x,omega):
        return np.real(integrand(x,omega))
    def imag_func(x,omega):
        return np.imag(integrand(x,omega))
    a = 0.05 #Lower bound
    b = 20.0 #Upper bound
    real_integral = integrate.quad(real_func, a, b, args=(omega))
    imag_integral = integrate.quad(imag_func, a, b, args=(omega))
    return real_integral[0] + 1j*imag_integral[0]

vintegral = np.vectorize(integral)

omega = np.linspace(-30, 30, 601)
I = integral(omega)

plt.plot(omega, I.real, omega, I.imag)
plt.show()
</code>```

Edit

I found an analytical representation of the integrand in question, defined below. It still seems that there are some numerical difficulties with it though as NaNs are returned when I try to integrate.

```<code>    def pyfunc(z):
        return np.sqrt(z**2-1)

    def func(theta):
        t1 = 1/np.sqrt(1+np.tan(theta)**2)
        t2 = -1/np.sqrt(1+1/np.tan(theta)**2)
        return t1, t2

    def integrand(s, omega):
            sigma = np.pi/(np.pi+2)
            xs = np.exp(-np.pi*s/(2*sigma))
            x1 = -2*sigma/np.pi*(np.log(xs/(1+np.sqrt(1-xs**2)))+ np.sqrt(1-xs**2))
            x2 = 1-2*sigma/np.pi*(1-xs)
            zeta = x2+x1*1j
            Vc = 1/(2*sigma)
            theta =  -1*np.arcsin(np.exp(-np.pi/(2.0*sigma)*s))
            t1, t2 = func(theta)
            return np.real((t1-1j*t2)/pyfunc(zeta))*np.exp(1j*omega*s/Vc)

    def integral(omega):
        def real_func(x,omega):
            return np.real(integrand(x,omega))
        def imag_func(x,omega):
            return np.imag(integrand(x,omega))
        a = 0
        b = np.inf
        real_integral = integrate.quad(real_func, a, b, args=(omega))
        imag_integral = integrate.quad(imag_func, a, b, args=(omega))
        return real_integral[0] + 1j*imag_integral[0]

    vintegral = np.vectorize(integral)
</code>```


问:
Method with low memory requirement for large-scale eigenvalue problem
I am working on the flow stability problem. In this work the main complication is solving generalized eigenvalue problem for a large scale Non-Hermitian matrix. I need only one eigenvalue (most left on real axis), so I used shift-invert transformation. I used SLEPc library to solve this problem, especially Krylov-Shur method with MUMPS LU decomposition. It works fine for small problem (matrix size up to 500k). But real cases give me matrix size about 15M, so I just do not have such amount of memory to solve these cases with MUMPS. I tried to find an iterative method for my problem (method with small memory requirement), but they do not work for me. Generalized Davidson and Jacobi-Davidson are not converged. Also I tried iterative linear solver instead of direct MUMPS for Krylov-Shur method, but they also are not converged. I concentrated on one relatively small problem (about 70k matrix size) that solved fine using MUMPS and cannot find any iterative method to solve it.

Is it possible that iterative eigenvalue solver (or solver with small memory requirement) does not exist for my matrix? I will very appreciate if someone gives me an advice what I can do more to solve my problem.

ContentLicense: CC BY-SA 3.0


问:
Fortran library for minimization or maximization of functions (optimization problem)
I'm looking for a library or individual package which provides convenient optimization subroutine for Fortran, any suggestions?

I'm testing the codes from Numerical Recipe, downhill simplex method and direction-set method, but the input-output interface is not convenient for chi-squared function, because we generally need to find the minimum of chi-squared function within a cube-like parameter space. Thus, I need a more convenient utility especially for least chi-squared problem.


问:
Time-stepper approach to eigenvalue problem
For a linear system

$$
M \dot{u} = Au \qquad \textrm{or} \qquad \dot{u} = L u
$$

The generalized eigenvalue problem is
$$
A e = \lambda M e
$$

We can use the time-stepper approach which essentially computes the eigenvalues of the operator

$$
E = \exp(L T)
$$

If $\mu$ is an eigenvalue of $E$ then it is related to the original eigenvalue by

$$
\mu = \exp(\lambda T)
$$
and $T$ has to be chosen appropriately. If I have a code $E_{\Delta} \approx E$ which solves the initial value problem, then I can couple this to ARPACK to compute the eigenvalues. I would use the standard eigenvalue solver of ARPACK on the operator $E_\Delta$.

But I would also like to use the inner product induced by $M$ in the Arnoldi process so that the eigenvectors have unit norm in the inner product of $M$. How can I achieve this with ARPACK ?


问:
Comparison between Voronoi and Delaunay 2D quadrature methods
This question is a search for further answers from a question on maths.stackexchange.com.

I've inherited some numerical quadrature code that is designed to integrate sparse 2D data. The quadrature is not being done in any of what I would consider the "normal" ways, but seems to be using an ad-hoc method based on Voronoi cell areas. However the data has some unusual properties that may make this sensible.

The samples fall roughly on a grid such that there is either one or zero sample points inside each grid cell. (But randomly positioned inside the cell, not at the cell center).

The integration then proceeds by adding "ghost" points to the exterior of the region (to bound the Voronoi regions at the non-ghost points), generating a Delaunay triangulation of the domain and then assigning a weight for each sample point equal to the area of the corresponding Voronoi region.

I can see that this scheme tries to minimises the error for C1 functions.

However having gone to the trouble of triangulating the domain, I would have expected the integration to be based upon the neighbouring linear elements, which should give exact answers for linear functions, and thus I'd expect an error estimate (at least for C2 functions) of one order of $h$ higher.

However some numerical experiments I performed suggested that for smooth data both methods had similar convergence rates. (Details can be found in my answer to the original math.SE question)

It was noted in one of the answer to a original question that one possible reason for preferring the Voronoi based scheme is continuity of the weights with changes in the sample locations, (which fails for the Delaunay triangular weights as the mesh can change discontinuously with sample locations)

So my questions are:

Are there known advantages to using the Voronoi method over the Delaunay style method, or vice-versa? Can these be quantified?

Why am I seeing roughly the same convergence rates from the two methods?


问:
How to compute frank copula and its derivative accurately?
I need to fit a model using MLE with Frank copula by linking two discrete univate distribution function $u = F(x)$ and $v = F(y)$ together, and the joint distribution function is

$$
\Phi(x,y) = C(F(x),G(y))
$$

where

$$
C(u,v) = \frac{1}{\theta} \log\!\left( 1+\frac{(\exp(-\theta u)-1)(\exp(-\theta v)-1)}{\exp(-\theta)-1} \right)
$$

and end up with the copula density as

$$\begin{align}
c(F(x),G(y)) = &~C(F(x),G(y)) - C(F(x-1),G(y)) - C(F(x),G(y-1)) \\
& + C(F(x-1),G(y-1))
\end{align}$$

Therefore, I have to evaluate $\log c(u,v)$ and ultimately I have to evaluate $C$ and its first derivatives $C_u$, $C_v$ and $C_{\theta}$.

I couldn't get it work with my optimization algorithm (I am using fmin_bfgs and fmin_l_bfgs_b from scipy).

Most of the time, fmin_bfgs try to evaluate some very big point and make my code overflow (get -inf) and end with incorrect result.

I tried to use fmin_l_bfgs_b so that I can put bounds on the variable to prevent overflow, but I end up to only make it works (not overflow) in a too small interval, $[-2,2]$ or $[-1,1]$ in some cases. The resulting point will be on the bound and it must not be the correct solution.

I guess the problem is that $\log c$, $C_u$, $C_v$ and $C_{\theta}$ involve too much exp and log and I lost too much precision similar to the case like $\exp(x)-1$ and $\log(1+x)$.

Are there any better method to compute $C$, $C_u$, $C_v$ and $C_{\theta}$ accurately? (or compute $\log c$ and its derivatives directly?)

I am computing them using:

```<code>def frank_copula(u,v,t):
    if t == 0.:
        return u*v
    else:
        x1 = expm1(-u*t)
        y1 = expm1(-v*t)
        z1 = expm1(-t)
        return log1p(x1*y1/z1)/(-t)

def frank_copula_u(u,v,t):
    if t == 0.:
        return v
    else:
        x = exp(t*u)
        y = exp(t*v)
        y1 = expm1(t*v)
        z = exp(t)
        return z*y1/(z*(x+y1) - x*y)

def frank_copula_v(u,v,t):
    if t == 0.:
        return u
    else:
        x = exp(t*u)
        x1 = expm1(t*u)
        y = exp(t*v)
        z = exp(t)
        return z*x1/(z*(x1+y) - x*y)

def frank_copula_t(u,v,t):
    if t == 0.: # central difference here
        h = 1e-8
        return (frank_copula(u,v,t+h) - frank_copula(u,v,t-h))/(2.*h)
    else:
        x = exp(-t*u)
        y = exp(-t*v)
        z = exp(-t)
        x1 = expm1(-t*u)
        y1 = expm1(-t*v)
        z1 = expm1(-t)
        tmp1 = -z1*(x1*y*v + y1*x*u) + x1*y1*z
        tmp2 = z1*z1 + x1*y1*z1
        tmp3 = log1p(x1*y1/z1)/(-t)
        return (tmp1/tmp2 + tmp3)/(-t)
</code>```

EDIT

I fixed the dependency parameter $\theta$ first, and then find out the parameter of $F(x)$ and $G(y)$ together by maximizing $\log L(p) = \sum \log c_{\theta}(F(x|p),G(y|p))$. They involve the parameter $p$ of dimension 2 to 5.

Here is the output of how fmin_l_bfgs_b algorithm evaluate my (minus) log-likelihood function using my workaround, set the bounds to $[-0.2,1.2]$ and $[-1,1]$:

```<code>L: array([ 0.,  0.,  0.,  0.,  0.]) => 769.446206148
dL: array([ 0.,  0.,  0.,  0.,  0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
        -52.57668841])
L: array([ 1.2,  1.2,  1.2,  1. ,  1. ]) => 8019.20848421
dL: array([ 1.2,  1.2,  1.2,  1. ,  1. ]) => array([ 10020.77504911,   6095.3498353 ,   8824.13266013,   6763.34035373,
         2021.60815303])
L: array([ 0.31272257,  0.31272257,  0.31272257,  0.26060214,  0.26060214]) => 638.623452278
dL: array([ 0.31272257,  0.31272257,  0.31272257,  0.26060214,  0.26060214]) => array([-77.88688381,  54.50079579,  11.01938969,  98.73509615,  78.99317652])
L: array([ 0.41609651,  0.21822845,  0.28865935,  0.12751238,  0.13439021]) => 612.115679161
dL: array([ 0.41609651,  0.21822845,  0.28865935,  0.12751238,  0.13439021]) => array([-92.53682091,  12.87793807, -22.45693087,  26.74828638,  25.65511475])
L: array([  7.05680743e-01,   1.28761254e-01,   3.21734086e-01,
        -3.25757720e-02,   4.93345936e-04]) => 590.768526126
dL: array([  7.05680743e-01,   1.28761254e-01,   3.21734086e-01,
        -3.25757720e-02,   4.93345936e-04]) => array([-15.61044704,  16.3196277 ,  16.76613757,  10.17553022,   1.73461386])
L: array([ 0.81933183,  0.05230577,  0.26232909, -0.07694669, -0.01632839]) => 587.63001378
dL: array([ 0.81933183,  0.05230577,  0.26232909, -0.07694669, -0.01632839]) => array([-10.36049322,   5.82433346,  10.12673749,  -7.60549293,  -3.70091331])
L: array([ 0.93612938, -0.03119036,  0.15964292, -0.05604783, -0.02092306]) => 585.655142497
dL: array([ 0.93612938, -0.03119036,  0.15964292, -0.05604783, -0.02092306]) => array([-1.60415082,  1.31457884,  5.95207277, -9.31825269, -5.19784526])
L: array([ 0.99490642, -0.09547364,  0.05431211,  0.01089551,  0.01278564]) => 584.829795091
dL: array([ 0.99490642, -0.09547364,  0.05431211,  0.01089551,  0.01278564]) => array([-0.44946283, -0.58748774, -1.01623299, -2.83092327,  1.18001099])
L: array([ 0.99380363, -0.10055615,  0.0527463 ,  0.03080186, -0.00557434]) => 584.840769319
dL: array([ 0.99380363, -0.10055615,  0.0527463 ,  0.03080186, -0.00557434]) => array([ 0.29741289,  0.67052685, -0.11044264,  1.88188906, -3.31504582])
L: array([ 0.99442739, -0.09768137,  0.05363196,  0.01954239,  0.00481047]) => 584.813957801
dL: array([ 0.99442739, -0.09768137,  0.05363196,  0.01954239,  0.00481047]) => array([-0.13879393, -0.04846326, -0.63399355, -0.79472897, -0.7858597 ])
L: array([ 0.99245896, -0.09809096,  0.05491007,  0.02184427,  0.00701989]) => 584.811245584
dL: array([ 0.99245896, -0.09809096,  0.05491007,  0.02184427,  0.00701989]) => array([ 0.16175504,  0.3075338 , -0.23277139, -0.09032161, -0.17666573])
L: array([ 0.99043781, -0.09901399,  0.056742  ,  0.02308671,  0.00811308]) => 584.810330203
dL: array([ 0.99043781, -0.09901399,  0.056742  ,  0.02308671,  0.00811308]) => array([ 0.08694667,  0.38532343, -0.11389679,  0.14966267,  0.06054983])
L: array([ 0.98934696, -0.10209532,  0.05898625,  0.02414812,  0.00863789]) => 584.809284488
dL: array([ 0.98934696, -0.10209532,  0.05898625,  0.02414812,  0.00863789]) => array([ 0.04661414,  0.30695628, -0.01908606,  0.20837613,  0.14123964])
L: array([ 0.9893386 , -0.10784629,  0.06164861,  0.02504006,  0.00849419]) => 584.808428091
dL: array([ 0.9893386 , -0.10784629,  0.06164861,  0.02504006,  0.00849419]) => array([ 0.00282033,  0.01810351,  0.00812252,  0.0266529 ,  0.02348458])
L: array([ 0.98941539, -0.10801042,  0.06168044,  0.02499464,  0.00841819]) => 584.80842517
dL: array([ 0.98941539, -0.10801042,  0.06168044,  0.02499464,  0.00841819]) => array([-0.00429622, -0.00085152, -0.00194202,  0.00036292,  0.0010925 ])
</code>```

and here is the output if I use my original bounds $[-5,5]$ (this bound works if I treat $F(x)$ and $G(y)$ as independent):

(Remarks: If I treat them as independent, I just need to evaluate the density function $f$,$g$ and they are much more stable and therefore I guess my problem is due to inaccuracy)

```<code>L: array([ 0.,  0.,  0.,  0.,  0.]) => 769.446206148
dL: array([ 0.,  0.,  0.,  0.,  0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
        -52.57668841])
L: array([ 5.,  5.,  5.,  5.,  5.]) => inf
dL: array([ 5.,  5.,  5.,  5.,  5.]) => array([ nan,  nan,  nan,  nan,  nan])
L: array([ 0.,  0.,  0.,  0.,  0.]) => 769.446206148
dL: array([ 0.,  0.,  0.,  0.,  0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
        -52.57668841])
</code>```

and sorry I just found that fmin_bfgs works and will converge to the same point got by using my workaround method. But it will evaluate some big point and get -inf and nan:

```<code>dL: array([ 0.,  0.,  0.,  0.,  0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
        -52.57668841])
L: array([ 0.,  0.,  0.,  0.,  0.]) => 769.446206148
L: array([ 338.93276474,  114.42079214,  216.55003298,  137.12475499,
         52.57668841]) => inf
dL: array([ 338.93276474,  114.42079214,  216.55003298,  137.12475499,
         52.57668841]) => array([ nan,  nan,  nan,  nan,  nan])
L: array([  3.38932765e-06,   1.14420792e-06,   2.16550033e-06,
         1.37124755e-06,   5.25766884e-07]) => 769.444241862
dL: array([  3.38932765e-06,   1.14420792e-06,   2.16550033e-06,
         1.37124755e-06,   5.25766884e-07]) => array([-338.93116156, -114.42016384, -216.54896297, -137.12398571,
        -52.57656892])
L: array([ 0.00238973,  0.00080675,  0.00152684,  0.00096683,  0.00037071]) => 768.063616272
dL: array([ 0.00238973,  0.00080675,  0.00152684,  0.00096683,  0.00037071]) => array([-337.80294205, -113.97709305, -215.79496551, -136.58082878,
        -52.49109471])
L: array([ 169.46757724,   57.21079945,  108.27577991,   68.56286091,
         26.28852956]) => inf
dL: array([ 169.46757724,   57.21079945,  108.27577991,   68.56286091,
         26.28852956]) => array([ nan,  nan,  nan,  nan,  nan])
L: array([ 0.00238973,  0.00080675,  0.00152684,  0.00096683,  0.00037071]) => 768.063616272
dL: array([ 0.00238973,  0.00080675,  0.00152684,  0.00096683,  0.00037071]) => array([-337.80294205, -113.97709305, -215.79496551, -136.58082878,
        -52.49109471])
L: array([ 338.93276474,  114.42079214,  216.55003298,  137.12475499,
         52.57668841]) => inf
L: array([ 169.46638237,   57.21039607,  108.27501649,   68.56237749,
         26.28834421]) => inf
L: array([ 84.73319119,  28.60519803,  54.13750824,  34.28118875,  13.1441721 ]) => inf
L: array([ 42.36659559,  14.30259902,  27.06875412,  17.14059437,   6.57208605]) => inf
L: array([ 21.1832978 ,   7.15129951,  13.53437706,   8.57029719,   3.28604303]) => inf
L: array([ 10.5916489 ,   3.57564975,   6.76718853,   4.28514859,   1.64302151]) => inf
L: array([ 5.29582445,  1.78782488,  3.38359427,  2.1425743 ,  0.82151076]) => inf
L: array([ 2.64791222,  0.89391244,  1.69179713,  1.07128715,  0.41075538]) => inf
L: array([ 1.32395611,  0.44695622,  0.84589857,  0.53564357,  0.20537769]) => 2010.68946691
L: array([ 0.25288777,  0.08537274,  0.16157439,  0.10231284,  0.03922902]) => 651.514679659
dL: array([ 0.25288777,  0.08537274,  0.16157439,  0.10231284,  0.03922902]) => array([-214.23525183,  -56.37608777, -123.86215026,  -59.8381013 ,
        -29.35530006])
L: array([ 2.45189744, -1.10216675, -0.00860584, -2.25562141, -0.0889141 ]) => 1140.89767258
dL: array([ 2.45189744, -1.10216675, -0.00860584, -2.25562141, -0.0889141 ]) => array([ 200.97116275, -107.57429828,   77.9175913 , -284.63382774,
        222.78613298])
L: array([ 0.60625088, -0.10545527,  0.1342278 , -0.27658817,  0.01863746]) => 641.856514926
dL: array([ 0.60625088, -0.10545527,  0.1342278 , -0.27658817,  0.01863746]) => array([-179.74804571,  -91.13312507, -130.64600443, -147.53451709,
        -47.29257281])
L: array([ 0.94280274, -0.44136287, -0.34555214, -0.04312928,  0.3963997 ]) => 650.978041317
dL: array([ 0.94280274, -0.44136287, -0.34555214, -0.04312928,  0.3963997 ]) => array([-147.33106574,  -91.10716372, -144.5277363 , -121.87325508,
         59.49497911])
L: array([ 0.7272913 , -0.22626399, -0.03832441, -0.19262499,  0.15449914]) => 638.059222532
dL: array([ 0.7272913 , -0.22626399, -0.03832441, -0.19262499,  0.15449914]) => array([-170.21539216,  -91.09157378, -134.59864658, -140.28045749,
        -16.11667578])
...
L: array([ 0.98943719, -0.10800294,  0.06167076,  0.02498195,  0.00840848]) => 584.808425121
dL: array([ 0.98943719, -0.10800294,  0.06167076,  0.02498195,  0.00840848]) => array([ -6.05248630e-08,   2.30028962e-08,   3.00105580e-08,
        -1.12310441e-07,   1.26875060e-08])
</code>```

There may be a possibility that my gradient function is computed wrongly, therefore, I check the my gradient function againist finite difference:

```<code>def approx_gradient(f):
    def wrapped_f(x0):
        n = len(x0)
        G = np.empty((n,),dtype=np.float64)
        h = np.sqrt(np.finfo(np.float64).eps)

        g1 = f(x0)
        for i in xrange(n):
            x1 = np.array(x0,dtype=np.float64)
            x1[i] += h
            g2 = f(x1)
            G[i] = (g2 - g1) / h
        return G
    return wrapped_f
</code>```

Print the error as:

```<code># r1 is gradient computed by code, r2 is gradient computed by above finite-difference
d = np.log10(np.mean(np.abs(r1-r2)))
print("dL: %s => %s, error: %s\n" % (repr(x0),repr(r1),str(d)))
</code>```

and the output:

```<code>dL: array([ 0.,  0.,  0.,  0.,  0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
        -52.57668841]), error: -1.12866497845
dL: array([ 338.93276474,  114.42079214,  216.55003298,  137.12475499,
         52.57668841]) => array([ nan,  nan,  nan,  nan,  nan]), error: nan
dL: array([  3.38932765e-06,   1.14420792e-06,   2.16550033e-06,
         1.37124755e-06,   5.25766884e-07]) => array([-338.93116156, -114.42016384, -216.54896297, -137.12398571,
        -52.57656892]), error: -1.3444523876
dL: array([ 0.00238973,  0.00080675,  0.00152684,  0.00096683,  0.00037071]) => array([-337.80294205, -113.97709305, -215.79496551, -136.58082878,
        -52.49109471]), error: -1.76425541273
dL: array([ 169.46757724,   57.21079945,  108.27577991,   68.56286091,
         26.28852956]) => array([ nan,  nan,  nan,  nan,  nan]), error: nan
dL: array([ 0.00238973,  0.00080675,  0.00152684,  0.00096683,  0.00037071]) => array([-337.80294205, -113.97709305, -215.79496551, -136.58082878,
        -52.49109471]), error: -1.76425541273
dL: array([ 0.25288777,  0.08537274,  0.16157439,  0.10231284,  0.03922902]) => array([-214.23525183,  -56.37608777, -123.86215026,  -59.8381013 ,
        -29.35530006]), error: -4.13744558124
dL: array([ 2.45189744, -1.10216675, -0.00860584, -2.25562141, -0.0889141 ]) => array([ 200.97116275, -107.57429828,   77.9175913 , -284.63382774,
        222.78613298]), error: 4.55716763394
dL: array([ 0.60625088, -0.10545527,  0.1342278 , -0.27658817,  0.01863746]) => array([-179.74804571,  -91.13312507, -130.64600443, -147.53451709,
        -47.29257281]), error: -3.6510702955
dL: array([ 0.94280274, -0.44136287, -0.34555214, -0.04312928,  0.3963997 ]) => array([-147.33106574,  -91.10716372, -144.5277363 , -121.87325508,
         59.49497911]), error: -2.92020205735
dL: array([ 0.7272913 , -0.22626399, -0.03832441, -0.19262499,  0.15449914]) => array([-170.21539216,  -91.09157378, -134.59864658, -140.28045749,
        -16.11667578]), error: -3.52122115041
dL: array([ 0.87382717, -0.37020575, -0.21805711, -0.03796395,  0.05122782]) => array([-164.76068544,  -94.42547784, -141.97757919, -126.19087232,
        -39.8567946 ]), error: -3.19928906949
dL: array([ 0.87817125, -0.41141175, -0.15890429, -0.01557913,  0.05239629]) => array([-153.02101444,  -89.44357123, -130.9150735 , -117.22771105,
        -36.90142624]), error: -3.0131329785
dL: array([ 0.89554758, -0.57623572,  0.07770701,  0.07396015,  0.05707019]) => array([-98.94618778, -66.39110534, -79.97680149, -75.88852736, -21.65497929]), error: -3.5340970988
dL: array([ 1.31009037,  1.3391422 , -0.80615157, -0.3413817 , -0.09094827]) => array([  88.75622027,  103.66551146,   28.52931607,   59.92216422,
         34.82483184]), error: -3.86721832568
dL: array([ 1.03325948,  0.06005647, -0.2159125 , -0.0640172 ,  0.00789817]) => array([-57.66593306, -27.77908949, -59.01722568, -47.64785975, -13.70147329]), error: -3.21832537521
dL: array([  9.95932857e-01,  -1.40930257e-01,   1.39440628e-01,
         4.05132617e-02,   5.47288338e-04]) => array([ 27.98975787,  12.98469404,  25.68770109,  19.28943637,   5.48862252]), error: -4.19602531575
dL: array([ 0.98974304, -0.10584977,  0.05390259,  0.02224462,  0.00674051]) => array([-3.07471251, -1.62491769, -2.83340373, -2.41904387, -1.23502546]), error: -3.84209274725
dL: array([ 0.98982872, -0.10776008,  0.0609875 ,  0.0246589 ,  0.00848128]) => array([-0.07129277, -0.05980753, -0.10964847, -0.12157649,  0.00033618]), error: -3.57935797704
dL: array([ 0.98941444, -0.10804747,  0.06173506,  0.02499305,  0.00836133]) => array([-0.00140251, -0.00150752,  0.00292727,  0.00098813, -0.01230324]), error: -4.19996237842
dL: array([ 0.98944092, -0.10800269,  0.06166992,  0.0249806 ,  0.00840967]) => array([ 0.00141414,  0.000485  ,  0.00088848,  0.00042683,  0.00068224]), error: -4.05494671838
dL: array([ 0.98943748, -0.10800296,  0.06167061,  0.02498174,  0.00840839]) => array([  2.32462795e-05,  -1.47758964e-05,  -5.04953544e-06,
        -4.29179409e-05,  -1.69174238e-05]), error: -4.30601047846
dL: array([ 0.9894372 , -0.10800295,  0.06167076,  0.02498195,  0.00840848]) => array([  3.89230032e-07,  -8.59830965e-07,  -1.92379457e-07,
        -1.03506122e-06,  -1.88641715e-06]), error: -4.04379091415
dL: array([ 0.98943719, -0.10800294,  0.06167076,  0.02498195,  0.00840848]) => array([ -6.05248630e-08,   2.30028962e-08,   3.00105580e-08,
        -1.12310441e-07,   1.26875060e-08]), error: -3.67971485116
</code>```

OK. I think my updated question would be:

I still think my gradient does not compute accurately in the box $[-5,5]^n$. I will approximate fisher information using my gradient function, I couldn't get a semi-positive definite information matrix.

Are there any way to avoid BFGS to evaluate some "big" points?

Are there any way to compute the copula density by computing $f$ and $g$ only for discrete margin?


问:
Compute drag force on a mesh in SU2 CFD (beginner)
I am a complete beginner in computational fluid dynamics. I have a 3D mesh in Blender representing the object I would like to do my tests on. I can export it in (I think) any useful format for 3D. I would like to know where I can start to make this kind of simulation with SU2 CFD:

External flow

Air at standard conditions

Air flows against my mesh at some speed (say, 10 m/s)

I would like to know what the drag force on my object will be

Is it possible in some easy way? I tried to have a look around but couldn't find any tutorial or example.

Thank you


问:
Package for integration over non-rectangular region
I want to compute the expected value of a multivariate function f(x) wrt to dirichlet distribution. My problem is "penta-nomial" (i.e 5 variables) so calculating the explicit form of the expected value seems unreasonable. Is there a way to numerically integrate it? Is there some package (e.g in R) that can let me do that? I realize I can solve it using first principle but I want the computation to be really fast.

```<code>f(x) = \sum_{0,4}(q_i*log(n/q_i))
</code>```


问:
Why is my lower convex hull extraction algorithm not working?
Recently, I wrote an algorithm to obtain a delaunay triangulation of a random point set in $I=[-10,10]$x$[-10,10] \subset R^2$ by projecting these points onto the 3 dimensional paraboloid $z=x^2+y^2$, computing the convex hull, and then projecting the lower hull back onto $R^2$.  I have verified that my convex hull algorithm is producing the correct result.  Thus, I'm sure the problem is with my approach to obtaining the lower hull.

My lower hull extraction algorithm is similar to the incremental construction algorithm for generating convex hulls.  That is, if we chose a point at negative infinity and shined a light from this point up toward the hull, we should be able to see the entire lower hull and be able to extract it completely.  Computationally, of course, we cannot really use negative infinity as a number.  But there must be some optimal point under the hull where the entire lower hull should be visible. (note: visibility here is determined by the signed orientation).

In my problem, since the point set is close to the origin (x,y)=(0,0), I chose to look for an optimal viewpoint along the negative z-axis.  I figured that since the points on the lower hull are on the paraboloid $z=x^2 + y^2$, we can exploit the fact that the paraboloid is monotonically increasing as the (x,y) input coordinate moves in the radial direction away from the origin (0,0).  Thus, the point in the point set whose x,y coordinate is furthest from the origin (0,0) will not only have the highest z-coordinate on the paraboloid, but also have the steepest tangent plane.  Furthermore, where this tangent plane intersects the negative z-axis should be the optimal point which can see the entire lower hull.

When I experimented with this algorithm, it doesn't seem to work perfectly.  That is, it is able to extract almost all the facets of the lower hull.  It only seems to run into problems extracting facets whose projected delaunay triangles (onto $R^2$) are extremely thin and located on the 2D hull (that is, the outer boundary triangles).  In these cases, i have had to lower my initial viewpoint by as much as $10^5$ or more, depending on the thinness of the triangle on the outer boundary of the delaunay triangulation in 2D.

I know that there are other approaches to obtaining the lower hull, but I'm really more interested in understanding where the flaw is in my original algorithm.  Since the paraboloid is monotonically concave up, shouldn't the optimal viewpoint be located at a point on the z-axis where the tangent plane from $(x^*,y^*,z^*)$ hits the z-axis, where $(x^*,y^*,z^*)$ is the point projected from the given point set (in $R^2$) onto the paraboloid (in $R^3$ whose (x,y) coordinate is furthest from the origin?


问:
Produce one smooth curve on one triangle mesh
I hope to get one smooth curve on one triangle mesh. I get one path on the mesh at first. The path consists of vertices of the mesh. I can see the path from the image below.

Each one green dot indicates one vertex on the mesh. The method that I am using is to fit one smooth curve using the vertices and project the smooth curve on the mesh. To project the smooth curve, I find the closest point on the mesh for each vertex of the smooth curve.

However the curve is not on the mesh because it intersects with the mesh. I realize that the vertices need to move along the edge of the triangle shared with the previous vertex.

In my view, there are two kinds of methods to produce the smooth curve. The first kind is to move the points of the curve along the mesh to make the curve smooth. The second kind is to fit the points to obtain one smooth curve and project the smooth curve on the mesh,but the projection should be complex.

Can you give me some advices or links to websites or papers? Thanks in advance.


问:
total memory usage of MPI shared memory
I am trying to use the MPI share memory feature. I have several SMP nodes, and each of them has four cores. I need an array of size N for each node that should be accessed by all four cores in each node. My plan is to construct a shared window of size N/4 using MPI_Win_allocate_shared, and I expect that memory usage of each node would be N. In the example below, N is 4X10^9 bytes, but the memory usage of each node is not 4GB but 16GB. Am I missing something?</p>

```<code>#include <iostream>
#include <mpi.h>

int main(int argc, char** argv) {
   MPI_Init(&argc, &argv);

   int rank_all;
   int rank_sm;
   int size_sm;

   // all communicator
   MPI_Comm comm_sm;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank_all);

   // shared memory communicator
   MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &comm_sm);
   MPI_Comm_rank(comm_sm, &rank_sm);
   MPI_Comm_size(comm_sm, &size_sm);

   std::size_t local_window_count(1000000000);

   char* base_ptr;
   MPI_Win win_sm;
   int disp_unit(sizeof(char));
   MPI_Win_allocate_shared(local_window_count * disp_unit, disp_unit, MPI_INFO_NULL, comm_sm, &base_ptr, &win_sm);

   // write
   char buffer;
   if (rank_sm == 0) {
      buffer = 'A';
   }
   else if (rank_sm == 1) {
      buffer = 'C';
   }
   else if (rank_sm == 2) {
      buffer = 'G';
   }
   else {
      buffer = 'T';
   }

   MPI_Win_fence(0, win_sm);

   for (std::size_t it = 0; it < local_window_count; it++) {
      base_ptr[it] = buffer;
   }

   MPI_Win_fence(0, win_sm);

   // read
   long long int index_start(-1 * rank_sm * local_window_count);
   long long int index_end((size_sm - rank_sm) * local_window_count - 1);

   for (long long int it_rel = index_start; it_rel < index_end; it_rel++) {
      buffer = base_ptr[it_rel];

      if (it_rel == index_start) {
         std::cout << rank_sm << " start: " << buffer << std::endl;
      }
      else if (it_rel == (index_end - 1)) {
         std::cout << rank_sm << " end: " << buffer << std::endl;
      }
   }

   MPI_Finalize();

   return 0;
}
</code>```

Title: total memory usage of MPI shared memory


问:
Solvers for stiff initial value ODEs with sparse Jacobian
What ODE solvers are optimized for solving stiff systems with sparse Jacobian? Such systems appear, for instance, when a parabolic PDE is discretized in space using typical finite difference or finite element methods.</p>


问:
Solving a nonlinear equation with a Markov process and RVs
Assume that we have the following equation and the following assumption. The scope is to solve for some particular variables expressed later.

Update

$$E_{t}\left[ b(A_{t+1})^{1-\gamma} *R_{t+1}^{-\gamma}*[r_k(A_{t+1},a)-r_h(A_{t+1})+\eta)] \right]=0$$

where $R=(1-\theta_1-\theta_2)*r^f + \theta_1*r_k(A,a)+\theta_2 *[r_h(A)+\eta]$

Also consider that the number $b(A)$ is pinned down recursively by the following equation:

$$b_t = G[E_t(b_{t+1}*R_{t+1}^{1-\gamma})]$$, where $G(\cdot)$ is a known function

Main Assumptions

$A_t$ is a random variable that follows a Markov process, and we know its value at time $t$, but not at time $t+1$.

$a$ is an $IID$ random variable lognormally distributed with some mean $\mu_a$ and variance $\sigma^2_a$.

$\eta$ is another $IID$ random variable normally distributed with some  mean $\mu_{\eta}$ and variance $\sigma^2_\eta$.

The RVs are independent of each other. That is, $a$ is independent of $\eta$, while these two can be correlated with $A$, thought are still independently drawn from $A$.

$b(A)$ is a number that depends on $A$ and I do not know their functional relationship.

$r^f$ is a known number.

by $r_k(A,a) $ and respectively $r_h(A)$ I mean a functional relationship (which I know) between these two that depend on $A, a$  and $A$ respectively.

To sum up, the expectations operator is about the unknowns in the next period $A,a$ and $\eta$. To fix the ideas. Let's suppose a simple case, where $a,\eta$ take two values and $A$ follows a two-state Markov process and also take two values. Let's also suppose that at time $t$ we are in one of these states. In that case, the Expectations operator is defined over the next possible state for $A$, and the possible values for the $IID$ RVs.

THE COMPUTATIONAL TASK: Solve for $\theta_1$ and $\theta_2$

Apparently, the solution for each $\theta$ will depend on the possible values of $A$, so I expect to be a vector.

Any guidance on how to solve this equation will be appreciated.

ContentLicense: CC BY-SA 4.0


问:
Maximize sum of Rayleigh quotients
I want to maximize the sum of Rayleigh quotients:

$$\max_x\sum_{i=1}^n\frac{x^\top A_i x}{x^\top B_i x}$$

where $A_i$ and $B_i$ is positive definite. I've found a similar question here:
minimization problem: sum of Rayleigh quotients. But it only deals with the case $n=2$, and therefore it seems that the methods there are not applicable for my problem.

Is there any suggestion for this problem? Does gradient descent/ascent work for this problem?

Thanks!

Title: Maximize sum of Rayleigh quotients


问:
Euler Equation Eigensystem with Gravity in the Energy Flux
I am modifying a conservative form of the Euler equations with gravity in the energy flux (see previous question: Energy Conservation in Conservation Laws with Source Terms) for use in a Riemann solver.  I'm seeing some unphysical/unstable results despite total energy conservation to roundoff.

The conservation law equation with source term (a.k.a balance law)
$\frac{\partial \vec{Q}}{\partial t} + \nabla \cdot \vec{F} = \vec{S}  $, is solved, where the variables are defined in $\vec{Q}$,the flux, $\vec{F_z}$ in the $z$-direction, and the source $\vec{S}$ of modified equations are:

$ \vec{Q}= \left[\begin{matrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ E_{total} \end{matrix}\right] \quad \quad \vec{F_z}= \left[\begin{matrix} \rho w \\ \rho uw \\ \rho vw \\ \rho w^2+p \\ (E_{total} + p)w \end{matrix}\right] \quad \quad \vec{S}= \left[\begin{matrix} 0 \\ 0 \\ 0 \\ -\rho g \\ 0 \end{matrix}\right]$

where $\rho$ is the density, $u$,$v$, and $w$ are velocities in the $x$, $y$ and $z$ directions respectively.  The total energy is defined as $E_{total} = \frac{p}{\gamma-1} + \frac{1}{2}\rho q^2 + \rho\psi$, with $q^2=u^2+v^2+w^2$, and $\psi=gz$, with $g$ being the acceleration due to gravity in the $z$-direction, and $\gamma$ is the ratio of specific heats.  The pressure $p$ can be defined from the total energy, $p=(\gamma-1)(E_{total} - \frac{1}{2}\rho q^2 - \rho\psi)$.

I find the flux Jacobian, $A = \frac{\partial \vec{F_z}}{\partial \vec{Q}}$,

$ A = \left[\begin{matrix}
0 & 0 & 0 & 1 & 0 \\
-uw & w & 0 & u & 0 \\
-vw & 0 & w & v & 0 \\
(\gamma-1)(\frac{1}{2} q^2 - \psi) - w^2 & -(\gamma-1)u & - (\gamma-1)v & (3-\gamma)w& (\gamma-1) \\
[(\gamma-1)(\frac{1}{2} q^2 - \psi) - h']w & -(\gamma-1)uw & -(\gamma-1)vw &  h'-(\gamma-1)w^2 & \gamma w
 \end{matrix}\right]$

where $h' = h + \psi$, with $h = \frac{a^2}{\gamma-1} + \frac{1}{2}q^2$, and $a^2 = \frac{\gamma p}{\rho}$.

This systems yields eigenvalues $\lambda =[u-a, u,u,u,u+a] $, and right eigenvectors:

$R = \left[ \begin{matrix}
1    & 0 & 0 & 1 & 1 \\
u    & 1 & 0 & u & u \\
v    & 0 & 1 & v & v \\
w-a  & 0 & 0 & w & w+a \\
h'-ua& u & v & \frac{1}{2}q^2+\psi & h'+ua
\end{matrix}\right]$

and left eigenvectors:

$L = R^{-1} = \frac{\gamma-1}{2a^2}\left[ \begin{matrix}
(\frac{1}{2}q^2-\psi)+\frac{wa}{\gamma-1} & -u & -v & -w-\frac{a}{\gamma-1} & 1 \\
\frac{2ua^2}{\gamma-1} & -\frac{2a^2}{\gamma-1} & 0 & 0 & 0 \\
-\frac{2va^2}{\gamma-1} & 0 & \frac{2a^2}{\gamma-1} & 0 & 0 \\
2h'-2q^2 & 2u & 2v & 2w & -2 \\
(\frac{1}{2}q^2-\psi) - \frac{wa}{\gamma-1} & -u & -v & -w + \frac{a}{\gamma-1} & 1
\end{matrix}\right]$

Does anyone see any problems with this formulation?  Does anybody have experience including gravity as a source term for the momentum equation, and as a flux for the energy equation?  Or does anyone see problems with the calculation of the eigensystem?

Thanks for the help.


问:
Mixed DG for Poisson with mixed BC's
I am trying to find a good reference on a proper weak formulation for mixed DG (Raviart Thomas and DG) formulation for a Poisson equation with mixed boundary conditions. Can anyone suggest a good reference?

Since what follows is lengthy, I am happy to offer a modicum of recompense (of your choice) for good answers.

Basically I am trying to solve is:
$$
\begin{align}
\mathbf{u} + \mathbf{k} \cdot \nabla p &= 0 \quad in \quad\Omega \\
\nabla \cdot \mathbf{u} &= q  \quad in \quad \Omega \\
p &= g  \quad on \quad \partial  \Omega_D\\
\mathbf{u} \cdot \mathbf{n} &= f \quad on \quad \partial \Omega_N
\end{align}
$$

So I cleaned things up a bit
$$
\begin{align}
\mathbf{k}^{-1} \cdot \mathbf{u} +  \nabla p &= 0 \quad in \quad\Omega \\
-\nabla \cdot \mathbf{u} &= -q  \quad in \quad \Omega \\
p &= g  \quad on \quad \partial  \Omega_D\\
\mathbf{u} \cdot \mathbf{n} &= f \quad on \quad \partial \Omega_N
\end{align}
$$

multiplied the equations with the test functions $\mathbf{v}$ and $w$ and integrating by parts:
$$
\begin{align}
\int_{\Omega} (\mathbf{k}^{-1} \cdot \mathbf{u}) \mathbf{v}\mathrm{d}x + \int_{\Omega_i} \{p\} \mathbf{v} \cdot \mathbf{n} \mathrm{d}s + \int_{\Omega_D} g \mathbf{v} \cdot \mathbf{n} \mathrm{d}s + \int_{\Omega_N} p \mathbf{v} \cdot \mathbf{n} \mathrm{d}s - \int_{\Omega} p \nabla \cdot \mathbf{v} \mathrm{d}x= 0
\end{align}
$$
and
$$
- \int_{\Omega_D} (\mathbf{u} \cdot \mathbf{n}) w \mathrm{d}s - \int_{\Omega_N} f w \mathrm{d}s + \int_{\Omega} ( \mathbf{u} \cdot \nabla w) \mathrm{d}x = - \int_{\Omega} qw \mathrm{d}x
$$

where $\{p\}=0.5(p^++p^-)$ . Since $p$ is discontinuous, I need a trace term (I think). What about $\mathbf{u}$ ? since it is a Raviart Thomas element, do I need to do anything for it? Should I change $\mathbf{u}$ to be DG (of order degree($p$)+1) as well? Does it matter?

I get the following weak formulations:
$$
\begin{align}
A(\{\mathbf{u},p\},\{\mathbf{v},w\}) &= (\mathbf{k}^{-1} \cdot \mathbf{u},\mathbf{v})_{\Omega} + (\{p\},\mathbf{v}\cdot \mathbf{n})_{\Omega_i}+ (p,\mathbf{v}\cdot \mathbf{n})_{\Omega_N} - (p, \nabla \cdot \mathbf{v})_{\Omega} \\
 &- (\mathbf{u} \cdot \mathbf{n}, w)_{\Omega_D} + (\mathbf{u}, \nabla w)_{\Omega} \\
F(\{\mathbf{v},w \}) &= -(g,\mathbf{v} \cdot \mathbf{n})_{\Omega_D}+(f,w)_{\Omega_N}- (q,w)_{\Omega}
\end{align}
$$

I perform surface integrals where necessary, and all boundary conditions are imposed weakly. Would this be a stable formulation? If I do not do integration by parts for the second equation $-{\textrm{div}}\ {\mathbf u} = -q$, would imposing $\mathbf{u}\cdot \mathbf{n} = f$ on $\partial \Omega_N$ strongly be better/worse or is it irrelevant?

Thank you very much in advance.

Title: Mixed DG for Poisson with mixed BC's


问:
numerical solver for stochastic optimal control problems
can any one recommend numerical solver (c/c++ library preferred) for stochastic optimal control problems? For deterministic optimal control I found something like that: http://abs-5.me.washington.edu/ocp/index.html Now looking for something similar for stochastic OCP.</p>


问:
4th order tensor
I'm new with FEniCS and Python and I'm stuck with this issue:
is there a way to write a 4th order tensor in an easy way to implement?

I have to compute the following stiffnes tensor:

$A_{ijkl}= \lambda \delta_{ij}\delta_{kl} + \mu (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$

Title: 4th order tensor


问:
Numerical integration and filtering of acceleration experimental data
I have a vector containing acceleration measurements and the corresponding vector of times in which measurements are taken.

To obtain velocity and displacement I used the cumtrapz() function already implemented in MATLAB. This should be fine, as I read here Numerical integration of non-uniform acceleration samples and here https://stackoverflow.com/questions/9881430/numerical-integration-using-simpsons-rule-on-discrete-data.

My doubts are related to the disturbances introduced by the integration algorithm as I found in a book a chapter in which the author says that the trapezium rule introduces low-frequency disturbances.

I tried to remove these disturbances and unwanted 0-frequency components (the acceleration was not 0-mean therefore I had a linear trend in the velocity) applying a 2nd order, highpass, butterworht filter with cut-off frequency at 5 Hz (in both directions to null the phase shift) to the velocity and the displacement but I don't know if this is enough or too much because I need very accurate values for velocity and displacement since I have to use them to plot a surface in the state-space (EDIT: phase-space).

How can I choose the correct cut-off frequency in order not to loose any important information knowing that the experimental structure is excited in a range of frequencies from 5 to 100 Hz?

EDIT
Since I have to integrate noisy data, is it advisable to use an higher order method if my goal is accuracy?


问:
Eigenvalue analysis of preconditioned partial differential operator
today, I encountered a confused problem by accident, but I have no ideas to deal with it fully. The question can be described as follows: for example, when we need to use FDM/FEM to discrete the partial differential operator with variable coefficients,
$\mathcal{L} = d_1(x,y)\frac{\partial^2 u(x,y)}{\partial x^2} + d_2(x,y)\frac{\partial^2 u(x,y)}{\partial y^2} + v(x,y)u(x,y)$, where the coefficients $d_1(x,y),d_2(x,y),v(x,y)\geq 0$. For simplicity, we use the following notations to explain the discretization by FDM:
\begin{equation}
d_1(x,y),\ d_2(x,y),\ v(x,y) \rightharpoonup D_1, D_2, V\ \ (\mathrm{three\ diagonal\ matrices,\ whose\ entries\ are\ all\ nonnegative})\\
\frac{\partial^2 u(x,y)}{\partial x^2}\rightharpoonup L_1 (\mathrm{discretized\ matrix})\\
\frac{\partial^2 u(x,y)}{\partial y^2}\rightharpoonup L_2 (\mathrm{discretized\ matrix})\\
\end{equation}
so we can obtain the discretized matrix of $\mathcal{L}$ as follows,
$A:= D_1L_1 + D_2L_2 + VI$ ($I$ is the identity matrix). Maybe, we want to construct a preconditioner as $P:=d_1L_1 + d_2L_2+\hat{v}I$, where $d_1,d_2,\hat{v}$ are all the mean value of diagonal entries (for example, if $D_1 = diag(s_1,s_2,\ldots,s_n)$, then $d_1 = \frac{s_1 + s_2+\ldots + s_n}{n}$). in fact, this preconditioner can be regard as the discretized matrix of $\mathcal{L} = d_1\frac{\partial^2 u(x,y)}{\partial x^2} + d_2\frac{\partial^2 u(x,y)}{\partial y^2} + \hat{v}u(x,y)$, so I want to know that:

(1) Is this idea reasonable?

(2) if so, how to analyze the spectrum of preconditioned matrix $P^{-1}A$, for the sake of simplicity, take $v(x,y) = 0,$, I can write
\begin{equation}
(d_1L_1 + d_2L_2)^{-1}(D_1L_1 + D_2 L_2) = I + (d_1L_1 + d_2L_2)^{-1}[(D_1 - d_1I)L_1 + (D_2-d_2 I) L_2],
\end{equation}
it means that there are some eigenvalue closed to 1 and this can translate a rapid convergence of Krylov subspace methods? However, this analysis seems to be very rough. So are there some more detailed spectral analysis of preconditioned matrix ? Maybe, some one can provide some references to explain this idea ?

(3) Maybe, some other boundary condition can be also considered.

I have do a small test problem: the spectral distributions of original and preconditioned matrix, please see the following figures:


问:
Two Dimensional Nozzle Flow
I've been looking at the book "Computational Fluid Dynamics - The Basics With Applications", by John D. Anderson, and his solution for the one-dimensional subsonic-supersonic nozzle flow using the MacCormack scheme. I was wondering if there are any book chapters or academic papers out there explaining, in the level of detail given by Anderson, how to solve two dimensional nozzle flow. I found the following AIAA paper, but cannot access it since I'm unfortunately not a member: http://arc.aiaa.org/doi/abs/10.2514/3.61365

Any help would be greatly appreciated.

Title: Two Dimensional Nozzle Flow


问:
Iterative algorithm prove precise conditions for convergence
Question: Consider the iterative improvement algorithm below.

Starting with $Az_i = r_i$ and $(A + E)\hat{z}_i = r_i$ derive a formula showing how the absolute error in the $(i + 1)^{st}$ iterate $\Vert\hat{x}_{i+1} - x\Vert$ is bound by a multiple of the absolute error in the $i$-th iterate $\Vert\hat{x}_{i} - x\Vert$. Argue that the magnitude of this multiple is dependent on the condition of A, and is less than one (hence convergence) if $A$ is well-conditioned.

My work:

I know the following:

$$\frac{1}{\kappa(A)} \frac{\Vert r_i\Vert}{\Vert b\Vert} \le \frac{\Vert\hat{x}_{i} - x\Vert}{\Vert x\Vert} \le \kappa(A) \frac{\Vert r_i\Vert}{\Vert b\Vert}$$
being $\kappa(A)$ the condition number of $A$.

Now do:
$$
\Vert\hat{x}_{i+1} - x\Vert = \Vert(\hat{x}_{i} + z_i) - x\Vert \le \Vert\hat{x}_{i} - x\Vert + \Vert z_i\Vert \le \frac{\Vert\hat{x}_{i} - x\Vert}{\Vert x\Vert} \le \kappa(A) \frac{\Vert r_i\Vert\, \Vert x\Vert}{\Vert b\Vert} + \Vert A^{-1}r_i\Vert$$

$$= \Vert A\Vert\, \Vert A^{-1}\Vert \frac{\Vert r_i\Vert\, \Vert x\Vert}{\Vert b\Vert} + \Vert A^{-1}r_i\Vert \le \Vert A\Vert \Vert A^{-1}\Vert \frac{\Vert r_i\Vert\, \Vert x\Vert}{\Vert b\Vert} + \Vert A^{-1}\Vert\, \Vert r_i\Vert$$

$$= \Vert A^{-1}\Vert\, \Vert r_i\Vert \left(\Vert A\Vert \frac{\Vert x\Vert}{\Vert b\Vert}\right)$$

What I did is not correct because I need to bound the absolute error of $||\hat{x}_{i+1} - x||$ by a multiple of  $||\hat{x}_{i} - x||$ which I am not sure how to do.


问:
What is the relation between Kruskal tensor and CP decomposition?
In Matlab Tensor Toolbox there is a tensor type called "Kruskal tensors", I found its form is similar to the CP decomposition.

Wikipedia mentioned:

"As such, many of the methods have been independently invented
  several times, often with subtle variations, leading to a confusing
  literature. Abstract and general mathematical theorems are rare
  (though see Kruskal[1] with regard to the CP decomposition)"

In a paper on [A concise proof of Kruskal's theorem on tensor decomposition]， it mentioned:

"...In these fields, the expression of a tensor as a sum of rank-1
  tensors is often referred to as the Candecomp or Parafac
  decomposition."

So, what is the relation or difference between Kruskal tensor decomposition and CP decomposition?

Are they actually the same thing but Kruskal decomposition is devised by a different person in different field for third-order tensor?


问:
Poisson equation with pure Neumann boundary conditions (using FEM)
I am trying to derive the correct variational form for the Poisson equation with pure Neumann boundary conditions, and an additional contraint $\int_{\Omega} u \, {\rm d} x = 0$, as described in this link.

As I understand it, the problem can be defined as:

Find $(u, c) \in V \times \mathbb{R}$ such that

$F(u, c) = \frac1{2} \int \nabla u \cdot \nabla u - \int f u -
 \int_{\partial \Omega}g u + c \int u$

is minimized wrt $u$ and $c$

To minimize with respect to $u$, I started writing

$\left. \frac{\partial}{\partial \epsilon}F(u + \epsilon v, c)\right|_{\epsilon=0} = 0$

which led to

$\int \nabla u \cdot \nabla v + \int c v = \int f v + \int_{\partial \Omega} g v $.

With respect to $c$ this gives $\int u = 0$.

Just adding these two equalities leads to the variational form mentioned in the link (where $d \int u = 0$ for all $d \in \mathbb{R}$ ensures that the two equations remain true separatly). Is this correct so far?

Now, in the finite element description, we have one more unknown, $c$, i.e. the Lagrange multiplier, no? To go to a set of equations, the variational form of the equation must be true for all $N$ testfunctions (shapefunctions) $v$. Am I right that this will lead to a $N \times (N + 1)$ matrix equation? Am I also correct that, since the expression should be valid for all $d \in \mathbb{R}$ we can pick a random value in the actual evaluation of these forms?


问:
Plotting renormalization group flow diagram from recursion relation
How do I plot a flow diagram (given certain initial conditions, the trajectory in y-x space) for a RG (renormalization group) flow given by the following recursion relations?

\begin{eqnarray}
x^{\prime}&=&\frac{x(1+y)^2}{(x+y)(1+xy)} \\
 y^{\prime}&=&\frac{y(x+y)}{1+xy}  \\
\end{eqnarray}

I have explored several options such as RSolve, Equation trekker package on Mathematica, but none of these can be used to plot a flow diagram for a continuous transformation.


问:
Piecewise linear approximation of an experimental data curve. Knots position constrained
I need to fit a curve, obtained from experimental data, with a piecewise linear model (4 knots and therefore 3 lines).

I tried using the MATLAB function

```<code>coeff = lsqnonlin(@fun,X0,[],[],[],x,y);
</code>```

where $fun$ is a function defined as:

```<code>[f] = fun (X,x,y)
A = X(1);
.
.
.
F = X(6);

if x(i) <= kn1
    f(i) = A+B*x(i)-y(i);
elseif x(i) > kn1 && x(i) <= kn2
    f(i) = C+D*x(i)-y(i);
else
    f(i) = E+F*x(i)-y(i);
end
</code>```

In this way I should find the coefficients $A,B,C,D,E,F$ which minimize the difference between the model and the real curve (represented by $y$). $X$ is the coefficients vector.

The main problem I have is that the knots $kn1$ and $kn2$ cannot be placed anywhere: usually MATLAB places the knots at the same distance one from the other, but in my case I need that the knots are placed in the exact point where the real curve changes its slope (of course the change in the slope does not occur suddenly, it is gradual since it is a curved line but I need to approximate it).

To do that I thought about some constrained optimization to be run before the fitting, in order to find those 2 points before computing the function $fun$. The MATLAB function I could use to do that is

```<code>x_opt = fmincon ( fun2, x02, A2, b, Aeq, beq, lb, ub )
</code>```

but I do not know how to express the constraints related to the slope changes.

Is there any other way I could to this?


问:
Frozen coefficient method (von Neumann stability analysis)
Earlier it was considered that frozen coefficients method for Neumann stability analysis for finite difference scheme is more heuristic than rigorous. But I have read some information in a book by Gustafsson,Kreiss,Oliger that there are some conditions, which give us possibility to use frozen coefficients+Neumann analysis at the level of full rigor. As a book didn't give any reference, could you help me to find these conditions?</p>


问:
smallest eigenvalues for linear elasticity
I want to compute a few tens of the smallest eigenvalues of a linear system which is a discretization of a linear elasticity. In the presence of additional constraints like Dirichlet boundary conditions, does it help to use the rigid body displacements to deflate the linear system? In that case they might only be approximations of zero eigenmodes, right? Thank you</p>


问:
Optimizing rank computation for very large sparse matrices
I have a sparse matrix such as

```<code>A =

   (1,1)        1
   (3,1)        1
   (1,2)        1
   (2,2)        1
   (1,3)        1
   (3,3)        1
   (4,3)        1
   (4,4)        1
</code>```

The full matrix of A can see look like as following:

```<code>full(A) =

     1     1     1     0
     0     1     0     0
     1     0     1     0
     0     0     1     1
</code>```

I want to find the rank of matrix A by fast way (because my matrix can extend to $10000 \times 20000$). I try to do it by two ways but it give the different result

Convert to full matrix and find rank using

```<code>rank(full(A)) = 3
</code>```

Find the rank using sprank

```<code>sprank(A) = 4
</code>```

The true answer must be 3 that means using first way. However, it take long time to find the rank,especially matrix with large size. I know the reason why the second way give 4 because sprank only tells you how many rows/columns of your matrix have non-zero elements, while rank is reporting the actual rank of the matrix which indicates how many rows of your matrix are linearly independent. sprank(A) is 4 but rank(A) is only 3 because you can write the third row as a linear combination of the other rows, specifically A(2,:) - A(1,:).

My problem is that how to find the rank of a sparse matrix with lowest time consumption

Update: I tried to use some way. However, it reported larger time consumption comparison with rank function. Could you suggest to me other way?

```<code>%% Create random matrix
 G = sparse(randi(2,1000,1000))-1;
 A=sparse(G) %% Because my input matrix is sparse matrix
%% Measure performance
>> tic; rank(full(A)); toc
Elapsed time is 0.710750 seconds.
>> tic; svds(A); toc
Elapsed time is 1.130674 seconds.
>> tic; eigs(A); toc
Warning: Only 3 of the 6 requested eigenvalues converged.
> In eigs>processEUPDinfo at 1472
  In eigs at 365
Elapsed time is 4.894653 seconds.
</code>```


问:
Eigenvalues and Timestep restriction Follow up
This is a follow-up question to the previous questions I had on eigenvalues. Please let me know if I should edit the previous question itself for asking this.

If the eigenvalues of a matrix generated from a lower order coarse mesh discretization are negative, can we comment anything on the eigenvalues of the matrix that we get by refining the mesh or using a higher order polynomial for discretization of the same pde?


问:
How to tell if symplectic integrator is suitable for my problem?
This question follows my earlier question How can I call the Boost C++ odeint Runge-Kutta integrator for a system of ODEs?

I'm still considering the system:
$$ \ddot x = - \frac{Ax}{||x||^3} $$

My intention was to use a classical Runge-Kutta 45 method to solve ODEs of my system. However, I have seen recommendations for using symplectic integrators, because they do a better job at conserving the total energy of the system.

How do I know if a symplectic integrator is a better choice for my system? And either way, what are the downsides of using symplectic integrators over, say, RK45 ?


问:
How to solve an ode with stochastic time-dependent input
I am trying to repeat an example I found in a paper.

I have to solve this ODE:

$25 \ddot{x} + 15 \dot{x} + 330000 x = p(t)$

where $p(t)$ is a white noise sequence band-limited into the 10-25 Hz range.

Then it says the system is simulated using a Runge-Kutta procedure. The sampling frequency is set to 1000 Hz and a Gaussian white noise is added to the data in such a way that the noise contributes to 5% of the signal r.m.s value (how do I use this last info?).

The main problem is related to the band-limited white noise. I followed the instructions I found here https://dsp.stackexchange.com/questions/9654/how-to-generate-band-limited-gaussian-white-noise and I wrote the following code:

```<code>% White noise excitation
% design FIR filter to filter noise in the range [10-25] Hz
Fs = 1000; % sampling frequency

% Time infos (for white noise)
T = 1/Fs;
tspan = 0:T:4; % I saw from the plots in the paper that the simulation last 4 seconds
tstart = tspan(1);
tend = tspan (end);

b = fir1(64, [10/(Fs/2) 25/(Fs/2)]);
% generate Gaussian (normally-distributed) white noise
noise = randn(length(tspan), 1);
% apply filter to yield bandlimited noise
p = filter(b,1,noise);
</code>```

Now I have to define the function for the ode, but I do not know how to give it the $p$ (white noise):

I tried coding this way

```<code>odefun = @(t,u,m,k1,c,p)[u(2); 1/m*(-k1*u(1)-c*u(2)+p)];
q0 = [0;0]; % Initial conditions
[t,q,te,qe,ie] = ode45(@(t,u)odefun(t,u,m,k1,c,p),tspan,q0(:,1),options);
</code>```

But I get this error message:

@(T,U)ODEFUN(T,U,M,K1,C,NB) returns a vector of length 4002, but the length of initial
conditions vector is 2. The vector returned by @(T,U)ODEFUN(T,U,M,K1,C,NB) and the
initial conditions vector must have the same number of elements.

I already had a look at the following questions but I do not understand where the problem could be.

https://stackoverflow.com/questions/14343563/how-solve-a-system-of-ordinary-differntial-equation-with-time-dependent-paramete?lq=1

https://stackoverflow.com/questions/19732442/solving-an-ode-when-the-function-is-given-as-discrete-values-matlab?lq=1

https://stackoverflow.com/questions/19732442/solving-an-ode-when-the-function-is-given-as-discrete-values-matlab

EDIT

I might have found a way to do it:

```<code>function [y] = p(t)
    b = fir1(64, [10/(Fs/2) 25/(Fs/2)]);
    % generate Gaussian (normally-distributed) white noise
    noise = randn(length(t), 1);
    % apply filter to yield bandlimited noise
    y = filter(b,1,noise);
end

odefun = @(t,u,m,k1,c,p)[u(2); 1/m*(-k1*u(1)-c*u(2)+p(t))];

[t,q,te,qe,ie] = ode45(@(t,u)odefun2(t,u,m,k2,c,@p),tspan,q0(:,1),options);
</code>```

Does it look fine to you?


问:
how to pick a submatrix of a sparse matrix quickly(matlab)
I have a large sparse matrix A and i want to pick the submatrix of it to do block jacobi iterations.

For the blocks, i have get a matrix Q that contains the index of nonzero entries in its Jth column define the nodes of Jth block.

my code is:

```<code>for i=1:J
lg=find(Q(level).agg(:,i));
D(i).matrix=inv((sparse(A(lg,lg))));
end
D_inverse=sparse(blkdiag(D(1:J).matrix));
</code>```

and I found that the operations A(lg,lg) cost much time, because A is sparse format,

Is there any other ways to do this purpose (building D_inverse)quickly and efficently?

Thank you


问:
Checking if convex polytope is nonempty
I am currently running a linear program with MATLAB to determine, by the exitflag of linprog, if two rotated and shifted hypercubes have nonempty intersection. I wondered if this is a waste of time, or if linprog does nothing more than checking feasibility before starting the optimization.</p>

Title: Checking if convex polytope is nonempty


问:
Is it normal to expect the error of simulation of a damped harmonic oscillator to decrease as the damping factor decreases?
I am simulating a damped harmonic oscillator using the RK4 method of numerical integration. I am comparing the simulated results with the analytical ones (for the free evolution case) and obtaining the error as a function of step-size, for different values of the quality factor Q.

$$\ddot{X} + \frac{1}{Q}\dot{X} + X = 0$$
The analytical solution is given by:
$$x(\tau) = e^{-\tau / 2Q} \left[c_1 \cos\left(\tau \sqrt{1 - \frac{1}{4Q^2}}\right) + c_2\sin\left(\tau \sqrt{1 - \frac{1}{4Q^2}}\right)\right]$$
where $c_2$ and $c_2$ can be found by using the boundary conditions: $x(0) = 1$ and $\dot{x}(0) = 0$.

For RK4 method, we obtain two first order equations which we can use in the standard RK4 formula:
Eqn1:
$$\dot{X} = P \Rightarrow \dot{P} = \ddot{X}  $$
Substitute into original equation to obtain:
$$\Rightarrow \dot{P} ~+ \frac{1}{Q}P + X = 0$$
Which leads to Eqn2:
$$\dot{P} = -\left(\frac{P}{Q}+X\right)$$

So these equations allow us to obtain the next point in $X$ and $dX$ for a particular step-size, $h$, using $$y_{n+1} = y_n + h\left(\frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6}\right)$$

Error Investigation

I am interested in the relative error for different step-sizes and Quality factors.

$$Err_{relative} = \frac{\Sigma_i (x_i - X_i)^2}{\Sigma_i x_i^2}$$

We can logically expect the error to increase as the step-size increases, which is what my graphs show. I run this simulation for different values of Q.

I would expect the error to decrease as Q increases, but instead the trajectory is just shifted upwards, i.e: the error is increased but the behaviour with respect to the step size is what we should expect.

It might be hard to see on the graph but those trajectories are for Q = 5 and Q = 95. The trajectory for Q = 5 is lower than the other one. I would have expected it to be the other way round.

Assuming I haven't goofed up my code, are my assumptions correct or is this normal behaviour?


问:
Coding a convex problem in CVX
I am new to CVX and am trying to simulate this convex problem I found in a paper.

$$\min_{\gamma,\mathbf{mu},\mathbf{G},\mathbf{\Omega},t} \text{Tr}(\mathbf{G}\mathbf{C}\mathbf{G}^H)+t \\ s.t. -t-\frac{\gamma}{2}\sum\limits_{i=1}^n\log(1-\frac{2\mu_i}{\gamma})-\gamma \log(\delta)\leq0\\\mu_1\geq\mu_2\geq\ldots\geq\mu_n\geq0\\2\mu_1<\gamma\\ \text{S}_k(C^{1/2} \Lambda C^{1/2})\leq \sum\limits_{i=1}^{k} \mu_i \quad \forall k = 0,1,\dots,n-1\\
\text{Tr}(\Omega C)\leq \sum\limits_{i=1}^{n}\mu_i\\\Omega\geq \Lambda$$

where $$\Lambda = (\mathbf{G}\mathbf{H}-\mathbf{I})^H(\mathbf{G}\mathbf{H}-\mathbf{I})$$ and $S_k(\mathbf{A})$ is the sum of the largest $k$ eigenvalues of $\mathbf{A}$. In the problem, the scalar $\delta$, matrices $\mathbf{H}$ and $\mathbf{C}$ are assumed known. $\mathbf{C}$ is symmetric positive definite. I am finding it difficult to reformulate the first constraint with the log summation.

Any help will be greatly appreciated.

Title: Coding a convex problem in CVX


问:
Can you provide me with the source code of DMRG writen by C++?
I want to study the Density matrix renormalization Group(DMRG). Can some one provide me with the souce code which is writen by C++ for study purpose only? The code is not necessiary very professional, but works.


问:
In-place QR update: deleting a column
Background

I'm trying to do an update to a "thin" QR decomposition ($A = QR$, where $Q$ is $\mathbf{R}^{m,n}$, the first few columns (up to the matrix rank) of an orthogonal matrix and $R$ is upper triangular and square and $\mathbf{R}^{n,n}$).  The matrix itself is a full column rank "tall" matrix ($m \ge n$), and I'm periodically adding and removing columns from $A$ as part of a larger pivoting optimization algorithm.

I'm doing the decomposition in-place inside $A$ using Householder transformations, with the Householder vectors stored in the lower triangular portion of the matrix and $R$ stored in the upper triangular portion of $A$.  I need to store the diagonals for $R$ in another extra array.  So I only need $\mathcal{O}(n)$ more memory beyond the space of the matrix being decomposed itself.

Requirements

At each step of the pivoting algorithm, I need to calculate

$$R^T(Rx - Q^T f)$$
$$Q Q^T f $$

The $Q^T f$ term in both equations can be calculated without needing to store $Q$ directly just by applying each Householder transformation to $f$ as they come in from newly added columns.  That is, at each step $f_{n} := H_n f_{n-1}$ for some Householder transformation $H_n$.  Unfortunately I don't see a way to update $Q Q^T f$ in the same way without explicitly storing $Q$ somewhere, so my guess is I still need to hold on to $Q$. ($Q Q^T f$ is basically the projection of $f$ on to the column space of $A$).

Updating

When deleting a column, all the literature I've found recommends using a series of Givens rotations for the "downstream" columns to push them from Hessenberg form back to triangular form.  However, I don't have anywhere to store these so I can calculate $Q Q^T f$ when I need it.  For the column being removed I need to keep it available to be re-added in the future, so I can't use it for scratch space.  Ideally I'd like to do the decomposition using only $\mathcal{O}(n)$ additional memory beyond the space provided by $A$ in the first place.

My question

Is there some other way to push the downstream columns back to triangular format that wouldn't use more memory?  I've tried to think of a way to update the Householder transformations, but I haven't been able to figure out a way that isn't $\mathcal{O}(n^3)$.  I'd like the update itself to be $\mathcal{O}(n^2)$ as it is with the Givens rotations.

Alternatively, is there a way of updating $Q Q^T f$ at each stage in the pivoting algorithm when a column is added or removed from $A$ without explicitly storing $Q$?

Title: In-place QR update: deleting a column


问:
I am looking for a complex sparse matrix EigenVector solver for GPGPU; preferably CUDA
So far the closest I've found is ViennaCL, which has a Lanczos implementation for Eigenvalues.

It is not clear that EigenVectors are produced by this library.

Does anyone here know whether ViennaCL solves for EigenVectors? And how to get them?

Or of any other sparse matrix EigenVector solver for GPGPU?

ETA:From what I have read, Lanczos is only an algorithm for converting from a general matrix to tridiagonal form. The QR algorithm (or other) is needed for extracting EigenValues/Vectors from the tridiagonal matrix.

Assuming I can upgrade the ViennaCL library to handle complex matricies, does anyone have a recommendation for a GPU based QR algorithm implementation?


问:
ARPACK gives different answers from Matlab and NAG
I'm playing with ARPACK. I looked into the examples they provide, zndrv4.f illustrating the usage of the routine znaupd, in the directory of ARPACK/EXAMPLES/COMPLEX/. I also came cross NAG Fortran Library. In NAG, there are some linear problem solvers F12***. The F12*** routines in NAG are equivalent to the znaupd in ARPACK. So I want to check if they will yield the same results.

I first looked into the example provided in the user guide of F12ARF at http://www.nag.co.uk/numeric/fl/nagdoc_fl22/pdf/F12/f12arf.pdf for example. In the end, it yields the results of

```<code>509.9390
380.9092
659.1558
271.9412
</code>```

around the shift=500. I solved the same generalized eigenvalue problem in Matlab. Matlab gave me the same results.

But when I used znaupd from ARPACK to solve the same problem, I obtained different answers. The 4 eigenvalues are now

```<code>Eigenvalue               Residual
501.65650188259684       2.01604D-01
480.15153312181440       2.10853D-01
526.52596256924164       1.91850D-01
461.99019999608828       2.19336D-01
</code>```

The routines in NAG and ARPACK both use SHIFTED INVERSE mode and solve the same generalized problem.
I'm not sure what was wrong, but you can see the residual is too big (usually it should be around 1e-15).
I attached my scripts of the ARPACK zndrv4.f (it's basically the same as the example file provided by ARPACK, I just need to change a little bit the matrices around line 174 to be the same as that in NAG.) and the Matlab file zndrv4.m in the compressed file below.

https://www.dropbox.com/s/mcq2jkmrgk2839u/zndrv4.zip?dl=0

Then later I found that actually the M matrix in F12ARF is normalized (entries divided by a six). I followed this and got the above wrong results in ARPACK. Now if I didn't divide the entries by six in zndrv4.f, the ARPACK script gives me the correct answer (the same as Matlab).

```<code>Eigenvalue               Residual
4.90161D+02              5.20906D-16
5.50553D+02              8.16401D-16
4.33527D+02              2.29794D-15
6.14759D+02              4.93630D-15
</code>```

I'm confused by this. It seems to say that the ARPACK routine is not robust. Just scaling the matrix with a constant 6 will lead to such a big change in residual? Thank you very much.


问:
Finding a scalar field in order to generate a solenoidal vector from a given vector
I am working on generating a (complex) solenoidal vector field $\mathbf{A}$ from a prescribed (complex) vector $\mathbf{a}$ and the gradient of a scalar, $b$, such that
$$\mathbf{A} = \mathbf{a} + \nabla b.$$
The condition that I require is then
$$\nabla\cdot\mathbf{A} = \nabla\cdot\mathbf{a} + \nabla^2 b = 0 \qquad \Longrightarrow \qquad \nabla^2 b = - B,$$
where $B = \nabla\cdot\mathbf{a}$.

To complicate this slightly, I want to decompose $\mathbf{A},\mathbf{a}, B$ and $b$ into
$$\mathbf{A} = \hat{\mathbf{A}} \exp[i(\omega t - \mathbf{k}\cdot\mathbf{x})] = \hat{\mathbf{A}} \exp[i\phi],\\
\mathbf{a} = \hat{\mathbf{a}} \exp[i\phi],\\
B = \hat{B} \exp[i\phi],\\
b = \hat{b} \exp[i\phi].$$

With this decomposition, $B = (\nabla\cdot\hat{\mathbf{a}} - i\mathbf{k}\cdot\hat{\mathbf{a}}) \exp[i\phi]$, and the Poisson equation above becomes
$$(\dagger)\qquad\nabla^2 b = (\nabla^2 \hat{b} - 2i\mathbf{k}\cdot\nabla b - \lvert\mathbf{k}\rvert^2 b) \exp[i\phi] = -\hat{B} \exp[i\phi].$$

This is then solved for $\hat{b}$ iteratively using the following solver scheme (given for 1-dimension, but it is easily extended to higher dimensions):
$$(\ddagger)\qquad \hat{b}_j = \frac{\Delta x^2 \hat{B}_j + \hat{b}_{j+1} + \hat{b}_{j-1} - ik\Delta x(\hat{b}_{j+1} - \hat{b}_{j-1})}{2 + k^2 \Delta x^2},$$
where $j \in \{2,\ldots,N-1\}$ indicates the discrete grid location, the position $x_j = j \Delta x$, and $\Delta x$ is the grid spacing. The boundary $j=1,N$ is fixed at zero. $\hat{B}_j$ is evaluated as $\hat{B}_j = (\hat{a}_{j+1} - \hat{a}_{j-1})/2\Delta x -ik\hat{a}_j$.

The numerical solution to $\hat{b}$ produced satisfies $$\nabla^2 \hat{b} - 2i\mathbf{k}\cdot\nabla b - \lvert\mathbf{k}\rvert^2 b = -\hat{B}$$ very well, however if I construct
$$\hat{A}_j = \hat{a}_j + \frac{\hat{b}_{j+1} - \hat{b}_{j-1}}{2\Delta{x}} - ik\hat{b}_j$$
and evaluate
$$\nabla \cdot \mathbf{A} \simeq \Big(\frac{\hat{A}_{j+1} - \hat{A}_{j-1}}{2\Delta x} - ik \hat{A}_j\Big) \exp[i\phi],$$
then this is not only non-zero but also large. Putting the expression for $\hat{A}_j$ in terms of $\hat{a}$ and $\hat{b}$, one finds that the above equation reads
$$(*)\qquad \nabla \cdot \mathbf{A} \simeq \Big( \frac{\hat{a}_{j+1} - \hat{a}_{j-1}}{2\Delta x} - ik \hat{a}_j + \frac{\hat{b}_{j+2} - 2\hat{b}_{j} + \hat{b}_{j-2}}{4\Delta x^2} - 2ik\frac{\hat{b}_{j+1} - \hat{b}_{j-1}}{2\Delta x} - k^2\hat{b}_j \Big) \exp[i\phi]\\
= \Big( \hat{B} + \nabla_2^2\hat{b} - 2ik \frac{\partial\hat{b}}{\partial x} - 4k^2 \hat{b} \Big) \exp[i\phi].$$
This is just the 1-dimensional version of $(\dagger)$, but with the Laplacian operator evaluated as
$$\nabla_2^2 \hat{b} = \frac{\hat{b}_{j+2} - 2\hat{b}_j + \hat{b}_{j-2}}{4\Delta x^2},$$ suggesting that it is this difference which is responsible for $\nabla \cdot \mathbf{A} \neq 0$.

To overcome this, I attempted to invert ($*$) rather than use ($\ddagger$), this time with boundary conditions $\hat{b}_j = 0$ for $j = 1, 2, N-1, N$, but the (iterative) solution to the Poisson equation diverges.

I was wondering what is responsible for this? Could my setting the two boundary cells to zero cause such a problem? The solver scheme for point $j$ involves $j+2, j+1, j-1, j-2$. Given the non-solenoidal vector $\mathbf{a}$, is there a better way to find a scalar $b$ such that $\nabla\cdot(\mathbf{a} + \nabla b) = 0$?

Thank you for any help or suggestions!


问:
How big a matrix can we row reduce in reasonable time?
I have very large matrices that I would like to row reduce (I need to keep track of the steps and find a basis of the kernel/image, not just find the rank).

The good news is that I work mod 2 and the matrices are very sparse. The problem is their size, right now I have a beast that is roughly 3 millions x 1 million, but I expect much bigger, probably billions. I might have a thousand of those matrices that I would like to reduce eventually.

I found the linbox project that might be helpful, but before I wanted to know if anybody has a feeling whether or not this is realistic ? I could be using the Grid of my university which has "the combined processing power of 4,568 cores: 1,346 Intel cores, 3,222 AMD cores, with over 13.5TB of RAM and 1.2PB of disk space.". What do you think ?

Thanks!


问:
Tracking the speed of 2D oscillations on a lattice
I wrote a Monte Carlo simulation of the 2D Lotka-Volterra model on a discrete lattice (with periodic boundary conditions). A video that I produced (which images the system after some number of monte carlo steps) can be found here. I am wondering, what is the best way to track the spatial "speed" of the oscillations of the two species across the lattice from one video frame to the next?

My model is done with the following rates, where $A$ denotes a predator node on the lattice, $B$ denotes a prey node on the lattice, and $\emptyset$ denotes an empty space on the lattice:

$
A + B \xrightarrow{\lambda} A+A,\;(predator\;consuming\; prey)\\
B + \emptyset \xrightarrow{\sigma} B+B,\;(prey\;growing\;into\;empty\;adjacent\;spaces)\\
A \xrightarrow{\mu} \emptyset,\;(predators\;dying\;off)\\
$

For each Monte Carlo step I choose one space on the lattice randomly, do a cumsum(each neighbor's rate./sum(neighbors' rates)), and then generate a random number to decide which of the three reactions above will be done. In my video, $\lambda=1.6$, $\sigma=50$, and $\mu=0.1$.

I am not sure if I should try to track a certain shape in each frame and measure its deformation, or if I should try to use some other sort of spatial metric. To be clear, the data for each video frame is stored in an N by N matrix (with values 0 for empty, 1 for predator, and 2 for prey), so I am not limited to an analysis of the video I provided, it is only for explanation. My source code for this project can be found here.

What is the typical way of approaching this sort of problem, perhaps someone can provide a reference?

Thanks.


问:
Data format for Cartesian mesh with cut cells
Is there a mesh data format which is particularly well suited for Cartesian meshes with Dey-Mittra style cut-cells? Ideally, one which does not require converting the cut-cells into a more general polyhedron.

Non-proprietary formats with open-source implementations are also preferable.


问:
What is a reasonable manufactured solution to test finite difference method?
What is a reasonable manufactured solution to test the following equation against its finite difference approximation? I want it to look like a Cosine function about $0$, rotated about $Z$ axis (basically not varying about $\theta$) such that $f$  turns out to be a not so bad function in terms of $r,\theta$

$$u_{rr} \;+\; \frac{1}{r} u_r \;+\; \frac{1}{r^2} u_{\theta\theta}\;+\; Cu \;=\; f(r,\theta)$$

Also, it should avoid the discontinuity at $r = 0$


问:
Time-stable spectral decomposition algorithm
Consider an $n \times n$ real, time-dependent matrix $A(t)$ such that $A(t) = A(t)^{T} > 0$ and $A(t)$ is continuous on $[a,b]$. Then it is posible to specify a matrix $S(t) \in SO(n)$ such that $S(t)A(t)S^{T}(t) = \Lambda(t)$, where $\Lambda$ is a diagonal matrix. For instance, we can obtain such $S(t)$ using Jacobi rotation method.

I'm looking for time-continuous on $[a,b]$ matrices $S(t)$ and $\Lambda(t)$. Lets consider a net
$$
  a = t_0 < t_1 <  \ldots < t_n = b
$$
Assume that $S(t_k)$, $\Lambda(t_k)$ are computed. How to receive than $S(t_{k+1})$ and $\Lambda(t_{k+1})$ nearest possible to $S(t_{k})$ and $\Lambda(t_k)$? Which algorithm should I use?


问:
infinite but non-periodic space with PMLs in Comsol
I have a problem with implementing proper boundary conditions on sides of my simulation (in RF or WaveOptics module of Comsol). I want to obtain an infinite but not periodic space (beacuse then I will put something inside). For example, it can be a glass-vacuum interface with incident plain waves. For such a simulation one can check if Fresnel equations work as they should and if not - the simulation is not properly set. I was trying with PMLs around the simulation area and in/out ports from the bottom/top side. But it doesn't work at all. Should I set another type of boundaries on the sides except PMLs? </p>


问:
Finite Difference in Polar Co-ordinates
Background

Please note that I am duplicating the question on scicomp. I have already asked this in math.

I am trying to come up with a scheme in Polar Co-ordinates for the following PDE:

PDE I am trying to solve
$$  u_t - (u_{rr} + \frac{1}{r} * u_r + \frac{1}{\theta} * u_{\theta\theta} + bu) = f(r,\theta, t)$$

I am extending the ideas from the book [1] which is already given for the PDE:

PDE in [1]

$$u_{rr} + \frac{1}{r} * u_r + \frac{1}{\theta} * u_{\theta\theta} = f(r,\theta, t)$$

The book discusses a particular scheme in pages 333-334 to solve the above PDE. I am listing the scheme below for convenience and its notation.

Scheme proposed for the above PDE in [1]

$$\frac{1}{r_i} (r_{i+1/2} \frac{u_{i+1j} - u_{ij}}{\Delta r} - r_{i-1/2}\frac{u_{ij} - u_{i-1j}}{\Delta r}) \frac{1}{\Delta r} + \frac{1}{r_i^2} \frac{u_{ij+1} - 2u_{ij} + u_{ij-1}}{\Delta \theta^2} = f_{ij}$$

where $u_{ij}$ and $f_{ij}$ are the grid functions at $(r_i, \theta_j) = (i\Delta r, j\Delta \theta)$

Further, to avoid difficulties arising at the origin because of the $\frac{1}{r}$ terms in the scheme, the book goes on to integrate the PDE over a small disc (D) of radius $\epsilon$ near the origin $r=0$, applies the Gauss-Divergence Theorem to get the following result:

$$\iint_D f \;r\;dr\;d\theta =  \iint_D (u_{rr} + \frac{1}{r} * u_r + \frac{1}{\theta} * u_{\theta\theta}) \; \;r\;dr\;d\theta = \iint_D \nabla. (\nabla u) \; \;r\;dr\;d\theta  = \int_{0}^{2\pi} \frac{\partial u}{\partial r} \;\epsilon\;d\theta$$

These are later approximated as follows:

$$f(0) * \epsilon^2 * \pi = \sum_{j=1}^{J} \frac{u_{1j}-u_0}{\Delta r}  \epsilon\; \Delta \theta$$

From which it solves for $u_0$.

Note: $J\Delta \theta = 2 \pi$ and $I\Delta r = R$ where R: radius of the domain.

Question

i) How do I extend the idea for the terms $u_t$, $bu$ terms above when computing an approximation for $u_0$ for my PDE?

I propose to treat $bu$ term where $b$ is a constant the same way as $f$ is treated in the approximation above. But it is unclear to me as to what I need to do for $u_t$.

Any help is greatly appreciated. Even pointers to read a particular book/topic.

Thanks!

References:
[1] Elliptic Partial Differential Equations and Difference Schemes
Finite Difference Schemes and Partial Differential Equations, Second Edition. 2004, 311-338

Title: Finite Difference in Polar Co-ordinates


问:
Speeding up the classical Jacobi method using Scheduled Relaxation method?
There has been quite a flutter recently in the iterative world about an algorithm that speeds up the classical Jacobi method by as much as 200 times using a scheduled relaxation method where a combination of over and under relaxation is employed. The results as shown in this paper are remarkable and are indicative of a resuscitation of Jacobi method as a method for practical purposes. </p>

<p>However, when I implemented it I cannot get it to work. I was wondering if anybody else in this community has tried implementing this method and if so was the attempt successful. It would be interesting to hear if there is anything special that needs to be taken care of when implementing this algorithm. </p>


问:
2nd order differences for modelling a physically discrete phenomenon?
Consider the diffusion of some unknown $R$ on a physically discrete, 1-D ring:

```<code>    x---------x---------x---------x
n = 0         1         2         3
</code>```

We let grid point $n=0$ be the same as grid point $n=3$, giving us our ring. We are given a diffusion coefficient $D$, which is still multiplied by some to some notion of the "Laplacian" (which in 1-D, is just the second order difference) of $R$ in order to tell us how $R$ moves diffuses between neighbouring nodes

We can use techniques such as 2nd order central differences to calculate some notion of second order difference of $R$ on these physically discrete points, if we consider the numerical discretization of the grid to be the same as the physical discretization.

What if the physical discretization is not uniform?

```<code>    x----------------x--x-------x
n = 0                1  2       3
</code>```

What notion of "second order difference" in $R$ might I now consider given that:

I can't simply consider the physical discretization to be the same as the numerical discretization, without complicating things (non-uniform numerical grid spacing);

but more importantly, I don't have any continuous physical situation to go back in order to derive a different second order differencing scheme using Taylor series, so I can't use some finer uniform numerical grid, which also happens to coincide with the physical points -- it would make no sense to have a numerical grid point that doesn't coincide with a physical grid point.

Is my only solution to make some sort of assumption regarding the "continuous" nature of $R$ between the grid points? Maybe I assume that $R$ varies linearly between the physically discrete points? Then the problem I outline in point 2. isn't there, and I can define a finer numerical grid? This is the solution advocated in a similar question here: How can I numerically differentiate an unevenly sampled function?

What other ideas might I consider?


问:
Absorbing BC's / PML on a graph
The wave equation,
$$\ddot{u} = c^2 \Delta u,$$
 can be generalized to abstract graphs by using the negative graph Laplacian in place of the physical Laplacian.

Is there a graph-theoretic analog of Perfectly Matched Layers (PML), where "outward traveling waves" exit the graph without reflection?

It is not clear to me what the definition of "outward traveling waves" should be in the graph context, but it does seem intuitively reasonable to consider nonresonant waves propagating through a small region of interest (subgraph) embedded in a much larger sparse graph.

Title: Absorbing BC's / PML on a graph


问:
Stationary 2D/3D Navier-Stokes source code
Trying to solve stationary Navier-Stokes problem for incompressible laminar Newtonian fluid.

I've found a couple solutions for instationary Navier-Stokes equations (like FeniCS examples or CFD Python), but nothing on stationary.

As far as I understand methods are totally so this code couldn't be reused properly. Does anybody know any software which may help?


问:
Find the solution of linear equation using Wiedemann/ Krylov method

Let given

$M =$

```<code> 1     0     1
 0     1     1
 1     1     1
</code>```

and
$b =$

```<code> 1
 0
 1
</code>```

How to find the solution $x_3$ where $x=${$x_1,x_2,x_3$}?

Solution: Based on Wiedemann algorithm we has $u_i=[0; 0; 1]$.
I am getting confused in the Step2. How to find $M_s$?
Is it just the $3^th$ columns of $M_t$ as $M_s=[1 1 1]$ or it is created from Krylov sequence? Please help me find M_s and minimal polynomial?

Matlab code:

```<code>%% Input Generation Matrix M, b and index ith
%% Output single symbol x
M=[1 0 1;0 1 1;1 1 1];
b=[1; 0; 1];
index=2;
[K N]=size(M);
Mt=M';
Imatrix=eye(K,K);
x=[];
for index=1:K
u_i=Imatrix(:,index);
se_Krylov=[];
se_Krylov(:,1)=u_i;
    for i=2:(2*K) % due to i=1 is u_i
       se_Krylov(:,i)=mod((Mt^i)*u_i,2);
    end
%% M_s is the operator Mt retricted to S
M_s=se_Krylov(index,:);

%% Compute the polynominal using Berlekmap Massey
[f, LCP] = Berlekamp_Massey2(M_s);
d=size(f,2)-1; %deg of f: x^d+x^(d-1)+....

x_comma=zeros(K,1);
    for i=d:-1:1
        x_comma=mod(x_comma+f(d).*(Mt)^(i-1)*u_i,2);
    end
x_single=mod(x_comma'*b,2) ;%Inner product
x(index)=x_single;
end
</code>```


问:
Fast algorithms to solve Markov Decision Processes
In my master thesis I used an Algorithm called Approximative Dynamic Programming [1] to solve equations of the form

$$
\max_{\pi}\mathbb{E}^{\pi}\left\{\sum_{t=0}^{T}\gamma^tC_t^{\pi}(S_t,A_t^{\pi}(S_t))\right\}.
$$

It uses Monte-Carlo sampling and an approximation $\overline{V}_t$ of the value function to come around the curse of dimensionality. As a decision function serves a convex linear program

$$
\hat{v}_t^n=\max_{a_t\in A_t^n}\left(C_t(S_t^n,a_t)+\gamma \overline{V}_{t+1}^{n-1}(S^M(S_t^n,a_t))\right),
$$

and the update is done by using a stepsize $\alpha$

$$
\overline{V}_t^n(S_t^n)=(1-\alpha_{n-1})\overline{V}_t^{n-1}(S_t^n)+\alpha_{n-1}\hat{v}_t^n.
$$
I wondered now if they exists similar algorithms, despite the many flavours of this algorithm, to deliver high quality solutions in reasonable time or if the used principle is the only one to solve such kind of problems?

Cheers,
Reza

[1] POWELL, W. : Approximate Dynamic Programming: Solving the curses of dimensionality. Bd. 703. Wiley-Blackwell, 2007


问:
Is it possible to predict the null space of a structure from contributing elements null spaces?
I am trying to solve an almost incompressible problem with heterogeneous properties by domain decomposition. Solution with CG converges slowly or divergerces completely. My problem becomes ill-conditioned (condition number in the order of 10^5) due to the poissons ratio approaching to 0.5 and stiffness differences in between materials. I detected that in element stiffness matrix, several additional zero energy modes (hourglass modes etc) are formed in addition to standard rigid body modes.

I am planning to prevent these modes by MatNullSpace object in Petsc. Since the computation of nullspaces with SVD or QR is expensive, I am trying to find a cheaper way of predicting the nullspace of floating subdomains.

Is it possible to predict nullspaces for subdomain from contributing finite element's null spaces?

I read a broad sentence in a publication says that "The rigid body modes of a collection of elements is equal to the assembly of the rigid body modes of the individual elements taking into account the multiplicity of those degrees of freedom that lie in multiple neighboring elements."

does that really make sense?


问:
Understanding how Numpy does SVD

Possible Duplicate:

  Understanding how Numpy does SVD

I have been using different methods to calculate both the rank of a matrix and the solution of a matrix system of equations. I came across the function linalg.svd, which I now learned is in fact a LAPACK routine. Comparing this to my own effort of solving the system with Gaussian Elimination, it appears to be both faster and more precise. I'm trying to understand why.

As far as I know, the linalg.svd function uses a QR algorithm to calculate the eigenvalues of my matrix. I know how this works mathematically, but I don't know how Numpy manages to do it so quickly and without losing much precision.

I'm obviously not saying my script is in any way close to what the LAPACK authors do, but I'm trying to understand the fundamentals behind it.

So my question: How does the numpy.svd function work, and more specifically, how does it manage to do it fast and accurately (compared to gaussian elimination)?

Thank you for your time.

Title: Understanding how Numpy does SVD


问:
Arbitrary Choosing of the Solution Domain - Navier Stokes and Manufactured Solutions
I want to verify a finite-volume solver (SIMPLE-Algorithm) for the incompressible Navier-Stokes equations by using a manufactured solution. I use Dirichlet boundary conditions for the velocity at all boundaries. The manufactured solution for the velocities I use, is constructed as the curl of a vector field and thereby fulfills the continuity equation at any point.

Depending on how I choose the domain on which to solve the equations I run into convergence problems: The residual of the pressure correction equation stagnates after an initial decrease.

I am approximating the mass fluxes through boundary faces using the midpoint rule and apparently this leads to a net loss/gain of mass considering the entire problem domain.

Is this the normal behavior? Are there any remedies other than choosing a manufactured solution that identically vanishes at the boundaries? Would it be an option to provide the exact mass flux instead of the approximate (not sure if this is allowed, since the Dirichlet boundary condition only fixes the velocity)?


问:
Reference request for numerical variational method
I have a variational problem where the unknown function is a periodic path $\gamma:[0,1)\to\mathbb{R}^2$, and the functional is
$$ \int_0^1\left( \tfrac12\|\dot\gamma(s)\|^2 + \mathcal{F}[\gamma]\right)\,ds. $$
Here $\mathcal F$ depends on values of $\gamma$ at different values of $s$, so it's not just a function of $\gamma(s)$, like so:
$$ \mathcal{F}[\gamma](s) = \sum_{0<j<n}\frac{1}{\|\gamma(s+j/n)-\gamma(s)\|}. $$

My current approach is to expand $\gamma$ into its Fourier components $e^{2\pi \mathrm{i} k s}$ to ensure periodicity, and use a Newton-type method to minimize the functional as a function of the Fourier coefficients. I can calculate the functional and its gradient efficiently for any given set of Fourier coefficients, in part because the integral can be done with the trapezoid rule, which converges quickly for periodic functions.

One of the problems I have is that it is difficult to explore all the different minima the functional might have, because apart from some linear constraints I can impose on the Fourier coefficients, I have very little control over which local minimum I get.

Is there a good graduate/research-level book on numerical variational methods that would explain the different things I should try? Survey articles?

Edit: the exact functional I'm using comes from the n-body choreography problem, but I'm interested more in the general ideas here than in this specific application.


问:
Perron-Frobenius theorem on general real symmetric matrices
From the Perron-Frobenius theorem, it might be concluded that the spectral radius is the largest eigenvalue for positive matrices, ie, matrices with strictly positive entries. In other words, the largest magnitude eigenvalue of a symmetric, real, positive matrix, is a positive eigenvalue.

However, I'm interested if the theorem is applicable to an arbitrary real symmetric matrix $R$ with no zero entries. Namely, the particular interest is whether the [power iteration] method2, might be used to extract the eigenvector and the associated largest positive eigenvalue of an input matrix.

Obviously, I'm concerned with the possibility of largest magnitude eigenvalue being negative. Is the Power Iteration inappropriate in this case? If not in general, does stating that the eigenvector of all ones, $1_n$, corresponding to the eigenvalue of $0$ spans the null-space of $R$?


问:
Get a single value from a global Vec in petsc
I don't understand the concept with local and global vec's in petsc when it come to DMDAVecGetArray and VecGetValues. I'm creating some DMDA Vec in parallel which works well. Now I want to access certain non local values of my Vec without receiving the whole global array (will not fit into my main memory) on each processor.

Is there something I'm getting wrong? In the documentation of VecGetValues it is written "Gets values from certain locations of a vector. Currently can only get values on the same processor".

thx


问:
Search Direction in Conjugate Gradient
Could you help me with a Conjugate Gradient question? In using CG to solve Ax=b, why is the search direction $p_{k+1}$ in CG chosen as a linear combination of the residual $r_k$ and previous direction $p_k$? I know the set of search directions are supposed to be conjugate with respect to A. One way to generate such sets is using the eigenvectors of A, but this is very computationally expensive. Another way is to use modified Gram Schmidt, but this requires so much storage. Hence, we choose $p_{k+1}=-r_k+\beta p_k$... but why this choice?</p>

Title: Search Direction in Conjugate Gradient


问:
Thin plate stiffness: analytical formula to validate FEM model
I tried to compute analitically the stiffness of a cantilever thin plate (shown in picture). The plate is also homogeneous and isotropic. The aim is to compare the result I obtain with the result I get when running the FEM analysis.

$a = 0.127; b = 0.0508; h = 0.0012; E = 7*10^{10}; \nu = 0.3$

This is my theoretical derivation:

$D \nabla^4 w(x,y,t) + \rho h \ddot{w}(x,y,t) = 0$

where $w$ is the coordinate in the $z$ direction and $D = \dfrac{E h^3}{12(1-\nu^2)}$ is the flexural rigidity, $\rho$ is the density and $h$ the thickness.

$w(x,y,t) = \vec{N}(x,y) \vec{u}(t)$ approximation with a trial set

I am interested in the first bending mode only, therefore the trial set $\vec{N}(x,y)$ has only one element:

$N_i(x,y) = X_r(x)Y_s(y) = X_1(x)Y_1(y)$

The boundary conditions $\vec{N}(x,y)$ has to satisfy are

$w(0,y) = 0$ constraint on the translation

$w_{/x}(0,y) = 0$ constraint on the rotation ($_{/x}$ means derivation with respect to $x$)

Therefore $X_1(x) = \biggl(\dfrac{x}{a} \biggr)^2$ and $Y_1(y) = \biggl(\dfrac{y}{b} \biggr)^0 = 1$

The stiffnes $K$ can be found solving:

$K = \int\int_{A} D [N_{/xx}(N_{/xx}+\nu N_{/yy})+N_{/yy}(\nu N_{/xx}+N_{/yy})+2 (1-\nu)N_{/xy}N_{/xy}]dA $

The only terms different from $0$ are $N_{/x}$ and $N_{/xx}$. After the integration on the area $A$ of the plate I obtain

$K = \dfrac{4Db}{a^3} = 1098.8$

Then I run the FEM analysis and I post-process the results: knowing the time history of the acceleration and the input force applied to the plate ($p = \sin(\omega t)$ applied at the bottom tip corner of the plate, $\omega = 1^{st}$ natural frequency of the plate) I can compute:

$f_i(x,\dot{x}) = p_i-m\ddot{x}_i$

Then I can extract the curve $f(x)$ setting $\dot{x} = 0$. Fitting this curve with a line and computing the slope of this line I get $K = 3618.0$ which is different from the analytical value.

What am I missing/doing wrong?


问:
Convergence of a DASPK depending on DAE formulation
I have a system of non-linear DAE and I noticed that the system does not converge if some of the equations are not differentiated. For example, if the control volume equation is represented as this:
$V_{total}  = V_1 + V_2$
The solver won't convege.

However if I differentiate the equation like this:

$ 0 = m_1(\partial v_1/\partial p)(dp/dt)+\dot{m}_2v_2+ m_2(\partial v_2/\partial p)(dp/dt)+(\dot{m}_2v_2) $

Then it does converge.
I'm using DASPK on octave.

Background:
The equation is part of another five equations that are used to solve the dynamics of a steam accumulator. In such a system, a metal vessel is almost half filled with liquid water and steam is occupying the top other half of the vessel. The system operates under varying pressures and the two phases (1: liquid and 2: steam) are not always in equilibrium.

The full system of equations is as follows:
the Unknowns are
$M_{1},M_{2},M_{1r},M_{2r},p,\dot{m}_r$
But notice that $\dot{m}_r$ is actually calculated readily from 4

$\frac{d}{dt}[M_1+M_2] = \dot{m}_f - \dot{m}_s$

$\frac{d}{dt}[M_1h_1+M_2h_2-pV_t] = \dot{m}_fh_f - \dot{m}_sh_s$

$V_{t}  = V_1 + V_2$

$\frac{d}{dt}[M_{1r}+M_{2r}] = \dot{m}_{dc} - \dot{m}_r$

$\frac{d}{dt}[M_{1r}h_1+M_{1r}h_2-pV_r] = \dot{m}_{dc}h_1 - \dot{m}_sh_s\alpha_r$

$V_{r}  = V_{1r} + V_{2r}$


问:
How to determine the minimum number of multiplication needed for a specific expression?
Is there any algorithm to determine the minimum number of multiplication(division) of a specific expression? and the optimal expression form for implementation?

For example, given values of $\cos\alpha, \sin\alpha$, $\cos\beta, \sin\beta$ and $\cos\theta,\sin\theta$, what is the minimum number of multiplication in order to obtain numerical matrix:

How to obtain the final optimal expression (codes) for implementation?


问:
Library for calculating determinants with Kronecker products
I need to calculate a determinant consisting of vectors, using the Kronecker product as product. As an example I would need to be able to calculate:

$\left| \begin{array}{cc}
\left(\begin{array}{c}1\\0\end{array}\right) & 0 \\
0 & \left(\begin{array}{c}0\\1\end{array}\right) \end{array} \right|=\left(\begin{array}{c}0\\1\\0\\0\end{array}\right)$

I know of libraries containing routines for tensor network calculations, but those libraries to my knowledge do not implement determinants. Doing this right will require a good amount of work, so I wanted to ask beforehand whether anybody knows of any code that implements this already.

This is pretty specialized and therefore a long shot, but I wanted to ask nonetheless before spending days implementing it myself. Any language would be acceptable (be it Mathematica or Fortran77 or anything else).

Thanks


问:
3D animation of a vector field (ParaView)
I am trying to create an animation of a 3d vector field in ParaView.

All my data is in three matrices (201 x 201 x 201) corresponding to vectors components in x,y,z directions. And it's calculated in MATLAB.
If I do:

```<code>B = cat(4,Bx,By,Bz);
B = permute(B,[4 1 2 3]);
fid = fopen('Bvec.raw','w'); fwrite(fid,B,'float'); fclose(fid);
</code>```

I can load Bvec.raw into ParaView and visualise the vector field at one instant of time.

Now, how would I go about making an animation if I had say 100 (Bx, By, Bz) triples at different instants of time? Without creating 100 different raw files and loading them into paraview by hand.

ALSO, I'm not set on using paraview, it's just the first software I tried, so if you think there's an easier/better/ more competent way of going about it I'd appreciate suggestions.