问:
Avoid matrix multiplication in algebraic multigrid method
Currently when I try to solve a linear algebra system of the form of $A x =b$ I use the algebraic multigrid method. The algebraic multigrid method uses a Galerkin product to form the coarse grid matrix ($A_c = RA_fP$). Where $R$ is the restriction operator, $P$ the prolongation operator, $A_c$ the coarse grid matrix and $A_f$ the fine grid matrix.

In the current state of my implementation I perform the matrix multiplication to form the coarse grid. This step is quite expensive and I was asking myself if there is a way to use the restriction and prolongation operators directly in combination with the fine grid matrix instead of calculating the coarse grid matrix.

The first problem I was facing is to form the strong connection matrix from the fine grid matrix $A_f$, the restriction operator $R$ and the prolongation operator $P$.

Is there any available algebraic multigrid package that is not calculating the coarse grid matrix but using the operators?


问:
Trying to understand splitting-based iterative method for 2D Laplace problem
I am trying to understand the theory behind a splitting based iterative method which uses the incomplete Cholesky factorization. Before giving the specific details, let me first give the problem background:

I have discretized a linear system $Au = f$ for the following PDE problem:

$$-\frac{\partial u}{\partial x^{2}} -\frac{\partial u}{\partial y^{2}} = f $$

If for convenience we eliminate the boundary conditions and have the usual lexicographic ordering while using finite differences, then $A$ will have the same shape as described on page 4 here. Note that it is symmetric and positive definite. We will assume that it is of size $n \times n$.

Since it is a large sparse matrix, direct solvers are not a good method. So, I want to employ splitting based iterative solvers, i.e. we let $A = M - N$ which defines the following recursion:

$$u_{k+1} = M^{-1}Nu_{k} + M^{-1}f$$
where the starting vector is the zero vector.

My first question: If we let $L$ denote the incomplete Cholesky factor with zero fill-in, then since $L \cdot L^{T}$ is an approximation of A, can we take $M = L \cdot L^{T}$ in the above method?

Second question: If the previous answer is yes, how can I theoretically prove that the above method converges? I know that I need to prove that the spectral radius(i.e. maximum absolute eigenvalue) $\rho (B_{M}) < 1$, where $B_{M} = I - M^{-1}A$ is the iteration matrix. Especially since it is not clear what form $M^{-1}$ takes here, I am unsure of how to proceed.

The third question is closely related to the previous one. Given that the stopping criterion is $\frac{\left\lVert f - Au_{k}\right\rVert _{2}}{\left\lVert f\right\rVert _{2}} \leq 10^{-10}$, where $f = -12x^{2}y^{5} - 20x^{4}y^{3}$, I am also interested in quantifying the required amount of iterations for convergence. All I know is that the smaller $\rho (B_{M})$ is the faster it converges, but how to translate this into a minimum number of iterations?

The final question has to do with the amount of computational work per iteration, if we look at the above recursion then it is clear that at the very least a linear system solve with $M$ is required, but it is unclear to me what amount of work this requires.

Edit: I have indeed looked into the alternative direct methods suggested in the comments, but for me it's important to gain a good grasp for relatively simpler methods like the one described here. Especially for comparative purposes I'd like to be able to answer the above questions, so that I can better consider its weaknesses and see what kind of alternative methods work best. I still do appreciate the suggestions though!

ContentLicense: CC BY-SA 4.0


问:
Can Representation Theory be studied computationally / numerically?
Can a subfield such as the representation theory of Lie algebras be studied computationally / numerically -- is there an interplay between the abstract and the concrete?  I would be grateful for an example, e.g. a reference to a journal paper.</p>


问:
Is it possble to do this complex symbolic calculation with Matlab?
Sorry it's bit abrupt, but recently I am caught up in some symbolic calcualtion which is tedious and almost impossible with mere human hands, so just wondering is it possible to solve the double integral of what's written below using Matlab or some other platform?

$$\int_{0}^{a}dk\frac{1}{k}\int_{|k-q|}^{|k+q|}dk'\frac{1}{k'}(k^2+q^2-k'^2)(k'^2+q^2-k^2)(k'^2-(k+q)^2)(k'^2-(k-q)^2)[\frac{1}{v_F(k-k')-w}+\frac{1}{v_F(k-k')+w}].$$

And to clarify, all the variables are positive real numbers. The first integral is about $k'$ with upper and lower limit being $|k+q|$ and |k-q| respectively, and second integral is about $k$ from $0$ to $a$. The symbols other than $k$ and $k'$ can all be seen as real positive constants.

OwnerDisplayName: Xinxin Peng


问:
Why the two Gram-Schmidt algorithms produce different results for qr factorization?
For the qr factorization using classic Gram-Schmidt algorithm, I found the 2 different implementations below. The first one uses the for loop to compute the upper triangular matrix R, but the second one uses the matrix-vector multiplication. Since they are mathematically equivalent, but they produce different results. I cannot find why it is so, can anyone give me some indications? Furthermore, I often find that whenever I understand an algorithm theoretically, if I implemented it, I will have some trouble because of different numerical results. This is why I always ask questions about how, actually, to implement the algorithm with specific code. Because even if I know the algorithm, I still cannot generate the correct results. Thanks.

%   test the classical Gram-Schmidt algorithm
clc;clear;format compact;
A = hilb(7);
%%  method 1
R(1,1) = norm(A(:,1));
Q(:,1) = A(:,1)/R(1,1);
m = size(A,2);
for j=2:m
    for i=1:j-1
        R(i,j) = Q(:,i)'*A(:,j);%   the difference between 2 methods from below
    end
    q_hat = A(:,j)-Q(:,1:j-1)*R(1:j-1,j);
    R(j,j) = norm(q_hat);
    Q(:,j) = q_hat/R(j,j);
end
%%  method 2
R1(1,1) = norm(A(:,1));
Q1(:,1)=A(:,1)/R1(1,1);
for j=2:m
    R1(1:j-1,j) = Q1(:,1:j-1)'*A(:,j);  %   the difference between 2 methods from above
    q_hat = A(:,j)-Q1(:,1:j-1)*R1(1:j-1,j);
    R1(j,j) = norm(q_hat);
    Q1(:,j) = q_hat/R(j,j);
