问:
Python package to calculate static force and moments of rigid body
Which Python package is suited for solving problems of the following type?

Given the rigid body depicted in violet in the following sketch

I would like to do the following:

Define cartesian coordinate systems $\text{cs0}(x_0, y_0, z_0)$, $\text{cs1}(x_1, y_1, z_1)$ and $\text{cs2}(x_2, y_2, z_2)$.

Specify forces and moments acting on the center of coordinate systems.

Calculate forces and moments acting on fixed clamp at origin of $\text{cs0}$.

In order to preserve flexibility and intuitive usage, object orientated syntax would be preferred.

Pseudocode:

```<code>...
cs0 = initCartesianCoordinatesystem()
cs1 = cartesianCoordinatesystem( cs0, tranlationX = 0, translationY = 20,\
    translationZ = 0, rotationX = 0, rotationY = 0, rotationZ = 0)
cs2 = cartesianCoordinatesystem( cs1, translationX = 10, translationY = 0,\
    translationZ = 0, rotationX = 0, rotationY = 135, rotationZ = 0)

body = rigidBody()
body.addForce( cs1, (1, 0, 0) )
body.addMoment( cs1, (0, 0, 1) )
body.addMoment( cs2, (0, 2, 0) )
force, moment = body.staticEquilibrium( cs0 )
</code>```

Leading to
$$ force = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}_{\text{cs0}} $$
$$ moment = \begin{pmatrix} 0 \\ 2+20 \\ 1 \end{pmatrix}_{\text{cs0}} $$

The following features would be helpful:

Visualization

Units


问:
Name for vectors in a Krylov space but not the preceding one
It seems to me that a useful concept to define when studying Krylov subspace methods is the idea of a vector that belongs to a Krylov subspace $\mathcal{K}_{n+1}(A,b)$ but not to the preceding one $\mathcal{K}_{n}(A,b)$, i.e., they can be written as a polynomial $p(A)b$ of degree exactly $n$, with nonsingular leading term.

Is there a standard name or notation for vectors that belong to $\mathcal{K}_{n+1}(A,b) \setminus \mathcal{K}_{n}(A,b)$? Or, how can I describe or denote them? Something like "$v$ is a proper vector in $\mathcal{K}_{n+1}(A,b)$" sounds good, but I don't want to invent a new term when there already exists a perfectly cromulent word for it1.

1: sorry, I had to.

ContentLicense: CC BY-SA 3.0


问:
Eigenvalue ODE in Spherical Coordinates--Numerical
I wish to solve an eigenvalue problem:

$$\nabla^{2}f=Ef $$

If I assume spherical symmetry $f(r,\theta,\phi)=f(r)$, I can reduce the problem to 1D:

$$(\frac{2}{r}\frac{d}{dr}+\frac{d^{2}}{dr^{2}})f=Ef$$

My boundary conditions are $f(r=0)=0$ and $f(r=L)=0$, where $L$ is the size of the discrete grid.

I am attempting to solve this problem using a orthogonal basis set expansion, but the issue I am experiencing is the divergence of $r$ at $r=0$. Is there a better method for solving this type of problem? How can I accommodate the singularity?


问:
Pros and cons of optimization vs. variational calculus, re: nonlinear elasticity
I'm trying to solve a problem of finding the displacements of an elastic material subjected to external forces.  Those external forces are themselves a nonlinear function of the material displacements.

The first step is to define a functional for the energy in terms of the displacements.  The solution is then the displacement which minimizes that functional over sufficiently smooth functions.

From this functional, Euler-Lagrange equations can be derived, which must be satisfied by the solution.  Those equations are a system of elliptic PDEs.

To solve this problem in practice, it is necessary to discretize.  This can be done directly for the functional, so that the integral of the Lagrangian over the spatial domain becomes a sum over a spatial grid, with the appropriate translation of spatial derivatives to finite differences.  The problem is to minimize the sum, which is generally non-convex.  The alternative is to discretize the Euler-Lagrange equations, which results in a large nonlinear system of equations.  The functional contains terms which are fourth order polynomials of the displacements, while the Euler-Lagrange equations have terms which are third order polynomials of the displacements.  There are no discontinuities in the problem.

My question is about the relative advantages and disadvantages of the two approaches.  I have been trying quite hard to solve the second problem, i.e. solve the large nonlinear system, and I just can't seem to make progress.  I have an analytic expression for the Jacobian, and I can ensure that I have a good initial guess in the sense that it is "close" to the solution as measured by the norm of the difference, but I have no guarantee that, for example, Newton's method will converge from this initial guess.  I've attempted some simple globalization strategies, such as a parametrized continuation, pseudo-transient continuation, and some ad hoc regularization.  The problem with Newton's method without a continuation is that it seems to always find a point at which the Jacobian becomes singular.  The parametrized continuation seems to avoid this, but the Newton iterations stagnate.  Quasi-Newton methods like L-BFGS work very well up to a point, and then fail catastrophically, with massive increase in the residual norm in a single iteration.

I have not attempted recasting the problem as a minimization for two reasons - first, I don't understand what advantages that would have over solving the nonlinear PDEs, and second, I would have to do a lot of calculations arising from discretizing at an early step.  I realize this question may be a bit vague but I'd really appreciate any guidance - thanks!

Edit: I'll try to include some more information.  One advantage I can see of discretizing the energy functional and approaching the problem as a minimization is that the Hessian would be symmetric.  This contrasts with the discretization of the Euler-Lagrange equations, for which the Jacobian is not symmetric.  In both approaches the Hessian/Jacobian is not positive definite.  It should be possible to calculate a value which can be added to the diagonal to make a positive definite matrix, as discussed in Nocedal and Wright, but I'm unsure how to actually implement that.  I suppose the idea is to put in a positive definite approximation to the Jacobian without altering the function itself.

Edit 2: Another thought about why I seem to be having so much trouble is that because my external data comes on a Cartesian grid, I am using finite difference instead of finite element.  I realize that the vast majority of work in elasticity uses finite element, but I have yet to see an explicit discussion of the advantages in terms of convergence.  I use centered difference operators on a 27 point stencil.  The boundary conditions are homogeneous Dirichlet.


问:
Is there any ML solutions for permutations?
In Machine Learning field, is there any work done, or research, or info in general about finding optimal permutations?

For eg. to find a number of optimal arrangements of sliced shapes on 2D sheet of material - so that there is low number of leftovers?

I assume it can be done by brute-forcing permutations, but what about neural networks or Machine Learning in general?


问:
Low memory algorithm for matrix diagonalisation
I'm trying to find the largest eigenvalues of very large $N \times N$ matrices  ($N = 10^{10}$ and larger). The matrices are not sparse but the multiplication operation is fast.

For now, I'm using methods without explicit matrix storage, from the simple power method to Krylov-Schur (using petsc/slepc in C). But all these methods require to store at least a vector of size $N$.

Is there a diagonalisation algorithm which does not need as much memory?


问:
Crank-Nicolson integrator: oscillations with complex matrix
I'm working on a Real-Time TDDFT implementation and I am currently comparing different propagation schemes for the propagation of the Kohn-Sham wave function,
$$
\phi(t+\Delta t) = \hat{\mathcal{U}}\phi(t).
$$
An external field is incorporated by a vector potential and the kinetic energy matrix element is modified such that the hamilton matrix gets a complex part:
$$
\mathcal{H}_{ij} = \mathcal{H}_{ij}^{(R)} + i \mathcal{Z}_{ij}
$$
where $\mathcal{H}^{(R)}$ denotes all the real parts in the hamilton matrix and $\mathcal{Z}$ is a matrix that incorporates the imaginary part of the terms arising due to the vector potential. I think details are not that important here.
Among other things, I use the instantaneous total energy $E_\text{tot}(t)$ to compare different propagators$^*$.

Most implemented integration schemes are exponential-based, i.e.
$$
\mathcal{U}=\exp\left(f(\mathcal{H})\right),
$$
e.g. Exponential Midpoint with $f(\mathcal{H})=-i\mathcal{S}^{-1}\mathcal{H}(t+\Delta t/2)$ or Enforced Time-Reversal Symmetry with $f(\mathcal{H})=-i\Delta t/2(\mathcal{S}^{-1}\mathcal{H}(t+\Delta t)+\mathcal{S}^{-1}\mathcal{H}(t))$ ($\mathcal{S}$ is the overlap matrix). These methods mostly compute the same results but when I use the Crank-Nicolson propagator,
$$
\mathcal{U}=\frac{\mathcal{S}-i\mathcal{H}(t)\Delta t/2}{\mathcal{S}+i\mathcal{H}(t)\Delta t/2},
$$
it only reproduces the results when $\mathcal{Z}=0$, else I get (a while after the pulse ends, constant) amplitude oscillations that even persist after a pulse is switched off. Reducing the step size does not change anything.

So my question is what the reason for this could be. I don't know too much about details of applying these propagation schemes to the given kind of problem.

$^*$ $E_\text{tot}(t) = \sum_i f_i \epsilon_i(t) - F[\rho(t)] + E_{nn}$ where $\epsilon_i$ are the eigenvalues of the time-dependent hamiltonian, $F$ denotes some functionals incorporating XC and electrostatic corrections and $E_{nn}$ is the internuclear repulsion energy.


问:
FiPy with derivative source terms
I have a coupled nonlinear PDE system in 1 spatial dimension, in which I want to solve using FiPy. The dependent variables are $n$ and $T$:
\begin{align}
\frac{\partial n}{\partial t} \,&=\, D\,\frac{\partial^2 n}{\partial x^2} \\
\frac{\partial T}{\partial t} \,&=\, \frac{D}{\zeta}\,\frac{\partial^2 T}{\partial x^2} \,+\, \left(\frac{\zeta + 1}{\zeta}\right) \frac{D}{n} \, \frac{\partial n}{\partial x} \, \frac{\partial T}{\partial x}
\end{align}

As of now, both $D$ and $\zeta$ are simple numerical parameters.

I am having serious trouble representing the second term in the $\frac{\partial T}{\partial t}$ equation, as there seems to be no nice way of writing it. Should it be written as a Convection term, or as a Source term?

This is the code so far, omitting some things:

```<code>L = 5.0
nx = 100
mesh = Grid1D(nx=nx, dx=L/nx)
x = mesh.cellCenters[0]

density = CellVariable(name=r"$n$", mesh=mesh)     # The n dependent variable
temperature = CellVariable(name=r"$T$", mesh=mesh) # The T dependent variable

D = 5.0
zeta = 0.5

density_equation = TransientTerm(var=density) == \
                    DiffusionTerm(coeff=D, var=density)

temp_equation = TransientTerm(var=temperature) == \
                 DiffusionTerm(coeff= D / zeta, var=temperature) \
                          + S_T    # <- What should be filled in for this?

semifull_eq = density_equation & temp_equation

if __name__ == '__main__':
    viewer = Viewer((density, temperature))

timeStepDuration = dx**2 / (2*diffusivity)
for steps in range(100):
    semifull_eq.solve(dt=timeStepDuration)
    if __name__ == '__main__':
        viewer.plot()
</code>```

I have left out the boundary and initial conditions for brevity.

Title: FiPy with derivative source terms


问:
Newton - Raphson method : maxima of function in 2 variables
I am computing the maximum of a function (with two-variables) using Newton-Raphson method. The function is : $e^{-(x \ - x_0)^2 - (y \ - y_0)^2}$, whose maxima exists at $(x_0,y_0)$.

The Jacobian will be a function of the second derivatives (including mixed) of this function. It can be written as : $J = e^{-(x \ - x_o)^2 - (y \ - y_o)^2} \
\begin{Bmatrix}
-2 \ + \ 4 \times (x \ - \ x_o)^2 && 4 \times (x \ - \ x_o) \ (y \ - \ y_o) \\
4 \times (x \ - \ x_o) \ (y \ - \ y_o) && -2 \ + \ 4 \times (y \ - \ y_o)^2
\end{Bmatrix}$

The above expression analytically calculates the second derivatives to evaluate the Jacobian. The Newton-Raphson can now be formulated as :

$$
\vec{x}^{\, n+1} \ = \ \vec{x}^{\, n} \ - \ J^{-1} \ f(\vec{x}^{\,n})
$$

where $\vec{x}$ includes $x$ and $y$.

An alternative way to evaluate the Jacobian (second derivatives) includes using Finite Differences like : $J_{i,j} (\vec{x}) \ = \ \frac{f_i(\vec{x} \ + \ h_j \vec{e_j}) \ - \ f_i(\vec{x})}{h_j}$, where  $h_j$ is approximately taken to be $x_j \ \times \sqrt{\epsilon}$, with $\epsilon$ being the computer precision (generally taken as $10^{-16}$). This can be roughly called as Secant Method (though not exactly, as it means).

Practically, analytical way of calculating the Jacobian should yield faster convergence (in a few iterations) in comparison to the aternative way described above. However, I seem to get the converse result (for the guess values ($ x_0$ and $y_0$ respectively) taken to be 3.6 and 4.6 in the code below !). In other words, I seem to get the converged solution (for the case of Analytical Jacobi) after 26 iterations, whereas 15 for the other case.

Here is the code, where I have implemented this.

