问:
Parallel In Time with Multigrid
I am trying to solve the linear finite element equation $M\ddot{u}+Ku=F(t)$, where $M$ is the mass matrix ,$K$ the stiffness matrix and $F(t)$ the external load vector, parallel in time using XBraid package.Using Newmark's method for time integration we take three equations:

$(K+\frac{M}{dt^2β})u_{n+1}=F_{n+1}+\frac{M}{dt^2β}(u_n+dt\dot{u}_n+(0.5-β)dt^2\ddot{u}_n)$

$\ddot{u}_{n+1}=\frac{1}{βdt^2}(u_{n+1}-u_n-dt\dot{u}_n-dt^2(0.5-β)\ddot{u}_n)$

$\dot{u}_{n+1}=\dot{u}_{n}+dt((1-γ)\ddot{u}_{n}+γ\ddot{u}_{n+1})$

We arrange these three unknowns in a block vector $w_{n+1}=[u_{n+1},\ddot{u}_{n+1},\dot{u}_{n+1}]$ and with some further effort we take something like $$ Aw=f $$ where $w=[w_0\;w_1 \dots w_{n+1}]^T$, $f=[f_0\;f_1 \dots f_{n+1}]^T$ and

$$A=\begin{bmatrix}
    I         \\
    Φ && I  \\
     & &Φ&& \ddots \\
      &&&&\ddots    \\
       &&&&&Φ  &&&I \\
    \end{bmatrix}
$$

with $Φ$ be the time step matrix which take the solution from one time to another.

I tried to implement this idea in Xbraid,however I noticed some strange results. For example I noticed that for very small time-points(50 time points) the residual provided by Xbraid grows at each iteration and at some point decreases rapidly and the method converges giving me the same results as the serial problem. For larger time-steps (over 1000 for example) the residual grows rapidly (takes very large values) and after some point it start to decrease very slowly reaching a point where it stucks and no furhter degradation takes place, leading to non-convergence. Of course this behaviour is not normal in both cases, but I don't know what can cause this behaviour. I provide the basic core of Xbraid implementation in case someone familiar with the software can suggest a solution.

```<code>   /*Time-step Function. Newmark Method is used to take the solution from t->t+1*/

     void Newmark_BetaMethod(double t,double dt,double gamma,double beta,int NN,int NEy,int *dbc,double
         **K_global,double **M_global,double **Matrix_global,double *Fext,double *RHS,Solver *solver,double
                *tempDispl,double *tempVeloc,double *tempAccel,double *displ,double *veloc,double *accel){

    /*create Matrix_global*/
    for(int i=0;i<2*NN;i++){
        for(int j=0;j<2*NN;j++){
            Matrix_global[i][j]=K_global[i][j]+(M_global[i][j])/(beta*dt*dt);
        }
    }

    /*hold values at time step t*/
    for(int i=0;i<2*NN;i++){
        tempDispl[i]=displ[i];
        tempVeloc[i]=veloc[i];
        tempAccel[i]=accel[i];
    }

    /*calculate RHS*/
    for(int i=0;i<2*NN;i++){
        RHS[i]=0.0;
    }

    for(int i=0;i<2*NN;i++){
        for(int j=0;j<2*NN;j++){
            RHS[i]+=(M_global[i][j]*(displ[j]+dt*veloc[j]+(0.5-beta)*(dt*dt)*accel[j]))/(beta*dt*dt);
        }
        RHS[i]+=t*Fext[i];
    }

    /*zero the rows where homogeneous dirichlet_bc refer to*/
    BoundaryCondForRHS(NN,NEy,dbc,RHS);

    /*Find displacement at t+1*/
    solver->ConjugateGradient(Matrix_global,RHS,displ);

    /*Update acceleration*/
    for(int i=0;i<2*NN;i++){
        accel[i]= (displ[i]-tempDispl[i]-dt*tempVeloc[i]-(dt*dt)*(0.5-beta)*tempAccel[i])/(beta*dt*dt);
    }

    /*Update velocity*/
    for(int i=0;i<2*NN;i++){
        veloc[i]= tempVeloc[i]+dt*((1-gamma)*tempAccel[i]+gamma*accel[i]);
    }
}

//include all time-dependent variables->create a block vector
class BraidVector{

    public:
        double *displ; //a displacements vector
        double *veloc; //velocities vector
        double *accel; //acceleration
        BraidVector(double *u,double *v,double *acc);
        virtual ~BraidVector() {};
};
BraidVector::BraidVector(double *u,double *v,double *acc){
    displ=u;
    veloc=v;
    accel=acc;
}

//include all time-independent information here
class MyBraidApp:public BraidApp
{

    protected:

    public:

        MyBraidApp(MPI_Comm comm_t_, int rank_,double **Matrix_global_,double **M_global_,double **K_global_,double *Fext_,double *RHS_,double gamma_,double beta_,double *tempD,double *tempV,double *tempA,int NN,int NEy_,int *dbc,Solver *solver_ ,double tstart_ =0.0,double tstop_=1.0,int ntime_=100);

        int rank;
        double **Matrix_global; //This matrix pointer holds the sum of K+M/beta*dt*dt
        double **M_global;      //global mass matrix
        double **K_global;      //global stiffness matrix
        double *Fext;           //global total load vector
        double *RHS;            //temporary vector that holds t*Fext+(M/beta*dt*dt)*{displ+dt*velocity+(0.5-beta)*dt*dt*acceleration}
        double beta;            //Newmark parameter
        double gamma;           //Newmark parameter
        double *tempDispl;      //temporary vector holds old values of displacements
        double *tempVeloc;
        double *tempAccel;
        int nodes;              //total number of nodes in spatial mesh
        int NEy;                //number of elemenets in y->helps to find the size of boundary nodes
        int *displ_bc;          //holds the dirichlet boundary nodes where u=0;
        Solver *solver;         //pointer to a solver object (needed for solving tthe linear system at each step)

        virtual ~MyBraidApp();

        virtual int Init(double t,braid_Vector *w_ptr);
        virtual int Step(braid_Vector w_,braid_Vector wstop, braid_Vector fstop,BraidStepStatus &pstatus);
        virtual int Clone(braid_Vector w_,braid_Vector *s_ptr);
        virtual int Sum(double alpha,braid_Vector x,double beta,braid_Vector y);
        virtual int BufSize(int *size_ptr,BraidBufferStatus &status);
        virtual int BufPack(braid_Vector w_,void *buffer,BraidBufferStatus &status);
        virtual int BufUnpack(void *buffer,braid_Vector *w_ptr,BraidBufferStatus &status);
        virtual int SpatialNorm(braid_Vector w_,double *norm_ptr);
        virtual int Access(braid_Vector w_,BraidAccessStatus &astatus);
        virtual int Free(braid_Vector w_);

        // Not needed in this example
        virtual int Residual(braid_Vector     w_,
                        braid_Vector     r_,
                        BraidStepStatus &pstatus) { return 0; }

        // Not needed in this example
        virtual int Coarsen(braid_Vector   fw_,
                       braid_Vector  *cw_ptr,
                       BraidCoarsenRefStatus &status) { return 0; }

        // Not needed in this example
        virtual int Refine(braid_Vector   cw_,
                      braid_Vector  *fw_ptr,
                      BraidCoarsenRefStatus &status)  { return 0; }
};

/*Construct the braid app*/
MyBraidApp::MyBraidApp(MPI_Comm comm_t_, int rank_,double **Matrix_global_,double **M_global_,double **K_global_,double *Fext_,double *RHS_,double gamma_,double beta_,double *tempD,double *tempV,double *tempA,int NN,int NEy_,int *dbc,Solver *solver_  ,double tstart_ ,double tstop_,int ntime_)
    :BraidApp(comm_t_, tstart_, tstop_, ntime_)
{
    rank=rank_;
    Matrix_global=Matrix_global_;
    M_global=M_global_;
    K_global=K_global_;
    Fext=Fext_;
    RHS=RHS_;
    beta=beta_;
    gamma=gamma_;
    nodes=NN;
    tempDispl=tempD;
    tempAccel=tempA;
    tempVeloc=tempV;
    NEy=NEy_;
    displ_bc=dbc;
    solver=solver_;

}

MyBraidApp::~MyBraidApp(){
    delete []tempDispl;
    delete []tempVeloc;
    delete []tempAccel;
    delete []Fext;
    delete []RHS;
    delete []displ_bc;
    int size=2*nodes;
    for(int i=0;i<size;i++){
       delete M_global[i];
       delete Matrix_global[i];
       delete K_global[i];
    }
    delete [] M_global;
    delete [] K_global;
    delete [] Matrix_global;
}

/*Initialize the block vector w*/
int MyBraidApp::Init(double t,braid_Vector *w_ptr){
    int size=2*nodes;
    double *displ=new double[size];
    double *veloc=new double[size];
    double *accel=new double[size];
    BraidVector *w=new BraidVector(displ,veloc,accel);
    /*w->displ=new double[size];
    w->veloc=new double[size];
    w->accel=new double[size];
    w->size=size;*/

    //Initial conditions at t=0
    if(t==tstart){
        for(int i=0;i<size;i++){
            w->displ[i]=0.0;
            w->veloc[i]=0.0;
            w->accel[i]=0.0;

        }
    }
    else{
        //for t>0 use random values as initial guess
        for(int i=0;i<size;i++){
            w->displ[i]=1e-2;
            w->veloc[i]=1e-2;
            w->accel[i]=1e-2;
        }
    }
    *w_ptr=(braid_Vector) w;
    return 0;
}

/*Takes a time step using Newmark method*/
int MyBraidApp::Step(braid_Vector w_,braid_Vector wstop, braid_Vector fstop,BraidStepStatus &pstatus){

    BraidVector *w=(BraidVector*)w_;
    double t1,t2;
    // Get time step information->Return tstart and tstop
    pstatus.GetTstartTstop(&t1, &t2);
    double dt=t2-t1; //calculate dt
    /*Call Newmark method to advance in time*/
    Newmark_BetaMethod(t2,dt,gamma,beta,nodes,NEy,displ_bc,K_global,M_global,Matrix_global,Fext,RHS,solver,tempDispl,tempVeloc,tempAccel,w->displ,w->veloc,w->accel);

    return 0;
}

/*Create a new clone vector s of block vector w*/
int MyBraidApp::Clone(braid_Vector w_,braid_Vector *s_ptr){
    int size=2*nodes;
    BraidVector *w = (BraidVector*) w_;
    double *displ=new double[size];
    double *veloc=new double[size];
    double *accel=new double[size];
    BraidVector *s=new BraidVector(displ,veloc,accel);
    /*s->displ=new double[size];
    s->veloc=new double[size];
    s->accel=new double[size];
    s->size=size;*/
    for(int i=0;i<size;i++){
        s->displ[i]=w->displ[i];
        s->veloc[i]=w->veloc[i];
        s->accel[i]=w->accel[i];
    }
    *s_ptr=(braid_Vector) s;
    return 0;
}

int MyBraidApp::Sum(double alpha,braid_Vector x,double beta,braid_Vector y){
    int size=2*nodes;
    BraidVector *vec1 = (BraidVector*) x;
    BraidVector *vec2 = (BraidVector*) y;
    for(int i=0;i<size;i++){
        vec2->displ[i]=alpha*(vec1->displ[i])+beta*(vec2->displ[i]);
        vec2->veloc[i]=alpha*(vec1->veloc[i])+beta*(vec2->veloc[i]);
        vec2->accel[i]=alpha*(vec1->accel[i])+beta*(vec2->accel[i]);
    }

    return 0;

}

/*Create the buffer's size that holds the data of braid vector*/
int MyBraidApp::BufSize(int *size_ptr,BraidBufferStatus &status){
    int buffsize=3*(2*nodes);  //buffsize holds the size of the block vector w =size of displ+ size of veloc +size of accel
    *size_ptr=buffsize*sizeof(double);
    return 0;
}

/*pack displacement,velocity and acceleration vector into a buffer->This buffer is send via proccesses*/
int MyBraidApp::BufPack(braid_Vector w_,void *buffer,BraidBufferStatus &status){
    BraidVector *w=(BraidVector*)w_;
    double *dbuffer=(double *)buffer;
    int size=2*nodes;
    for(int i=0;i<size;i++){
        dbuffer[i]=w->displ[i];
        dbuffer[size+i]=w->veloc[i];
        dbuffer[2*size+i]=w->accel[i];
    }
    status.SetSize((3*size)*sizeof(double));
    return 0;
}

/*Create a new braid vector,unpack the data, and initialize it*/
int MyBraidApp::BufUnpack(void *buffer,braid_Vector *w_ptr,BraidBufferStatus &status){

    int size=2*nodes;
    double *dbuffer = (double *) buffer;
    double *displ=new double[size];
    double *veloc=new double[size];
    double *accel=new double[size];
    BraidVector *w=new BraidVector(displ,veloc,accel);
    /*w->displ=new double[size];
    w->veloc=new double[size];
    w->accel=new double[size];*/
    for(int i=0;i<size;i++){
        w->displ[i]=dbuffer[i];
        w->veloc[i]=dbuffer[size+i];
        w->accel[i]=dbuffer[2*size+i];
    }

    *w_ptr = (braid_Vector) w;

    return 0;
}

/*Euclidean norm is used for convergece*/
int MyBraidApp::SpatialNorm(braid_Vector w_,double *norm_ptr){

    int size=2*nodes;
    BraidVector *w=(BraidVector*)w_;
    double norm=0.0;
    for(int i=0;i<size;i++){
        norm+=(w->displ)[i]*(w->displ)[i]+(w->veloc)[i]*(w->veloc)[i]+(w->accel)[i]*(w->accel)[i];
    }

    *norm_ptr=sqrt(norm);
    return 0;
}

/*Clean up*/
int MyBraidApp::Free(braid_Vector w_){

    BraidVector *w=(BraidVector*)w_;
    delete []w->displ;
    delete []w->veloc;
    delete []w->accel;
    delete w;

    return 0;
}

int MyBraidApp::Access(braid_Vector w_,BraidAccessStatus &astatus){
    char filename[255];
    FILE *file;
    BraidVector *w = (BraidVector*) w_;

    int size=2*nodes;
    int done,iter,index,level;
    double t;
    astatus.GetTILD(&t, &iter, &level, &done);
    astatus.GetTIndex(&index);

    // Print information to file
    sprintf(filename, "%s.%07d.%05d", "Plane stress", index, rank);
    file = fopen(filename, "w");
    //fprintf(file, "%d\n",    Npoints);
    fprintf(file, "%.14e\n", tstart );
    fprintf(file, "%.14e\n", tstop );
    //fprintf(file, "%d\n",    size );

    fprintf(file, "%s\n" , "*Displacements*");
    for (int i=0;i<size;i++){
      fprintf(file, "%.14e\n", w->displ[i]);
    }

    fprintf(file, "%s\n" , "*Velocities*");
    for (int i=0;i<size;i++){
      fprintf(file, "%.14e\n", w->veloc[i]);
    }

    fprintf(file, "%s\n" , "*Accelerations*");
    for (int i=0;i<size;i++){
      fprintf(file, "%.14e\n", w->accel[i]);
    }

    fflush(file);
    fclose(file);
    return 0;
}
</code>```

