Title: Calculating the exponential of a complex matrix

Body:
I am trying to calculate the exponent of a 3 x 3 matrix using the formula

$\sum_{i=0}^\infty\frac {A^n}{n!}$

I believe that my error may lay in the scalar division with a factorial or the member function performing the actual exponential calls. After a couple of iterations, I see my output values explode for the simple case of having the matrix in question filled with 1 + i elements.

Output
    This program will calculat

```<code>> Blockquotee the exp of a matrix A
Calculating power of matrix
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
Performing scalar division
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)

Performing Summation
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
1

Calculating power of matrix
(0,6)(0,6)(0,6)
(0,6)(0,6)(0,6)
(0,6)(0,6)(0,6)
Performing scalar division
(0,6)(0,6)(0,6)
(0,6)(0,6)(0,6)
(0,6)(0,6)(0,6)

Performing Summation
(1,7)(1,7)(1,7)
(1,7)(1,7)(1,7)
(1,7)(1,7)(1,7)
2

Calculating power of matrix
(-18,18)(-18,18)(-18,18)
(-18,18)(-18,18)(-18,18)
(-18,18)(-18,18)(-18,18)
Performing scalar division
(-18,18)(-18,18)(-18,18)
(-18,18)(-18,18)(-18,18)
(-18,18)(-18,18)(-18,18)

Performing Summation
(-17,25)(-17,25)(-17,25)
(-17,25)(-17,25)(-17,25)
(-17,25)(-17,25)(-17,25)
3

Calculating power of matrix
(-108,0)(-108,0)(-108,0)
(-108,0)(-108,0)(-108,0)
(-108,0)(-108,0)(-108,0)
Performing scalar division
(-108,0)(-108,0)(-108,0)
(-108,0)(-108,0)(-108,0)
(-108,0)(-108,0)(-108,0)

Performing Summation
(-125,25)(-125,25)(-125,25)
(-125,25)(-125,25)(-125,25)
(-125,25)(-125,25)(-125,25)
4

Calculating power of matrix
(-324,-324)(-324,-324)(-324,-324)
(-324,-324)(-324,-324)(-324,-324)
(-324,-324)(-324,-324)(-324,-324)
Performing scalar division
(-324,-324)(-324,-324)(-324,-324)
(-324,-324)(-324,-324)(-324,-324)
(-324,-324)(-324,-324)(-324,-324)

Performing Summation
(-449,-299)(-449,-299)(-449,-299)
</code>```

(-449,-299)(-449,-299)(-449,-299)
(-449,-299)(-449,-299)(-449,-299)
5

```<code>Calculating power of matrix
(0,-1944)(0,-1944)(0,-1944)
(0,-1944)(0,-1944)(0,-1944)
(0,-1944)(0,-1944)(0,-1944)
Performing scalar division
(0,-1944)(0,-1944)(0,-1944)
(0,-1944)(0,-1944)(0,-1944)
(0,-1944)(0,-1944)(0,-1944)

Performing Summation
(-449,-2243)(-449,-2243)(-449,-2243)
(-449,-2243)(-449,-2243)(-449,-2243)
(-449,-2243)(-449,-2243)(-449,-2243)
6

Calculating power of matrix
(5832,-5832)(5832,-5832)(5832,-5832)
(5832,-5832)(5832,-5832)(5832,-5832)
(5832,-5832)(5832,-5832)(5832,-5832)
Performing scalar division
(5832,-5832)(5832,-5832)(5832,-5832)
(5832,-5832)(5832,-5832)(5832,-5832)
(5832,-5832)(5832,-5832)(5832,-5832)

Performing Summation
(5383,-8075)(5383,-8075)(5383,-8075)
(5383,-8075)(5383,-8075)(5383,-8075)
(5383,-8075)(5383,-8075)(5383,-8075)
7

Calculating power of matrix
(34992,0)(34992,0)(34992,0)
(34992,0)(34992,0)(34992,0)
(34992,0)(34992,0)(34992,0)
Performing scalar division
(34992,0)(34992,0)(34992,0)
(34992,0)(34992,0)(34992,0)
(34992,0)(34992,0)(34992,0)

Performing Summation
(40375,-8075)(40375,-8075)(40375,-8075)
(40375,-8075)(40375,-8075)(40375,-8075)
(40375,-8075)(40375,-8075)(40375,-8075)
8

Calculating power of matrix
(104976,104976)(104976,104976)(104976,104976)
(104976,104976)(104976,104976)(104976,104976)
(104976,104976)(104976,104976)(104976,104976)
Performing scalar division
(104976,104976)(104976,104976)(104976,104976)
(104976,104976)(104976,104976)(104976,104976)
(104976,104976)(104976,104976)(104976,104976)

Performing Summation
(145351,96901)(145351,96901)(145351,96901)
(145351,96901)(145351,96901)(145351,96901)
(145351,96901)(145351,96901)(145351,96901)
9
</code>```

```<code>#pragma once
#include <iostream>
#include <complex>
#include <cmath>
#include <cassert>

using namespace std;

class ComplexMatrix
{
    private:
    complex<long double>** Arr;
    int mi = 3;
    int mj = 3;
public:
    ComplexMatrix();
    //~ComplexMatrix();
    ComplexMatrix(int i, int j);
    ComplexMatrix(const ComplexMatrix&);
    void Initialise(complex<long double>);
    void DisplayMatrix();
    void DeleteMatrix();
    void EnterComplexMatrix(int, int);

    ComplexMatrix matrixPower(ComplexMatrix&, int);
    ComplexMatrix ScalarDivision_Fac(ComplexMatrix& , complex<long
    double>);
    ComplexMatrix MatrixAddition(ComplexMatrix&, ComplexMatrix&);

    ComplexMatrix matrixSq(ComplexMatrix&);
    ComplexMatrix Multiply(ComplexMatrix, ComplexMatrix);

    ComplexMatrix matrixe_A(ComplexMatrix&, int);

    ComplexMatrix operator*(const ComplexMatrix&);
    ComplexMatrix operator=(const ComplexMatrix&);

    friend ComplexMatrix operator+(const ComplexMatrix&, const ComplexMatrix&);
};

//Constructor to initialise a default 3 x 3 complex matrix with 0s for real and imaginary values
ComplexMatrix::ComplexMatrix()
{
    Arr = new complex<long double> * [mi];

    for (int x = 0; x < mi; x++)
    {
        Arr[x] = new complex<long double>  [mj];
    }

    for (int x1 = 0; x1 < mi; x1++)
    {
        for (int y1 = 0; y1 < mj; y1++)
        {
            Arr[x1][y1] = (0.0, 0.0);
        }
    }
}

//Constructor to initialise a user defined complex matrix of a particular size with 0s for real
and imaginary values
ComplexMatrix::ComplexMatrix(int i, int j)
{
mi = i;
mj = j;
Arr = new complex<long double> * [i];

for (int x = 0; x < i; x++)
{
    Arr[x] = new complex<long double>[j];
}
for (int x1 = 0; x1 < i; x1++)
{
    for (int y1 = 0; y1 < j; y1++)
    {
        Arr[x1][y1] = (0.0, 0.0);
    }
    }
}

//Copy Constructor
ComplexMatrix::ComplexMatrix(const ComplexMatrix& CM)
//lets only have ARR mean one thing to make it easier to read, understand,
and to avoid this->
clutter.
{
    Arr = new complex<long double> * [mi];

    for (int x = 0; x < mi; x++)
    {
        Arr[x] = new complex<long double>[mj];
    }

for (int x1 = 0; x1 < mi; x1++)
{
    for (int y1 = 0; y1 < mj; y1++)
    {
        Arr[x1][y1] = CM.Arr[x1][y1];
    }
    }
}
/*
ComplexMatrix::~ComplexMatrix()
{
    for (int x = 0; x < mi; x++)
    {
        delete[] Arr[x];
    }
    delete[] Arr;
}
*/
//Initialise matrix elements to a particular value
void ComplexMatrix::Initialise(complex<long double> x)
{
    for (int i = 0; i < mi; i++)
    {
        for (int j = 0; j < mj; j++)
        {
            Arr[i][j] = x;
        }
    }
}

//Display the matrix member function
void ComplexMatrix::DisplayMatrix()
{
    for (int x = 0; x < mi; x++)
    {
        for (int y = 0; y < mj; y++)
        {
            cout << Arr[x][y];
        }
        cout << endl;
    }
}

//Delete the memory allocated by the matrix member function
void ComplexMatrix::DeleteMatrix()
{
    for (int x = 0; x < mi; x++)
    {
        delete[] Arr[x];
    }
    delete[] Arr;
}

//Enter complex matrix elements member function
void ComplexMatrix::EnterComplexMatrix(int i, int j)
{
    double real, img;
complex < long double> temp = (0.0, 0.0);
cout << "Your matrix will have " << i * j << " elements" << endl;

//Prompt for user input and assign values for real and imaginary values
for (int x = 0; x < i; x++)
{
    for (int y = 0; y < j; y++)
    {
        cout << "Enter the details for the real part of element" << "[" <<
 x << "]" << "[" << y
        << "]" << endl;
        cin >> real;
        cout << "Enter the details for the real part of element" << "[" <<
 x << "]" << "[" << y
        << "]" << endl;
        cin >> img;
        temp = (real, img);
        Arr[x][y] = temp;
        }
    }
}

ComplexMatrix ComplexMatrix::Multiply(ComplexMatrix x, ComplexMatrix y)
{
    ComplexMatrix z(3, 3);

    for (int x1 = 0; x1 < 3; ++x1)
    {
        for (int y1 = 0; y1 < 3; ++y1)
        {
            for (int z1 = 0; z1 < 3; ++z1)
            {
                Arr[x1][y1] += x.Arr[x1][z1] * y.Arr[z1][y1];
            }
        }
    }
    return z;
}

ComplexMatrix ComplexMatrix::ScalarDivision_Fac(ComplexMatrix& x,
complex<long double> n)
{
ComplexMatrix newCompArr(3, 3);
complex <long double> fac = 0.0;
int n1 = static_cast <int>(n.real());
n1 = static_cast <int>(n1);
complex <long double> i1;

for (int i = 1; i < n1; i++)
{
    i1 = i;
    fac = fac * i1;
}

for (int x1 = 0; x1 < mi; x1++)
{
    for (int y1 = 0; y1 < mj; y1++)
    {
        newCompArr.Arr[x1][y1] = x.Arr[x1][y1] / fac;
    }
}

return newCompArr;
}

ComplexMatrix ComplexMatrix::matrixSq(ComplexMatrix& x)
{
    ComplexMatrix result(mi, mj);
    result = x * x;
    return result;
}

ComplexMatrix ComplexMatrix::matrixPower(ComplexMatrix& a, int n)
{
    ComplexMatrix result(mi, mj);
    ComplexMatrix temp(mi, mj);

    temp = a;

    if (n % 2 == 0)
    {
        for (int i = 1; i < n / 2; i++)
        {
            result = temp * a;
            temp = result;
        }
        result = temp;
        result = result.matrixSq(result);
    }
    else
    {
        for (int j = 0; j < (n - 1) ; j++)
        {
            result = temp * a;
            temp = result;
        }
        result = temp;
    }
    return result;
}

ComplexMatrix ComplexMatrix::matrixe_A(ComplexMatrix& A, int n)
{
    ComplexMatrix expA(mi, mj);
    ComplexMatrix sum(mi, mj);
    sum.Initialise({ 0.0, 0.0 });
    ComplexMatrix A_n(mi, mj);
    ComplexMatrix A_n_div_n(mi, mj);
    ComplexMatrix temp(mi, mj);
    ComplexMatrix zero(mi, mj);
    zero.Initialise({ 0.0, 0.0 });
    complex <long double> j;

    for (int i = 1; i < n; i++)
    {
        cout << "Calculating power of matrix" << endl;
        A_n = A.matrixPower(A, i);
        A_n.DisplayMatrix();

        A_n_div_n = A_n;
        cout << "Performing scalar division" << endl;
        A_n_div_n.ScalarDivision(A_n_div_n, i);
        A_n_div_n.DisplayMatrix();
        cout << endl;

        temp = zero;
        temp = A_n_div_n;
        cout << "Performing Summation" << endl;
        sum = sum + temp;
        sum.DisplayMatrix();

        cout << i << endl;
        cout << endl;

        temp = zero;
        A_n = A.matrixPower(A, i);
    }
    return sum;
}

ComplexMatrix ComplexMatrix::operator*(const ComplexMatrix& CompArr)
{
    ComplexMatrix newCompArr(3, 3);

    for (int x1 = 0; x1 < 3; ++x1)
    {
    for (int y1 = 0; y1 < 3; ++y1)
        {
            newCompArr.Arr[x1][y1] = {0.0, 0.0};
            for (int z1 = 0; z1 < 3; ++z1)
            {
                newCompArr.Arr[x1][y1] += Arr[x1][z1] * CompArr.Arr[z1][y1];
            }
        }
    }
    return newCompArr;
}

ComplexMatrix ComplexMatrix::operator=(const ComplexMatrix& CM)
{
    mi = 3;
    mj = 3;

    //ComplexMatrix Arr(3, 3);

    for (int x1 = 0; x1 < mi; x1++)
    {
        for (int y1 = 0; y1 < mj; y1++)
        {
            Arr[x1][y1] = CM.Arr[x1][y1];
        }
    }
    //return Arr;
    return *this;
}

ComplexMatrix operator+(const ComplexMatrix &x, const ComplexMatrix &y)
{
    int mi = 3;
    int mj = 3;
    ComplexMatrix newCompArr(3, 3);
    for (int x1 = 0; x1 < mi; x1++)
    {
        for (int y1 = 0; y1 < mj; y1++)
        {
            newCompArr.Arr[x1][y1] = {  (x.Arr[x1][y1].real() + y.Arr[x1][y1].real()) ,
                                    (x.Arr[x1][y1].imag() + y.Arr[x1][y1].imag()) };
        }
    }
    return newCompArr;
}

#include <iostream>
#include "Header.h"
#include <complex>
#include <cmath>

using namespace std;

int main()
{
    std::cout << "This program will calculate the exp of a matrix A\n";
   complex<long double> x = {1, 1};
    complex<long double> y = {2, 2};
    complex<long double> z = { 0.0, 0.0 };

    ComplexMatrix z1(3, 3);
    ComplexMatrix z2(3, 3);
    ComplexMatrix z3(3, 3);

    z2.Initialise(x);
    z3 = z2.matrixe_A(z2, 10);
    //z3.DisplayMatrix();

    z1.DeleteMatrix();
    z2.DeleteMatrix();
}
</code>```