end
%%  compare
norm(A-Q*R)
norm(Q'*Q-eye(m))
norm(A-Q1*R1)
norm(Q1'*Q1-eye(m))
```

<p>My results are as follows:</p>

<pre class="lang-matlab prettyprint-override"><code>ans =
     0
ans =
    0.3369
ans =
   8.9238e-10
ans =
    0.1890
</code>```


问:
Hit-n-Run Monte Carlo on convex polytope
So, I'm currently trying to implement a MCMC to uniformly sampling hyper-points from the polytope defined as $\mathbb{K}=\{x\in\mathbb{R}^{n}\;\;\text{s.t.}\;\; A\,x=b \}$ in the specific case where, given a generic linear transformation $A\in\mathbb{R}^{m\times n}$, $b\equiv0\in\mathbb{R}^m$ and the boundary conditions $0\leq x\leq1$ hold.

Now, although I was able to succesfully perform the simulation (I am using Julia), there are some things I am not quite sure about:

given the fact that polytopes of this this kind tend to have star-like shapes in higher dimensions, two pre-processing steps are required, before launching the simulation:

the first regards the so-called blocked-flux adjustment which consists, quote from the text of the exercise, find fluxes $i$ such that $\max_{x\in K} x_i = \min_{x\in K} x_i = z_i$ and remove such variables from the system, adjusting the vector $b$. Can anyone please explain to me what the heck this means?

the second consists in finding an optimal inner point in $\mathbb{K}$ as a starting point of the chain, which intuitevely enough has to be located far from the vertexes of the polytope. The text tells me this: it can be done e.g. by computing $\frac{1}{2n}\sum_{i=1}^n (x^{\min,i} +  x^{\max,i})$ where $x^{\min,i} \in \arg\min_{x\in K} x_i$ and $x^{\max,i} \in \arg\max_{x\in K} x_i$. Here I just do not understand the notation: I suppose I should compute a weighted average of the midpoints of each edge of the polytope but I can't see how this is related to the above formulation.

as any coherent MCMC, the walk in the state space has to be s.t. it satisfies the detailed balance, which, for the present case, the text tells me that it should be the target distribution $$p(x) \propto \delta^m(Sx-b)\prod_i \theta(u_i-x_i)\theta(x_i-l_i) $$
again, I have no idea how this is obtained nor how to compute them.


问:
Solving saddle point problem having non-invertible top-left block with a PETSc nested matrix
My system is a symmetric FE problem with lagrange multipliers:

$Z=\begin{pmatrix}A & C^T \\ C & 0\end{pmatrix}$

The matrix $A$ is positive semi-definite, non-invertible. The whole matrix is invertible.

I am working on the development of a finite element software where I want to solve such problems with PETSc in parallel. At this time I would like the $Z$ matrix to be a MATNEST in which the $A$ and $C$ matrices are assembled independently.

With the representative simple code at the end of this message (1D laplacian with Lagrange multipliers), and using a PCFIELDSPLIT preconditioner, the solver diverges. If I artificially force the matrix $A$ to be invertible (adding the command-line option -force_invertible), then I get the right solution.

So, what can I do to solve such a system with a nested matrix?

Another question: is it possible to use direct methods (Mumps) with nested matrices in parallel?

```<code>#include "petscsys.h" /* framework routines */
#include "petscvec.h" /* vectors */
#include "petscmat.h" /* matrices */
#include "petscksp.h"

#include <vector>
#include <string>
#include <iostream>
#include <numeric>

// Try to solve with Petsc a simple 1D laplacian problem (or a series of springs)
// o-////-o-////-o-////-o-////-o-////-o-////-o-////-o-////-o ...
//
// The boundary conditions (Dirichlet) are imposed using Lagrange multipliers.
// The system to be solved is of the form:
//
//      Z     x  =  y
//    --^--  -^-   -^-
//    [A Ct] [u] = [b]
//    [C 0 ] [l]   [d]
//
// A is positive semi-definite (1 eigenvalue is 0)
// Z is indefinite invertible

static char help[] = "Saddle point problem: scalar 1D laplacian with lagrange multipliers.\n\n";
int main(int argc,char **argv)
{

  // Initialization ------------------------------------------------ //
  MPI_Init(NULL, NULL);
  PetscErrorCode ierr;
  ierr = PetscInitialize(&argc,&argv,NULL,help);CHKERRQ(ierr);

  int rank, nproc;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  MPI_Comm_size(PETSC_COMM_WORLD, &nproc);

  PetscViewer viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD);
  PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE);

  PetscBool forceInvertible=PETSC_FALSE;
  PetscOptionsGetBool(NULL,NULL, "-force_invertible", &forceInvertible, NULL);

  // Problem definition and dof spliting among processes ----------- //
  double k=1.;

  std::vector<int> globalDofs={0, 1, 2, 3, 4, 5, 6, 7, 8};

  std::vector<std::pair<int, double>> globalDirichlet={{0, 0.},
                                                       {8, 1.}}; // first: global dof number, second: imposed value

  auto start=[&](int rank){
    int q=(globalDofs.size()-1)/nproc;
    int r=(globalDofs.size()-1)%nproc;
    return q*rank + std::min(rank, r);
  };

  std::vector<int> dofs;
  for (int i=start(rank); i<=start(rank+1); ++i) {
    dofs.push_back(globalDofs[i]);
  }

  std::vector<int> isLocal(dofs.size(), 1);
  if (rank!=0) isLocal[0]=0;
  int nLocalDofs=std::accumulate(isLocal.begin(), isLocal.end(), 0);

  std::vector<int> dirichletIdx;
  for (int j=0; j<globalDirichlet.size(); ++j) {
    for (int i=0; i<dofs.size(); ++i) {
      if (isLocal[i] && globalDirichlet[j].first==dofs[i]) {
        dirichletIdx.push_back(j);
      }
    }
  }

  for (int p=0; p<nproc; ++p) {
    MPI_Barrier(PETSC_COMM_WORLD);
    if (rank==p) {
      std::cout << "Rank: " << rank << std::endl;

      std::cout << "  dofs   : ";
      for (int d: dofs) {
        std::cout << d << " ";
      }
      std::cout << std::endl;

      std::cout << "  isLocal: ";
      for (int d: isLocal) {
        std::cout << d << " ";
      }
      std::cout << std::endl;

      std::cout << "  nLocalDofs: " << nLocalDofs;
      std::cout << std::endl;

      std::cout << "  dirichlet: ";
      for (int i: dirichletIdx) {
        std::cout << "{"
                  << globalDirichlet[i].first << ", "
                  << globalDirichlet[i].second
                  << "} ";
      }
      std::cout << std::endl;

    }
  }

  // Matrix A ------------------------------------------------------ //
  Mat A;
  MatCreate(PETSC_COMM_WORLD, &A);
  MatSetType(A, MATMPIAIJ);
  MatSetSizes(A, nLocalDofs, nLocalDofs, PETSC_DETERMINE, PETSC_DETERMINE);

  MatMPIAIJSetPreallocation(A, 5, NULL, 1, NULL);

  auto setValue=[&](int i, int j, double v){
    if (forceInvertible && (i==2 || j==2)) {
      MatSetValue(A, 2, 2, 1., ADD_VALUES);
    }
    else {
      MatSetValue(A, i, j, v, ADD_VALUES);
    }
  };

  for (int i=0; i<dofs.size()-1; ++i) {
    // k * [ -1  1 ]
    //     [  1 -1 ]
    setValue(dofs[i]  , dofs[i]  , -k);
    setValue(dofs[i+1], dofs[i+1], -k);
    setValue(dofs[i]  , dofs[i+1],  k);
    setValue(dofs[i+1], dofs[i]  ,  k);

  }

  MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

  MatView(A, viewer);

  // Matrix C ------------------------------------------------------ //
  Mat C;
  MatCreate(PETSC_COMM_WORLD, &C);
  MatSetType(C, MATMPIAIJ);

  MatSetSizes(C, dirichletIdx.size(), nLocalDofs, PETSC_DETERMINE, PETSC_DETERMINE);
  MatMPIAIJSetPreallocation(C, 1, NULL, 0, NULL);
  for (int i: dirichletIdx) {
    MatSetValue(C, i, globalDirichlet[i].first, 1., ADD_VALUES);
  }

  MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);

  // Zero Matrix --------------------------------------------------- //
  Mat N;
  MatCreate(PETSC_COMM_WORLD, &N);
  MatSetType(N, MATMPIAIJ);
  MatSetSizes(N, dirichletIdx.size(), dirichletIdx.size(), PETSC_DETERMINE, PETSC_DETERMINE);

  MatMPIAIJSetPreallocation(N, 0, NULL, 0, NULL);
  MatAssemblyBegin(N, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(N, MAT_FINAL_ASSEMBLY);

  // Matrix Z ------------------------------------------------------ //
  Mat Z;
  Mat subZ[4];

  subZ[0]=A;
  MatCreateTranspose(C, &subZ[1]);
  subZ[2]=C;
  subZ[3]=N;

  MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, subZ, &Z);
  MatView(Z, viewer);

  // Vector b ------------------------------------------------------ //
  Vec b;
  VecCreate(PETSC_COMM_WORLD, &b);
  VecSetType(b, VECMPI);
  VecSetSizes(b, nLocalDofs, PETSC_DECIDE);

  VecAssemblyBegin(b);
  VecAssemblyEnd(b);

  // Vector d ------------------------------------------------------ //
  Vec d;
  VecCreate(PETSC_COMM_WORLD, &d);
  VecSetType(d, VECMPI);
  VecSetSizes(d, dirichletIdx.size(), PETSC_DECIDE);

  for (int i: dirichletIdx) {
    VecSetValue(d, i, globalDirichlet[i].second, ADD_VALUES);
  }

  VecAssemblyBegin(d);
  VecAssemblyEnd(d);

  // Vector y ------------------------------------------------------ //
  Vec y;
  Vec suby[2];

  suby[0]=b;
  suby[1]=d;

  VecCreateNest(PETSC_COMM_WORLD, 2, NULL, suby, &y);

  VecView(y, viewer);

  // KSP ----------------------------------------------------------- //

  KSP ksp;
  KSPCreate(PETSC_COMM_WORLD,&ksp);
  KSPSetOperators(ksp, Z, Z);
  KSPSetFromOptions(ksp);

  PC pc;
  KSPGetPC(ksp,&pc);
  PCSetType(pc, PCFIELDSPLIT);
  PCFieldSplitSetDetectSaddlePoint(pc, PETSC_TRUE);
  IS isg[2];
  MatNestGetISs(Z, isg, NULL);
  PCFieldSplitSetIS(pc, "u", isg[0]);
  PCFieldSplitSetIS(pc, "l", isg[1]);

  // PC pc;
  // KSPSetType(ksp, KSPPREONLY);
  // KSPGetPC(ksp,&pc);
  // PCSetType(pc, PCLU);
  // PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);

  Vec x;
  ierr = VecDuplicate(y, &x);
  KSPSolve(ksp, y, x);

  // View ---------------------------------------------------------- //
  MatView(Z, viewer);
  VecView(y, viewer);
  KSPView(ksp, viewer);
  VecView(x, viewer);
  // --------------------------------------------------------------- //

  ierr = PetscFinalize();

  MPI_Finalize();
  return ierr;
}
</code>```


问:
Explanation of subspace strategy regarding CG described in Golub's book
I was wondering about the last paragraph in Matrix Computations (4th edition) by Golub, Chapter 11 (11.3.3), specifically his explanation of subspace strategy for Conjugate Gradient.

Note that in Algorithm 11.3.1, the step-k optimization is over the affine space $x_k + span\{\nabla\phi(x_k)\}$.

Given that $Ax_0 \approx b$, our plan is to produce a nested sequence of subspaces

$\space S1 \subset S_2 \subset S_3 \subset \cdots $

that satisfy $dim(S_k) =k$ and to solve the problem
  $\underset{x\in x_0 + S_k}{min} \phi(x)$
  each step along the way.

If $x_k$ is the step-k minimizer, then because of the nesting we have $\phi(x_1)\geq \phi(x_2) \geq \cdots \phi(x_n)=\phi(x_*)$. since $S_n = \mathbb{R}^n $, we ultimately obtain $x_* = A^{-1}b$. Even though this is a finite-step solution framework, it may not be attractive if $n$ is extremely large. The challenge is to find a subspace sequence that promotes rapid decrease in the value of $\phi$, for then we may be able to terminate the iteration long before $k$ equals $n$.

With this goal in mind we note that at $x_k$ the function $\phi$ decreases most rapidly in the direction of the negative gradient. Thus, it makes sense to choose $S_{k+1}$ so that it includes $x_k$ and the gradient $g_k = \nabla\phi(x_k) = Ax_k -b$. This strategy guarantees that $x_{k+1}$ is at least as good as a steepest descent update:

$\underset{x\in x_0 + S_{k+1}}{min} \phi(x_{k+1} \leq \underset{\mu \in \mathbb R}{min}\space \phi(x_k -\mu g_k )$

If $x_0$ is an initial guess and we define $g_0 = Ax_0 -b$, then since $\nabla \phi(x_k) \in span \{x_k,Ax_k\}$ it follows that the only way to satisfy this requirement is to set:

$S_k = \mathcal{K}(A,g_0,k) = span\{g_0, Ag_0, A^2g_0,\dots, A^{k-1}g_0\}$

Why does Golub write that $\nabla \phi(x_k)\in span(x_k, Ax_k)$?

I don't understand how this is possible since $\nabla\phi(x_k) = Ax_k -b$, how do you get $b$ from $span(x_k, Ax_k)$?

I think there's something trivial that I'm missing here...


问:
2-dimensional Gauss-Hermite quadrature in R
A similar question was asked here and the given answer is perfect for a unidimensional integration.

I need to make bidimensional integration in R with a Gauss-Hermite quadrature:
$$\int_{R^2} h(p1,p2) \phi(p1,p2) dp1 dp2$$
with $p1$ and $p2$ two parameters that follow a multivariate normal distribution of density $\phi(.)$ such as:
$$p1,p2 \sim \mathcal{N_2} (\mu, \Sigma) $$
with $\mu = \left(\matrix{0\\ 1} \right)$ and $\Sigma = \left[\matrix{0.6&0.5\\0.5&0.25} \right]$.

I tried some packages in R but the results I got were very different from one another. I also tried the solution here but I'm not sure of the results I got.

I am looking for a simple solution, preferably based on gauss.quad.

Not being an expert, I have no idea of when to divide by $\pi$, if I should make a rescaling in my case and so on.

Thanks in advance for any help !


问:
Correct approach for thermal finite element simulation of layered assembly
I would like to optimise the heat transfer on a PCB. Several dies are on the top and cooling air is going through the fins in heat sink on the bottom. The assembly consists of several layers like shown in the simplified figure below.

I would like to investigate the influence of for instance several thermal greases on the maximum temperature of the chips: I want to vary the thermal coefficient and the thickness of the thermal grease layer.

The “problem”/”issue” is that the dimensions in the Z direction (micrometer) are much smaller than X and Y directions (~10 mm). This makes the meshing process complicated using “naive” finite element simulations.

Is there a better approach to analyse this type of layered assembly using finite element simulation?


问:
How to minimize a integral function using a constant step gradient method in Python?
I am developing a practical work of the following system of ode

\begin{align}x'(t) &= k_1h(t) - (k_2+k_3)x(t)\\
y'(t) &= k_3x(t)\end{align}

and $z(t) = (1-k_4)(x(t)+y(t))+k_4h(t)$, where $h(t) = 2t+1$

We now consider the following cost function:

$$J(u) = \frac{1}{2}\int_{t_0}^{t_f}|z_u(t)-\bar{z}(t)|^2w(t)dt+\frac{\lambda}{2}|u-u_{ap}|^2.$$
Where $u=(k_1,k_2,k_3,k_4)$, the function $w(t)$ modulates the weight that is given to each value taken from the $\bar{z}$ data ($\bar{z}$ is an interpolated function of the data). If we consider that all the values have the same importance we can only take $w(t) = 1.$ The second term present in $J(u)$ is interesting when we have an idea of the values of the parameters to optimize. It therefore allows you to stay neighborhood to these values.

In order to be able to minimize $J(u)$ we introduce the Lagrangian:

$$L(y,u,p)=\frac{1}{2}\int_{t_0}^{t_f}|C(Y(t),u)-\bar{z}(t)|^2w(t)dt+\frac{\lambda}{2}|u-u_{ap}|^2+\int_{t_0}^{t_f}(Y'(t)-f(Y(t),u))p(t)dt,$$

where $Y(t)=(x(t),y(t))$ and $C(Y(t),u)=z_u(t)=(1-k_4)(x(t)+y(t))+k_4h(t).$

So
$$J(u) = L(Y_u,u,p).$$

When $D_uJ(u) = 0$ leads to the resolution of the following adjoint equation

\begin{align}-p'(t) &= (D_yf)^*p(t)-(C(Y(t),u)-\bar{z}(t))w(t)(D_yC)\\
p(T)&=0 \end{align}

I want to minimize $J(u)$ using a constant $\rho$-step gradient method:

$$u^{k+1}=u^k-\rho D_uJ(u^k).$$

So $D_uJ(u) = D_uL(y,u,p) = \lambda(u-u_{ap})\cdot\bar{u}+\int_0^T(C(Y(t),u)-\bar{z}(t)w(t))D_uC^*-D_uf^* p(t))\cdot \bar{u}.$

Where $\bar{u}$ is a vector of $\mathbb{R}^4$ and \begin{equation*}
D_yf =
\begin{pmatrix}
-(k_2+k_3) & 0 \\
k_3& 0
\end{pmatrix},
\end{equation*}
\begin{equation*}
D_uf =
\begin{pmatrix}
h(t) & -x(t) & -x(t) & 0 \\
0 & 0 & x(t) & 0
\end{pmatrix},
\end{equation*}
\begin{equation*}
D_yC =
\begin{pmatrix}
1-k_4 & 1-k_4
\end{pmatrix},
\end{equation*}
\begin{equation*}
D_uC =
\begin{pmatrix}
0 & 0 & -(x(t)+y(t))+h(t)
\end{pmatrix}
\end{equation*}

I have tried in several ways to "build" the integral function but as $C(Y(u),t)$ and $\bar{z}$ are interpolated functions and in addition the matrices $D_uJ(u),\dots$ shown also contain these interpolated functions, how can I implement?

Here my code:

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

# Creation of system of ode
def myfct_R(Y, t, params):
    k1, k2, k3, k4 = params
    x, y = Y
    h = 3*t+2 # pour tester la convergence des schémas
    return np.array([k1*h-(k2+k3)*x, k3*x])
# Resolution of system of ode
def z_fct(t,y0,params):
    k1, k2, k3, k4 = params
    sol_fct = odeint(myfct_R, y0, t,args=(params,))
    return (1-k4)*(sol_fct[:,0]+sol_fct[:,1])+k4*h(t), Y

params = [0.01,0.02,0.03,0.1]
Nt = 2**14 +1
y0 = [0.0,0.0]
t = np.linspace(0, tf, Nt) # tableau de t = 0->tmax avec N valeurs
zfunct,Z = z_fct(t,y0,params)

#---------------Adjoint equation-----------------

#------------****************------------

Cyu = interp1d(t,zfunct)
def Padj_fct(p, t,  params):
    p1, p2 = p # p(t) = (p1(t),p2(t))
    k1,k2,k3,k4 = params
    w = 2*t**2 -2*t -1
    Cyui = Cyu(t)#interp1d(t,zfunct)

    Dyf = np.array([[-(k2+k3),0],[k3,0]])
    Dyfa = np.transpose(Dyf)
    DyC = np.array([1-k4,1-k4])
    zb = h(t)
    H = Cyui - zb
    return Dyfa.dot(p)-(H*w)*np.transpose(DyC)

def sol_adj(t,y0,params):
    k1, k2, k3, k4 = params
    Y = odeint(Padj_fct, y0, t,args=(params,))

    return Y

padj = sol_adj(t,y0,params)
</code>```

And from here I don't know how to proceed with the calculation of the gradient method, I would also like to know if the Padj_fct function is good to numerically approximate the attached equation presented above.

Note: $DyC^* = DyC^T$ is the transpose matrix


问:
Conjugacy in Non-linear Conjugate Gradient Descent
In linear conjugate gradient method, our goal is to solve the system of linear equations $$Ax = b$$ where A is a symmetric positive definite matrix, and that is equivalent to finding the minimizer of the convex quadratic form : $$φ(x) = \frac12 x^T A x - x^T b$$

In the non linear CG method, our goal is finding the minimum of $f(x)$, a general continuous function for which the gradient $f'$ can be computed.

My question is, in linear CG, we say that the directions $p_{k}$ are conjugate with recpect to the matrix $A$ if $$p_{i}^TAp_{j} = 0$$ for all $i\neq j$, where $A$ is the matrix that appears in the system we would like to solve.

However, in nonlinear, such a matrix doesn't exist, does it? In that case, how do we define the conjugacy of our search directions or what do we mean by keep saying that the directions are conjugate?

Thanks in advance.


问:
Velocity self-advection stability problem with QUICK upwind scheme
I'm trying to implement a finite-difference fluid solver (incompressible, inviscous) using the QUICK scheme for advection. Mostly I've been following the Norris thesis chapter 2 (https://ses.library.usyd.edu.au/handle/2123/376) for the math.

It got advection of scalar fields to be stable and non-oscillating by using Adams-Bashforth (AB2) time integration and applying the ULTRA flux-limiter of Leonard.

The code breaks down, though, when applying this advection scheme to the velocity field itself. In particular at the domain boundaries (using open/zero-pressure boundary conditions) oscillations appear. In 1D the oscillation amplitude seems to stay limited, but in 2D diverging modes appear at the domain corners.

I tried a variety of methods of boundary cell handling without success.

Is it necessary to modify the QUICK/AB2 scheme as such for velocity self-advection?

I can post visualizations if helpful.


问:
Efficient inversion of multidimensional non linear function
I have a function $f:x\mapsto \vec{y}$ with $x \in [0,1]$ and $\vec{y}=(y_1,...,y_n) \in \mathbb{R}^n$. $n$ is a small integer e.g. $n=8$. Each of the component functions in $y_i(x)$ "oscilate" up and down (roughly between -1.5 and 1.5). A snippet of two such component functions looks like this,

I know that $f$ is invertible because it was designed to be bjective (using a numerical method). In other words, the vectors $\vec{y}$ are unique.

Edit:
Another way to view this goes as follows. $\vec{y}$ is a set of codewords which encode each position $x$. The generating function $f(x)$ is designed so that these code words are unique.

I am searching for an efficient numerical inversion method for $f$. I.e., I am looking for $f^{-1}: \vec{y} \mapsto x$.

So far, I have tried a neural network (using sklearn.neural_network.MLPRegressor) with one hidden layer with 200 neurons. This works fine, but I am looking for a solution for which the evaluation of $f^{-1}$ is computationally cheaper. A solution which gives me a small number of possible candidates would be acceptable, too. How can I solve this problem?


问:
Good 3D surface fits for multiscale oscillatory surfaces
I have a 3D surface in $x$, $y$, and $z$. where $z$ is a function of $x$ and $y$ and my points are on a structured grid in $x$ and $y$. My function $z$ is highly oscillatory and irregular with multiscale phenomena. I'd like a fitting tool that will capture the large scale phenomena and provide a differentiable closed function without getting derailed by the small scale phenomena. Are there any good 3D data fitting toolboxes or methods?


问:
Defining appropriate test function spaces for the finite element solution of Euler's fluid equations
I have the following coupled equations for the conservation of mass and momentum of a compressible fluid :

\begin{equation}
\rho_t + (\rho u)_z = 0,
\end{equation}

$$
(\rho u)_t + (\rho u^2)_z + \tfrac{1}{k}\rho_z = \left(\frac{dp}{dz}\right)_f + \rho g
$$
over the domain $[0, 2L]\times[0, \infty)$ for density and velocity $\rho = \rho(z, t)$, $u = u(z, t)$ respectively. Note that
$$
g = g(z) = \begin{cases}
9.81 & \text{if } 0 \le z \le L
\\
-9.81 & \text{if } L < z \le 2L
\end{cases}
$$
and that $\left(\frac{dp}{dz}\right)_f$ is a frictional pressure term that is non-linear in $\rho$ and $u$.

The weak variational problem is:

$\textit{find $u\in V_1, \rho\in V_2$ such that}$
  \begin{equation}
\int_0^{2L} v\rho_t - v_z\rho u dz = 0
\end{equation}
  $$
\int_0^{2L}q(\rho u)_t - q_z\rho u^2 - \tfrac{1}{k}q_z\rho - q\left(\frac{dp}{dz}\right)_f - q\rho gdz = 0
$$
  $\textit{for all test functions $v, q$ in spaces $V_1$, $V_2$ respectively.}$

I am then discretising these equations in time using an implicit backward Euler discretisation for implementation in the finite element solver software Firedrake (this is a standard approach when using Firedrake).

The problem I'm having is that the Firedrake solver does not converge with the inclusion of $\left(\frac{dp}{dz}\right)_f$, and gives the wrong solution when $\left(\frac{dp}{dz}\right)_f$ is excluded.

So in order to narrow down the causes of this divergence, my question is: are there specific spaces for $V_1$ and $V_2$ that should be used for this kind of variational problem?


问:
sequence\ flow analysis with score- beginner data scientist
I'm trying to create a model that will give me the best parts that lead to a maximized score out of a sequence.

my data(spark rdd) looks like this:

("dan", "john", "john", "guy", 45)

("john", "dan", "guy", 60)

("guy", "guy", "guy", "guy", "guy", "guy", "guy", "guy" 8)

("guy", "john", "john", "guy", "mark", 0)

the number of names is known, and finite.

the length of the longest sequence is also known and finite.

the sequence is by order of appearance, with the score being the outcome.

i also have the data of the specific time stamp when each occurrence happens.

I'm trying to find the best combination/ the best parts that will maximize my final score,(so i can say "john" is crucial to getting a high score, but 'guy" is not very effective)

any help will be extremely appreciated,
thanks in advance!


问:
Interior point of convex polytope
Suppose the convex polytope is the set of feasible solutions $\mathbf{x}\in\mathbb{R}^n$ for the linear system $\mathbf{A}\mathbf{x}=\mathbf{b}\,,\; \mathbf{A}\in\mathbb{R}^{m\times n}$ subject to some (linear) constraints $\mathbf{0}\leq\mathbf{x}\leq\mathbf{v}$; how is it possible to select an interior point of said polytope, i.e. it doesn't lie on an edge or vertex of the polytope?

Does the algorithm requires the knowledge of what the vertices of the polytope are?

Title: Interior point of convex polytope


问:
Solve system of PDES
I need to solve the following system of partial differential equations

\begin{align}
\psi_{x} &= -i a \psi - b \phi,\quad &(1)\\
\phi_{x} &= b \psi + i a \phi,\quad &(2)\\
\psi_{t} &= - ic \psi + d \phi, \quad &(3)\\
\phi_{t} &= - d \psi + i c \phi, \quad &(4)
\end{align}

where $\psi=\psi(x,t),\; \phi=\phi(x,t)$ and $\psi_{x} $ means derivative w.r.t $x$ and $a,\;b,\;c,\;d$ are constants. How to solve (1)-(4)?

Title: Solve system of PDES


问:
Polynomial rooting - fast root finding
I need to solve a rooting problem of a polynomial (the order of which is 2(N-1), where N=48).

So far I'm using Python Numpy algorithm (that relies on computing the eigenvalues of the companion matrix), but this is too slow solution for me. The polynomial can be written as:
$$\sum_{l=-N+1}^{N-1}D_{l}z^{l}$$
where $D_{l}$ is the sum of the elements along an $N \times N$ matrix:
$$D_{l}=\sum_{n-m=l}D_{mn}$$
In the attached image you can see all the roots (94 roots, half inside the unitary circle, half outside it); any algorithm that allows the calculation of M≪N closest(<1) to the unit circle roots?

I add here some further figures with hopes to improve the understanding of my question. I have done some tests plotting the roots locus (see pictures) concerning only some $D_{l}$ but this does not make so sense, as Federico Poloni says in the comments:

roots are a global thing that depends on the whole polynomial (or
  matrix)

However, it is interesting to note that some roots (those closer to zero) keep a fixed position on the plane (if we do a comparison between solving the all polynomial and reduced version of it). All this in view to better understand the whole polynomial (or matrix) structure and possibly to find a faster way to calculate only desired roots.

Title: Polynomial rooting - fast root finding


问:
Pressure boundary conditions in Stokes Equation in 2D
I am solving the steady-state incompressible Stokes equations in 2D:

\begin{equation}
\frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} = 0,
\end{equation}

\begin{equation}
\mu\left[\frac{\partial^2 u_x}{\partial x^2} + \frac{\partial^2 u_x}{\partial y^2}\right] - \frac{\partial p}{\partial x} = 0,
\end{equation}
and
\begin{equation}
\mu\left[\frac{\partial^2 u_y}{\partial x^2} + \frac{\partial^2 u_y}{\partial y^2}\right] - \frac{\partial p}{\partial y} = 0.
\end{equation}

I'm using finite differences with a staggered grid approach (velocities defined at the faces and pressure at the cell centers).

I'm struggling to find how to implement pressure boundary conditions in the left and bottom faces. My pressure stencil is an "L" stencil: I use $p_{i,j+1}$, $p_{ij}$, and $p_{i+1,j}$. Hence, at the left and bottom-left volumes, I don't use the pressure value in the "ghost cell". Any clue on how to implement pressure BCs in this setting? I would like, for instance, to validate my code with a parallel plate simulation (prescribed pressure at left and right faces, and no-slip at top and bottom faces).


问:
Solving a huge least squares system of equations when I can only evaluate Ax
I have a situation where I can generate a system of $M$ linear equations for $N$ variables ($N \ne M$). Implicitly this is of the form $Ax=b$ with $A \in \mathbb{R}^{M \times N}$, although I never explicitly construct $A$ (even as a sparse matrix), but I can evaluate $Ax$ implicitly for any vector $x \in \mathbb{R}^N$ and calculate $b$ (but this is quite a non-trivial process to do this). I wanted to use LSMR or LSQR (https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.lsqr.html#scipy.sparse.linalg.lsqr), however they require definining a linear operator that can evaluate both $Ax$ and $A^Ty$ for any $x \in \mathbb{R}^N$, $y \in \mathbb{R}^M$. Evaluating $A^T y$ in my case is not that easy I don't think, so are there any approaches I can use that don't require doing this? It feels like it should be possible as I am still fully evaluating all of the equations that I want to approximately hold, but perhaps it isn't!


问:
Proving convexity of Frobenius norm and correlation function formulations of an optimization problem
I have been working on formulating my requirements in the form of an optimization problem in a multi-output regression setting.
Firstly, I would like to make the variables I used in the problem and their dimensions clear.

There are $Q$ outputs, $N$ samples, and $P$ covariates (features).

$\hat{Y}\in\mathbb R^{Q\times N}$ is the prediction matrix,

$$\hat{Y}=WX$$

where $W\in\mathbb R^{Q\times P}$ is the parameter matrix. These are the parameters we estimate in the regression.

$X\in\mathbb R^{P\times N}$ is the data matrix

My requirement:

I want to preserve the correlation structure in multi-output regression, i.e., the correlation structure of the predictions and the correlation structure of the data should remain close. I had two formulations for it.
I am assuming that the $Y$'s are standardized, Hence $YY^{T}$ is equal to the correlation matrix of the data.

Formulation 1: Minimize the Frobenius norm of the difference between the correlation matrices of the data and the predictions

$$W^{*}=\underset{w}{\mathrm{argmin}}\frac{1}{\mathcal{N}}\|YY^{T}-\hat{Y}\hat{Y}^{T}\|_{F}^{2}$$

Formulation 2: Maximize the correlation between the correlation matrices of the data and the predictions

$$W^{*}=\underset{w}{\mathrm{argmax}} \operatorname{cor}\left(Y{Y}^{T}, \hat{Y}\hat{Y}^{T}\right)$$

where $\operatorname{cor}$ is the correlation function.

Now I have to prove the convexity of the above formulations:

My attempts: (I am giving a very brief description of what I have done and I can elaborate if it helps, kindly ask in the comments section)

I didn't try to generalize the proof for matrices of any size. I started with a matrix $W$ of very small dimensions ($2\times 1$) and $X$ of dimension($1\times 2$) and then tried to prove that the Hessian matrix for both above formulations is PSD (positive semi-definite). The result tends to depend on the values of the matrix $X$ and $Y$.

I am stuck with generalizing this proof. Kindly give your valuable suggestions to prove the convexity/nonconvexity of both the above formulations.

UPDATE: In Formulation 2 I was wondering if I could use the positive semi definiteness properties of the matrices $YY$$^{T}$,$XX$$^{T}$ so I tried the following:

Note: i am replacing the feature correlation matrix $XX$$^{T}$ with $B$ and the output correlation matrix
$YY$$^{T}$ with $C$
\begin{align}
W=\begin{bmatrix}
 \alpha_{11}&\alpha_{12}\\
 \alpha_{21}&\alpha_{22}\\
\end{bmatrix}
\end{align}

\begin{align}
B=\begin{bmatrix}
 b_{11}&b_{12}\\
 b_{21}&b_{22}\\
 \end{bmatrix}
\end{align}
\begin{align}
C=\begin{bmatrix}
 c_{11}&c_{12}\\
 c_{21}&c_{22}\\
 \end{bmatrix}
\end{align}
$\hat{Y}\hat{Y}^{T}=WXX^{T}W^{T}=WBW^{T}$

now i am rewriting $\operatorname{cor}\left(Y{Y}^{T}, \hat{Y}\hat{Y}^{T}\right)$ as $Vec(C)\cdot\!Vec(WBW^{T})$

now this is my new loss function :
$L(\alpha)=Vec(C)\cdot\!Vec(WBW^{T})$.(as $W$ is a function of $\alpha$ and . here is the dot product between the two vectors)

I have calculated the Hessian of the above loss function $L(\alpha)$ for dimensions of $W$($2\times2$),$C$($2\times2$),$B$($2\times2$) and ultimately the Hessian turns out to be $2*C\otimes B$ .Where $\otimes$ is the Kronecker product between two matrices. As both $B$ and $C$ are PSD matrices the Kronecker product of them should also be a PSD matrix (please find the link to the proof below) which means my loss function $L(\alpha)$ should be convex. But however, this may not hold for larger examples as I didn't generalize the proof.

Look at Corollary 2.4 for proof


问:
SHREC 2010 Descriptors
I will appreciate if I may find someone how can clarify for me the part regarding the quality of feature descriptor, shown in the figure below:

and this screenshot is from the article: SHREC

All my problem is that I'm not able to understand the situation clearly, appreciate simple words hoping to understand.
Further indications, what I'm trying to do is mentioned in this question, but I wasn't able to get a fair answer for it. : Evaluation

So, it's about having same descriptor but for two models, main and scaled(or rotated, translated...) and trying to check how much my descriptor is good and well.

Any other recommendations regarding such situation are much appreciated,I'm implementing my work in python.

Title: SHREC 2010 Descriptors


问:
Are there unproblematic max constraints when modelling problems as Linear Programs?
Suppose we have a linear objective function that we want to maximize.
All variables are from the set of reals.
We have a constraint of the form:

$$\max(x_1,x_2) + \max(x_3,x_4)\leq c\,, \text{ with } c\in\mathbb{R}$$

Can we transform this such that we obtain a linear program?
I have read that constraints of the form

$$\max(x_1,x_2)\geq x_3$$

are problematic because this transforms into

$x_1\geq x_3$ OR
$x_2\geq x_3$

Back to my constraint, I would do the following transformation:

$m_1 + m_2\leq c$

$\max(x_1,x_2)\leq m_1$

$\max(x_3,x_4)\leq m_2$

which is transformed into

$m_1 + m_2\leq c$

$x_1\leq m_1$

$x_2\leq m_1$

$x_3\leq m_2$

$x_4\leq m_2$

which looks plausible to me as there is no OR introduced and both have equivalent meanings to me. The blowup is just polynomial in the original size of the input.

Could you please say if this transformation is correct or not with a proof? I am very interested in understanding how such proofs work. I know that linear constraints split the convex feasible set into a convex feasible and unfeasible set.
The constraints that I introduced are linear which should make the feasible set convex.
But how does this relate to my original problem which had a maximization?
Will the solution of the obtained linear program be the optimal solution of the original optimization problem?


问:
Normalized Coordinate system, Hermitian Cubic Shape function
How can I convert a local shape function given in terms of normalized/natural coordinates to global shape function in terms of x and y coordinate system in Hermitian cubic shape function.
Please explain using a 1 Dimensional problem.
What is a general expression for Hermitian cubic shape function.</p>


问:
Shooting Method with RK4 and Thermal Radiation
I am attempting to numerically solve the following problem. I decompose it into a system of two first order ODEs and then solve via the shooting method. I use the fourth order Runge-Kutta (RK4) method to solve each iteration of the shooting method.

$$k \frac{d^2T}{dx^2}=h(T-T_{\inf}) + \sigma(T^4-T^4_{surr})\\
T_{inf} = T_{surr} = 200 \, K\\
k = 1 \quad h = 0.05 \, \text{m}^{-2}\quad \sigma = 2.7*10^{-9} \, \text{K}^{-3}\,\text{m}^{-2}\quad T_L = 300 \, K \quad T_R = 400 \, K\\
L = 10 \, \text{m}$$

If I set $\sigma = 0$ then the solution converges for RK4 in just 3 iterations. I have checked it against the analytical solution for when $\sigma=0$.

However, if I keep the radiation term ($\sigma\neq0$) then RK4 quickly diverges during the first iteration of the shooting method. About half-way thru the solution domain it explodes.

Is there something special or unique about the following system that I am not aware? Some special case or sensitivity to initial conditions?

The system of first order ODEs that I am applying RK4 to is:
$$
\frac{dz}{dx} = w\\
\frac{dw}{dx} = h(z-T_{inf}) + \sigma (z^4-T^4_{inf})\\
\text{where} \quad z(0) = T_L\\
\text{and} \quad w(0) = 1 \; \text{(Guess 1)}
$$


问:
Is the given equation for Fractional Fourier Transform wrong?
In order to compute $$ \hat{f}_k = \sum_{m = 0}^{M - 1} e^{-2\pi i k m \theta} f_m, \ \ k = 0, ..., M - 1$$

for any $\theta$, my book states that this can be done using fractional Fourier transform as follows.

So, as you can see, they define the $c$-vector, then the $g$-vector, and then the final equation is a 2-times application of the FFT and 1-time application of the inverse FFT.

Is this equation correct? I have implemented it in code using some simple sequence $f_m$, and it does not seem to work.

The resulting vector $\hat{f}_k$ does indeed work ... but not for all $k$! It is wrong for the first few $k$ ($k = 0, 1, ...$) and only equals the correct values for the bigger $k$ ($k = M - 1, M - 2, ...$). So I am not sure what is going wrong here since the equation clearly is doing half the job but ...

I used R with its FFT-function.


问:
Fastest way to find roots of second order polynomial upto single decimal point?
What is fasted way to find roots of second order polynomial up to single decimal point using a program.</p>


问:
When does reduced integration lead to artificial zero energy modes in stiffness matrix?
This question relates to the topic of locking free finite element development.

In the case of application of reduced integration to global stiffness matrix for the Timoshenko beam element with quadratic shape function of field variables, I get an integrand up to 4th degree in case I calculate this integrand with 2 Gauss quadrature points should I expect additional zero eigenvalues to the resultant matrix?

To generalize, under what condition do I get additional zero eigenvalues when I integrate an $n$ degree polynomial function with reduced integration of Gauss quadratures?


问:
Unsteady diffusion equation with spatial finite elements and Forward Euler in time
I have solved the unsteady diffusion equation using piece-wise linear Finite elements(triangles) for spatial discretisation and Forward Euler for temporal discretisation. I have the following mesh convergence results. I'm not sure how to interpret them. The error was calculated using an analytical solution at time = 2 seconds. Time step, dt of 0.04 was used.

How do I verify the order of convergence, does the slope have to be 2 for second order accuracy in space?

Do we have anything like Courant number ($\frac{\Delta t}{\Delta h^2}\leq1/2$) to restrict the time step?


问:
exploding gradients in gradient descent procedure of multi-output ridge regression
Multi-output ridge regression:

$$W^{*}=\underset{W}{\arg \min } \frac{1}{\mathcal{N}}\|Y-WX\|_{F}^{2}+\lambda\|W\|_{F}^{2}$$

There are $Q$ outputs, $N$ samples, and $P$ covariates (features).

$\hat{Y}\in\mathbb R^{Q\times N}$ is the prediction matrix,

$$\hat{Y}=WX$$

where $W\in\mathbb R^{Q\times P}$ is the parameter matrix. These are the parameters we estimate in the regression.

$X\in\mathbb R^{P\times N}$ is the data matrix.

I have been trying to implement a multi-output ridge regression algorithm. I have derived the gradient descent update for the parameter $W$.The evaluation metric I am using is RRMSE(Relative root mean square error). I have observed that my RRMSE values are going very high for certain values of hyperparameters $\lambda$ (L2-regularization constant) and $\eta$ which is the learning rate in gradient descent.      My $\eta$ remains constant throughout the gradient descent procedure.

$$\operatorname{RRMSE}=\sqrt{\frac{\sum_{\left(\mathbf{x}_{i}, \mathbf{y}_{i}\right) \in D_{\text {test}}}\left(\hat{\mathbf{y}}_{i}-\mathbf{y}_{i}\right)^{2}}{\sum_{\left(\mathbf{x}_{i}, \mathbf{y}_{i}\right) \in D_{\text {test}}}\left(\hat{Y}-\mathbf{y}_{i}\right)^{2}}}$$

where $\left(\mathbf{x}_{i}, \mathbf{y}_{i}\right)$ is the $i$ th sample $\mathbf{x}_{i}$ with ground truth target $\mathbf{y}_{i}, \hat{\mathbf{y}}_{i}$ is the prediction of $\mathbf{y}_{i}$ and $\hat{Y}$ is the average of the targets over the training set $D_{\text {train}} .$ A lower RRMSE indicates better performance.

My attempts: I have read about the impact of un-normalized data on gradients but in my case the data is normalized.

I want some help with the two doubts I have right now.

$1.$What are the theoretically possible reasons for gradients to explode in gradient descent?

$2.$ Are there better ways to choose $\eta$ in order to avoid the exploding gradients problem.

Please provide your inputs on the above two doubts.


问:
Advantage of fractional Fourier transform over multiscale wavelet?
What could be the arguments of using fractional Fourier transform instead of multiscale wavelet for data analysis ?

Optimization of the good time-frequency domain parameter? good in the sens of best time-frequency domains that minimize spectral entropy of the data. But often, best entropy is achieved in a full 100% frequency domain.

Due to the physic world behavior, particle-path is locally described by fractional Fourier in N-slit problem context. Usefull for describing binary brownian data or differential binary brownian data ? ( sign(cumsum(randn(n,1))) or abs(diff(sign(cumsum(randn(n,1))))) )

Wavelet could be thinking as just a fast and efficient dyadic scheme of a particular sparse fractional transform ?

Continuous transformation from time to frequency: What would be the purpose of this transformation? Maybe to bring up patterns?

I really appreciate any strong argument about advantage of fractional Fourier Transform in comparison to wavelet.

(Our universe choosed fractional transform coding scheme for particles, there must be a reason of efficiency somewhere...)


问:
Time-dependent Schrodinger equation implementation in FEniCS
For our Bachelors thesis we're trying to solve the Schrodinger equation $i\partial_tu = -\nabla^2u+Vu$ in FEniCS. Given the domain $[-5, 5]^2$ with an initial value of $u_0(x, y)=e^{(-2(x^2+y^2))}$ and with $V(x,y)=\frac{x^2}{4}+9y^2$ at time $T=2$ the solution $|u|^2$ should look like the following image given by our supervisor.

This is however the result we get:

The equation is dicretized with Crank-Nicolson as follows:
(where $u=p+qi$)

For the real part:

$\frac{\Delta t}{2}(\nabla^2p^{n+1}-Vp^{n+1})-q^{n+1}=-q^n + \frac{\Delta t}{2}(Vp^n-\nabla^2p^n)$

for the imaginary part:

$\frac{\Delta t}{2}(\nabla^2q^{n+1}-Vq^{n+1})+p^{n+1}=-p^n + \frac{\Delta t}{2}(\nabla^2q^n-Vq^n)$

In the code below, $u[0]$ is $p$ and $u[1]$ is $q$

We're having trouble creating the correct residual when solving for $u$.

In FEniCS the equation is written in weak form and when multiplying with the test function we change sign of $dot(grad(), grad())$ to accommodate for multiplication with $i$. This is were we probably are at fault because we do not know how to construct the residual. When changing the sign of the gradient for the imaginary case the solution does not converge.

```<code>T = 2.0            # final time
num_steps = 200  # number of time steps
dt = T / num_steps # time step size

mesh = RectangleMesh(Point(-5, -5), Point(5, 5), 55, 55)
V = VectorFunctionSpace(mesh, 'Lagrange', 1, dim=2)
# Define initial value
u_0 = Expression(  ( 'exp(-2*(pow(x[0], 2)+pow(x[1], 2)))', '0'),
degree=1)
u_n = interpolate(u_0, V)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V,(0, 0), boundary)

u = Function(V)
v = TestFunction(V)
f = Expression(('(1/4)*pow(x[0], 2) + 9*pow(x[1], 2)', '0'), degree=2)

t=0

#Residual
ReVL = -u[1]*v[0]*dx + (dt/2)*(-dot(grad(u[0]), grad(v[0])) -
f[0]*u[0]*v[0])*dx
ReHL = -u_n[1]*v[0]*dx + (dt/2)*(f[0]*u_n[0]*v[0] + dot(grad(u_n[0]),
grad(v[0])))*dx
ImVL = u[0]*v[1]*dx + (dt/2)*(-dot(grad(u[1]), grad(v[1])) -
f[0]*u[1]*v[1])*dx
ImHL = u_n[0]*v[1]*dx + dt/2*(f[0]*u_n[1]*v[1] + dot(grad(u_n[1]),
grad(v[1])))*dx

FReal = ReVL - ReHL
FIm = ImVL - ImHL

Fny = FReal + FIm

J = derivative(Fny, u)
for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solve(Fny == 0, u, bc, J=J)
    # Update previous solution
    u_n.assign(u)
</code>```

We've also tried our luck doing as the example found in this link without success:
Gross-Pitaevskii in FEniCS

What are we doing wrong here?


问:
Comparison of convection time - theoretical value vs computed
This is a follow up to my previous post here,

I'm solving for convection in 1D

$$\frac{\partial C}{\partial t} = - v\frac{\partial C}{\partial x}$$

The discretization of the above equation is presented using a backward difference formula,
$$\frac{d C}{dt} = -v \frac{C_{i} - C_{i-1}}{\Delta x}$$

The boundary condition is imposed at the left end of the domain (i.e inlet)

$$\frac{d C}{dt}_{i=1} = -v \frac{C_{i} - C_{0}}{\Delta x}$$

and $C_0$ is specified
for observing the transient change in C over a time span
[0 tend].

Similar to the exponential kinetics

$C(t) = C_{\infty} (1-\exp(-\frac{t}{\tau}))$

presented for describing mass transport with diffusion mechanism, I'd to know if a similar expression can be derived for describing convection.


问:
Incorporating radiation boundary condition at the edge in finite difference
I am trying to solve the 2-d heat equation on a rectangle using finite difference method. I am confused as to how to incorporate non linear radiation boundary condition at the edge.

$-k\frac{\partial T}{\partial n}=\sigma\epsilon(T^4-T^4_{\infty})+h(T-T_{\infty})$

Is the following implementation reasonable?

$-k\frac{T^{n+1}_{i}-T^{n+1}_{i-1}}{\Delta x}=\sigma\epsilon((T^n_{i})^4-T^4_{\infty})+h(T^{n+1}_{i}-T_{\infty})$.

My logic is that taking sufficiently small time steps will give accurate results using this formulation.


问:
Why it doesn't matter to use velocity gradient or shear rate tensor to calculate wall shear stress?
The wall shear stress is defined based on Matyka et. al. eq. A.4:

$$\vec{\tau} = 2 \mu (\mathbf{S} \cdot \vec{n} - (\vec{n} \cdot \mathbf{S} \cdot \vec{n})\vec{n})$$

Where $\mathbf{S} = \frac{1}{2} (\nabla \mathbf{u} + (\nabla \mathbf{u})^{T})$ is shear rate tensor and $\mu$ is kinematic viscosity. There is one more popular definition of wall shear stress as:

$$\vec{\tau^{*}} = \mu \nabla \mathbf{u} \cdot \vec{n}$$

In the first look, these two definitions are really different, but I could simplify the Matyka's equation by remembering that wall shear stress is calculated at the wall and at the wall ($\partial \Omega$) we have no-slip boundary condition: $\mathbf{u} = 0$:

$$\vec{\tau} = \mu (\nabla \mathbf{u} \cdot \vec{n} + (\nabla \mathbf{u})^{T} \cdot \vec{n} - (\vec{n} \cdot (\nabla \mathbf{u} + (\nabla \mathbf{u})^{T}) \cdot \vec{n})\vec{n})$$

Now, we have:

$$(\nabla \mathbf{u})^{T} \cdot \vec{n} = \frac{\partial u_{j}}{\partial x_{i}} n_{j} \mathbf{e}_{i}$$

Here Einstein summation is used and $\mathbf{e}_{i}$ are unit vectors.

But:

$$\frac{\partial}{\partial x_{i}} (u_{j}n_{j}) = \frac{\partial u_{j}}{\partial x_{i}} n_{j} + \frac{\partial n_{j}}{\partial x_{i}} u_{j}$$

or:

$$(\nabla \mathbf{u})^{T} \cdot \vec{n} = \nabla (\mathbf{u} \cdot \vec{n}) - (\nabla \vec{n})^{T} \cdot \mathbf{u} = \nabla (\mathbf{u} \cdot \vec{n})$$

Because $\mathbf{u} = 0$ at the $\partial \Omega$ wall.

Also:

$$\vec{n} \cdot (\nabla \mathbf{u} + (\nabla \mathbf{u})^{T}) \cdot \vec{n} = \vec{n} \cdot \nabla \mathbf{u} \cdot \vec{n} + \vec{n} \cdot (\nabla \mathbf{u})^{T} \cdot \vec{n} = \vec{n} \cdot \nabla \mathbf{u} \cdot \vec{n} + \nabla (\mathbf{u} \cdot \vec{n}) \cdot \vec{n}$$

But:

$$\vec{n} \cdot \nabla \mathbf{u} \cdot \vec{n} = \frac{\partial u_{\alpha}}{\partial x_{\beta}} n_{\alpha} n_{\beta} = (\frac{\partial}{\partial x_{\beta}} (u_{\alpha} n_{\alpha}) - \frac{\partial n_{\alpha}}{\partial x_{\beta}} u_{\alpha})n_{\beta} = \vec{n} \cdot \nabla (\mathbf{u} \cdot \vec{n}) - (\vec{n} \cdot \nabla \vec{n}) \cdot \mathbf{u}$$

Again the last term is zero because of no-slip boundary condition so:

$$\vec{n} \cdot (\nabla \mathbf{u} + (\nabla \mathbf{u})^{T}) \cdot \vec{n} = 2 \nabla (\mathbf{u} \cdot \vec{n}) \cdot \vec{n}$$

So finally:

$$\vec{\tau} = \mu \nabla \mathbf{u} \cdot \vec{n} + \mu \nabla (\mathbf{u} \cdot \vec{n}) -2 \mu (\nabla (\mathbf{u} \cdot \vec{n})\cdot \vec{n})\vec{n}$$

So basically the difference of $\vec{\tau}$ and $\vec{\tau^{*}}$ is just:

$$\vec{\tau} - \vec{\tau^{*}} = \mu \nabla (\mathbf{u} \cdot \vec{n}) -2 \mu (\nabla (\mathbf{u} \cdot \vec{n})\cdot \vec{n})\vec{n}$$

I calculated the wall shear stress with both $\vec{\tau}$ and $\vec{\tau^{*}}$ and they are the same up to fourth digit, which means basically the term $\mu \nabla (\mathbf{u} \cdot \vec{n}) -2 \mu (\nabla (\mathbf{u} \cdot \vec{n})\cdot \vec{n})\vec{n}$ must be close to zero. Is there any reason that the gradient of normal component of velocity should be close to zero? I really appreciate some suggestions or insight here.


问:
Implementing Level Sets
I am trying to implement a basic level set program for fluid simulation in C++ and visualize it using SFML, so far I have the following simple program for propagating a circular curve:

https://gist.github.com/BhayanakMoth/da8250c849aeead098a61afbf655fb76

However, I am not sure if this the correct way.

Since I am working with a 128×128 grid, the resolution of the grid might not be enough and as a result tracking the curve at phi(x,y) = 0 is not possible in my current code. To deal with that I have tried defining the boundary like |phi(x,y)| < C, but I am not getting a continuous curve all the time unless the C is 0.6-0.9 which would probably be inaccurate. Is there another way to implement Level Set or how can I improve upon this?

Title: Implementing Level Sets


Title: Arnoldi Decomposition Algorithm

Body:
I try to get into GMRES via Arnoldi-Decomposition. For my understanding, I Implemented the Arnoldi-Decomposition in python.

```<code>import numpy as np
n = 5
A = np.random.rand(n,n)
v = np.random.rand(n)
M = np.zeros_like(A)
M[:,0] = v
for i in range(1,5):
   M[:,i] = np.matmul(A,M[:,i-1])

d = np.linalg.matrix_rank(M)

V = np.zeros((n,d))
H = np.zeros((d,d))

V[:,0] = v / np.linalg.norm(v)

for j in range(d-1):
    vj = V[:,j]

    Avj = np.matmul(A, vj)

    H[:j+1,j] =  np.matmul(Avj, V[:,:j+1])

    w = Avj - np.matmul(V[:,:j+1] ,H[:j+1,j])
    H[j+1,j] = np.linalg.norm(w)
    V[:,j+1]  = w/  H[j+1,j]
</code>```

I expected the result to be like $AV_k=H_kV_k$, i.e. to get an arnoldi-decomposition. But when I counter-checked the results, it seems to be off by one iteration:

```<code>>>np.matmul(A,V)
array([[ 1.74946273,  0.34244799, -0.03368428, -0.52017016, -0.36756412],
       [ 1.01841093,  0.09805327, -0.02055688, -0.82253407, -0.37112225],
       [ 0.81156945,  0.3374741 , -0.11957267,  0.17757239,  0.23021054],
       [ 0.99533205, -0.32793397, -0.12764172,  0.09559959,  0.14891387],
       [ 1.35827635, -0.02004282,  0.24724491, -0.07114606,  0.32279871]])
>>np.matmul(V,H)
array([[ 1.74946273,  0.34244799, -0.03368428, -0.52017016,  0.        ],
       [ 1.01841093,  0.09805327, -0.02055688, -0.82253407,  0.        ],
       [ 0.81156945,  0.3374741 , -0.11957267,  0.17757239,  0.        ],
       [ 0.99533205, -0.32793397, -0.12764172,  0.09559959,  0.        ],
       [ 1.35827635, -0.02004282,  0.24724491, -0.07114606,  0.        ]])
</code>```

So it seems to be that $AV_k = V_{k+1}H_{k+1}$. In my lecture notes, I found this relationship without a proof. So, can somebody show me a proof of it?


问:
Are linear, CTCS codes always stable?
I would like to solve some equations which basically look like this
$$\frac{\partial u}{\partial x}=F\left(v,\frac{\partial v}{\partial y},\frac{\partial^2 v}{\partial y^2}\right),$$
$$\frac{\partial v}{\partial x}=G\left(u,\frac{\partial u}{\partial y},\frac{\partial^2 u}{\partial y^2}\right).$$
My current plan is to solve this by replacing every derivative with a centred-difference equation and so the numerical scheme looks like this
$$u_{i+1,j}=u_{i,j} + \Delta x\,F\left(v_{i+1/2,j},v_{i+1/2,j+1},v_{i+1/2,j-1}\right)$$
$$v_{i+3/2,j}=v_{i+1/2,j} + \Delta x\,G\left(u_{i+1/2,j},u_{i+1/2,j+1},u_{i+1/2,j-1}\right),$$
where I define each $v_{i+1/2,j}$ to be located halfway between $u_{i,j}$ and $u_{i+1,j}$ in the $x$-direction. I initialisise $u$ and $v$ at $x=x_{min}$ and then iterate along in $x$. Is this scheme stable? Do you know if a scheme for solving such a problem is documented anywhere. I have coded it up and it seems to be stable for certain values of $\Delta x$ and $\Delta y$. However, I do not know if there is a fundamental problem with my scheme or if I have simply made a silly coding error somewhere.

Edit:

I have linked a pdf showing the expressions for F and G. They are highlighted in red text on page 3.

https://www.dropbox.com/s/43q5kkkqrrkbvrq/Normal_mode_code.pdf?dl=0

Title: Are linear, CTCS codes always stable?


问:
Residual of Poisson equation with periodic boundaries
I am trying to write a multigrid solver for Poisson's equation, $-\Delta u=f$, on the unit square, $\Omega=(0,1)^2$ with periodic boundaries. My primary source has been Multigrid by Trottenberg, Oosterlee and Schuller, 2001. Based on their discussion I'm using the standard 2nd order central difference for the derivatives:
$$-\Delta_h u(x_i,y_j) \approx \frac{1}{h^2} \left[ 4u(x_i,y_j) -u(x_{i+1},y_j)  - u(x_{i-1},y_j)-u(x_i,y_{j+1})-u(x_i,y_{j-1}) \right]$$
where
$$x_i=\frac{i}{h} \text{ } i =1,2,\dots ,n \text{ and } y_j=\frac{j}{h} \text{ } j=1,2,\dots , n$$
with the following boundary conditions:
$$ u(x_0, y_j) = u(x_n, y_j)$$
$$ u(x_1, y_j) = u(x_{n+1}, y_j) $$
$$ u(x_i, y_0) = u(x_i,y_n)$$
$$ u(x_i, y_1) = u(x_i,y_{n+1})$$

For existence, I'm enforcing the compatibility condition by altering the right hand side as follows:
$$f_h(x_i,y_j)=f(x_i,y_j)-h^2\sum_{i,j}f(x_i,y_j).$$
For uniqueness, I'm enforcing the global constraint:
$$\sum_{i,j}u_h(x_i,y_j)=0.$$

Now, here is where my question comes in:

Before writing the multigrid solver, I wanted to test my implementation of the smoother (in this case I'm using Gauss-Seidel, although eventually I plan to switch to Red-Black Gauss Seidel). This gives me the following iterative scheme:
$$u^{m+1}_{i,j}=\frac{h^2}{4}f_h(x_i,y_j)+\frac{1}{4}\left[ u^m_{i+1,j}+ u^{m+1}_{i-1,j} +u^m_{i,j+1}+u^{m+1}_{i,j-1} \right] \text{ }i,j=1,2,\dots n$$
where I have used the shorthand notation: $u_{i,j}=u(x_i,y_j)$. My first test case was:
$$u(x,y)=\cos(2\pi x)\cos(2 \pi y),$$
giving us the source term:
$$f(x,y)=8\pi^2\cos(2\pi x)\cos(2 \pi y).$$

The solution is periodic on $(0,1)^2$ and also satisfies the continuous compatibility condition:
$$\int_0^1\int_0^1 f(x,y) dx dy=0.$$
<s>However, it doesn't satisfy the discrete compatibility condition</s> It does saitsisfy the discrete compatibility condition (see edit):
$$ \sum_{i,j}f(x_i,y_j)=0,$$
and so I have imposed the condition above so that $\sum_{i,j}f_h(x_i,y_j)=0$.
What is strange to me is that when I solve the system $-\Delta_h u_{i,j}=f_h(x_i,y_j)$, <s> the residual does not converge to zero</s> (see edit). After noticing this, I removed the enforcement of the discrete compatibility condition and solved the system $-\Delta_h u_{i,j}=f(x_i,y_j)$. In this case the residual converges to zero. Since the correction to the RHS is $O(h^2)$ in both cases my numerical solution converges to the analytic solution as $n\to\infty$.

Since I plan to use a Multigrid method, I really want the residual to go to zero. For a practical problem I will only have data on the source term at discrete points, and thus, I will only know if the discrete compatibility condition is satisfied. Now for the question, based on my experiments it seems like my residual, $R_{i,j}$, should be:
$$R_{i,j}=f(x_i,y_j)+\Delta_h u_{i,j}$$
where $f(x_i,y_j)$ is the source term without the compatibility correction.

This confuses me when I look at the mathematics, though. I am solving the system $-\Delta_h u_{i,j}=f_h(x_i,y_j)$ and thus, to me it seems like I should be checking how well my iterative scheme satisfies this equation. In other words, it seems like the residual should be:
$$\bar{R}_{i,j}=f_h(x_i,y_j)+\Delta_h u_{i,j}.$$

Can anyone find the flaw in my reasoning?

EDIT:
Problem solved. I was evaluating the compatibility condition incorrectly. The correct compatibility condition is:
$$ \sum_{i=1}^n \sum_{j=1}^n f(x_i, y_j)=0 $$
NOT
$$ \sum_{i=0}^n \sum_{j=0}^n f(x_i, y_j)=0. $$
This counts the boundary terms twice.

The correct residual does converge to zero and is given by:
$$\bar{R}_{i,j}=f_h(x_i,y_j)+\Delta_h u_{i,j}.$$


问:
Numerical errors due to terms of the form $\frac{1}{r}$ (r goes to 0 at the boundary) while using finite difference method
I am trying to solve a system of differential equations using finite difference method.

There are few terms of the form $\frac{A(r)}{r}$, both $A(r)$ and r go to zero at the boundary. Analytically this term is well defined but numerically these terms are leading to error which very quickly break the simulation.

I am using leapfrog method if that matters. I read somewhere that leapfrog method has dissipation inbuilt so it should be stable, but it's clearly not the case.

Any advice on how to proceed?
Should I change my algorithm? Or use higher order schemes?

Thanks in advance


问:
Evaluating integral $F(r_2) = \frac{1}{r_2^n} \int_0^{r_2} r_1^n f(r_1)\ \mathrm{d}r_1$ without growing instability
I have the following expression to be numerically integrated in a vector-based library (e.g. numpy, MATLAB, etc),
$$
F(r_2) = \frac{1}{r_2^n} \int_0^{r_2} r_1^n f(r_1)\ \mathrm{d}r_1,
$$
where $n$ is an integer that can have value up to $40$ and $r_1$ and $r_2$ can have values of many orders of magnitude (e.g. $r_2 \in [10^{-10},10^4]$).

The problem is, the numerical integration method that I use introduces a slight numerical instability which makes the integration result $\lim_{r_2\rightarrow 0}\int_0^{r_2}r_1^n f(r_1)\mathrm{d}r_1$ not exactly $0$, but something around $10^{-13}$. So when the integration result is divided by $r_2^n\sim10^{-10n}$ where $n \geq 3$, the result of $F(r_2\rightarrow 0)$ becomes really big.

I am looking for a way to allow the multiplication of $r_1^n/r_2^n$ takes place before the integration, so the numerical instability does not grow, or anything that prevents the numerical instability from the integration from growing really big.

The question is: how can I integrate the expression above numerically using numpy, MATLAB, or any vector-based library without growing the numerical instability?

P.S: I avoid for-loops if possible because I program in Python and for-loops can be really slow.


问:
Decoupling of the fluid viscosity and incompressibility constraint in a partitioned scheme
I am new here and this is my very first question: I hope I respect all the criteria and rules.

I am just getting started with partitioned approaches for solving FSI problems. I am interested in a Chorin-Temam projection scheme for the fluid, with semi implicit coupling. This is the paper I am referring to: https://arxiv.org/abs/1711.10829

I am having some doubts regarding the following fact: if, for visualization purposes, I save the fluid velocity $\hat{\textbf{u}_f}$ that I compute in the explicit step, then what I am visualizing is not the physical velocity, correct? Because $\hat{\textbf{u}_f}$ does not satisfy the incompressibility constraint.

Sorry if this is an easy question, but I have just started to study partitioned approaches and would like to understand also the very basic aspects about them.


问:
Parallel plate capacitor as simple superposition of Coulomb-fields
I am trying to calculate and visualize the electric field inside and outside of a parallel plate capacitor by assuming a uniform distribution of point charges on each plate and by adding up the Coulomb-fields directly (I want this approach for educational purposes) with the following Python program.

Mathematically speaking, the electric field $\vec{E}(X)$ of a point charge is given by Coulombs law, i.e.

$$
\vec{E}(X) = \frac{1}{4\pi\epsilon_0} \frac{q}{r^3} \cdot \vec{r}
$$

where $X$ is the point where you measure the field strength and $\vec{r}$ is the vector pointing from the location $P$ of the point charge to $X$. $r$ denotes the absolute value of $\vec{r}$. $q$ denotes the charge of the point charge and $\epsilon_0$ is a constant.

Now I assume a uniform distribution of two parallel plates with distance $d$ between the plates. On one plate with positive charges on the other with negative ones.

The total electric field at point $X$, is then given by simple superposition of the Coulomb fields, i.e.

$$
E_{\mathrm{total}} = \vec{E}(X) = \frac{q}{4\pi\epsilon_0} \sum_i \frac{1}{r^3} \cdot \vec{r} \propto  q \sum_i \frac{1}{r^3} \cdot \vec{r}
$$

Since my program should only demonstrate the principle I set $\frac{1}{4\pi\epsilon_0}$ to $1$ and choose $q=1$ for the positive plate and $q=-1$ for the negative one.

To translate the notation to my program I denote the position of a point charge by the vector $\vec{a}$ with coordinates $a_0,a_1,a_2$, and introduce a cartesian coordinate system with coordinates $x$ (to the right), $y$ (to the top) and $z$ (into the plane of view).

Then I can write the field of one charge in componentes as:

$$
\begin{align}
E_x &= q\cdot(x-a_0)/r^3 \\
E_y &= q \cdot (y-a_1)/r^3
\end{align}
$$

where

$$
r^3 = ((x-a_0)^2 + (x-a_1)^2 + (z-a_2))^{1.5}
$$

I look only at the plane $z = 0$ and since I draw in 2D (projection to the plane $z = 0$) I omit $E_z$ (but I am not sure, if this introduces a mistake).

Then I use a loop to equally distribute the charges on my planes and calculate the total field as described above.

# Imports
import matplotlib
from numpy import *
from pylab import *
import matplotlib.colors as colors

def Cap(N = 10,plus=1,minus=1,app=''):
    # Aspect ratio
    aspekt='auto'

    ## Field between the plates
    xMin,xMax = -0.15,0.15
    yMin,yMax = -10,10
    zMin,zMax = -10,10

    ## Field outside
    # xMin,xMax = -3,3
    # yMin,yMax = -3,3

    ## Plate distance

    d = 0.1

    # Case of plates very far apart
    # d = 50
    # xMin,xMax=-d*1.25,d*1.25
    # yMin,yMax=-d*1.25,d*1.25
    # aspekt = 'equal'

    ## Grids
    xG,yG = meshgrid(linspace(xMin,xMax,18),linspace(yMin,yMax,18))
    xGM,yGM = meshgrid(linspace(xMin,xMax,500),linspace(yMin,yMax,500))

    # Class that represents a point charge
    class charge:
        def __init__(self, q, pos):
            self.q=q
            self.pos=pos

    # Coulomb Law (I only look a the plane z = 0)
    def E_point_charge(q, a, x, y):
        z = 0
        r3 = ((x-a[0])**2+(y-a[1])**2 + (z-a[2])**2)**(1.5)
        return q*(x-a[0])/r3,q*(y-a[1])/r3

    # Total electric field
    def E_total(x, y, charges):
        Ex, Ey=0, 0
        for C in charges:
            E=E_point_charge(C.q, C.pos, x, y)
            Ex=Ex+E[0]
            Ey=Ey+E[1]
            #print(Ex)
        return [Ex, Ey]

    ## Initialize

    close('all')
    fig = figure()

    chargesPlus=[]
    chargesMinus=[]

    ## Generate charges (1/N**2 to have for every N the same total charge on each plate)
    for i in arange(yMin*0.9,yMax*0.9,0.9*(yMax-yMin)/N):
        for j in arange(zMin*0.9,zMax*0.9,0.9*(zMax-zMin)/N):
            if plus != 0:
                chargesPlus.append(charge(1/N**2,[d/2,i,j]))
                print("Charge + at %s,%s" %(i,j))
            if minus != 0:
                chargesMinus.append(charge(-1/N**2,[-d/2,i,j]))
                print("Charge - at %s,%s" %(i,j))
    charges = chargesPlus + chargesMinus

    #
    E_totalX,E_totalY = E_total(xG,yG,charges)
    EE = sqrt(E_totalX**2 + E_totalY**2)
    E_DirX,E_DirY= E_totalX/EE,E_totalY/EE

    #
    E_totalXM,E_totalYM = E_total(xGM,yGM,charges)
    EEM = sqrt(E_totalXM**2 + E_totalYM**2)

    # Plots
    ax = subplot(111)
    fig.add_subplot(ax)

    colormap='Spectral_r'

    # Colormap
    I = ax.imshow(EEM,extent=[np.min(xGM),np.max(xGM),np.min(yGM),np.max(yGM)],cmap=colormap, norm=colors.LogNorm(vmin=EEM.min(),vmax=EEM.max()),alpha=0.5,origin='lower',interpolation='bilinear')

    # Vector field
    Q = ax.quiver(xG,yG,E_DirX,E_DirY,color='gray',alpha=0.7,width=0.002,scale_units='inches',scale=4)

    # Colorbar
    fig.colorbar(I,extend='max',orientation='horizontal')

    # Plot point charges
    for C in charges:
        if C.q>0:
            plot(C.pos[0], C.pos[1], 'bo', ms=10*sqrt(C.q))
        if C.q<0:
            plot(C.pos[0], C.pos[1], 'ro', ms=10*sqrt(-C.q))

    # Axis sessings
    xlabel('$x$')
    ylabel('$y$')
    ax.set_aspect(aspekt)

    ax.set_xlim(xMin, xMax)
    ax.set_ylim(yMin, yMax)

    savefig("Capacitor_N_%05d_%app.pdf" %(N,app),dpi=fig.dpi,papersize='a4')
    show()

Cap(200)
```

<p>However the computation time is very long, so is there better approach to do this under the constraint to use the direct summation of Coulomb fields (because of educational purposes) and to use Python as a tool.</p>

<p>The most time consuming step seems to be the calculation of the total field for the colormap plot.</p>

<p><strong>Edit</strong></p>

<p>I tried coarser grid for the colormap and then tried to use the <code>timeit</code> module to find the most time consuming part, but this doesn't seem to report reliable times. Maybe because of asynchronous processes. To show the problem I did a screencast for <span class="math-container">$N = 70$</span> here: https://www.pastefile.com/mBh38p</p>

<p>As you can see the plotting itself seems to be the problem. (I used a old i7-3520M CPU in my screencast)</p>

<p>The modified code looks like:</p>

```#!/usr/bin/env python3

# Imports
import matplotlib
from numpy import *
from pylab import *
import matplotlib.colors as colors

from timeit import default_timer as timer
from datetime import timedelta

def Cap(N = 10,plus=1,minus=1,app=''):
    print('### N = %s ###' %N)

    # Aspect ratio
    aspekt='auto'

    ## Field between the plates
    xMin,xMax = -0.15,0.15
    yMin,yMax = -10,10
    zMin,zMax = -10,10

    ## Field outside
    # xMin,xMax = -3,3
    # yMin,yMax = -3,3

    ## Plate distance

    d = 0.1

    # Case of plates very far apart
    # d = 50
    # xMin,xMax=-d*1.25,d*1.25
    # yMin,yMax=-d*1.25,d*1.25
    # aspekt = 'equal'

    ## Grids
    xG,yG = meshgrid(linspace(xMin,xMax,18),linspace(yMin,yMax,18))
    xGM,yGM = meshgrid(linspace(xMin,xMax,40),linspace(yMin,yMax,40))

    # Class that represents a point charge
    class charge:
        def __init__(self, q, pos):
            self.q=q
            self.pos=pos

    # Coulomb Law (I only look a the plane z = 0)
    def E_point_charge(q, a, x, y):
        z = 0
        r3 = ((x-a[0])**2+(y-a[1])**2 + (z-a[2])**2)**(1.5)
        return q*(x-a[0])/r3,q*(y-a[1])/r3

    # Total electric field
    def E_total(x, y, charges):
        Ex, Ey=0, 0
        for C in charges:
            E=E_point_charge(C.q, C.pos, x, y)
            Ex=Ex+E[0]
            Ey=Ey+E[1]
            #print(Ex)
        return [Ex, Ey]

    ## Initialize

    close('all')
    fig = figure()

    charges = []

    ## Generate charges (1/N**2 to have for every N the same total charge on each plate)
    print("Generate charges")
    start = timer()
    n = 0
    for i in arange(yMin*0.9,yMax*0.9,0.9*(yMax-yMin)/N):
        for j in arange(zMin*0.9,zMax*0.9,0.9*(zMax-zMin)/N):
            if plus != 0:
                charges.append(charge(1/N**2,[d/2,i,j]))
                #print("Charge + at %s,%s" %(i,j))
                #print(n)
                n += 1
            if minus != 0:
                charges.append(charge(-1/N**2,[-d/2,i,j]))
                #print("Charge - at %s,%s" %(i,j))
                #print(n)
                n +=1
    end=timer()
    print(timedelta(seconds=end-start))

    print("Calculating E_total for Quiver Plot")

    start=timer()
    E_totalX,E_totalY = E_total(xG,yG,charges)
    EE = sqrt(E_totalX**2 + E_totalY**2)
    E_DirX,E_DirY= E_totalX/EE,E_totalY/EE
    end=timer()
    print(timedelta(seconds=end-start))

    print("Calculating E_total for colormap")
    start=timer()
    E_totalXM,E_totalYM = E_total(xGM,yGM,charges)
    EEM = sqrt(E_totalXM**2 + E_totalYM**2)
    end=timer()
    print(timedelta(seconds=end-start))

    # Plots
    ax = subplot(111)
    fig.add_subplot(ax)

    #colormap='Spectral_r'
    colormap='coolwarm'

    print("Making colormap")
    start=timer()
    # Colormap
    I = ax.imshow(EEM,extent=[np.min(xGM),np.max(xGM),np.min(yGM),np.max(yGM)],cmap=colormap, norm=colors.LogNorm(vmin=EEM.min(),vmax=EEM.max()),alpha=0.5,origin='lower',interpolation='bilinear')

    # Colorbar
    fig.colorbar(I,extend='max',orientation='horizontal')
    end=timer()
    print(timedelta(seconds=end-start))

    print("Making quiver plot")
    start=timer()
    # # Vector field
    #Q = ax.quiver(xG,yG,E_DirX,E_DirY,color='gray',alpha=0.7,width=0.002,scale_units='inches',scale=4)
    Q = ax.quiver(xG,yG,E_DirX,E_DirY)
    end=timer()
    print(timedelta(seconds=end-start))

    ax.set_title('N = %s' %N)

    # Plot point charges
    for C in charges:
        if C.q>0:
            plot(C.pos[0], C.pos[1], 'bo', ms=10*sqrt(C.q))
        if C.q<0:
            plot(C.pos[0], C.pos[1], 'ro', ms=10*sqrt(-C.q))

    # Axis sessings
    xlabel('$x$')
    ylabel('$y$')
    ax.set_aspect(aspekt)

    ax.set_xlim(xMin, xMax)
    ax.set_ylim(yMin, yMax)

    print("Savefig")
    savefig("Capacitor_N_%05d_%app.pdf" %(N,app),dpi=fig.dpi,papersize='a4')

    print("Rendering")
    show()

Cap(70)

# N = arange(30,300,30)

# for i in N:
#     Cap(i)
```


问:
Optimization of centers and radii of circles under the non-collision constraint
I want to optimize a function w.r.t. $n$ circles parametrized by centers and radii:

$$\min_{C, R} f(C,R)$$

where $C\in\mathbb R^{n \times 2}$ and $R\in\mathbb R^n$. For example,

```<code>C[1,1]: x-coordinate of the first circle
C[1,2]: y-coordinate of the first circle
R[1]  : radius of the first circle
</code>```

While optimizing, I want to impose a constraint that those circles do not collide with each other. Anyone to help with this?


问:
Alternatives to breadth-first-search in 3D grid cluster detection?
I've got a question about a good way to find the quickest algorithm for my problem:

problem:
I've got a 3D cubical grid containing voxels that are either 1 or 0. It is stored as a flattened array. If multiple 1-voxels are neighboring each other, they form a cluster together.
I want to detect clusters smaller than a given size (n).

current solution (slow):

For each 1-voxel, I'm doing a breadth first search to find how many neighbors it has, and how many neighbors the neighbors have, etc, until reaching a certain amount of neighbors (n).

this is horribly slow when setting n to my preferred amount n>50, because the total cube edge size is 64

Any different ways to do this that are not breadth first search? Everything is open: could even be something from computer vision like 3d pooling or whatever, but I'm a little stuck. Or am I better off trying to optimize my current approach? If interested, I posted my c# code in code review

Thanks a bunch in advance!


问:
Solving a nonlinear problem with a very small components with finite element method
In solving nonlinear hyperelastic solid mechanics problems, to converge to the correct solution we need to do step-by-step loading which makes the deformation at each step very small (for my particular problem $\mathcal{O}(e^{-11})$). Moreover, due to the material properties, the nonlinear weak form and consequently the linear system that should be solved at each step of the nonlinear solver have very small components. In my experience, this makes the iterative solvers very inefficient especially considering that many components are close to machine epsilon precision and the relative norm is too small and is not a good measure for the convergence. I would appreciate if you can help me to know

1) is what I am observing true or I am making a mistake somewhere?

2) what is the smart way to deal with such matrices that the weak form and consequently the resulted linear system have very small components?