```<code>#include<iostream>
#include<cmath>
#include<fstream>
#include<iomanip>

using namespace std;

double eps = 1e-16;

void user_input(double x0, double y0)
{
        cout << "Enter the value of x0 - to find the maximum of the exponential function\n" << endl;
        cin >> x0;
        cout << "Enter the value of y0 - to find the maximum of the exponential function\n" << endl;
        cin >> y0;

}

//evaluates the exponential at x0 and y0
double exp_fn(double x, double y, double x0, double y0)
{
    double val;
    val = exp(-(x - x0)*(x - x0) - (y - y0)*(y - y0));
    //cout << val << endl;

    return val;
}

//x- gradient of the given function
double f_x(double x, double y, double x0, double y0)
{

  double value;
  value = -2*exp_fn(x,y,x0,y0)*(x - x0);

  return (value);
}

//y- gradient of the given function
double f_y(double x, double y, double x0, double y0)
{

  double value;
  value = -2*exp_fn(x,y,x0,y0)*(y - y0);

  return (value);
}

//first component in Jacobian
double f_xx(double x, double y, double x0, double y0)
{
        double value;
        value = exp_fn(x,y,x0,y0)*( -2 + 4*(x - x0)*(x - x0));

        return value;
}

//second component in Jacobian
double f_xy(double x, double y, double x0, double y0) //returns the same value as f_yx
{
        double value;
        value = exp_fn(x,y,x0,y0)*(x - x0)*(y - y0)*4;

        return value;
}

//last component in Jacobian
double f_yy(double x, double y, double x0, double y0)
{
        double value;
        value = exp_fn(x,y,x0,y0)*( -2 + 4*(y - y0)*(y - y0));

        return value;
}

//equivalent to f_xx, but using finite difference, not analytical
double f_xx_sec(double x, double y, double x0, double y0)
{
  double value, h_x;
  h_x = x*sqrt(eps);

  value = (f_x(x + h_x, y, x0, y0) - f_x(x, y, x0, y0))/h_x;

  return value;
}

//equivalent to f_xy, but using finite difference, not analytical
double f_xy_sec(double x, double y, double x0, double y0)
{
  double value, h_y;
  h_y = y*sqrt(eps);

  value = (f_x(x, y + h_y, x0, y0) - f_x(x, y, x0, y0))/h_y;

  return value;
}

//equivalent to f_yx, but using finite difference, not analytical
double f_yx_sec(double x, double y, double x0, double y0)
{
  double value, h_x;
  h_x = x*sqrt(eps);

  value = (f_y(x + h_x, y, x0, y0) - f_y(x, y, x0, y0))/h_x;

  return value;
}

//equivalent to f_yy, but using finite difference, not analytical
double f_yy_sec(double x, double y, double x0, double y0)
{
  double value, h_y;
  h_y = y*sqrt(eps);

  value = (f_y(x, y + h_y, x0, y0) - f_y(x, y, x0, y0))/h_y;

  return value;
}

//determinant of the matrix
double det(double a, double b, double c, double d)
{
    double det_val;
    det_val = (a*d - b*c);
    return det_val;

}

int main(int argv, char** argc)
{

    int user = 0, iter = 1; //user = 0 determines the default x0 and y0 values to be used.
    //Else, for user to input these values, equate user = 1

    int secant_met = 1; //1 causes secant method of Jacobian computation. 0 uses the conventional way

    double x0 = 4, y0 = 7, x_old =3.6, y_old =6.7, error = 1., magnitude; //x0, y0 - guess  values

    double x_new, y_new, mult_x, mult_y, lambda = 1.0; // lambda refers to the relaxation factor

    ofstream fout;
    fout.open("convergence_data.dat");
    fout << "Error" << "\t" << "X_updated" << "\t" << "Y_updated" << "\t" << "Iteration" << endl;

    if (user == 1)
      user_input(x0,y0);

    while(error > 1e-8)
    {

      if (secant_met ==1) //using the secant way of evaluating Jacobian - Finite differences
        {
          mult_x = f_x(x_old, y_old, x0, y0)*f_yy_sec(x_old, y_old, x0, y0) - f_y(x_old, y_old, x0, y0)*f_xy_sec(x_old, y_old, x0, y0);
          mult_y = f_y(x_old, y_old, x0, y0)*f_xx_sec(x_old, y_old, x0, y0) - f_x(x_old, y_old, x0, y0)*f_yx_sec(x_old, y_old, x0, y0);

          x_new = x_old - lambda*mult_x/det(f_xx_sec(x_old, y_old, x0, y0), f_xy_sec(x_old, y_old, x0, y0), f_yx_sec(x_old, y_old, x0, y0), f_yy_sec(x_old, y_old, x0, y0)); //Generalized newton's method
          y_new = y_old - lambda*mult_y/det(f_xx_sec(x_old, y_old, x0, y0), f_xy_sec(x_old, y_old, x0, y0), f_yx_sec(x_old, y_old, x0, y0), f_yy_sec(x_old, y_old, x0, y0));
        }

      else // implementing conventional way of computing Jacobian
        {
          mult_x = f_x(x_old, y_old, x0, y0)*f_yy(x_old, y_old, x0, y0) - f_y(x_old, y_old, x0, y0)*f_xy(x_old, y_old, x0, y0);
          mult_y = f_y(x_old, y_old, x0, y0)*f_xx(x_old, y_old, x0, y0) - f_x(x_old, y_old, x0, y0)*f_xy(x_old, y_old, x0, y0);

          cout << setprecision(7) << f_xx(x_old, y_old, x0, y0) << "\t" << setprecision(7) << f_xy(x_old, y_old, x0, y0) << "\t" << setprecision(7)<< f_yy(x_old, y_old, x0, y0) << endl;

            x_new = x_old - lambda*mult_x/det(f_xx(x_old, y_old, x0, y0), f_xy(x_old, y_old, x0, y0), f_xy(x_old, y_old, x0, y0), f_yy(x_old, y_old, x0, y0)); //Generalized newton's method
            y_new = y_old - lambda*mult_y/det(f_xx(x_old, y_old, x0, y0), f_xy(x_old, y_old, x0, y0), f_xy(x_old, y_old, x0, y0), f_yy(x_old, y_old, x0, y0));
        }

        magnitude = fabs( (x_new - x_old)*(x_new - x_old) + (y_new - y_old)*(y_new - y_old));
        error = sqrt(magnitude); //calculating the error for convergence

        x_old = x_new;
        y_old = y_new;

        cout << fixed;
        cout << "Error = " <<  error << "\t" << "Iteration" << "\t" <<  iter << endl << "X updated" << "\t" << x_new << "\t" << "Y_updated" << "\t" << y_new << endl;
        fout << fixed;
        fout << error << "\t" << x_new << "\t" << y_new << "\t" << iter << endl;

        iter++;
    }

    fout.close();
    return 0;
}
</code>```

Can someone give me an insight as to where I am going wrong in the code, wherein I have implemented both these methods ? Any kind of help would be deeply appreciated.

EDIT:

This is the plot of convergence error vs iterations which I seem to get when using both the methods, as described above.


问:
Second-order PDE with seven variables
I need to solve the following partial differential equation in seven variables with four boundary conditions. I don't think Mathematica has the capacity to solve this differential equation.

Do you know of scientific packages which can be used to solve this equation and plot graphs of the dependent variable against each of the independent variables?

Title: Second-order PDE with seven variables


问:
Computational approaches to analyse numerical solutions
I apologise beforehand if the question isn't well defined or is too broad, I'm just having some difficulty finding the information I need. I also apologise if this should have been asked on stack overflow.

I have computed some numerical solutions using spectral methods to some different hydrodynamics problems that I'm interested in. The solutions are given as time series arrays. For example, at each time $t = 1, \dots, n,$ the solution for the velocity field $v$ is an $m \times m$ array.
I would like to analyse my system to determine as much information about it as possible. For example, in my 2D vortex simulations, I would like to numerically determine the dominant behaviour of the system, or numerically find/predict any splitting or merging of vortices, or perhaps find relationships between the different variables of the system etc. The problem is, I don't actually know of many computational methods that I can use to determine these kinds of behaviour. Some things I do know we can use (and have implemented) are

Autocorrelation functions (in both time and space).

Singular value decomposition (to determine dominant behaviour).

Analysing the spectra of the solution (radially averaged bins, FFT).

Basic statistical analysis, analysing fluctuation data.

Coupling tracer particles/dyes into the DE system to observe phenomena.

Visualisation of the numerical solutions.

Obviously, comparing the numerical solution with the analytical solution is also a useful technique, though most of the time this isn't possible in hydrodynamics.

I was hoping that there were more numerical techniques available that I could learn and use to study my system. So I was wondering if someone could let me know of any papers or books that detail approaches to analysing numerical solutions, or could list down some methods that they know of. Any help would be greatly appreciated.

EDIT 02/12

I was given a review paper today by Jiang et al. which outlines nine numerical methods for detecing vortices, and also a paper by Jeong et al.. It is mostly these types of algorithms/computational methods that I am interested in using.


问:
Vectorised root finding in Python
I have an array of size (254, 80) which I am trying to use Scipy's fsolve on. I have found that the speed of using fsolve on a vector is quicker than it is in a for loop but only for vectors upto about 100 values long. After this, the speed quickly drops off and becomes very slow, sometimes completely stopping.

I'm currently looping through one dimension of the array and using a vectorised fsolve on the smaller dimension but it's still taking longer than I would expect/like.

Does anyone have a good work around for this or a know of a similar function which will be happy handling a vector of a larger size? I have tried Scipy's Newton Krylov function which, upon testing, is really fast for simple functions but it's slower than fsolve on my function. - Will this be faster if I use the analytical jacobian? If so, how do I implement it?

I wondered if it is a problem with my initial solutions but the difference is in the range of ~0.001 - 0.8 which doesn't seem too bad. Although perhaps it is this?

Here's the current code:

```<code>for i in range(array.shape[0]):
    f = lambda y: a[i] - m[i]*y - md[i]*(( y**4 + 2*(y**2)*np.cos(Thetas[i,:]) )**0.25)
    ystar[i,:] = fsolve(f, y0[i])
</code>```

(Variables a, m and md are all the same size i.e. (254,80))

Digging in to this further, I have found that a function such as

```<code>f = lambda y: y*np.tanh(y) - a0/(m**2)
</code>```

is faster to solve than

```<code>f = lambda y: (m**2)y*np.tanh(y) - a0
</code>```

where m and a0 are large 2D np arrays.

Can anyone explain why this is?

Title: Vectorised root finding in Python


问:
What are the numerical properties to consider between Augmented Lagrangian and the Penalty Method?
I'm interested in (locally) minimizing a smooth nonconvex objective function:
$$
f(\textbf{x}_1, \textbf{y}_1,\cdots, \textbf{x}_n, \textbf{y}_n)=\sum_{i=1}^ng(\textbf{x}_i, \textbf{y}_i)
$$
Subject to $\textbf{y}_{i+1}=\textbf{h}(\textbf{x}_i, \textbf{y}_i)$ for each $i$.

What are the pros and cons of using the penalty method vs Augemented Lagrangian for this problem?

Penalty method. Until convergence:

Minimize wrt $\textbf{x}_i, \textbf{y}_i$ the objective $f+\sum_i\nu_i\|\textbf{y}_{i+1}-\textbf{h}(\textbf{x}_i, \textbf{y}_i)\|$ with gradient descent.

Update  $\nu_i\gets\nu_i+\alpha\|\textbf{y}_{i+1}-\textbf{h}(\textbf{x}_i, \textbf{y}_i)\|^2$

Augmented Lagrangian method. Until convergence:

Minimize wrt $\textbf{x}_i, \textbf{y}_i$ the objective $f+\sum_i\boldsymbol\lambda_i\cdot(\textbf{y}_{i+1}-\textbf{h}(\textbf{x}_i, \textbf{y}_i))+\sum_i\nu_i\|\textbf{y}_{i+1}-\textbf{h}(\textbf{x}_i, \textbf{y}_i)\|^2$ with gradient descent.

Update $\boldsymbol \lambda_i\gets \boldsymbol\lambda_i+\alpha(\textbf{y}_{i+1}-\textbf{h}(\textbf{x}_i, \textbf{y}_i))$ and $\nu_i\gets\nu_i+\alpha\|\textbf{y}_{i+1}-\textbf{h}(\textbf{x}_i, \textbf{y}_i)\|^2$


问:
Reference request: free/open source books on intro numerical PDE
I'm interested in finding a free online textbook that is an introductory textbook on numerical PDE for engineers. It should cover the basics of finite differences, finite elements, and Fourier series methods for very simple PDE cases, mainly the 1D heat and 1D wave equations, and should not be concerned with rigor, proofs of convergence, etc.

The purpose of this request is to find a decent online resource for students to use as an alternative to an official course textbook. The students are undergraduate engineers who have seen multivariate calculus and linear algebra but no other mathematics.

Does anybody have suggestions?