问:
Well posedness dirichlet BC in elasticity
I always thought that a Dirichlet BC need to be prescribed on at least one part of the domain in order for linear elasticity to be well-posed. Is this correct?

But I was recently told by a colleague that this is not true and that the problem can be well-posed by specifying 2 tractions.

ContentLicense: CC BY-SA 4.0


问:
Assume $AX = C$. How to determine which entry of $BX - D$ is non-negative?
Let $A,B$ be $n \times n$ matrices and $C,D$ be $n \times 1$ matrices. Moreover, all entries of $A,B,C,D$ are non-negative. Assume that there is a unique matrix $X$ that solves $AX = C$.

My goal is not to calculate $X$, but to determine which entry of $BX - D$ is non-negative. In solving my problem, I need to repeat this procedure several times.

I would like to ask if there is an efficient method (or references) to do so. Thank you so much!


问:
Multidimensional center-of-mass coordinates for PDE solutions
In this question I'm looking for a discussion on coordinate transformations. Say I have a PDEs of the form
$$
\frac{\partial}{\partial t} \, V(t,x) \ = \ \Bigg\{r(x,t) - \sum_{m=1}^D \mu_m(t,x)\frac{\partial}{\partial x_m} - \frac{1}{2} \sum_{m=1}^D\sum_{n=1}^D\Gamma_{mn}(t,x) \frac{\partial^2}{\partial x_m \partial x_n}\,\Bigg\} \; V(t,x) \,,
$$
and I want to transform the coordinates $\{x_i\}_{i=1 \ldots D}$ to what I'd call generalized center-of-mass coordinates. That is, in the new coordinates $\{y_j\}$, the first dimension should be $y_1 = \frac{1}{D}\sum_i x_i$.

As far as I remember, there is a Jacobi matrix $J$ which connects the two coordinate systems,
$$
\begin{pmatrix}y_1 \\ \vdots \\ y_D\end{pmatrix} = J \begin{pmatrix}x_1 \\ \vdots \\ x_D\end{pmatrix}
$$
For center-of-mass coordinates, the first row of $J$ should consist of constants, $J_{1i} = \frac{1}{D}(1,\ldots,1)$. With this, there are $D-1$ rows left to choose, with the constraint that $J$ should have full rank and usually also be orthogonal. Moreover, I would choose the Jacobian as a matrix with constant entries (unlike in spherical coordinates, say, where the Jacobian depends on the position).

Questions:

What are standard approaches for generalized center-of-mass coordinates?

Do these approaches fit in the simple framework above? If yes, what is the strategy for assembling the other coordinates $y_2, \ldots, y_D$ beside the center of mass $y_1$?

Is there a systematic way to find a Jacobi matrix such that the PDE above remains "as easy as possible"?

As a quick idea, one could simply choose the rows $2$ to $D$ with random constants, and then apply Gram-Schmidt to obtain an orthogonal matrix -- any comments on that?

Please feel free to answer only parts of these questions and maybe also to provide an alternative view.


问:
CUDA and implicit array expansion
I am retrying this question since I didn’t do so properly last time:

I have implement a simply c++ CUDA example that takes A [m x 1] and B [1 x N] and adds them element wise with what I am used to being called “implicit” array expansion (MATLAB term). I have two questions:

Why do I not run into errors with two threads simultaneously reading from the same element of a or b?

Is this the “proper” way to perform these types of element wise operations where “element add” represent any number of custom functions that take any permutations of row and column matrixes and calculate on all permutations of combinations?

```<code>#include <stdio.h>

#define Row  500
#define Col  60000

__global__ void elementadd( int * a, int * b, int * c )
{

  int ro = blockIdx.x * blockDim.x + threadIdx.x;
  int co = blockIdx.y * blockDim.y + threadIdx.y;

  if (ro < Row && co < Col)
  {
    c[ro*Col + co] = a[ro] + b[co];
  }
}

int main()
{
  int *a, *b, *c_cpu, *c_gpu;

  float s1 = (float)(Row) * (float)(Col);
  float size1 = s1 *  sizeof (int); // Number of bytes of an N x N matrix
  int sizeR = Row * sizeof(int);
  int sizeC = Col * sizeof(int);

  // Allocate memory
  cudaMallocManaged (&a, sizeR);
  cudaMallocManaged (&b, sizeC);
  cudaMallocManaged (&c_cpu, size1);
  cudaMallocManaged (&c_gpu, size1);

//initialize 2d
for(int i = 0; i < Row; ++i){
  a[i] = i*2;
}

for(int i = 0; i < Col; ++i){
b[i] = i+2;
}

for(int ro = 0; ro < Row; ++ro){
for(int co = 0; co < Col; ++co){
c_gpu[ro*Col + co] = 0;
c_gpu[ro*Col + co] = 0;
}
}

  dim3 threads_per_block (32, 32, 1); // A 16 x 16 block threads
  dim3 number_of_blocks ((Row / threads_per_block.x) + 1, (Col / threads_per_block.y) + 1, 1);

  elementadd<<< number_of_blocks, threads_per_block >>> ( a, b, c_gpu );

  cudaDeviceSynchronize(); // Wait for the GPU to finish before proceeding

for(int ro = 0; ro < Row; ++ro){
for(int co = 0; co < Col; ++co){
printf("entry[%d][%d] is %d a:%d b:%d \n",ro,co,c_gpu[ro*Col + co],a[ro],b[co]);
}
}

  // Free all our allocated memory
  cudaFree(a); cudaFree(b);
  cudaFree( c_cpu ); cudaFree( c_gpu );
}
<span class="math-container">```
```

Title: CUDA and implicit array expansion


问:
Is there some algorithm to find the shortest path in the context of genetics and breeding?
In genetics and breeding, we typically have two parent genotype (may or may not be the same) which can produce a set of offspring with certain probability (assuming simple Mendelian inheritance).

I am looking for a way to determine given a set of parent, what is the shortest path to a particular offspring genotype (and yes, I am talking about breeding flowers in Animal Crossing!).

This seems to be very similar to the shortest path problem, but now the edges is given by two nodes interacting with each other. Is there a way to achieve this?

I was directed to here from stack overflow. My original post was at stackoverflow.

Edit: A genotype is a word of four groups of two chars. These groups can be xx, Xx or XX where x = {r,b,y,w}. For example; RRbbyyWw is a red flower.

An offspring of two parents inherits a random combination of letters. For example, if parent A has genotype RRbbyyWw and parent B has genotype rrBBYyww then the offspring might be RRBbyyWw.

If, by luck, the offspring has inherited rRBbyyWw (rR is not allowed) then it is reordered so it becomes RrBbyyWw.

Now, the question is; given genotypes of a set of parents and a target genotype, is there an algorithm which can generate the shortest path in terms of pairing parents?

Edit2: I am aware of the data mining done by others (e.g. here and here.) Therefore, I am not looking for a particular path to get to certain result, rather a method to find a (shortest) path.


问:
Moller Plesset Theory- MP2 correction
I have a molecule LiBH2 with C2V point group. I used 3-21G basis set and obtained RHF coefficients and symmetry of orbitals. I wish to find the set of orbitals (a,b,i,j)( where a and b are occupied and i and j are virtual oritals) that contribute to MP2 correlation energy. That is, in the process I have to find the non vanishing integrals. I tried finding and tried seeing what combination of MOs gives me A1 symmetry, because that will not vanish. However I am getting 156 such combinations which should not be the case. Can someone please help? I have been trying for really long now.

I have attached the output file for reference.

Thanks a lot in advance. Also, you feel this question is off-topic then please let me know where should I post this cause I have already tried it on chemistry.stackexchange.com

Title: Moller Plesset Theory- MP2 correction


问:
Create random number generator from a top hat distribution
How would I use a random number generator to give random numbers from a top hat distribution in Fortran?
I've seen lots of information about how to make a Gaussian random number generator, but not one with a top hat distribution.</p>


问:
To use the confluent hypergeometric function or not to?
I am numerically computing the following integral as a function of positive $k$.
$$I(k) := \int_0^\infty x^b(k+x)^{a-1} e^{-x} dx \tag1$$
It is shown in math.stackexchange.com that this can be expressed as the confluent hypergeometric function of the second kind
$$I(k) = k^{a+b} \Gamma(b+1)\, U(b+1, a+b+1, k)  \tag2$$
under the condition $k>0, \Re(b)>-1$.

For $k\to 0^+$, $I(k)\to \Gamma(a+b)$, while $U(b+1,a+b+1,k)\sim\frac1{k^{a+b}}$. So there is a problem of multiplying a small number with a large number to get a finite number, not a numerically ideal situation. If I use the Python package Scipy, is it more accurate to evaluate the integral $(1)$ directly with Scipy.integrate.quad, or use the confluent hypergeometric function module scipy.special and compute $U$ via expression $(2)$? Or should I divide the axis of $k$ into two parts $[0, K]$ and $[K,\infty)$ and evaluate $I(k)$ on the left interval via $(1)$ and the right interval via $(2)$? If so, how should I choose the best $K$?


问:
Name of method for quickly sampling many small-probability booleans
I have a list of probabilities $p_k$, each corresponding to a Bernoulli distribution. These values of $p_k$ are all small (let's say 0.1% or less). I want to prepare a data structure that allows quickly sampling each of the Bernoulli distributions and returning the sparse set of them that returned 1 instead of 0.

Naively, the expected runtime cost of this sampling task would be $O(N)$ where $N$ is the number of probabilities. I want the expected runtime of producing another set of sparse samples to instead be like $\tilde{O}(pN)$ where $p$ is the maximum probability $p_k$.

I know that this sort of thing can be done, I just can't remember the name of the method. The method involves converting the probabilities into spans within an exponential distribution, creating a search tree over the spans, then using the tree to quickly identify which span a sample from that exponential distribution landed inside (which tells you the next value in the sparse sample), then recursing on the spans past the one that was landed in.

(An alternative method would be to group the probabilities into groups of size $1/p$, precomputing the probability distribution for the number of True's that appear when sampling all elements within the group, then using alias sampling to sample how many samples-without-repetition to extract [also using alias sampling].)

It would be particularly useful if a method achieving this kind of speedup was available in a common python library such as scipy.


问:
Simulating the response of nonlinear system with stiff differential equations
I want to simulate the response of a nonlinear system given in the following form:

$$ \dot{x_1} = f_1(\bar{x_1})+g_1(\bar{x_1})x_2, \ x_1(0) = 0.2 $$
$$ \dot{x_2} = f_2(\bar{x_2})+g_2(\bar{x_2})x_3, \ x_2(0) = 0.5 $$
$$ \dot{x_2} = f_3(\bar{x_3})+g_3(\bar{x_3})u, \ x_3(0) = 0.5 $$
$$ y = x_1 $$

where $u$ is the control input to the system, $y$ is the outputof the systen and the functions $f_i$ and $g_i$, $i = 1,2,3$ are known and nonlinear. In order to come up with the control input at each time step some computations have to be carried out, which are the following:

$$ \epsilon_1 = S_1^{-1}(\frac{x(1)-y_d}{p_1(t)}) $$
$$ n_1 = k_1[\frac{1}{p_1(t)\frac{dS_1}{d\epsilon_1}(\epsilon_1)}+p_1(t)\frac{dS_1}{d\epsilon_1}(\epsilon_1)]\epsilon_1 $$
$$ \dot{ζ_1} = \frac{\epsilon_1n_1}{p_1(t)\frac{dS_1}{d\epsilon_1}(\epsilon_1)} $$
$$ a_1 = N(ζ_1)n_1 $$

$$ \epsilon_2 = S_2^{-1}(\frac{x(2)-a_1}{p_2(t)}) $$
$$ n_2 = k_2[\frac{1}{p_2(t)\frac{dS_2}{d\epsilon_2}(\epsilon_2)}+p_2(t)\frac{dS_2}{d\epsilon_2}(\epsilon_2)]\epsilon_2 $$
$$ \dot{ζ_2} = \frac{\epsilon_2n_2}{p_2(t)\frac{dS_2}{d\epsilon_2}(\epsilon_2)} $$
$$ a_2 = N(ζ_2)n_2 $$

$$ \epsilon_3 = S_3^{-1}(\frac{x(3)-a_2}{p_3(t)}) $$
$$ n_3 = k_3[\frac{1}{p_3(t)\frac{dS_3}{d\epsilon_3}(\epsilon_3)}+p_3(t)\frac{dS_3}{d\epsilon_3}(\epsilon_3)]\epsilon_3 $$
$$ \dot{ζ_3} = \frac{\epsilon_3n_3}{p_3(t)\frac{dS_3}{d\epsilon_3}(\epsilon_3)} $$
$$ u = N(ζ_3)n_3 $$

The functions $S_i(\epsilon_i)$, $i=1,2,3$ are some smooth bijective functions, which are known and so their inverse functions and their derivatives can be computed. Same goes for the functions $p_i(t)$, $i=1,2,3$ which are some performance functions. Finally, the function $N(ζ)$ is an even Nussbaum function, also known and $y_d(t)$ is a desired trajectory for the output of the system. The numerical solution of the closed loop (with $u$ as the input) is stiff and so I will use the ode15s MATLAB solver for the simulation. Now, my problem is how can I compute the values of the $ζ_i$, $i=1,2,3$ variables in order to compute the value of $u$. Some pseudocode I have thought off is:

```<code>function dx = odefun(t,x)

    compute ε1, p1, n1, ζ1
    compute a1

    compute ε2, p2, n2, ζ2
    compute a2

    compute ε3, p3, n3, ζ3
    compute u

    dx = [x1_dot ; x2_dot ; x3_dot]