What I am doing now to solve this nonlinear system is to reduce the absolute convergence norm criteria to a very small number $\mathcal{O}(e^{-13})$ when the absolute norm for 0 iteration is $\mathcal{O}(e^{-9})$. And use a relatively small relaxation value.

Thank you, everybody, for your help.


问:
How to achieve (approx) unit scaling of a non-linear diffusion (heat) equation with a wildly varying diffusion coefficient?
I have numerical issues with a poorly scaled one-dimensional non-linear diffusion equation in physical co-ordinates

$$ \frac{\partial{u}}{\partial{t}}(x,t) = \frac{\partial}{\partial{x}}\left(D(u)  \frac{\partial{u}(x,t)}{\partial{x}}\right)$$

$$x \in (0,L),\  t \in (0,t_\text{f}) $$

with initial conditions
$$u(x,0) = u_0(x)$$ and Double-Neumann Boundary Conditions
$$D(u)\frac{\partial{u}}{\partial{x}}_{\Large\mid \normalsize{x=0}} = f(t); \quad D(u)\frac{\partial{u}}{\partial{x}}_{\Large\mid \normalsize{x=L}} = 0$$.

In numerically solving this PDE, I end up with convergence errors due to poor scaling of the problem with $u \approx \mathcal{O}(10^3)$, $x \approx \mathcal{O}(10^-6)$, and the most problematic of them all $D(u) \approx \mathcal{O}(10^{-12})\text{ to }\mathcal{O}(10^{-16}) $.