问:
Odd-even decoupling at faces of cells
I am currently solving PDEs using the finite volume method. The surface integrals of the equations that I am solving involves computing face gradients. The current algorithm that we use to compute face gradients is
$$\nabla\phi_{face}=\frac{1}{2}(\nabla\phi_L + \nabla\phi_R) - \boldsymbol{\hat{d}}(\frac{1}{2}(\nabla\phi_L + \nabla\phi_R)\cdot\boldsymbol{\hat{d}})+(\phi_L-\phi_R)\frac{\boldsymbol{\hat{d}}}{|\boldsymbol{d}|}$$
The LHS is the face gradient of some quantity $\phi$. $\nabla\phi_L$ and $\nabla\phi_R$ are the cell gradients at the left and right neighboring cells of the face, respectively. $\boldsymbol{\hat{d}}$ and $\boldsymbol{d}$ are the unit distance and distance vector between the left and right cell centers, respectively. The RHS consists of 3 terms. The latter 2 terms are the "odd-even decoupling" correction terms that are applied to avoid numerical oscillations that would otherwise arise if the computation of the face gradient was a simple averaging (i.e., just the first term on the RHS). </p>

<p>I have some questions regarding odd-even decoupling:
(1) Is odd-even decoupling also necessary for computing non-gradient values at the face? For example, if instead of computing the temperature gradient at the face, I instead wanted the temperature at that face. (2)
Is there a variant of this equation for non-gradient terms? For example, what if I just want to compute a face value? This equation can almost be duplicated for a plain $\phi$, non-gradient, term if it weren't for the 3rd term on the RHS. </p>

Title: Odd-even decoupling at faces of cells


问:
Does radioisotope decay affect computer systems?
It is known that there exist a number of radioisotopes of elements commonly used in computer systems. Is the decay of these materials known to affect device performance over time, or is its impact so small as to be immeasurable?

This question occurred to me last night as a friend "cleaned the N64 cartridge" by blowing on the contacts. I've also heard that cosmic ray collisions regularly affect satellite computation. There are a number of "material parameters" to computational science we, as developers and computer scientists, rarely consider.


问:
Rigid and flexible ZIF-8 framework
What are the differences in parameters in the force field for MD simulation study of RIGID and FLEXIBLE ZIF-8 framework?Is the electrostatic interaction term present in the rigid one?

Title: Rigid and flexible ZIF-8 framework


问:
BFGS convergence problem
I would like first to state that it is beyond my capability to identify whether this is a BFGS issue or a R package problem.

I've been doing some mixed logit regression using the R package "mlogit". In total I have 32 random parameters and one fixed parameter to be estimated. Correlations are allowed among random parameters.
With BFGS method, the code stopped after 13 iteration. Judging from the Chi-square value, it is far from converge. With BHHH method, however, it usually converge (although very slowly). Is this a package problem? or is it BHHH method more appropriate than BFGS in my context?

I doubt if it is related to the step size, but it's only my doubt.
The code and data I've used are here.

Title: BFGS convergence problem


问:
Band structure of nonlinear Schrodinger equation with one dimensional potential
I have a nonlinear Schrodinger equation which reads:

$$ \frac{1}{2} \frac{d^2u}{dx^2}+ |u|^2u + V(x)u = -i \frac{du}{dz},$$

where $V(x)=\cos(wx)+ i a \sin(wx)$ and $w$, $a$ are numbers.

How to solve its band structure, meaning its dispersion relation?

The plane wave expansion method is used for a real potential. Can it be applied for this complex potential also?

What are other methods, if any, in Matlab or Mathematica or other software?


问:
Nonlinear global optimization algorithm that can use dynamic programming
I've asked this question on stackoverflow 2 weeks ago, but, judging by zero response, that probably was the wrong forum. Therefore copying it here:

Let F0,...,Fn be functions with parameters p0,...,pn and G be a composite function. Say, for n=3 G(x)=F3(p3,F2(p2,F1(p1,F0(p0,x)))).

I'd like to globally optimize the parameters p0,...,pn with making as few function calls as possible. There are several ways to optimize:

1) Ignore that G is a composition of Fis, treat it as a blackbox with parameters p0,...,pn, and use any generic global optimizer.

2) Optimize hierarchically; for clarity let's assume n=3:

Loop A: Freeze parameters p0,p1,p2 and optimize p3 in isolation. because G depends only on p3 and the output of prior functions that have frozen params I only have to evaluate F3; the outputs of other Fs don't change between iterations.

Loop B: Freeze parameters p0,p1, change p2, and repeat the above loop A. Repeat for optimizing p2.

Loop C: Freeze p0, change p1, repeat loop B, thus optimizing p1.

Loop D: Obvious.

This "dynamic programming" approach would lead to optimizing G by, hopefully, making fewer iterations. I have a couple of questions though:

First, I'm not sure that such hierarchical optimization approach would lead to solution of the same quality as the generic "blackbox" approach. Has any research have been done regarding the quality of results produced by hierarchical versus blackbox optimization?

Second, are there any optimization packages available out there supporting hierarchical optimization like that? I'm aware of a few available generic blackbox optimizers (NLOpt, etc), but not generic hierarchical optimizers.


问:
What is Chebfun `eigs` doing
What is this doing? Looks like the original eigenvalue problem is converted into generalized eigenvalue problems with different dimensions of collocation points. Can someone explain more about this?

eigs.m

```<code>% Linear combination coefficients for convergence test. The convergence of the
% combination is the same as the worst constituent function. The nontrivial
% coefficents are to make accidental cancellations extremely unlikely.
coeff = 1./(2*(1:k)');
for dim = [dimVals NaN]
    [V, D, P] = getEigenvalues(discA, discB, k, sigma);

    % Combine the eigenfunctions into a composite.
    v = V*coeff(1:size(V,2));

    % Convert the different components into cells
    u = partition(discA, P*v);

    % Test the happiness of the function pieces:
    vscale = zeros(1, sum(isFun));   % intrinsic scaling only
    [isDone, cutoff] = testConvergence(discA, u(isFun), vscale, prefs);

    if ( all(isDone) )
        break
    elseif ( ~isnan(dim) )
        % Update the discretiztion dimension on unhappy pieces:
        discA.dimension(~isDone) = dim;
    end
end
</code>```

eigs.m/getEigenvalues

```<code>function [V, D, P] = getEigenvalues(discA, discB, k, sigma)
% Formulate the discrete problem and solve for the eigenvalues

    % Discretize the LHS operator (incl. constraints/continuity):
    [PA, P, C, ignored, PS] = matrix(discA);

    % Discretize the RHS operator, or use identity.
    if ( ~isempty(discB) )
        % TODO: This is untidy. Can we make a method to do this? NH Apr 2014.
        discB.dimension = discA.dimension;
        PB = matrix(discB);
        % Project RHS matrix and prepend rows for the LHS constraints.
        PB = [ zeros(size(C)) ; PB ];
    else
        PB = [ zeros(size(C)) ; PS ];
    end

    % Compute eigenvalues.
    if ( length(PA) <= 2000 )
        [V, D] = eig(full(PA), full(PB));
        lam = diag(D);

        % Find the ones we're looking for.
        lam = deflate(lam, size(C,1));
        idx = nearest(lam, sigma);
        idx = filter(idx, P*V, k, discA);

        % Extract them:
        V = V(:,idx);
        D = D(idx,idx);

    else
        % TODO: Experimental.
        [V, D] = eigs(PA, PB, k, sigma);
    end

end
</code>```

For the full code see https://github.com/chebfun/chebfun/blob/development/%40linop/eigs.m

Title: What is Chebfun `eigs` doing


问:
Understanding Boundary Condition in FEM
I am trying to understand Dirichlet and Neumann boundary conditions in FEM and I wanted to know if my inference is correct. To articulate my understanding, lets consider a simple case of TE and TM propagation in rectangular waveguides.

In the case of TE ($E_z$=0 and $H_z$=$\phi$) and in the case of TM ($H_z$=0 and $E_z$=$\phi$). Consequently, in TE we solve for $H_z$ using Helmholtz equation with Neumann boundary condition on all sides and in the case of TM we solve for $E_z$ with Dirichlet boundary condition on all sides.

The particular method of FEM and assembly procedure I am using is given by Clayton R Paul and Sadiku wherein the nodal value vector is bifurcated into two sub vectors namely the free node vector($ \phi_{f} $) and the prescribed node vector($ \phi_{p}$). The prescribed node vector is filled by applying Dirichlet boundary conditions to exterior nodes or by assigning excitation to nodes. The free node vector is unknown and has to be computed.

My question is particularly about Neumann boundary condition. My understanding is that since Neumann is implicit and a natural boundary condition in FEM, the node vector($ \phi$) is entirely comprised of free node vector implying that every node in the domain is free and unknown. Is my understanding correct?

Title: Understanding Boundary Condition in FEM


问:
Vortex Lattice Method: better basis than horse-shoe vortex?
I read several introductory texts about potential flow and vortex lattice method. Basically, it is fitting of some velocity field described by conditions on velocity at some control points using basis functions placed in other points.

horse-shoe vortex:

In order to fullfill Helmholdz theorem that vortex filaments $\Gamma$ are closed loops (does not end in the air) horse-shoe vortex is used as basis function. But it is not very nice because:

It is quite computationally costly to evaluate resulting velocity field using Biot–Savart law

It imposes that vortex fillament is stright line behind the wing while in reality vortex fillaments curve and merge.

Delocalized trailing vortex:

Since we don't know how exactly filament looks like far behind the wing it would be better not to impose particular shape, but rather assume its position is uncertain (=delocalized) => vorticity $\Gamma$ smoothly blurred behind the wing. Do we really need strictly irotational vortex which concentrate vorticity in singular line?

Brainstorming:

So I'm thinking if it would be possible to use somethink like rotational dipole which is velocity field generated by one infinitisimal element of vortex filament $d\Gamma$. I'm not exaclty sure about formula but something like:

$\vec \phi(\vec r)=(v_x,v_y,v_z) = (0, z/|\vec r|^3, y/|\vec r|^3)$

This function is both faster to evaluate than horse-shoe and smoothly decays behind the wing. But I'm not sure how well it fulfills condition for proper potential flow basis function. Maybe with little modifications the errors (residual divergence and rotation) can be minimized ?

Maybe residual of such simple basisfunctions can be fitted to capture viscosity (like Lamb–Oseen vortex )? Possibly somebody already have done this?


问:
Structural mechanics traction boundary condition question
In structural mechanics, are the boundary conditions "free surface," "Traction free", "stress free" all equivalent Neumann boundary conditions? </p>


问:
Obtaining the lagrangian multipliers in an optimization problem
Suppose we have this simple optimization problem
\begin{align*}
    \underset{x\in V}{\text{min}} &~ f(x) \\
    \text{s.t.}& ~x \leq \beta
\end{align*}

Using slack variables
\begin{align*}
    \underset{x\in V}{\text{min}} &~ f(x) \\
    \text{s.t.}&~ x - \beta + s \leq 0 \\
               & s \geq 0
\end{align*}

Lagrangian
\begin{equation*}
    \mathcal{L} \left(x, \lambda, s, \mu \right) = f(x) + \langle \lambda, x -\beta + s \rangle_{L_2} +\langle \mu, s \rangle_{L_2}
\end{equation*}
KKT conditions $\forall \delta \lambda, \delta x, \delta \mu, \delta s$
\begin{align*}
    \mathcal{L}_x (\delta x) &= \langle f'(x), \delta x \rangle_{L_2} + \langle \lambda, \delta x \rangle_{L_2} = 0\\
    \mathcal{L}_{\lambda} (\delta \lambda) &= \langle x - \beta + s, \delta \lambda \rangle_{L_2} = 0 \\
    \mathcal{L}_{s} (\delta s) &= \langle \lambda , \delta s \rangle_{L_2} + \langle \mu, \delta s \rangle_{L_2}= 0 \\
    s \geq 0 ~ \mu \leq 0 & ~\mu s= 0 ~\text{a.e.}
\end{align*}
Now assume we use an optimization algorithm that replaces the function $f(x)$ with an approximation easier to solve. It also replaces the box constraint $\beta$ with a value closer to the current optimization iteration.  One such method is MMA. Because of this substitution, the KKT conditions are different throughout the optimization. After all, the box constraint is different and their lagrange multiplier will be different as well. Is there any way to recover the lagrange multipliers of the original problem during the optimization?


问:
How do 'virtual kinematics/functions' play a role during deriving weak form formulations for physical problems?
I wanna ask a question that confuses me quite a long time. I saw many guys, in the context of computational mechanics,  they seemed to choose the virtual functions or kinematics in a way that some kinematical quantities are kept/fixed zero, but  they ask the remaining  of  them   to act as roles during the derivation of  a weak form equation. In the end, we can generally get a system of weak equations. I am all the time wondering whether we have a chance to change the nature of the original physical problem of a strong form.  Who can shed lights on my doubts in a physcial or mathematical view?

Edit:
Say, u=u1+u2+u3 for whatever reason. Its variation is delta u1+ delta u2 +delta u3. The weak form may be integral sigma(u): sym gradient(delta u1) =0, integral sigma(u): sym gradient(delta u2) =0, integral sigma(u): sym gradient(delta u3) =0.


问:
Complex Integral Equation Solution in MATLAB
I need to solve an integral equation in the form:

$$A(z)+\int\limits^{z_2}_{z_1}B(z') \frac{z^N}{z^N-z'^N} \frac{e^{i\beta}}{|z|}\mathrm{d}z'=0 $$