end
</code>```

I know I can compute the values of the variables $\epsilon_i, n_i, p_i$ by using some function handles but know comes my problem regarding the variables $ζ_i$ due to the fact that they are computed though the differential equations. Is there a way to include them in the ode solver and solve them all at once in the proper way ?


问:
Two-dimensional modelling of a flate-plate reactor

I am trying to simulate the unsteady pollutant concentration along the reactor by solving the 2nd order PDE below with the stated BC and IC. which method is appropriate to solve convection-diffusion-reaction equation (unsteady state) using MATLAB. I solved the steady state using the function pdepe.
Note: I have no knowledge of finite volume/element methods to solving PDE's. I can't seem to recognize why the upwind scheme and second order central difference were used to discretize the convection and diffusion terms respectively. However,I have some knowledge about finite difference methods. Any sources I could read to help better my understand in order to handle the problem

```<code>% solve Ct=DCzz-vxCx

clear all

% parameters

D=9.1e-10; % Diffusion Coefficient for phenol in wastewater

vx=(0.05/10/3600)/(0.01*30*0.01)

% Domain and step

Lx=60

Lz=0.01

Nx=100

Nz=100

Nt=1000

dx=Lx/(Nx-1)

dz=Lz/(Nz-1)

% Satisfy Courant Number

C=0.05

ux=1

uz=1

dt=C/((ux/dx)+(uz/dx))

% Field Variables

Cn=ones(Nz,Nx)

x=linspace(0,Lx,Nx)

z=linspace(0,Lz,Nz)

[X,Z]=meshgrid(x,z)

% Initial conditions

Cn(:,:)=1

t=0

% loop

for n=1:Nt

Cc=Cn

for i=2:Nx-1
    for j=2:Nz-1
        Cn(j,i)=Cc(j,i)+(D*dt/(dz*dz))(Cc(j+1,i)-2*Cc(j,i)+Cc(j-1,i))-(vx*dt/dx)(Cc(j,i)-Cc(j,i-1))
    end
end

% Boundary conditions

Cn(1,:)=1

Cn(end,:)=1

Cn(:,1)=1

Cn(:, end)=Cn(:, end-1)

% Visualize

mesh(x,z,Cn); axis[0 Lx 0 Lz 0 Lz 0 1]

pause(0.01)

end
</code>```


问:
Least-squares fit of explicit parabolic sheet to data points
For a given set of data points

$$\{(x_i, y_i, z_i)\}$$

there exists some

$$f_{ABC}(x,y)=Ax^2+Bxy+Cy^2$$

that minimizes

$$\sum_i(f_{ABC}(x_i,y_i)-z_i)^2$$

$A$, $B$, and $C$ can be found quickly by directly solving/minimizing

$$(\begin{bmatrix}x_1^2&x_1y_1&y_2^2\\x_1^2&x_2y_2&y_2^2\\x_1^2&x_3y_3&y_3^2\\\vdots&\vdots&\vdots\end{bmatrix}\begin{bmatrix}A\\B\\C\end{bmatrix}-\begin{bmatrix}z_1\\z_2\\z_3\\\vdots\end{bmatrix})^2$$

However, suppose we want to reduce the model's degrees of freedom from 3 to 2 by constraining $f_{ABC}$ to be a parabolic sheet. I think this would mean that $\begin{bmatrix}A&B&C\end{bmatrix}$ has to be of the form $\begin{bmatrix}U^2&2UV&V^2\end{bmatrix}$ so that $f_{ABC}(x,y)=(Ux+Vy)^2$.

It seems this is no longer solvable directly as a linear system.

What would be a performant/fast way to solve this constrained version of the problem?

Background

The above question is a simpler version of my actual problem, which is to fit an explicit bicubic surface patch $g(x,y)=Ax^2+Bxy+Cy^2+D+Ex^3+Fy^3+Gx^2y+Hxy^2+Ix+Jy$ that is constrained to be quadratic (instead of cubic) along some (any) direction.

Thoughts

Thinking further, my best idea currently is to

Directly fit an unconstrained quadratic surface i.e. find a solution in the span of $\begin{bmatrix}x_i^2&x_iy_i&y_i^2&1&x_i&y_i\end{bmatrix}$. Call this the "quadratic solution".

Orthogonalize the remaining basis vectors $\begin{bmatrix}x_i^3&y_i^3&x_i^2y_i&x_iy_i^2\end{bmatrix}$ w.r.t. the quadratic solution.

Somehow solve a constrained minimization within this reduced basis that is orthogonal to the quadratic solution, though I don't really know how to go about this.

I think the eventual cubic terms should be of the form $Ex^3+Fy^3+Gx^2y+Hxy^2=(Ux+Vy)^3$ so the constraint is $\begin{bmatrix}E&F&G&H\end{bmatrix}=\begin{bmatrix}U^3&3U^2V&3UV^2&V^3\end{bmatrix}$

How to do this constrained optimization?


问:
fixed point iteration on DD method
I have to solve the the problem $u_t+\Delta^2u=f(u)$, where $f(u)$ is non-linear, using domain-decomposition method.

My approach is first using fixed point iteration on mixed form i.e to say $u^{k+1}_t+\Delta w^{k+1}=f(u^k)$ and $w^{k+1}=\Delta u^{k+1}$, then using usual domain decomposition procedure.

Another approach: Use fixed point iteration on the sub-problem i.e to say first use domain decomposition then for each sub-problem use fixed point iteration.

1) which approach is good?

2) Is the selection of the operator $T=$$\begin{bmatrix}(\partial_t)^{-1}+(\Delta)^{-1}\\(\Delta)^{-1}+I\end{bmatrix}$ is correct for first approach.

Thanks in advance.

Title: fixed point iteration on DD method


问:
How to construct a Fortin Operator for Crouzeix-Raviart Element?
I want to prove the LBB condition for the Stokes Equations discretised by the Crouzeix-Raviart element.

The continuous Stokes Equation in the weak formulation is

Find $u \in H_0^1(\Omega, \mathbb{R}^d)$ and $p \in L_0^2(\Omega)$ with

\begin{align} \int_\Omega \nabla u : \nabla v dx - \int_\Omega
 p \text{ div} (v) dx &= \int_{\Omega} fv dx \qquad \ \forall v \in H_0^1(\Omega, \mathbb{R}^d)\\
 \int_\Omega q \text{ div} (u) dx&= 0 \qquad \qquad \quad \forall q \in
 L_0^2(\Omega) \end{align}

It has a unique solution in the spaces $V = H_0^1(\Omega, \mathbb{R}^d)$ and $Q= L_0^2(\Omega)$.

For the discrete spaces I have

\begin{align}
V_h &= \{ v_h \in H_0^1(\Omega, \mathbb{R}^d) : v_h|_T\in P^3(T),v_h|_E\in P^2(E),\forall T \in T_h\forall E\subset \partial T\}\\
Q_h &= \{q_h\in L_0^1(\Omega):q_h|_T\in P^1(T),\forall T \in T_h \}
\end{align}

To prove solvability, I want to prove the discrete LBB condition

\begin{align}
\sup\limits_{v_h\in V_h} \frac{\int_\Omega \text{div}(v_h)q_h dx}{\|v_h\|_{V}} \geq \beta_1 \|q_h\|_{Q} \qquad \forall q_h \in Q_h
\end{align}

For the continuous spaces $V,Q$ the LBB condition holds.

Assume there exists a Fortin operator $\Pi_h:V \rightarrow V_h$ with
\begin{align}
\|\Pi_h v \|_V &\preccurlyeq \|v\|_V \qquad \quad \forall v\in V \\
b(\Pi_h v, q_h) &= b(v, q_h) \qquad \ \forall q_h \in Q_h
\end{align}
Then the continuous LBB implies the discrete LBB.

Can you please tell me how I can construct such an Fortin Operator and then prove its two properties.

Especially it would be very helpful to see the full scaling argument that proves its continuity.

I would be thankful for any help.

I have also posted this question on Math-Stack exchange, but so far without success:
https://math.stackexchange.com/questions/3710994/fortin-operator-for-crouzeix-raviart-element


问:
Comparing custom linear regression solver to SciPy equivalent in Python
From a given data set, I set out to complete a task which is below

Fit the data of the previous exercise to fit Eq. (8.18) using the SciPy function scipy.optimize.curve_fit(). Plot the data as symbols and the fit as a line on linear and on log-log axes in two separate plots in the same figure window. Compare your results to those of the previous exercise.

Here is the data I had.

I attempted the problem with the following Python script.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(K, t, p):
    return (K * np.power(t, p))

def func1(c, m, x):
    return (m * x + c)

def LineFitWt(x, y, dy):
    """Fit to straight line.
    Inputs: x, y, and dy (y-uncertainty) arrays.
    Ouputs: slope and y-intercept of best fit to data.
    """
    dy2 = dy ** 2
    norm = (1. / dy2).sum()
    xhat = (x / dy2).sum() / norm
    yhat = (y / dy2).sum() / norm
    slope = ((x - xhat) * y / dy2).sum() / ((x - xhat) * x / dy2).sum()
    yint = yhat - slope * xhat
    dy2_slope = 1. / ((x - xhat) * x / dy2).sum()
    dy2_yint = dy2_slope * (x * x / dy2).sum() / norm
    return slope, yint, np.sqrt(dy2_slope), np.sqrt(dy2_yint)

def redchisq(x, y, dy, slope, yint):
    chisq = (((y - yint - slope * x) / dy) ** 2).sum()
    return chisq / float(x.size - 2)

# Extract data from the text file for t, s and ds
t, s, ds = np.loadtxt("growthdata.txt", skiprows=3, unpack=True)

print ("time =", t)
print("size = ", s)
print ("Uncertainty =", ds)

# Set X and Y for the relevant axis
X = np.log(t)
Y = np.log(s)
dY = np.log(ds)

# Call LineFitWt function to calculate the gradient(m), y-intercept(c) and uncertainties in these(dm, dc)
m, c, dm, dc = LineFitWt(X, Y, dY)
rchisq = redchisq(X, Y, dY, m, c)
SciPy_Fit = curve_fit(func, t, s, p0=([0.54, 5.83]))

# Assign values for scipy fit
p_s = SciPy_Fit[0][0]
K_s = np.exp(SciPy_Fit[0][1])

# Calculate straight line properties for scipy parameters
m_s = p_s
c_s = np.log(K_s)

print ("m = ", m)
print ("c = ", c)
print ("dm = ",dm)
print ("dc = ", dc)
print ("Reduced chi square = ", rchisq)
print ("SciPy_Fit values  = ", SciPy_Fit)
print ("SciPy Fit p = ", p_s)
print ("SciPy Fit K = ", K_s)

# Calculate the values for p and K
p = m
K = np.exp(c)

print ("p = ", p)
print ("K = ", K)

# Calculate values for custom fit points (y = mx + c)
Xext = 0.05*(X.max()-X.min())
Xfit = (np.array([X.min()-Xext, X.max()+Xext]))
Yfit= (c+m*Xfit)

# Calculate points for log log graph using Custom fit
Y_custom = func(K, t, p)

# Calculate values for SciPy y = mx + c fit
X_scipy = X
Y_scipy = (m_s * X + c_s)

# Calculate values for log-log plot using SciPy parameters
Y_scipy_loglog = func(K_s, t, p_s)

# Assign a figure object to plot on.
plt.figure()
plt.subplot(2, 1, 1)
plt.plot(X, Y,"x", label="Data")
plt.errorbar(X, Y, yerr=dY, zorder=-1, label="Unc in s")
plt.plot(Xfit, Yfit,"+--", label="Custom Fit")
plt.plot(X_scipy, Y_scipy, "D--", label="SciPy Fit")
plt.text(-1, -2, 'Custom fit m={0:0.4f}, c={1:0.4f}'.format(m,c))
plt.text(-1, -5, 'SciPy Fit m={0:0.4f}, c={1:0.4f}'.format(m_s,c_s))
plt.xlabel("ln s")
plt.ylabel("ln r")
plt.legend()
plt.plot()

plt.subplot(2,1,2)
plt.loglog(t, s, "x", label="Data")
plt.errorbar(t, s, yerr=ds, label="Uncertainty in s", zorder=-1)
plt.loglog(t, Y_custom, "--", label="Custom fit")
plt.loglog(t, Y_scipy_loglog, "D--", label="SciPy Fit")
plt.xlabel("t")
plt.ylabel("s")
plt.text(0.1, 1, 'Custom fit K={0:0.4f}, p={1:0.4f}'.format(K, p))
plt.text(0.1, 10, 'SciPy Fit K={0:0.4f}, p={1:0.4f}'.format(K_s, p_s))
plt.legend()
plt.tight_layout()
plt.show()
```