By applying standard scaling techniques $\bar{t} = \frac{t}{t_\text{c}}$, $\bar{x} = \frac{x}{x_\text{c}}$, $\bar{u} = \frac{u}{u_\text{c}}$, $\bar{D} = \frac{D}{D_\text{c}}$, and setting $x_\text{c} = L$, $D_\text{c} = D_{\max}$, $t_\text{c} = \frac{L^2}{D_{\max}}$ and $u_\text{c} = {u_{\max}}$, I was able to arrive at a better-scaled non-dimensionalised version of the PDE as

$$ \frac{\partial{\bar{u}}}{\partial{\bar{t}}}(\bar{x},\bar{t}) = \frac{\partial}{\partial{\bar{x}}}\left(\bar{D}(\bar{u})  \frac{\partial{\bar{u}}(\bar{x},\bar{t})}{\partial{\bar{x}}}\right)$$

in the unit interval $\bar{x} \in (0,1)$ with $\bar{u} \in (0,1)$. However, due to the 4 orders of magnitude variation in $\bar{D}$ as well as the poor resolution of its values, it doesn't look like the linear scaling $\bar{D} = \frac{D}{D_\text{max}}$ is the correct one to use.

Can someone please help to scale this PDE for better numerical conditioning?


问:
how to construct chebyshev differentiation matrix of even/odd polynomial
I know there is efficient way to construct chebyshev differentiation matrix $D_N$. The idea is to fit the function $f(x)$ with a polynomial of $p_N(x)$ of order $N$, and use the derivative $p_N'(x)$ to approxiate $f'(x)$ at chebyshev points $x_0,..,x_N$.