where $A(z)$ distribution is known and we are solving for $B(z')$ distribution and $z,z'$ are complex variables. I also know the value of $B(z_2)$.

Can anyone help me how I can solve this equation using MATLAB (preferably numerically)?


问:
How to apply an integrated constrain condition in FEM?
I'm running some simulation using FEM. In my model I need to apply a constraint condition to the governing equation.  My governing equation similar to the diffusion equation as below:

$$\frac{\partial c}{\partial t}=\nabla(D\nabla c)\quad\text{on}\quad\Omega$$

and the flux is defined as follow:

$$-D\nabla c\cdot n=J$$

And the constraint condition is an integral. For short, I need to constrain the total flux alone the specific surface equals to a constant, which means:

$$\int_{\partial\Omega_{1}}(-D\nabla c\cdot n)d\Gamma=\text{Constant}$$
where $$\partial\Omega_{1}$$ is part of my domain's surface, not the total surface.

I can write out the weak form for the governing equation(based on N-R method) as follow:
$$R_{c}^{I}=\int_{\Omega}\dot{c}N^{I}d\Omega+\int_{\Omega}D\nabla c\nabla N^{I}d\Omega-\int_{\partial\Omega}D\nabla c\cdot n N^{I}d\Gamma$$
and the related stiffness matrix:
$$K_{c\dot{c}}^{IJ}=-\frac{\partial R_{c}^{I}}{\partial\dot{c}^{J}}=-\int_{\Omega}N^{J}N^{I}d\Omega$$
$$K_{cc}^{IJ}=-\frac{\partial R_{c}^{I}}{\partial c^{J}}=-\int_{\Omega}D\nabla N^{J}\nabla N^{I}d\Omega+\int_{\partial\Omega_{1}}D\nabla N\cdot n N^{I}d\Omega$$

How can I apply $$\int_{\partial\Omega_{1}}(-D\nabla c\cdot n)d\Gamma=\text{Constant}$$ to my FEM system?

Someone told me that Lagrange multiplier can solve this problem, but I still don't know how to do it. I'm a newbie of FEM, so any details are helpful for me.


问:
Dealing with non-physical negative ODE solutions using ODEPACK
Hi and thank you all again. I am solving the reaction-diffusion-advection equation as follows

$ \partial_{t} n\left(t,z\right) = -\partial_{z}\left(\Phi\left(t,z\right)\right)  + p\left(t,z\right) -n\left(t,z\right)l\left(t,z\right)$

The terms $p\left(t,z\right)$ and $l\left(t,z\right)$ are computed according to big sets of chemical, and photo-chemical reactions. These reactions are characterized by rates which span huge different scales. Hence the system is stiff. I clarify this point because otherwise a simple Euler-Implicit scheme is perfectly valid for any time step if only transport $\Phi\left(t,z\right)$ is present. Nevertheless the solution absolutely fails involving chemistry if I cannot handle the truncation error. Since I need to integrate for great amounts of times (up to millions of years) to achieve the stationary state, hence I need an adaptive time step scheme. I solve the stationary state by pure time evolution because of further research requirements.

I found the best solution was using ODEPACK. Taking the spatial discretization we end up with a system of ordinary differential equations in which each dependent variable is the density of a certain chemical species at a certain height.

https://computation.llnl.gov/casc/odepack/

Here I attach a publication which basically solves the same kind of system as mine (planetary atmosphere), and uses ODEPACK.

In order to use ODEPACK within my fortran90 code I work with a modern fortran ODEPACK interface to call the old fashion fortran77 ODEPACK routines. In particular the one known as DLSODAR.

I faced the problem that sometimes DLSODAR yields negative densities which does not make sense in this physical context.

My solution was simple. The integration step is made through calling another routine within DLSODAR code, named DSTODA. I introduce some code to just detect if any dependent variable becomes negative

```<code>C-----------------------------------------------------------------------
C   CALL DSTODA(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPRJA,DSOLSY)
C-----------------------------------------------------------------------
  DO 295 I = 1,N
  295    YOLD(I) = Y(I)
  TOLD = TN
  HOLD = H
  SIGNVAL = .TRUE.
  CALL DSTODA (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LEWT),
 1   RWORK(LSAVF), RWORK(LACOR), RWORK(LWM), IWORK(LIWM),
 2   F, JAC, DPRJA, DSOLSY)

  CALL POSITIVENESS (SIGNVAL,NEQ,Y) ! (30 - 11 - 2017)
  IF(SIGNVAL.EQV..FALSE.) KFLAG = - 3

  KGO = 1 - KFLAG
  GO TO (300, 530, 540, 400), KGO ! IF KFLAG = -3, solution is negative and DLSODAR returns
                                  ! from the last non negative solution to repeat if desired
                                  ! with a smaller time step H. (14 - 12 - 2017)

C-----------------------------------------------------------------------
</code>```

CALL POSITIVENESS(SIGNVAL,NEQ,Y) statement basically looks for each Y(I) possible and when it is less than zero, the routine returns the last Y with all values greater than zero, altogether with the attempted time step. Both stored before DSTODA call in YOLD HOLD respectively.

After that returned unfinished solution, I reduce the last attempted time step by a factor $1/10$ and proceed from the last step. Furthermore, when from DSTODA routine the time derivative and jacobians are called, if some density is negative, those auxiliar routines do not change the previous computed time derivative and jacobian. I found this important because sometimes DSTODA may yield a positive integrated solution but not correct. Hence not changing the time derivative in any case will change the solution's sign. I think this is due to the fact that sometimes the time derivative or jacobian are computed with intermediate Y solutions that not being positive may get a positive one.

--------------------------------------------------------------------------

It works, except for a certain kind of reactions. The physical origin of the equation does not care because as other reactions is characterized by a forward and reverse reaction rates. Anyway the forward mechanism has $k_{f} \approx {10^{-24},10^{-12}} cm^{3} s^{-1}$ and reverse rate $k_{r} \approx {10^{-38},10^{-27}} s^{-1}$. I will give the particular example as an independent simulation.

The reaction is as follows, considering we have D = 500 height layers (index j) and nElems = 3 (index i) chemical species

y(1+(j-1)*nElems) + y(2+(j-1)*nElems) <---> y(3+(j-1)*nElems)

And the code

```<code>use dlsodar_f95

implicit none
integer i,j,jj,u, D, nElems
integer lrw, liw, MU, ML
double precision p(500),kf(500),kr(500)
double precision n(3,500),sol_f(3*500),dummy1(7),dummy2(2)
double precision sol(3*500)
double precision t, xend, dt, dti, tf
double precision Rtol, Atol(3*500)
character(len=1024) forward, reverse

D = 500
nElems = 3

! Importing reaction constants
open(unit=10,file='reaction.dat')
  read(10,*) p(:)
  read(10,*) forward, kf(:)
  read(10,*) reverse, kr(:)
close(10)

! Importing intial water density profile
n(:,:) = 0.0d0
open(unit=10,file='inProf.dat')
  read(10,*)
  do i = 1,D
        read(10,*) dummy1(:), n(3,i), n(1,i), dummy2(:), n(2,i)
  end do
close(10)

! Solving using DLSODAR
t = 0.0d0
dt = 0.0d0
Rtol = 1e-1
Atol(:) = 1e-2
MU = 3
ML = 3
lrw = 22+(10 + (2*(MU)+ ML))*nElems*D
liw = 20+nElems*D
tf = 1e1
dti = 1d-15 ! dti is the trial initial time step
u = 0
do i = 1,D
  do j=1,nElems
    u = u + 1
    sol(u) = n(j,i)
  end do
end do
do while (t.lt.tf)
      xend = tf
  call ODE_DLSODAR(t,xend,sol,derivs,jac,
 *                     DLSODAR_dummy_Constraint,0,
 *                     xend,sol_f,dti,dt,nElems,nElems,liw,lrw,
 *                     relative_tolerance=RTOL,
 *                     absolute_tolerances=ATOL,
 *                     max_steps=int(1e9), ! If max_steps is too big, maybe there is an allocation error
 *                     store_steps=.false.,verbose_level=1)
  !
  ! If solution is negative this loop will repeat from last non negative
  ! iteration but with a smaller time step
  !
  t = xend
  if (t.lt.tf) then

    dti =  dt/10.
    u = 0
    do i = 1,D
      do j=1,nElems
        u = u + 1
        sol(u) = sol_f(u)
      end do
    end do
  end if
end do

u = 0
do i = 1,D
  do j=1,nElems
    u = u + 1
    n(j,i) = sol_f(u)
  end do
end do
open(unit=10,file='finProf.dat')
  do i = 1,D
        write(10,*) p(i), n(1,i), n(2,i), n(3,i),
 *                  kf(i),kr(i)
  end do
close(10)

contains
c --------------------------------------------------------------

subroutine jac(NEQ, X, Y, ML, MU, PD, NROWPD)
  implicit none
  ! Arguments
  integer, intent(in) :: NEQ,ML,MU,NROWPD
  double precision, intent(in) :: X
  double precision, dimension(NEQ), intent(in) :: Y
  double precision, dimension(NROWPD,NEQ), intent(out) :: PD
  ! Local variables
  integer i, j,u
  double precision kf(500),kr(500)
  character(len=1024) forward, reverse

  do u = 1, NEQ
    if (Y(u).NE.Y(u)) then
      print*, "yields NAN"
      stop
    end if
    if (Y(u).lt.0.0d0) then
      return
    end if
  end do

      PD(:,:) = 0.0d0

! Importing reaction constants
open(unit=10,file='reaction.dat')
  read(10,*)
  read(10,*) forward, kf(:)
  read(10,*) reverse, kr(:)
close(10)

  print*, "dentro jac:",x

  do i=1,500
      PD(2,3+(i-1)*3) =  kr(i)
      PD(3,2+(i-1)*3) = -kf(i)*Y((i-1)*3+1)
      PD(3,3+(i-1)*3) =  kr(i)
      PD(4,1+(i-1)*3) = -kf(i)*Y((i-1)*3+2)
      PD(4,2+(i-1)*3) = -kf(i)*Y((i-1)*3+1)
      PD(4,3+(i-1)*3) = -kr(i)
      PD(5,1+(i-1)*3) = -kf(i)*Y((i-1)*3+2)
      PD(5,2+(i-1)*3) =  kf(i)*Y((i-1)*3+1)
      PD(6,1+(i-1)*3) =  kf(i)*Y((i-1)*3+2)
  end do

  return
end subroutine jac
c --------------------------------------------------------------

subroutine derivs(NEQ,X,Y,DYDX)
  implicit none
  integer, intent(in) :: NEQ! number of variables
  double precision, intent(in) :: X! independent variable, time
  double precision, dimension(NEQ), intent(in) :: Y ! dependent variables. Unknowns algorithm is solving
  double precision, dimension(NEQ), intent(out) :: DYDX ! time derivative
  integer u ! running through nElems*D positions
  integer i, j, counter
  double precision dummy(2) ! dummy variable
  double precision kf(500),kr(500)
  character(len=1024) forward, reverse

  do u = 1, NEQ
    if (Y(u).NE.Y(u)) then
      print*, "yields NAN"
      stop
    end if
    if (Y(u).lt.0.0d0) then
      return
    end if
  end do

      DYDX(:) = 0.0D0

! Importing reaction constants
open(unit=10,file='reaction.dat')
  read(10,*)
  read(10,*) forward, kf(:)
  read(10,*) reverse, kr(:)
close(10)

  print*, "dentro dydx:",x

  do i=1,500
      DYDX((i-1)*3+1) = -kf(i)*Y((i-1)*3+1)*Y((i-1)*3+2)
 *                           + kr(i)*Y((i-1)*3+3)
      DYDX((i-1)*3+2) = -kf(i)*Y((i-1)*3+1)*Y((i-1)*3+2)
 *                           + kr(i)*Y((i-1)*3+3)
      DYDX((i-1)*3+3) =  kf(i)*Y((i-1)*3+1)*Y((i-1)*3+2)
 *                           - kr(i)*Y((i-1)*3+3)
  end do

  return
end subroutine derivs
c --------------------------------------------------------------

end
</code>```

And part of the output is as follows (I removed from the source code some prints* for the sake of being clearer)

```<code> dentro dydx:   4.1576048498046223
 dentro dydx:   4.1577169032470946
 dentro dydx:   4.1577169032470946
 dentro dydx:   4.1578289566895670
 dentro dydx:   4.1578289566895670
 dentro dydx:   4.2698823991657884
 dentro dydx:   4.2698823991657884
 dentro dydx:   4.1858423173086221
 dentro dydx:   4.1858423173086221
 Y(I) :  -1.4252481539273311E-013
 Where I :          25
 Ha salido negativo
 el xend alcanzado es :   4.1578289566895670
 pasos temporales, dti :   1.1205344247195033E-004 , dt :  0.11205344247622159
 dentro dydx:   4.1578289566895670
 dentro dydx:   4.1690343009371889
 dentro dydx:   4.1690343009371889
 Y(I) :  -8.0043582938772017E-007
 Where I :           4
 Ha salido negativo
 el xend alcanzado es :   4.1578289566895670
 pasos temporales, dti :   1.1205344247622160E-002 , dt :   1.1205344247622160E-002
 dentro dydx:   4.1578289566895670
 dentro dydx:   4.1589494911143294
 dentro dydx:   4.1589494911143294
 Y(I) :  -1.0391305382066468E-010
 Where I :           4
 Ha salido negativo
 el xend alcanzado es :   4.1578289566895670
 pasos temporales, dti :   1.1205344247622159E-003 , dt :   1.1205344247622159E-003
 dentro dydx:   4.1578289566895670
 dentro dydx:   4.1579410101320429
 dentro dydx:   4.1579410101320429
 dentro dydx:   4.1580530635745188
 dentro dydx:   4.1580530635745188
 dentro dydx:   4.2701065061024490
 dentro dydx:   4.2701065061024490
 dentro dydx:   4.1860664242065013
 dentro dydx:   4.1860664242065013
 Y(I) :  -1.0170965354084183E-013
 Where I :          25
 Ha salido negativo
 el xend alcanzado es :   4.1580530635745188
 pasos temporales, dti :   1.1205344247622159E-004 , dt :  0.11205344252793006
</code>```

Basically the time step is stuck between 1e-3 s and 1e-1 s. It is not worth to wait more than a day just to see if ODEPACK finally decides to increase dt. I tried changing tolerances too. But the biggest issue is that the solution is wrong too. As an example a t = 100 s I display the simulation output

In the red solution (corresponding to Y(3)) those peaks are wrong. The initial condition is taken from a well-mixed hydrostatic atmosphere with a given composition.

As an extension, if we let DLSODAR to run on its own, without any modifications at al, we have the following output, for t = 1e12 s

Where the white discontinuity is how gnuplot deals when we have negative values on a logarithmic scale. I should not accept negative solutions at any time value. This output only takes a few seconds to be computed.

The final (hopefully correct) stationary state solution if anyway we let DLSODAR to keep computing no matter negative numbers looks like

Any suggestions when ODEPACK yields negative non physical solutions?. Is this a known problem common to occur?. I feel the problem should be so easy to solve, but my mind cannot see it clearer.

Again, thank you so much.


问:
Is reduced stiffness matrix positive definite too?
The stiffness matrix $K$ in a finite element analysis is a symmetric positive definite matrix. When we introduce essential boundary conditions, we remove rows and columns associated with prescribed displacements, for example. And we obtain a linear system of equations $\bar{K}\bar{X}=\bar{F}$ in which $\bar{K}$ is the reduced stiffness matrix and is obtained by removing rows and columns of stiffness matrix $K$ associated with essential boundary conditions. My project is about damage in cyclic loading (low-cycle, $10^3$). This is a non-linear analysis because of plasticity and damage involved in. So, the run-time of the code that I am writing will be too large. I want to use Conjugate Gradient method for solving the systems of linear equations during the analysis. But the method requires that the reduced stiffness matrix $\bar{K}$ is symmetric and positive definite. The symmetry of $\bar{K}$ is trivial, but I do not know whether it is always positive definite or not.</p>


问:
Multithreading determinant computation of a Big Matrix
I've been trying to design an algorithm for computing determinant of Big square Matrices (N <= 1000). I am allowed to use multithreading, but still, I can't design the algorithm. The assignment states, that I should use Gaussian elimination method. Yet, there are no input data, so the testing data is on me. For big matrices (N ~1000) it's a binary matrix, otherwise - matrix of integers. I wonder, if there exists a relatively time efficient ( less than 10 seconds ) multithread algorithm for this task. Would be glad to hear opinions.

EDIT:: I came up with idea of using partial pivoting to find a row with maximum at current iterating position, swapping it with current row (such heuristic keeps precision the best). And then, the parallel part, is to do row adding in multithread ( I used 8 threads, since it gave the best results on quad-core) if the matrix is big enough.
Such method gave more precise results, and efficiency raised with matrix size relatively to one thread method.
It actually turned out to be quite easy to implement.


问:
Find a vector B that minimizes |W-A*B|
I want to find a candidate vector $B$ that

$$\min|(W - A_i * B_i)|$$
$$ a_i > 0,\ A_i=\{a_0,...,a_i\},\ B_i=\{-1,0,1\}^i$$

For example, given

$$W = 0.6,\quad A_4 = [0.1, 0.2, 0.4, 0.7] $$

one of answer will be

$$B_4 = [1, -1, 0, 1] $$

I only get that by searching an entire tree (DFS, BFS).
I tried with a ternary tree, but 0 value in vector $B$ shuffles the ascendant order.

Dynamic-programming doesn't seem to work, divide-conquered also look at the entire tree.

What would be an efficient approach to the problem?

ContentLicense: CC BY-SA 4.0

Title: Find a vector B that minimizes |W-A*B|


问:
When numerically computing eigenstates during a coupled-mode-space NEGF calculation, do phases matter?
The coupled mode space NEGF method for computing transistor characteristics involves expanding the electronic wavefunction in a mode space basis

$$\Psi(x,y,z) = \sum_n\phi_n(x)\xi_n(y,z;x)$$

where $\xi_n(y,z;x)$ is the nth mode at position x (the x dimension being the direction of transport).

The method involves computing coupling between modes along the transport direction. E.g. coupling coefficients like

$$b_{mn}(x) = \langle\xi_m|\frac{\partial}{\partial x}|\xi_n\rangle$$

The numerical algorithm will divide the transistor into $X$ slices $\Delta x$ apart, and compute the modes on each slice. The issue is my numerical solver will sometimes compute $\xi_n$ as the nth eigenvalue, and sometimes compute $-\xi_n$, so phases might not be consistent across each slice. This changing phase will result in a large $\frac{\partial}{\partial x}|\xi_n\rangle$

My question is: Do the phases need to be consistent across each slice? If they do, how do NEGF algorithms ensure they are consistent?


问:
Solving 2D Schrodinger Equation with ARPACK: Can I ensure all eigenvectors have the same phase?
I use arpack to solve the 2D Schrodinger, and eigenvalue problem of the form

$$Hx = \epsilon x$$

on a uniform grid. All eigenvectors are real in my case.

Arpack doesn't normalise the eigenvectors, and phases between the eigenvectors will differ (some are multiplied by -1, some aren't).

