问:
Inverse Newton Method for optimization: is this the correct algorithm?
I am trying to implement the algorithm in this article. I have already asked a question before about it here, and I am trying to figure out what I am doing wrong. This time, it's this section of the article that I am struggling with (bold part is mine).

[...] it can be shown that the vector $c_r$ that satisfied (29) can be estimated by minimizing the function $$\chi(c_r) = \frac{1}{2} c_r'[G(c_r) + \alpha I]c_r - c_r'\tilde{m}_r$$ where $I \in \mathcal{R}^{(s_1s_2)\times(s_1s_2)}$ is an identity matrix. The optimization for $c_r$ can be performed by using the inverse Newton method since the first and second derivative of $\chi(c_r)$ are available in closed form as $$\nabla \chi(c_r) = (G(c_r) + \alpha I)c_r - \tilde{m}_r$$ $$\nabla \nabla \chi(c_r) = G(c_r) + \alpha I$$

I looked for this inverse Newton method and didn't find anything with that exact name (not really my area of expertise). However, I did find Newton's method for higher dimensions in Wikipedia. Is that the algorithm that I should use? If not, could you give me some pointer to what I could use? I implemented it in Matlab, but my final result is still not even close to what I should have. I simulated some data, transformed it, added some noise, then tried to invert it back. I should get something like what I started with, but it's totally different and nonsensical. So I am trying to rule out possible failure areas. Here is the matlab code for all the functions, and the formula for $G(c_r)$. If anyone could double check my work, it would mean a lot. I think it looks much more difficult than it really is, but I can't wrap my head around what I'm doing wrong.

Formula for $G(c_r)$:

$$
G(c_r) =
\tilde{K}_0 \begin{bmatrix}
H(\tilde{K}'_{0,1}c_r) & 0 & \cdots & 0 \\
0 & H(\tilde{K}'_{0,2}c_r) & \cdots & 0 \\
\vdots & \vdots &  & \vdots \\
0 & 0 & \cdots & H(\tilde{K}'_{0,N_xN_y}c_r)
\end{bmatrix} \tilde{K}_0'
$$

Here, $\tilde{K}_{0,i}$ is the ith column of $\tilde{K}_0$ and $H$ is the Heavyside function.

(Inverse?) Newton method: $x_{k+1} = x_k - \gamma [f''(x_k)]^{-1} f' (x_k)$

cr = ones(s1*s2, 1);  % Initial Guess
alpha = 1;  % Initial Guess
gamma = 1; % Step size.
prev_chi_val = 1E8;  % Just to check when it's converging
epsilon = 0.001;
maxiter = 50;

for num = 1:maxiter
   G_calc = G(cr, K0til, Nx, Ny);  % Needed for all 3 functions. Defined at the end
   chi_val = chi(cr, G_calc, alpha, mtil);  % Defined at the end
   D_chi = der1_chi(cr, G_calc, alpha, mtil); % Defined at the end
   DD_chi = der2_chi(G_calc, alpha); % Defined at the end
   cr = cr - gamma .* inv(DD_chi) * D_chi;  % Performing newton step
   fprintf('Iter %d: chi= %f, delta= %f\n', num, chi_val, prev_chi_val - chi_val);
   if abs(prev_chi_val - chi_val) < epsilon
       break
   end
   prev_chi_val = chi_val;  % Update goodness parameter
end
fprintf('Done\n')
```

<p>Function to calculate $G(c_r)$</p>

```function out = G(cr, K0til, Nx, Ny)
    temp = zeros(Nx * Ny, Nx * Ny);
    for ii = 1:Nx*Ny
        temp(ii, ii) = heaviside(K0til(1:end, ii)' * cr);
    end
    out = K0til * temp * K0til';
end
```

<p>Function to calculate $\chi(c_r)$</p>

```function out = chi(cr, G_calc, alpha, mtil)
    out = 0.5 * cr' * (G_calc + alpha .* eye(size(G_calc))) * cr - cr' * mtil;
end
```

<p>Function to calculate $\nabla \chi(c_r)$</p>

```function out = der1_chi(cr, G_calc, alpha, mtil)
    out = (G_calc + alpha .* eye(size(G_calc))) * cr - mtil;
end
```

<p>Function to calculate $\nabla \nabla\chi(c_r)$</p>

```function out = der2_chi(G_calc, alpha)
    out = G_calc + alpha .* eye(size(G_calc));
end
```

Sizes of the matrices. These are all either 2D matrices ($s_1s_2$ means multiplying these two numbers; they're the number of singular values in a previous step) or vectors, and one, $\chi$, is a number.

$\tilde{m}_r$: ($s_1s_2 \times 1$)

$c_r$ : ($s_1s_2 \times 1$)

$G$ : ($s_1s_2 \times s_1s_2$)

$\chi$ : ($1\times1$)

$\nabla \chi$ : ($s_1s_2 \times 1$)

$\nabla \nabla \chi$ :  ($s_1s_2 \times s_1s_2$)

$\tilde{K}_0$ : ($s_1s_2 \times N_xN_y$)|


问:
First-principles benchmark of CFD solver
A decade ago I saw a ~20 coupled state + random number problem that was used as a benchmark for meteorological CFD tools.  It had a fractal dimension (of the chaotic attractor) around 2.3, I think.  I saw matlab code for it.

I'm looking for the method now, and can't remember what it was called.  To be clear it was entirely synthetic.  It was not flow in 2d or 3d.  It wasn't a cylinder, cavity, or other plausible mechanical system.  It was something like 20 nonlinear coupled states, and a random number generator.

Can you help me find this method?  It would be able to differentiate between ensemble Kalman filter or 4dvar in terms of performance, I think.

ContentLicense: CC BY-SA 4.0


问:
Optimality conditions for optimal control: BVP - DAE
I am solving an optimal control problem of the form
$$
\min_u \qquad\int_0^T \langle u(t), u(t) \rangle \, \mathrm{d}t \\
s.t. \quad \dot{x} = \tilde{f}(x) + u, \quad x(0)=x_0 \\
\qquad  \tilde{\Phi}(x(T)) = 0\}
$$
reformulated as
$$
\min_{\xi(T)} \qquad \xi(T) \\
s.t. \quad \dot{z} = f(z,u), \quad z(0) = z_0 \\
\qquad  \Phi(z(T)) = 0
$$
where the dual variable on the dynamics is denoted by $\lambda$ and the dual variable on the terminal constraint is denoted by $\mu$. Also, we define $z := (x^{\top}, \xi)^{\top}$, $f(z,u) := (f(x,u)^{\top}, \langle u, u \rangle)^{\top}$, and $\Phi(z(T)) := \tilde{\Phi}(x(T))$.

I would like to use a BVP-DAE solver to solve directly the optimality conditions
$$
\dot{\lambda} = -\frac{df}{dz}, \quad \lambda(T) = \frac{d\xi(T)}{dz} + \mu \frac{d\Phi(z(T))}{dz} \\
\dot{z} = f(z,p), \quad z(0) = z_0 \\
\Phi(z(T)) = 0
$$
and was wondering what numerical software I could use. Ideally it would be a collocation-based DAE solver. I've looked into Sundials a bit, but was wondering what I should use. I've looked at a few previous questions but am looking for existing software. Thanks!


问:
Why this LJ molecular dynamics result doesn't converge?
I am doing a molecular dynamics simulation of Leonard Jones 6-12 potential. But instead of converging to a particular value. It always stays between -5.58 to -5.62. The standard value is -5.517. The problem is that however long it runs (I have tested up to 30000 steps) it does not converge to even 1 digit after the decimal point. It changes between -5.58 to -5.62 but never reaches -5.517. My timestep is 0.005 (as given in reference value). All my other parameters are the same as the value in reference. I am using a velocity-Verlet algorithm. The standard reference is available here.

I have used velocity rescaling and that's why kinetic energy remains almost constant, at least up to 2 digits after decimal always. Total energy change is quite small after a few steps, the change is mainly due to fluctuation in potential energy i.e. in the range -5.58 to -5.62.

The code to reproduce the error is below.

import numpy as np
import matplotlib.pyplot as plt
from numba import jit
import os
import sys

# All globals are compile time constants
# recompilation needed if you change this values
# Better way: hand a tuple of all needed vars to the functions
# params=(NSteps,deltat,temp,DumpFreq,epsilon,DIM,N,density,Rcutoff)

# Setting up the simulation
NSteps =100000 # Number of steps
deltat = 0.005 # Time step in reduced time units
temp   = 0.851# #Reduced temperature
DumpFreq = 100 # Save the position to file every DumpFreq steps
epsilon = 1.0 # LJ parameter for the energy between particles
DIM     =3
N       =500
density =0.776
Rcutoff =3

params=(NSteps,deltat,temp,DumpFreq,epsilon,DIM,N,density,Rcutoff)

#----------------------Function Definitions---------------------

#------------------Initialise Configuration--------

error_model=True
#If it is really required to search for division by zeros (additional cost)?

@jit(nopython=True,error_model="numpy")
def initialise_config(N,DIM,density):
    velocity = (np.random.randn(N,DIM)-0.5)

    # Set initial momentum to zero
    COM_V = np.sum(velocity)/N     #Center of mass velocity
    velocity = velocity - COM_V    # Fix any center-of-mass drift

    # Calculate initial kinetic energy
    k_energy=0
    for i in range (N):
        k_energy+=np.dot(velocity[i],velocity[i])
    vscale=np.sqrt(DIM*temp/k_energy)
    velocity*=vscale

    #wrong array initialization (use tuple)
    #Initialize with zeroes
    coords = np.zeros((N,DIM))

    # Get the cooresponding box size
    L = (N/density)**(1.0/DIM)

    """ Find the lowest perfect cube greater than or equal to the number of
     particles"""
    nCube = 2

    while (nCube**3 < N):
        nCube = nCube + 1

    # Assign particle positions
    ip=-1
    x=0
    y=0
    z=0

    for i in range(0,nCube):
        for j in range(0,nCube):
            for k in range(0,nCube):
                if(ip<N):
                    x=(i+0.5)*(L/nCube)
                    y=(j+0.5)*(L/nCube)
                    z=(k+0.5)*(L/nCube)
                    coords[ip]=np.array([x,y,z])
                    ip=ip+1
                else:
                    break
    return (coords/L),velocity,L

@jit(nopython=True)
def wrap(pos,L):
    '''Apply perodic boundary conditions.'''

    #correct array indexing
    for i in range (len(pos)):
        for k in range(DIM):
                if (pos[i,k]>0.5):
                    pos[i,k]=pos[i,k]-1
                if (pos[i,k]<-0.5):
                    pos[i,k]=pos[i,k]+1

    return (pos)

@jit(nopython=True,error_model="numpy")
def LJ_Forces(pos,acc,epsilon,L,DIM,N):
    # Compute forces on positions using the Lennard-Jones potential
    # Uses double nested loop which is slow O(N^2) time unsuitable for large systems
    Sij = np.zeros(DIM) # Box scaled units
    Rij = np.zeros(DIM) # Real space units

    #Set all variables to zero
    ene_pot = np.zeros(N)
    acc = acc*0
    virial=0.0

    # Loop over all pairs of particles
    for i in range(N-1):
        for j in range(i+1,N): #i+1 to N ensures we do not double count
            Sij = pos[i]-pos[j] # Distance in box scaled units
            for l in range(DIM): # Periodic interactions
                if (np.abs(Sij[l])>0.5):
                    Sij[l] = Sij[l] - np.copysign(1.0,Sij[l]) # If distance is greater than 0.5  (scaled units) then subtract 0.5 to find periodic interaction distance.

            Rij   = L*Sij # Scale the box to the real units in this case reduced LJ units
            Rsqij = np.dot(Rij,Rij) # Calculate the square of the distance

            if(Rsqij < Rcutoff**2):
                # Calculate LJ potential inside cutoff
                # We calculate parts of the LJ potential at a time to improve the efficieny of the computation (most important for compiled code)
                rm2      = 1.0/Rsqij # 1/r^2
                rm6      = rm2**3
                forcefact=(rm2**4)*(rm6-0.5) # 1/r^6
                phi      =4*(rm6**2-rm6)

                ene_pot[i] = ene_pot[i]+0.5*phi # Accumulate energy

                ene_pot[j] = ene_pot[j]+0.5*phi # Accumulate energy

                virial     = virial-forcefact*Rsqij # Virial is needed to calculate the pressure
                acc[i]     = acc[i]+forcefact*Sij # Accumulate forces
                acc[j]     = acc[j]-forcefact*Sij # (Fji=-Fij)
    #If you want to get get the best performance, sum directly in the loop intead of
    #summing at the end np.sum(ene_pot)
    return 48*acc, np.sum(ene_pot)/N, -virial/DIM # return the acceleration vector, potential energy and virial coefficient

@jit(nopython=True,error_model="numpy")
def Calculate_Temperature(vel,L,DIM,N):

    ene_kin = 0.0

    for i in range(N):
        real_vel = L*vel[i]
        ene_kin = ene_kin + 0.5*np.dot(real_vel,real_vel)

    ene_kin_aver = 1.0*ene_kin/N
    temperature = 2.0*ene_kin_aver/DIM

    return ene_kin_aver,temperature

# Main MD loop
@jit(nopython=True,error_model="numpy")
def main(params):
    NSteps,deltat,temp,DumpFreq,epsilon,DIM,N,density,Rcutoff=params
    # Vectors to store parameter values at each step

    ene_kin_aver = np.ones(NSteps)
    ene_pot_aver = np.ones(NSteps)
    temperature = np.ones(NSteps)
    virial = np.ones(NSteps)
    pressure = np.ones(NSteps)

    pos,vel,L = initialise_config(N,DIM,density)
    acc = (np.random.randn(N,DIM)-0.5)
    volume=L**3

    # Open file which we will save the outputs to
    # Unsupported operations have to be in an objectmode block
    # or simply write the outputs at the end in a pure Python function
    """
    if os.path.exists('energy2'):
        os.remove('energy2')
    f = open('traj.xyz', 'w')
    """
    for k in range(NSteps):

        # Refold positions according to periodic boundary conditions
        pos=wrap(pos,L)

        # r(t+dt) modify positions according to velocity and acceleration
        pos = pos + deltat*vel + 0.5*(deltat**2.0)*acc # Step 1

        # Calculate temperature
        ene_kin_aver[k],temperature[k] = Calculate_Temperature(vel,L,DIM,N)

        # Rescale velocities and take half step
        chi = np.sqrt(temp/temperature[k])
        vel = chi*vel + 0.5*deltat*acc # v(t+dt/2) Step 2

        # Compute forces a(t+dt),ene_pot,virial
        acc, ene_pot_aver[k], virial[k] = LJ_Forces(pos,acc,epsilon,L,DIM,N) # Step 3

        # Complete the velocity step
        vel = vel + 0.5*deltat*acc # v(t+dt/2) Step 4

        # Calculate temperature
        ene_kin_aver[k],temperature[k] = Calculate_Temperature(vel,L,DIM,N)

        # Calculate pressure
        pressure[k]= density*temperature[k] + virial[k]/volume

        # Print output to file every DumpFreq number of steps
        """
        if(k%DumpFreq==0): # The % symbol is the modulus so if the Step is a whole multiple of DumpFreq then print the values

            f.write("%s\n" %(N)) # Write the number of particles to file
            # Write all of the quantities at this step to the file
            f.write("Energy %s, Temperature %.5f\n" %(ene_kin_aver[k]+ene_pot_aver[k],temperature[k]))
            for n in range(N): # Write the positions to file
                f.write("X"+" ")
                for l in range(DIM):
                    f.write(str(pos[n][l]*L)+" ")
                f.write("\n")

        #Simple prints without formating are supported
        """
        if (k%150==0):
            print("step: ",k,"KE: ",ene_kin_aver[k],"PE :",ene_pot_aver[k],"\n Total Energy: ",ene_kin_aver[k]+ene_pot_aver[k])
            #print("\rStep: {0} KE: {1}   PE: {2} Energy:  {3}".format(k, ene_kin_aver[k], ene_pot_aver[k],ene_kin_aver[k]+ene_pot_aver[k]))
            #sys.stdout.write("\rStep: {0} KE: {1}   PE: {2} Energy:  {3}".format(k, ene_kin_aver[k], ene_pot_aver[k],ene_kin_aver[k]+ene_pot_aver[k]))
            #sys.stdout.flush()

    return ene_kin_aver, ene_pot_aver, temperature, pressure, pos

#------------------------------------------------------
ene_kin_aver, ene_pot_aver, temperature, pressure, pos = main(params)

# Plot all of the quantities
def plot():
    plt.figure(figsize=[7,12])
    plt.rc('xtick', labelsize=15)
    plt.rc('ytick', labelsize=15)
    plt.subplot(4, 1, 1)
    plt.plot(ene_kin_aver,'k-')
    plt.ylabel(r"$E_{K}", fontsize=20)
plt.subplot(4, 1, 2)
plt.plot(ene_pot_aver,'k-')
plt.ylabel(r"$E_{P}$", fontsize=20)
plt.subplot(4, 1, 3)
plt.plot(temperature,'k-')
plt.ylabel(r"$T$", fontsize=20)
plt.subplot(4, 1, 4)
plt.plot(pressure,'k-')
plt.ylabel(r"$P$", fontsize=20)
    plt.show()

plot()

```

I have tried with different timesteps but no result. I think I have used reduced units everywhere. Epsilon is taken as 1 and sigma is never used (always taken as 1). Box size is the length of the cube grid required (L) to make a system of required reduced density (density as a parameter in the code). It is calculated by the system and no. of particles and density (reduced) is taken as input.

I would also like to know how to calculate the radial distribution function using this code. I know its definition but can't understand how to implement it in code. As far as I can understand, my origin will be (0,0,0) as I have at first in initialise_config function made all the coordinates with respect to the center of mass. But how to implement a function to efficiently calculate the number of particles at a distance $r$? Making a histogram corresponding to the radial distance of each LJ particle from origin may work but I am not sure how it will produce the negative parts in the $\rho$ axis in the following plot.

.

Any help in implementing this will be much appreciated. I don't have any clue on how to sort these 2 problems.


问:
Petsc matrix with DM objects takes too much memory
I have been learning petsc recently and made a simple solver for stokes equation on 3D staggered grid.
My code is running on server with 200GB ram, but if I run a solver on 256^3 grid, memory used by the solver is over 200GB.

Checking with the -log_view flag, I've figured out most of the memories were used by the matrix created by

```<code>DMCreateMatrix(Mat &A)
</code>```

I understand that the DM structure helps the MPI communications and therefore may need extra memories, but it seems too much in my opinion. Is there any way to reduce the memory usage by Matrix created from DM objects?


问:
Spectral Clustering by Andrew Ng paper: theorem proof question
I recently read the paper "On Spectral Clustering: Analysis and an Algorithm" by Ng et al.

Much of the paper centers on Theorem 2 and equation 8.
To me, it appears there is no given or referenced proof of this Theorem.
I am hoping that someone can point me in the right direction to understanding the derivation of this result.

Perhaps a longer version of the paper (I am unable to find one so far), or is the result 'obvious' and I am missing something?
Any help is greatly appreciated!


问:
What is the state-of-the-art in parallel sparse matrix and dense vector multiplication?
For sure, there has been many highly optimized library on this. But I am working on a matrix-free context since the problem size does not allow explicit storage of sparse matrix elements. I'd love have a general knowledge on how people handle this sort of problem, especially on how to alleviate the memory bottleneck.

=========================update========================

As suggested by other commentators, I shall further elaborate on my situation.

My current implementation is below, which is just the basic one showing up in various textbooks.

```<code>void SpMVMtply(vector<det>& basis, ColVec& x, ColVec& b)
{

    /*
         "basis" is the orthonormal basis that spans the physical space
         "x" ,"b" are just those in Ax = b
    */
    int len = x.size();

    #pragma omp parallel
    {
        #pragma omp for nowait
        for (int i=0; i<len; i++)
        {
            vector<det> pool_det;
            vector<cdouble> pool_amp;

            /*
               this function generates the non-zero off diagonal elements
               with respect to the given basis vector. pool_det is the
               collection of generated basis vectors, and pool_map collection
               of the amplitude.
            */

            OffDiagGen(basis[i], pool_det, pool_amp);

            // this function generates the diagonal element
            cdouble sum = x(i)* DiagCal(basis[i]);

            /*
                Below, basis_pos is an unordered_map with basis vector as key
                and that basis vector's position in "basis" array as value.

            */
            for (int j=0; j< pool_det.size(); j++)
                sum += x(basis_pos[pool_det[j]]) * conj(pool_amp[j]);

            b(i) = sum ;
        }
    }
}
</code>```

The size of vector could reach 1 billion, so the explicit storage of sparse matrix is impossible. A rough profiling shows that the non-sequential accessing of memory in x(basis_pos[pool_det[j]]) is the most time-consuming part. On a Xeon Gold 6242 with 12 cores in use, it takes 6.2 seconds to such part for a vector with 3 million elements, while it takes only 1.5 seconds when that part is replaced by a constant. That's why I call it a memory bottleneck. My main concern is that if such bottleneck is inevitable, since even when the sparse matrix is explicitly stored, such non-sequential memory access also shows up. Based on my knowledge I haven't figured out a way to avoid it.


问:
Does BFGS preserve the bandedness of the inverse hessian?
In the BFGS method we perform iterations by calculating an approximation $\boldsymbol{H}_k$ to the inverse Hessian $\boldsymbol{H}$ of the objective function. This method belongs to a family of methods in which $\boldsymbol{H}_k$ is guaranteed to preserve the positive definite nature of $\boldsymbol{H}$ for all $k$, provided the line search is exact.

Now, suppose the Hessian is a banded matrix, can we assure that $\boldsymbol{H}_k$ will also be banded? If not - is there any other quasi-Newton method that guarantees this?


问:
Tackling multiscale problem in numerical simulation
In a dusty plasma system there are more than one component with different masses, i.e, electrons, ions,neutrals and dust grains. Accordingly, there are more than one temporal and spatial scales involved. I am aware that such a system is difficult to simulate, because, due to very different scales involved it may require a very large no. sample to be collected and this may not be feasible  with the available computational resources. But, such multicomponent simulations are done, for example, using Particle In Cell method. How this problem of multiple scales is overcome in PIC simulation ?</p>


问:
Fastest way to solve linear programming with 6 complex inequalities and 5 nonnegative variables
I have a program where I need to solve a linear programming problem in a fast loop. The language I'm using is Java and any kind of bindings to other languages are not acceptable. Libraries might be acceptable, given permissive license, but I generally prefer own code to libraries.

The problem is very simple. It has five variables, six inequalities of the type expression >= 0 and all of the five variables obviously need to be non-negative too.

Specifically, the problem is to minimize e, given:

```<code>a + b - c*A - B >= 0
C - c >= 0
D + d - c*E - b >= 0
a + b + c - F >= 0
e - a*G >= 0
e - d*H >= 0

a >= 0
b >= 0
c >= 0
d >= 0
e >= 0
</code>```

where a, b, c, d and e are variables and A ... H are constants. Some of the constants (A, B, C, D, and F) change rapidly and some (E, G and H) stay constant for the duration of the application execution.

The problem is that I have to find an optimal solution preferably over 10 million times per second in a single modern CPU core (yes, I'm using multiple cores, each of them solving the same linear programming problem repeatedly; the aggregate throughput needs to be in excess of 100 million solutions per second).

What would be the fastest approach for solving such a small linear programming problem? I see that there are quite many algorithms capable of solving even complex problems, but are they beneficial for such simple problems?

I see that there are 462 ways to pick 5 equations out of 11, so the approach of solving a different linear system of equations and checking whether the non-picked conditions hold, would not be fastest, right?


问:
How to determine the finite difference coefficient matrix in 2D with periodic BC?
I'm solving a PDE in matlab using ode15s, and since the spatial dimension is 2, and number of variables grow large very quickly, I need to supply the structure of the jacobian (called JPattern in Matlab) by indicating the position of nonzero elements of the jacobian. The jacobian 'pattern' for each point in 2D (indicated as c) is as follow:

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

To make myself more clear, if it is just the laplacian, then the pattern is

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

It is easy to turn the laplacian to the finite difference matrix, which is just a banded  matrix with five diagonals grouped together, with another five diagonals separated by the number of rows (or columns depending on how the points are ordered). With periodic B.C., there are extra points in the 'corners' of each sub matrix inside. What about the one for my problem? Are there ways to figure out the finite difference matrix given a pattern like the above?


问:
How is System.Decimal represented in memory bits?
I am trying to look at how different floating points are stored in memory.

Firstly I looked at the System.Double (accessible by keyword Double in vb.net) which I think I understand. It is stored as follows:

$\pm (I+a)\times 2^b$

where

$\pm$ requires one bit.

$a \in [0,\sum_{k=1}^{52}2^{-k}]$ and consumes $52$ bits.

$b \in [-2^{10}+2, 2^{10}-1]$ and consumes the remaining $11$ bits.

So in total, this representation requires $64$ bits.

Note $11$ bits can represent a total of $2^{11}=2048$ distinct numbers while the set of possibilities of $b$ are $|[-2^{10}+2, 2^{10}-1]|=2046$. The remaining two binary representations of $All[0]:="00000000000"$ and $All[1]:="11111111111"$ are omitted to cover special cases of $\pm 0.$, $\pm \infty$ and different types of not a number, $NaN$.

$I$ is $1$ except for special cases of $b=All[0]$ where it is $0$ so doesn't consume any additional bit.

Like if $b=All[0]$ and $a$ also has all its bits $0$ then the resulting number represent $\pm0$.

Similarly, if $b=All[1]$ and $a$ has all its bits $0$ then the resulting number represents $\pm\infty$.

Similarly, if $b=All[1]$ and $a$ has some particular bits 1 then the resulting numbers represents various types of $NaN$.

Similarly, if $b=All[0]$ and $a$ has some particular bits 1 then the resulting $b$ is replaced by its minimum value $-2^{10}+2$ and since $I=0$, in this case, it allows for more closer numbers to $0.$ be represented (Subnormal numbers).

This is how I understand the entire Double (64bit floating point). Now when trying to understand the System.Decimal in the same spirit I am unable to decipher all the details.

I understand that decimal has the following structure:

$\pm a \times {10}^{-b}$

where $a$ is a 96 bit number so $a \in [0,2^{96}-1]$ while $b \in [0,28]$ is a scaling factor that decide where to put the decimal point. Thus it can store up to $28$ decimal places.

The documentation states that this number requires $128$ bits. Which means using the 96 bit for the $a$ and one bit for the sign, this scaling factor is using the remaining $128-96-1=31$ bits. Can somebody explain how this scaling factor is using the $31$ bits? It's presumingly has a base two representation of the scaling factor or something. Can someone who knows explain how the $31$ bits are used in this Decimal data type?

https://docs.microsoft.com/en-us/dotnet/visual-basic/language-reference/data-types/decimal-data-type


问:
How to solve odd-order differential equations in FEM? Petrov-Galerkin?
I've recently learned about using weighted residuals with the Galerkin method to numerically approximate even-order differential equations (for linear elements, I'm still a beginner). It seems for odd-order differential equations you cannot use integration by parts to reduce the order of the third-order differential operator.

I've read that the Petrov-Galerkin method can approximate these solutions, but I cannot find cogent material on it or examples of it in action. Can someone please explain?


问:
Lattice spring models vs. finite element models
I am a beginning graduate student in the field of continuum mechanics. It is my understanding that most problems in this field are numerically solved via finite element methods (FEM). However, I have also seen  papers where people are using lattice spring models (LSM) as their chosen numerical method.

I am not able to understand the difference between the two methods at a fundamental level as they both use finite elements for discretization. Can somebody please help me compare/contrast the two methods? I will also appreciate if somebody can point me to places where such comparisons might have been made. I tried to search for it, but couldn't find anything of much relevance on the internet.


问:
Choice of method for Stochastic Differential Equation Solution
I have to simulate the following differential equation:

$$\dot{x}=p$$

$$\dot{p}=-x-x^3-\gamma p+F\cos\omega t+ \delta F(t)$$

where

$$\langle F(t)F(t')\rangle=2\gamma k_BT\delta(t-t')$$

$$\langle F(t)\rangle=0$$

$F=0.006$, $\omega=1.016$, $k_BT=0.008$, $\gamma=0.001$

I want to simulate for a long time of $T=800$ and want a decently accurate solution. There are many methods to choose from like Euler-Maruyama, Milstein, Runge-Kutta, etc. I was wondering which method to choose.

What are the qualifying criteria between all these methods?

When is the improved accuracy useless due to the presence of a stochastic term?


问:
Modelling of Stefan Maxwell equation
I am trying to solve Maxwell Stefan's equation over a membrane to get the transient mole fraction distribution over the membrane thickness 'z'. But somehow I am not able to code it using ODE45, more likely I am not able to write the system to solve using ODE45. It will be really great if someone can help me with the primary syntaxes and function. The equation I am trying to solve is

$\frac{dy_{H_2}}{dz}={\frac{1}{C \times D_{H_2,{H_2O}}}}\left[y_{H_2}(N_{H_2}+N_{H_2O})-N_{H_2}\right]$

where C is concentration, $D_{H_2,H_2O}$ is the binary diffusion coefficient, N is molar flux and y is mole fraction

Thank you in advance.

Title: Modelling of Stefan Maxwell equation


问:
Boundary Conditions for Continuum Mechanics
If I'm given some sort of shape and there is a displacement given but there are no external forces acting on the body, do I still need to write down the traction free boundary conditions or will there only be displacement boundary conditions?</p>


问:
von Neumann analysis: computation of maximum value of amplification factor
In this question I adress the stability property of a numerical scheme (a FDM described below) to the problem

$$
\begin{align}
u_{t} &= \alpha u_{xx}+\beta u_{xxxx}, & &x\in\mathbb{R},\;t>0, \\[3mm]
u(x,0) &= \sin x, & &x\in\mathbb{R}, \\[3mm]
u(x+2\pi,t)&=u(x,t), & &x\in\mathbb{R},\;t>0.
\end{align}
$$

where $\beta \leq0$ and $\alpha >0$. It can be shown for this setting that the problem is well-posed. I discretize the problem by a central finite difference in space and forward Euler in time:

$$
\begin{eqnarray}
    u_{j}^{n+1} &=& \underbrace{\frac{\Delta t\beta}{\Delta x^{4}}}_{\equiv b_{2}}u_{j+2}(t) + \underbrace{\left( \frac{\Delta t\alpha}{\Delta x^{2}} - \frac{4\Delta t\beta}{\Delta x^{4}} \right)}_{\equiv b_{1}} u_{j+1}(t) \\
    &+& \underbrace{\left(1 + \frac{\Delta t\beta}{\Delta x^{4}} -\frac{\Delta t\alpha}{\Delta x^{2}}\right)}_{\equiv b_{0}}u_{j}(t) + \underbrace{\left( \frac{\Delta t\alpha}{\Delta x^{2}} - \frac{4\Delta t\beta}{\Delta x^{4}} \right)}_{\equiv b_{-1}}u_{j-1}(t) \nonumber \\
    &+& \underbrace{\frac{\Delta t\beta}{\Delta x^{4}}}_{\equiv b_{-2}} u_{j-2}(t) \nonumber.
\end{eqnarray}
$$

To determine stability of the scheme we want

$$
\max_{k}\vert g_{k}\vert = \max_{k}\vert \sum_{l=-2}^{2}b_{l}\exp(\text{i}kl\Delta x)\vert \leq 1+\theta\Delta t,\quad \text{for some $\theta>0$}.
$$

My problem is that I need a little help with the algebra, which gets messy. I obtain

$$
\begin{align}
&\vert g_{k}\vert = \left\vert b_{2}e^{2\text{i}k\Delta x} + b_{1}e^{\text{i}k\Delta x} + b_{0}e^{0\text{i}k\Delta x} + b_{-1}e^{-\text{i}k\Delta x} + b_{-2}e^{-2\text{i}k\Delta x} \right\vert = \\[3mm]
&= \left\vert 2b_{2}\cos(2k\Delta x) + 2b_{1}\cos(k\Delta x) + b_{0} \right\vert.
\end{align}
$$

From here I am wondering, is there anymore simplifications I can make before taking an upper bound to this result? Can I simply continue by making

$$
\begin{align}
&\left\vert 2b_{2}\cos(2k\Delta x) + 2b_{1}\cos(k\Delta x) + b_{0} \right\vert \leq 2\vert b_{2}\vert\vert\cos(2k\Delta x)\vert + 2b_{1}\vert\cos(k\Delta x)\vert +\vert b_{0}\vert \leq \\[3mm]
&\leq -2\frac{\Delta t\beta}{\Delta x^{4}}\vert\cos(2k\Delta x)\vert +2\left( \frac{\Delta t\alpha}{\Delta x^{2}} - \frac{4\Delta t\beta}{\Delta x^{4}} \right)\vert\cos(k\Delta x)\vert + 1 -\frac{\Delta t\beta}{\Delta x^{4}}+\frac{\Delta t\alpha}{\Delta x^{2}}
\end{align}
$$

and obtain a reasonably tight upper bound?


问:
What is the algorithm to convert an adjacency matrix to a block diagonal structure?
I'm a little at loss as to what my supervisor meant by the "hierarchical block decomposition" of a matrix, but the goal is to put a sparse symmetric adjacency matrix into a block diagonal structure to highlight the structure of the clusters present in the data.

(The data captures whether the two proteins are adjacent in a molecule: if they are, the corresponding adjacency matrix has one on their row/column, otherwise zero. I sum a few of these adjacency matrices(one for each unique molecule) and divide the result by the number of molecules(to simulate the mean).

I'm wondering how to get to the block diagonal structure from there. People over here suggest scipy's reverse CM and seaborn's heatmap should be able to do that ostensibly, but i had no luck with either so far.
Should i go with SVD? What am i missng?
This is what I got. Thank you kindly!

```<code>def construct_adjacency_matrix(clusters=dict, nomenclature_namespace=dict):

   nbrpairs  = tree2tuplearr(clusters['nbrtree'])
    dim       = len(nomenclature_namespace.items())
    keys      = list(nomenclature_namespace.keys())
    substrate = np.zeros((dim, dim))

    for pair in nbrpairs:
        if (pair[0] not in keys or pair[1] not in keys):
            pass
        else:
            substrate[nomenclature_namespace[pair[0]],
                      nomenclature_namespace[pair[1]]] = 1
    return substrate

def get_simplemean(targetgroup=str):
    targetspath = './../clusterdata/targetgroups/{}/'.format(targetgroup)
    batch       = os.listdir(targetspath)
    adjmats     = [construct_adjacency_matrix(openjson(targetspath+x), namespace) for x in batch]
    # laplacians  = []
    reduced = functools.reduce(lambda x, y: np.add(x,y), adjmats)
    reduced = np.divide(reduced, len(adjmats))
    plt.matshow(reduced)
    seaborn.clustermap(reduced)
    plt.show()
</code>```

There is obviously certain structure in the data and in fact I already have the graphs of clusters, but how do i put into the form that resembles this below? I realize that this picture has quite a few more datapoints.


问:
The relation between PDE order and discretization order
In Jasak's Ph.D. thesis (2000), a notion is given about discretization of a transport equation:
For good accuracy, it is necessary for the order of the discretization to
be equal to or higher than the order of the equation that is being discretized.

Question:

Why?


问:
Optimizing vectors with equal elements
I am trying to distribute power across different devices, so that the sum is as equal as possible to the power setpoint. At the same time, the sum of power per phase must not exceed the power of the fuse and must also be between min and max power (depending on the device). So far I have reduce the problem to the following:
$$
\sum_i^n \sum^z P_i= P_{setpoint}
$$
$$
\sum_i^n P_{x i}\leq P_{xfuse}
$$
$$
\sum_i^n P_{y i} \leq P_{yfuse}
$$
$$
\sum_i^n P_{z i} \leq P_{zfuse}
$$
$$
P_{min} \leq P_i \leq P_{max}
$$
Where I get stuck is at defining devices, since some use one phase, others three, so their power vector would look like: $P_i = (0, 0, P_z)$ or $P_i = (P_x, 0, 0)$ (for one phase devices) or $P_j = (P_x, P_y, P_z)$ but $P_x = P_y = P_z$ (for three phase devices). Also, power per device is defined as $ P_{device} = Px + Py + Pz$.

I think this could be easily solved with linear optimization, but I don't know how to enforce the eqality of the elements for three phase devices. Also, how would I ensure that the sum of power of all devices is as equal as possible to the setpoint?

I searched for existing answers, but I haven't found anything helpful(I might be using the wrong terms).

Edit

As said in the comments, my optimization goal was not defined well. What I really want is to  maximize the power of each device with above constraints. In other words, I want all devices to still use as much power as possible, but their power still sums to $P_{setpoint}$.

Title: Optimizing vectors with equal elements


问:
How does Block-diagonal Schur factorization in Matlab works?
Can anyone explain how does the bdschur command of Matlab works? And how should the CONDMAX and BLKS be specified? It would be better if there is an associated example.

Source: https://www.mathworks.com/help/control/ref/bdschur.html


问:
Poisson equation with FFT and normalization
I'm trying to understand how to solve Poisson equation with FFT. Say, if we have the simplest periodic example
$$u_{xx}=-4\pi\cos(x)$$

The solution then should be
$$u=4\pi\cos(x)$$
I really get the cosine, but with a different magnitude. So what I`m missing, and how to normalize solution with FFT properly?

My code is following:

```<code>clear all;
Xmin=-8*pi;
Xmax=8*pi;
Nsteps=800;
h=(Xmax-Xmin)/Nsteps;
x=[Xmin:h:Xmax]+10e-6;
f=-4*pi*cos(x);

F=fft(f);
k=2*pi./x;
FI=-F./k.^2;
fi=real(ifft(FI));

fiReal=4*pi*cos(x);

%%
figure,
subplot(1, 2, 1);
plot(x,fi) ; title('FFT solution')
subplot(1, 2, 2);
plot(x,fiReal); title('Analytical solution')
</code>```


问:
Is there a subfield within computational science research that's done on pencil and paper?
Is there an area of computational science research that can be done on pencil and paper (with results written up for a journal format later on)?

I'm wondering if there is abstract proof-based linear algebra or abstract algebraic objects being developed for computational science.  Seems like it would be really fun if such an area of research existed.


问:
Numerical method for harmonic oscillator with jumping constant
Let $k_1 \neq k_2$ be positive reals, $t_0 > 0$ and consider the following Cauchy problem in $[0,+\infty)$:
\begin{cases}
y(t) + k(t)y''(t) = 0 \newline
y(0) = 1/\sqrt{k_1} \newline
y'(0) = 0,
\end{cases}
where
\begin{equation}
k(t) =
\begin{cases}
k_1 \hspace{1 cm} 0 \leq t \leq t_0 \newline
k_2 \hspace{1 cm} t > t_0
\end{cases}
\end{equation}

It is clear that there exists a unique $C^1$ solution $y$ in $[0,+\infty)$, obtained by gluing two smooth solutions of the equation in $[0,t_0]$ and $[t_0,+\infty)$. An explicit computation gives:
\begin{equation}
y(t) =
\begin{cases}
\dfrac{1}{\sqrt{k_1}} \mathrm{cos} \bigg(\dfrac{t}{\sqrt{k_1}} \bigg) \hspace{1 cm} 0 \leq t \leq t_0; \newline
A \ \mathrm{cos}\bigg( \dfrac{t}{\sqrt{k_2}} + \phi \bigg) \hspace{1 cm} t > t_0,
\end{cases}
\end{equation}
where $A$ and $\phi$ are obtained imposing the $C^1$ condition. With further computations one discovers that the amplitude is constant, i.e. $A = 1/\sqrt{k_1}$, if and only if $t_0 = n \pi \sqrt{k_1}$ for some $n \in \mathbb N$.

I am searching for a good numerical method for the problem in the constant-amplitude case. I have tried several classical schemes, but I discovered that when $k_2 >> k_1$ the second amplitude also grows bigger. Any suggestion?


问:
Bifurcation points on homotopy path by numerical continuation?
I am trying to implement an algorithm that finds (possibly) all solutions of a system of nonlinear polynomial equations $$F(X) = 0$$ I thought about using the (convex) homotopy $$H(X,T) = TF(X) + (1-T)G(X_0)$$ with $X_0$ a known solution to $G(X_0)=0$. By increasing $T$ iteratively from $0$ to $1$, the sought solution of $F(X)$ will be found.

During the hike on that path, I expect to encounter points of different bifurcation types. I hoped to use some of the many well-implemented numerical continuation algorithms that are available on the internet and would allow me to detect those points with their belonging paths.

I'm currently stuck at getting the numerical continuation running on the homotopy. To be more precise, the paths of solutions during numerical continuation is far from being a solution to $H(X,1)=0$. I do have a solid background in programming but only mediocre knowledge of mathematics. Therefore I hoped someone could review my approach and answer the following questions:

Is numerical continuation applicable to the homotopy method? I see that there are some differences compared to the continuation of an ODE (i.e. in the predictor/corrector scheme) but isn't the approach the same?

If not, can one detect the above-mentioned bifurcation points differently?


问:
Random Orthogonal Matrix Generation
This post is inspired by N. Higham post "What is Random Orthogonal matrix?".

In this post, N. Higham links to the two papers:

G. W. Stewart, The efficient generation of random orthogonal matrices with an application to condition estimators, SIAM J. Numer. Anal. 17(3), 403–409, 1980.

T. W. Anderson, I. Olkin and L. G. Underhill, Generation of random orthogonal matrices, SIAM J. Sci. Statist. Comput. 8(4), 625–629, 1987.

that illustrate $\mathcal O\left(\frac{4n^3}{3}\right)$ algorithms to construct a Haar distributed random orthogonal matrix based on Householder matrices and Givens rotations, respectively.

Is $\mathcal O(N^3)$ a proven complexity bound?

Is there any progress in the generation of random orthogonal matrices that are faster than $\mathcal O(N^3)$ or, at least, $\mathcal O\left(\frac{4n^3}{3}\right)$?

In this post on C++ generation of Random orthogonal matrices, I found a link to the paper preprint by F. Mezzardi which was an interesting read but was focused on other aspects of the random orthogonal matrix generations (different distributions, following of the generated matrix to the given distribution). In this question, I am more interested in speed and complexity.

Title: Random Orthogonal Matrix Generation


问:
Numerical flow visualization in 2D for a moving boundary,
I have a rigid body that moves according to a set of governing ODEs, and I'd like to numerically visualize the vortices that are shed by this object.  How could I proceed?

I've been reading up on the Navier-Stokes equation and have a better grasp of it now, and I've watched some tutorials on implementing a 2D solver.  But, I still have some confusion:

Do I prescribe the movement of the rigid body, when trying to solve the N-S equations?

It would seem strange to do so, because the movement of the rigid body came from dynamical solutions to the governing ODEs.


问:
Is there a published RQ decomposition column-major algorthm?
I am refactoring an existing algorithm where where a RQ decomposition (as opposed to the more common QR) would be rather useful.

Most common books on the subject (e.g. Golub and Van Loan) discuss QR extensively and provide algorithms for common implementations (e.g. Algorithm 5.2.1 for Householder or Algorithm 5.2.5 for Modified Gram Schmidt).

Unfortunately, it does not discuss RQ, except for a small note regarding LAPACK routines. Wikipedia suggests that the differences are trivial, where "RQ decomposition is Gram–Schmidt orthogonalization of rows of A, started from the last row".

Since the matrices I am working on are column-major, this is not especially attractive.

Are there well-known column-major implementations of RQ, especially using MGS or Householder?  Or even if there is some blocking algorithm that is efficient on column data?

Other details:
This particular problem does not require pivoting, and I only need R - I do not need Q.


问:
Fastest matrix library for Android (with GPU is possible)
I was working on an Android app that requires some linear algebra with matrices. The matrices will be somewhat medium-sized as they are not too small or too big. I was originally using jBlas because of the benchmarks I have seen online but after visiting multiple forums, they have all said that jBlas would be nowhere near as fast as what the benchmarks originally state. I have also heard about the Snapdragon Math Library to help improve the speed but I have not heard too much about it. From the general community, what linear-algebra libraries do you recommend that will be the fastest to use on Android? All help is greatly appreciated!</p>


问:
Implementation of Z^2 error estimator in Abaqus for adaptive mesh refinement
Currently, I am working on a remeshing routine for my simulations (Abaqus 6.14-1) using python scripts. The simulation deals with the Brinell indentation test and as the remeshing software I use Gmsh 4.5.2 (although the newest probably also works fine).

The main script is given by

```<code>main.py:
createIndentationTest.py
runJob.py                    #(only one small step)
odbExport.py
writeGeo.py                  #(.geo files are Gmsh readable geometry files)
remeshGeo.py                 #(here occurs the change from quadrangular elements to triangular ones)
importInAbaqus.py
runJob.py                    #(only one small step)
while totalTime-stepTime>0:
    odbExport.py
    writeGeo.py
    remeshGeo.py
    importInAbaqus.py
    runJob.py                #(only one small step)
    totalTime-=stepTime
</code>```

As an example, my last routine created 6 separate jobs, where after every step the mesh was exported to Gmsh, remeshed and imported back in Abaqus to continue the simulation.
The first job starts with quadrangular elements (element size of 0.25) for one step. Then the remaining 5 jobs use triangular elements with a custom density function, as can be seen in the following figure (The reason for switching from quadrangular elements to triangular ones during the remeshing process is that it was easier for me to work with in Gmsh).

Now, part of my thesis is to determine the error during the whole simulation and remeshing process and I have been given the task to use the $Z^2$ error estimator.
When compared to my reference simulation, where I ran a whole simulation from start to finish with quadrangular elements with the element size of 0.5, I can see there is a slight discrepancy when analyzing the stresses, as shown in Figure 2.

For now, I assume that the results on the left-hand side are more accurate, therefore I can assume that there is an error (or deviation) in my simulation results on the right-hand side where the remeshing was implemented.
For the past 2 weeks, I have studied some works on the $Z^2$ error estimator, but I fail to see how to implement this in my script. Therefore, I wanted to ask, if anyone knows how I can apply the error estimator in my python code.

I know that I have to export some data from the COORD fieldOutput and use these to solve this equation

$F(\underline{a})=\sum_{i=0}^n\,\left(\sigma_h(x_i,y_i)-\sigma_p^*(x_i,y_i)\right)^2=\sum_{i=0}^n\,\left(\sigma_h(x_i,y_i)-\underline{P}(x_i,y_i)\,\underline{a}\right)^2$

where $\underline{a}$ is a set of unknowns, $\sigma_h$ is the approximate solution to $\sigma$, $\underline{P}$ is a set of appropriate polynomial terms and $\sigma_p^*=\underline{P}\,\underline{a}$.

Another question I have stems from the fact that I change element types during the simulation process. To estimate the error, the $Z^2$ method uses so-called superconvergent points. But I have no clue if there will arise a problem when using these since I change my element type once during the whole simulation process.

I hope, someone can help me with this.


问:
Why does the naive barycentric hodgestar fail?
The discrete exterior calculus is defined first using circumcentric dual cells, because the primal and dual edges are orthogonal and thus the dual cells are convex. This leads to a diagonal hodge star for 1-forms that is defined as $\star_1 = \frac{|\star\sigma^1|}{|\sigma^1|}$ for a primal edge $\sigma^1$ and a dual edge $\star\sigma^1$, what has many computational advantages.

For barycentric cells, the operators are more complicated, constructions for $\star_1$ are found for example in this paper by Mohamed, Hirani and Samtaney. Their hodge star has no sparse inverse and while they suggest how a sparse inverse could be constructed, it would not satisfy the identity $\star_1\star_1^{-1} = -1$.

The reason for the more complicated constructions is, that the dual edges are no longer orthogonal and are kinked at their intersection with the primal edges.

(Image by Jiansong Deng)

I am now wondering, why you cannot construct a diagonal hodge star by intergrating over the kinked lines. The relevant identities are

$$
\langle\omega^1, \sigma^1\rangle = \langle\star\omega^1, \star\sigma^1\rangle \\
\langle\text{d}\omega^1, \sigma^2\rangle = \langle\omega^1, \partial\sigma^2\rangle \\
\langle\text{d}(\star\omega^1), \star\sigma^0\rangle = \langle(\star\omega^1), \partial(\star\sigma^0)\rangle
$$

For a circumcentric dual 2D mesh $\star\omega^1$ represents are straight line that connects the circumcenter of two triangles and intersects the center of the primal edge.

In a barycentric mesh it is the combination of two line segments. The first identity should hold by definition, when we construct the hodge star by integrating over both segments. And I think Stokes' theorem should hold when integrating over the whole boundary of a barycentric 2-cell.

Which part am I missing, why you cannot construct a diagonal definition that allows for non-convex boundaries composed of split edges with a kink in the middle?


问:
Solving a spectral system by reducing it to a single frequency - Feasability of approach?
I'm trying to solve the linear non-paraxial pulse propagation equation

$$\partial_z\hat{E}=ik_z\hat{E}$$
for a field defined as $$E=E(r, t, z)$$
The equation given above uses $$\hat{E}=\hat{E}(k_\perp, \omega, z)$$. Therefore, I have to use a FFT twice, once, for going from $r$ to $k_\perp$, and once for going from $t$ to $\omega$. This also means that I have to limit the used frequencies, after a FFT going from $t$ to $\omega$ will only support a certain amount of frequencies, depending on the resolution of $t$ and the limits of $t$, i.e. $t_{min}$ and $t_{max}$.

To properly resolve longer pulses at short wavelengths, this makes rather high resolutions necessary to capture all involved frequencies, and thereby increasing computation time. Therefore I had the idea to set $t$ to $0$, and $\omega$ to $\omega_0$ and thereby effectively reducing the 3d-system $(x, y, t)$ to a 2d-system $(x, y)$. When transforming $E(k_\perp, t, z)$ to $E(k_\perp, \omega, z)$ this would require a FFT over a single number (after $t$ only contains one value), and therefore $E(k_\perp, t, z)\equiv E(k_\perp, \omega, z)$.

Does that make sense, or am I introducing additional errors here? Or are there other methods to reduce the amount of variables resulting from the 3d-system, to speed up computation?


问:
How does the "Stable Fluids" algorithm by Jos Stam relate to the SIMPLE and PISO algorithms?
The "Stable Fluids" paper (*) by Jos Stam starts by acknowledging that "Our method cannot be found in the computational fluids literature, since it is custom made for computer graphics applications. The model would not be accurate enough for most engineering applications."

Nonetheless, it is based on the Navier-Stokes equations, solved in an Eulerian structured grid (well, with a semi-Lagrangian scheme for tracking properties that travel with the fluid).

When putting side-to-side his algorithm with the description of the SIMPLE/PISO algorithms as formulated in the CFD literature, I'm not able to find their common points or their differences.

I mean: can the algorithm by Jos Stam be understood as an approximation/modification of either SIMPLE or PISO? Or what should you modify in his algorithm so that it becomes SIMPLE/PISO?

Thanks a lot!

(*) Note: The link to the paper in PDF format has been obtained from the official Jos Stam site, so I believe it is legit to put it here.


问:
Why is scoring a separate problem even after docking is done?
With regards to protein-ligand interactions, we talk about two concepts, docking and scoring. I understand that docking means to find the best orientation of a ligand at the active site of a protein, and scoring is to compare several docked ligands to see which ligand best binds to the protein.

To quote this paper,

"While there is now a general consensus that most of the popular
  docking algorithms perform fairly well in generating sound poses,
  scoring functions most often fail to adequately evaluate the binding
  affinity."

What does it mean to be docked, if not being the pose with the best binding affinity? If not, how is docking done without calculation of binding affinity?


问:
Choosing the pivot for the rotation matrix in similarity transformation
I have arrived at an equation in the similarity transformation -

$M_r$ = $T. M_{r-1}. T^t$ ,where $T$ is the rotation matrix and $M_r$ ,$M_{r-1}$ are similar matrices. My aim is to find the rotation matrix $T$ with a pivot element $(a,b)$ such that all elements $m_{ij} \in M$ will be as follows :

$m_{ij} \neq 0$ , if $|i-j| = 1$

$m_{ij} \neq 0$ , if $i=j$ , i.e, it lies on the principal diagonal

I found the similarity transformation of a matrix in an application for Microwave filter design by Richard J. Cameron, "General coupling matrix synthesis methods for Chebyshev filtering functions" where they used similarity transformation using rotation matrix for annihilation of elements. But the concept of setting the pivot is not very clear to me.

Can someone explain how the pivot is set to eliminate a particular element of $M$, for instance?


问:
Solving nonlinear pendulum using Runge-Kutta 4 for smaller steps
I am trying to solve nonlinear pendulum using 4th order Runge-Kutta method for limits between a=0.0 to b=110 seconds and simulated the results to observe the pendulum movement. But when I increase the number of steps N in RK4 method, the pendulum isn't decaying, it is oscillating forever, but for smaller value of N, I can observe the decaying effects.

from math import sin, pi
import numpy as np
import matplotlib.pyplot as plt

# constants
G = 9.51        # acceleration due to gravity
L = 0.1         # length of the pendulum
T = 0.5         # time period of the pendulum
ANGLE = 90.0    # initial angle of the pendulum

def f(r, t):
    """Vectorized simultaneous ODEs functions."""
    theta = r[0]
    omega = r[1]
    ftheta = omega
    fomega = - (G / L) * sin(theta)

    return np.array([ftheta, fomega], dtype=float)

a = 0.0
b = 110.0
N = 2000
h = (b - a) / N

tpoints = np.arange(a, b, h)
thetapoints = []
omegapoints = []

theta0 = (pi / 180) * ANGLE
omega0 = (2 * pi) / T
r = np.array([theta0, omega0], dtype=float)
for t in tpoints:
    thetapoints.append(r[0])
    omegapoints.append(r[1])
    k1 = h * f(r, t)
    k2 = h * f(r + 0.5 * k1, t + 0.5 * h)
    k3 = h * f(r + 0.5 * k2, t + 0.5 * h)
    k4 = h * f(r + k3, t + h)
    r += (k1 + 2 * k2 + 2 * k3 + k4) / 6

ypoints = -L * np.cos(thetapoints)
xpoints = L * np.sin(thetapoints)

plt.plot(tpoints, omegapoints, 'k', lw=0.8)
plt.xlabel('t')
plt.ylabel('omega')
plt.savefig('2000_steps.jpeg', dpi=300, format='jpeg')
plt.show()

```

Output of omega vs time:

N = 1000

N = 2000

I strongly believe this is due to the approximation

$$
\sin \theta \approx \theta
$$

for smaller values of $\theta$ in the ordinary differential equation.


问:
Convection equation expanded in Legendre polynomials
In some areas of physics, for example radiation transport in plasmas, it is beneficial to re-express the convection equation in terms of a Legendre polynomial expansion. This results in a set of coupled equations for the expansion coefficients.

So, if the standard convection equation for a function f is:
$\partial f/\partial t=-v_x \partial f/\partial x$, after expansion, the first two of the coupled equations are

$\partial f_0/\partial t=-a v\partial f_1/\partial x$

and

$\partial f_1/\partial t=-b v\partial f_0/\partial x$

where $a$ and $b$ are constants and $v$ is the magnitude of $v_x$.

Since these equations describe convection, their numerical solution suffers from the usual problems of fluid convection, for example spurious short scale oscillations, problems at sharp interfaces, false extrema etc. These problems have been solved for the standard convection equation (e.g. WENO, Van-Leer-type flux limiters, artificial diffusion/viscosity, Riemann solvers etc). Are there equivalents of these for the coupled equations?

One approach I have tried is to re-express these coupled equations in terms of new function g and h like this:

$\partial g/\partial t=-c v\partial g/\partial x$

and

$\partial h/\partial t=-k v\partial h/\partial x$

which are in the usual convection form, but for other reasons I need to work directly with the functions $f_0$ and $f_1$ so this doesn't help.


问:
Using Spherical Bessel Functions in Python
My question relates to using spherical Bessel functions in Python.

If my ODE contains a spherical Bessel function of the form $$j_\ell(tx)$$ and similarly $$y_\ell(tx)$$

for given values of $t$ and $x$.

I used the following:

```<code>from scipy.special import spherical_jn(n, z, derivative)
from scipy.special import spherical_yn(n, z, derivative)
</code>```

I thought I could maybe define jn as follows:

```<code>def Jn(tx,n):
    return (sqrt(pi/(2*(t*x))*jn(n+0.5,t*x))
</code>```

And then just input it into my ODE.

Any help here is appreciated. I can post the ODE if required.


问:
Unable to achieve semi-linear running time in computation of continuant
I am trying to compute the continuant of a list of numbers $a_0, a_1,...,a_n$, defined by the recursion relation: $K_{n+1} = a_{n+1} K_n + K_{n-1}$ and $K_0 = 1$ (see Wikipedia).

I am trying to use divide-and-conquer as suggested by Knuth, 1971, Page 271-272.

However, my implementation:

```<code>mpz_class continuant(size_t s, size_t t, ConvergentList& convergents) {
    if (s == t) return convergents[s];
    if (t - s <= 2) {  // 2 can be adjusted
        mpz_class k_n = convergents[s++], k_n1 = 1, temp;
        while (s != t) {
            temp = k_n1;
            k_n1 = k_n;
            k_n = convergents[s++] * k_n + temp;
        }
        return k_n * convergents[s] + k_n1;
    }
    else {
        auto mid = s + (t - s) / 2;
        return continuant(s, mid, convergents) * continuant(mid + 1, t, convergents)
        + continuant(s, mid - 1, convergents) * continuant(mid + 2, t, convergents);
    }
}
</code>```

runs in $O(n^2)$ time instead of semi-linear time.


问:
Applying boundary Conditions on FEM
I have a partial differential equatons as shown below.
$$\dfrac{d}{dx}((1+x)\dfrac{du(x)}{dx})=0$$
With the following boundary conditions.
$$u(0)=0, u(3)=10$$
To solve it using FEM, I multiplied the PDE with residual function then converted it to its weak form.

Then constructed to global matrices using three elements within the global boundary.

Remember that, to solve the matrix system we need to solve following equation.
$$[A][u]=[b]+[D]$$ $[D]$ corresponds to the matrix who includes the information about boundary conditions.

Since we do not have a constant on the right hand side on the first equation, we do not have the $[b]$ matrix.

I constructed $[D]$ matrix, as shown below.
$$D=\left(\begin{array}{c}(1+x)B \\0 \\0 \\-(1+x)B\end{array}\right)$$

After this point what should I do? How should $1+x$ term affect the imposing Dirichlet BCs to this system?

Title: Applying boundary Conditions on FEM


问:
Scipy.integrate.odeint is returning curves with almost the same frequency for different damping ratios, shouldn't they be different?
I am trying to solve the ODE for a harmonic oscillator using Scipy's odeint solver for different dampening factors.

I'm using the following code, based off of this example:

```<code>from scipy.integrate import odeint

m = 1  # kg
k = 10  # N/m
omega = np.sqrt(k / m)

def eps(c):
  return (c / (2 * mass * np.sqrt(kspring/mass)))

def calc_deri(y, t, eps, omega):
    return (y[1], -eps * omega * y[1] - omega **2 * y[0])

time_vec = np.linspace(0, 8, 80)
yinit = (1, 0)

plt.figure(figsize=(10, 7))
for c in [1, 2, 3, 4, 5, 6]:
  eps_c = eps(c)
  yarr = odeint(calc_deri, yinit, time_vec, args=(eps_c, omega))
  plt.plot(time_vec, yarr[:, 0])

plt.show()
</code>```

My problem is that no matter what combination of k,m, and c I use, I always get a set of curves which have very close oscillation frequencies (I can't tell from the plot whether there is no difference at all or whether the difference is so small as to not show on the graph).

Something like this:

But I am hoping to get something that looks more like this:

What am I doing wrong?


问:
Whi are chi-squared distributions in (C++) boost::random and C++ STL faster than in boost::math?
I am trying to generate random chi-squared numbers in C++, according to some degree of freedom (which can be a float). Several libraries can be used for that purpose, among which the C++11 Standard Library with the header <random>, or two Boost libraries, boost::mathor boost::random.

These libraries don't exactly have the same user interface, but more surprisingly, they definitely don't produce numbers at the same speed. I have performed a test with the next toy code:

```<code>  #include <chrono>
  #include <random>
  #include <boost/math/distributions/chi_squared.hpp>
  #include <boost/random/chi_squared_distribution.hpp>

  unsigned int N = 1e6;
  double df = 2.4;      //degree of freedom for chi-square

  //generator and distributions
  std::mt19937 gen;  //random number engine, boost::mt19937 doesn't change anything
  std::chi_squared_distribution<> chi2STL(df); //STL implementation
  boost::random::chi_squared_distribution<> chi2BoostRandom(df); //boost::random
  boost::math::chi_squared chi2BoostMath(df); //boost::math
  boost::uniform_real<> unif(0., 1.); //uniform distribution for boost::math
                                      //taking STL implementation does not change anything

  //time measure
  auto point1 = std::chrono::high_resolution_clock::now();

  //generate chi2 numbers with STL
  for(unsigned int i = 0; i < N; ++i)
  {
    chi2STL(gen);
  }

  auto point2 = std::chrono::high_resolution_clock::now();

  //generate chi2 numbers with boost::random
  for(unsigned int i = 0; i < N; ++i)
  {
    chi2BoostRandom(gen);
  }

  auto point3 = std::chrono::high_resolution_clock::now();

  //generate chi2 numbers with boost::math
  for(unsigned int i = 0; i < N; ++i)
  {
    boost::math::quantile(chi2BoostMath, unif(gen));
  }

  auto point4 = std::chrono::high_resolution_clock::now();

  //get time
  std::cerr << "Time STL loop: "   << std::chrono::duration_cast<std::chrono::microseconds>(point2 - point1).count()/1000000. << " s" << std::endl;
  std::cerr << "Time Boost Random loop: " << std::chrono::duration_cast<std::chrono::microseconds>(point3 - point2).count()/1000000. << " s" << std::endl;
  std::cerr << "Time Boost Math loop: "   << std::chrono::duration_cast<std::chrono::microseconds>(point4 - point3).count()/1000000. << " s" << std::endl;
</code>```

With the next output:

```<code>Time STL loop: 0.474826 s
Time Boost Random loop: 0.374868 s
Time Boost Math loop: 5.02392 s
</code>```

And it turns out that the boost::math method is way slower than the two other ones. I investigated a little more this surprising point, using a normal law instead of a chi-square law. In that case all three methods produce numbers at the same speed, which makes me think that the boost::quantile method is not the time-consuming one.

By the way, I checked the produced numbers, which all seem to be consistent with the chi-square law, at least for my needs.

It seems that the method implemented in boost::math is just less efficient than the two other ones. (Sadly, it is the only one for which I found implementation details online: it makes use of numerical calculations of the incomplete gamma function).

Is such a simplistic interpretation correct ? And in that case, shouldn't the boost::math implementation be deprecated ?


问:
How can I optimize that loop?
I need to populate a matrix $A_{kl}$, where

$$ k = (m-1)J+n$$
$$ l = (p-1)J+q$$

And

$$m,p = 1, 2, ..., I$$
$$n,q = 1, 2, ..., J$$

Its components are (mnpq). For populate it, I'm using a expensive 4 Do loop

```<code>Do[
 Do[
  Do[
   Do[
    Print[m, n, p, q];
    k = (m - 1) nC + n;
    l = (p - 1) nC + q;
    If[k <= l, A[[k, l]] = cf[Nfunc, xi, yi, wix, wiy, m, n, p, q],
     0];
    , {q, 1, J, 1}]
   , {p, 1, I, 1}]
  , {n, 1, J, 1}]
 , {m, 1, I, 1}]
</code>```

Knowing that $A_{kl}$ for a $I=J=2$, its components are (mnpq)

$$ \begin{bmatrix}
(1111) & (1112) & (1121) & (1122)\\
   & (1212) & (1221) & (1222)\\
 symm. &        & (2121) & (2122)\\
   &        &        & (2222)\\
\end{bmatrix} $$

Does anyone know a more efficient way to populate it?

Title: How can I optimize that loop?


问:
Solving a complex ODE with large number of variables (>1e6 variables) - best practise?
I have to solve a non-linear ODE of the shape
$$\partial_zA=f(A)$$
with $f$ a non-linear function and $A$ a matrix/vector with >1e6 variables (i.e. $A$ is a matrix with >1000x1000 entries). For each calculation of the right hand side of the system (i.e. $f(A)$) I have to calculate roughly five matrix-matrix multiplications $B\cdot A$ with $B$ a constant matrix and roughly five FFTs over each row of $A$. $A$ is using complex numbers.

To solve this system I intended to use a solver such as odeint from boost. I also intended to run the calculation on a GPU, to use the increased speed of matrix-matrix-multiplications and batch FFTs. I had to find out, though, that odeint does not support thrust::complex as value type. Therefore, my current approach is to copy $A$ to the GPU in each step, and return $f(A)$ afterwards to the solver. This is rather costly, and therefore I would like to avoid it.

Thus, are there other strategies (apart from writing my own solver) which could help solving this PDE in reasonable time? Should I still try to go for the GPU, or rather try things like MPI (or a hybrid approach)?
Moreover, I noticed that for several sets of parameters the equation behaves rather stiff. Are there options for solving the system in that case?


问:
Stress visualization for structured and unstrcutured mesh
I put together a FEM solver to do some structural analysis. At the moment I am trying to add stress visualization feature to my code. The post-processing is done with Paraview and the results with which I am trying to verify my code are from the following problem:

The structure is two square plates welded on one edge together to form the figure shown in the images below. The structure is clamped at two points as shown by red circles and a uniformely distributed load is applied on the upper edge of the structure as shown below.
The simulation was run once using structured mesh and second time using unstructured mesh the results using structured mesh are somehow reasonable even though not completely satisfactory but the results from unstructured mesh are absolutely off base!

At a first glimpse what could be done wrong in transition from structured mesh to unstructured mesh (from coding point of view) that results in such discrepancy?


问:
Advice for a topic in a seminar
I am a master student in Mathematics, and I have to prepare a seminar for a course in mathematical methods for applied sciences. I have a good background in numerical analysis for ODEs, PDEs and hence I also would like to run some simulations in Python. I also have some experience with standard numerical methods for SDEs.

The topic of such a talk could be

a chapter from classical Murray's books, so maybe some biological models, but that is of course not mandatory.

I also thought to search for some articles like this one from SIAM, but honestly I'm pretty open to other possibilities.

something that can be modeled by some ODEs that are also "hard" to solve (stiffness) or something that has some properties and one wants the numerical solution to satisfy such a property.

Of course, I do not want the level to be too low, but nor to high, since I'm still a student.

Title: Advice for a topic in a seminar


问:
Finding isocontour lines, populating points on said lines (in MATLAB)
I'm trying to extract data from a potential field. I have an example code of a potential field, using meshgrid in MATLAB. I've pasted the code below.

What I'd like to do next is 'convert' this field into a 2D graph, collect/graph the isocontours, and place data-points along the contour lines, ie. I want to collect a list of XY points along each contour (evenly spaced, probably) making the data easy to transfer or manipulate in other software etc.

ultimately without revealing too much about my research think I need to build up all the data points this way for multiple slices of 2D so that I have a densely populated 3D field of such points, eg. a constellation of spheres whose size depends on the potential of the isocontour it is on.

eg if this was my isoplot, I'd want dots/circles populated along the isocontour lines, with smaller tightly packed dots where there is the steeper gradient near the charges/obstacles, and larger dots in the periphery, filling up the space.

```<code>close all
clear all
clc

xf = 0:0.1:12;
yf = 0:0.1:12;
[X Y] = meshgrid(xf,yf);

Goal=[10;10];
Obs1=[3;3];
Obs2=[9;9];
KG = 10;
Ko = 20;

rG = sqrt((Goal(1)-X).^2+(Goal(2)-Y).^2);
VG = KG*rG; %PF of Goal

ro1 = sqrt((Obs1(1)-X).^2+(Obs1(2)-Y).^2);
Vo1 = Ko./ro1; %PF of Obstacle 1

ro2 = sqrt((Obs2(1)-X).^2+(Obs2(2)-Y).^2);
Vo2 = Ko./ro2; %PF of Obstacle 2

figure
mesh(xf,yf,VG+Vo1+Vo2)
axis([0 12 0 12 0 200])
xlabel('x')
ylabel('y')
zlabel('V')

H4_PotentialNavigation

H4_PlatoonRobots

function [t,x]=H4_PotentialNavigation()

    clc

    time = [0 300];

    x0 = [0;0;0];

    Goal = [10;10];
    Obs1 = [3;3];
    Obs2 = [9;9];
    KG = 30;
    Ko = 30;

    [t x] = ode23(@vehicle,time,x0)
    figure
    plot(x(:,1),x(:,2),'r',Goal(1),Goal(2),'o',Obs1(1),Obs1(2),'x',Obs2(1),Obs2(2),'x')
    axis([0 12 0 12])

    function dx=vehicle(t,x)

        rG = sqrt((Goal(1)-x(1))^2+(Goal(2)-x(2))^2);
        FGx = KG*(Goal(1)-x(1))/rG;
        FGy = KG*(Goal(2)-x(2))/rG;

        ro1 = sqrt((Obs1(1)-x(1))^2+(Obs1(2)-x(2))^2);
        Fo1x = -Ko*(Obs1(1)-x(1))/ro1^3;
        Fo1y = -Ko*(Obs1(2)-x(2))/ro1^3;

        ro2 = sqrt((Obs2(1)-x(1))^2+(Obs2(2)-x(2))^2);
        Fo2x = -Ko*(Obs2(1)-x(1))/ro2^3;
        Fo2y = -Ko*(Obs2(2)-x(2))/ro2^3;

        Fx = (FGx+Fo1x+Fo2x);
        Fy = (FGy+Fo1y+Fo2y);

        alpha = atan(Fy/Fx);

        v = 1;
        L = 2;

        if rG < 0.05 % Vehicle @ goal - stop
            v = 0;
        end

        K = 2;
        ph = K*(alpha-x(3));
        dx = [v*cos(ph)*cos(x(3));
              v*cos(ph)*sin(x(3));
              v*sin(ph)/L];
    end
end

function [t,x] = H4_PlatoonRobots()

clc

time = [0 300];
x0 = [0;0;0;0;0;0;0;0;0;0;0;0];

Goal = [10;10];
Obs1 = [3;3];
Obs2 = [9;9];
KG = 30;
Ko = 30;

[t x] = ode23(@vehicles,time,x0)
figure
plot (x(:,1),x(:,3),'b',x(:,5),x(:,7),'r',x(:,9),x(:,11),'g',...
    Goal(1),Goal(2),'o',Obs1(1),Obs1(2),'x',Obs2(1),Obs2(2),'x')
axis([0 12 0 12])

    function dx = vehicles(t,x)
        Goal = [10;10];
        Obs1=[3;3];
        Obs2=[9;9];
        KG=20;
        Ko=30;
        Kij=20;
        KLi=100;
        KF=10;
        rD=0.5; % Separation from Leader

        rG = sqrt((Goal(1)-x(1))^2+(Goal(2)-x(3))^2); % Leader and Goal
        FGx = KG*(Goal(1)-x(1))/max([rG 0.01]); % Division by Zero catch
        FGy = KG*(Goal(2)-x(3))/max([rG 0.01]);

        ro1 = sqrt((Obs1(1)-x(1))^2+(Obs1(2)-x(3))^2);
        Fo1x = -Ko*(Obs1(1)-x(1))/ro1^3;
        Fo1y = -Ko*(Obs1(2)-x(2))/ro1^3;

        ro2 = sqrt((Obs2(1)-x(1))^2+(Obs2(2)-x(3))^2);
        Fo2x = -Ko*(Obs2(1)-x(1))/ro2^3;
        Fo2y = -Ko*(Obs2(2)-x(2))/ro2^3;

        FrLx = -KF*x(2);
        FrLy = -KF*x(4);

        Fx = (FGx + Fo1x + Fo2x + FrLx);
        Fy = (FGy + Fo1y + Fo2y + FrLy);

        xL = [0 1; 0 0]*[x(1);x(2)]+[0;1]*Fx; % Leader Dynamics
        yL = [0 1; 0 0]*[x(3);x(4)]+[0;1]*Fy;

        r1o1 = sqrt((Obs1(1)-x(5))^2+(Obs1(2)-x(7))^2); % Follower 1 and Obs1
        F1o1x = -Ko*(Obs1(1)-x(5))/r1o1^3;
        F1o1y = -Ko*(Obs1(2)-x(7))/r1o1^3;

        r1o2 = sqrt((Obs2(1)-x(5))^2+(Obs2(2)-x(7))^2); % Follower 1 and Obs2
        F1o2x = -Ko*(Obs2(1)-x(5))/r1o2^3;
        F1o2y = -Ko*(Obs2(2)-x(7))/r1o2^3;

        rL1 = sqrt((x(1)-x(5))^2+(x(3)-x(7))^2); % Follower 1 and Leader
        FL1x = KLi *(rL1-rD)*(x(1)-x(5))/max([rL1 0.01]);
        FL1y = KLi *(rL1-rD)*(x(3)-x(7))/max([rL1 0.01]);

        r21 = sqrt((x(9)-x(5))^2+(x(11)-x(7))^2); % Follower 1 and Follower 2
        F21x = -2* Kij *(x(9)-x(5))/max([r21^4 0.01]);
        F21y = -2* Kij *(x(11)-x(7))/max([r21^4 0.01]);

        Fr1x = -KF*x(6);
        Fr1y = -KF*x(8);

        F1x = (F1o1x+F1o2x+FL1x+F21x+Fr1x);
        F1y = (F1o1y+F1o2y+FL1y+F21y+Fr1y);

        x1 = [0 1; 0 0]*[x(5);x(6)]+[0;1]*F1x; % Follower 1 Dynamics
        y1 = [0 1; 0 0]*[x(7);x(8)]+[0;1]*F1y;

        r2o1 = sqrt((Obs1(1)-x(9))^2+(Obs1(2)-x(11))^2); % Follower 2 and Obs1
        F2o1x = -Ko*(Obs1(1)-x(9))/r2o1^3;
        F2o1y = -Ko*(Obs1(2)-x(11))/r2o1^3;

        r2o2 = sqrt((Obs2(1)-x(9))^2+(Obs2(2)-x(11))^2); % Follower 2 and Obs2
        F2o2x = -Ko*(Obs2(1)-x(9))/r2o2^3;
        F2o2y = -Ko*(Obs2(2)-x(11))/r2o2^3;

        rL2 = sqrt((x(1)-x(9))^2+(x(3)-x(11))^2); % Follower 2 and Leader
        FL2x = KLi *(rL2-rD)*(x(1)-x(9))/max([rL2 0.01]);
        FL2y = KLi *(rL2-rD)*(x(3)-x(11))/max([rL2 0.01]);

        r12 = sqrt((x(5)-x(9))^2+(x(7)-x(11))^2); % Follower 2 and Follower 1
        F12x = -2* Kij *(x(5)-x(9))/max([r12^4 0.01]);
        F12y = -2* Kij *(x(7)-x(11))/max([r12^4 0.01]);

        Fr2x = -KF*x(10);
        Fr2y = -KF*x(12);

        F2x = (F2o1x+F2o2x+FL2x+F12x+Fr2x);
        F2y = (F2o1y+F2o2y+FL2y+F12y+Fr2y);

        x2 = [0 1; 0 0]*[x(9);x(10)]+[0;1]*F2x; % Follower 2 Dynamics
        y2 = [0 1; 0 0]*[x(11);x(12)]+[0;1]*F2y;

        dx = [xL(1);xL(2);yL(1);yL(2);x1(1);x1(2);y1(1);y1(2);x2(1);x2(2);...
            y2(1);y2(2)];

        if x(2) == 0 && x(4) == 0 % avoid local minima
            w = 0.001*randn(2,1);
            dx = dx + [0;w(1);0;w(2);0;0;0;0;0;0;0;0];
        end

        if r12 == 0 % Seperate followers from initial condition
            w = 0.001*randn(4,1);
            dx = dx + [0;0;0;0;0;w(1);0;w(2);0;w(3);0;w(4)];
        end
    end
end
</code>```


问:
Add orthogonality constraint to PDE over entire domain finite element
I have a nonlinear differential equation for $u:\mathbb{R}^2\rightarrow\mathbb{R}^3$ that I can express in the form:
$ \nabla (g(u)_{ij} \nabla(u)^j) = 0 $
 which is similar to the "Nonlinear Poisson" example given in the Fenics tutorial except with the nonlinearity $g$ now a matrix.

However I would like to add an additional constraint of the form:

$u_x \cdot u_y = 0$.

I suspect there is a way to do this using Lagrange Multipliers, but I am not sure how to express this.


问:
Stokes problem with imposed acceleration on boundaries (projection scheme)
I am trying to solve FSI problems with finite elements and using a projection scheme (I am taking as reference the review of Guermond:

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

These methods are based on the Helmoltz decomposition as I understand:
$$[L^2(\Omega)]^d = H \oplus \nabla H^1 (\Omega),$$
where
$$H=\{ v\in [L^2(\Omega)]^d; \nabla\cdot v = 0; v\cdot n|_{\Gamma} = 0\}.$$
So to solve the Navier-Stokes equations
$$\frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot \nabla)\mathbf{u} = -\nu \Delta\mathbf{u} + \nabla p,\\
\nabla\cdot \mathbf{u} = 0,$$
first one solves a problem for velocity only, where the pressure is approximated in some way, and then solves a problem for the pressure to enforce the incompressibility constraint
$$ \frac{\partial\mathbf{u}}{\partial t} -\nabla p = 0, \qquad \frac{\partial p}{\partial n}\bigg|_{\Gamma} = 0,$$
where the boundary condition should come from the Helmoltz decomposition.
My doubt is that I usually see this framework for standard Navier-Stokes problems, where boundaries are fixed and velocity is imposed with a Dirichlet condition. If that is the case, taking the product of the NS equations with the normal to the boundary $n|_{\Gamma}$, one sees that the condition respected by pressure at the boundary is
$$ \frac{\partial p}{\partial n}\bigg|_{\Gamma} = \nu \Delta \mathbf{u} \cdot n|_{\Gamma}.$$

Now, if I want to apply this to a FSI problem, one of the boundaries is moving and therefore the velocity is assigned but changing in time, that is an acceleration is imposed. My question is, does the classic framework of pressure projection still work? And how does one derive rigorously the two subproblems to be solved (velocity and pressure), and more importantly, which are the boundary conditions to be imposed on pressure and velocity at the boundaries and at the FSI interfaces?
To me, it doesn't look like that the Helmoltz decomposition in its basic form works if the velocity on Dirichlet boundaries is time dependent.
Intuitively, if again I take the product of NS equations with the boundary normal $n|_{\Gamma}$, the $\frac{\partial \mathbf{u}}{\partial t}$ does not cancel if the boundary is accelerating, therefore the pressure satisfies now
$$ \frac{\partial p}{\partial n}\bigg|_{\Gamma} = \left(-\frac{\partial \mathbf{u}}{\partial t} + \nu\Delta\mathbf{u} \right)\cdot n|_{\Gamma}.$$
But do I simply add the acceleration as boundary condition to the pressure problem? And how do I rigorously derive it? Also, is there a good reference on this that I am missing out?


问:
Difficulty to prove the coercivity of $a(u, v)=\int_{\Omega} \nabla u \cdot \nabla v$
I'm stuck at proving the coercivity of the bilinear functional in the variational formulation of the problem:

\begin{array}{c}
-\nabla^{2} u=f \quad \text { in } \Omega \\
u=g_{D} \text { on } \partial \Omega_{D} \quad \text { and } \quad \frac{\partial u}{\partial n}=g_{N} \text { on } \partial \Omega_{N}
\end{array}

where $\partial \Omega_{D} \cup \partial \Omega_{N}=\partial \Omega$ and $\partial \Omega_{D}$ and $\partial \Omega_{N}$ are distinct.

The variational formulation is :

Find $u \in \mathcal{H}_{E}^{1}$ such that
$$
\int_{\Omega} \nabla u \cdot \nabla v=\int_{\Omega} v f+\int_{\partial \Omega_{N}} v g_{N} \quad \text { for all } v \in \mathcal{H}_{E_{0}}^{1}
$$
where:
\begin{aligned}
\mathcal{H}_{E}^{1} &:=\left\{u \in \mathcal{H}^{1}(\Omega) | u=g_{D} \text { on } \partial \Omega_{D}\right\} \\
\mathcal{H}_{E_{0}}^{1} &:=\left\{v \in \mathcal{H}^{1}(\Omega) | v=0 \text { on } \partial \Omega_{D}\right\}
\end{aligned}

The proof by contradiction holds for $a(u, v)=\int_{\Omega} \nabla u \cdot \nabla v+\int_{\partial \Omega} uv$, but here, our bilinear form doesn't contain the term $\int_{\partial \Omega}  uv$

Thank you for your help


问:
Discretisation by FVM with triangulare mesh
I work on electromagnetic modeling, and I use the finite volume method as a tool, so I want to develop code by the modified finite volume method (with triangular mesh), my problem studied is in a permeable environment (ferromagnetic materials).

in first i try to discetize  $-\mathrm{div}(\mathrm{grad}(A)=\mu_0(J + \mathrm{rot}(M))$ by fvm with rectangular mesh, (see picture)

My question Now

I search to determine the magnetic vector potential $A$ at the centers of gravity of the triangles,

my problem is in how to discretize this equation above by FVM triagnular like in rectangular mesh

$$-\left(\frac{\partial^2 A}{\partial x^2} + \frac{\partial^2 A}{\partial y^2}\right)  = \mu_0 \left(J_S + \frac{\partial M_y}{\partial x} - \frac{\partial M_x}{\partial y}\right)\, .$$


问:
Finding points inside cells of power (generalized Voronoi) diagram
Suppose we have a set of points $p_1,\ldots,p_n\in\mathbb R^d$ as well as a set of weights $w_1,\ldots,w_n\in\mathbb R$.  Recall that the power cell associated to the pair $(p_k,w_k)$ is given by:
$$\mathcal C_k:=\{x\in\mathbb R^d:\|x-p_k\|_2^2-w_k\leq\|x-p_\ell\|_2^2-w_\ell\ \forall \ell\in\{1,\ldots,n\}\}\subseteq\mathbb R^d$$
As a special case, the $\mathcal C_k$'s become Voronoi cells when all the weights equal zero.

Assuming $\mathcal C_k\neq\emptyset$ for all $k$, is there an easily-implemented, moderately efficient (at least for low dimensions $d$) numerical algorithm for finding a set of points $x_1,\ldots,x_n\in\mathbb R^d$ with $x_k\in\mathcal C_k$ for all $k\in\{1,\ldots,n\}$?

I don't care which point we find, so long as it is inside the power cell.  The obvious algorithm "construct the power diagram" is hard to implement---especially for $d>3$---and potentially wasteful.  My hope is that for low $d$ maybe there's a simpler way to solve this problem that scales roughly in $n$.

If it helps, we can possibly assume we have a set of points $\tilde x_1,\ldots,\tilde x_n\in\mathbb R^d$ such that each $x_k$ isn't too far away from $\mathcal C_k$, i.e., try to make a local update.


问:
Why are "instructions retired" more stable than "cycles" or "task-clock" when using "perf"?
I tested the perf tool in Ubuntu 18.04 on a simple benchmark in our compiler (parsing some file). I ran perf several times and put here the slowest:

```<code>$ sudo perf stat ./parse
Construct
Parse
Parsing: 93ms
Counting: 14ms
Total: 107ms
Count: 450009
String size (bytes):      2250042
Allocator usage (bytes): 48400928

 Performance counter stats for './parse':

        112.329506      task-clock (msec)         #    0.999 CPUs utilized
                 0      context-switches          #    0.000 K/sec
                 0      cpu-migrations            #    0.000 K/sec
            14,560      page-faults               #    0.130 M/sec
       418,991,316      cycles                    #    3.730 GHz
       843,385,199      instructions              #    2.01  insn per cycle
       139,033,655      branches                  # 1237.730 M/sec
         1,033,015      branch-misses             #    0.74% of all branches

       0.112438670 seconds time elapsed
</code>```

and fastest run:

```<code>$ sudo perf stat ./parse
Construct
Parse
Parsing: 80ms
Counting: 12ms
Total: 93ms
Count: 450009
String size (bytes):      2250042
Allocator usage (bytes): 48400928

 Performance counter stats for './parse':

         97.922823      task-clock (msec)         #    0.998 CPUs utilized
                 0      context-switches          #    0.000 K/sec
                 0      cpu-migrations            #    0.000 K/sec
            14,559      page-faults               #    0.149 M/sec
       376,416,206      cycles                    #    3.844 GHz
       843,259,218      instructions              #    2.24  insn per cycle
       138,973,382      branches                  # 1419.213 M/sec
         1,031,451      branch-misses             #    0.74% of all branches

       0.098073332 seconds time elapsed
</code>```

When comparing the various statistics, here are the percentage improvements from the slowest to the fastest run:

```<code>task-clock      12.83% faster
cycles          10.16% less
frequency        3.06% faster
instructions     0.015% less
insn per cycle  11.44% more
branches         0.043% less
branch-misses    0.15% less
</code>```

Why is instructions 3 orders of magnitude more stable than the total time (task-clock) or cycles?

One explanation could be that my CPU (Intel(R) Xeon(R) CPU E3-1505M v6 @ 3.00GHz) executes roughly the same number of instructions each time, but depending on various factors (temperature, past history, other programs running in linux at the same time, etc.) it can sometimes execute ~10% more instructions per cycle by hiding latency / executing instructions in parallel, which means it needs ~10% less cycles to execute all of them, and given the fact that the frequency is roughly the same (within ~3%) this translates to ~10% less time. Is this a correct explanation?


问:
How to numerically evaluate this double Integral?
I want to evaluate the following integral:

$$\int_{0}^{60} \ \left(\int_{0}^{2z} 0.5\cdot t \left(\mathrm{erf}(t-a) -1 \right)J_{0}(qt)\mathrm{d}t \right)^2 \mathrm{exp}\left(-\frac{(z-a)^2}{2s^2}\right)\mathrm{d}z $$

Where a, q and s are constants.
Without the square outside the first integral, this can be easily evaluated using dblquad from scipy as:

```<code>import math
from scipy.integrate import dblquad
from scipy.special import erf, jv

def h(t, z):
    return f(t) * g(z)

def f(t):
    return 0.5 * t * (erf(t - a) - 1) * jv(0, q * t)

def g(z):
    return math.exp(-((z - a) ** 2)/(2 * (s ** 2)))

if __name__ == '__main__':
    a, q, s = 0, 1, 2  # set the constants
    result, abserr = dblquad(h, 0, 60, lambda z: 0, lambda z: 2 * z)
    print(f'result: {result}, abserr: {abserr}')

</code>```

I looked about the square of integrals, and as shown here square of integral is equivalent to double integral:

$$\left(\int_a^bf(x)\text{d}x\right)^2 = \int_a^b \int_a^b f(x) f(y) \text{d}x\text{d}y$$

This makes the problem more complicated as now I have to evaluate triple integral which makes evaluation slower. Is there any better way to tackle this problem? Any help would be highly appreciated.


问:
Conjugate Gradient for nonlinear equation system
Is it possible to apply adaptions of the conjugate gradient algorithm i.e. Fletcher-Reeves, Polak-Ribere or others to systems of nonlinear equations? How should the equation system be adjusted so one can use those algorithms?

I am asking this because in all of the literature, the pseudocode says to compute the conjugated directions by gradients, i.e.:

$p_{k+1} = -\nabla f_{k+1} + \beta_{k+1} p_k $
with
$\beta_{k+1} = \frac{\nabla f_{k+1}^T \nabla f_{k+1}}{\nabla f_{k}^T \nabla f_{k}}$

But in an system of equations, $\nabla f$ results in the Jacobian. What do I miss?


问:
When to discretize nonlinear Poisson Equation
I am trying to solve a nonlinear poisson equation of the form:

$u_{xx} + f(u_y)u_{yy} = 0$.

In trying to get a handle on this problem, it seems like there are two approaches. I could either discretize this directly, and end up with a nonlinear system of equations for $u_{i,j}$ or I could apply newtons method first and then discretize. Which is preferable and why?


问:
Division and Modulus for Large DIVISORS
Excel and also Excel VBA have no built-in support for arbitrary precision arithmetic. There are a few very large add-ins that can be installed to do these sorts of calculations where the operands are hundreds to thousands of digits long.

For recreation I am attempting to code in VBA similar functions without creating a big-integer class.

I am not a mathematician. But I have discovered that it was fairly easy to implement functions for Addition, Subtraction, and Multiplication.

For example, here is my function for addition for arbitrarily large operands:

```<code>Function AddStrings<span class="math-container">$(ByVal a$</span>, ByVal b$)

    Dim carry&, i&, k&, ka&, kb&, p&, pLen&, s$, z&, subt
    Const size& = 28

    ka = Len(a): kb = Len(b)
    Select Case kb
        Case Is < ka: k = ka: If k > kb Then b = String<span class="math-container">$(k - kb, "0") & b
        Case Else:    k = kb: If k > ka Then a = String$</span>(k - ka, "0") & a
    End Select
    s = Space$(k + 2)

    p = k + 1
    For i = -Int(-k / size) To 1 Step -1
        p = p - size
        pLen = size: If p < 1 Then p = 1: pLen = k Mod size
        subt = CDec(Mid<span class="math-container">$(a, p, pLen)) + Mid$</span>(b, p, pLen)
        If carry Then subt = subt + carry
        z = Len(subt)
        Select Case z
            Case pLen:      carry = 0: Mid<span class="math-container">$(s, p, pLen) = subt
            Case Is > pLen: carry = 1: Mid$</span>(s, p, pLen) = Mid<span class="math-container">$(subt, 2)
            Case Else:      carry = 0: Mid$</span>(s, p, pLen) = String<span class="math-container">$(pLen - z, "0") & subt
        End Select
    Next
    If carry Then s = "1" & s
    AddStrings = Trim$</span>(s)

End Function
</code>```

This function seems to work perfectly regardless of the size of operands I give it. It breaks the operands down into 28-digit chunks and iteratively processes them.

My multiplication routine works perfectly as well. It uses 14-digit chunks.

However, Division and Mod have me scratching my head.

My strategy that worked well for the Addition, Subtraction, and Multiplication was segmenting the operand strings into a number of digits that would fit into VBA native data types and doing the calculations piecemeal. This worked very well and I can now do those calculations for string values hundreds of digits long.

But this does NOT work for division or modulus.

By tinkering I have discovered that for modulus I can have the dividend any length and get the correct result through iteration of the native Mod function, as long as the divisor fits into a native data type.

For example this function works for HUGE dividends but only for divisors up to eight digits in length:

```<code>Function SuperMod(ByVal a<span class="math-container">$, b$</span>)
    Do While Len(a) > 9
        a = Mid<span class="math-container">$(a, 1, 9) Mod b & Mid$</span>(a, 10)
    Loop
    SuperMod = a Mod b
End Function
</code>```

For example:

```<code>MsgBox SuperMod("654654321902010548705495321500805040401121210900506078798456849816241968462130000051204789540453745898347543753489751261610093245892737467263468234623689898054456361110010151021101827", 99999999)
</code>```

...correctly displays:  84245606

But attempting even one more digit in the divisor results in an error because the native Mod function cannot handle divisors that large.

My goal is to figure out a way where the divisor can be any arbitrary length, just like the dividend currently can be.

I was able to figure out a different way where the divisor can be up to 15-digits with an arbitrarily sized dividend by replacing a Mod b with a - Int(a / b) * b. But this fails for divisors larger than 15 digits because the native division and Int() functions fail at that point:

```<code>Function SuperMod2(ByVal a<span class="math-container">$, b$</span>)
    Do While Len(a) > 15
        a = Mid<span class="math-container">$(a, 1, 15) - Int(Mid$</span>(a, 1, 15) / b) * b & Mid$(a, 16)
    Loop
    SuperMod2 = a - Int(a / b) * b
End Function
</code>```

This version allows for:

```<code>MsgBox SuperMod2("654654321902010548705495321500805040401121210900506078798456849816241968462130000051204789540453745898347543753489751261610093245892737467263468234623689898054456361110010151021101827", 999999999999999)
</code>```

...correctly displays:  194745570864677

VBA does support a Decimal (28-digits of precision) data subtype for Variant variables, but for some reason, that does not help in this case. But even if it did, it would just extend by roughly twice the size of divisors this method can handle.

What sort of math manipulation can I employ to allow for divisors of arbitrary size both for modulus and for division?

I do not see how to chunk-ify both operand strings like I was able to do for Addition, Subtraction, and Multiplication to iteratively process the results.

Everything I find by searching (even here on Stack) refers to large dividends. As you can see, I have that figured out. My problem is arbitrarily large DIVISORS.

I would very much appreciate anyone taking the time to point me in the right direction.

IS CHUNKING THE DIVISOR OPERAND TO ITERATE DIVISION POSSIBLE? As shown and outlined, that technique works fine for Addition, Subtraction, and Multiplication. I cannot see how to make it work for the DIVISOR operand in Division or Mod.

IF NOT, WHAT OTHER TECHNIQUE CAN I PERUSE FOR ARBITRARILY SIZED DIVISORS?

Title: Division and Modulus for Large DIVISORS


问:
Plot of ratio of two integrals:
Consider the following integrals
$$
I_1(x) =\int_0^\infty\mathrm dy \frac{F(x + \mathrm iy) − F(x −\mathrm iy)}{\mathrm e^{2πy}-1},
$$

And

$$I_2(x) =\int_1^x F(t)dt$$

Where, $ F(z) = \sin^2[π\Gamma(z)/(2z)]$

I want to plot (with considerable accuracy)
$$R(x)=I_1(x)/I_2(x)$$

Also , see if it's apparent from plot : $R(x)\rightarrow0$ as $x\rightarrow\infty$

I want growth rate of local minimas and Maximas of $I_1(x)$ as $x\rightarrow\infty$

See this MSE post for more details .

Title: Plot of ratio of two integrals:


问:
Circumferencial waves on a cylinder/sphere
I was wondering how we can introduce $e^{ik.x}$ terms associated with circumferencially propagating waves? In this case $\hat{e}_\theta$ is the direction of wave propagation. However, I was not able to express the exponent  of exponential using a dot product because of polar coordinates. How can we introduce this dot product? Can spherical harmonics or Fourier transform help here?


问:
Unsupervised event prediction with time series data
I have to solve a problem using unsupervised time series forecast. I have a dataset where three columns are given,
1. Endtime: when the row was generated ( it generates a row after each 54-55 seconds)
2. Min: min value of a wave
3. Max: max value of a wave.

Here I need to predict when we would observe an abrupt change (increase) in amplitude of the wave using unsupervised method. The increase may occur at the start or end or in middle. It's not fixed. ( Below is an example snippet of data and plot for your reference).

Please let me know if anyone has done similar projects. It would be a great help if you could give me some hints.


问:
Does order of data points matter for approximating AUC with unit time steps?
I have a time series of data where the increment is every minute. In order to approximate AUC, I just compute the sum of the data values, since they would all be multiplied by 1 anyway per Riemann sum approximation of the integral. </p>

<p>I'm fairly confident that this means the data can be scrambled and still produce the correct AUC, but I was hoping to check here in case I'm missing anything silly. </p>


问:
Solve 1D saturation equation using FEM
I want to solve the following 1D equation to model the time evolution of saturation in a porous medium using finite elements (I'm a beginner).
The permeability $K$ and the pressure $p$ depend on the saturation $s$.
The initial condition in the field $\Omega = (0,l)$ is given as: $$s(t=0)=0.1, \text{ in } \Omega $$
$$ s[0] = 0.5; \text{left boundary}$$

$$\frac{\partial s}{\partial t} = - \frac{1}{\phi} \nabla \cdot u \\
u = \frac{K(s)}{\mu} \nabla p(s)$$

For the time-development, I discretized the time derivative using finite differences (explicit Euler).
I'm aware that this is not optimal, because a linear system needs to be solved to obtain the new timestep, but I think that it should work for a preliminary solution.

I used linear shape functions and the standard Galerkin approach (shape function = weight function).
My approximate solution is:
$$ s = \sum_i \tilde s_i N_i$$
Where $N$ is the linear shape function. For the discrete weak formulation, with the exponent $n$ denoting the time, I get:

$$ MS^{n+1} = MS^{n} + \frac{\Delta t}{\phi}\frac{K(S^n)}{\mu}Dp^n $$
With the mass matrix:
$$ m_{ij} = \int_\Omega N_iN_jd\Omega$$
and the stiffness matrix: $$d_{ij} = \int_\Omega \nabla N_i \nabla N_j d\Omega$$.

However, the approach seems to be faulty because the solution doesn't converge and it oscillates strongly at the saturation jump.
I had a look at different problems like the heat equation but found it to be different from my system. Especially because the saturation dependencies of $K$ and $p$ are strongly non-linear (exponential functions).

Could it be that the Standard Galerkin approach is not applicable?
I would be very grateful for any hints.

PS: I could successfully simulate the case using a particle-based simulation (SPH).

Title: Solve 1D saturation equation using FEM


问:
Volumetric Mesh Cutting Algorithms For Simulating Sculpture
I am interested in whether there an available algorithms to give me insights on how to cut a volumetric mesh?</p>


问:
COMSOL Circularl polarization
I'm having some problems trying to implement circularly polarized light in COMSOL Muliphysics.

For a isotropic homogenous media, I've obtained without problems the TE and TM reflectance curves. Following the Fresnel Equation application tutorial, I've defined a periodic port for incident light with Ex=1 Ey=i Ez=0 and output periodic port with Ex=1 Ey=i Ez=0. Then, I just insert Floquet periodic boundary conditions over the x and y planes of my rectangular structure. The results are somehow so goofy.

What do you think? is there any missing condition I've forgotten?

Title: COMSOL Circularl polarization


问:
Triangle on top of diagonal least squares
I need to solve many least squares problems with the following matrices:
$$
\pmatrix{ R \\ D_i }
$$
where $R$ is upper triangular and $D_i$ is diagonal. $R$ is the same for all the problems, while $D_i$ changes for each problem. Is there some way to factor $R$ to make it efficient to solve these problems?

A Givens approach doesn't seem to work, since it will introduce many non-zero entries into the $D_i$ part of the matrix.

Right now I am doing a full QR factorization of each matrix, which doesn't account for the sparse structure of the problem and is inefficient.