<p>I then obtained the following graphs<img src="https://i.stack.imgur.com/CvKu0.jpg" alt="Matplotlib plots"></p>

<p>Here are the printed outputs.</p>

```time = [ 0.12  0.18  0.42  0.9   2.1   6.   18.   42.  ]
size =  [ 115.  130.  202.  335.  510.  890. 1700. 2600.]
Uncertainty = [10. 12. 14. 18. 20. 30. 40. 50.]
m =  0.5419106669494728
c =  5.837596806432137
dm =  0.5518341454711854
dc =  1.0280703524098318
Reduced chi square =  0.0002705479986781181
SciPy_Fit values  =  (array([  0.67373601, -10.71638416]), array([[-3.48433902e+11,  1.40195650e+13], [ 1.40195650e+13, -5.64090353e+14]]))
SciPy Fit p =  0.6737360069383264
SciPy Fit K =  2.217856767946089e-05
p =  0.5419106669494728
K =  342.9541642911326
```

What am I missing here? Why is the SciPy fit so far from the data and the custom fit?


问:
Why is LAPACK (seemingly) suboptimal for packed and banded eigenvalue problems?
Based on this LAPACK routines list, it looks like there is no relatively robust representation (RRR) driver routine for either packed or banded symmetric eigenvalue problems. According to the relevant LAPACK documentation though, it appears that the RRR driver is the fastest algorithm in most cases and uses the least workspace. Furthermore, based on the method used by LAPACK for symmetric eigenvalue problems, LAPACK tridiagonalizes the given matrix for all cases (full, packed, banded) and then runs an appropriate routine (RRR or otherwise) on the symmetric tridiagonal matrix thus yielded.

So this raises the question: if RRR is typically the most efficient routine, and it can (seemingly?) work for packed and banded matrices (since it works on symmetric tridiagonal matrices), why is there no RRR driver routine for packed or banded matrices?

(And if there is no good reason, would it be potentially more efficient for a user to tridiagonalize themselves using the computational routine *btrd and then call the driver routine stevr on the symmetric tridiagonal matrix thus yielded?)


问:
Compute the sum of probabilities when they are given as logits
Say I have a set of numerous probabilities given by their logarithm : $\{\ln p_i, 1 \leq i \leq N\}$.

I want to compute $\sum p_i$, if possible without exponentiating $\ln p_i$, since some of those probabilities are really small and I would suffer a dramatic loss of precision by doing so.

Do you know of any clever trick ?

Edit 14/06

I compute the probabilities along a probability tree whose depth can go s high as $D = 10 000$. This tree is typically sparse, but I don't have a better upper bound than $N = 2^D$

Some of these probabilities get very low (e.g. $10^{-200}$). I work in python. I haven't witnessed firsthand a precision loss, but I suspected it might happen and decided to reach out to more kowledgeable people.

My current implementation keeps only the 1000 highest log probabilities, exponentiates and sums them.

Computing the log of the sum is perfect for my application, since I have a natural baseline for the log probs. I'd gladly mark this as an approved answer.


问:
What is the meaning of the Helmholtz wave equation?
I am trying to build understanding on the Helmholtz wave equation $\Delta p + k^2 p = 0$, where $p$ is the deviation from ambient pressure and $k$ the wave number, in order to use it in numerical simulations. However, I'm not sure how to interpret the results. Take as an example the figure 5 on this work

More precisely, my questions are:

Since we are not in time-domain anymore, does that mean the output will be a stationary wave?

Is it possible to draw an analogy between the frequency-domain form and its time-domain form? For example, I know what happens when a propagating wave hits an obstacle and is scattered in time-domain, but in frequency-domain will we have the same behaviour?

Does the real and imaginary parts have a physical meaning?

I'd be very thankful if someone could answer at least part of my questions. Thanks a lot!


问:
A parallelized GMRES solver?
My application calls for solving a dense, 40,000 x 40,000, ill-conditioned linear system. The native SciPy GMRES solver with preconditioning has worked well for my application and solving a single system is within the ability of my local hardware. I am now looking to perform many such calculations and would like to move the calculations to cluster to improve speed. My understanding is that the native GMRES solver is default configured for use on a single core. Is there a straightforward way to incorporate the solver with Pool? Or does anyone know of an efficient parallelized GMRES solver with a Python wrapper?</p>

Title: A parallelized GMRES solver?


问:
At what l/d ratio will a frame element start to behave as a shell element?
I'm working in ETABS. There are few columns of dimension 300mmx1400mm. The height of building is 36.6 meter above ground level and the dimension of building is 26mx68m. I'm getting the time period of this model as 3.08 which is way too high. Hence i thought of changing the above mentioned columns to shell elements which will be more stiffer. Is this procedure correct?</p>


问:
scipy's solve_ivp returns erroneous results for a stiff differential equation
I'm using scipy's solve_ivp for solving a stiff differential equation. I'm using method BDF for solving the same.
I have already used MATLAB's ode23s and I'm getting correct results in MATLAB.
However, solve_ivp returns the wrong results. The problem might be in computing X.
My code is :

```<code>X0=np.array([0,0,0,0,0,0])
swd=0
start =timer()
for k in range(0,100):

    tspan=[k*(1e-4), (k+.7)*(1e-4)]
    tspan1=[(k+.7)*(1e-4), (k+1)*(1e-4)]

    swm=1
    sol=solve_ivp(conv,tspan,X0,method='BDF')
    TT=sol.t
    YY=sol.y

    swm=0
    sol1=solve_ivp(conv,tspan1,YY[:,-1],method='BDF')
    TT1=sol1.t
    YY1=sol1.y
    X0=YY1[:,-1]
    plt.plot(TT,YY[0,:],'r')
    plt.plot(TT,YY[1,:],'b')

end=timer()
print(end-start)
</code>```

for the function defined as

```<code>def conv(t,X):

    global swd
        if (swd==0) or (swd==1 and X[0]>0):

            Xdot = A11.dot(X).reshape(6,1) + B11.dot(U)
            swd=1
        else:

            Xdot = A10.dot(X).reshape(6,1) + B10.dot(U)
            swd=0
    return np.squeeze(np.asarray(Xdot))
</code>```

The matrices A's are 66* , B's are 62* resp. X is 61*. One of the X values is computed incorrectly and then the whole plot is returned incorrect. Just the plot for the first cycle is returned correctly.
Details: The given code is of a boost converter. swm and swd are the switching states of the MOSFET and diode respectively.
As far as I could debug, the diode switches correctly for the first cycle. But for the second and other cycles, the diode doesn't switch off and just remains on.


问:
Combining many probabilities, modifying, seeking general formula
CONTEXT

I need to combine the probability of occurrence of many thousands of events for millions of individuals (trees) in an agent-based/individual-based simulation model developed in NetLogo (agent-based simulation) and Go (to parallelize probability calculations). For those interested, the probabilities are for infection events that result from multiple other events, e.g., exposure, transmission, each with their own probabilities, for individual trees in a forested landscape.

Each individual can potentially be infected by thousands of other individuals in its surroundings. The full model includes millions of individuals, and so I am designing the model to calculate, modify, and update probabilities as quickly and efficiently as possible as conditions in the simulation change.

Based on the probability math that, if  and  are events that are independent and not mutually exclusive and occur with different probabilities  and , with the overall probability  calculated as:

or

In addition, each event probability  is the product of two other probabilities  and , where  is constant across events but  differs between events).

Therefore, the overall equation used in my model to calculate each individual tree's probability of infection  by any other individual (i.e., combining probability of many such infection events) is as follows:

or

or

Over time in the simulation, the value of  for a given individual tree will change (depending on changing state variables of the individual trees) while the values of  remain constant. Therefore, as the value of  changes between time steps, the value of  will need to be modified to reflect this change.

However, given the amount of values involved in the calculation (thousands of calculations for each of millions of individual trees), it is not practical to store each value of  for each individual (tree) once the (maybe not possible given the constraints of the software and computing resources I am using). In addition, I want to limit the amount of infomation that passes back and forth between NetLogo (which stores the states and probabilities for each individual) and Go (where probabilities are calculated but not stored) to avoid exceeding memory constraints?

Rather than recalculate the overall probability from scratch, which requires lots of time and computing resources, I want to be able to modify the overall probability
for a given individual tree based on the change in the value of .

QUESTION

What is a generalized form for calculating the probability that at least one event occurs given N independent events with different probabilities given the criteria above?

How can I efficiently implement the modification of each individual's probability of infection without having to re-calculate from scratch?


问:
Optimization for sampling multiple points of maximized minimum distance
I'm trying to find a way to sample new points that have maximum minimum-distance (maximin distance).
The current situation is where there are ns number of pre-existing sample points. I want N number of new sample points that are filling the space as much as possible (hence, the maximization of minimum distance).

The most conventional way would be to sample single points in a loop (performing N loops).
However, what I want is to find multiple points that have distances maximized for each other as well as from the existing points.
Currently, I'm using a population-based global optimizer where a variable is set as follows:

When sampling N points with nx number of variables, a single population would be a vector with a length of (N x nx). The objective function to be maximized would be the total sum of minimum distances between the N points and also the distance between N and ns (existing) points.

The problem with this method is that not only the problem is highly non-linear but also the size of the number of variables for the optimization problem (N x nx) increases rapidly as the number of variables for the sample space (in other words, dimension) nx increases.

I was wondering is there any other way than to vectorize multiple samples?
And also, is there a global optimization algorithm that gives multiple individuals from a population as its optimum result (not a multi-objective optimization as the problem in my case has a single objective: distance)?


问:
Non-negative Least Squares to perform Inverse Laplace with weights
I'm trying to perform the inverse Laplace transform of a (noisy) dataset $y_i$ using Tikhonov regularization:

$$\min \sum_{i=1}^{N} \left(\int_0^\infty e^{-s_i t} f(t) \, dt - y_i \right)^2 - \lambda^2 \int_0^\infty f(t)^2 dt$$

where I know that $f(t) \geq 0$ for $t\geq 0$. I therefore convert that to a Non-negative least squares problem following the answer here with

$$\begin{aligned}
A_{ij} &= e^{-s_i t_j} \sqrt{\Delta t} \\
b_i&=y_i
\end{aligned}
$$
where $\mathbf{A}$ is an $N \times L$ matrix.

Analytically we can calculate that
$$ M(s_i, s_j) = \int_0^\infty e^{-s_i t}e^{-s_j t} dt = \frac{1}{s_i+s_j} \iff \mathbf{M} = \mathbf{A} \mathbf{A}^T$$
and numerically I get this up to a precision of $10^{-3}$ if I choose my $t$-axis with $L=1000$ points evenly spaced in log-space between $10^{-10}$ and $10^{10}$ (which is much better than with linear spacing).

My problem is this: I would like to include a diagonal weights matrix $W_{ii} = 1/\sigma_{y_i}^2$ as

$$\min \sum_{ij} \left(\int_0^\infty e^{-s_i t} f(t) \, dt - y_i \right) W_{ij} \left(\int_0^\infty e^{-s_j t} f(t) \, dt - y_j \right) - \lambda^2 \int_0^\infty f(t)^2 dt$$

and thought that I could just update $\mathbf{A}$ and $\mathbf{b}$ to
$$\begin{aligned}
\mathbf{A}' \to \mathbf{W}^{1/2} \mathbf{A} \\
\mathbf{b}' \to \mathbf{W}^{1/2} \mathbf{b}
\end{aligned}
$$
as described e.g. here.

However, if I now calculate $W_{ij} M(s_i, s_j)\longleftrightarrow \mathbf{A}' (\mathbf{A}')^T$ analytically and numerically I get different answers. This is the upper 4x4 part for analytical and numerical respectively:

```<code>array([[2.82777399e-06, 3.92720333e-06, 4.59359056e-06, 5.08745371e-06],
   [3.92720333e-06, 6.13584813e-06, 7.65547638e-06, 8.83179899e-06],
   [4.59359056e-06, 7.65547638e-06, 9.94943930e-06, 1.18061968e-05],
   [5.08745371e-06, 8.83179899e-06, 1.18061968e-05, 1.43013249e-05]])
</code>```

and

```<code> array([[2.82752363e-06, 1.88501575e-06, 1.41376182e-06, 1.13100945e-06],
   [8.18040652e-06, 6.13530489e-06, 4.90824391e-06, 4.09020326e-06],
   [1.49228376e-05, 1.19382701e-05, 9.94855841e-06, 8.52733578e-06],
   [2.28800939e-05, 1.90667449e-05, 1.63429242e-05, 1.43000587e-05]])
</code>```