I wonder if we know f(x) is an even function on $[-1,1]$, and we want to approximate $f(x)$ with a polynomial of only even terms $p_{2N}(x) = a_0 + a_2 x^2 + a_4 x^4 + ... a_{2N}x^{2N}$, is there an efficient way to construct the differentiation matrix $D_{2N}$?


问:
Plot FDE Amplification Factor in Complex Plane
I am reviewing stability of finite difference methods. I am using the von Neumann method to study a linear PDE and an explicit finite difference scheme. The goal is to develop an expression for an amplification factor $G$.

I've read that if $G$ is a complex function it can be plotted in the complex plane to more easily determine where the solution of the finite difference equation will be stable. For example, consider if $G$ was

$$
G = (1-2d) + 2d\cos\theta
$$

where $d$ is determined from my grid and time step sizes.

What is the template or formulas used to re-write this in terms of a circle (or some other shape, whatever it is) in the complex plane? Once I know that, I can plot it and hopefully tell where the stable region is.

My reference indicates this example of $G$ is an oscillation centered at $(1-2d,0i)$ with an amplitude of $2d$. It isn't clear to me how that was determined or how you could use that to plot it.


问:
Stress during unloading in FEM simulation of an elasto-plastic truss element
I am simulating the behaviour of a single bar composed of two node points by fixing one point and applying a force on the oder node with $F_x = 1$ and fixing $u_y=0$. I am assuming the material to behave elasto plastic with isotropic hardening.