Normalising the eigenvectors after the arpack run is easy, but I am not sure how to ensure all eigenvectors have the same phase. Is there some way to compare the phase between two orthogonal eigenvectors (say, by checking the sign near a boundary)?


问:
functional second derivative
I'm trying to build a numerical solution for a parameter estimation problem of reaction-diffusion equation, using the adjoint method.

To summarize it, I'm trying to minimize the function
$$ G=\frac{1}{2}\int_\Omega{(u(T,x)-u_f)^2d\Omega} + \frac{1}{2}\int_\Omega{D^2d\Omega}+ \frac{1}{2}\int_\Omega{k^2d\Omega} $$

subject to

\begin{align}
F = \partial_t u - \nabla(D\nabla u) - ku(1-u) = 0 \\
u(0,x) = u_0 \\
u(t,x) = 0 ~\forall x\in\Gamma :=\partial\Omega
\end{align}

where $u$ is the tumor density function, $\Omega$ is the space domain of the problem and $D$, $k$ are the unknown diffusion and proliferation terms (to optimize in order to find the minimal $G$).

Thank to GoHookie's help, I have now a proper derivation of the Lagrangian
$$ L = G - \lambda F $$ giving me an expression of the total derivative of $G$ w.r.t the unknown parameters $D$ and $k$

$$
\frac{\partial L}{\partial D}(\delta D) = \int_\Omega\int_t^T{(\frac{D}{T}-\nabla \lambda\nabla u)\delta Dd\Omega dt} \\
\frac{\partial L}{\partial k}(\delta k) = \int_\Omega\int_t^T{(\frac{k}{T}-\lambda u(1-u))\delta k d\Omega dt}
$$

Since the problem is non-linear, I have to use an iterative optimization methods, which requires me also to provide an expression of the Hessian of $G$.

Using again an adjoint method, leads to calculate the various second derivatives of $L$ w.r.t $u, D$ and $k$.

I think I can manage it except for the case of
$$ \frac{\partial^2 L}{\partial u\partial D}(\delta u,\delta D) = -\int_\Omega\int_t^T{(\nabla \lambda\nabla \delta u)\delta Dd\Omega dt} $$

since then, I think I need to factorize this expression both for $\delta u$ and $\delta D$ (i.e. transform the expression in something like $\int_\Omega\int_0^T{\left( J(u,\lambda,\nabla u,\nabla \lambda)\delta u\delta \lambda\right) d\Omega dt}$ ) but none of the divergence theorem, Green theorem or integration by part methods seems to help me doing this factorization.

Any help ?

Thx

Title: functional second derivative


问:
Obtainting KKT for QSDP for the trace inequality constraint
I am working on developing my own solver(for implementation on hardware), based on IPM for following problem:
\begin{equation}
\begin{split}
\min_{X} \;  \frac{1}{2}&\|X\|_F^2 + trace(CX)\\
\text{s.t. } & trace(AX)\leq b\\
& X\succeq 0\\
\end{split}
\end{equation}
As a reference I already know how to obtain KKT equations for
\begin{equation}
\begin{split}
\min_{X} \;  \frac{1}{2}&\|X\|_F^2 + trace(CX)\\
\text{s.t. } & trace(AX) = b\\
& X\succeq 0\\
\end{split}
\end{equation}
Dual of this problem will be
\begin{equation}
\begin{split}
\max_{X,y,S}\;& -\frac{1}{2}\|X\|_F^2+by+\beta\log\det(S)\\
\text{s.t. }\;& A^Ty-X+S=C,\; S\succeq 0
\end{split}
\end{equation}
KKT system will be
\begin{equation}
\begin{array}{ccccccc}
-X & + & A^Ty & + & S & = C, & S\succeq 0\\
trace(AX) &  &  &  &  & = b, & X\succeq 0\\
&  & XS  &  &  & = 0, &\\
\end{array}
\end{equation}
Later I sove it with Newton-Rhapson method, have no problems there.
My question, therefore, be, how to obtain KKT system, when $trace(AX) \leq b$?

All responses or references are appreciated.


问:
Analysis and numerical methods for PDEs arising in industrial problems
Assume some basic knowledge of numerical methods for PDEs, acquired through A. Quarteroni's Numerical models for differential problems.

I'm looking for a reference to get started on the analysis of PDEs and modeling for industrial problems and related computational approaches (with sample Matlab code available to illustrate them).

Are there monographs that deal with concrete complicated (possibly messy)  industrial problems (elasticity, fluid mechanics, continuum mechanics, chemical reactions, etc.) involving PDEs (i.e. PDE models that arise in real-life "industrial" settings)?

I'm looking for rigorous analysis and numerical analysis. Also, I'd prefer if such books included fairly complete source codes of implementations of the proposed numerical algorithms in MATLAB.


问:
scaling in discretized PDE system
I want to solve the following system via Matlab
$\Omega=(0,1)^2$

$$\Delta y=\frac{1}{\alpha} p$$
$$ -\Delta p= y -1 $$
$$p|_{\partial \Omega}=0,~y|_{\partial \Omega}=0$$

using

```<code>    n=50;
    alpha=1;
    A=gallery('poisson',n);

    K=[A -1/alpha*speye(n^2);
       speye(n^2) A];
    b=[zeros(n^2,1) -ones(n^2,1)];
    solution=K\b;
    %embedding with boundaries...
</code>```

To my surprise this leads to (I think) a correct answer, however I usually have this scaling factor $\frac{1}{h^2}=(n+1)^2$ in it. Taking $-\Delta y =1$ with zero boundaries as an example, this leads to $$\frac{1}{h^2}Ay=1_{n\times1}$$

Is my code incorrect, or does this factor vanish somehow?

Title: scaling in discretized PDE system


问:
Why do higher order finite elements (Q2) do not perform well for large Peclet number flows, as compared to Q1 finite elments?
I a solving the 2d steady state convection-diffusion problem on the famous flow around a cylinder in rectangular domain benchmark. My numerical results show that with Q1 finite elements, the solver is converging for Pe=5000, whereas with Q2  discretization, the solver fails to converge. What can be the reason for poor performance in case of Q2 FEM discretization?  </p>


问:
(Approximate) Incremental Projection Method for Navier-Stokes equations
I am trying to implement an incremental projection method for the 2D incompressible Navier-Stokes.  The type of projection method I am trying is

$$
\frac{u^{*} - u^{n}}{dt} = - \nabla p^{n} - u \cdot \nabla u + \nu \nabla^{2} u\\
u^{*}|_{\partial D} = u^{n+1}|_{\partial D}\\
\nabla^{2} \phi^{n+1} = \frac{\rho}{dt} \nabla \cdot u^{*}\\
\nabla \phi^{n+1}|_{\partial D} = \frac{1}{dt} \left(u^{*} - u^{n+1}\right)\\
\frac{u^{n+1} - u^{*}}{dt} = - \nabla \phi^{n+1}\\
p^{n+1} = p^{n} + \phi^{n+1} (5)
$$

In terms of the velocity-pressure arrangement I am using a collocated formulation using the Rhie-Chow interpolation to have stable coupling.  I however have run into a problem in terms of the practical implementation of the pressure update.  The way I update pressure is as follows