So only the diagonal is the same, but the off-diagonal terms are wrong and what's worse, the result is no longer symmetric. So what am I doing wrong here? Any help is much appreciated!


问:
Numerical scheme to calculate the normal mode of a set of hyperbolic PDEs?
I would like to solve the linearised, ideal, MHD equations, where the gas pressure is zero.
$$\frac{\partial u_x}{\partial t}=v_A^2(x,z)\left[\nabla_{||}b_x - \frac{\partial b_{||}}{\partial x}\right],$$
$$\frac{\partial u_\perp}{\partial t}=v_A^2(x,z)\left[\nabla_{||}b_\perp - \nabla_\perp b_{||} \right],$$
$$\frac{\partial b_x}{\partial t}=\nabla_{||}u_x,$$
$$\frac{\partial b_\perp}{\partial t}=\nabla_{||}u_\perp,$$
$$\frac{\partial b_{||}}{\partial t}=-\left[\frac{\partial u_x}{\partial x}+\nabla_\perp u_\perp\right],$$
where
$$\nabla_\perp = \cos\alpha\frac{\partial}{\partial y} - \sin\alpha\frac{\partial }{\partial z}$$
$$\nabla_{||} = \sin\alpha\frac{\partial}{\partial y} + \cos\alpha\frac{\partial }{\partial z}$$
and $v_A(x,z)$ denotes the background Alfven speed.
The magnetic field perturbations $b_x$, $b_\perp$, $b_{||}$ are dimensionless because they have been normalised by $B_0$ and the background field is given by
$$\boldsymbol{B}_0=B_0(\sin\alpha\,\boldsymbol{\hat{y}}+\cos\alpha\,\boldsymbol{\hat{z}}).$$

I would like to assume a $y$ and $t$ dependence of the form $\exp[i(k_y y + \omega t)]$ as this means we don't need to integrate in $y$ and $t$ which saves time.

Do you know a numerical scheme I could use to solve these equations? The boundary conditions I had in mind are periodic in $z$ and imposing a driver at $x=x_{min}$ and $x=x_{max}$, where $v_A(x,z)$, $k_y$ and $\omega$ are prescribed beforehand.


问:
Origin of phrase `computational microscope'
I have heard the term 'computational microscope' used to describe the practice of molecular simulation (in the context e.g. computational chemistry, materials science) and its use as a numerical tool for predicting scientific properties of materials, etc. I have also heard 'numerical microscope' in the same context.

Is there a well-known origin for this terminology / expression? References would be ideal.


问:
Effect of reducing flux consistency order at boundary on convergence order
Consider the 1D nonstationary convection-diffusion PDE
$$
\begin{alignat}{2}
\partial_t u &= -a \partial_x u + D \partial_{xx}u, &\qquad x \in (0,1), t \in (0,T), \\
f(t) &= \left.\left( a u - D\partial_x u \right) \right|_{x=0} & \qquad t \in (0,T), \\
0 &= \left. -D\partial_x u \right|_{x=1} & \qquad t \in (0,T), \\
u_0 &= \left. u \right|_{t=0} & \qquad x \in (0,1),
\end{alignat}
$$
where $a, D > 0$ are constant and $f,u_0$ are given functions. The problem is convection-dominated (i.e., $a \gg D$).

A cell-centered finite volume method with uniform grid is used to discretize the PDE in space. In the interior, central differences for the diffusive flux and a second-order accurate convective flux (e.g., QUICK) can be used. Crank-Nicolson is used in time to obtain an overall second-order convergent scheme.

The boundary conditions eliminate the diffusive fluxes in the discretization at the boundaries.
For the right faces of the first and last cell, however, there are not enough neighbouring cells to apply the second-order convective flux (i.e., the QUICK scheme requires the left and right cells of a face and one additional upstream cell).

Now, for these two faces, an upwind flux is used, which is only first-order accurate. Theoretically, the whole method is, hence, reduced to first-order consistency.

In computational experiments, second-order convergence is observed, though. This suggests the reasoning that an $O(h)$ error committed in a region of size $O(h)$ does itself result in an $O(h^2)$ error, which does not impede overall $O(h^2)$ convergence.

Is this reasoning correct? How can it be proved / are there any references?

According to (Randall LeVeque,  Finite Difference Methods for Ordinary and Partial Differential Equations , Sec. 2.12), this does not hold for Poisson's equation with central and one-sided finite differences. So it may not hold in general (i.e., general equation type, general discretization method).


问:
Solving Stokes Equations in 3D - Do I need to treat pressure-velocity coupling iteratively?
I need to develop a code to solve Stokes Equations in 3D in cubic geometries (structured grid, uniform mesh spacing).

My code needs to take a pressure gradient in one direction as a BC (pinlet=p1, poutlet=p0) and no-slip BC's in the other 4 faces of the cubic geometry.

I started to develop this code in Python, using a staggered grid FV approach; pressure at the center of the FV cells, ux velocities at the faces normal to the x direction and so on.

My code gives strange results, which has made me wonder if I am considering the correct approach to solve such system.

I am currently assembling a system of equations that considers [p, ux, uy] as unknowns and solving it in "one go". That is, I build a matrix and solve the linear system, which gives me the solution. The degrees of freedom are coupled through Stokes equations and continuity equation, but the coupling is "solved" by assembling the coefficients to the matrix.

However, whenever I research numerical codes and books on Navier-Stokes equations, they usually give a couple of strategies to treat the pressure-velocity coupling, typically "guess and correct" approaches, where you guess a pressure field, and correct it iteratively to satisfy continuity (e.g. SIMPLE, SIMPLEC, PISO). These approaches are usually formulated in the Navier-Stokes context, but Stokes equations have the same type of pressure-velocity coupling; they just don't have the nonlinear inertia term. So I assume these strategies would be applicable to Stokes equations as well.

So my question to the community is: Is there something "wrong" with my approach to solve Stokes equations in one go, implicitly, by solving one linear system? Or do I need to consider iterative methods to treat the pressure-velocity coupling?

Thank you!
Rafael March.


问:
Entropy of a spatially and temporally varying fluid system?
I'm trying to analyze the change in system thermodynamic entropy of an ideal gas system.

For reference, I'm analyzing solution techniques for solving the 1D compressible Euler equations:
$$
\partial_t \rho + \partial_x (\rho u_x) = 0\\
\partial_t (\rho u_x) + \partial_x (\rho u_x^2 + P) = 0\\
\partial_t e + \partial_x ((e + P) u_x) = 0\\
e = \frac{P}{\gamma - 1} + \frac{1}{2} \rho u_x^2\\
\partial_t (\rho s) + \partial_x  (\rho s u_x) = \dot{S}_{gen}
$$

The first test problem I want to analyze is the classical shock tube configuration where the jump discontinuity is located at $x=0$:
$$
\rho_L = 1\\
P_L = 1\\
\rho_R = 4\\
P_R = 4\\
u_x(t=0) = 0\\
\gamma = \frac{5}{3}
$$

To get an idea of what I should expect, I computed the analytical solution for this problem (using Toro's Riemann Solvers book), and computed the entropy as
$$
\frac{\rho s - S_0}{c_v} = \ln\left(\frac{P}{P_0} \left(\frac{\rho_0}{\rho}\right)^\gamma\right)
$$
where my understanding is that $\rho_0$, $P_0$, and $S_0$ are user-chosen reference values. I decided to chose $\rho_{0} = 4$, $P_0 = 4$, and $\frac{S_0}{c_v} = 1$, though this doesn't necessarily effect the final result. Here's the entropy density at two times $t=0$ and $t=1$:

As a sanity check I defined the domain $x \in [-4, 4]$ and tried to integrate $\rho s$ over the domain and compared this to the initial system entropy:

$$
\left.\int^4_{-4} \rho s dx\right|_{t=1} - \left.\int^4_{-4} \rho s dx\right|_{t=0} = -0.477216
$$

This would be acceptable if there was mass flow out of the domain boundaries, but I chose the domain such that the shock and rarefaction fronts don't reach the domain boundaries. This seems to imply that
$$
\int^1_0 \int^4_{-4} \dot{S}_{gen} dx dt = -0.477216
$$
which is only possible if $\dot{S}_{gen} < 0$ somewhere in the domain, a clear violation of the 2nd law of Thermodynamics. However, my understanding is that the Euler equations should obey the 2nd law.

My question is how do I resolve this issue? How should you use the 2nd law to analyze the entropy flow/generation of a spatially and temporally varying fluid field (both analytically, as well as for CFD simulations)?

Here's a plot of $u_x$ at $t=1$ for reference:


问:
Is the matrix exponential and the Jordan canonical form actually useful for solving differential equations?
All of my yearlong graduate-level Linear Algebra course notes from my professor—an algebraist/representation theorist—shows his love for the exponential map $e^A$ and the Jordan canonical form—and one main reason given by him is that they help solve linear differential equations.

But how true is this in practice?

I have little scientific computing experience.

The thing is, in my research interest, fluid dynamics, not only are the evolution equations nonlinear, often with nonconstant coefficients. Do the matrix exponential and the Jordan form help with modeling and simulation in fluid dynamics problems?


问:
Numerical calculation of the Berry connection
I'm doing some numerical calculations involving Hermitian matrices, and derivatives of the eigenvectors.
Essentially, I have an n x n, Hermitian matrix H(x), which is dependent on some continuous parameter (x). I then need to calculate the Berry connection for this system. So suppose $U^\dagger(x) H(x) U(x) = \epsilon_n \delta_{nm}$, I'm essentially computing $i U^\dagger(x) \partial_x U(x) = \hat{A}$.

To do this, I tried to evaluate $H(x), H(x+\Delta x)$, and obtained $U(x),U(x+\Delta x)$. I therefore approximated $\hat{A} \approx i U^\dagger(x) (U(x+\Delta x )- U(x))/(\Delta x)$. This is carried out in LAPACK, using HEEVD (CHEEVD), since H is Hermitian, and parametrized with double precision.

However, since the diagonalization at every x carries with it a random phase (apparently), the object $\hat{A}$ contains errors, and no longer represents just the derivative of the rotation matrix $U$. In fact, it even turns out to be non-Hermitian at some x points.

I imagine the solution has to be some sort of uniform gauge condition for all $x$, ensuring the smoothness of the connection. Is there anyway to add this, or implement this within LAPACK? I have so far been unsuccessful.

Thanks!

P.S.
Naturally, the Berry connection is not gauge invariant, so it won't directly enter any physical observable. The issue is, however, that errors in the evaluation of $\hat{A}$ propagate directly to observable quantities, like the berry curvature for example. Quantities $i[A,A]$ suddenly contain disallowed imaginary parts, coming in from errors in $\hat{A}$.


问:
Plotting a Magnetic Field in Spherical Coordinates in Python
I am modeling a Helmholtz Coil as two dipoles from far away and I want to plot the magnetic field.

$$\mathbf{B}(\mathbf{r}) = \frac{\mu_0 |\mathbf{m}|}{4\pi r^3}\left(2\cos\theta\,\hat{\mathbf{r}} + \sin\theta\,\hat{\boldsymbol{\theta}}\right)\, .$$

I am struggling to plot the magnetic field in spherical coordinates in Python. Can somebody show me an easy way to plot 3D functions in spherical coordinates?


问:
Compute Nullspace of Sparse Matrix
I am computing the nullspace of a sparse rectangular $m$ x $n$ matrix $A$, where $m$ << $n$. I do this by computing the QR decomposition of $A^T$ and extract the $n-m$ right-most columns of the resulting $Q$. Symbolically:
$$A^T=QR, Q=[Q_1 Q_2] , Q_2 = Nullspace(A)$$

Using this technique, I get the results I want but I would like it to be faster.

Here is the sparsity pattern of $A^T$ as generated by Mathematica (2,337 x 36,090; 42,900 nonzeros):

and the resulting nullspace (36,090 x 33,753; 8,167,054 nonzeros)

The sparsity pattern arises out of the constraints on the contact of two subdomains in a finite element problem. The coefficients on the left are positive whereas the coefficients on the right are negative. This is due to a sign convention. I ordered the node indices of each subdomain using a nodal connectivity matrix from the finite element mesh and fed that into METIS (which made my conjugate gradient solves MUCH faster!).

I primarily use Eigen and Intel MKL for my matrix algebra, but also use SuiteSparse for their much-faster QR decomposition routine. Although not important for this question, I am including my procedure for interfacing between Eigen and SuiteSparse, just in case it does matter. Here is how I compute the nullspace (kernel) of $A$:

```<code>using namespace Eigen;
typedef SparseMatrix<double,ColMajor,int> Sparse_t;

cholmod_common* cc, Common;
cc = &Common;
cholmod_l_start(cc);

Sparse_t AT; //assume A^T is already defined
//unfortunately necessary copy with SuiteSparse_long indices -- minimal impact on performance
SparseMatrix<double, ColMajor, SuiteSparse_long> tmp(AT);

const int rows = AT.rows();
const int cols = AT.cols();
const int ordering = SPQR_ORDERING_BEST;
const double tol = SPQR_NO_TOL;
const SuiteSparse_long econ = rows;
const int getCTX = 1;