I need to mention that I ran Xbraid sequentially as mentioned in the manual and the results were exactly the same as my serial code.

Title: Parallel In Time with Multigrid


问:
Optimize linear equation using inner products and subject to L1 norm
I have a linear system of the form $A x = b$ where $A$ and $b$ are known, $A$ is "square", and $\lvert b \rvert_1 = \lvert x \rvert_1 = 1$. Unfortunately, I am working in a framework that only gives me access to a dot product, not the individual vector or operator elements. The only way I have found to solve this problem is by optimizing the quantity:

$$ d^2 = \lvert b - A x \rvert_2 = \langle b, b \rangle - 2 \langle b, A x \rangle + \langle A x, A x \rangle $$

subject to the normalization $\lvert x \rvert_1 = 1$. This tends to work, but it becomes computationally challenging to compute the last term because of the double operator. (Forming the vector $A x$ is actually very difficult, so I have to compute the inner product "all at once" to use this toolset: $\langle A x, A x \rangle = x^T A^T A x$.)

Now, if this were an $L^2$ constrained problem, I could use ${d^\prime}^2 = 1 - \langle b, A x \rangle$ as the "overlap" between states $b$ and $A x$, which makes more sense to me coming from a quantum mechanics perspective, yet, because $x$ is only $L^1$ normalized, I would still have to ensure that $b$ and $A x$ are $L^2$ normalized before optimizing this way -- which takes me back to computing $\langle A x, A x \rangle$ at every step, which is, again, too expensive. The operator $A$ does guarantee that the $L^1$ norm is preserved, but the $L^2$ norm is, in general, not.

I can't think of another approach to solve this equation using only inner products. Have I missed something obvious? Or is there literature out there about solving this problem? Google has yet to be helpful on this issue, because "$L^1$ optimization" is a common search and does not seem to be related to my question, since $A$ in that problem is apparently sparse. Please explain any answers at an advanced-undergrad-core-math level, because I am not a mathematician.

ContentLicense: CC BY-SA 4.0


问:
discretizing advection equation with variable wave speed + stability
I currently have a code that solves $u_t+ cu_x=0$ with periodic boundary conditions, and constant $c$ (using an upwind method). I'm wondering how I would alter this code to solve something of the form
$$u_t +c(x) u_x=0.$$ I'm thinking it will look of like

$$u_i^{j+1}=u_i^j-c_i \Delta t \left(\frac{u_i^j-u_{i-1}^j}{\Delta x} \right),$$
where $c_i=c(x_i)$?

I know that generally, the upwind scheme is stable for $c \frac{\Delta t}{\Delta x} \leq 1$. How would this change with a variable $c$? Am I better off using another method?


问:
Error in Monte Carlo integration
I am looking for a concise description to help me understand the error for Monte Carlo Integration using Uniform Sampling and Importance Sampling

For Importance Sampling I have that the error is just the square root of the variance...

$$\sqrt{\frac{1}{N}<f(Xi)^2>-<\frac{1}{N}f(Xi)>^2)}$$

for Uniform Sampling what would the error be similar such as

$$\frac{b-a}{\sqrt{N}}\sqrt{\frac{1}{N}<f(Xi)^2>-<\frac{1}{N}f(Xi)>^2)}$$

where $a$ and $b$ are the limits of the integration??

This to me isn't right as my understanding was that an importance sampling would decrease the variance of the MC integration.

I can't find any clear explanation of the difference online so any help is appreciated.

Title: Error in Monte Carlo integration


问:
Multigrid Reduction In Time Convergence
I am trying to solve a 2D dynamic linear elasticity model parallel in time using Xbraid. The spatial domain is [0,1]x[0,1] and time domain [0,1]. For time integration I am using a backward Euler scheme and FEM for spatial descretization. I tried to run a simulation for a small spatial problem(4x4 quadrilateal elements) and 10000 time-steps and I noticed a weird behaviour of Xbraid's residual(a flactuation). At the beginning the residual decreases, for example starts from 3.24e-01 decreases for the next 30 iterations until 2.45e-04 .Then for the next 40 iterations there is a small increase until 2.09e-03 and the decreases again. Then the residual is decreasing normally until convergence. What may cause this behaviour?Is this normal?
I want to note that I am using a standard coarsening and a FCF relaxation at all levels.V-cycle is used.

A photo of the residual is uploaded.

Title: Multigrid Reduction In Time Convergence


问:
3d surface image for given datasets

I have some data that looks like as follows

x-cor  y-cor  z-cor

0.02  0.0251  0.02

0.01  0.0257  0.02

0.014  0.02   0.03

0.023  0.003  0.05

0.013  0.002  0.77

....................

....................

I would like to draw an image like the above. Can you please suggest a way to start that in Matlab/python?

image source: https://www.pnas.org/content/pnas/109/37/14761.full.pdf

Thanks a lot in advance!

Title: 3d surface image for given datasets


问:
Algorithm for computing inner products multiple times
I am taking a computational linear algebra course and i got stuck during a homework problem concerning the computation of inner products. I am supposed to compute the inner product:$$\mathrm{a}_{\mathrm{k}}=\frac{\left(\mathrm{g}_{\mathrm{k}} \mid \mathrm{f}\right)}{\left(\mathrm{g}_{\mathrm{k}} \mid \mathrm{g}_{\mathrm{k}}\right)}$$ for several values of k.

The inner product is defined as $$(U \mid V)=\int_{-1}^{1} U(x) V(x) d x$$The function $$ f = f(x)=\frac{1}{1+25 x^{2}}$$ and $g_{k}$ are the Legendre polynomials in the interval [-1, 1].

Is there any built-in command in NumPy that could be useful? I tried symbolic manipulation using SymPy but it did not go well. After that, i wrote the basic loop together with the relevant matrices, but i still cannot compute those inner products.

Any help is appreciated. Thanks in advance!


问:
Explicit DG time step restriction for compressible Navier-Stokes equations
Hesthaven's book 1 mentions the following time step restriction for Navier-Stokes equations (see (7.32) in 2008 edition)
$$
\Delta t \approx \frac{h}{N^2} \frac{C}{|u| + |c| + \frac
{N^2 \mu}{h}}
$$
($h$-element size, $N$-degree of polynomial basis, $C$-constant)

Isn't this equation dimensionally incorrect? $\mu/h$ doesn't have dimensions of $u$. I think $\mu$ should be replaced by the kinematic viscosity $\nu$. However, there is nothing mentioned about this in "corrections" either (2,3).

Can anyone help me with giving the correct form of this?

Edit

$\Delta t$ - time step

$h$ - characteristic element size (usually radius of inscribed or circumscribed circle)

$N$ - degree of polynomial approximation

$C$ - constant, $O(1)$

$u$ - fluid velocity in $x$ direction

$c$ - speed of sound

$\mu$ - dynamic viscosity

$\nu$ - kinematic viscosity


问:
Integration of a discretized field in a cylindrical coordinate system
I would like to integrate a discretized field in cylindrical coordinates, given as A(r, z), with z being spaced regularly (spacing distance dz), and r being spaced non-regularly. Moreover, r[0]>0 due to the chosen discretization. For integration, I would like to use the trapezoidal approach, and start with the z-component, to integrate everything into a vector aux:

```<code>for (int i = 0; i < Nr; ++i) {
    for (int j = 0; j < Nt; ++j) {
        double I = A(i, j);
        if ((j == 0) || (j == Nt - 1))
            I *= 0.5;
        aux(i) += I;
    }
    aux(i) *= dz;
}
</code>```

For the integration along the r-axis I then have to consider the different spacing, resulting in

```<code>double integral = 0;
for (int i = Nr; i > 0; --i) {
    integral += (r(i) * aux(i) + r(i - 1) * aux(i - 1)) / 2
            * (r(i) - r(i-1));
}
</code>```

Finally, after r[0] > 0 I also have to take that into account by adding

```<code>integral += 0.5 * aux(0) * r(0) * r(0);
</code>```

and add the missing $2\pi$:

```<code>integral *= 2 * M_PI;
</code>```

Is that approach correct, or did I miss a factor here somewhere?


问:
Finite volume reconstruction techniques
For a cell-centered finite-volume calculation, with a nontrivial stencil (either using an unstructured grid or a strongly non-uniform structured grid), what are the main techniques for reconstruction of function values and derivatives on cell faces? What are their advantages and disadvantages, e.g., compared to the polynomial interpolation? Is there a book or a recent review article covering this?</p>

Title: Finite volume reconstruction techniques


问:
Construction of Prolongation and Restriction Operator for Geometric Multigrid (2D-FEM): Resulting in a Decreasing Solution
Consider the following problem,
$$ -\Delta u(x) = f(x), \qquad x \in \Omega \\ u(x) = 0,\qquad x \in \partial \Omega$$
with $\Omega = [0,1]\times [0,1]$ being the domain and $\partial \Omega$ being the boundary. The source term $f(x) \in L^2(\Omega)$ is a smooth source term. The exact solution is known.

I have generated the two following triangulations of the domain, dimensioned $81\times 81$ and $289 \times 289$.

<img src="https://i.stack.imgur.com/zoayW.png" alt="$81\times 81$ nodes" /><img src="https://i.stack.imgur.com/5BMdx.png" alt="$289 \times 289$ nodes" />

Using regular P1-FEM, I have constructed a linear system for each of these two meshes,
$$ A_1 u_1 = f_1, \qquad A_2 u_2 = f_2.$$

The direct solutions of the two systems are working as intended, with the finer mesh having a relative error (of $L^2$-norm) of $0.0077$ compared to the known exact solution and the coarser mesh being a bit less accurate.

Now, rather than using a direct solver, I am attempting to construct a two-grid solution for this problem. I am using a V-cycle method, with the following structure:

Start with the fine grid, and perform 3 iterations of the dampened Jacobi method (weight $\omega = 2/3$), using zero as initial guess. The damped Jacobi method I have implemented has the following structure
$$ u^{n+1} = u^n + \omega D^{-1} (f - Au^n),$$ where $D$ is the diagonal of $A$. This seem to work as intended, as the $L^2$-norm  of the relaxed solution is lower than the $L^2$-norm for the initial guess. We can denote this relaxed solution by $u_2$.

Compute the residual of the relaxed solution by $r_2 = A_2 u_2 - f_2$.