```<code>Pnew(1:nx,1:ny) = P(1:nx,1:ny) + Phi(1:nx,1:ny)
P(1:nx,1:ny) = Pnew(1:nx,1:ny)
</code>```

I am testing my solver on the lid driven cavity.  My velocity profile is correct, however my pressure profile is completely wrong as it contains wiggles.  Here is a picture:

I have tested my solver on a pressure-free projection method(where there is no $\nabla p^{n}$ in the intermediate veloctiy step and it works fine.

EDIT: (1/17/18)

I forgot to mention but I am using 2nd-order central discretization for both convection and diffusion.  My integration technique is currently forward Euler for debugging purposes.


问:
Finding integer/lattice points (coordinates) inside a polytope/polyhedra?
I am using Python but I wouldn't mind changing language. All I have gotten from my research are tools to count the number of (lattice) points inside a region given the equations for the planes that enclose it. Other tools are made to optimize a given function inside the polytope (linear programming).

How about finding the lattice points alone? For example, a function of the kind

```<code>latticePoints( 'x < 5 & x > 0' ) = [ 1, 2, 3, 4]
</code>```

Plus I am looking for something to work in the multivariate scenario (constrains on x, y, z, ...).

I am currently trying to solve this using ppl.


问:
Stability of different quadrature rules in 1st-kind Volterra integral equation
I am dealing with a integral equation
$$ f'(t) = -\int_0^t K(s) f(t-s)\quad t\in [0,t_\max] \tag{1}$$
in which $f(t)$ and $f'(t)$ are known, well-behaved functions of $t$ and $K(t)$ is the unknown. In addition, the physics of the problem dictate that $f(t)$ and $K(t)$ are even about $t=0$. I have discretized this equation on a uniform grid with $n$ intervals of width $\Delta t$. That is, I map
$$
\begin{align}
t &\mapsto i \Delta t \\
t_\max &\mapsto n\Delta t \\
\int_0^t ds &\mapsto \Delta t\sum_{j=0}^i w^{(i+1)}_j
\end{align}
$$
where $w^{(i+1)}$ are the weights for an $(i+1)$-point quadrature rule.
With this, I arrive at a linear system for $K_1,\dots,K_n$:
$$
f_i' + \Delta t w_0^{(i+1)}f_i K_0 = -\Delta t \sum_{j=1}^i w_j^{(i+1)} f_{i-j} K_j\quad i=1,\dots,n \tag{2}
$$
with $K_0 = K(0)$ assumed known. In order to accommodate quadrature rules that require an odd number of points (e.g., Simpson's), I use the evenness of $f$ and $K$ to produce the equivalent discretization
$$
f_i' + \frac{\Delta t}{2} w_i^{(2i+1)}f_i K_0 = -\Delta t \sum_{j=1}^i w_{i-j}^{(2i+1)} f_{i-j} K_j\quad i=1,\dots,n \tag{3}
$$
based on extending the integral to the range $[-t,t]$.

I am testing my numerical solutions by comparing with analytic solutions to Eq. (1). Namely,
$$
f(t) = \cos{(\omega t)}
\implies K(t) = \omega^2
$$
 and
$$
f(t) = e^{-t/2} \left[\cos{(\alpha t/2)} + \frac{\sin{(\alpha t/2)}}{\alpha}\right]
\implies K(t) = \frac{1+\alpha^2}{4}e^{-|t|}
$$
In both cases, I find that the trapezoidal rule produces a numerical solution $K(t)$ that is in excellent agreement with the analytic one. On the other hand Simpson's rule appears to be unstable, producing numerical solutions that oscillate and rapidly grow in amplitude at $t$ increases. A sample for each test case is shown below:

Evidently, there is something about Simpson's rule that makes the linear system in Eq. (3) ill-conditioned. Writing Eq. (3) in standard form, I have
$$
b_i = \sum_{j}A_{ij}x_j
$$
with
$$
b_i=f'_i+\frac{\Delta t}{2}w^{(2i+1)}_if_iK_0\\
A_{ij}=-\Delta t w_{i-j}^{(2i+1)} f_{i-j}\\
x_j=K_j
$$
Since the coefficient matrix $A_{ij}$ is lower-triangular, its eigenvalues are its diagonal entries:
$$\lambda_i = A_{ii} = -\Delta t w_0^{(2i+1)} f_0\quad i=1,\dots,n$$
Whatever the quadrature rule, $w_0^{(2i+1)}$ is the same for all $i\ge 1$. For the trapezoidal rule $w_0=\frac{1}{2}$, and for Simpson's rule $w_0=\frac{1}{3}$. Then in both cases, there is only one unique eigenvalue of the matrix $A$.

Questions

Can I use this information about the eigenvalues to ascertain why Simpson's rule is unstable, while trapezoidal is stable?

Is there some analytic feature of Eq. (1) that could have led me to know a priori that Simpson's rule is unstable?

If so, is there a more general result (or even rules of thumb) for 1st-kind Volterra equations that says certain types quadrature rules are to be preferred?


问:
Can Mallat's pyramidal algorithm be extended to non-power of 2 input sizes?
Mallat's pyramidal algorithm for the discrete wavelet transform operates on power-of-2 vector lengths. Can it be extended to work on inputs of any size without resorting to zero padding?</p>


问:
How to construct Diffusion Wavelet Packets?
I understand the idea of constructing Low-pass and High-pass filters as a projection on the Numerical Range and Numerical Kernel of dyadic powers of a diffusion operator  in the work Diffusion Wavelet Packets.

This can be represented as a tree and is implemented here:
 https://github.com/aweinstein/dw/blob/master/Wavelets/DWPTree.m

I don't understand how to split each $W_{j}$. To my understanding I could make another QR decomposition, as is done for the $V_{j}$ subspace at level $J$.

The authors say (in the end of section 4 of the above mentioned article) to choose each $W_{j}$'s children as having the same dimension.

The problem is: they don't explain why and how to do this splitting, so I started to read the code on github.

In line 404, using a handle to function DefaultSplitFcn, I found a method nowhere explained in articles or lectures about diffusion wavelets.

Using the researches' words (line 503, DefaultSplitFcn function):

Averages the endpoints of the approximate frequency range of the node

My questions are: why can I do this?

Thanks in advance.


问:
What could be causing multi-dimensional numerical integration inconsistency?
I'm trying to numerically integrate a multi-dimensional expression. The integrand is complicated; for example this is the integrand for $N=4$:
$$\begin{aligned}&x_1^6x_2^5x_3^3x_4^2(x_1-x_1x_2)(x_1-x_1x_2x_3)(x_1-x_1x_2x_3x_4)\times{}\\&{}\times(x_1x_2-x_1x_2x_3)(x_1x_2-x_1x_2x_3x_4)(x_1x_2x_3-x_1x_2x_3x_4)\end{aligned}$$

And I'm most interested in expressions that have even larger $N$. The integration limits in my trial case are $x_1$ from 0 to 10, and everything else from 0 to 1.

To numerically integrate this I'm using Mathematica while double checking with the VEGAS integrator in C++ (via the GNU scientific library). The problem is the two integrators match for lower dimensions but not for higher ones. The agreement is great for $N =2, 3$, but for $N=4, 5$ there's a considerable discrepancy of about a factor of 5 that gets worse with higher $N$. Agreement doesn't fail completely; the results are still of the same order of magnitude, but a factor of 5 is too large for comfort especially with both integrators reporting error bars much lower than this.

I am unfamiliar with the internals of both integrators and how numerical integration in higher dimensions works. Is there a (non-bug) reason for why the two integrators give the same results at low dimensions but not at higher dimensions? If the only possible explanation is "a bug in your code" that's fine, but if there's some other possibility I'd like to check that out too.


问:
What is the complexity of calculating K-th largest real part eigenvalue of non-normal sparse matrix
I just need to calculate the largest real part of eigenvalues of a Jacobian which is highly non-normal and singular. Most of the eigenvalues are negative, and some of them are positive but near to zero relatively. I want to find the largest one except the zeroes. And I do not need the eigenvectors.

However, there are several zero eigenvalues. So I want to know if it is worthy to do some transformation to cancel the known zero eigenmode (similar transform to the null space of the modes). Or I can choose to calculate the first few eigenvalues directly.

Here is the sample matlab code to test how much is the cost to calculate first K eigenvalue with largest real part of a random matrix.

```<code>N=3000; % N hundreds to thousand
M=rand(N,N)+eye(N); % shift to avoid zero
K=8; % K is less than 10
t = zeros(K,1);
for k=1:K
    f=@()(eigs(M,k,"largestreal"));
    t(k) = timeit(f);
end
figure
plot(t)
ylabel("time(s)")
xlabel("k")
</code>```

However, the plot output is wavy and non-monotonic.

Is there any analysis on this problem? I read the ARPACK's document, there is no detail of how the it is implemented.


问:
Software for easy spherical coordinate plotting from data file
Does anyone have any suggestions for good open source plotting software which has good graphical spherical coordinate plotting (from a data file). I have a three column data file where column 1 is $\theta$, column 2 is $\phi$ and column three is $f(\theta, \phi)$. I would like to plot this in spherical coordinates $(\theta, \phi, \rho = f(\theta, \phi))$, where points of equal $\rho$ are represented by equivalent colours. Does anyone have any recommendations regarding this type of plotting?

Also, is it possible to easily represent the plot with different colours as a function of the spherical coordinate $\rho$ (hence distance from the origin)?

By the way I used Python to produce these data files.


问:
Methods and tools to solve the two-temperature model (TTM)
I would like to model heat diffusion at the gold / water interface after excitation of the metal surface by an ultrafast laser pulse (ca. 80 fs).

An appropriate model to start with would be the "two temperature model" (TTM) or "parabolic two step" (PTS). This model considers the electronic and lattice systems as separate with electronic and lattice temperature $T_e$ and $T_l$, respectively. The electronic and lattice systems are coupled through an electron-phonon coupling constant $G$ and the electronic system is excited by a source term $S$ representing the incoming laser pulse. The model can be written as two coupled PDEs [from Smith (2001) Appl. Phys. Lett. 78:1240]:
$$
C_e(T_e) \frac{\partial T_e}{\partial t} = \frac{\partial}{\partial x} \left( k_e(T_e, T_l) \frac{\partial T_e}{\partial x} \right) - G \left[ T_e - T_l \right] + S \\
C_l \frac{\partial T_l}{\partial t} = G \left[ T_e - T_l \right],
$$
where $C_e$, $C_l$ are electron and lattice heat capacities, respectively, and $k_e$ is the thermal conductivity.

I am specifically interested in determining the temperatures of the system for the first 100 ps.

The question I would like to ask is: What are the different methods and tools available to solve the problem, with an emphasis on entry-level methodologies?


问:
1D irregular interpolation with $O(1)$ or $O(\log(N))$ evaluation
In Floater's paper on barycentric rational interpolation, he shows that a stable interpolant using irregularly space points can be evaluated in $O(N)$ operations.

For equally spaces samples, cubic b-splines can be used to generate an interpolator that requires $O(1)$ evaluation, and using the Catmull-Rom curve, $O(\log(N))$ operations can interpolate the curve, but the parameterization cannot be supplied.

Is there an algorithm which creates a $C^1$ interpolant from irregular samples (with user-supplied parameterization) which can be evaluated $O(1)$ or $O(\log(N))$ time?


问:
Backward stable algorithm to get orthogonal projection onto the column space of a matrix
I have to find the orthogonal projection of a vector $b$ onto the matrix $A$ of size $m \times n$.

In my application, I don't have the luxury of calculating the QR factorization. All I have are routines that give me the least squares solutions with minimum norm, i.e., ${A^\dagger b}$. I can, of course, use matrix multiplication. There are two ways to get $\mbox{proj}_{A}(b)$, $A (A^{\dagger}b)$ or $(A^T)^{\dagger}(A^Tb)$.

What would be the better strategy for following cases in terms of backward stability?

$m>n$ and A is full rank.

$m>n$ and A is not full rank.

$n>m$ and A is not full rank.

$n=1$.

For the last part, this suggests that there is no backward stable algorithm on the account of the solution having an outer product. I am not convinced, though, because we never explicitly calculate the outer product.


问:
Parallel compact schemes using the Parallal Diagonal Dominant (PDD) algorithm
I would like to use the PDD algorithm developed by Sun to solve tridiagonal matrices in parallel for the following compact finite difference scheme:

$
\begin{align}
\dfrac{1}{4}f^{'}_{i-1} + f^{'}_i + \dfrac{1}{4}f^{'}_{i+1} &= \dfrac{3}{4} \dfrac{-f_{i-1} + f_{i+1}}{dx} \; \; \; \; \;\;\;\;\;\;\;\;\;\; \text{(Interior Points)}
\\
f^{'}_0 + 3f^{'}_1 &= \dfrac{-17 f_0 + 9 f_1 + 9 f_2 - f_3}{6 \Delta x} \; \; \; \text{(Boundary points)}
\end{align}
$

The PDD algorithm consists in splitting the diagonal matrix into a block diagonal, $\tilde{A}$, and offset matrix, $\Delta A$. The following example is for a 8 grid points domain divided into 2 sub-domains (processors):

$
\begin{align}
A &= \begin{bmatrix}
    1       & 3 & & &  \\
    \alpha       & 1 & \alpha \\
    &\alpha       & 1 & \alpha \\
    &&\alpha       & 1 & \alpha \\
    &&&\alpha       & 1 & \alpha \\
    &&&&\alpha       & 1 & \alpha \\
    &&&&& \alpha & 1 & \alpha \\
    &&&&&& 3 & 1
     \end{bmatrix} \\
&= \begin{bmatrix}
    1       & 3 & & &  \\
    \alpha       & 1 & \alpha \\
    &\alpha       & 1 & \alpha \\
    &&\alpha       & 1 &  \\
    &&&        & 1 & \alpha \\
    &&&&\alpha       & 1 & \alpha \\
    &&&&& \alpha & 1 & \alpha \\
    &&&&&& 3 & 1
     \end{bmatrix} +
    \begin{bmatrix}
    0& & & &  \\
    &  0    &  &  \\
    &&   0    &  &  \\
    &&&  0       & \alpha \\
    &&&  \alpha      & 0 &  \\
    &&&&&   0   &   \\
    &&&&&& 0 &    \\
    &&&&&&&  0
    \end{bmatrix} \\
\end{align} \\
\Delta A = [\alpha \, e_4, \alpha \, e_3][e^T_3, e^T_4] = VE^T
$

Where $e_i$ represents a column vector with its ith elements being 1 and zero otherwise for $ 0 < i < n-1$. With this splitting, we have:

$
x^{-1} = A^{-1} d = (\tilde{A} + VE^T)^{-1} d = \tilde{A}^{-1} d - \tilde{A}^{-1}(1+E^T\tilde{A}^{-1}V)^{-1}E^T\tilde{A}^{-1}d
$

This generally leads to have a pentadiagonal reduced matrix solved in parallel. However, the author shows that for a sufficiently large number of grid points per processors, the block-diagonal terms are below machine precision and the matrix reduces into 2x2 independent block problems. This is even more true for diagonally dominant matrices.

The algorithm consists of the following steps:

Solve independently on each $i$ processor: $A_i[\tilde{x}^i, v^i, w^i] = [d^i, \alpha e_0, \alpha e_{m-1}]$

Send $\tilde{x}^i_0, v_0^i$ from the $i^{\text{th}}$ node to the $(i-1)^{\text{th}})$ node.

Solve:
$
\begin{bmatrix}
1 & w_{m-1}^i \\ v_0^{i+1} & 1
\end{bmatrix}
\begin{bmatrix}
y_{2i} \\ y_{2i+1}
\end{bmatrix} =
\begin{bmatrix} \tilde{x}^i_{m-1} \\ \tilde{x}_0^{i+1}\end{bmatrix}
$

Send $y_{2i+1}$ from the $i^{\text{th}}$ to the $(i+1)^{\text{th}}$ node.

$x^i = \tilde{x}^i - \begin{bmatrix}v^i, w^i\end{bmatrix} \begin{bmatrix} y_{2i-1} \\ y_{2i} \end{bmatrix}$

I wrote a MPI fortran code which follows the algorithm. I tested the code with explicit finite difference and also in serial to make sure the parallelism and compact scheme were implemented correctly. In parallel, I am getting large errors right at the nodes boundary. I have verified that my implementation follows the algorithm perfectly. Anyone know where the problem may lie?

Here is the standalone PDD algorithm for the compact scheme:

```<code>subroutine tridiag_parallel(nPoints, dx, phi, dphi, src, des)

    integer,  intent(in) :: nPoints, src, des
    real(dp), intent(in) :: dx
    real(dp), dimension(-buf:), intent(in)  :: phi
    real(dp), dimension(   0:), intent(out) :: dphi
    integer  :: ii, min_x, max_x
    integer  :: recv_rqst, send_rqst
    real(dp) :: a_im, c_ip1m, x_recv, v_recv
    real(dp) :: y_2im1, y_2i, y_2ip1
    real(dp), dimension(0:nPoints-1) :: Ai_low, Ai_mid, Ai_upp, e_m
    real(dp), dimension(0:nPoints-1) :: x_tilde, v, w, rhs, rhs_v, rhs_w

    Ai_low = alpha; Ai_mid = 1._dp; Ai_upp = alpha
    y_2i = 0._dp; y_2ip1 = 0._dp; y_2im1 = 0._dp; x_recv = 0._dp; v_recv = 0._dp

    ! Build RHS
    min_x = 0; max_x = nPoints-1
    if(src == mpi_proc_null) then
        min_x = 1; Ai_upp(0) = 3._dp
        rhs(0) = 1._dp/(6._dp*dx) * (-17._dp*phi(0) + 9._dp*(phi(1)+phi(2)) - phi(3))
    end if

    if(des == mpi_proc_null) then
        max_x = nPoints-2; Ai_low(nPoints-1) = 3._dp

        rhs(nPoints-1) = 1._dp/(6._dp*dx) * (phi(nPoints-4) - 9._dp*(phi(nPoints-3)+phi(nPoints-2)) &
                                        + 17._dp*phi(nPoints-1))
    end if

    do ii=min_x,max_x
        rhs(ii) = 3._dp/(4._dp*dx) * (-phi(ii-1) + phi(ii+1))
    end do

    ! Build a_im, c_im, e_m, rhs_v, rhs_w
    a_im = alpha; c_ip1m = alpha; e_m = 0._dp
    if(src == mpi_proc_null) a_im   = 0._dp
    if(des == mpi_proc_null) c_ip1m = 0._dp

    rhs_v = 0._dp; rhs_v(0)         = a_im
    rhs_w = 0._dp; rhs_w(nPoints-1) = c_ip1m

    ! Solve for x_tilde, v and w
    call tridiag(nPoints, Ai_low, Ai_mid, Ai_upp, rhs,   x_tilde)
    call tridiag(nPoints, Ai_low, Ai_mid, Ai_upp, rhs_v, v)
    call tridiag(nPoints, Ai_low, Ai_mid, Ai_upp, rhs_w, w)

    ! First communication: send tilde{x_0}, v_0 to previous processor
    call mpi_isend(x_tilde(0), 1, mpi_dp, src, 200, mpi_comm_world, send_rqst, ierror)
    call mpi_irecv(x_recv,     1, mpi_dp, des, 200, mpi_comm_world, recv_rqst, ierror)
    call mpi_wait(send_rqst, stat, ierror); call mpi_wait(recv_rqst, stat, ierror)

    call mpi_isend(v(0),   1, mpi_dp, src, 300, mpi_comm_world, send_rqst, ierror)
    call mpi_irecv(v_recv, 1, mpi_dp, des, 300, mpi_comm_world, recv_rqst, ierror)
    call mpi_wait(send_rqst, stat, ierror); call mpi_wait(recv_rqst, stat, ierror)

    ! Solve for y_2i, y_2ip1
    y_2i   = (x_tilde(nPoints-1) - w(nPoints-1)*x_recv) / (1._dp - w(nPoints-1)*v_recv)
    y_2ip1 = x_recv - v_recv * y_2i

    ! Second communication: send y_2i to next processor
    call mpi_isend(y_2i  , 1, mpi_dp, des, 400, mpi_comm_world, send_rqst, ierror)
    call mpi_irecv(y_2im1, 1, mpi_dp, src, 400, mpi_comm_world, recv_rqst, ierror)
    call mpi_wait(send_rqst, stat, ierror); call mpi_wait(recv_rqst, stat, ierror)

    ! Derivative
    do ii=0,nPoints-1
        dphi(ii) = x_tilde(ii) - (v(ii)*y_2im1 + w(ii)*y_2i)
    end do
end subroutine tridiag_parallel
</code>```

Below is the actual program (not sure if necessary but here it is anyway):

```<code>program main
    use, intrinsic :: iso_fortran_env, dp => real64
    use mpi, mpi_dp => mpi_double_precision, stat => mpi_status_ignore
    implicit none

    ! Grid
    integer  :: nx, grid_nx
    integer, parameter :: buf = 2
    real(dp) :: dx, proc_xmin
    real(dp), parameter :: grid_xmin = -1._dp, grid_xmax = 1._dp
    real(dp), allocatable, dimension(:) :: x, x_grid

    ! MPI
    integer :: nProcs, myrank, ierror, Xsrc, Xdes
    logical :: master_proc = .false.

    ! Compact
    real(dp), parameter :: alpha = 1._dp/4._dp

    ! Subscripts => s: serial, p: parallel, e: exact
    real(dp), allocatable, dimension(:) :: u_s, u_p, dudx_s, dudx_p, dudx_e
    integer :: ii, nn, send_rqst, recv_rqst

    ! Init mpi
    call mpi_init(ierror)
    call mpi_comm_size(mpi_comm_world, nProcs, ierror)
    call mpi_comm_rank(mpi_comm_world, myrank, ierror)
    if(myrank == 0) master_proc = .true.

    call init_grid

    ! Serial
    if(master_proc) then
        Xsrc = mpi_proc_null; Xdes = mpi_proc_null
        call setup_case(grid_nx, x_grid, u_s, dudx_s)
        call tridiag_parallel(grid_nx, dx, u_s, dudx_s, Xsrc, Xdes)
        !call FD2_x(grid_nx, dx, u_s, dudx_s, Xsrc, Xdes)
    end if

    ! Parallel
    call init_blocks
    call setup_case(nx, x, u_p, dudx_p)
    call tridiag_parallel(nx, dx, u_p, dudx_p, Xsrc, Xdes)
    !call FD2_x(nx, dx, u_p, dudx_p, Xsrc, Xdes)
    call compute_exact

    ! Compare
    call compare(dudx_s, dudx_p)
    call compute_errors(dudx_e, dudx_p)

    ! Finalize
    call mpi_finalize(ierror)
contains
subroutine init_grid
    ! Read in grid size
    if(master_proc) then
        print'(A)', "Enter the total number of grid points"
        read(*,*) grid_nx
        do nn=1,nProcs-1
            call mpi_isend(grid_nx, 1, mpi_int, nn, 100, mpi_comm_world, send_rqst, ierror)
            call mpi_wait(send_rqst, stat, ierror)
        end do
    else
        call mpi_irecv(grid_nx, 1, mpi_int, 0, 100, mpi_comm_world, recv_rqst, ierror)
        call mpi_wait(recv_rqst, stat, ierror)
    end if

    if(mod(grid_nx, nProcs) /= 0) then
        print'(A)', "The number of points is not divisible by the number of processors."
        print'(A)', "Stopping the simulation"
        call mpi_abort(mpi_comm_world, 10, ierror)
    end if

    dx = (grid_xmax - grid_xmin) / real(grid_nx - 1, dp)
    nx = grid_nx / nProcs
    proc_xmin = grid_xmin + mod(myrank, nProcs) * dx * nx

    allocate(x(0:nx-1), x_grid(0:grid_nx-1))
    x = [ (proc_xmin + ii*dx, ii=0,nx-1) ]
    x_grid = [ (grid_xmin + ii*dx, ii=0,grid_nx-1) ]
end subroutine init_grid

subroutine init_blocks
    Xsrc = myrank-1; Xdes = myrank+1
    if(mod(myrank,   nProcs) == 0) Xsrc = mpi_proc_null
    if(mod(myrank+1, nProcs) == 0) Xdes = mpi_proc_null
end subroutine

subroutine apply_bc(u)
    real(dp), dimension(-2:), intent(inout) :: u
    call mpi_sendrecv(u(0),  2, mpi_dp, Xsrc, 100,  &
                      u(nx), 2, mpi_dp, Xdes, 100, mpi_comm_world, stat, ierror)
    call mpi_sendrecv(u(nx-2), 2, mpi_dp, Xdes, 200, &
                      u(-2),   2, mpi_dp, Xsrc, 200, mpi_comm_world, stat, ierror)
end subroutine apply_bc

subroutine setup_case(nx, x, u, dudx)
    integer,  intent(in) :: nx
    real(dp), dimension(0:), intent(in) :: x
    real(dp), allocatable, dimension(:), intent(out) :: u, dudx

    allocate(u(-2:nx+1), dudx(0:nx-1)); u = 0._dp
    u(0:nx-1) = [(exp(-10._dp*x(ii)**2), ii=0,nx-1)]
    call apply_bc(u)
end subroutine setup_case

subroutine compute_exact
    allocate(dudx_e(0:nx-1))

    dudx_e = [(-20._dp*x(ii)*exp(-10._dp*x(ii)**2), ii=0,nx-1)]
end subroutine compute_exact

! 1st order derivative, 2nd order explicit central difference
subroutine FD2_x(nx, dx, phi, dphi, Xsrc, Xdes)
    integer,  intent(in) :: nx, Xsrc, Xdes
    real(dp), intent(in) :: dx
    real(dp), dimension(-2:), intent(in)    :: phi
    real(dp), dimension( 0:), intent(inout) :: dphi
    integer :: nn, min_x, max_x

    min_x = 0; max_x = nx-1

    ! BC
    if(Xsrc == mpi_proc_null) then
        min_x = 1
        dphi(0) = (phi(1)-phi(0)) / dx
    end if
    if(Xdes == mpi_proc_null) then
        max_x = nx-2
        dphi(nx-1) = (-phi(nx-1)+phi(nx-2)) / dx
    end if

    ! Interior points
    do nn=min_x,max_x
        dphi(nn) = (phi(nn+1)-phi(nn-1)) / (2._dp*dx)
    end do
end subroutine FD2_x

subroutine tridiag(n, a, b, c, r, u)
    integer, intent(in) :: n
    real(dp), dimension(0:), intent(in)  :: a, b, c, r
    real(dp), dimension(0:), intent(out) :: u
    real(dp), dimension(0:n-1) :: gam
    integer  :: j
    real(dp) :: bet

    bet  = b(0)
    u(0) = r(0)/bet

    do j = 1,n-1
        gam(j) = c(j-1) / bet
        bet    = b(j) - a(j) * gam(j)

        u(j) = (r(j) - a(j) * u(j-1))/bet
    end do

    do j =n-2,0,-1
        u(j) = u(j) - gam(j+1) * u(j+1)
    end do
end subroutine tridiag

subroutine compare(serial, parallel)
    real(dp), dimension(0:), intent(in) :: serial, parallel
    integer  :: proc_id
    real(dp) :: Ltwo, Linf
    real(dp), dimension(0:nx-1,0:nProcs-1) :: U_all
    real(dp), dimension(0:grid_nx-1) :: grid_parallel

    if(master_proc) then
        U_all(:,0) = parallel(:)
        do proc_id=1,nProcs-1
            call mpi_recv(U_all(:,proc_id), nx, mpi_dp, proc_id, proc_id+100, mpi_comm_world,   &
                            stat, ierror)
        end do
    else
        call mpi_send(parallel(0:nx-1), nx, mpi_dp, 0, myrank+100, mpi_comm_world, ierror)
    end if
    call mpi_barrier(mpi_comm_world, ierror)

    if(master_proc) then
        grid_parallel = [(U_all(mod(ii,nx), int(ii/nx)), ii=0,grid_nx-1)]

        Ltwo = sqrt(sum((serial-grid_parallel)**2)/real(grid_nx,dp))
        Linf = maxval(abs(serial-grid_parallel))

        print'(A)',"Serial vs Parallel:"
        print'(2(A,ES15.6))', "Ltwo =", Ltwo, "; Linf =", Linf

        open(unit=100,file='postProcessing/compare.dat')
        do ii=0,grid_nx-1
            write(100,*) x_grid(ii), serial(ii), grid_parallel(ii), abs(grid_parallel(ii)-serial(ii))
        end do
        close(100)
    end if
end subroutine compare
subroutine compute_errors(exact, finite)
    real(dp), dimension(0:), intent(in) :: exact, finite
    real(dp) :: m_sum, m_max, Ltwo, Linf

    m_sum = sum((exact-finite)**2)
    call mpi_reduce(m_sum, Ltwo, 1, mpi_dp, mpi_sum, 0, mpi_comm_world, ierror)
    Ltwo = sqrt(Ltwo/real(nx,dp))

    m_max = maxval(abs(exact-finite))
    call mpi_reduce(m_max, Linf, 1, mpi_dp, mpi_max, 0, mpi_comm_world, ierror)

    if(master_proc) then

        print'(A)', "Finite difference vs Exact:"
        print'(2(A,ES15.6))', "Ltwo =", Ltwo, "; Linf =", Linf
    end if
end subroutine compute_errors
end program
</code>```


问:
Is there a software package for Quantum Monte Carlo estimation of the exchange correlation functional?
Quantum Monte Carlo (QMC) calculations historically have been used to parameterise the exchange-correlation functionals for Density Functional Theory (DFT). For example, this article explains a way to do that: Quantum Monte Carlo investigations of density functional theory of the strongly inhomogeneous electron gas.

However, I cannot seem to find any software package that allows one to do this for arbitrary input potentials. Is there some QMC program package that gives this possibility, or are they all in-house software?

I found that the QMCPACK has a vhxc potential option, which might be what I am looking for, but it is undocumented. Thanks!


问:
Convergence of Gauss quadrature for a discontinuous function
Is there a known error estimate for Gaussian quadrature when applied to a discontinuous function?

For simple one-dimensional experiments, the error appears to be bounded by $C h$, where $C$ is some constant and $h = \max |x_{i+1}-x_i|$ is the maximum spacing between two nodes, but I haven't found a clean way to show this.


问:
Speed and accuracy of Strassen vs Winograd matrix multiplication algorithms
I am doing work which requires as fast matrix multiplication as possible and just want to double-check with this community that the Winograd variant of Strassen's MM algorithm is the fastest practical (so no Coppersmith-Winograd) algorithm.

I am doing a lot of processing with the data before, between, and after each MM to the point where using Mathematica or Matlab would be a hindrance.

Also, I am curious if anyone has a good idea for the error in the regular Strassen vs Winograd variant?  In "Exploiting Parallelism in Matrix Computation Kernels for SMP systems" D'Alberto et al. briefly mention Strassen as being more accurate but this seems counter-intuitive since Winograd has less operations overall.

Edit: We are using matrices up to size 2^16 x 2^16 ~ 4 billion doubles, so a sub-cubic algorithm is definitely faster than naive.

Edit 2: On the accuracy of Strassen vs Winograd, if anyone is interested.  In "Accuracy and Stability of Numerical Algorithms" Higham has an in-depth analysis of the error of the two algorithms and shows Strassen has slower error growth w/respect to size of matrices.  Also of note, Strassen more error-prone (against itself) for matrices with all positive entries.

Thanks


问:
Sequential Quadratic Programming for Quadratically Constrained Quadratic Programs
A standard Quadratically Constrained Quadratic Program (QCQP) is of the form:

$$ \underset{x}{minimize} \frac{1}{2}x^TP_{0}x +  q_{0}^{T}x
$$
$$
   subject \; to \quad \frac{1}{2}x^TP_{i}x +  q_{i}^{T}x + r_{i}\leq 0\quad i \in \{1,2,\ldots,m\}
$$
$$
\quad Ax = b
$$

It seems to me that this type of program is specifically well suited for a Sequential Quadratic Programming (SQP) solution, especially when the $P_{i}$ are positive semi-definite. However, I have only seen SDP relaxations in literature.

Are there any commerical software packages that solve these types of problems with SQP? What other methods are actually used in practice?


问:
What is a performant clustering algorithm for approx 10,000 vectors of approx 30 dimension?
I have a set of real-valed vectors, for example

$S = \{v_1, v_2, ..., v_k\}$

$v_i = \begin{pmatrix} age_i \\ height_i \\ weight_i \\ ... \end{pmatrix}$

or whatever. Each vector has on the order of 30 components, and $k\approx 10,000$. The vectors are not normalized.

Two objectives:

count how many clusters are naturally formed by the data set,

given arbitrary $N$, choose $N$ representatives from the set (or construct representative vectors) which span the clusters.

There are lots of clustering algorithms to choose from. I can be flexible on what the definition of clusters entails, but I do want "reasonable" behaviour if $N$ is different from the number of clusters which naturally exist:

if N is too small, the biggest, densest clusters should be represented first,

if N is too big the extra representatives should sit near the big clusters rather than anywhere else, and the extra representatives should not all congregate in the same place.

Simplicity and performance are the priority. The definition of a cluster is negotiable. It would be nice if the algorithm could be tuned to specify that certain components are "important" in defining a cluster, but that's just a wish.

What would be a fast, simple clustering method for this kind of data?


问:
Eigenvalue problem (LAPACK)
I am working on a project in numerical analysis which I have to program in C (using Lapack and Blas). Matrix is given which is tridiagonal and "almost" symmetric (one element is to be changed to make it symmetric). In order to calculate certain eigenvalues and eigenvectors I am planning to use dsyevx() routine which works only for symmetric matrices. Is there any routine for orthogonal similarity transformation to make it symmetric, or should I opt for different routine for calculating eigenvalues of general $N\times N$ matrix?

Title: Eigenvalue problem (LAPACK)


问:
How can I solve the wave equation for a circular rod in cylindrical coordinates using finite differences?
I have a problem with the stability of finite difference method for the wave equation in cylindrical coordinates.

the equation is:
$$
\frac{\partial^2 \omega_n}{\partial r^2}+\frac{1}{r}\frac{\partial \omega_n}{\partial r}+\frac{\partial^2\omega_n}{\partial z^2}-\frac{\omega_n}{r^2}=\frac{1}{c^2}\frac{\partial^2\omega_n}{\partial t^2}
$$

and when I run it in Matlab, I get the following result:

Does anyone know where the problem is?

I put the code here so that if anyone knows where the problem is, help me out.

```<code>%domain
clear;
Lx=10;
Ly=10;
dx=0.1;
dy=dx;
nx=fix(Lx/dx);
ny=fix(Ly/dy);
x=linspace(0,Lx,nx);
y=linspace(0,Ly,ny);
%time
T=10;
%varables
wn=zeros(nx,ny);
wnm1=wn; %w at time n-1
wnp1=wn; %w at tome n+1
%parameters
CFL=0.5; %CFL=c.dt/dx
c=1;
dt=(CFL*dx/c);

%dt=1e-6;

%%initial conditions
%%time stepping loop
t=0;
while(t<T)
    %reflecting boundary conditions
    wn(:,[1 end])=0;
    wn([1 end],:)=0;

    %solution
    t=t+dt;
    wnm1=wn; wn=wnp1; %save current and previous arrays

    %sourse
    wn(50,50)=dt^2*20*sin(30*pi*t/20);
    for i=2:nx-1 for j=2:ny-1

           wnp1(i,j)=2*wn(i,j)-wnm1(i,j)...
               +CFL^2*(((wn(i+1,j)-2*wn(i,j)+wn(i-1,j))/dx^2+(1/(2*i*dx^2))*(wn(i+1,j)-wn(i-1,j))+(wn(i,j+1)-2*wn(i,j)+wn(i,j-1))/dy^2-wn(i,j)/(i*dx)^2));

        end,end
    %chek convergence

    %visualize at selected steps
    clf;
    subplot(2,1,1);
    imagesc(x,y,wn'); colorbar;caxis([-0.02 0.02])
    title(sprintf('t=%.2f',t));
    subplot(2,1,2);
    mesh(x,y,wn'); colorbar;caxis([-0.02 0.02])
    axis([0 Lx 0 Ly -0.05 0.05]);
    shg; pause(0.01);
end
</code>```


问:
Meaning of "points where the first- and second-order optimality conditions fail"
I am reading a paper that summarizes a set of simulations. Basically, the authors are trying to minimize some function using different optimization algorithms. They conclude: "Our findings point to instances of convergence at points where the first- and second-order optimality conditions fail."

What does this mean in plain English -- that parameter updates have essentially stopped even though the algorithm has not identified an optimal point? My understanding is:

First-order conditions: Ensure that a critical point has been found -- maximum, minimum, or inflection.

Second order conditions: Determine that this point is of the nature we want (e.g. that it is indeed a minimum).

Is this correct? My knowledge of optimization is a bit rusty, so I just wanted to confirm that I understand.


问:
surface triangulation to boundary representation
Are there any open source tools available which can take a surface triangulation of mesh volumes and convert them to a boundary representation?

I am interested in refining multiple region tetrahedral meshes with boundary conditions applied to some of the surfaces.  I have attempted to read a mesh into gmsh, but it is unclear what steps would be necessary to convert it to a boundary representation.


问:
Validating a Markov chain and Bagging to find overall model
let´s say I estimated a Markov chain by splitting my data into traning and test data, and estimating the transition matrix for the traning data yields
$$ P=
    \begin{bmatrix}
    p11 & p12 & p13 \\
    p21 & p22 & p23 \\
    p31 & p32 & p33 \\
    \end{bmatrix}
$$

i.e, the transitions between states 1, 2 and 3. It is, of course, possible to obtain the n-step ahead matrix from this.

Question 1: How do I evaluate the predictive capacity of this model on the test data? I know I should evaluate the MSE of this matrix on the transitions that actually occurred on the test, but how to do it? The caret package in R has some code, but only for other types of predictors (Knn, naive Bayes, etc, in any case, nothing whose output is a matrix). An easy approach would be to obtain the observed transitions matrix for the test data, divide by rows, and then get the element-wise difference with P, but thats probably not the best approach...

Question2 (Related to question 1): I am unable to deploy a model in my entire dataset, due to data size. Is it possible to evaluate several Markov models on samples and then take some kind of aggregation of these to better reflect the overall data? Kind of what bagging is about, but once again, no examples exist for Markov chains, only decision trees, etc.


问:
Nonlinear Sylvester-Like Equation
Maybe you can point me to some results already developed for this.

I have to solve for $X$ the following "Sylvester-like" equation:

$$ AX - XB = F(X)$$

where $A\in\mathbb{R}^{a\times n}$, $B\in\mathbb{R}^{p\times b}$, $X\in\mathbb{R}^{n\times p}$ and $F(X)$ is a nonlinear function of the elements of X.

I know I can solve this using generic numerical algorithms to solve for the roots of $ AX - XB - F(X) = 0$. Nevertheless, I would like to know if there is any algorithm or procedure that actually exploits the "Sylvester-like" structure of the equation.

Thanks in advance!

Title: Nonlinear Sylvester-Like Equation