cholmod_sparse AT_cm = viewAsCholmod(tmp);
cholmod_sparse* I_cm = cholmod_l_speye(rows, rows, CHOLMOD_REAL, cc);
cholmod_sparse* Q_cm = NULL;

//perform QR factorization and request (Q^T * I)^T = Q
SuiteSparseQR<double>(ordering, tol, econ, getCTX, &AT_cm, I_cm, NULL, &Q_cm, NULL, NULL, NULL, NULL, NULL, NULL, cc);
Sparse_t kernel = viewAsEigen<Sparse_t::Scalar, Sparse_t::Options, Sparse_t::Index>(*Q_cm).rightCols(rows - cols);
cholmod_l_free_sparse(&I_cm, &cc);
cholmod_l_free_sparse(&Q_cm, &cc);
</code>```

From my research on this website and other sources, the QR decomposition seems to be the best way to compute the nullspace of a sparse matrix. Furthermore, SuiteSparse's SPQR appears to be the fastest and most versatile routine for sparse matrices. I have compiled everything on Windows and linked SuiteSparse against Intel MKL (probably not relevant in this problem) and Intel TBB, and I do achieve limited multithreading while using SPQR.

I think my code is optimized for my platform of choice (Windows) on a six-core CPU laptop. I didn't have much luck getting the GPU accelerated features of SPQR to work, but that could be my weak GPU. However, my real question is if there is a better way of approaching this problem. Is there a linear algebra trick I am missing? Maybe some matrix decomposition I can perform on $A^T$ before performing the QR?

Thanks! (I wonder how long it will take someone to delete the pleasantries)

Update: I used umfpack to calculate the LU decomposition first, and then calculate the QR decomposition. $$P_rRAP_c=LU , U^T=QR , Q = [Q_1Q_2] , Nullspace(A) = P_cQ_2$$ I have not posted the additional code because it is a confusing c-style mess! The nullspace is computed in a fraction of the time (including the LU decomposition), but the resulting sparsity pattern is quite bizarre which appears to slow down subsequent operations. My numerical results appear to be unaffected. If $A$ is very sparse, there is negligible improvement in my overall finite element calculations. As $A$ becomes less sparse (more constraints), there is a very significant increase in performance. Maybe a different node ordering strategy would fix this.

Compared to the QR decomposition, the LU decomposition is instantaneous. Is there any way to compute the nullspace directly from a decomposition in the form of $P_rRAP_c=LU$?

Title: Compute Nullspace of Sparse Matrix


问:
Why is my numerical solution to a set of ODEs infinite?
I am trying to solve the following linear PDEs
$$\frac{\partial u_x}{\partial x}=-[i\omega b_{||}+\nabla_\perp u_\perp],$$
$$\frac{\partial b_{||}}{\partial x}=-\frac{i}{\omega}\mathcal{L}u_x,$$
$$\mathcal{L}u_\perp = i\omega\nabla_\perp b_{||},$$
for $u_x=\texttt{ux}$, $b_{||}=\texttt{b_par}$ and $u_\perp=\texttt{u_perp}$ which are functions of $x$ and $z$.
Where
$$\nabla_\perp = ik_\perp - \sin\alpha\frac{\partial }{\partial z},$$
$$\mathcal{L} = \cos^2\alpha\frac{\partial^2}{\partial z^2}+\frac{\omega^2}{v_A^2(x,z)}.$$

To solve this I am trying to use the method of lines. I discretize the $z$-derivatives which gives set of ODEs. The code is for the case where $v_A(x,z)=v_{A0}$. The attached graph shows that the solutions go to infinity. Does this mean the scheme is numerically unstable? Note that if I change the line of code which reads

```<code>    dux_dx = -(1j * omega * b_par + 1j * k_perp * u_perp - \
                np.sin(alpha) / (2 * dz) * (np.roll(u_perp,-1) - np.roll(u_perp,1)))
</code>```

to

```<code>    dux_dx = -(1j * omega * b_par_exact(x,z) + nabla_perp0 * u_perp_exact(x,z))
</code>```

then the code works perfectly. This suggests the matrix inversion is causing the solution to go to infinity. Do you know how to fix this problem?

Code:

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

def ux_exact(x,z):
    return ux0 * np.exp(1j * (kx * x + kz * z))

def b_par_exact(x,z):   # Note that b_par denotes b_{||}
    return b_par0 * np.exp(1j * (kx * x + kz * z))

def u_perp_exact(x,z): # Note that perp is short for perpendicular
    return u_perp0 * np.exp(1j * (kx * x + kz * z))

def calc_u_perp(x, b_par):
    # Calculate u_perp by inverting the matrix ab,
    # Note that A is nearly tridiagonal and so a computationally efficent algorithm can be used
    # See https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)
    [ab, u, v] = calc_A_u_v(x)
    d = vector_d(b_par)

    y = solve_banded((1,1), ab, d)
    q = solve_banded((1,1), ab, u)

    vTy = v[0] * y[0] + v[-1] * y[-1]
    vTq = v[0] * q[0] + v[-1] * q[-1]
    u_perp = y - vTy / (1 + vTq) * q

    return u_perp

def calc_A_u_v(x):
    a = np.full(nz, np.cos(alpha) ** 2 / dz ** 2, dtype=complex) # Off diagonal elements
    b = -2 * a + omega ** 2 / vA0 ** 2                           # Diagonal elements
    b[0]  += b[1]
    b[-1] += a[0] ** 2 / b[1]

    ab = np.zeros((3, nz), dtype=complex)
    ab[0,:] = a
    ab[1,:] = b
    ab[2,:] = a

    u = np.zeros(nz, dtype=complex)
    v = np.zeros(nz, dtype=complex)
    u[0] = -b[1]
    u[-1] = a[0]
    v[0] = 1
    v[-1] = -a[0] / b[1]

    return [ab, u, v]

def vector_d(b_par):
    return -omega * (k_perp * b_par + \
            1j * np.sin(alpha) / (2 * dz) * (np.roll(b_par,-1) - np.roll(b_par,1)))

def dUdx(x, U):
    ux = U[0:nz]
    b_par = U[nz:]
    u_perp = calc_u_perp(x,b_par)
    # dux_dx = -(1j * omega * b_par_exact(x,z) + nabla_perp0 * u_perp_exact(x,z))
    dux_dx = -(1j * omega * b_par + 1j * k_perp * u_perp - \
                np.sin(alpha) / (2 * dz) * (np.roll(u_perp,-1) - np.roll(u_perp,1)))
    db_par_dx = -1j / omega * ( \
                np.cos(alpha) ** 2 / dz ** 2 * (np.roll(ux,-1) - 2 * ux  + np.roll(ux,1))  + \
                omega ** 2 / vA0 ** 2 * ux)
    return np.concatenate((dux_dx, db_par_dx))

Lz = 1
kx = 1
kz = np.pi / Lz

alpha = 0.25 * np.pi
k_perp = 0.5 * np.pi
vA0 = 1
omega = vA0 * np.sqrt(kx ** 2 + k_perp ** 2 + kz ** 2 - 2 * kz * k_perp * np.sin(alpha) + 0j)

nx = 256
x_min = 0
x_max = 6
dx = (x_max - x_min) / (nx - 1)
x = np.linspace(x_min, x_max, nx)

nz = 128
z_min = 0
z_max =  2 * Lz
dz = (z_max - z_min) / nz
z = np.linspace(z_min + dz / 2, z_max - dz / 2, nz)

X, Z = np.meshgrid(x, z)

nabla_perp0 = 1j * (k_perp - kz * np.sin(alpha))
nabla_par0  = 1j * kz * np.cos(alpha)
L0          = nabla_par0 ** 2 + omega ** 2 / vA0 ** 2

ux0     = 1
b_par0  = -L0 / (omega * kx) * ux0
u_perp0 = -1j * kx * nabla_perp0 / (L0 + nabla_perp0 ** 2) * ux0

ux_x_min    = ux_exact(x_min, z)
b_par_x_min = b_par_exact(x_min, z)
U_x_min = np.concatenate((ux_x_min, b_par_x_min))
sol = solve_ivp(dUdx, [x_min, x_max], U_x_min, t_eval=x)

ux = sol.y[0:nz]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Z, np.real(ux))
ax.set_xlabel('x')
ax.set_ylabel('z')
ax.set_title(r'Re(<span class="math-container">$u_x$</span>)')
plt.savefig('Figures/real_part_ux.png')
</code>```

Ouputted figure:


问:
Some formulations of domain coupling lead to saddle point problems. Is this merely an artifact?
Background I wanted to learn how to couple FEM and BEM (for the Poisson equation), because I wanted to better understand how open boundary conditions look like. Therefore I worked through the relevant sections of “Principles of Boundary Element Methods” by Martin Costabe, summer course lecture, 1986 (28 pages). The last section provides a formulation of the coupling, which leads to a symmetric saddle point problem. Of course, I would prefer a symmetric positive definite problem instead.

Even without BEM, some formulations of domain coupling lead to saddle point problems. Consider two domains $\Omega_1$ and $\Omega_2$ with common boundary $\Gamma$, and setup Poisson equations with homogeneous Neumann boundary conditions on $\Gamma$ in both domains. To couple the two domains, require that the potentials agree on the boundary, and add the transpose of this restriction by a Lagrange multiplier (to account for the fact that the normal derivative on both sides of $\Gamma$ is no longer restricted to 0, but must only agree with the normal derivative on the other side). The block structure of the matrix of the resulting linear system is

$\pmatrix{A_1 & 0 & B_1 \\
0 & A_2 & B_2 \\
C_1 & C_2 & 0}$

Here $C_1 = B_1^T$ and $C_2 = B_2^T$.

In this case, $C_1$ and $C_2$ have a simple structure, like $C_1=\pmatrix{0 & I}$ and $C_2 = \pmatrix{-I & 0}$ if we order the nodes suitably. Therefore, it is often possible to eliminate the constraint and the Lagrange multiplier explicitly. My confusion is that resulting system seems to be easier to solve than directly trying to solve the initial saddle point problem.

What is different when I try to couple FEM and BEM is that the structure feels more complicated to me, and I know longer know how to explicitly eliminate the Lagrange multiplier. So now I try to understand why certain approaches to couple domains lead to saddle point problems in the first place. One of my guesses is that the Lagrange multiplier which occur in this formulation effectively correspond to the normal derivative of the potential on the inner interface. But why should a formulation where both the potential and (directional) derivatives of the potential occur lead to saddle point problems? (And is there a simple why to avoid saddle point problems, or at least a simple way to solve the resulting saddle point problems without sacrificing performance of iterative solvers like CG?)


问:
Cell-centered DG extension to the two-point flux approximation scheme
A current problem that I am working on requires me to compute the solution from the heat diffusion evolution on a discontinuous function. More precisely - I have a Delaunay triangulation and within each triangle a polynomial function is defined such that it can be discontinuous at the edges. For a constant function within each triangle the two-point flux approximation scheme is applicable (precisely because I have a Delaunay triangulation).

Is there a DG extension/analogue to such a cell-centered scheme? I would like to be able to use a similar cell-centered scheme on cells with higher polynomial degrees.

To put it simply, I am looking for a (simpler) method to perform homogeneous diffusion on a mesh where discontinuous polynomial functions are defined within each cell (preferably one that extends to higher degree polynomials unlike the FVM scheme that I mentioned).


问:
How property-invariance is imposed to neural nets?
I was wondering how specific symmetries or constraints such as property-invariance transformation are imposed on any (deep) neural net when they are trained.

I'll appreciate it if anyone can aware me of the relevant researches in this direction.


问:
Solver for large dense BVP system in python
I have a large system of boundary value problems of the form
$$ \frac{d^2 y }{dt^2} = C(t) y + b(t), $$
where the variable $y$ is a vector that has anywhere from 50 to around 500 components, $C$ is a dense matrix whose entries come from data from a simulation, so there is no closed form for the entries, and $b$ also comes from simulation data.

My boundary conditions are $$y'(0)=y'(T)=0,$$ and I have on the order of 100 or 1000 time steps.

Is there a good python package/environment that could be used to solve this efficiently (or at all)? I would also be willing to use a solver in a different language if necessary.

I have tried to use the solver in Scipy, which works well when the number of components of $y$ is less than 100 but has memory problems if there are more. I have also tried to use Fenics to solve this, but Fenics cannot create a tensor space large enough to deal with $C$.


问:
Is Romberg integration method implemented as weighted function values numerically correct?
I have to integrate expression f(x) * g(x) for many different functions f but just one g.

I want to integrate it as sum of weighted values of f(x) * g(x) instead of calculating the table.  Note that in Python I may write:

```<code>sum(w[i] * f(x[i]) * g(x[i]) for i in range(2 ** k + 1))
</code>```

as:

```<code>wg = [w[i] * g(x[i]) for i in range(2 ** k + 1)]
sum(wg[i] * f(x[i]) for i in range(2 ** k + 1))
</code>```

where w[i] are weights of function values used by the Romberg method which may be calculated like:

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