Restrict the residual $r_2$ to the coarse mesh. Here, I am using an injection restriction operator, where I simply remove the finer nodes and keep the nodes which are part of both meshes. The general structure of this operator is
$$ R = \begin{bmatrix}1 & 0 & 0 & \cdots & & & 0\\ 0 & 1 & 0 & \cdots & & & 0\\ 0 & 0 & \ddots & \cdots & & & 0 \\ 0 & 0 & 0 & 1 & 0 & \cdots & 0\end{bmatrix}, \qquad (81 \times 289)$$
The implementation of this operator seem to work as I intend it to. The dimensions become correct, and the "finer" nodes are removed from the residual. We can denote the restricted residual $r_2$ as $Rr_2$.

Solve the coarse system directly. Using Matlab notation and denoting the solution as $e$, $e = A_1 \setminus (Rr_2)$. This results in a $L^2$ error which is larger than the relaxed solution $u_2$, which seems odd.

Apply a prolongation operator to the solution $e$. The prolongation operator that I am using is a linear interpolation. The nodes which are shared between the coarse grid and the fine grid are left unchanged, while the finer nodes are computed as the average of the two closest coarse nodes (on the same edge). This results in the following structure of the prolongation operator, $$ I = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & \ddots & \cdots & 0 \\ 0 & 0 & \cdots & & 1 \\ 0.5 & 0.5 & 0 & \cdots & 0 \\ 0 & 0 & 0.5 & 0.5 & 0\\
\vdots & \vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0.5 & 0 & 0.5\end{bmatrix}, \qquad (289 \times 81).$$ Again, the dimension of this prolongation operator seem correct, and the vectors seem to update as intended. We can denote $u_1 = Ie$ as the prolonged solution $e$ onto the finer grid.

I update the original relaxed solution $u_2$ as $u_2 = u_2 - Ie$.

Lastly, I apply the same dampened Jacobi iteration (same as step 1), with the solution $u_2$ as the initial guess. The final solution $u$ is obtained.

I think that this is a reasonable approach, and the same principle works for my test program (1D-FDM). Now, the issue that I am facing is that the solution after performing these steps is much lower than the corresponding reference solution.

The following two Figures illustrate the issue well,

The same tendencies are shown in both solution, but the two-grid solution is much smaller than the direct solution.

Comparing the $L^2$-errors of each solution, the following result is shown:

Two-grid solution: 0.0079995 (rel. error 0.24)

Direct solution: 0.000256 (rel. error 0.0077)

This is a huge difference. Adding more levels, or increasing the number of cycles (the principle above is only one cycle) makes the issue worse.

My question is, what am I doing wrong? This is a very odd behavior, but the implementation of the restriction and prolongation operators agree with what I have come up with using pen and paper.

My first thought is that there is some kind of scaling factor missing, since both solution have the same tendencies but have different sizes. I would assume that it is missing in the restriction operator, but adding any kind of factor there seems wrong, since the row-sum should be 1.

I am running out of ideas right now, and would appreciate any suggestion very much!
Thanks in advance!

Edit: If I multiply the restriction operator by a scalar, say $4$, the $L^2$-error of the two-grid solution becomes $0.000612$ (rel.error 0.0184), which is significantly better. But I really cannot understand why this would be a better approach? Could the scalar be related to the number of levels, the ratio of the number of nodes between the levels?

Edit 2: These are the current results, for different attempts on the restriction operator:

Restriction operator based only on injection:

$L^2$ error: 0.0327

Same as (1) but with a factor 4 in front:

$L^2$ error: 0.0022

Restriction operator as transpose of prolongation:

$L^2$ error: 0.0093

Same as (3) but use injection on boundary nodes:

$L^2$ error: 0.0054


问:
Finding a CFD paper with extra degree of freedom variable in mass conservation
I am trying to find a paper that I saw about a year ago. I am not sure of the actual date of the paper.

I believe it was a finite difference CFD paper. The interesting part of the paper was the addition of a single degree of freedom variable to all of the mass conservation equations. This residual mass conservation variable was included to stabilize the iteration of their scheme.

Alternatively, I would be interested in suggestions on handling convergence issues that may be caused by replacing one of the mass conservation equations with a pressure average along the outlet of the fluid domain.

I am also considering switching to the pressure Poisson equations with a projection method to be selected later.


问:
In sights into why higher order finite differencing method leads faster to instability
I was playing around with numerically solving the 1D wave equation with density and stiffness varying with position using central differencing methods and noticed that for certain discretization steps sizes $\Delta t$ and $\Delta x$ the solution could become unstable. Therefore, I thought maybe if I increased the order of the differencing method I could get away with a courser grid, but the opposite seems to be the case.

To make it easier to analyse I simplified the problem to constant material parameters, so the PDE considered can be reduced to

$$
\frac{\partial^2}{\partial t^2} u(t,x) = c^2 \frac{\partial^2}{\partial x^2} u(t,x), \tag{1}
$$

with $c$ the wave speed and boundary conditions such that $u(t,x+L) = u(t,x)$, mainly because this was easiest to implement with both differencing methods. For discretization I used a uniform grid such that

$$
u^n_k = u(n\,\Delta t, k\,\Delta x). \tag{2}
$$

For the integration with respect to time I used Verlet integration, so

$$
\frac{\partial^2}{\partial t^2} u(t,x) = \frac{u^{n+1}_k - 2\,u^n_k + u^{n-1}_k}{\Delta t^2} + O(\Delta t^2). \tag{3}
$$

My initial approach used the same differencing method as $(3)$ for the spacial partial derivatives,

$$
\frac{\partial^2}{\partial x^2} u(t,x) = \frac{u^n_{k+1} - 2\,u^n_k + u^n_{k-1}}{\Delta x^2} + O(\Delta x^2). \tag{4}
$$

Combining $(3)$ and $(4)$ therefore yields

$$
u^{n+1}_k = 2\,u^n_k - u^{n-1}_k + v^2 (u^n_{k+1} - 2\,u^n_k + u^n_{k-1}), \tag{5}
$$

with

$$
v = \frac{c}{\Delta x / \Delta t},
$$

which can be seen as the wave speed normalized with respect to the discretization.

For $v>1$ the system becomes unstable, which also makes sense to me, since waves should propagate more than $\Delta x$ in one time step, but $(5)$ only considers neigboring grid positions.

For my second approach used a higher order differencing method for the spacial partial derivatives,

$$
\frac{\partial^2}{\partial x^2} u(t,x) = \frac{-u^n_{k+2} + 16\,u^n_{k+1} - 30\,u^n_k + 16\,u^n_{k-1} - u^n_{k-2}}{12\,\Delta x^2} + O(\Delta x^4). \tag{6}
$$

Combining $(3)$ and $(6)$ therefore yields

$$
u^{n+1}_k = 2\,u^n_k - u^{n-1}_k + v^2 \frac{-u^n_{k+2} + 16\,u^n_{k+1} - 30\,u^n_k + 16\,u^n_{k-1} - u^n_{k-2}}{12}. \tag{7}
$$

My initial guess would be that $(7)$ would be able to handle values for $v$ up to two, since it now sees two neigboring grid positions, so it can "see" waves coming from twice the distance. However, the opposite seems to be the case. Namely, now the system becomes unstable for $v > \sqrt{3/4} \approx 0.866$. And it is not obvious to me why $\sqrt{3/4}$, so I wonder if there is an intuitive explanation for this?

I determined the stability of each system by formulating it as $x[n+1] = A\,x[n]$, with

$$
x[n] = \begin{bmatrix}u^{n-1}_1 & \cdots & u^{n-1}_N & u^n_1 & \cdots & u^n_N\end{bmatrix}^\top
$$

which is stable if $A$ has all it eigenvalues inside the unit disk.


问:
If dot product is commutative, why does MATLAB give different answers?
Why does the dot() function in MATLAB return different expressions based on the order in which I pass vectors?


问:
Minimal covering of rectangle with fixed, overlapping rectangles
A finite set $R$ of fixed, axis aligned 2-D rectangles $r_i=\left\{x_{0i},y_{0i},W_i,H_i \right\}$ is given. These rectangles are potentially overlapping. Given a new axis aligned rectangle $t$, I would like to find the subset $S \subset R$ with minimum number of elements $M$ such that $t \subset \underset{i=1,\ldots,M}{\bigcup} s_i$. Is this a known computational geometry problem?


问:
Efficient computation of marginalized multivariate normal likelihood
In general,if we know that the marginal Gaussian distribution for some variable $\textbf{x}$ and a conditional Gaussian distribution for some $\textbf{y}|\textbf{x}$ of the forms:
$$p(\textbf{x}) = \mathcal{N}(\textbf{x}|\boldsymbol{\mu},\Lambda^{-1})$$

$$p(\textbf{y}|\textbf{x}) = \mathcal{N}(\textbf{y}|A\textbf{x}+\textbf{b},L^{-1})$$
then the marginal distribution of $\textbf{y}$ and the conditional distribution of $\textbf{x}|\textbf{y}$ are given by:

$$ p(\textbf{y}) = \mathcal{N}(\textbf{y}|A\boldsymbol{\mu}+\textbf{b},L^{-1}+A\Lambda^{-1}A^{T})$$
$$p(\textbf{x}|\textbf{y}) = \mathcal{N}(\textbf{x}| \Sigma \left\{ A^{T} L ( \textbf{y} - \textbf{b} ) + \Lambda \boldsymbol{\mu} \right\} , \Sigma)$$

where :
$\Sigma = (\Lambda + A^{T}LA)^{-1}$

In my situation, I have a prior and likelihood: $$\boldsymbol{\beta} \sim \mathcal{N}(0,I_p\sigma^2_\beta)$$
$$\hat{\boldsymbol{\beta}} | \boldsymbol{\beta} \sim \mathcal{N}(\hat{\textbf{S}}\hat{\textbf{R}}\hat{\textbf{S}}^{-1}\boldsymbol{\beta},\hat{\textbf{S}}\hat{\textbf{R}}\hat{\textbf{S}}).$$
Where $\hat{\textbf{S}}$ is a diagonal matrix and $\hat{\textbf{R}}$ is a symmetric correlation matrix. By making substitutions I arrive at the following marginalized likelihood :

$$\hat{\boldsymbol{\beta}}|\sigma_\beta^2 \sim \mathcal{N}(0,\sigma_\beta^2\hat{\textbf{S}}\hat{\textbf{R}}\hat{\textbf{S}}^{-2}\hat{\textbf{R}}\hat{\textbf{S}}+\hat{\textbf{S}}\hat{\textbf{R}}\hat{\textbf{S}})$$

My question is this, if I want to find the value of $\sigma_\beta^2$ that maximizes the marginalized likelihood, do I have to recompute the cholesky of that entire covariance matrix from scratch every time I try out a new value of $\sigma_\beta^2$, or do I have other options?

If $\hat{\textbf{S}}=\text{diag}(1)$, then we can simply take  the eigenvalue decomposition of $\textbf{R}$:
$$\hat{\boldsymbol{\beta}}|\sigma_\beta^2 \sim \mathcal{N}(0,\sigma_\beta^2\textbf{R}^2+\textbf{R}) = \mathcal{N}(0,\sigma_\beta^2\textbf{Q}\textbf{D}^2\textbf{Q}+\textbf{Q}\textbf{D}\textbf{Q}) \\ \rightarrow \textbf{Q}^T\hat{\boldsymbol{\beta}}|\sigma_\beta^2 \sim \mathcal{N}(0,\sigma_\beta^2\textbf{D}^2+\textbf{D})$$
Where $\textbf{Q}$ are the eigenvectors of $\textbf{R}$ and $\textbf{D}$ is the diagonal matrix of eigenvalues.  That's really easy to compute for different values of $\sigma_\beta^2$.  With $\hat{\textbf{S}}\neq\text{diag}(1)$ None of that works.  My intuition is that I'm trying to get a cheap full rank update of an eigenvalue decomposition, which is in general a no-go, but I would love to hear thoughts.


问:
Formula for overdetermined logical matrix pseudoinverse not requiring SVD?
In https://commons.wikimedia.org/wiki/File:YI_%3D_PI.png, you will find a formula-based solution for an overdetermined logical matrix pseudoinverse. This simple formula gives the same result as the Moore-Penrose method. Does anyone know of any other overdetermined logical matrices that can be solved without using SVD methods?</p>


问:
Phase portrait of non-linear system of ode (Triple Galaxy System)
Context: Hello, I am an undergrad student of physics self-studying Python Programming. I am trying to find the value of H_lam for which limit cycles corresponding to the triple galaxy system governed by the equations
given here occur

Code

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