My loading simulation for $F_x = 0$ to $F_x = 1$ with increment $1/10$ seems to be okay, my unloading process with same increment size is somehow weird. Maybe I am wrong but the stress should be $0$ for $F_x = 0$ again, right?.


问:
Unconditionally stable numerical method for 1st order non-linear coupled ODEs?
I am attempting to numerically solve the following system of ODEs:

$$\begin{gather}\frac{dT_1}{dt} = f_1(T_1,T_2), \quad T_1(t=0)=T_{1,0} \\[3pt]
\frac{dT_2}{dt} = f_2(T_1,T_2), \quad T_2(t=0)=T_{2,0}\end{gather}$$

The functions $f_1$ and $f_2$ have $T$ and $T^4$ terms and therefore are non-linear. I tried the forward Euler method but this becomes unstable even for very small $\Delta t$.

Is there any unconditionally stable method of numerical integration to solve this system?

OwnerDisplayName: cryonole


问:
Convergent Finite Difference Scheme for Parabolic Equation
Consider the PDE $$u_t = b_{11}u_{xx} + 2b_{12}u_{xy} + b_{22}u_{yy},$$
where $b_{11}, b_{22} > 0$, and $b_{12}^2 < b_1b_2$.
In Strikwerda's book, the ADI scehme \begin{align*}
\left( 1 - \frac{k}{2}  b_{11} \delta_x^2\right)v^{n+1/2} &= \left( 1 + \frac{k}{2}  b_{22} \delta_y^2\right)v^{n} + kb_{12} \delta_{0x} \delta_{0y} v^n \\
\left( 1 - \frac{k}{2}  b_{22} \delta_x^2\right)v^{n+1} &= \left( 1 + \frac{k}{2}  b_{11} \delta_x^2\right)v^{n+1/2} + kb_{12} \delta_{0x} \delta_{0y} v^{n+1/2} \\
\end{align*}
is stated to be convergent. I was able to show that it was consistent, but I am struggling with stability. Making the substitution $v^{n} = g^n e^{il \theta}e^{i m \phi}$, and $v^{n+1/2} = \tilde{g}g^n e^{il \theta}e^{i m \phi}$,
\begin{align*}
\left(1 + 2\mu b_{11} \sin^2{ (\theta/2
)} \right) \tilde{g} &= \left(1 - \mu( 2b_{22} \sin^2{ (\theta/2
)} + b_{12}\sin{\theta}\sin{\phi}) \right) \\
\left(1 + 2\mu b_{22} \sin^2{ (\theta/2
)} \right) g &= \left(1 - \mu( 2b_{11} \sin^2{ (\phi/2
)} + b_{12}\sin{\theta}\sin{\phi}) \right)\tilde{g}, \\
\end{align*}
where $\mu = \frac{k}{h^2}$. This implies that $$g = \frac{\left(1 - \mu( 2b_{11} \sin^2{ (\theta/2
)} + b_{12}\sin{\theta}\sin{\phi}) \right)\left(1 - \mu( 2b_{22} \sin^2{ (\phi/2
)} + b_{12}\sin{\theta}\sin{\phi}) \right)}{\left(1 + 2\mu b_{11} \sin^2{ (\theta/2
)} \right)\left(1 + 2\mu b_{22} \sin^2{ (\phi/2
)} \right)}$$
Ultimately, I want to find some restriction on $\mu$ such that $\left| g \right| \leq 1 $. Since most ADI schemes are unconditionally stable, I expect that this fraction is automatically bounded above by $1$, but I'm not sure how to show that. Any help would be much appreciated.