w = romb(np.eye(2 ** k + 1))
</code>```

Is such implementation of Romberg method safe?

(The question has been also asked at SO: https://stackoverflow.com/questions/62659798/is-romberg-integration-method-implemented-as-weighted-function-values-numericall)


问:
Euler Explict Method for solving the PDE for option prices under the Schwartz mean reverting model. Numerical Finance
I have to solve the following PDE for a Call option :

$\partial_tV + \{ \alpha - (\mu - \lambda/ \alpha -log(S))\}S\partial_SV + 1/2 \sigma^{2}S^{2}\partial_{S}^{2}V - rV = 0$

Where V(S,t) is the price of the option as a function of stock price S and time t, $\alpha$ the strength of mean reversion, $\mu$ mean rate of return, $r$ is the risk-free interest rate, and $\sigma$  is the volatility of the stock, $\lambda$ Market price of Risk.
I have solved the PDE with finite difference for S and, as required for the exercise, Euler Explict for the time (I upload the Matlab code here with all the data).
My question is about the boundary condictions. For the Black-Scholes equation is easy:

Left end S = 0: So payoff(S) = payoff(0) = max(0-K,0) = 0, hence $V(0,t) = 0 \forall t\ge 0$.

Right end S = M: Where M is 'large' with respect to K, payoff(S) = max(S-K,0) = S - K, so the value of an option with payoff equal to S - K is just $V(S,t) = S - Ke^{-rT}$, $\forall t\ge 0$.

The left end condiction remains the same, but i don't understand how to change the right end.

```<code>% Generate a Path for the Schwartz single factor model.
% Call
% We start from : dS = ...., then X = log(S).
clear all
close all

%% Parameters
T = 1;           % End time
K = 1.5;         % Strike price
alpha = 10;     % strength of mean reversion
mu = 1;         % long term mean
sigma = 0.7;    % volatility
r = 0.1;        % risk-free interest rate
lambda = mu - r;% Market Price of Risk
lambdaHat = lambda/alpha;
Nt = 5000;      % number of time steps

%% Eu for PDE
ds = 0.01;
dt = T/Nt;
t = 0:dt:T;
Ns = round(3*(K/ds));   % Number of "price" step
S = ds*(0:(Ns));
S1 = S(2:(end-1));       % all but endpoints
V = max(S-K,0);
hold on
plot(S,V,'o-b')

%% PDE solution
Ws = alpha*(mu - lambdaHat - log(S1)).*S1;
Zs = 0.5*sigma^2*S1.^2;
Vnew = zeros(Ns+1,1);
Vold = V;
for k = 2:Nt
    F = zeros(Ns+1,1)';
    dV = (Vold(3:end) - Vold(1:(end-2)))/(2*ds);
    ddV = (Vold(3:end) - 2*Vold(2:(end-1)) + Vold(1:(end-2)))/(ds^2);
    F(2:(end-1)) = Ws.*dV + Zs.*ddV - r*Vold(2:(end-1));
    Vnew = Vold + dt*F;
    %% Boundary cond
    Vnew(1) = 0; % left end condition
    Vnew(end) = S(end) - K*exp(-r*k*dt); % right end condition
    % Update solution
    Vold = Vnew;
end
hold on
plot(S,Vnew,'-r')
</code>```


问:
Surface mesh from labeled 3D points
I'm trying to figure out how to create a surface mesh from a set of labeled 3D points. The 3D object could be something like part of a cave system or asteroid where there would be parts of the surface not visible from looking on the outside. Each point has a radius and not all points are on the surface. The output should be a Nx3 matrix (coordinate matrix) and another matrix to represent the triangles of the mesh (as in the Libigl tutorial https://libigl.github.io/tutorial/). What algorithms are there to find the points on the surface and what algorithms can create meshes from 3D points with complicated pockets? A example 2D slice is below.

Title: Surface mesh from labeled 3D points


问:
Adaptive quadrature methods for Fourier Optics
In Fourier Optics one often needs to compute approximations to bivariate integrals like
$$
\int_{-\frac{l}{2}}^{\frac{l}{2}}\int_{-\frac{l}{2}}^{\frac{l}{2}} {\rm e}^{i\phi(\xi,\eta)}\mathrm{exp}\left[i\frac{\pi}{\lambda d}\left((x-\xi)^2-(y-\eta)^2\right)\right]{\rm d}\xi{\rm d}\eta
$$
These integrals generally reduce to a form that, after discretization, can be solved by a 2-D DFT (see [1, Chap. 5]). There are some technical issues however that have to be taken care of, like aliasing, etc. I was wondering if there are any adaptive quadrature algorithms that would work more like a "black box"?

[1] Goodman, J. W. (2017). Introduction to Fourier Optics.


问:
Computation of a functional for large values
Consider the following function :

$$f(x) = \sin^2(\frac{π\Gamma(x)}{2x})$$

Now consider the following functional :

$$I(x)=\int_0^\infty \frac{f(x + iy) − f(x − iy)}{e^{2πy}-1} dy$$

I need values for of $I(x)$ for large $x$ at least upto 10^2.

I tried it on Mathematica but I couldn't .

Here are some calculated values for small X's :img1 img2

Also if possible , what is the apparent asymptotic of minimas and Maximas of the functional?


问:
Benchmark instances for directed 3-Cycle cover
The directed 3-Cycle cover asks for a vertex-covering set of oriented cycles with at least three vertices per cycle such that every vertex is covered by exactly one cycle.

I have scrutinzed the internet for some time for freely availabe benchmark problems but couldn't find anything useful.

So my question is where I can find such instances with known optimal solutions that were hard to find, resp., how such instances can be created.

The reason for asking is that I found a new algorithm for the problem and want to get a feeling about its quality and robustness.


问:
RK4-method starts oscillating above certain input parameters
I am trying to solve an equation of the following type
$$\partial_zE(z)=-c_0J$$
with $$J=c_1\beta E^3(z)$$
using the boost::odeint-framework and a fixed time stepper, with $c_0$, $c_1$ and $\beta$ constant values.

The right hand side of the function is calculated by the following procedure (depending on the step size $\Delta z$):

Multiplication of $E$ with a vector $\exp(i\cdot k_z\cdot\Delta z)$ to move the field from $z=0$ to $z=\Delta z$

Calculation of $J$, with the equation given above

Multiplication of $J$ with a vector $\exp(-i\cdot k_z\cdot\Delta z)$ to move the result back from $z=\Delta z$ to $z=0$

Return the value to the solver

This means that for a common RK4-approach I have to call that function four times, once for $z=0$, twice for $z=0.5\Delta z$ and once for $z=\Delta z$. This approach works fine for $\max(E)\ll E_{lim}$, but for $\max(E)\approx E_{lim}$ I have to reduce my step size $\Delta z$ significantly, else my resulting values for $E$ will explode within a few steps. $E_{lim}$ can be estimated by using $J_{lim}\approx 10^{13}$ (which might depend on the values for $c_0$ and $c_1$). If $J$ stays below that threshold, even large steps are not an issue, but the closer I get to that limit, the smaller I have to choose my steps (which in turn increase my calculation time).

Why does my equation behave like that? When calculating the solution for a similar type equation
$$f'(t)=-h\cdot f^3$$
on WolframAlpha I get
$$f(t)=\pm\frac{1}{\sqrt{c_0+2ht}}$$

To address the comment from @MaximUmansky:

I am trying to solve a reduced version of the UPPE, written as
$$\partial_z\hat{E}(\omega, z)=ik_z\hat{E}(\omega, z)-\frac{\omega}{2\varepsilon_0 c^2k_z}\hat{J}(\omega, z)$$
with $J$ given by the non-linear $K$-photon-absorption ($K=2$):
$$\begin{split}J&=\frac{2\varepsilon_0n_0c\sigma_2I^2\hbar\omega_0\varrho_{nt}}{I}E\\
&=\varepsilon_0^2n_0^2c^2\sigma_2\hbar\omega_0E^3\end{split}$$
with $\sigma_2$ the absorption cross section which can be reduced to the equation given above:
$$J=c_1\beta E^3(z)$$
For solving the equation I can split it into the linear and the non-linear part, and start by solving the non-linear part using a RK4-solver, and applying the linear part afterwards. Implementation of the right-hand side function is similar as described above, except by the addition of an FFT for transferring $\hat{E}\rightarrow E$ and back, before and after the calculation of $J$.

Still, as soon as $J$ reaches a value $J\geq J_{lim}$ (which occurs for either large values of $E$ or larger values for $\sigma_2$, i.e. a large absorption), the problems described above occur, and the closer $J$ gets to $J_{lim}$, the smaller I have to choose my step size, to prevent $J$ becoming larger than $J_{lim}$, which is a behavior I still do not understand.


问:
Choosing good modelling method for solving Boltzmann equation
I'm writing a solver for Boltzmann Equations (BE) including a force term in rarefied plasma, for my PhD. The aim is to see if an instability occurs inside an electric streamer (theoretically it should, basing on the analysis of the dispersion relation).

I used finite differences method with MPI parallelization, but the model is unstable before the instability shows. Unfortunately, the methods that I used are slow. I'm thinking of changing the paradigm and turning to Particle-In-Cell (PIC) model. Anyway there still are other methods like Lattice Boltzmann Method (LBM) that might be used without reducing particles to macro-particles as in PIC while improving the speed.

Have you got any recommendations for me where to turn with my investigations? Maybe it would be useful to go few steps higher, trying to see if hydrodynamic approach would do? However, the analysis linked is based purely on kinetic theory. I note that I got the access to huge computing resources.

The set of BE being solved is:

On the right hand side there is a BGK collision term ( for $\alpha$ type of particles, $\Phi_\alpha$ is a Maxwellian distribution function shifted by the velocity of a streamer), r is displacement vector, v is velocity vector, E is a constant field that acts on the whole fluid, but E’ and B’ are locally-variable fields that are related to the small wavy plasma disturbance calculated from Maxwell equations in each step.


问:
Preconditioning the $[1 \quad-2 \quad 1]$ Finite Difference matrix
Let $A$ be the well known tridiagonal matrix coming from the 1D Finite difference discretization of the Laplacian, with stencil $\frac{[1 \quad-2 \quad 1]}{h^2}$.

The system $Ax = b$ is very large, so an iterative method must be used. Also, it is symmetric and positive definite, so I could employ PCG to solve it. The big question is: what preconditioner should I set? I've searched on the web but, surprisingly, I couldn't find anything.

Any reference, MatLab/C++ code or anything in between would be highly appreciated.

Best regards


问:
Convergence of Truncated Newton for non-convex Hessian
I was wondering if anyone could enlighten me about the convergence properties of the truncated newton method in case of a non-positive definite hessian $\nabla^2 f = H$. From the Book 'Numerical Optimization' from Nocedal & Wright, the Algorithm consists of two loops, one inner to compute the (Newton-)search direction $p_k$ by using the conjugate gradient method (CG) and one outer loop to compute the stepwise $\alpha_k$.

The algorithm roughly looks like this:

Outer Loop For $k = 0,1,2,...$:

Given some initial $x_0$, $r_0 = \nabla f_k$ and $p_0 = -\nabla f_k$

Define stopping tolerance of CG $\eta_k=min$ $(0.5,\sqrt \Vert \nabla f \Vert)$ $\Vert \nabla f \Vert$ for superlinear convergence.

Inner Loop For $j = 0,1,2,...$:

Compute $p_k$ by solving $Hp_k$ = $-\nabla f$ using the CG.

If $\Vert r_j \Vert \leq \eta_k$, stopping criteria is fulfilled and $p_k$ is found.

If during CG the Hessian loses positive definiteness $p_j^T H p_j \leq 0$ than use last valid computed direction $p_{j-1}$ where  $p_j^T H p_j > 0$. If that happens during first iteration $j=0$, use the gradient as new direction $p_k = -\nabla f$.

Apply some line search to determine $\alpha_k$ and set $x_{k+1} = x_k + \alpha_k p_k$

Now as mentioned, since  $\underset{k \rightarrow \infty}{lim}\eta_k = 0$ the convergence of truncated newton method is q-superlinear if $H$ is positive definite. But if the inner loop gets disrupted by $H$ becoming non-positive definite, CG will not fullfill the stopping criteria. Especially if non-convexity gets detected in the first iteration and $p_k = -\nabla f$, truncated newton method basically becomes method of greatest descent which has only linear convergence properties. Am I therefore right to assume that truncated newton only has superlinear convergence if $H$ stays positive definite for all times ($k=0,1,2,...$)?


问:
Discretization formula for a system of two differential equations. "Solution to one of these is the initial condition of the other". In which sense?

Consider the following stochastic differential equation
\begin{equation}
dy=\left(A-\left(A+B\right)y\right)dt+C\sqrt{y\left(1-y\right)}dW\tag{1}
\end{equation}
where $A$, $B$ and $C$ are parameters and $dW$ is a Wiener increment.

Equation $(1)$ will be our point of reference in what follows.

Now, first let us consider a "method" for equation $\left(1 \right)$ which can be described by the following one-step discretization scheme:
\begin{equation}
y_{n+1}=y_n+\left(A-\left(A+B\right)y_n\right)\Delta t +C\sqrt{y_n\left(1-y_n\right)}\Delta W_n + D\left(y_n\right)\left(y_n-y_{n+1}\right)\tag{2}
\end{equation}
where $\Delta t$ is the length of the time discretization interval, $\Delta W_n$ is a Wiener increment and $D(y_n)$ is the system of control functions and takes the form
$$
D(y_n)=d^0(y_n)\Delta t + d^1\left(y_n\right)|\Delta W_n|
$$
with
$$
d^1(y)=
\begin{cases}
C\sqrt{\frac{1-\varepsilon}{\varepsilon}}\hspace{0.5cm}\text{if }y<\varepsilon\\
C\sqrt{\frac{1-y}{y}}\hspace{0.5cm}\text{if }\varepsilon\le y<\frac{1}{2}\\
C\sqrt{\frac{y}{1-y}}\hspace{0.5cm}\text{if }\frac{1}{2}\le y\le 1-\varepsilon\\
C\sqrt{\frac{1-\varepsilon}{\varepsilon}}\hspace{0.5cm}\text{if }y>1-\varepsilon
\end{cases}
$$
At this point, let us consider a "method" which decomposes $\left(1\right)$ into two equations. Specifically, the first equation is a stochastic one, that consists of the diffusion term of $\left(1\right)$ only (see eqtn $\left(3\right)$), while the second one is an ordinary differential equation (see eqtn $\left(4\right)$) that consists of the drift part of $\left(1\right)$. We have:

$\begin{equation}
dy_1=C\sqrt{y_1\left(1-y_1\right)}dW\tag{3}
\end{equation}$
$\begin{equation}
dy_2=\left(A-\left(A+B\right)y_2\right)dt\tag{4}
\end{equation}$

This last method approximates the solution to $\left(3\right)$ at each time step using $\left(2\right)$ (and numerical solution to $\left(3\right)$ is used as the initial condition in $\left(4\right)$), while $\left(4\right)$ can be solved using the Euler method. Thus, such a method can be described by the following one step discretization formula:
$$
y_{n+1}=y_n+\left(A-\left(A+B\right)y_n\right)\Delta t + \dfrac{C\sqrt{y_n\left(1-y_n\right)}\Delta W_n}{1+d^1\left(y_n\right)|\Delta W_n|}\left(1-\left(A+B\right)\Delta t\right)\tag{5}
$$

My doubts:

I cannot understand in which way the last method approximates solution to $\left(3\right)$ at each time step using $\left(2\right)$. Could you please explicit such an approximation? How is it obtained by means of $\left(2\right)$?

In which sense numerical solution to $\left(3\right)$ is used as the initial condition in $\left(4\right)$? Which is such an initial condition?

Could you please explicit the way in which solution to $\left(3\right)$ and solution to $\left(4\right)$ are combined so as to obtain discretization formula $\left(5\right)$?


问:
OpenMP inconsistently segfaulting
I parallelized my code with openMP, and now have this bug in my code that's really odd. If it's the first time I run on a computer, it segfaults, but if I run the program again it runs fine. I've never seen anything like this and was curious if anyone has seen this before or has good resources or hints on how to fix this. It's a very simple loop that it segfaults on, and I parallelized every loop with default none. The variable declaration followed by the loop is below. I had previously posted a stripped down version of the loop for readability, but it was requested that I fill it in for the sake of more clarity :

```<code> real(dp)    :: q, oneOverq, L, D, Pressure, u, v, En, rho, oneOverrho
 real(dp)    :: normal(Mesh%numDims), wL(Mesh%numFields), cd_dql(Mesh%numFields, Mesh%numEdges)
 real(dp)    :: norm_dx1(Mesh%numDims, Mesh%numDims), norm_dx2(Mesh%numDims, Mesh%numDims)
 real(dp)    :: L_dx(Mesh%numDims, Mesh%numNodes), D_dx(Mesh%numDims, Mesh%numNodes)
 real(dp)    :: p_dql(Mesh%numFields), L_dql(Mesh%numFields, Mesh%numEdges), Res(Mesh%numFields, Mesh%numTri)
 INTEGER(i4) :: i, j, k, m, nodecount, cellNum, edgeNum, Node1, Node2, node, cell, cells(3)

   !<span class="math-container">$OMP PARALLEL DO Reduction(+:L, D, cL_dql, L_dx, cD_dql, D_Dx) default(none) shared(Mesh, W, grad, gamma, dRdX, Jac_2nd) &
   !$</span>OMP private(cellnum, edgenum, node1, node2, rho, oneoverrho, u, v, En, pressure, normal, norm_dx1, norm_dx2, p_dql, wL, cell, node, cells)
   !DIR$ IVDEP:LOOP
   do cellnum = 1, Mesh%numTri
       edgeNum = Mesh%airfoilEdges(i)
       Node1   = Mesh%edgelist(edgenum,1)
       Node2   = Mesh%edgelist(edgenum,2)
       cellNum = Mesh%edgeList(edgenum,3)
       call cell_stencil(Mesh, cellnum, cells)

       wL = W(:, cellnum)

       call linear_limited_reconstruction(cellnum, edgenum, Mesh, grad(:,:,cellNum), wL)

       rho = wL(1)
       oneOverrho = 1.0_dp/rho
       u   = wL(2)*oneOverrho
       v   = wL(3)*oneOverrho
       En  = wL(4)*oneOverrho

       call get_pressure(rho, En, u, v, pressure)

       p_dql(1) = (gamma - 1._dp)*half*( u**2 + v**2 )
       p_dql(2) = (gamma - 1._dp)*-u
       p_dql(3) = (gamma - 1._dp)*-v
       p_dql(4) = gamma-1._dp

       call norm_dx(node1, node2, Mesh, normal, norm_dx1, norm_dx2)

       L     = L + Pressure*normal(2)
       D     = D + Pressure*normal(1)

       do j = 1, 4
           if (j == 1) then
               cell = cellnum
           else
               cell = cells(j-1)
           end if
           if (cell /= 0) then
               cL_dql(:,cell) = cL_dql(:,cell) + matmul(P_dql(:),Jac_2nd%wL_rc_dq(:,:,j,edgenum))*normal(2)
               cD_dql(:,cell) = cD_dql(:,cell) + matmul(P_dql(:),Jac_2nd%wL_rc_dq(:,:,j,edgenum))*normal(1)
           end if
       end do

       L_dx(:, node1)   = L_dx(:, node1) + Pressure*norm_dx1(2, :)
       L_dx(:, node2)   = L_dx(:, node2) + Pressure*norm_dx2(2, :)
       D_dx(:, node1)   = D_dx(:, node1) + Pressure*norm_dx1(1, :)
       D_dx(:, node2)   = D_dx(:, node2) + Pressure*norm_dx2(1, :)

       do j = 1, 6
           node = Mesh%NodeStencil(j, cellnum)
           if (node /= 0) then
               do k = 1, Mesh%numFields
                   L_dx(:,node) = L_dx(:,node) + P_dql(k)*dRdX%wL_rc_dx(k,:,j,edgenum)*normal(2)
                   D_dx(:,node) = D_dx(:,node) + P_dql(k)*dRdX%wL_rc_dx(k,:,j,edgenum)*normal(1)
               end do
           end if
       end do

   end do
   !$OMP END PARALLEL DO