# Define function
def TGS(t, z):
    x = np.empty(3); u = np.empty(3)
    y = np.empty(3); v = np.empty(3)
    x[0]= z[8]; x[1]= z[0]; x[2]= z[4]
    y[0]= z[9]; y[1]= z[1]; y[2]= z[5]
    u[0]= z[10]; u[1]= z[2]; u[2]= z[6]
    v[0]= z[11]; v[1]= z[3]; v[2]= z[7]
    G = 6.67e-17; m= [2e+43,2e+43,2e+43]
    R = []
    for i in [0,1,2]:
        R1=[]
        for j in [0,1,2]:
            R1.append(np.sqrt((x[j]-x[i])**2+(y[j]-y[i])**2))
        R.append(R1)
    return [u[1], sum((G*m[j]*(x[j]-x[1])/((R[1][j])**3)) for j in [0,2])+(H_lam**2)*x[1],\
            u[2], sum((G*m[j]*(x[j]-x[2])/((R[2][j])**3)) for j in [0,1])+(H_lam**2)*x[2],\
            v[1], sum((G*m[j]*(y[j]-y[1])/((R[1][j])**3)) for j in [0,2])+(H_lam**2)*y[1],\
            v[2], sum((G*m[j]*(y[j]-y[2])/((R[2][j])**3)) for j in [0,1])+(H_lam**2)*y[2]]

# Set limits for evaluation
a, b = 0,1000

# Define H_lams and set the stage for solving the Diff eq
H_lams = [61]
t = np.linspace(a, b, 50000)

# Solve the Diff eq
for H_lam in H_lams:
    sol = solve_ivp(tbs, [a, b], [0.15, -0.4, 0, -90, 0, 0.8, 0, 0], t_eval=t)
    plt.plot(sol.y[0],sol.y[2], "-")