问:
Linearize non-linear PDE with BCs to hyperbolic problem: How does linearization affect BCs?
I am working with the Shallow Water equations that is a system of non-linear PDEs that simulate water waves propagation on some domain, in my cases the $x$-axis. I have here a GIF showing the results (is it possible to insert the GIF into the post? Because I have not figured out how). I will try to keep it short and what I am about to do is linearize the system below and then solve the linearized hyperbolic system analytically (my problem is how I should deal with the boundary conditions for the characteristic variables I obtain). So we have

$$
\begin{cases}
h_{t} + (hv)_{x} = 0, \\[2mm]
(hv)_{t} + (hv^{2} + \frac{1}{2}gh^{2})_{x} = 0, \\[2mm]
\text{on} \quad (x,t)\in [0,L]\times[0,\infty)
\end{cases}
\quad \iff
\begin{pmatrix}h \\ hv\end{pmatrix}_{t} + \begin{pmatrix}hv \\ hv^{2} + \frac{1}{2}gh^{2}\end{pmatrix}_{x} = u_{t} + f_{x}(u) = 0.
$$

The initial and boundary conditions are given by

$$
\begin{cases}
\text{Initial conditions:} \quad h(x,0) = H + \varepsilon e^{-(x-L/2)^{2}/w^{2}} \quad \text{and} \quad v(x,0) = 0, \\[2mm]
\text{Boundary conditions:} \quad v(0) = 0 \quad \text{and} \quad v(L) = 0.
\end{cases}
$$

Now if a linearize this system around some constant state, say $u_{0} = (h_{0},h_{0}v_{0})$, then I get a hyperbolic system. If we denote

$$
q = u - u_{0}
$$

the perturbation from this constant state, then we may linearize the system by the following approach:

$$
f(u) = \begin{pmatrix}hv \\ hv^{2} + \frac{1}{2}gh^{2}\end{pmatrix} = \left.\begin{cases} u_{1} &= h \\ u_{2} &= hv  \end{cases}\right\} = \begin{pmatrix}u_{2} \\ \frac{u_{2}^{2}}{u_{1}} + \frac{1}{2}gu_{1}^{2}\end{pmatrix}
\quad \implies
f'(u) = \begin{pmatrix}0 & 1\\ - \frac{u_{2}^{2}}{u_{1}^{2}} + gu_{1} & 2\frac{u_{2}}{u_{1}}\end{pmatrix}.
$$

If we evaluate $f'(u)$ in $u_{0}$ we obtain the linearized system matrix $A$,

$$
A = f'(u)\vert_{u=u_{0}} = \left.\begin{pmatrix}0 & 1\\ - \frac{u_{2}^{2}}{u_{1}^{2}} + gu_{1} & 2\frac{u_{2}}{u_{1}}\end{pmatrix}\right\vert_{u=u_{0}} = \begin{pmatrix}0 & 1\\ - v_{0}^{2} + gh_{0} & 2v_{0}\end{pmatrix}.
$$