</code>```

I can put in the full code up in the code block, but for the sake of readability I thought it was best not to. If anyone has any advice I'd greatly appreciate it.

Thanks!

Title: OpenMP inconsistently segfaulting


问:
What is the generalization of the resource allocation problem I'm dealing with here?
I'm dealing with a problem as follows:

I have a finite set of money  to spend over  different raffles, and I need to spend approximately to my budget, with the goal of maximizing my probability of my total winnings across all raffles being a constant . Each raffle has a different prize amount , and there is a certain amount of time  between now and when the raffle winner is drawn, and at each interval a random number of people  enter each raffle(which I have perfect information of), and a random number of people refund their raffle ticket . So, the probability I win each raffle can go up or down each temporal period. I start the scenario with a differing number of raffle tickets in each game, and I have a defined value function for each raffle so I know the increase I get in probability of winning the raffle based on the amount of money I spend buying raffle tickets for it.

I've looked into multi-armed bandit problems and knapsack problems but I don't feel the problem perfectly fits into either bucket. Are there multiple problems at hand here, and that's why it doesn't appear to fit into any of those categories?


问:
FEM Meshing artifact at nodes with fewer neighbors
I wrote a 2D-FEM solver to solve some diffusion process and wanted to verify my code with a test problem.
The input was $f(x,y) = x^2+y^2$ and I applied the stiffness matrix on it to get $\Delta f = 4$.
The lumped mass matrix $M$ and stiffness matrix $D$ are assembled from an analytical solution of triangles.

Unfortunately, I get too high values at nodes with fewer neighbor elements (5 instead of six) than average (see below), whereas elements with more neighbors underestimate the analytical value.
I also tried a finer resolution, but this didn't help.

<img src="https://i.stack.imgur.com/qsQM4.png" alt="solution of $\Delta f = 4$ for $f(x,y) = x^2+y^2$ " />

Note that I manually added the value 4 at the boundary nodes, since the corresponding entries in the stiffness matrix are 0.

Do you have an idea, what the problem might be?
In my test simulations for the diffusion process, I couldn't find this meshing artifact. Only in the test problem.


问:
Generally quadratic convex problem with one non-convex term
How would you approach a standard convex quadratic problem  with convex constraints but one non-convex term ? Say $|x|^{0.4}$.

$$\min_x \frac{1}{2} x^{T}Qx + g^Tx + c^T \mathrm{sign}(x) |x|^{0.4}  $$
subject to $$Ax\leq b$$

Where $|\cdot|, \mathrm{sign}()$ are element-wise.

Is there any other way to solve the problem within convex framework ? What`s the best way to approximate if the dimensions are large ?

This is similar to Reformulate a strictly convex QP problem containing absolute value term but here there is a power term.


问:
implementation of shell elements in a topology optimization algorithm
I am working on developing a topology optimization solver based on the finite element method and I want to add a triangular shell element in it. I used the classic finite element method but I didn’t separate in-plane and off-plane strains. My strain-displacement matrix consists of the following:

$$
B_i=\left(
\begin{array}{cccccc}
\frac{\partial N_i}{\partial x} & 0 & 0 & 0 & z\frac{\partial N_i}{\partial x} & 0\\
0 & \frac{\partial N_i}{\partial y} & 0 & -z\frac{\partial N_i}{\partial y} & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0\\
\frac{\partial N_i}{\partial y} & \frac{\partial N_i}{\partial x} & 0 & -z\frac{\partial N_i}{\partial x} & z\frac{\partial N_i}{\partial y} & 0\\
0 & 0 & \frac{\partial N_i}{\partial y} & 0 & \kappa N_i & 0\\
0 & 0 & \frac{\partial N_i}{\partial x} & -\kappa N_i & 0 & 0\\
\end{array}\right)
$$

And the shape function is:

$$
N(\xi,\eta)=\begin{cases}
N_1=1-\xi-\eta\\
N_2=\xi\\
N_3=\eta
\end{cases}
$$

I tried to validate my implementation with an example a cantilever beam: fixed on the left side and a concentrated force applied in the lower right corner and here is my Von Mises stresses distribution:

And my topology optimization is (optimal criteria used):

The results I obtain from Abaqus with the same properties, for Von Mises stresses distribution are depicted as follows:

The results that I get from my code are different from the results of Abaqus and I am not sure why. But here are their solutions and modifications I tried without any tangible difference in the obtained results:

I used a transformation matrix otherwise the calculation doesn’t run

The transformation matrix applied only in the calculation of the stiffness matrix (its dimension is $18\times 18$)

I tried to formulate the shape function in $(x,y)$ basis instead of $(\xi,\eta)$ but there was no difference

What do you think might be the source of error in my implementation? If you need more information about the way I implemented my code please tell me to clarify.


问:
How to speed up the Mixed-Integer Quadratic Program process?
Currently, I am solving a problem in the format:

M is an integer as well. The problem that troubles me is that X is a vector in {0,1} with a size of 7000. I use the solver in https://github.com/oxfordcontrol/miosqp. It seems they are using the branch-and-bound algorithm. And they claim the speed is comparable to Gurobi. Is there a way to speed this calculation process up as it already took me 10 hours and I still can not get a solution. The parameters I can find to set the solver are the "integer feasibility tolerance", the "maximum number of iterations" and whether to use the "depth-first" rule or the "depth-first until the first incumbent and then best bound" rule.

I am totally new to the Mixed-Integer Quadratic Program so I don't know what these parameters mean. Is there a way to speed it up if I don't require a very accurate solution? Of course, I still don't want to solution to be very bad. But I can tolerate the solver to get a solution within an hour.


问:
Explanation of Givens rotation in Jacobi Rotation SVD
I'm trying to implement Singular Value Decomposition (homework of sorts) via the Jacobi Rotation method (more info here, pages 11 and 12).

I am stuck at the bullet saying (sorry for the picture, but I'm not too familiar with writing formulas... yet):

I know how to find the sine and cosine which rotate a vector (say, as explained in Wikipedia), but I don't understand the process above:

how should I find s and c when B is a 2x2 matrix is involved (instead of a 2x1 vector, as in Wikipedia)?

how should I find the s-hat and c-hat in the right matrix?


问:
Derivative-free ill-conditioned non-linear least squares
I am looking for a package which can solve (non-linear) least squares problems without the use of derivatives (because of an expensive model), but which also deals with ill-conditioning well (such as when the Jacobian is sparse).

For example, the DFO-LS (Derivative-Free Optimizer for Least-Squares) solver documentation makes no mention of conditioning, whereas other solvers which use methods such as truncated SVD or Tikhonov regularisation never include good derivative free methods.

If none exist, is it maybe the case that the trust-region framework can already deal well with an ill-conditioned Jacobian and so something like DFO-LS is still worth trying, or is there no good way of solving such a case?