# Make a little extra room for legend
plt.xlim([-1e+10,1e+10])
plt.axis("equal")
plt.legend([f"<span class="math-container">$\H_lam={m}$</span>" for m in H_lams])
plt.axes().set_aspect(1)
</code>```

My thoughts and the problem: x, y are the x and y coordinates while X and Y are the velocities along the x-axis and y-axis respectively. Note that I wanted to keep the first galaxy at origin with zero velocity. The masses for each body, in this case, are equal but I would like to alter them later.

Now, this code that I have written does give me a plot, but I am not sure if it is truly the phase portrait I am looking for. Can someone please suggest a way to verify it?

Moreover, is there a way to make this code more efficient?

Edit: The plot I am getting is


问:
Is it possible to shift the numerical results in order to compare it with the analytical solution?
I'm using lattice Boltzmann method (LBM) to do simulation of flow in a pipe. The diameter of the pipe is 0.02 m and its height is 0.1 m. I put a parabolic velocity profile at the inlet:

$$\mathbf{u}(\mathbf{r},t) = \bar{u}(t) \Bigg(2 - \frac{2|\mathbf{r}-\mathbf{r}_{C}|^{2}}{R^{2}}\Bigg)$$

Where $R$ is the radius of the pipe, $\mathbf{r}_{C}$ is the position of the center of inlet circle, which in my case is: $\mathbf{r}_{C} = (0,0,0)$ and the axis of the pipe is oriented at the z-direction. Finally for $u_{max}(t)$, I have:

$$\bar{u}(t) = \frac{u_{avg}}{2} + \frac{u_{amp}}{2}\sin{(\omega t)}$$

Where $u_{avg} = 0.000124289772727275$ $\frac{\text{m}}{\text{s}}$ and $u_{amp} = 0.00004143$ $\frac{\text{m}}{\text{s}}$ and $\omega = 2 \pi$. This funny numbers for $u_{avg}$ and $u_{amp}$ are chosen because of keeping the Reynolds number at 1. In fact, for density and viscosity, I have: $\rho = 1056$ $\frac{\text{kg}}{\text{m}^{3}}$ and $\mu = 0.0035$ $\text{Pa} \cdot \text{s}$:

$$\text{Re} = \frac{2 \rho R u_{max}}{\mu} = \frac{2 \times 1056 \times 0.01 \times 0.000165719696969700}{0.0035} = 1$$

Note that: $u_{max} = u_{avg} + u_{amp} = 0.000165719696969700$ $\frac{\text{m}}{\text{s}}$.

Because, this is a pipe, I expect to see the same $\bar{u}$ at the outlet after simulation. I run this simulation for 12 seconds in reality, which is equal to 12 cycles of this sinusoidal waveform. My Mach number (i.e. lattice Mach number that has nothing to do with the physical Mach number) and relaxation time as well as my spatio-temporal resolutions are:

$$\text{Mach} = 0.002$$

$$\tau = 0.85$$

$$\Delta x = 0.2 \ \text{mm}$$

$$\Delta t = 1.39 \ \text{ms}$$

Now this the comparison of $\bar{u}(t)$ measured at the outlet versus time:

The blue curve in the above plot is just the analytical $\bar{u}(t)$ defined above and the orange one is the one that I measured at the outlet. You see that there is a time lag or phase shift between the $\bar{u}(t)$ that I measure at the outlet and the analytical solution. In fact, the reason seems to be because of my initial condition that I start from zero velocity everywhere and my simulation needs sometimes to sync outlet with the inlet. In fact, the time that is needed to  transport information from inlet to outlet could be estimated as:

$$t^{*} = \frac{L}{c_{s}} = \frac{0.1}{\frac{1}{\sqrt{3}}\frac{\Delta x}{\Delta t}} = 1.20 \ \text{s}$$

As we could see from the plot the initial part of the simulation results that orange curve is still equal to zero is around 1 s, which confirms our above estimation for the time that is needed to carry the information from inlet to outlet. In fact, because of artificial compressibility of the LBM framework, there is a 1.20 time lag between inlet and outlet. I know, by decreasing the Mach number, I could decrease the time lag but this makes my simulations pretty expensive for my actual application. Now, if I look at the end of the signal in the above plot and try to find one cycle of my simulation result and compare it with the one cycle of the analytical solution, I have this plot:

In fact, I just matched the position of the last peak in my full signal visualization and then plot just the last cycle in the above chart. As you see, if we put the time lag aside, the amplitude and shape of the signal extracted from LBM matches the analytical solution quite well (analytical solution: blue curve and numerical results: orange dashed curve).

So finally my question: Am I allowed to match the last the peak of these two signals and claim that: if we put the time lag aside, which is purely an artifact of the LBM, the two signals match quite well and my framework works fine?


问:
Optimize speed for calculating the approximate inverse of a large matrix
I am searching for a faster method to calculate an approximate inverse of a large matrix (up to 32000x32000) resulting from a discrete non-linear system of partial differential equations. I'm using C++ together with the Eigen library (MKL activated). Unfortunately the resulting matrix is not sparse as it can be seen in:

The blue entries denotes negative and the orange positive numbers.
Currently I'm using a block matrix inversion scheme and then a partialPivLu() decomposition from Eigen to calculate the necessary inverses. The matrix has no specific properties like definiteness. Is there the possibility of a faster algorithm for a matrix with such a structure than to calculate the inverse or rather to solve the linear system of equations $Ax=b$ with the standard implementation of a Housholder or Lu decomposition?


问:
Resampling of values between body fitted and cartesian grids
Assuming there exist two 2D grids on the same domain, one of them is a cartesian one while the other is curvilinear (body fitted regular grid, with curved coordinate lines).

I am looking for a way to transfer values (interpolate) between these grids. Assuming the body fitted grid is given in the form of coordinates of the cell centers in the cartesian system doing the interpolation from the cartesian to the body fitted grid is relatively easy, as simple bicubic etc. interpolation can be used on the cell center coordinates to find the corresponding cartesian coordinates and get the values from the grid with these.

On the other hand for the other direction, body fitted to cartesian, this cannot be achieved so easily, as it would essentially require the inversion of the body fitted to cartesian diffeomorphism (mapping).

I was wondering if someone may have some suggestions for this problem. The numerical inversion of the mapping can be carried out, however I am a little worried about the compound effect of the inaccuracies during this process with all the interpolation involved.

ContentLicense: CC BY-SA 3.0


问:
Hardware supporting floats with fraction beyond 64 bit
Is there any computation accelerator (like a GPGPU) available, that natively (this means in hardware, not emulated by a library) supports arithmetics using floating point numbers with a fractional part longer than 64 bit (like e.g. IEEE 754 quadruple-precision)?</p>


问:
How to curve-fit the lower envelope of random sequence?
I'm more or less familiar with procedures and methods to fit a curve to experimental data, and I have done this many times using Matlab. However this time I have a problem that I'm not sure how to solve and I'm stuck.

I have some experimental data which I have plotted in the figure below (blue curve).

The data seem to be lowerbounded by some kind of smooth logarithmic-like curve like the red line shown in the figure (I have plotted this red curve manually using Paint just to illustrate what I mean).

I have several mathematical expressions that could be good candidates to fit the red line, but I'm not sure how to fit them in the way shown in the figure. Standard data fitting procedures use the whole data set and provide a curve somewhere in the middle of the blue curve, but I need something like the red curve.

I work with Matlab and I have tried to use the envelope function (lower envelope), the absolute value of the Hilbert transform and the convex hull (convhull) but nothing works.

Any ideas or suggestions'


问:
How can I improve the accuracy of the calculation of the magnetic field in Gmsh/GetDP?
I need to calculate the magnetic field along a straight line in proximity of an array of 6 magnets. I used the tutorial files "magnets" included in Gmsh and I slightly modified the file in order to describe the position of the magnets in my desired configuration.

I set the mesh size with a rather small value:

Mesh size at infinity: 0.01;
Mesh size on magnets:   0.001

and I use the "H" formulation. However, the resulting calculated magnetic field is very "fuzzy" as you can see in the figure below. Different colors are different components of the magnetic field. This fuzziness is clearly not physical and has a numeric origin. Am I missing something?

TL/DR:
How can I improve the accuracy of the numerical calculation and avoid numerical noise when calculating the magnetic field of an array of magnets?


问:
updating the matrix Adjugate/Cofactor
I would like to calculate the Adjugate matrix of a given matrix $A$, and its updates in the diagonal:
$B=A-\lambda I$, where $I$ is the identity matrix, $\lambda$ is a scalar. To this end, I am using the algorithm explained here, using a decomposition $A=XDY$. Unfortunately, the matrix $A$ or $B$ could be singular.

Is it possible to obtain a suitable decomposition such, that for the decomposition of $B$ I could save computer time by reusing $X$, $D$, or $Y$ from the decomposition of $A$?

Title: updating the matrix Adjugate/Cofactor


问:
How to make a directed graph symmetric?
Say I have a directed graph given as an adjacency matrix $A$ in CSR format represented by the arrays ia (row indexes) and ja (column indexes).
In my application the graph serves as the underlying spatial connectivity for an RBF-FD discretization of a PDE. The graph originates from finding the $N$-nearest neighbor stencil for each node $\boldsymbol{x}_i$ of a $d$-dimensional point cloud.

For some purposes it is benefical to use symmetric stencils $S$, meaning that $\boldsymbol{x}_j \in S(\boldsymbol{x}_i)$ implies $\boldsymbol{x}_i \in S(\boldsymbol{x}_j)$ for all $i,j$ (or at least for interior nodes). This corresponds to making the directed graph symmetric (undirected). One example of where this is needed is in the METIS partitioning library, which expects graphs to be symmetric. Similarly, the Reverse Cuthill McKee ordering algorithm also expects a symmetric matrix.

Question: Given an adjacency matrix in CSR format as arrays ia and ja, how can I find the symmetric graph adjacency matrix arrays ias, jas?

I have noted in Scipy, the way they achieve this is by forming the matrix $A + A^T$ (see Scipy source here). Note that I don't want to compute the actual matrix values, but only the structure of the resulting matrix.

Is there any other way to achieve this, apart from the naive way of of initializing $M$ empty adjacency lists (for all $M$ nodes of my point cloud), iterating through all the nodes, and pushing each connection $(i,j)$ into the right list? It seems kind of unneccesary to process each connection twice.

Title: How to make a directed graph symmetric?


问:
Does incomplete LU preconditioning improve the asymptotic scaling of Krylov subspace methods?
It is well known that unpreconditioned Krylov subspace methods applied to the finite-difference-discretised Poisson equation with $n$ grid points per direction require $O(n \, |\log(\varepsilon)|)$ iterations to reduce the error in the initial guess by a factor $\varepsilon \in (0,1)$. This is typically a prohibitively large number of iterations; hence it is often suggested to reduce the number of iterations through incomplete LU (ILU) preconditioning.

My question is: does ILU reduce the $\text{# iterations} = O(n \,|\log(\varepsilon)|)$  estimate asymptotically (e.g. to $\text{# iterations} = O(\sqrt{n} \,|\log(\varepsilon)|)$), or does it merely reduce the prefactor implied by the big-O notation?

My own experiments seem to suggest that ILU preconditioning only improves the prefactor, but maybe I am not using the method properly?


问:
performance comparison between PETSc and SLATE
We want to start a new project to solve a large-scale inverse problem (O(10^6) number of parameters) to invert for subsurface wave speeds. We will use FEM to solve forward and adjoint PDEs. In our project, we will need to have both sparse and dense solvers. Previously in similar projects, we were using PETSc library but today when I was searching for new libraries, SLATE library came to my attention (link). I was just wondering if there is any comparison between these two libraries from the performance point of view. Perhaps to make this question richer, I know Trilinos is another competitive library, but based on my personal experience and research that I did, I found PETSc is faster while Trilinos is more reliable. Please comment on this if you have any objections.</p>


问:
How to derive the adjoint sensitivity equations for a least squares objective function gradient
The Problem

I would like to determine the gradient of a least  squares objective function which depends on a vector of 40 parameters $p$, and the solution of a system of 32 differential equations. In specific, let the aforementioned system be:

$$\dot{y} = f(t,y(t,p),p) \tag{1}$$

then the least squares objective function is

$$G(t,p) = \frac{1}{2}\sum^n_{i=1}||y(t_i,p) - \hat{y}(t_i)||^2\tag{2}$$

where $\hat{y}$ is some measured data, and $||\cdot||$ is the L-2 norm.

My Attempt

For a small number of parameters I can use the typical approach employing the sensitivity equations which yields the following for the gradient of the objective function:

$$\frac{d}{dp}G(t,p) = \sum^n_{i=1}\left(y(t_i,p) - \hat{y}(t_i)\right)\frac{\partial y(t_i,p)}{\partial p} \tag{3}$$

where $\frac{\partial y(t,p)}{\partial p}$ is determined by solving the sensitivity equations simultaneously with $(1)$, which is one extra system of 32 differential equations for every value of $p$. For a small number of parameters $p$ the equation systems are kept reasonably small. For my case however with 40 parameters, each with a system of 32 equations to solve, there is a total of 1312 to solve which makes it very demanding on the solver.

I think that the solution for when the sensitivities of many parameters need to be calculated is to use the adjoint sensitivity equations where according to this website, they are as follows:

For an objective function $$G(t,p) = \int^T_0 g(t,p)dt \tag{5}$$ we have to solve a single extra differential equation for $\lambda^*$, namely: $$\frac{d\lambda^*}{dt} = \frac{\partial g}{\partial y} - \lambda^* J_y, \ \lambda^* = 0 \tag{6}$$ Then the gradient of the objective function is $$\frac{d}{dp}G(t,p) = -\int_{0}^{T}\lambda^{\star}\frac{\partial f}{\partial p}dt-\lambda^{T}(0)\frac{\partial y}{\partial p} \tag{7}$$

The Question

What is the function $g$ for my case where I have the least squares objective function?

How do I determine the value of $\frac{\partial y}{\partial p}$? Isn't the point of this approach to avoid the calculation of this value so that I don't have to use the sensitivity equations?

For the case of the least squares objective function in eqn. $(2)$, how is the gradient formulated using the adjoint method shown in equation $(7)$?

EDIT: After much searching I found the answer to the first question in this paper. Because equation $(5)$ is inheritly continuous, the least squares formulation has to be translated to such a format using the Diract Delta function $\delta(t)$: $$G(t,p) = \int^T_0g(t,p)dt=\int^T_0\sum^n_{i=1}\sum^n_{j=1}\frac{(y_j(t_i,p)-\hat{y}_j(t_i))^2}{2}\delta(t_i-t)dt \tag{8}$$ As a result equation $(6)$ becomes a discontinuous differential equation: $$\frac{d\lambda^*}{dt}=\sum^n_{i=1}\left(y(t_i,p) - \hat{y}(t_i)\right)^2\delta(t-t_i) -\lambda^*J_y \tag{9}$$


问:
Finding block structure of a tensor
Are there any well-known algorithms for partitioning a dense tensor into block-sparse form?

In other words, I need to find a set of non-overlapping blocks that contain all non-zero entries of the tensor, while trying to minimize both the number of blocks and the the number of zeroes in each block.

The tensor in question is 3D if that helps.

Title: Finding block structure of a tensor


问:
Why is perfect sampling not used in large-scale lattice model simulations?
The statistical physics literature is replete with papers describing simulations of lattice models, such as the Ising model. Typically, these are done through Monte Carlo methods, such as the Metropolis algorithm. A major concern in these studies is whether the algorithm has reached equilibrium. If the dynamics are not run for a long enough time, they will not produce draws from the desired equilibrium distribution. Further, in many cases, the algorithm is never actually at equilibrium, but rather approximates it to a finer and finer resolution as the simulation proceeds. Hence, the distance to equilibrium is one key source of error in these simulation methods.

There is, however, a Monte Carlo method for sampling exactly from the equilibrium distribution. In general it is extremely computationally expensive, and therefore inferior to standard methods, but there are many cases where the specific structure of the problem can be leveraged to make it practical. The Ising model is one such case, and the algorithm proceeds as follows. Start one chain in the state with all spins up, and one chain in the state with all spins down. Couple the chains by choosing the same proposal spin for each step of the dynamics for each chain. Run the Metropolis dynamics until the chains coalesce at some time $T$. The value of the chain at the random time $T$ will be distributed exactly as the equilibrium distribution.

This is known as coupling from the past and is originally due to Propp and Wilson, in the paper

Exact Sampling with Coupled Markov Chains and Applications to Statistical Mechanics, published 1996. Further, it is known that the coalescence time is on average approximately the same as the mixing time for the chain (in the Ising case).

Naively, I would think perfect sampling in this way might be superior to traditional algorithms, such as Metropolis. At the cost of running two chains for the mixing time, I get a perfect i.i.d. sample from the target distribution. In comparison, I need to run the chain for the mixing time in the Metropolis algorithm to get what is merely an approximation to an i.i.d. sample (with each roughly independent sample requiring the mixing time again). Thus, for the cost of running twice as many chains, I am able to remove to the equilibration error completely. It's not obvious to me this is better, given a fixed budget of Metropolis runs (trading eliminating the error for halving the number of samples), but it seems like it might be. At the very least, I hoped to find some discussion of this method in the statistical physics simulation literature.

However, I did not. I looked in what I was told was the standard reference, A Guide to Monte Carlo Simulations by Binder and Landau, and could not find any mention of coupling from the past or other perfect sampling methods. Further, among the citations to the Propp and Wilson paper, I could not find any statistical physics papers.

I wonder then, why is the physics community not interested in perfect sampling? For most problems, there is not any monotonicity structure as in the Ising model, so it is clearly inferior to usual methods. But, there are also many cases where this kind of structure can be exploited. My guess is that the method is known to researchers but not used because it is inferior to the usual methods in some clear way, which I am not able to discern.

So, to conclude: is perfect sampling not used by physicists because there it is obviously worse than other Monte Carlo methods? If so, why is it worse?

(If the answer is model-dependent, let us restrict ourselves to the Ising and related – e.g. Edwards–Anderson – models, since those are the ones I am interested in at the moment.)


问:
Probability Density Function of DNS velocity field
I am currently working on a DNS turbulent solver and I would like to compare my IHT simulations to papers. Those papers show the Probability Density Function of the velocity field:

I would like to reproduce this figure from data. I have access to all the velocity fields at any point during the simulation. I don't understand how to compute this pdf from a 3D velocity field or how to take into account all my data.


问:
What's the best way to implement a least-squares estimation of a motor system in MATLAB?
Basically, I'm trying to use Least-Squares to estimate the parameters of a DC motor.

My system can be modeled by the following matrix equation:

$$\begin{bmatrix}V_{input}(t)\\0\end{bmatrix}=\begin{bmatrix}\frac{\text{d}i}{\text{d}t}(t) && I(t) && \omega(t) && 0 && 0 \\ 0 && 0 && -I(t) && \frac{\text{d}\omega}{\text{d}t}(t) && \omega(t) \end{bmatrix} \times \begin{bmatrix}L\\R\\k_t\\J\\B\end{bmatrix}$$

From what I understand (admittedly probably not enough), I can estimate what the constants are for the vector $\begin{bmatrix}L,R,k_t,J,B\end{bmatrix}^T$ if I happen to know the values of the time-dependant variables for a handful of data points. If I take the equation to be of the form $Y=M\times K$, then the least-squares estimate at a particular discrete moment should be given by $\hat K=(M^TM)^{-1}M^TY$.

I have a motor set up to be driven by a sinusoidal voltage signal, offset so that the current doesn't reverse (for simplicity). I've measured the angular rate, the applied voltage, and the applied current. To get the necessary derivatives from the noisy signals, I built a sinusoid-fitting script in MATLAB which minimizes the average absolute difference between the fit signal and the noisy signal across all the measured data points. I'm no statistician, but the resulting fit eyeballed as "okay" so I moved on.

Larger version

Using a smooth fit curve for each of my measured parameters, I calculated the necessary derivatives and plugged everything into MATLAB. (I'm aware that K = inv(M.'*M)*(M.'*Y) isn't the most computationally-friendly, but I need to stick to the experimental procedure I listed in my write-up.)

Larger version

My output for $\begin{bmatrix}L,R,k_t,J,B\end{bmatrix}^T$ is absolute garbage though, and none of the values MATLAB spits out seem plausible. There are absurdly large values, very small values, and lots of negative values. (Physical properties shouldn't be negative.)

Larger version

I had planned to estimate $\begin{bmatrix}L,R,k_t,J,B\end{bmatrix}^T$ for each of my data vectors and then average those results, but I clearly won't get anything meaningful when my results are already useless.

Would anyone here happen to have an idea where I might have gone wrong?

I've stared at my code until I went cross-eyed without making too much progress, so I suspect that it might be a case of "garbage-in leads to garbage-out." My slightly-educated-guesses for what's happening are:

I've overlooked something fundamental in the math.

There's a data mismatch somewhere.

I'm going to be lost within the floating-point errors. (Symbolic MATLAB says the determinant for $M^TM \equiv 0$, but I know people use this kind of math fairly frequently. Therefore, I'm left to conclude that the calculation works but is ridiculously sensitive to changes in the measured values or something.)

Does anyone have any thoughts? I'm honestly not sure what I'm missing. (This is based upon a modification from this paper.)

Note: this question was reposted to SciComp.SE upon the suggestion that Physic.SE might not be as well-suited for it.

Edit:
Using the pseudo-inverse and letting MATLAB do the $(M^TM)^{-1)M^T$ part more efficiently has helped, but my output still looks like garbage.

Larger version


问:
Given an unpivoted form of Aasen's algorithm, how does one add pivoting?
I've implemented the version of Aasen's algorithm described in the book Matrix Computations 4th Edition. The version there doesn't have pivoting. The book's description of how to add pivoting is a bit sketchy for me. Can someone illustrate with actual code how to add pivoting?

I should add that my implementation has support for Hermitian matrices, unlike the one in the book.

from sympy import *

def aasen(A):
    """Aasen without pivoting, but with support for Hermitian matrices."""
    n = A.rows
    alpha = zeros(n,1)
    beta = zeros(n,1)
    L = eye(n)
    h = zeros(n,1)
    v = zeros(n,1)
    for j in range(n):
        if j == 0:
            alpha[0] = A[0,0]
            v[1:n,0] = A[1:n,0]
        else:
            h[0] = beta[0].conjugate() * L[j,1].conjugate()
            for k in range(1,j):
                h[k] = beta[k-1] * L[j, k-1].conjugate() + alpha[k] * L[j,k].conjugate() + beta[k].conjugate() * L[j,k+1].conjugate()
            h[j] = A[j,j] - L[j,0:j].dot(h[0:j])
            alpha[j] = h[j] - beta[j-1] * L[j,j-1].conjugate()
            v[j+1:n,0] = A[j+1:n,j] - L[j+1:n,0:j+1] * Matrix(h[0:j+1])
        if j <= n-2:
            beta[j] = v[j+1]
        if j <= n - 3:
            L[j+2:n,j+1] = Matrix(v[j+2:n]) / v[j+1]
    T = zeros(n,n)
    for i in range(n):
        for j in range(n):
            if i == j+1:
                T[i,j] = beta[j]
            elif i == j:
                T[i,j] = alpha[j]
            elif j == i+1:
                T[i,j] = beta[i].conjugate()
    return L, T
```


问:
Minimax optimization with an oracle
I have an optimization problem of the following form:  $$\min_y\left[\max_x f(x,y)\right].$$  It is fairly straightforward to minimize $f(x,y)$ over $y$ with $x$ fixed, and similarly to maximize $f(x,y)$ over $x$ with $y$ fixed.  If it helps, my $x$ problem is concave---but while I can solve the $y$ problem in closed-form it is not convex.

An obvious candidate for an optimization algorithm for the minimax problem might look something like the following:
$$
\begin{array}{rl}
x_{k+1}&\gets\arg\max_x f(x,y_k)\\
y_{k+1}&\gets\arg\min_y f(x_{k+1},y)
\end{array}
$$
But, this algorithm doesn't have to converge.

Given oracles for optimizing in one variable at a time, can we construct an optimization algorithm for the minimax problem?