Furthermore if we choose $v_{0}=0$ we obtain the linearized system

$$
q_{t} + \begin{pmatrix}0 & 1\\ gh_{0} & 0\end{pmatrix}q_{x} = 0
$$

which can be seen to be hyperbolic. Now to my question: It is possible to transform this system and decouple it into the so called characteristic equations (these are two partial differential equations modeling advection). But what will the boundary conditions be for each of these equations? I have read in Leveque $\S$ 7.3 that these equations each have one boundary condition, but what are the domains of definition for each characteristic variable?Is it not still that $x\in[0,L]$? For me it seems that one bounary condition would not require each equation to stay in the interior.

Let me show in mathematical terms what I mean: After linearization we may transform the system by diagnoalizing $A$. We obtain

$$
A = \begin{pmatrix}\xi_{1} & \xi_{2}\end{pmatrix}\begin{pmatrix}\lambda_{1} & 0 \\ 0 & \lambda_{2}\end{pmatrix}\begin{pmatrix}\xi_{1} & \xi_{2}\end{pmatrix}^{-1} \quad \text{where} \quad \lambda_{1} = -\sqrt{gh_{0}},\; \xi_{1} = \begin{pmatrix}-1 \\ \sqrt{gh_{0}}\end{pmatrix},\; \lambda_{2} = \sqrt{gh_{0}},\; \xi_{2} = \begin{pmatrix}1 \\ \sqrt{gh_{0}}\end{pmatrix}.
$$

Now if we denote $R = \begin{pmatrix}\xi_{1} & \xi_{2}\end{pmatrix}$ and set $z = R^{-1}q$ we obtain the decoupled system

$$
z_{t} + \begin{pmatrix}-\sqrt{gh_{0}} & 0 \\ 0 & \sqrt{gh_{0}}\end{pmatrix}z_{x} = z_{t} + \Lambda z_{x} = 0.
$$

The initial conditions for this transformed system become

$$
\begin{cases}
\text{Initial conditions:} \quad z(x,0) = R^{-1}u(x,0) = R^{-1}\begin{pmatrix}h(x,0)\\ h(x,0)v(x,0)\end{pmatrix} = \frac{1}{2}\begin{pmatrix}-1\\ 1\end{pmatrix}h(x,0)
\end{cases}.
$$

But what happens to the boundary conditions? Are they given by

$$
\begin{cases}
\text{Boundary condition } x=0:\quad z(0,t) = R^{-1}u(0,t) = R^{-1}\begin{pmatrix}h(0,t)\\ h(0,t)v(0,t)\end{pmatrix} = \frac{1}{2}\begin{pmatrix}-1\\ 1\end{pmatrix}h(0,t), \\[2mm]
\text{Boundary condition } x=L:\quad z(L,t) = R^{-1}u(L,t) = R^{-1}\begin{pmatrix}h(L,t)\\ h(L,t)v(L,t)\end{pmatrix} = \frac{1}{2}\begin{pmatrix}-1\\ 1\end{pmatrix}h(L,t)\;?
\end{cases}
$$

Because there is no boundary value prescribed for $h(0,t)$ or $h(L,t)$. So how do I go about solving this analytically from here, how do I treat reflecting boundary conditions like this on a finite domain for a system of hyberbolic equations?

\Best regrads

EDIT
I show here the result I obtained for the linearized system:

As you can see the solution just flows through the boundaries. How do I prescribe values at the boundaries here exactly so that the waves reflect? I have found that for reflecting boundary conditions there are two relations (I do not know if this is mentioned in Leveque), namely

$$
\begin{align}
z_{2}(0,t) = -z_{1}(0,t) \quad \text{and} \quad z_{1}(L,t) = -z_{2}(L,t),
\end{align}
$$

but I do not know how these two relations impose reflectiveboundary conditions or how to do it.


问:
Automatically generate constraints for trajectory optimization
This is a follow up to my previous post here

I'm interested in performing trajectory optimization from the problem mentioned in abov link.

I want to supply the following as dynamical constraints to  MATLAB's fmincon.

$$ \frac{d \phi}{dt}  = -M^TDM\phi   + W \hspace{1cm}(1)$$

$\phi$ is a vector with variables [ A B C D E F].

I've been referring to this article that illustrates how trapezoidal collocation is done

For example, if $\phi$ is a scalar

$$\dot\phi = f \hspace{1cm} (2)$$

integrating the above from time $t_k$ to $t_{k+1}$ and approximating using trapezoidal quadrature gives,

$$ \phi_{k+1} = \phi_k + \frac{1}{2}h_k(f_{k+1} + f_k) \hspace{1cm} (3)$$

This will translate into a set of non-linear equality constraints (ref) for k = 0 to N-1.

I would like to know how to generate these constraints automatically. Is it possible to generate these constraints automatically and pass it as a handle to fmincon. Because it is possible to manually write for a single variable but for multivariable it gets complex.

I also looked at a repository available here that can be used in MATLAB. However, I am not sure if this can be used to generate constraints.

Also, I would like to understand how equation (3) should actually be written as an equality constraint.
Should one solve the ode system in (1) and use the values of $\phi$ at the collocation points to evaluate the RHS of equation (3)? What happens to RHS?

Any suggestions/explanations for the above questions will be really useful.

EDIT:
I tried following the example code given here to understand how ceq (equality constraint) has to be implemented

For example, in the code provided in the above link

```<code>function [ c, ceq ] = double_integrator_constraints( x )
</code>```

is given.

Could someone explain how the input argument, x, is computed?
Also, I don't understand how ydesiredend has to be mentioned in equality constraints.

In my problem, only the initial conditions are defined for the dynamical constraints and the terminal conditions are not available. In that case , I am not sure how the equality constraints have to be defined.

EDIT2:
Update: I could solve this problem using fminunc and lsqnonlin.

```<code>Dhat0 = %input vector
% fun   = @objfun;
% [Dhat,fval] = fminunc(fun, Dhat0)

%% lsqnonlin
Dhat = lsqnonlin(@(Dhat) objfun(Dhat),Dhat0)

function f = objfun(Dhat)

%% Integrator settings
tspan = %tspan
options = odeset('abstol', 1e-10, 'reltol', 1e-9);

%% generate exact solution
    phi0 = % initial condition vector
    [t, phi]  = ode15s(@(t,phi) exact(t,phi), tspan , phi0 ,options);

%% generate approximate solution

    [t, phi_tilde]  = ode15s(@(t,phi_tilde) approx(t,phi_tilde, Dhat), tspan , phi0 ,options);

%% objective function for fminunc
    % diff = (phi - phi_tilde).*(phi - phi_tilde);
    % f = sum(diff, 'all')

%% objective function for lsqnonlin
    f  = phi - phi_tilde
end
</code>```

I am still interested in understanding how the constrained optimization problem can be solved. I'm trying to set up the problem this way.

I would like to understand how noncol should be called from fmincon.

```<code>nonlcon = @defects;
Dhat= fmincon(@objfun,Dhat0,A,b,Aeq,beq,lb,ub,nonlcon)
</code>```

For my system,
A,b,Aeq,beq,lb,ub = []

But, I am not sure from where to pass these arguments for  defects(dt,x,f).

```<code>function [c ceq] = defects(dt,x,f)
    % ref: https://github.com/MatthewPeterKelly/OptimTraj
    % This function computes the defects that are used to enforce the
    % continuous dynamics of the system along the trajectory.
    %
    % INPUTS:
    %   dt = time step (scalar)
    %   x = [nState, nTime] = state at each grid-point along the trajectory
    %   f = [nState, nTime] = dynamics of the state along the trajectory
    %
    % OUTPUTS:
    %   defects = [nState, nTime-1] = error in dynamics along the trajectory
    %   defectsGrad = [nState, nTime-1, nDecVars] = gradient of defects

nTime = size(x,2);

idxLow = 1:(nTime-1);
idxUpp = 2:nTime;

xLow = x(:,idxLow);
xUpp = x(:,idxUpp);

fLow = f(:,idxLow);
fUpp = f(:,idxUpp);

% This is the key line:  (Trapazoid Rule)
defects = xUpp-xLow - 0.5*dt*(fLow+fUpp);
ceq = reshape(defects,numel(defects),1);
c = [];
end
end
</code>```

Any suggestions will be really helpful


问:
Compute mass matrix in vibrations problem by using finite element method
I have to compute the mass matrix of a Hexahedral mesh.

There are 3 methods to compute mass matrix. I'm interested in one method which consists of dividing the mass of the element by the number of nodes and attributes this mass at each node for each degree of freedoms.

My problems is that: I don't know how to compute the volume of the element and then multiply it by the density to calculate mass.

In 2D, I have the solution, but not in 3D.

I'm open to any help on this subject due to the fact that I don't know any article on this subject.


问:
Lexicographically order matrix into a vector
I am trying to implement the algorithm contained in this article here. It is about solving a 2 and 2.5D Fredholm integral, focused on bidimensional NMR experiments. I've made significant progress, despite this being my very first time doing any computational science. The output I get is clearly wrong, but seems shifted from what would be ideal. I am trying to stamp out problematic areas, and am wondering if I interpreted this passage wrong:

In principle, the 2-D problem in (2) can be reduced to a 1-D problem $m_r = K_0f_r + \epsilon_r \quad r=1,...,R$. Here, the vectors $m_r$, $f_r$ and $e_r$ are obtained by lexicographically ordering matrices $M_r$, $F_r$, $E_r$.

For reference, eq (2) that they mention is the Fredholm integral we are trying to solve, which relates the measured data $M_r$ and the joint probability distribution $\mathcal{F}_r$ of the relaxation times (x, y):

$$M_r(\tau_1,\tau_2) = \int \int k_1(x, \tau_1)k_2(y, \tau_2)\mathcal{F}_r(x,y)dxdy + \epsilon_r(\tau_1,\tau_2) \quad r=1,...,R$$

According to the text, the shape of $M_r$ is $N_1\times N_2$ and $m_r$ is $(N_1N_2\times 1)$. I thought I could use a simple reshape, but looking at definitions of lexicographic, there seems to be a sort involved. So, which one of these two python snippets is right? Or is it something else entirely?

m_r = np.sort(M_r, axis=None)[:, np.newaxis] or

m_r = M_r.reshape((N1*N2, 1))?

Later, they also state that

The matrix $F_r$ is estimated by reordering $f_r$ into the matrix notation.

So, I could do a simple reshape back. But if I have to unsort $F_r$, which has the shape $N_x \times N_y$, I am lost.


问:
Efficient way to find eigenvalues of complex symmetric matrix with real off-diagonal elements
My goal is to find all eigenvalues (and eigenvectors) in a given range of magnitudes of a complex symmetric matrix with real off-diagonal elements (only diagonal elements are complex). Currently I'm using PARPACK for non-hermitian matrices and discard eigenvalues outside of the window of interest, but this method is more general than I need. Moreover, PARPACK only lets you use a maximum of matrix_size / NCV processors (NCV is proportional to the desired number of eigenvalues), but even within this limitation it does not scale efficiently. Do you know any algorithm or library that can take advantage of these restrictions and do better than PARPACK? Or at least scale better than PARPACK?


问:
2D diffusion equation using Finite Volume Method
i am working on an assignment problem:
Consider a two-dimensional rectangular plate of dimension L = 1 m in the x direction and H = 2 m in the y direction. The plate material has constant thermal conductivity. The steady-state temperature distribution within this plate is to be determined for the following imposed boundary conditions: (i) y = 0, T = 100 ̊C, (ii) x = 0, T = 0 ̊C, (iii) y = H, T = 0 ̊C, and (iv) x = L, T = 0 ̊C. Choose a uniform grid size of 0.05 m in both directions. Solve the problem using the point-by-point Gauss-Seidel iterative method. Experiment with the
initial guess and comment on the number of iterations required for convergence in each case. Clearly explain your convergence criterion for the iterations and how it is implemented. Plot the temperature contours as the output.

I have solved this question in python and i am getting following results:
When initial guess = 0, No of iterations = 350
Now when i am taking initial guess less than 10 i get less no. of iterations but when i take initial guess to be larger than 10 i get larger number of iterations.
Can anyone explain. Please it is very urgent and important.


问:
Spectral solver on em-pic
I'm recently studying for the spectral solver to implement EM-PIC code. I read an article and have some questions.

Many PIC codes uses spectral solver to overcome numerical artifacts on FDTD.

In the high-order approximation of spatial derivative in PSTD and PSATD scheme,
$$
\partial_{x}E_{x}|^{n}_{i,j,k} = \Sigma_{l=0}^{p/2-1} c_{l,p}\frac{E_x|^{n}_{i+1/2+l,j,k}-E_x|^{n}_{i-1/2-l,j,k}}{\Delta x}
$$
after Fourier transform due to the efficiency of calculation,
$$
F[\partial_{x}E_{x}] = i[k_x]_p \hat{E_x}
\\
[k_x]_p = \Sigma_{l=0}^{p/2-1} c_{l,p}\frac{e^{ik_x(l+1/2)\Delta x}-e^{-ik_x(l+1/2)\Delta x}}{\Delta x}
$$
where $F$ represent the Fourier transform and $\hat{E_x} = F[E_x]$.

The term $[k_x]_p$ goes to $k_x$ as p goes to infinity. But commonly use values for $p$ in practice 32 or 64 (as the article said).

Is there any reason to use the $p$ value of 32 or 64 in practice?

If I let $p=\infty$, I don't need to use the above high-order approximation formula. Then, it becomes simpler and the computational cost will be reduced I think.

Title: Spectral solver on em-pic