The gradient descent-ascent (GDA) algorithm has a similar flavor, but for our problem finding the global optima for $x,y$ independently is possible so it doesn't seem necessary to just take small gradient steps.  If it helps, I'm happy to add a quadratic term to $f$ and assume I can optimize the objective $f(x,y)+\alpha\|y-y_0\|_2^2-\beta\|x-x_0\|_2^2$.

Title: Minimax optimization with an oracle


问:
Problem with recursive implementation of Subspace Iteration method in Numpy
I am having trouble with implementing the method of subspace iteration to find the eigenvalues and vectors of a random, symmetric matrix, A that is mxm with m = 10. The function that I have written to perform this task is

```<code>def subspaceiterate(A,V,v,j):
    if j == 0:
        return v
    else:
        v_jm1 = V[:,j-1]
        v_jm1 = np.reshape(v_jm1,(np.size(V,axis=0),1))
        v = v - np.matmul(v_jm1.T,np.matmul(A,v_jm1))
        j = j - 1
        return subspaceiterate(A,V,v,j)
</code>```

where A is my symmetric data matrix, V is an mxm matrix whose columns represent the current, calculated state of the eigenvectors, v is a particular mx1 column ofV that represents a particular eigenvector and j is the position of vin V.

As I understand it, the basic idea of Subspace Iteration is to simultaneously perform a Power Iteration on the different v's while subtracting away the components of v that are in the direction of the eigenvectors associated with eigenvalues larger than the one associated with the particular v in question. I have tried to design this function to operate in such a capacity but the results I get do not support the aforementioned idea. All the eigenvalues converge quickly (as they should), but do so to a value very close to the largest one instead of to the different, distinct eigenvalues (whose values I know because I have calculated them beforehand with np.linalg.eigs()).

The v's are normalized in the program before and after the function is called and the eigenvalues calculated with the Rayleigh-Taylor quotient v_j.T*A*v_j (implemented with np.matmul()). At this point I am lost and hoping someone brilliant can shine some light and help me out.  `


问:
Help with vectorized computation where vector multiplication would need to create an additional dimension
I want to use vectorized computation (with python and numpy) to do some calculation involving bessel functions.

As is often the case, I need to calculate the sums of bessel functions of (infinit) differnt order. I need to do that for a whole vector field, and would like to use vectorized computation. Actually, both, the bessel orders and the dimensions of the field would be handy to do vectorized, but they intersect. Here is a minimal, dummed down demo of the issue:

```<code>import numpy as np
from scipy.special import kn, iv

order = 20
size = 200

r = np.arange(-size, size, 1.0)
theta = np.arange(0, 2 * np.pi, 2 * np.pi / size)
z = np.arange(0, 1, 1 / size)

R, T, Z = np.meshgrid(r, theta, z)

m = np.arange(order)

eta = m * R
bessel_term = iv(m, eta) * kn(m, eta)

print(sum(bessel_term))
</code>```

What would need to happen here is, that of each value of R the bessel functions with the order-vector would need to be calculated. It would need to create an additional dimension in the vector, so to speak.

Instead I get the following error message which I understand, but want to avoid:

```<code>Traceback (most recent call last):
    eta = m * R
ValueError: operands could not be broadcast together with shapes (20,) (200,400,200)
</code>```

How can I do that? One way would be to iterate through the space dimensions manually (not vectorized) and use the vectorisation for the inner loop where the sums over the bessel functions are calculated. But to me it would be preferable if there was a way to indicate to numpy that at this point an additional dimention in the involved vectors is needed and collapse this dimension again when the infinite sum is calculated.

Background, the math:

In reality I want to calulate a magnetic field with three vector components in every point in space of my meshgrid. In every point I need to calculate an infinit sum of besselfunctions of infinit orders:

$$ B_r = j\frac{3 \mu_0 I_1}{4 \pi r} \sum_{m= ... -5, -2, 1, 4, 7...}^{\infty} 2 m I_m(\eta_m) K_m\left( \eta_m \frac{r}{a’} \right) + \frac{2 \pi r m q}{p} \left[ I_{m-1}(\eta_m) K_{m-1}\left(\eta_m \frac{r}{a’} \right) + I_{m+1}(\eta_m) K_{m+1}\left( \eta_m \frac{r}{a'} \right) \right] \exp[jm(\theta - 2 \pi z/p)]$$

with $r$, $z$ and $\theta$ as input parameters (space coordinates), $a',p, q$ as constant arguments and $\eta_m=|mq|$ as an abreviation.


问:
Upper bound on condition number in linear preconditioning
I'm studying iterative methods for solving linear system, and I find the following setting in Wikipedia:

Consider a matrix splitting $A = M-N$, where $A,M,N$ are all symmetric and positive definite matrices. Define iteration matrix $C=I-M^{-1}A$. We aim to show that the condition number $\kappa(I-C)=\kappa(M^{-1}A)$ is not too large, which ensures that iterative method converges fairly fast. Wikipedia gives the following bound on $\kappa(M^{-1}A)$:
$$\kappa(M^{-1}A) \leq \frac{1+\rho(C)}{1-\rho(C)}$$
where $\rho(C)$ denotes the spectral radius of matrix $C$. Here we are taking the condition number with respect to the Euclidean 2-norm. That is to say, we have $\kappa(M^{-1}A) = \frac{\sigma_{max}(M^{-1}A)}{\sigma_{min}(M^{-1}A)}$, where $\sigma_{max},\sigma_{min}$ denote the largest/smallest singular value of any matrix, respectively.

Any ideas on proving this upper bound? Any help/hint would be greatly appreciated!

(Link for reference: https://en.wikipedia.org/wiki/Preconditioner)


问:
Kronecker-factored least-squares?
Suppose $a_i,b_i$ are $d\times 1$ matrices, $y_i$ is scalar, and I need to find least-squares solution $w$ of the following system of $n$ equations:

$$y_i=(a_i^T\otimes b_i^T)w$$

Is there a specialized approach for this system?
Converting to normal equations and applying Cholesky to solve for $w$ takes takes $O(d^6)$ operations. This is a quite a bit expensive than verifying the solution at $O(d^2)$ operations. I was hoping for an approach that scales as $O(d^4)$ or $O(d^3)$. For my problem, $d\le n<d^2$

Title: Kronecker-factored least-squares?


问:
How applying the gradient descent method for solving a least square problem can remove the blur from an image?
I got an assignment where it asked to implement (in MATLAB) the gradient descent algorithm in order to resolve an ill posed least square problem:

$$
\min_u \Vert Gu - f \Vert
$$

where $u$ is the reconstructed image, $G$ is the blur applied to $u$ and $f$ is the blurred image.

I have successfully implemented it with MATLAB, it is effectively removing a portion of the applied blur, but I do not understand why applying a gradient descent method to this ill posed least square problem can effectively remove the blur from an image.


问:
Unsteady Stokes equations in ALE framework
I'm trying to solve Unsteady Stokes equations on a moving domain, using an ALE formulation, that is
$$\frac{\partial \mathbf{u}}{\partial t} - \mathbf{w}\cdot \nabla\mathbf{u} = \nu\Delta\mathbf{u} - \nabla p$$

where $\mathbf{w}$ is the velocity of the mesh.
I'm starting from the simple case of a rigid sphere translating with constant velocity, $\mathbf{W}$, in an unbounded viscous fluid. What I'm testing is that, if I solve the problem in a fixed control volume (the standard way, let's say), simply having an imposed velocity on the sphere in the middle of a large domain, I have that the force acting on it is given by the Stokes formula, that is $\mathbf{F} = 6\pi R \nu \mathbf{W}$ (normalized by density), once a stationary is reached.

Now, if I solve for the same thing on a control volume that moves following the sphere, I should have the same result, that is the same force acting on the sphere. I just plug the constant velocity $\mathbf{W}$ in place of $\mathbf{w}$ in the equation above:

$$\frac{\partial \mathbf{u}}{\partial t} - \mathbf{W}\cdot \nabla\mathbf{u} = \nu\Delta\mathbf{u} - \nabla p$$

and use the same boundary conditions: $\mathbf{u} = \mathbf{W}$ on the surface of the sphere, and $\mathbf{u} = \mathbf{0}$ on a boundary sufficiently far, or zero stress, should make no difference I think.

I solve the Unsteady Stokes equations problem with P2/P1 finite elements and a standard splitting technique of the Chorin kind: for reference I use the pressure correction projection method reported by Guermond and Shen

Guermond, J. L.; Minev, P.; Shen, Jie, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Eng. 195, No. 44-47, 6011-6045 (2006). ZBL1122.76072,

with BDF1 method for time stepping.

The problem is that the force results to be wrong: while in the fixed domain case it approaches the theoretical result, even though with some approximation, in the ALE case the error is much larger, of the order of 35%.

Changing the mesh or the time step seems to have no effect whatsoever on the final value, which makes me think that there must be something wrong in the setup, but I don't see what since I'm simply moving the points of the mesh, translating with a constant velocity (so in theory I don't even have to actually move them), and this shouldn't modify the physics.

Am I wrong? Are there any special precautions to be taken when the problem is that of Unsteady Stokes??


问:
Devising Convergent Numerical Scheme for PDE
I'm currently looking at the PDE
\begin{align*}
 u_t + \left[x(1-y) - (1-x)\right]u_x - (1-y) u_y + (z-xy) u_z  = (z-xy) u_{xy} - (1-x)u& \\
\end{align*}
with
\begin{align*}
u(x,y,z,0) = 1& \\
(x,y,z) \in [0,1]\times [0,1] \times [0,1]&,
\end{align*}
which I found via the conversion of a chemical master equation. The solution $u$ of this PDE is the probability generating function of the master equation. Is there any way to solve this numerically? It resembles advection-diffusion with a source, but it's not parabolic due to the mixed derivative term and variable coefficients. I don't believe the lack of boundary conditions is an issue, since if we just look at the advective terms, characteristics flow out of the boundary. Any help/suggestions would be greatly appreciated!


问:
A Question About a Claim from 1991 Computational EM paper about the Cancellation of certain Boundary Terms
Please let me know if this is not the appropriate site for this question. I found questions regarding EFIE/MFIE/CFIE on this site, so I thought my question might fit.

I am studying the paper by Putnam & Medgyesi on the CFIE method for BoR targets with junctions. In Section V (second paragraph) of the paper, the authors claim that by equating coefficients associated with half-triangle functions at the junction leads to a cancellation of boundary terms that one would get from integration-by-parts in the computation of the EFIE $L$ operator (see eqn. (5)) when the testing surface is open.

My issue is that I can't seem to reproduce this claim. I outline my approach below. Let us take eqn. (25) from the paper with $\beta = 0$. This formula is for the geometry in Fig (2). If I understand the discussion in Section V correctly, the boundary terms that arise in obtaining formulas for $\mathcal{L} (S_1^c , S_1^c ; R_1)$ and $\mathcal{L} (S_1^c , S_1^d ; R_1)$ should cancel if the expansion coefficients in the currents associated with the half-triangle functions are equated. After discretizing the generating curve of the BoR into $N_1 - 1 > 1$ straight line segments described by the interval $[t_{1k}^- , t_{1 , k+1}^-]$, $k \in \{ 1 , \cdots , N_1 - 1 \}$, (the subscript $1$ here indicates we are on the "outer" surface of the BoR) and parametrizing the surface by arclength along the generating curve $t$ and the polar angle $\varphi$, we expand the currents $\mathbf{J}^c_1$ according to

$$ \mathbf{J}_1^c \simeq \sum_{n \in \mathbb{Z}} \sum_{k = 1}^{N_1^c} \left( c_{1nk}^t \mathbf{J}_{1nk}^t + c_{1nk}^{\varphi} \mathbf{J}_{1nk}^{\varphi} \right) , $$

and likewise for $\mathbf{J}^d_1$. Above, $\mathbf{J}_{1nk}^t (t , \varphi)$ and $\mathbf{J}_{1nk}^{\varphi} (t , \varphi)$ are $\mathbf{u}_t$-directed and $\mathbf{u}_{\varphi}$-directed expansion functions which are both harmonic in $\varphi$ and use triangle functions in the $t$ variable, and $N_1^c$ is the number of segments making up the conducting part of the BoR outer surface (likewise, $N_1^d$). Above, $\mathbf{u}_t$ is the unit tangent vector along the BoR generating curve and $\mathbf{u}_{\varphi}$ is the usual polar unit vector. Equating the coefficients associated with half-triangle functions at the junctions means to me that $c_{1 n , N_1^c}^t = d_{1 n , 1}^t$ and $c_{1 n , N_1^c}^{\varphi} = d_{1 n , 1}^{\varphi}$. (This assumes the generating curve for the conducting part of the surface ends at the junction, and the dielectric part begins at the junction.)

The general formula for a matrix element is

$$ \mathcal{L}^{\alpha \beta}_{\ell k} (S_1^c, S_1^d ; R) = i k_0 \eta_0 \int_{S_1^c} ( \mathbf{J}_{1 n \ell}^{\alpha} )^* (\mathbf{r}) \cdot \int_{S_1^d} \{ \mathbf{J}_{1nk}^{\beta} (\mathbf{r}') G (\mathbf{r} - \mathbf{r}') + k_0^{-2} (\nabla G) (\mathbf{r} - \mathbf{r}') \mathrm{div}_S~{\mathbf{J}_{1nk}^{\beta}} (\mathbf{r}') \} d S' d S , $$

where $G (\mathbf{r}) = e^{- i k_0 |\mathbf{r}|} / (4 \pi |\mathbf{r}|)$, $\mathrm{div}_S$ is the surface divergence, and $\alpha , \beta \in \{ t , \varphi \}$. In the second term in the expression above, we integrate-by-parts to put the gradient on $G$ onto the testing function $(\mathbf{J}_{i n \ell}^{\alpha} )^* (\mathbf{r})$. Since $S_1^c$ is open, this will produce the boundary term

$$ \oint_{\partial S_1^c} [ ( \mathbf{J}_{1 n \ell}^{\alpha} )^* \cdot \mathbf{u}_t ] (\mathbf{r}) \left( \int_{S_1^d} G (\mathbf{r} - \mathbf{r}') \mathrm{div}_S~{\mathbf{J}_{1nk}^{\beta}} (\mathbf{r}') d S' \right) d s .$$

This boundary term is nonzero only when $\alpha = t$ and $\ell = N_1^c$. Based on the above formula, my understanding is that the authors are claiming the term

$$ \sum_{\beta \in \{t, \varphi\}} \oint_{\partial S_1^c} ( J_{1 n , N_1^c}^t )^* (\mathbf{r}) \left\lbrace \sum_{k = 1}^{N_1^c} c^{\beta}_{1nk} \left( \int_{S_1^c} G (\mathbf{r} - \mathbf{r}') \mathrm{div}_S~{\mathbf{J}_{1nk}^{\beta} (\mathbf{r}')} d S' \right) + \sum_{k = 1}^{N_1^d} d^{\beta}_{1nk} \left( \int_{S_1^d} G (\mathbf{r} - \mathbf{r}') \mathrm{div}_S~{\mathbf{J}_{1nk}^{\beta} (\mathbf{r}')} d S' \right) \right\rbrace d s $$

that would appear in equation (25) of the paper vanishes because $c_{1 n , N_1^c}^t = d_{1 n , 1}^t$ and $c_{1 n , N_1^c}^{\varphi} = d_{1 n , 1}^{\varphi}$. I am somewhat struggling to see this. The authors claim this is true because effectively $S_1^c \cup S_1^d$ form a closed surface. This suggests to me to try integrating-by-parts on the term

$$ \int_{S_1^c} G (\mathbf{r} - \mathbf{r}') \mathrm{div}_S~{\mathbf{J}_{1nk}^{\beta} (\mathbf{r}')} d S' = \oint_{\partial S_1^c} G (\mathbf{r} - \mathbf{r}') [ \mathbf{J}_{1 n k}^{\beta} \cdot \mathbf{u}_t ] (\mathbf{r}') d S' - \int_{S_1^c}  \nabla' G (\mathbf{r} - \mathbf{r}') \cdot \mathbf{J}_{1nk}^{\beta} (\mathbf{r}') d S' . $$

I can see how the boundary term for $\beta = t$ above would cancel out with a similar term in the similar formula for $S_1^d$ because $c_{1 n , N_1^c}^t = d_{1 n , 1}^t$. However, it is the lagging term

$$ \int_{S_1^c}  \nabla' G (\mathbf{r} - \mathbf{r}') \cdot \mathbf{J}_{1nk}^{\beta} (\mathbf{r}') d S' $$

that doesn't seem to disappear. This is making me quite confused about the authors' claim.

My question: Why do these boundary terms cancel? What am I missing? Is there anyone on this site who has worked on Method of Moments CEM codes or the EFIE and can shed light on the junction problem in general?


问:
"black box" preconditioner for shifted linear systems?
Does anyone know of any strategies for creating a preconditioner $P^{-1}_\sigma \approx (A+\sigma I)^{-1}$ given a preconditioner $P^{-1} \approx A^{-1}$, preferably such that the precomputation doesn't depend on the value of $\sigma$?


问:
Solving a linear system with block-banded matrix
It is required to solve a linear system $Ax=b$, where the matrix $A$ is symmetric, all the variables and coefficients are real.
The structure of $A$ is
\begin{equation}
A = \begin{pmatrix} A_{11} & A_{12} \\ A_{12}^T & A_{22} \end{pmatrix},
\end{equation}
where $A_{11}$, $A_{12}$ and $A_{22}$ are banded square blocks of dimension $N$, which may be large. $A_{12}$ is not symmetric.

I would like to avoid solving the full-matrix system, and to exploit the banded structure of the submatrices, without inverting any of the blocks.


问:
Book recommendation on numerical methods for solving Integro-Differential equations
I was wondering if anyone could recommend a good book or resource on numerical methods for solving integro-differential equations? Of course I am familiar with the methods for solving ODEs and PDEs numerically, so no problem there. But I don't have any frame of reference on the stability of numerical schemes when solving integro-differential equations, such as the Smoluchowski coagulation equation, or Volterra integro-differential equations. I can try to drop in some numerical integration routine into the ODE solver, but I am just not sure under what conditions this gives me a good solution versus accumulating errors over time.

I looked at one book by Lakshmikantham called Theory of Integro-Differential Equations, but I did not see any mention of numerical schemes. I will continue to look over that book, but if anyone knows of a good reference that would be very helpful.


问:
Can automatic differentiation be used on the parameters of an optimization problem?
If I wanted to perform an optimization using a Newton-based solver where the Hessian and gradient of a function are known analytically, and then use a package such as Adept to compute a Jacobian matrix ("reverse mode") with respect to the optimization parameters, would that be possible?  Would there be excessive memory usage if the optimization takes several thousand steps to converge?</p>


问:
Conservation of energy test for 2-body problem
I'm trying to implement a C++ code for the evaluation of the solution of an N-body system of ODE. I've started with a 2-body problem just to set the methods properly. I'm comparing the results of a 4 stages Runge-Kutta (RK4), a 4 steps Adams-Bashforth (AB4) and velocity Verlet method. The solutions I get have the correct behaviour for every methods I used, but once I check the conservation of energy for the Verlet method I obtain the following

As far as I understood from the article by Hairer, Lubich and Wanner, Stormer-Verlet method doesn't preserve the energy completely, and this should also apply to velocity Verlet. But I also know that Verlet should be particularly suited to this kind of problems, though I get these behaviours for the RK4 and AB4 energy conservation

which seems to me better than Verlet, they doesn't constantly grow with time at least. So is there something I'm missing about how Verlet method works or could it be just some code problem?

EDIT 1

I have some update and I answer here to Wolfgang Bangerth too.
Every picture above results from a time step of $h = 0.1$, using $1000000$ iterations.
I realized there was a problem in the code by the way. Now the energy for Verlet method is

which is the same as the previous one but it oscillates more. This oscillation reduce by decreasing the timestep (I tried with time step $h = 0.1, 0.05, 0.025$, related to a fixed number of iterations $N = 1000000$ and about $23, 12, 6$ orbits around the central body respectively) as well as with the code I used before and the slope is almost the same in both cases (about $3 \times 10^{-12}$). Both this picture and the previous one result from $h = 0.1$ for instance.

EDIT 2

Ok, now Verlet seems working properly. Even if the oscillation is clearly not centered on the initial energy and there's also a little modulation which seems to slowly grow with time. I also checked the angular momentum, which should be conserved even better than the energy by Verlet, and I think I can consider it as such (it actually oscillates a bit with a negative slope of order $10^{-21}$).


问:
Problem with solving coupled ODE and DAE equations with mass matrix (Error using daeic12 (line 77) This DAE appears to be of index greater than 1)
I am trying to solve 6 ODE equations coupled with 1 DAE one. The ODE equations have been discritized in space domain and ode15s MATLAB solver is used to solve the equations in time domain. I have defined a diagonal mass matrix including unity elements for ODEs and zero values for DAE equation. Assuming that "n" is the No. of nodes along the space, M has been defined as follows:

```<code>  O=ones(1,6*n);
  Ze=zeros(1,n);
  dia=horzcat(O,Ze);
  M=diag(dia);
</code>```

Here is the important parts of my main and function m.files:

```<code>.
.
.

y0(1:n)=zeros(n,1); %C0
y0(n+1:2*n)=zeros(n,1); %qc0
y0(2*n+1:3*n)=Tini*ones(n,1); %T0
y0(3*n+1:4*n)=Tini*ones(n,1); %Ts0
y0(4*n+1:5*n)=zeros(n,1); %N0
y0(5*n+1:6*n)=zeros(n,1); %qN0
y0(6*n+1:7*n)=zeros(n,1); %U0
t_int=1;
tf=3000;
nt=(tf/t_int);
tspan=(1:t_int:tf);
options=odeset('Mass',M);
[t,y]=ode15s(@adsor,tspan,y0,options);
</code>```

and

```<code>
function dydt=adsor(t,y)
.
.
.
dydt=zeros(7*n,1);
y(1,1)=Cin(1); %BC of C
y(n,1)=y(n-1,1); %BC of C
y(2*n+1,1)=Tg; %BC of Tg
y(3*n,1)=y(3*n-1,1); %BC of Tg
y(4*n+1,1)=Cin(2); %BC of N
y(5*n,1)=y(5*n-1,1); %BC of N
y(6*n+1,1)=Q/A; %BC of U
C(1:n,1)=y(1:n,1);
qc(1:n,1)=y(n+1:2*n,1);
T(1:n,1)=y(2*n+1:3*n,1);
Ts(1:n,1)=y(3*n+1:4*n,1);
N(1:n,1)=y(4*n+1:5*n,1);
qN(1:n,1)=y(5*n+1:6*n,1);
U(1:n,1)=y(6*n+1:7*n,1);

for i=2:n-1

    dqcdt(i,1)=K1*(qsc-qc(i));
    dCdt(i,1)=aEC(i).*C(i+1)+aWC(i).*C(i-1)-aPC(i).*C(i)+dqcdt(i-1));
    dqNdt(i,1)=K2*(qcn-qN(i));
    dNdt(i,1)=aEC(i).*N(i+1)+aWC(i).*N(i-1)-aPC(i).*N(i))+dqNdt(i-1);
    dTdt(i,1)=aET(i).*T(i+1)+aWT(i).*T(i-1)-aPT(i)+*T(i))+Ts(i);
    dTsdt(i,1)=T(i)-Ts(i)+dqcdt(i)+dqNdt(i);
    dUdt(i,1)=U(i)-U(i-1)+dqcdt(i)+dqNdt(i);

end
dydt(1:7*n,1)=[dCdt;0 ;dqcdt;0 ;dTdt;0 ;dTsdt;0 ;dNdt;0 ;dqNdt;0 ;dUdt;0];
</code>```

The equation located at the last line (i.e. dUdt(i,1)=...) is the DAE equation which is recognized by the defined mass matrix (please see "M" parameter defined before). I am not sure if this method is correct to specify the DAE equations! Anyway,when I run these codes, I face the following error:

Error using daeic12 (line 77)
This DAE appears to be of index greater than 1.

I checked the Jacobian and found that it is a n*n matrix at t0. There are two columns (the first and last ones) which are all zero and I think it makes the Jacobian matrix singular (as explained in https://ch.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html).

Could you please guide me how to overcome this error and fix my code? Any help is appreciated.

Best regards


问:
Continuous limit and nonlinear functional analysis
I have a kind of general question about approximation schemes in nonlinear functional analysis. Given a nonlinear map $\Phi$ from an open set (in an infinite dimensional Banach space) of functions to itself, I try to find fixed points of this map $\Phi(f)=f$, or equivalently to solve the roots of of $\Phi'=\Phi-{\bf 1}$ (where $\bf 1$ is identity functional ${\bf 1}(f)=f$) via some Newton method. (in my case $\Phi(f)(x) \sim \int f^n(y) K(x,y) dy$ for some integer $n>1$ and $K$ some integration kernel).

I am really a newbie in this and I have tried several types of discretization schemes, using polynomial interpolants and discretizing the functions $f$ on a grid and using piecewise-linear interpolants.

I am particularly interested in the theory that would tell me why or why not would the second option (with piecewise linear interpolants) converge well or not to the solutions of the continuous equation when I let the size of the grid tend to infinity, i.e. the spacing between two elements in the grid go to zero.

I empirically observe some erratic behaviour, and, while I have some physics intuition of why such limits can go bad sometimes, I would like to do some actual reading on this rather than doing the analysis myself like an amateur at this point, as I think this is some rather standard theory belonging maybe to Approximation theory or something. Google didn't quite help me on this.

(I originally asked on mathematics.stackexchange but discovered this scicomp and figured it might be a better place to ask).


问:
Open-source Hexahedral meshing tools
I'm looking for an open-source automatic meshing tool or library able to generate full or dominant hexahedral mesh starting from a provided arbitrary model i.e. a surface triangular mesh stored in STL file.

What I have found is snappyHexMesh embedded in OpenFOAM. Are there some other tools that fit my purpose?

Someone already asked a quite similar question here but it's a quite old answer. Has something new been released in the meanwhile?

Thank you

Title: Open-source Hexahedral meshing tools


问:
Methodology to develop parallel GPU simulations
I'm quite new to GPU or in general parallel programming. I am working on a Lattice-Boltzmann simulation that should be parallelized on GPU's. I wanted to asked more experienced people how to attack such a project. Would it be more beneficial to first implement all necessary features (basic algorithm with some extensions) and than to tackle the parallelization or should it be better done part by part. Speaking first the basic algorithm, than parallelizing this part. And then the same procedure for the extension?

Thank you in advance!


问:
Storing and retrieving two-dimensional and three-dimensional data
I work on computational geometry.
A huge number of two-dimensional and three-dimensional data are found in my project. Coordinates of polygon and polyhedrons vertices consisted of two-dimensional and three-dimensional data, respectively. Indeed, polyhedrons face data are also founded. Storing and retrieving them is nessecesery in some part of my code.

I want to know is there any efficient storing algorithm for my data which:

Require the least time and volume

Shapes can be easily retrieved

Please explain it in the simplest because I am new to this kind of algorithm.

I will appreciate any comments that give me a little help and forgive me for writing shortcomings.


问:
How to limit the output value of an Linear Congruential Generator, but still maintain an extremely long period?
Let us suppose that I have a full cycle LCG, following the standard equation

Hence, as it has a full cycle, it will have a period of m, and will output values ranging from 0 to  m. However, I want to limit the output values to somewhere between 0 to some arbitrary constant a, where a is less than m. All the while maintaining its original period length.

I don't have experience in this field, but I was thinking that wherever I want to use the output value, I just take mod a of the number, and be done with it? will this ruin the randomness or is a valid solution?

How do I go about doing this?


问:
Largest triangle that contains a point
Given the location of $n$ points on a 2D plane ($P_1, P_2, \ldots, P_n$); and the location of a special point $X$.

Find three points $P_i,P_j,P_k$ ($i \neq j \neq k$) such that point $X$ is inside the triangle $P_iP_jP_k$, and the minimum distance between point $X$ and the three sides of the triangle is the largest.

In other words, let $D_{xij}$ be closest distance between point $X$ and line $P_iP_j$, similarly for $D_{xjk}$ and $D_{xki}$. Find $i, j, k$ that maximizes $\min (D_{xij},D_{xjk},D_{xki})$.

The Naive solution runs in $\mathcal{O}(n^3)$. Goal: an algorithm that is better than $\mathcal{O}(n^3)$.

Title: Largest triangle that contains a point


问:
N-body correct scaling
I realized an usual way to scale an N-body problem for an N-body simulation is by choosing units such that gravitational constant $G = 1$, but I'm probably doing it the wrong way. Suppose I simply have a planet moving around the Sun (which I consider as the center of my reference frame). For example I choose

\begin{equation}
M_s = M_{sun} + M_1
\tag{1}
\end{equation}

\begin{equation}
D_s = D_1
\tag{2}
\end{equation}
as my mass and distance scaling units. Than I should choose $T_s$ in such a way that my new $G$, which I call $G'$, is equal $1$. So
\begin{equation}
G' = 1 = \frac{1}{G} \frac{D_s^{3}}{M_sT_s^2}
\tag{3}
\end{equation}
brings me to
\begin{equation}
T_s = \sqrt{\frac{D_s^3}{GM_s}}
\tag{4}
\end{equation}

After that, I divide masses and distances by $M_s$ and $D_s$, and velocities by $D_s/T_s$. This is what I understood from what I found about this scaling. Am i doing something wrong? I even tried different way to set $D_s$ and $M_s$, different starting units, but I always get wrong outputs.

Edit 1

This is the function I'm using (C++), with x as the vector containing velocity and position of the body moving aroun the Sun and y position vector of the Sun (considering the above example).

```<code>d = sqrt(pow((x[0] - y[0]),2) + pow((x[2] - y[1]),2) + pow((x[4] - y[2]),2));

//x system
p[0] = x[1];
p[1] = -G*m*(x[0] - y[0])/pow(d,3);
//y system
p[2] = x[3];
p[3] = -G*m*(x[2] - y[1])/pow(d,3);
//z system
p[4] = x[5];
p[5] = -G*m*(x[4] - y[2])/pow(d,3);
</code>```

The equations are specific for the choosen planet, so there's only the Sun mass on the right hand side. The $G$ is equal $1$ in this case.
The original dimensions are astronomic units for distances, kilograms for masses and days for times, so $[G] = \frac{AU^3}{kg\,days^2}$.

Title: N-body correct scaling


问:
Interpolation of a Concave Mesh
I would like to have an algorithm that interpolates the values attached to nodes in a concave mesh mesh.

To be me precise, assume we have a point cloud P (e.g. in 3 dimensions) and a list of edges E connecting it with form a concave mesh in the form $E=[[i_1,i_2],[j_,j_2],...]$, where $i,j$ are nodes.
Further assume that each point in P is associated with a value V (e.g. Temperature).

The algorithm I am looking for takes these objects together with an arbitrary point x in space and returns an interpolated value iff the point lies within the concave mesh and something unique if it doesn't (e.g. NaN). Hence if has a signature like:

```<code>f(P,E,V,x) -> float
</code>```

To illustrate this a bit. The point cloud might look like this with the hull of the edges shown in gray.

The algorithms I found mostly deal with convex point clouds exclusively.
In that case the Delaunay triangulation is unique (except in some edge cases) and the interpolation can be done e.g. with a linear barycentric interpolation.
However, I haven't found anything general about concave meshes. Yet I assume that this is something that has some general use as well and there must exist something already that I am not aware of.

Title: Interpolation of a Concave Mesh


问:
Parametric nonlinear programming
I believe, I have a parametric nonlinear optimization problem.
The non-convex constraints depend on some parameters, and I seek a solution that satisfies these constraints for all parameters in a specified continuous range.

Let $ {\bf A}(\vec{x},\vec{\theta}) := {\bf A} \in\mathbb{R}^{n\times n}$ be a matrix-valued function of $\vec{x}\in \mathbb{R}^n$ and $\vec{\theta}\in \mathbb{R}^m$ (here $m\ll n$).  Define
$ {\bf L}(\vec{x}):= {\bf L} \in\mathbb{R}^{n\times n}$ to be a matrix which is a function of $\vec{x}$, and $ \vec{F}(\vec{x}):= \vec{F} $ is a vector-valued function of $\vec{x}$.

Main problem:  The goal is to find any $\vec{x}$ such that
$$
{\bf A}(\vec{x},\vec{\theta})\vec{y} = \vec{F}(\vec{x}),
\\
{\bf L}(\vec{x})\vec{y}\ge \vec{0},
$$
holds for all $\vec{\theta}$ in some parameter space (ideally, $0<\vec{\theta}<\infty$ or $0<a\le \vec{\theta}\le b$ for scalars $a,b$).

For a fixed $\vec{\theta}$ I am able to solve the problem using nonlinear programming.  I have tried discretizing the parameter space and then solving the resulting coupled nonlinear programs, but this only leads to a solution that is valid for the discrete parameter values.

Is there a better way to think about this problem? Or any standard ways to solve it?  There is no objective function in the formulation I posed above, and I don't know if that is of any significance.

Title: Parametric nonlinear programming


问:
3D Matrix (Tensor) Operating On a 1D Vector
Say you have a tensor $T$ and the components are represented by a 3 by 3 by 3 matrix. And you want to use that tensor to map a vector $u$ into a new vector $s$, both of which are 3 by 1 column vectors.

How do I know how to convert the 3 by 3 by 3 matrix into a matrix that is compatible with multiplication with my vector? How do I know how to reshape the matrix? Which dimension to collapse?

I know that if I can collapse it somehow into a 3 by 3 by 1 matrix, multiplication with a 3 by 1 vector would yield another 3 by 1 vector (the mapped vector $s$ I think). But I have no intuition on how to perform that step.


问:
Integrating a wavelike equation with absorbing boundary conditions
I am trying to numerically solve the following equation:

$\frac{\partial^{2} \phi}{\partial t^{2}}-\frac{\partial^{2} \phi}{\partial x^{2}}+V(x) \phi(x, t)=0$ On some domain, with:
$\phi(x, 0) = I(x)$

$\frac{\partial \phi(x, 0)}{\partial t} = 0$

And with non-reflecting (absorbing) boundary conditions.

I am using finite difference method, so I get the following:

$\phi_{i}^{j+1}=-\phi_{i}^{j-1}+2 \phi_{i}^{j}+r^2\left(\phi_{i+1}^{j}-2 \phi_{i}^{j}+\phi_{i-1}^{j}\right) - (\Delta t)^2V_{i}\phi_{i}^{j}\,\,\, (1)$, where $j$ is used for time and $i$ for position, $r = \frac{\Delta t}{\Delta x}$ and $\Delta t$, $\Delta x$ are the time and space steps.

For $j = 0$, I use the $\frac{\partial \phi(x, 0)}{\partial t} = 0$ condition and the central difference approximation, $\left(\frac{\partial \phi}{\partial t}\right)_{j} = \frac{\phi_{j+1}-\phi_{j-1}}{2 \Delta x}$, to get:

$\phi_{i}^{1}= \phi_{i}^{0}+\frac{r^2}{2}\left(\phi_{i+1}^{0}-2 \phi_{i}^{0}+\phi_{i-1}^{0}\right) - \frac{(\Delta t)^2}{2}V_{i}\phi_{i}^{0} \,\,\, (1')$

To accommodate the boundary conditions I use the advection equation $\frac{\partial \phi}{\partial t}=- \frac{\partial \phi}{\partial x}$, again I use the central difference approximation, and I substitute into (1) to get:

$(1+r) \phi_{i, j+1}=2 r^{2} \phi_{i-1, j}+2\left(1-r^{2}\right) \phi_{i, j}-(1-r) \phi_{i, j-1}\,\,\, (2)$.

For $j = 0$ I use the same idea as before and I get:

$\phi_{i}^{1} = (1 - r^2) \phi_{i}^{0}+r^2\phi_{i-1}^{0}-\frac{(\Delta t)^2}{2}V_{i}\phi_{i}^{0}\,\,\, (2')$.

Now I am trying to implement this using python, and I use the following logic:

For the interior points I use equations (1), and (1'). For the boundary points I use equations (2), and (2'). But the results that I get are those for reflecting boundary conditions not absorbing ones.

I am using python, and if necessary I can post the code that I have written.

Any help or advice is very much appreciated. I also would appreciate any resources or references. Thanks


问:
"This DAE appears to be of index greater than 1" daeic12 (line76) error code
Hi I am trying to solve a set of pde converted into ODE and DAE using central finite difference method. I have used the MATLAB 'solve' command to determine the coefficients of fictitious nodes for boundary conditions. I am using the ode15s solver and getting "This DAE appears to be of index greater than 1" daeic12 (line76) error code. I have attached the code as a link below. The script file to run is main_code_live, which is a live script and explains how the pde have been converted into ode and dae along with the boundary conditions, rest all are necessary function files to run the main script. I have tried using a lot of suggestions found on Mathworks website but still not successful. Any help would be appreciated! https://drive.google.com/drive/folders/1emgOcdsnP81MP5ZBERefEJ0d4IP0qbSL?usp=sharing</p>


问:
How to convert this code to scan random datas instead of binned datas?
I am new to Fortran and I am trying to write a code using random data instead of binned data in x, y, z as shown in my sample code.

```<code>implicit real*8(a-h,o-z)
dimension rm(4),rp1(4),rip1(4),rp2(4),rip2(4),rp3(4),rip3(4),
d   rn(4),u1(4),u2(4),u3(4)

do ix= 1000,25000,1000
    x = ix/1000000.
    do iy = 1000,25000,1000
        y = iy/100000.
        do iz= 1,1000,25
            z = iz/10000.
            a = (x**2+y**2)/z
            b = x*y*z
            c = x*y**2+y*z**2+z*x**2
            fr = (a*b)/c
            if(fr.ge.0.05.and.fr.le.23)then
                write(40,*) x, y, x, fr
            else
            endif
        end do
    end do
end do

stop
</code>```

How to convert such code having binned data to a code using random draws.