问:
Curve fitting using a piecewise polynomial
I am trying to fit a piecewise polynomial function

Code:

```<code>import numpy as np
import scipy
from scipy.interpolate import UnivariateSpline, splrep
from scipy.optimize import curve_fit
from matplotlib import pyplot as plt

def piecewise_func(x, X, Y):
    """
    cond_l: condition list
    func_l: function list
    """
    spl = UnivariateSpline(X, Y, k=3, s=0.5)
    tck = (spl._data[8], spl._data[9], 3)  # tck = (knots, coefficients, degree)
    p = scipy.interpolate.PPoly.from_spline(tck)

    cond_l = []
    func_l = []
    for idx, i in enumerate(range(3, len(spl.get_knots()) + 3 - 1)):
        cond_l.append([(x >= p.x[i] & x < p.x[i + 1])])
        func_l.append([lambda x: p.c[3, i] + p.c[2, i] * x + p.c[1, i] * x ** 2 + p.c[0, i] * x ** 3])

    return np.piecewise(x, cond_l, func_l)

if __name__ == '__main__':

    xdata = [0.28190937,  0.63429607,  0.91620544,  1.68793236, 2.32350115,  2.95215219, 4.5,
         4.78103382,  7.2, 7.53430054,  8.03627018,  9., 9.86212529, 11.25951191, 11.62658532, 11.65598578, 13.90295926]
    ydata = [0.36273168,  0.81614628,  1.17887796,  1.4475374,   5.52692706,  2.17548169, 3.55313396,  3.80326533,  7.75556311,  8.30176616, 10.72117182, 11.2499386,
        11.72296513, 11.02146624, 14.51260631, 20.59365525, 21.77847853]

    spl = UnivariateSpline(xdata, ydata, k=3, s=1)

    plt.plot(xdata, ydata, '*')
    plt.plot(xdata, spl(xdata))
    plt.show()

    p, e = curve_fit(piecewise_func, xdata, ydata)
    # x_plot = np.linspace(0., 0.15, len(x))
    # plt.plot(x, y, "+")
    # plt.plot(x, (piecewise_func(x_plot, *p)), 'C3-', lw=3)
</code>```

I tried the UnivariateSpline function to interpolate, I see the following result

However, I don't want the polynomial curve to pass through all data points. I tried varying the smoothing factor but I am not able to obtain something like the one below.

Expected output:

I'm trying curve fitting (https://stackoverflow.com/questions/52483191/use-univariatespline-to-fit-data-tightly) to get the expected output and I have the following issues.

piecewise_func in the code posted returns the piecewise polynomial.
Passing this to curve_fit(piecewise_func, xdata, ydata) returns an error

Error:
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
ValueError: diff requires input that is at least one dimensional

I am not sure what is wrong.

Suggestions on how to get the expected fit will be
of great help.

Cross-posted: here


问:
Monte Carlo simulation for the quantum oscillator in the path integral approach

The theory

Consider a quantum harmonic oscillator described by the potential $V(q)=\frac{1}{2}m\omega^2 x^2$. In the path integral formulation, the partition function can be written as $$Z\propto\int Dx(\tau)\exp\left(-\frac{S_E}{\hbar}\right) $$ where $\tau$ is Euclidean time and $S_E$ is the Euclidean action for closed paths $$S_E=\int_0^{\beta \hbar}d\tau\left(\frac{1}{2}m\dot x^2(\tau)+\frac{1}{2}m\omega^2 x^2(\tau)\right).$$

To attack this problem numerically, we can make this quantity discrete by choosing a sequence of  $N$ points $x_1,...,x_N$ on the path $x(\tau)$ with distance $a$ between each other; rescaling these variables as $x_j=\sqrt{\frac{\hbar}{m\omega}}y_j$ we can write
$$\frac{S_E^{discrete}}{\hbar}=\sum_{j=1}^N\left(\frac{1}{2\eta}(y_{j+1}-y_j)^2+\frac{\eta}{2}y_j^2\right) $$ where $\eta=a\omega$. The continuum is then recovered in the limit $\eta\to 0$.
In this discrete form, the internal energy of the system transforms as follows: $$\frac{U^{cont}}{\hbar\omega}=\partial_\beta\log Z=\hbar\omega\left(\frac{1}{2}+\frac{1}{e^{N\eta}-1}\right) \  \Rightarrow \  \frac{U^{discrete}}{\hbar\omega}=\frac{1}{2\eta}-\frac{1}{2\eta^2}\langle \Delta y^2\rangle+\frac{1}{2}\langle y^2\rangle$$ where $y^2, \Delta y^2=(y_{j+1}-y_j)^2$ are first calculated as averages over a single path, and the averages $\langle y^2\rangle, \langle\Delta y^2\rangle$ are taken over many paths distributed according to the function $$p(y_1,...,y_N)=e^{-S_E^{discrete}/\hbar}.$$

The algorithm

In order to sample paths that follow this distribution, we can employ a MCMC method based on the Metropolis algorithm:

Initialize the discrete path in a cold configuration (all zeros) or hot configuration (random numbers)

For each Monte Carlo sweep, loop on all $j=1,...,N$ and propose to modify $y_j$ with a uniformly sampled $y_p\leftarrow [y_j-\delta,y_j+\delta]$ for some parameter $\delta$.

The update is accepted with probability $r=\min\left(1,e^{-\Delta S_{E}^{discrete}}\right)$, where

\begin{align*}
\Delta S_{E}^{discrete}&=(y^p)^2\left(\frac{\eta}{2}+\frac{1}{\eta}\right)-\frac{y^p}{\eta}(y_{j_0+1}+y_{j_0-1})-y_{j_0}^2\left(\frac{\eta}{2}+\frac{1}{\eta}\right)\\&-\frac{1}{2}y_{j_0}(y_{j_0+1}+y_{j_0-1}).
\end{align*}

The code

Here is my implementation of the algorithm in Python:

```<code>import numpy as np
import scipy as sp
import random as rand
import time as tm
import matplotlib.pyplot as plt

#STRAIGHT-LINE INITIALIZATION
def cold_path(N):

    return np.zeros(N)

#RANDOM PATH INITIALIZATION
def hot_path(N):

    return np.random.uniform(-1,1,N)

#MONTE CARLO SIMULATION
def Metropolis_HO(start,N,eta,delta,ntimes):

    #Set path at initial step
    if start=='cold':
        path=cold_path(N)

    elif start=='hot':
        path=hot_path(N)

    else:
        raise Exception('Choose either hot or cold starting path configuration.')

    #Initialize arrays of observables
    y_squared=np.zeros(ntimes)
    delta_y_squared=np.zeros(ntimes)

    #Useful constants
    c1=1./eta
    c2=1./eta+eta/2.

    #Iterate loop on all sites
    for i in range(ntimes):
        for j in range(N):

            #Set y as j-th point on path
            y=path[j]

            #Proposed modification
            y_p=rand.uniform(y-delta,y+delta)

            #Calculate accept probability
            force=path[(j+1)%N]+path[(j-1)%N]
            p_ratio=c1*y_p*force-c2*(y_p*y_p)-c1*y*force+c2*(y*y)

            #Accept-reject
            if rand.random()<np.exp(p_ratio):
                path[j]=y_p

        #Average of y^2 on the path
        temp=0.
        for k in range(N):
            temp+=path[k]**2
        y_squared[i]=temp/N
        #pos_squared[i]=np.average(path**2)

        #Average of Delta y^2 on the path
        temp=0.
        for k in range(N):
            temp+=(path[k]-path[(k+1)%N])**2
        delta_y_squared[i]=temp/N

    #Get rid of non-equilibrium states and decorrelate
    n_corr=1
    n_term=10000

    y_squared=y_squared[n_term:ntimes:n_corr]
    delta_y_squared=delta_y_squared[n_term:ntimes:n_corr]

    return y_squared, delta_y_squared

#MONTE CARLO INTERNAL ENERGY
def U(y_squared,delta_y_squared,eta):

    c1=1./(2*eta)
    c2=1./(2*eta**2)

    av1=np.average(y_squared)
    av2=np.average(delta_y_squared)

    return c1-c2*av2+av1/2.

#THEORETICAL EXPECTATION FOR THE INTERNAL ENERGY
def exact_U(N,eta):

    return 0.5+1./(np.exp(N*eta)-1)
</code>```

The issue

To test the code above, I have fixed $\eta$ and calculated the internal energy for various values of $N$, comparing it with the exact version. This is the code to reproduce my results, and the graph obtained:

```<code>###-PARAMETERS-###
start_path='cold'
ntimes=50000
eta=0.01
delta=np.sqrt(eta)

Nmin=1
Nmax=10

N_array=np.arange(Nmin,Nmax,1)
U_array=[]

start=tm.time()
for N in N_array:

    #Run Metropolis with each N and fixed eta
    y_squared, delta_y_squared=Metropolis_HO(start_path,N,eta,delta,ntimes)

    #Calculate and store the internal energy)
    U_array.append(U(y_squared,delta_y_squared,eta))

plt.plot(N_array,U_array,linestyle='None',marker='s',color='k')

N_array=np.arange(0.1,Nmax,0.1)
exact_U_array=[]
for N in N_array:

    #Calculate theoretical energy
    exact_U_array.append(exact_U(N,eta))

plt.plot(N_array,exact_U_array,'b')
plt.ylim(0,100)
plt.show()

end=tm.time()
print(f"Time taken: {end-start}.")
</code>```

At $\eta=0.01$, I should be expecting the data to fit the theoretical curve much better. How may I improve the code to get better results?

ContentLicense: CC BY-SA 4.0


问:
Optimising transition state geometries using gaussian 16
I'm attempting to carry out transition state optimisations in Gaussian and a lot of my calculations are coming up with errors, the main one being

no analytic 2nd derivatives for this method.

How can I get past this?
I am using the same input as for the local minima, but using the statement opt=(ts,noeigen,calcfc)
Am I perhaps missing something?

Any help would be greatly appreciated! Thank you!


问:
How to calculate Landau levels for a system e.g., Graphene numerically?
I am currently trying to calculate Landau levels for any system numerically in python. But, I am clueless about where to start rather, how to construct the Hamiltonian. Can anyone help?</p>


问:
Crank-Nicolson vs Spectral Methods for the TDSE
The time-dependent Schroedinger equation (TDSE) depends linearly on the system's initial state $\vert \psi(0) \rangle$, such that the solution can be generally written as

$$ \vert \psi(t) \rangle = \exp \left( \frac{i}{\hbar} \int_0^t \text{d}t' \widehat{H(t')} \right) \vert \psi(0) \rangle \quad , \quad (1) $$

where the [possibly time-dependent] Hamiltonian operator is given by $\widehat{H(t)}$.

A common technique to numerically solve the above equation is to use the Zassenhaus formula to approximately split the exponential for an infinitesimal time $\Delta t$. The error one gets depends on how this splitting is performed: simple schemes, such as the famous leapfrog method, have an [energy] error of $\mathcal{O}(\Delta t^2)$. Yoshida, Strang, Feit and others showed that smart splitting schemes can produce methods of arbitrary accuracy, albeit paying the price in computational cost. These splitting methods are known as "spectral methods". After choosing a splitting scheme the final wavefunction is obtained recursively and, for separable Hamiltonians expressible as a sum of kinetic plus potential energies, solved by moving back and forth in momentum and position representations by the repeated use of Fast Fourier Transforms (FFTs).

A much more general scheme is the Crank-Nicholson method, which in principle can be used to solve equations that are a lot more complicated than the TDSE and is a stable, also $\mathcal{O}(\Delta t^2)$ [in energy] method. In general, the solution to these equations cannot be written in the simple form (1), and spectral methods such as the ones described above are not available. In these cases it makes sense to think about Crank-Nicolson, but for the TDSE there is a plethora of super sophisticated/efficient spectral schemes available.

Why, then, do I keep stumbling into papers that employ Cranck-Nicholson to solve the TDSE?


问:
How to scale gradients in a gradient descent algorithm?
I am training a neural network with the multiobjective steepest gradient descent algorithm. The want to steer the direction of the gradient descent so that I land up at a point slightly above where I am landing now. How can I achieve this? Is there any mathematical way to scale the gradients? In case of the neural networks, I am only minimizing the loss function. I am trying to steer the direction of gradient descent by adding some random weights or biases to the loss. How can I mathematically scale the gradients to steer the direction?</p>


问:
SINDy Vs standard methods for system identification
I have been trying to understand the recently proposed Sparse Identification of Nonlinear Dynamics SINDy. Despite several attempts, I seem to fail to understand the difference between SINDy and the standard methods used for parameter identification in the field of system identification. The PNAS paper neither discusses nor compares with any existing methodologies in system identification. Both methods use least-squares to find the parameters. SINDy is claimed to be "data-driven" and based on Machine Learning; however, I do not see where ML comes into the picture. I very much appreciate any insights into this.</p>


问:
Non-Linear Distributed Delayed Kalman Filter
I have a system $\vec{x}_{i + 1} = \vec{x}_i + W_i$ where $W = N(\vec{\mu}, \Sigma)$. For some matrix $H_i$, let $y_i = H_i$ and let $z_i = y_i + R$. Where $R$ is some random variable. We are given $pdf(y_i | z_i)$. Using this, I can calculate the mean and the covariance of $y_i$ conditioned on the observation $z_i$ (I am using numerical integration or Monte Carlo methods). But due to computational performance reasons, the calculated value of the mean and the covariance is not available until time $i + d$.

The Question: What is a good heuristics to update the mean and the covariance of $x_i$?

I attempted to linearize the system by heuristically constructing a fake observation $z'_i$ and $R'_i = N(0, \Sigma_{R'})$ such that mean and covariance of $y_i$ conditioned on $z'_i$ matched the expected value and the covariance calculated in the non-linear system. And then I apply the fake observation using standard Kalman Filter logic. Unfortunately, I have noticed that there is "double counting" due to the delay. When $x_i$ moves significantly $E[y_i|z_i]$ is significantly different from $E[y_i]$. To account for this, my fake observation $z_i$ needs to be large. But since this observation is not accounted for until time $i + d$, at time $i+1$, I will have another extreme value of $z_{i + 1}$. This causes me to over shoot.

I have not used the standard particle filter, because my system is very high in dimension, and the variance of the particle filter estimator is way too high.

Is there a different approach I can try?


问:
How to choose MPI topology for uniform 3D Cartesian grid that is not cubical?
When writing an MPI-parallelised code, I need to split a global Cartesian domain of $(N_x, N_y, N_z)$ cells among $M$ MPI ranks.  For simplicity, we can assume that the decomposition is "simple", and can be expressed as the number of MPI ranks in each direction, $(M_x, M_y, M_z)$, such that $M_x M_y M_z = M$.  How can I choose $M_x$, $M_y$ and $M_z$ such that communication is minimised, where communication is via halos of the subdomains?

This seems to come down to minimising the surface area of the subdomains, which would coincide with having the subdomains as close to cubical as possible.  The MPI function MPI_Dims_create seems to be written for this purpose, but has no knowledge of the number of cells in each direction, so would only minimise communication for cubical domains and not, for example, for a long, thin, channel.

I am looking for a simple, general and possibly suboptimal solution to use as a starting point, but would also be interested in more advanced options.

I did come up with an exhaustive-search algorithm, but it seems quite complicated:

Find the set $P = (p_1, p_2, ..., p_R)$ of all prime factors of $M$

For each $R$-digit base-3 integer (easy to enumerate) $(d_1, d_2, ..., d_R)$, define a topology as $M_x = \Pi_j (pj | dj = 0), M_y = \Pi_j (pj | dj = 1), ... $ (i.e. find all the 3-tuples of products of prime factors)

Exhaustively search all these topologies to minimise the surface area $(A = 2 * (N_x/M_x  N_y/M_y + N_y/M_y  N_z/M_z + N_z/M_z  N_x/M_x))$ of the subdomain

There is also the related problem of what to do when the number of cells in a given direction is not a multiple of the number of ranks; I guess either you change the problem so that it is, or accept that some ranks will have one fewer layer of cells in that direction.


问:
approximate Poisson disk / blue noise sampling
Poisson disk or blue noise sampling consists of finding a random, maximal set of points $\{x_i\}$ in a domain $\Omega$ such that $\|x_i - x_j\| \ge r$ for some fixed radius $r$ and for all $i$, $j$.
There are several algorithms for efficient blue noise sampling, for example Bridson's algorithm.
Blue noise sampling is used in many fields, including computer graphics and mesh improvement.

Many papers use the empirical periodogram for evaluating the quality of the resulting blue noise samples.
This got me thinking:
Is there any published work on generating approximate blue noise samples through Fourier inversion of the ideal periodogram?
The resulting samples might have some pairs $x_i$, $x_j$ of points such that $\|x_i - x_j\| < r$.
But strict adherence to this exclusion condition isn't necessarily a hard requirement for a lot of the things people use blue noise for.
My hunch is that some kind of Fourier or spectral-based approach might (1) make better use of fast DFT algorithms, instead of spending all its time traversing a bunch of linked lists, and (2) be implementable in a single laptop screen of code at 12-point font.
Bridson's algorithm fails on both counts.

I found this paper which is kind of in the right direction -- they calculate that the power spectral density of blue noise is

$$P(k) = 1 - \rho\left(\frac{2\pi r}{k}\right)^{\frac{d}{2}}J_{d/2}(kr)$$

for samples of density $\rho$ in dimension $d$, where $J_\nu$ are the Bessel functions of the first kind.
But this is used more as a means to evaluate algorithms rather than the basis for an algorithm as such.


问:
Numerical precision on tricky coupled nonlinear boundary value problem on infinite interval
I am trying to solve with high precision the following coupled system $(f,h)$ on $[0,\infty]$:

$$-h''-\frac{1}{r}h'+\lambda_c h(f^2-1)+4\lambda_h h^3=0$$
$$-f''-\frac{1}{r}f'+\frac{1}{r^2}f+ f(-\lambda_f+2\lambda_c h^2)+\lambda_f f^3=0$$
where $\lambda_h, \lambda_f, \lambda_c $ are positive constants, and with boundary conditions $f(0)=0, f(\infty)=1, h'(0)=0, h(\infty)=0$.  The condition $h'(0)$ is necessary for $h$ to be finite at the origin.

The issue

The problem is tricky because of the singularity at $r=0$, and the boundary conditions at infinity which preclude a naive shooting approach. Further the mixed boundary conditions on $h$ preclude a naive relaxation method as $h(0)$ is not known.

In the following I'll fix $\lambda_h=\lambda_f=1$. I am having trouble getting scipy's solve_bvp to converge for this system for arbitrary  $\lambda_c$. Solving $f$ by itself by removing the coupling term with $h$ (no "backreaction") works fine. When solving the system with $h$ when $\lambda_c \sim 1 $, the solution converges to 0.1% accuracy, and the results look reasonable.

However, I am having trouble pushing to higher precisions, and the solution does not converge if $\lambda_c$ is too large, and has incorrect behavior for small $\lambda_c$ (setting $\lambda_c=0$, we should have $h=0$, however I have not been able to get that with precision higher than 0.1%).

code

The following snippet produces the below plot for $\lambda_c=1$ which converges and looks reasonable.

```<code>    def solve_h_f_sys(λ_c=1, λ_h=0, λ_f = 1, back_reaction=True, tol=1e-7,nb_rvals=500000,max_rval=100,max_nodes=1000000):

        x = np.logspace(-7,0,nb_rvals)*max_rval
        min_rval = x[0]
        br = int(back_reaction)

        def fun(r,y): # y is (f, rf', h, rh')
            z_h = λ_c*( (y[0])**2 -1) + 4*λ_h*(y[2]**2)

            z_f = λ_f*(y[0]**2-1) + 2*λ_c*(y[2]**2)*br

            return y[1]/r, y[0]/r + r*y[0]*z_f, y[3]/r, r*y[2]*z_h

        def bc(ya,yb):
            return ya[0]-ya[1], ((2/λ_f)**0.5)*yb[1] + max_rval*(yb[0]**2-1), ya[3]/min_rval, yb[2]+yb[3]

        y=[np.tanh(x), x/np.cosh(x)**2, 1/(1+x), -1/ (1+x)]

        #solve
        res=scipy.integrate.solve_bvp(fun, bc, x, y, tol=tol, max_nodes=max_nodes)
        print(res.message,"  ", len(res.x))

        f = res.sol(x)[0]
        h = res.sol(x)[2]
        rms_residuals = res.rms_residuals

    return x, h, f, rms_residuals

x, h, f, rms_residuals = solve_h_f_sys(λ_c=1, λ_h=1, λ_f = 1, back_reaction=True, tol=1e-3,nb_rvals=500000,max_rval=200,max_nodes=1000000)

plt.plot(x,f,label="f")
plt.plot(x,h,label="h")
plt.xlabel("r")
plt.legend()
</code>```

Setting $\lambda_c=1.2$ the solution does not converge and we get:

For small $\lambda_c=0$ I cannot get $h(0)=0$ to be true better than at the 0.1% level. Zoomed in plot for $\lambda_c=0$:

Asymptotic analysis

I find that at small $r$, $f(r) \sim \alpha r $ and $h(r) \sim a +br^2 +cr^4$ for constants $\alpha, a,b,c$. This seems to be respected by the numerical solution. At large $r$, if there is no coupling between $f$ and $h$ then we have $f(r)=1-C / r^2+\mathcal{O}\left(1 / r^4\right)$. Assuming $f$ has the same form when the coupling is turned on, I find the ansatz $h = K/r$ satisfies the first equation if $K^2 = (1+2C\lambda_c)/(4\lambda_h)$. Substituting in the second equation, this requires $C = (1+2\lambda_cK^2)/(2\lambda_f)$. There is a solution $(C,K)$ as long as $\lambda_c^2 < 2\lambda_h\lambda_f$. This suggests the ansatz $h \propto 1/r$ at large $r$ is maybe a poor approximation.

For now I the boundary conditions at infinity I am using are the ones suggested in this answer for  $f$:
$$f^{\prime}(r)^2=\frac{\lambda_{f}}{2}\left(1-f(r)^2\right)^2$$

These should be appropriate even when $f$ couples to $h$ because $h$ must go to zero at infinity so the backreaction term can be neglected. For $h$, I am imposing that $h \propto 1/r$. This is likely incorrect but better than imposing $h=0$ at infinity which does not give correct-looking solutions.


问:
How to take the Fourier transform of a Fibonacci chain in a Python script?
This may be an easy question to answer but I am really stuck.

In several topics (especially that of quasicrystals) the Fibonacci chain's Fourier transform and diffraction pattern is mentioned. Despite the chain not being periodic, its diffraction pattern shows sharp peaks (a hallmark of quasicrystals). See for example in Math SE

Does the cut-and-project method produce the Fibonacci chain?

I have gotten as far as generating the "canonical" Fibonacci chain using the classic substitution of starting with "L", and replacing "S" with "L" and "L with "LS" literally using strings, then converting the short space "S" to 1.0 and the long space "L" to $\tau = (1 + \sqrt{5}) / 2 = 2 \cos(36°)$ and taking the cumulative sum to produce atom positions in a 1D chain.

These atoms would then be point scatterers of an incident plane wave, so we can see them as identical point sources of waves in phase and use $x \sin(\theta)$ to assign a relative phase as seen from infinity at an angle $\theta$.

Below is shown the results for a chain after 20 iterations, producing 17,711 intervals, a long/short ratio of 1.6180340557 nearly converged to $\tau$ ~1.6180339887. The diffraction pattern (for an incident wavelength of unity) for the fibonacci atoms (blue) certainly shows many peaks between 0 and 30 degrees, and the red sequence (random pattern with the same "L" to "S" ratio) doesn't. The first 100 atoms for each are shown in the top plot with random on top.

But if I wanted to simply apply np.fft.fft() without calculating the physical diffraction pattern, how would I do that exactly?

Question: How to take the Fourier transform of a Fibonacci chain in a Python script?

```<code>import numpy as np
import matplotlib.pyplot as plt

s = 'L'
for i in range(20):
    wow = []
    for thing in s:
        if thing == 'L':
            wow.append('LS')
        else:
            wow.append('L')
    s = ''.join(wow)
    print('length: ', len(s), ' L/S ratio: ', s.count('L') / s.count('S'))

tau = (1 + 5**0.5) / 2

g = np.array([thing == 'L' for thing in s]).astype(float) * (tau-1) + 1.
g_rand = np.array([thing >= 1 / (1 + tau) for thing in np.random.random(len(s))]).astype(float) * (tau-1) + 1.

x = np.cumsum(g)
x_rand = np.cumsum(g_rand)

theta_deg = np.linspace(0, 30, 3000)
theta = np.radians(theta_deg)

i = np.complex(0, 1)
k = 2 * np.pi

amplitude = np.exp(-i * k * x * np.sin(theta[:, None])).sum(axis=1)
power = np.abs(amplitude)**2

amplitude_rand = np.exp(-i * k * x_rand * np.sin(theta[:, None])).sum(axis=1)
power_rand = np.abs(amplitude_rand)**2

if True:
    fig, (ax1, ax2) = plt.subplots(2, 1)
    n = 100
    ax1.plot(x_rand[:n], 0.1 * np.ones_like(x_rand[:n]), '.r')
    ax1.plot(x[:n], np.zeros_like(x[:n]), '.b')
    ax1.set_xlim(-2, n)
    ax2.plot(theta_deg, power_rand, '-r', linewidth=0.6)
    ax2.plot(theta_deg, power, '-b', linewidth=0.6)
    ax2.set_yscale('log')
    plt.show()
</code>```


问:
Can you left-precondition least squares?
Suppose I want to solve an overdetermined linear least squares problem

$$
x = \operatorname*{argmin}_{x\in\mathbb{R}^n} \| Ax - b\|^2
$$

where $A \in \mathbb{R}^{m\times n}$ has full column rank. I make no assumptions that this system is nearly consistent.

Suppose I know a have access to a matrix $P\in\mathbb{R}^{n\times m}$ which is a good approximation to the Moore–Penrose pseudoinverse $P \approx A^+$. I only have access to $P$ through  black-box subroutines to compute $y \mapsto Py$ and $z\mapsto P^*z$. Can I use $P$ to compute $x$ rapidly by some kind of iterative scheme?

My instinct is that this is not possible. All of the work I'm preconditioning for least squares I know of either use a right preconditioner (i.e., assume I can apply maps $R$ and $R^{-1}$ such that $AR^{-1}$ is well-conditioned) or use an (approximate) QR factorization of $A$. It makes sense why left preconditioning may not be possible for this kind of problem, but I haven't seen any authoritative source state that this type of left-preconditioning is not possible/not known to be possible.


问:
Difference between minimizing a function by gradient descent and by norm minimization?
I have working on 2 ways of training a neural network.  The first method uses gradient descent updates the model with Adam optimizer. the second method minimizes the norm of the gradient of the functions and updates the network. What is the difference between the two? does normal gradient descent guarantee minimum norm of the cost function? I am using binary cross entropy loss function. How does minimizing the norm of the gradient works?</p>


问:
Quadrature rules for products of 2D regions
I am interested in computing integrals of the form $\iint_{P\times P} Q(x_1,x_2,y_2,y_2) dxdy$ where $P$ is a polygon and $Q$ is a polynomial. The coordinates $(x_1,x_2)$ are in the plane of $P$. Of course $(y_1,y_2)$ represent again coordinates in the plane of $P$, possibly different from previous ones.

The above integral can be computed explicitly in terms of integrals of the form $I(p,q)=\int_P x_1^px_2^q$, but formulas quickly become really ugly (lots of binomial coefficients). There is an obvious strategy:

evaluate $I(p,q)$ numerically (exact quadrature rule) or analytically

find the exact decomposition of $\iint_{P\times P} Q(x_1,x_2,y_2,y_2) dxdy$ into sum of products of terms of the form $I(p,q)$

However, I'm interested to know if other strategies, more efficient and possibly more precise were developed.

Do you know any other efficient methods for computing $\iint_{P\times P} Q(x_1,x_2,y_2,y_2) dxdy$ where $P$ is a polygon and $Q$ is a polynomial?


问:
How to incorporate a continuity of 1st derivatives at an internal interface in a 2nd order Poisson equation in finite elements method?
How should I modify the following matrix formulation of a finite element (FE) approximation of Poisson equation at the (internal) interface denoted by the nodes $i-1$ to $i+1$ having a continuous 1st derivative $\epsilon_1 \nabla\cdot U_1 = \epsilon_2 \nabla\cdot U_2$ at the internal interface nodes (say from $i-1$ to $i+1$)?

$$\begin {bmatrix}
\epsilon_1 K_{11} & \epsilon_1 K_{12} & \epsilon_1 K_{13} & \dots & \epsilon_1 K_{1N}\\
\epsilon_1 K_{21} & \epsilon_1 K_{22} & \epsilon_1 K_{23} & \dots & \epsilon_1 K_{2N}\\
\vdots\\
\epsilon_1 K_{(i-1)1} & \epsilon_1 K_{(i-1)2} & \epsilon_1 K_{(i-1)3} & \dots & \epsilon_1 K_{(i-1)N}\\
\epsilon_1 K_{i1} & \epsilon_1 K_{i2} & \epsilon_1 K_{i3} & \dots & \epsilon_1 K_{iN}\\
\epsilon_2 K_{(i+1)1} & \epsilon_2 K_{(i+1)2} & \epsilon_2 K_{(i+1)3} & \dots & \epsilon_2 K_{(i+1)N}\\
\vdots\\
\epsilon_2 K_{N1} & \epsilon_2 K_{N2} & \epsilon_2 K_{N3} & \dots & \epsilon_2 K_{NN}
\end{bmatrix}
\begin {bmatrix}
u_1\\
u_2\\
\vdots\\
u_{i-1}\\
u_i\\
u_{i+1}\\
\vdots\\
u_N
\end{bmatrix}
=
\begin {bmatrix}
f_1\\
f_2\\
\vdots\\
f_{i-1}\\
f_i\\
f_{i+1}\\
\vdots\\
f_N
\end{bmatrix}$$


问:
What EXACTLY are the "outflow" and "pressure outlet" BCs in ANSYS Fluent?
I am just starting my MSc and we have been given an assignment in which we must model two-dimensional, laminar, incompressible Poiseuille flow in a channel. We are supposed to use a VELOCITY INLET inlet boundary condition and two different outlet boundary conditions - the OUTFLOW boundary condition and the PRESSURE OUTLET boundary condition. The problem is I don't understand exactly what these boundary conditions are. I am dealing with a simple rectangular domain with four sides (Inlet, outlet, and top/bottom wall). The flow moves from left to right. I will make a list below showing which boundary conditions I know and which I do not know.

VELOCITY

Inlet : Dirichlet, $\boldsymbol u=\boldsymbol u_{\text{inlet}}$

Outlet : ???

Top and Bottom Wall : Dirichlet, $\boldsymbol u=0$

PRESSURE

Inlet : ???

Outlet : ???

Top and Bottom Wall : Neumann, $\boldsymbol n\cdot\nabla p=0$

As you can see, I am unsure of the outlet conditions on $\boldsymbol u$, and the inlet AND outlet conditions on $p$. I don't know what these are when Fluent is using the OUTFLOW and PRESSURE OUTLET conditions. I had a look at this thread on CFD online, but it still isn't clear exactly what BCs Fluent is using. And the ANSYS official documentation is absolutely abysmal as well.

As I understand from the linked thread, the OUTFLOW condition assumes a Neumann condition on all flow variables except the pressure. So, what does it do with the pressure?!

And, it seems that the PRESSURE OUTLET boundary condition assumes a Dirichlet condition on the pressure. But what does it do with the velocity?!

Finally, in both of these cases, what boundary condition is the VELOCITY INLET choice assigning for the pressure at the inlet?

Would really appreciate some clarity here whatever rigorous mathematical language is needed to express it.

I asked this question on CFD Online but got no responses.


问:
What is the difference between pareto-stationary and pareto-optimal points?
I am training a neural network with an multi-objective optimization algorithm. The output that I get after training the network can be termed as pareto-stationary or pareto-optimal points?</p>


问:
Mass spectrometer atomic species identification
I'm starting to look at ways to analyze the output of a mass spectrometer. My first searches led to some white papers far from what I'm looking for (complex molecules vs simple atoms for me). With mass spectrometers being over a century old I'd thought it'd be easy to find algorithm to extract data from the spectra, but no luck so far.

So, I have a array of B/I points, I being the intensity which has some peaks I want to match to some atoms with various states of ionization. I also have a list of likely atoms and their possible charges (i.e. species), for instance:

```<code>A;Name;Mass;Abundancy;Z
1;H;1.007825037;99.985;1
2;H;2.014101787;0.015;1
3;He;3.016029297;0.00014;2
4;He;4.00260325;99.99986;2
6;Li;6.0151232;7.5;3
...about 300 species listed
</code>```

The data looks like this:

```<code>7.685483E+1     7.456893E-13
7.685368E+1     9.060985E-13
7.685262E+1     1.510505E-13
7.685148E+1     1.550802E-12
7.685043E+1     2.034891E-12
7.684936E+1     -2.060149E-14
7.684825E+1     1.334606E-12
7.684711E+1     1.191975E-12
...60000 lines like that
</code>```

Now I'm looking for a way to automatically find the species in the scan, and their relative quantities bearing in mind that:

some peaks are very large and reach saturation (can't match with a gaussian)

some peaks are very narrow (one single sample)

some are negative (apparently where the electron beam went through, I'm not clear on that; or when the calibration is off)

there's a wild variation in peak intensity, by several orders of magnitude (even looking at the graph in log scale it's hard to see some minor but present ones)

some species are far more likely than others (N vs U)

it all depends on the charge/mass ratio of the ions via the formula B=a+b*(A/q)^½ where a and B depends on the calibration and should be found by the algo based on some large and always present reference peaks.

Are there standard algorithms for that ? And if not where should I start (I'm not a physicist) ?


问:
Dispersion of Runge-Kutta methods when applied to systems of ODEs
I am interested in computing the dispersion / phase error(s ?) of an (explicit) Runge-Kutta method when applied to a linear system of ODEs

$$ u'(t) = A u(t). \tag{1} \label{1} $$

To begin, consider the linear test equation
$$u_t = \lambda u = (x + iy) u \tag{2}$$
and a $S$-stage Runge-Kutta method with stability polynomial $G_S = \sum_{i=0}^S \beta_i z^i = \sum_{i=0}^S \beta_i (\Delta t \lambda)^i = \sum_{i=0}^S \beta_i \big(\Delta t (x+iy) \big)^i = (P_S + iF_S)$.
Then, the next state $u_{n+1}$ can be computed as
$$u_{n+1} = G_S u_n = (P_S + iF_S) u_n. \tag{3}$$

This paper recalls a definition first presented in this publication:

The Runge-Kutta scheme defined through $G_S$ is dispersive of order $q$ if

$$ d(\Delta t, y, F_S,P_S) := \Delta t y - \tan^{-1} (F_S/P_S) = \mathcal{O}(\Delta t ^{q+1}) \tag{4} \label{4}.$$

Clearly, I can compute the dispersion $d$ for each eigenvalue $\lambda \in \sigma(A)$ \eqref{4}. Doing so, however, leaves me with a range of different values $d$. I would like to make the connection to an "actual observed dispersion".

For the case that \eqref{1} is the semi-discretization of e.g. the linear advection equation
$$u_t + a u_x = 0 \tag{5}$$
there is a very clear understanding what the dispersion of the Runge-Kutta scheme is:
It is the phase error of the true solution (orange) compared to the numerical one (blue):

Is there any known relation of the "observed" dispersion and the dispersion(s) computed for the eigenvalues of the system?


问:
How is the saddle point in a multi-objective optimization different frm a saddle point in a single objective optimization?
I am doing gradient descent with the multiple objectives for example f1 and f2 are two functions to be optimized. Multi-objective optimization outputs the pareto optimal points. How are the saddle points in multi-objective optimization different from the single objective optimization? Can we reach any saddle point through multi-objective optimization of 2 functions simultaneously? If yes then how is it different?</p>


问:
Can I use MOL to solve 2D steady state PDE in terms of r and z spatial coordinates?
recently I need to solve a 2D steady state PDE equation.
It’s not time dependent, and the only two independent variables are z and r direction.

So far for this solution, I was thinking using Method of Lines (MOL) to semi-discretize only the r radius direction and source term, leaving the spatial z coordinate and Temperature to be treated with ODE solver, ex:ODE45.

My question is: Is this the correct way to solve this 2D steady state PDE by MOL?
Since in general MOL is mostly dealing with one time dependent variable and one spatial coordinate in general PDEs.  Thus, I’m not sure if z spatial coordinate in this case, could be considered using implicit way to be solved (just like time step) ??

Any help or advice would be greatly appreciated.


问:
Optimizing noisy but deterministic one-variable functions within an interval
I am looking for advice on what numerical methods to consider to solve optimization problems like the following, ideally with as few evaluations of the objective function as possible.

I am looking to maximize a one-variable function $f(x)$ on an interval $[a, b]$.

The function is a black box and cannot be computed outside of this interval (other than assuming that $f(x) = f(b)$ for any $x > b$).

It is expected to generally have a "one hump" shape. However, it is also noisy so technically it has more than one local maximum. I am looking to get close to the global maximum, expected to be around the peak of the hump.

There tends to be a (noisy) plateau around the maximum.

The smoothened shape is not necessarily convex.

It may reach the maximum at an endpoint (in $a$ and $b$) and these cases should be handled well.

The function consists of steps and flat sections (essentially it is discrete valued as there is an underlying combinatorial structure). Therefore estimating the derivative is not easy.

It is computed deterministically (the "noise" is deterministic).

Here are four example shapes of the function, sampled at high-resolution, to give an idea of what it is like.


问:
Optimizing multivariable without known function (only input & output), but it's a polynomial
I have multiple set of input (5 parameters) and output from an unknown function. But I know this is a polynomial function.
What method can I use to optimize to find the polynomial variable that can fit those parameters? A lot of optimization method that I know need a known function.

I have tried optimizing the  polynomial function below using simulated annealing, but still fails.

```<code>y_pred = p1*p2*p3*p4*w1*x**5 + p1*p2*p3*p4*w2*x**4 + p1*p2*p3*p4*w3*x**3 + p1*p2*p3*p4*w4*x**2 + p1*p2*p3*p4*w5*x + p1*p2*p3*p4*w6
</code>```


问:
Software for Smith form or Hermite form of a sparse polynomial matrix
In a current research project, I have a number of matrices with coefficients in ℚ[] for which I want to understand how their rank depends on the value of the parameter .

These matrices are:

sparse (a typical example is a matrix of size 160x200 with less than 10 non-zero entries per row),

rather 'regular' (each entry is a polynomial of degree at most one), and relatively

close to full rank (in the above example of 160x200 rank is at least 150 no matter what the parameter is).

I would like to find software that can try to speed up calculations knowing at least some of these special features. So far I mostly tried Magma, and the calculations get endlessly long very fast for the sizes of matrices that I need. (I would also appreciate advice on Magma since perhaps I am forgetting some options that one should specify, I just tried the default versions of HermiteForm and ElementaryDivisors).


问:
where to use multiprocessing in multi-loop program
I am running a program that uses [a lot of] looping to perform a type of bootstrap analysis. The computational step - PCA - can take a long time (60+) seconds, which means that collectively this program could take almost a week to run serially. I have tried using multiprocessing on a high performance cluster - using concurrent.futures - but the analysis is just as slow.  I know that misusing multiprocessing can lead to same/worse performance than serial, so I'm wondering at what point in my code this happens.

My analysis program looks like this, serially:

def analysis(dataset, subset_sizes, n_iters):
    for size in subset_sizes:
        for _ in range(n_iters):
            subsampled_dataset = random_sample(dataset, size)
            computation(subsampled_dataset)

for dataset in datasets:
    for _ in range(n_trials):
        analysis(dataset, **kwargs)
```
<p>using multiprocessing, it looks like this:</p>
<pre class="lang-py prettyprint-override"><code>with concurrent.futures.ProcessPoolExecutor() as executor:

    futures = []
    for dataset in datasets: # here, datasets looks like [datasetN_trialN, ...] to use one loop
        futures.append(executor.submit(analysis, **kwargs)) # dataset is part of kwargs

    results = [f.result() for f in concurrent.futures.as_completed(futures)]
</code>```

Should I be calling executor.submit at a different point in my code?


问:
I'm working on a problem to proof the two-grid cycle, with ν1 pre-smoothing and ν2 post-smoothing steps

I'm facing trouble while proving the given formula. Could you please help me with this? Thanks :)


问:
Solving second order coupled differential equations in python
as I have to design a reactor and therefore have to get its length x, I have to solve the following differential equations:

$$D_{eg}\tfrac{d^2A_g}{dx^2}-u_g\tfrac{dA_g}{dx} = k_la_b\left(\tfrac{A_g}{H_A}-A_l\right)$$
$$\frac{}{}D_{el}\frac{d^2A_l}{dx^2}-u_l \tfrac{dA_l}{dx}+ k_la_b\left(\tfrac{A_g}{H_A}-A_l\right) = (1-\epsilon_B)\eta_ck_2A_sB_s$$
$$D_{el}\frac{d^2B_l}{dx^2}-u_l \tfrac{dB_l}{dx} = (1-\epsilon_B)\eta_cK_2A_sB_s$$

Only $\eta_c$, $A_s$ and $B_s$ are dependent on $A_G,A_L$ and $B_L$ and can be solved by the following circular arguments:
$$\eta_c = \tfrac{3A_sB_s}{\phi_T^2A_gB_l}(\phi_Tcoth(\phi_T)-1)$$
$$A_s = A_g\left ( f\left (1+\frac{\phi_{A,s}^2(\phi_Tcoth(\phi_T)-1)}{\phi_T^2Bi_{A,w}}  \right )^{-1}+(1-f)\left (1+\frac{\phi_{A,s}^2(\phi_Tcoth(\phi_T)-1)}{\phi_T^2Bi_{A,d}}^{-1}\right) \right )$$
$$B_s = B_l\left(1+\frac{\phi_{B,s}^2f(\phi_Tcoth(\phi_T)-1)}{\phi_T^2Bi_{B,w}}\right)^{-1}
$$
$$\phi_T = \phi_{A,s}(1-\tfrac{\phi_{B,s}^2}{3\phi_{A,s}^2}) \quad if \quad \phi_{A,s}\geq \phi_{B,s}$$
$$\phi_T = \phi_{B,s}(1-\tfrac{\phi_{A,s}^2}{3\phi_{B,s}^2}) \quad if \quad \phi_{A,s}\leq \phi_{B,s}$$
$$\phi_{A,s} = r_p\sqrt[2]{\tfrac{\alpha k_fB_s}{D_A}}$$
$$\phi_{B,s} = r_p\sqrt[2]{\tfrac{\alpha k_fA_s}{D_B}}$$

The Biot number is dependet on the length of the catalyst bed and must therefore be considered while integrating. The indexes only denouce the independet parameters that are being used:
$$Bi = \frac{\alpha_{Bi}x}{\lambda_{Bi}}$$

For the fractional wetting of A, this leads to:
$$Bi_{A,w} = \frac{\alpha_{A,w}x}{\lambda_{A,w}}$$

I tried to use the method of lines to linearize the differential equations and solve it stepwise. This leads to:

$$\frac{dA_{g,i}}{dx} = \frac{1}{u_g}\left (  D_{eg}\frac{A_g(y_{i+1},x)-2A_g(y_i,x)+A_g(y_{i-1},x)}{\Delta x^2}-k_la_b\left ( \frac{A_{g,i}}{H_A}-A_{l,i} \right )\right )$$
$$\frac{dA_{l,i}}{dx} = \frac{1}{u_l}\left (  D_{el}\frac{A_l(y_{i+1},x)-2A_l(y_i,x)+A_l(y_{i-1},x)}{\Delta x^2}+k_la_b\left ( \frac{A_{g,i}}{H_A}-A_{l,i} \right )-(1-\epsilon_B)\eta_ck_2A_sB_s\right )$$
$$\frac{dB_{l,i}}{dx} = \frac{1}{u_l}\left (  D_{el}\frac{B_l(y_{i+1},x)-2B_l(y_i,x)+B_l(y_{i-1},x)}{\Delta x^2}+(1-\epsilon_B)\eta_ck_2A_sB_s\right )$$

I now tried to solve it after this process: http://hplgit.github.io/prog4comp/doc/pub/p4c-sphinx-Python/._pylight006.html

I am now completely lost, as you can´t seem to pass several starting conditions into the function, the Biot numbers halt the prozess, as they are dependent on x.
The starting conditions at x=0 are: $A_G = A_{G,start}, A_L = 0, B_L = B_{L,start},A_{L,s} = 0, B_{L,s} = 0$. The concentrations for $A_G$ and $B_L$ can be chosen arbitrarily, at the moment.
My code until now is the following:

```<code>import scipy
import numpy as np
import mpmath

L = 2
N = 1000
x = np.array([0,L])
dx = x[1]-x[0]
D_eg = 2
D_el = 2
k_L = 3
a_B = 2
H_A = 10
u_G = 1
epsilon_B = 1
u_L = 1
K_2 = 20
f = 1
r_p = 1
k_r = 1
D_A = 1
D_B = 1
alpha_A_w = 1
L_B = 1
lambda_A_w = 1
alpha_A_d = 1
lambda_A_d = 1
alpha_B_w = 1
lambda_B_w = 1
length_mesh = 1000
A_G_Start = 10
B_L_Start = 10

A_G = np.zeros(length_mesh)
A_L = np.zeros(length_mesh)
B_L = np.zeros_like(A_G)

def odes(x,A_G, A_L, B_L):
    # print(A_G)

    N = len(A_G) - 1
    gas_A = np.zeros(N + 1)
    liquid_A = np.zeros_like(gas_A)
    liquid_B = np.zeros_like(gas_A)
    A_L_surface = np.zeros_like(gas_A)
    B_L_surface = np.zeros_like(gas_A)
    A_L_surface[0] = 0.001
    B_L_surface[0] = 0.0001
    for i in range(1, N):
        Thiele_A_ = Thiele_A(r_p,k_r,B_L_surface[i-1],D_A)

        Thiele_B_ = Thiele_B(r_p, k_r, A_L_surface[i - 1], D_B)

        Bischoff_ = Bischoff(Thiele_A_,Thiele_B_)

        Biot_A_w = Biot(alpha_A_w,x,lambda_A_w)
        Biot_A_d = Biot(alpha_A_d, x, lambda_A_d)
        Biot_B_w = Biot(alpha_B_w, x, lambda_B_w)
        A_L_surface[i] = A_L_s(Thiele_A_,Bischoff_,Biot_A_w,Biot_A_d,A_G[i-1])
        B_L_surface[i] = B_L_s(Thiele_B_,Bischoff_,Biot_B_w,B_L[i-1])

        eta = eta_c(A_L_surface[i],B_L_surface[i],A_G[i-1],B_L[i-1],Bischoff_)

        gas_A[i] = (D_eg / dx ** 2) * (A_G[i + 1] - 2 * A_G[i] + A_G[i - 1]) - k_L * a_B * (A_G[i] / H_A - A_L[i]) / u_G
        liquid_A[i] = ((D_el / dx ** 2) * (A_L[i + 1] - 2 * A_L[i] + A_L[i - 1]) + k_L * a_B * (
                    gas_A[i] / H_A - A_L[i]) - (1 - epsilon_B) * eta * K_2 * A_L_surface[i] * B_L_surface[i]) / u_L

        liquid_B[i] = ((D_el / dx ** 2) * (B_L[i + 1] - 2 * B_L[i] + B_L[i - 1]) - (
                    1 - epsilon_B) * eta * K_2 * A_L_surface[i] * B_L_surface[i]) / u_L
    return gas_A,liquid_A,liquid_B

def eta_c(A_L_s, B_L_s, A_G, B_L, Bischoff):
    return 3*A_L_s*B_L_s/(Bischoff**2*A_G*B_L)*(Bischoff*mpmath.coth(Bischoff)-1)

def Bischoff(Thiele_A, Thiele_B):
    if Thiele_A > Thiele_B:
        return Thiele_A * (1 - Thiele_B ** 2 / (3 * Thiele_A ** 2))
    elif Thiele_B >= Thiele_A:
        return Thiele_B * (1 - Thiele_A ** 2 / (3 * Thiele_B ** 2))

def Thiele_A(r_p, k_r, B_L_s, D_A):
    return r_p * (k_r * B_L_s / D_A) ** 0.5

def Thiele_B(r_p,k_r, A_L_s, D_B):
    return r_p / f * (k_r * A_L_s / D_B) ** 0.5

def A_L_s(Thiele_A, Bischoff, Biot_A_w, Biot_A_d, A_G):
    print( A_G * (
                f * (1 + (Thiele_A ** 2 * (Bischoff * mpmath.coth(Bischoff) - 1)) / (Bischoff ** 2 * Biot_A_w)) ** -1))
    return A_G * (
                f * (1 + (Thiele_A ** 2 * (Bischoff * mpmath.coth(Bischoff) - 1)) / (Bischoff ** 2 * Biot_A_w)) ** -1 +
                (1 - f) * (1 + (Thiele_A ** 2 * (Bischoff * mpmath.coth(Bischoff) - 1)) / (
                    Bischoff ** 2 * Biot_A_d)) ** -1)

def B_L_s(Thiele_B, Bischoff, Biot_B_w, B_L):
    return B_L * (
            (1 + (Thiele_B ** 2 * (Bischoff * mpmath.coth(Bischoff) - 1)) / (Bischoff ** 2 * Biot_B_w)) ** -1)

def Biot(alpha,L,lambdaa):
    return alpha*L/lambdaa

def ode_FE(f, U_0, dt, T):
    N_t = int(round(float(T)/dt))
    # Ensure that any list/tuple returned from f_ is wrapped as array
    f_ = lambda u, t: np.asarray(f(u, t))
    u = np.zeros((N_t+1, len(U_0)))
    t = np.linspace(0, N_t*dt, len(u))
    u[0] = U_0
    for n in range(N_t):
        u[n+1] = u[n] + dt*f_(u[n], t[n])
    return u, t

U0_AG = np.zeros(N+1)
U0_AL = np.zeros(N+1)
U0_BL = np.zeros(N+1)
A_G[0] = A_G_Start
A_L[0] = 0
B_L[0] = B_L_Start

x,A_G,A_L,B_L = ode_FE(odes,[A_G,A_L],dt,T=1.2)
</code>```

I would be incredibly happy and thankful, if you could help me.


问:
Implementations of the Eisenstat-trick for SSOR
Where can I find a source code for the SSOR method with Eisenstat's trick?

The original paper includes pseudo code but also seems to have minor typos. For that reason, I would be very happy to see an implementation somewhere. Unfortunately, I haven't any in the common libraries (perhaps the code is hidden to deep).


问:
How to save multiplication computation time between a dense vector and a not that sparse matrix?
I am trying to compute $\mathbf{X}\mathbf{u}$ for many times in my algorithm, where $\mathbf{X}\in \mathbb{R}^{n\times m}$ and $\mathbf{u} \in \mathbb{R}^{m}$. The problem is that, during the computational process of the algorithm, $\mathbf{X}$ is a constant, however, more and more elements in $\mathbf{X}$ can be changed to zero, which might help us to save computational time. $\mathbf{u}$ is just a variable. This process is kind similar to a process of changing a dense matrix to a band matrix and finally a sparse matrix. But the band matrix is not standard as it might look like a boomerang or something else.

I know that using a sparse matrix to represent the $\mathbf{X}$ can speed up the multiplication, but the advantage only appears when the sparsity is above 90% percent. It is not suitable for me because most of the time my $\mathbf{X}$ is not that sparse, and even when it reached 90% percent sparsity, the speed-up ratio of the Sparse Matrix multiplication in Scipy is not that oblivious.

I use NumPy in python to apply the multiplication $\mathbf{X}\mathbf{u}$, which is really fast. I wonder whether it is possible to be faster as we are sure that more and more elements in $K$ are zeros and unimportant. Can we speed up it by using space to exchange for computational time?


问:
Effect of a flipped image on projection matrix and the intrinsic/extrinsic calibration matrices
I already asked this question on StackExchange Mathematics, but it seems too domain-specific. Searching where to put Computer Vision related questions it was suggested to use Computational Science since there is no specific CV StackExchange branch.

With the direct linear transform, one can calibrate a camera and receive the projection matrix, which can be further decomposed in intrinsic and extrinsic camera parameters.

For this, at least 6 points need to be clicked on the image of which the 3D coordinates are known.

I wonder what is the effect of a flipped image (left to right) on the projection matrix or on the decomposed parameters?

I have tested this with an image. The resulting two p-matrices do not look similar to me. I hoped that only some signs changes but my intuition was wrong. When checking the decomposed matrices it seems like the intrinsics changes the sign for the x component and the skew of x:

$$
\begin{bmatrix}
-1 & -1 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
$$
which makes sense.
The rotation matrix seems to have negative signs and multiple places as well:
$$
\begin{bmatrix}
0 & -1 & -1 \\
-1 & -1 & 0 \\
0 & -1 & -1
\end{bmatrix}
$$
which makes sense to me when I draw the new axis of the camera, which is now no longer a right-handed coordinate system and seems to be mirrored.

In the end I failed to see a connection between the two translations T that I received from the two decomposed projection matrices. Otherwise, I could have combined K,R and T to receive the "for flipping corrected P" with
$$P_c = K_c [R_c|T_c]$$ with c indicating the correction.


问:
Kolmogorov $\epsilon$-entropy, $n$-width, and $\epsilon$-capacity and applications
What is the relationship between Kolmogorov $\epsilon$-entropy, Kolmogorov $n$-width, and Kolmogorov $\epsilon$-capacity of a set $M$ in a metric space $X$? (The $\epsilon$-capacity here is the largest number of points in $M$ whose distance from each other is $\geq \varepsilon$).

Are they are equivalent? If yes, in which sense (i.e. which bounds occur between these three quantities)?

Let us suppose that we have a continuous semigroup $(S_t)_{t\ge 0}$ acting on a Banach space $X$ (say, $X=L^p(\mathbb R)$) that gives solutions to a certain PDE. Suppose that we have a bound from above and from below of one of these quantities (say, the $\epsilon$-entropy) of $S_T(C)$ where $C$ is some subset of $X$ and for any $T>0$. Can we use this information to say anything about the convergence properties of a numerical approximation of the PDE?


问:
How to decouple nodes and cores with MPI?
I have code currently using multiprocessing with joblib to distribute computations from a loop over multiple cores within one cluster node, like so:

from joblib import Parallel, delayed
def function():
    return Parallel(n_jobs=28)(delayed(computation)() for _ in range(n_iter))
save(function()) # save results to some file
...
```
<p>I am limited to 28 cores on a cluster/node but I'd like to run the above task 4 times.  This might strain the multiprocessing load and hinder performance. I've looked into <code>mpi4py</code> which uses the MPI standard to concurrently run code on multiple nodes.</p>
<p>From what I've learned, I want my code, <code>script.py</code> to follow this pattern:</p>
<pre class="lang-py prettyprint-override"><code>from mpi4py import MPI
comm = MPI.COMM_WORLD
my_rank = comm.Get_rank() # number assigned to a node
# below, num_proc would be 5: node 0 receives data,
# and nodes 1-4 doing the work
num_proc = comm.Get_size()

if my_rank != 0:
    comm.send(function(), dest=0)
else:
    for proc_id in range(1, num_proc):
        result = comm.recv(source=proc_id)
        save(result)
</code>```

The script to submit a job with MPI using TORQUE has the general pattern:

#PBS -l nodes=5:ppn=28 # ppn is parts per node, or number of cores
...
mpirun -np $np python script.py
```

But from what I've seen1, $np here should be nodes * ppn.  However, in the loop in script.py above, I want num_proc to correspond to the number of nodes, and n_jobs in function's Parallel call to correspond to ppn.

How do I tell MPI to 'decouple' nodes and cores?  If I submit mpirun -np 5 python script.py, is joblib distributing my computation across 28 cores per node, and mpi4py distributing this function to 4 nodes?

I know there are various binding options with mpirun but based on the docs am not sure which one is appropriate.

References

1 see 'submitting a job' in this cluster user guide


问:
Calculating radial distribution function in a cylindrical pore
I have tried calculating radial distribution functions in a cylindrical pore geometry. The pore has z-axis periodicity and it does not have x and y periodicity. But I got some troubles on that. The main part of my code is like below.

```<code># Calculating distance between the referenced and selected particles.
for i in range(Nref):
    for j in range(Nsel):
        dz = zj - zi
        if abs(dz) > box[t,2]/2.0: dz = dz - box[t,2]*round(dz/box[t,2]) # periodic boundary condition
        dr = np.sqrt(dz**2)
        if dr < RMax:
            BinNumber = int(dr/binr)
            cnt[BinNumber] += 1

Volume = np.pi*pore_radius**2*box[t,2]
BulkDensity = Nsel/Volume # bulk density of selected particle
for i in range(NumBins):
    rIn = binr*i
    rOut = binr*(i+1)
    dN = cnt[i] # number of molecules in the bin
    dv = np.pi*pore_radius**2*(rOut - rIn)
    trdf = dN/dv/BulkDensity/Nref
    rdf[i,1] += trdf                          # rdf
rdf[:,1] /= Nsteps/dt
</code>```

Can we discuss which part is wrong?


问:
What are the most prominent use cases for Java in scientific computing?
What are the most prominent use cases for Java in scientific computing?

I've viewed the domain heavily C/C++, Fortran and possibly Python dominant, but what about Java?


问:
Algorithm for constructing a supercell using an arbitrary integer matrix
Software like Diamond and VESTA allow users to build a supercell by extending the structure along the basis vectors of the latttice. However, I'd like to do something more general and build a supercell from an arbitrary integer matrix transformation. For instance, if I want to convert a primitive hexagonal unit cell into an orthorhombic supercell, I can use the transformation matrix
$$M = \left[\begin{matrix}
 1 & 1 & 0 \\
-1 & 1 & 0 \\
 0 & 0 & 1 \\
\end{matrix}\right]$$
and get the supercell basis $S$ from the original basis $B$ from the matrix multiplication $S = BM$. I'll note I've been considering this problem for arbitrary dimensions, not just 2 or 3.

The strategies I've thought of are:

Using an LU decomposition of $M$ to generate translation vectors along which atoms can be copied (from the diagonal of the U component), then moving the atoms into the supercell basis $S$ by truncating the integer portion of the fractional coordinate relative to the supercell. This works in most cases, but some particular transformations cause different atoms to move into the same position in the supercell, leaving some sites empty. I've worked around this by padding with more translation vectors, but it's a suboptimal solution.

Enumerating the lattice points inside the parallelotope generated by the transformation matrix, and generating copies of the atomic positions with the addition of the lattice points to the atomic coordinates (fractional with respect to the original basis). The coordinates can then be reduced with left division by $M$. The issue here is enumerating points within the parallelotope generated by the transformation matrix. I know there are $\det(M)$ points inside the parallelotope, but I'm not sure how to traverse them. This seems similar to rasterization problems in computer graphics, and I'm wondering if a generalized form of Bresenham's algorithm could work here.

I assume the solution to this problem may also be related to how k-points are generated from an arbitrary mesh. I know I might be trying too hard to optimize this process, but I'm doing so because I want to understand the structure of the problem.


问:
Package Assignment problem to maximize profit
Problem description

We have a graph $G=(V,E)$. $V$ is the set of nodes. $c_{ij}$ is the profit of traveling through edge $(i,j)$. $T=\{1,2,3,...\}$ is the set of discrete time steps.

At each time step, we have a set of packages that we need to transport to their destination. $d^t_{i,j}$ is the maximum number of packages that needs to be transported from $i$ to $j$ at time $t$. $X^t_{i}$ is the total number of available trucks available at the factory $i$ at time $t$. $B^t$ is the set of available trucks at time $t$. Each truck $n$ is associated with a capacity $r_n$.

We assume that at each time step we don't need to assign all the packages to the trucks.

Objective

Assign the packages to trucks such that the profit is maximized

$$ $$

I would like to model the problem mentioned above as a linear programming optimization problem. I have limited experience in LP. Following is my current formulation of the problem. I tried to run a solver. But it can't solve it, even with a simple example.

Can anyone tell me what the problem with the formulation is? I would appreciate any feedback.

$$\max  \sum_{i,j \in V} \sum_{n \in B} x^t_{ij}c^t_{ij}z^t_{in}$$
$$0\le x^t_{ij} \le d^t_{ij} \quad \forall i,j \in V \quad (1)$$
$$0\le x^t_{ij} \le \sum_{n \in B}z^t_{in}r_n \quad \forall i,j \in V \quad (2)$$
$$\sum_{n \in B}z^t_{in} \le 1 \quad \forall i \in V \quad (3)$$
$$z^t_{in} \in \{0,1\}$$

Constraints:

$(1)$ represents that the number of packages to be transferred is not more than the available packages.

$(2)$ ensures that the number of packages doesn't exceed the capacity of the trucks.

$(3)$ to make sure each truck is selected once


问:
How to compare the bias between the two contour plots?
I wish to compare the value of the normalized concentration (c+), between two contour plots (one is experimental and another is a simulation). If I don't have a c+ value from the same point on both plots, how can I compare?  (A few data points and interpolation generate them).  I want to know which turbulence model performs better i.e which left-side plot is more similar to the right side plot.  The SKE2 is a modified k epsilon model. TIA

$$ figure 1: SKE2  vs Experimental$$

$$ figure 2:SKE vs Experimental $$


问:
Boundary conditions for compressible Euler equations
I want to solve the compressible 1D Euler equations numerically.

Theory says that for subsonic inflows, one can prescribe two variables, e.g. pressure and temperature. Density can then be computed with the equation of state, but how is the velocity calculated?

For subsonic outflows only one variable can be prescribed. Also here my question is: how to compute the remaining variables?


问:
Compute orthogonal complement using BLAS / LAPACK
Is there a fast method to compute an orthogonal complement of an arbitrary matrix $U\in\mathbb{R}^{m \times n}$ in BLAS / LAPACK?

Specifically, I want any matrix $V\in \mathbb{R}^{m \times (m - \text{rank}(M))}$ such that
$$
\text{rank}(U|V) = m.
$$

The singular value decomposition methods in LAPACK can do this if called correctly, but if $U\in\mathbb{R}^{m \times (m-1)}$ is full rank, this will take $O(m^3)$, when it should only take $O(m^2)$ (e.g., via Gram-Schmidt orthogonalization of a random vector).

Orthogonality of $V$ is a plus but not required.


问:
Calculate the time taken to run an algorithm on GPU
I want to calculate the total time taken for a fixed code run using an NVIDIA GPU (for instance, Tesla K40). The code has to run 1 million single-bit comparisons. All the comparisons are independent of each other. There is no memory access involved during the execution of the code since the complete data (2 bits per core) are stored in the local registers of the cores at the time of execution. So effectively, the code can run x number of comparison operations in parallel in a single instruction cycle (where x is the maximum number of parallel comparisons that can be executed on the mentioned GPU, I am not sure this will be the number of cores or the number of threads). I want to understand the process of calculating the following in this problem, using the datasheet of NVIDA Tesla K40:

x

time taken for each instruction cycle (for a single-bit comparison)

Total time taken for 1 million single-bit comparisons


问:
What's the best modern algorithm for recursive least squares?
Recursive least squares can be implemented using the Sherman–Morrison formula to avoid resolving, however, have better methods without $n^2$ cost been developed? I'm interested if there is a good randomized approach here.


问:
Difference between Numeric, Combinatorial, and Geometric Computing
In the paper [1], author has discussed a distinction between the 3 types of computations: numeric, combinatorial, and geometric. The author says that Geometric computation is one that has elements of both combinatorial and numeric computation. That is:

Geometric Computing = Numerical + Combinatorial Computing.

Numerical computing is exemplified by the problem of solving linear
systems of equations; graph searching is a typical problem of
combinatorial computing. The convex hull problem for a set of points
is a canonical example of geometric computing: the points are
specified by numbers but the convex hull is essentially a
combinatorial structure (a labeled graph). Linear programming counts
as geometric computing but solving linear systems does not qualify.

I am confused about the last sentence: "Linear Programming counts as geometric computing but solving linear systems does not qualify." Why is that so? We have a perfectly geometric representation of solving linear systems (through intersection of lines and planes). Why is Linear Programming geometric then?

So, my first question is: Why is LP geometric while solving linear systems is not?

Now, in book chapter [2], the author further clarifies:

The mere presence of combinatorial structures in a numerical
computation does not make a computation “geometric.” There must be
some nontrivial consistency condition holding between the numerical
data and the combinatorial data. Thus, we would not consider the
classical shortest-path problems on graphs to be geometric: the
numerical weights assigned to edges of the graphs are not restricted
by any consistency condition. Note that common restrictions on the
weights (positivity, integrality, etc.) are not consistency
restrictions. But the related Euclidean shortest-path problem (Chapter
31) is geometric.

So, this "consistency condition" seems to be the key distinction to determining if a problem is geometric or not. Could someone shed more light on these consistency conditions?

So, my second question is: What are these consistency conditions that makes a problem geometric?

References:

[1] Yap, Chee-Keng. “Towards Exact Geometric Computation".” Computational Geometry, 1997, 21.

[2] Sharma, Vikram, and Chee K Yap. “Chapter-45 ROBUST GEOMETRIC COMPUTATION,” Handbook of Discrete and Computational Geometry. Third edition. Boca Raton: CRC Press, 2017.


问:
spherical harmonics perturbation
Is there a way to use spherical harmonics to represent a perturbation on a sphere that has a golf ball dimple pattern? What I have write now looks very asymmetric:

With the following expression:
$${\displaystyle Y_{\ell }^{m}(\theta ,\varphi )={\sqrt {{\frac {(2\ell +1)}{4\pi }}{\frac {(\ell -m)!}{(\ell +m)!}}}}\,P_{\ell }^{m}(\cos {\theta })\,e^{im\varphi }} $$
where, $\ell=16$ and $m=4$.

Title: spherical harmonics perturbation


问:
Starting point to use Trilinos linear algebra packages?
I have experience with classical linear algebra packages in C++ like Eigen, Blaze, etc.
I have never wrote my own PetSC back-end solver but I used/modified several of them.

In a new project I would like to use Tpetra with CUDA acceleration (by kokkos?) for a dense linear algebra problem.

The initial version of my code would be just a single node (multi-core CPU + GPU)implementation.

I look at tutorials of Kokkos and Tpetra but couldn't find how to use them together with CUDA and OpenMP support.

Is there any tutorial, sample codes, suggestion to start use them with CUDA.


问:
Geant4, simulate photon emitted into lead
I have just learned Geant4 and recently I tried to simulate photons emitted into lead. After I collected the data and figured out the transmission of photons. I found my data is bigger than the transmission data from XCOM.

Here's my PhysicsLists

```<code>
void PhysicsList::ConstructParticle()
{

  G4BosonConstructor  pBosonConstructor;
  pBosonConstructor.ConstructParticle();

  G4LeptonConstructor pLeptonConstructor;
  pLeptonConstructor.ConstructParticle();

  G4MesonConstructor pMesonConstructor;
  pMesonConstructor.ConstructParticle();

  G4BaryonConstructor pBaryonConstructor;
  pBaryonConstructor.ConstructParticle();

  G4IonConstructor pIonConstructor;
  pIonConstructor.ConstructParticle();

  G4ShortLivedConstructor sLivedConstructor;
  sLivedConstructor.ConstructParticle();

}

#include "G4PhysicsListHelper.hh"

#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"

#include "G4eMultipleScattering.hh"

#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"

void PhysicsList::ConstructProcess()
{
  AddTransportation();
  G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();

  auto particleIterator=GetParticleIterator();
  particleIterator->reset();
  while( (*particleIterator)() ){
    G4ParticleDefinition* particle = particleIterator->value();
    G4String particleName = particle->GetParticleName();

    if (particleName == "gamma") {

      ph->RegisterProcess(new G4RayleighScattering, particle);
      ph->RegisterProcess(new G4PhotoElectricEffect, particle);
      G4ComptonScattering* cs   = new G4ComptonScattering;
      cs->SetEmModel(new G4KleinNishinaModel());
      ph->RegisterProcess(cs, particle);
      ph->RegisterProcess(new G4GammaConversion, particle);

    } else if (particleName == "e-") {

      ph->RegisterProcess(new G4eMultipleScattering, particle);
      ph->RegisterProcess(new G4eIonisation,         particle);
      ph->RegisterProcess(new G4eBremsstrahlung,     particle);

    } else if (particleName == "e+") {

      ph->RegisterProcess(new G4eMultipleScattering, particle);
      ph->RegisterProcess(new G4eIonisation,         particle);
      ph->RegisterProcess(new G4eBremsstrahlung,     particle);
      ph->RegisterProcess(new G4eplusAnnihilation,   particle);

    } else if (particleName == "mu+" ||
               particleName == "mu-"    ) {

      ph->RegisterProcess(new G4MuMultipleScattering(), particle);
      ph->RegisterProcess(new G4MuIonisation(), particle);
      ph->RegisterProcess(new G4MuBremsstrahlung(), particle);
      ph->RegisterProcess(new G4MuPairProduction(), particle);

    } else if( particleName == "proton" ||
               particleName == "pi-" ||
               particleName == "pi+"    ) {

      ph->RegisterProcess(new G4hMultipleScattering(), particle);
      ph->RegisterProcess(new G4hIonisation(), particle);
      ph->RegisterProcess(new G4hBremsstrahlung(), particle);
      ph->RegisterProcess(new G4hPairProduction(), particle);

    } else if( particleName == "alpha" ||
               particleName == "He3"    ) {

      ph->RegisterProcess(new G4hMultipleScattering(), particle);
      ph->RegisterProcess(new G4ionIonisation(), particle);
      ph->RegisterProcess(new G4NuclearStopping(), particle);

    } else if( particleName == "GenericIon" ) {

      ph->RegisterProcess(new G4hMultipleScattering(), particle);
      G4ionIonisation* ionIoni = new G4ionIonisation();
      ionIoni->SetEmModel(new G4IonParametrisedLossModel());
      ph->RegisterProcess(ionIoni, particle);
      ph->RegisterProcess(new G4NuclearStopping(), particle);

    } else if ((!particle->IsShortLived()) &&
               (particle->GetPDGCharge() != 0.0) &&
               (particle->GetParticleName() != "chargedgeantino")) {

      ph->RegisterProcess(new G4hMultipleScattering(), particle);
      ph->RegisterProcess(new G4hIonisation(), particle);
    }
  }
}

</code>```

And here's the chart of my data and data from XCOM.

Did I miss any physics processes? Or is there any ways to collect transmission energy? If someone could guide me to understanding how one does determine this I would be eternally grateful.


问:
How to generate coordinate points of a smallcircle on earth
I am looking up celestial navigation, and according to
https://youtu.be/-ARXW8InStY?t=3320
a specific sun angle reading (sun angle above the horizon) will be the same on a small-circle with the centre point where the sun is at its zenith on earth.

I would like to convince myself about that using Stellarium, but I would need an Ocatve code to generate longitudes and latitudes of points of a small-circle given a radius and a centre point.

Could someone help with that, or sketch how I could progress this, sadly I know almost nothing about spherical trigonometry.


问:
Instability due to different lumping dimensions in second-order linear system
I have a problem that consists of a linear hyperbolic differential equation with some boundary and initial conditions which leads after discretization in space using finite elements to the following second-order system:

\begin{equation}
\boldsymbol M \ddot {\boldsymbol d} + \boldsymbol C \dot {\boldsymbol d} + \boldsymbol {Kd} = \boldsymbol F \qquad (1)
\end{equation}

where

\begin{equation}
\boldsymbol M = \begin{bmatrix} \boldsymbol M_{11} & \boldsymbol 0 \\ \boldsymbol 0 & \boldsymbol 0 \end{bmatrix} \qquad
\boldsymbol C = \begin{bmatrix} \boldsymbol C_{11} & \boldsymbol 0 \\ \boldsymbol 0 & \boldsymbol C_{22} \end{bmatrix} \qquad
\boldsymbol K = \begin{bmatrix} \boldsymbol K_{11} & \boldsymbol K_{12} \\ \boldsymbol K_{12}^\top & \boldsymbol K_{22} \end{bmatrix}.
\end{equation}

All of the $\boldsymbol M_{ii}$, $\boldsymbol C_{ii}$ and $\boldsymbol K_{ij}$ are symmetric;
$\boldsymbol M_{11}$ is positive definite and has a mass-matrix-like structure $\int_{\Omega^e} N_aN_b d\Omega$ on the element level;
$\boldsymbol C_{11}$ is positive semi-definite and $\boldsymbol C_{22}$ is positive definite, and they both have a mass-matrix-like structure $\int_{\Gamma^e} N_aN_b d\Gamma$ on the element level; I'm not sure about the positive-definiteness of $\boldsymbol K$.

I wish to solve this system using the Newmark explicit method, where the solution stage is
\begin{equation}
(\boldsymbol M+\gamma\Delta t \boldsymbol C + \beta \Delta t^2\boldsymbol K)\boldsymbol a_{n+1} = \boldsymbol F_{n+1} - \boldsymbol C \tilde {\boldsymbol v}_{n+1} - \boldsymbol K \tilde {\boldsymbol u}_{n+1} \qquad (2)
\end{equation}
where $\beta=0$, $\gamma=0.5$, and the tildes indicate predicted values. The corresponding implicit scheme is stable (all the eigenvalues of the polynomial eigenvalue problem corresponding to (1) are in the left half plane). Still, the explicit scheme, with lumped $\boldsymbol M$ and $\boldsymbol C$, is not. Note that while $\boldsymbol M$ is singular, the term $\boldsymbol M+\gamma\Delta t \boldsymbol C$ is not, and allows lumping and inversion of the left-hand side of (2).

The finite element discretization is carried out by linear elements, row-sum lumping is applied here for both $\boldsymbol M$ and $\boldsymbol C$.
However, while the lumping of $\boldsymbol M$ is carried out inside the element $\Omega^e$, the lumping of $\boldsymbol C$ is on the element boundary $\Gamma^e$, so the lumped element matrices for $\boldsymbol M$ and $\boldsymbol C$ are different, as the latter does not involve the element dimension perpendicular to $\Gamma^e$. I suppose that this is what caused the instability.

Has anyone encountered instabilities caused by a dimensional mismatch in lumping?
What can be done to remedy the situation?


问:
Discrete laplacian 9 point
I am trying to write a code for 9 point discrete laplacian. I would like to write a matrix and solve the linear system $AU=F$ with gradient conjugate method.
I wrote the matrix this way

```<code>I=speye(M-1,M-1);
E=sparse(2:M-1,1:M-2,1,M-1,M-1);%matrice con 1 sulla sotto diagonale
T=E+E';
B=-1/6*(T+4*I);
D=1/3*(2*T-10*I);
A=(-kron(I,D)+kron(T,B));
</code>```

And I have a matrix $(M-1)^2(M-1)^2$. Then I wrote F like so

```<code>F=zeros(M-1);
for i=1:M-1
    for j=1:M-1
        F(i,j)=(dx^2)/12*(8*f(x(i+1),y(j+1))+f(x(i+2),y(j+1)+f(x(i),y(j+1))+...
            f(x(i+1),y(j+2))+f(x(i+1),y(j))));
    end
end
for i=1:M-1
    F(1,i)=F(1,i)+2*term(xmin,xmin+i*dx,scelta)/3+term(xmin,xmin+(i+1)*dx,scelta)/6+...
        term(xmin,xmin+(i-1)*dx,scelta)/6;
    F(i,1)=F(i,1)+2*term(xmin+i*dx,xmin,scelta)/3+term(xmin+(i+1)*dx,xmin,scelta)/6+...
        term(xmin+(i-1)*dx,xmin,scelta)/6;
    F(end,i)=F(end,i)+2*term(xmin+M*dx,xmin+i*dx,scelta)/3+term(xmin+M*dx,xmin+(i+1)*dx,scelta)/6+...
        term(xmin+M*dx,xmin+(i-1)*dx,scelta)/6;
    F(i,end)=F(i,end)+2*term(xmin+i*dx,xmin+M*dx,scelta)/3+term(xmin+(i+1)*dx,xmin+M*dx,scelta)/6+...
        term(xmin+(i-1)*dx,xmin+M*dx,scelta)/6;
end
F(1,1)=F(1,1)-term(xmin,xmin,scelta)/6;
F(1,M-1)=F(1,M-1)-term(xmin,xmin+M*dx,scelta)/6;
F(M-1,1)=F(M-1,1)-term(xmin+M*dx,xmin,scelta)/6;
F(M-1,M-1)=F(M-1,M-1)-term(xmin+M*dx,xmin+M*dx,scelta)/6;

F=reshape(F',[(M-1)^2,1]);
</code>```

where f is the source term.
I don't get why the error gets bigger when reducing by half the dx.

Title: Discrete laplacian 9 point


问:
How is the saddle point in a multi-objective gradient descent?
I am employing multi-objective gradient descent to solve a problem with the help of the neural network. I would like to know how does saddle point in a multi-objective optimization differs from a saddle point in a single-objective optimization. Does the definition of saddle point changes in a multi-objective optimization?</p>


问:
Solving compressible euler equations in non-conservative form
I am trying to solve the following compressible 1D system of equations in two non-conservative form PDEs for $\rho, Y$ and one ODE for $v$

\begin{align}
    \rho_t+(v+Q)\rho_x &= -\rho q(Y,T,\rho)-\nu\rho(Y^s-Y)H(Y-Y^s) \\
    Y_t+(v+Q)Y_x &= -\nu(Y^s-Y)H(Y-Y^s) \\
    v_t &= -\frac{1}{(\int_0^1\rho dx)}\Big[\int_0^1Q_t\rho dx+\int_0^1(v+Q)q\rho dx-\frac{1}{3}\frac{1}{Fr^2}\int_0^1(\rho-1) dx\\ &-\nu\int_0^1(v+Q)\rho(Y^s-Y)H(Y-Y^s) dx+\xi(v+Q)|v+Q|\Big] \nonumber\\
    T &= \frac{1}{\rho R(y)}
\end{align}
where $\rho,T, u$ are density, temperature and velocity of a flow. And $q$ is a known function of $\rho,T, Y$ and Y is the vapor content of the flow. And $v(0)=v0$
\begin{equation}
    u(x,t) = v(t)+\int_0^x q(Y(y,t),T(y,t),\rho(y,t))dy:=v(t)+Q[q](x,t)
\end{equation}
$Y^s=Y^s(T)$ is a given saturation function which depends on $T$.
And $Fr^2, \xi, \nu$ are scalers. And $H(Y-Y^s)$ is the heaviside function which takes the value 1 if $Y>Y^s$ otherwise its $0$.

The above system is subject to the following initial conditions, for $f$ is given
\begin{align}
        \rho(x,0) = f(x),\quad Y(x,0)=Y_0,\quad v(0)=v_0
\end{align}
And following inflow boundary conditions
\begin{align}
    \rho(0,t) &= 1 \quad \text{ if } u(0,t)>0 \\
    \rho(1,t) &= \rho_{top} \quad \text{ if } u(1,t)<0\\
    Y(0,t) &= 1 \quad \text{ if } u(0,t)>0\\
    Y(1,t) &= Y_{top} \quad \text{ if } u(1,t)<0
\end{align}

I've tried to find the solution to the problem by solving two PDEs by the explicit upwind method and ODE by Euler's method. My problem is the density $\rho$ seems to be blow up during my upwind simulations. Have I done some mistake in writing the upwing method (I will show that part of the code here.)

```<code>for n=2:n_t
    for j=2:n_x
        if v(n) > 0
            rh(1,n) = rh_bot;
            Y(1,n) = y_bot;
            T(1,n) = T_bot;

            % running from upwind direction
            rh(j,n) = rh(j,n-1)-(v(n)+Q(j,n))*dt/dx*(rh(j,n-1)-rh(j-1,n-1))-dt*rh(j,n-1)*q(j,n-1)-dt*conden(Y(j,n-1),Ys(j,n-1),nu)*rh(j,n-1);
            Y(j,n) = Y(j,n-1)-(v(n)+Q(j,n))*dt/dx*(Y(j,n-1)-Y(j-1,n-1))-dt*conden(Y(j,n-1),Ys(j,n-1),nu)*(1-Y(j,n-1));

            R = cal_R(Ra,Rw,Y(j,n));
            T(j,n) = 1/(rh(j,n)*R);
            Ys(j,n) = antonie(T(j,n),Tr,A,B);
            q(j,n) = q_heat_ener_sat(nu,omg,zeta,gammr,cvw,cva,Rw,Ra,eta*E_sol,Y(j,n),T(j,n),rh(j,n),Ys(j,n),T_wall(j));

        else
            rh(n_x,n) = rh_top;
            Y(n_x,n) = y_top;
            T(n_x,n) = T_top;

            % running from down wind direction
            rh(n_x-j+1,n) = rh(n_x-j+1,n-1)-(v(n)+Q(n_x-j+1,n))*dt/dx*(rh(n_x-j+2,n-1)-rh(n_x-j+1,n-1))-dt*rh(n_x-j+1,n-1)*q(n_x-j+1,n-1)-dt*conden(Y(n_x-j+1,n-1),Ys(n_x-j+1,n-1),nu)*rh(n_x-j+1,n-1); % modified equation with extra q term
            Y(n_x-j+1,n) = Y(n_x-j+1,n-1)-(v(n)+Q(n_x-j+1,n))*dt/dx*(Y(n_x-j+2,n-1)-Y(n_x-j+1,n-1))-dt*conden(Y(n_x-j+1,n-1),Ys(n_x-j+1,n-1),nu)*(1-Y(n_x-j+1,n-1));

            R = cal_R(Ra,Rw,Y(n_x-j+1,n));
            T(n_x-j+1,n) = 1/(rh(n_x-j+1,n)*R);
            Ys(n_x-j+1,n) = antonie(T(n_x-j+1,n),Tr,A,B);
            q(n_x-j+1,n) = q_heat_ener_sat(nu,omg,zeta,gammr,cvw,cva,Rw,Ra,eta*E_sol,Y(n_x-j+1,n),T(n_x-j+1,n),rh(n_x-j+1,n),Ys(n_x-j+1,n),T_wall(n_x-j+1));
        end
    end
end
</code>```


问:
Linear PDE solution with constraints
Consider the following linear PDE:

$$\nabla_q V(q) - M_d(q)M^{-1}(q)\nabla_q V_d(q) = 0,$$

where $V(q)$ and $M(q)$ are known and $M_d(q)$ is a grey box function (e.x., $M_d(q)$ is fitted using a neural network). The goal is to find $V_d(q)$.

I am able to find a solution for $V_d(q)$ using physics informed neural networks. But I need to add some constraints on $V_d(q)$, e.g., $V_d(q)$ is quasi-convex with a global minimum at $q^*$. Instead of using a simple MLP for $V_d(q)$, I also tried Input Convex Neural Networks, Input Invex Neural Networks, and hand-crafted loss for the constraints. In all of the cases, the neural network approach does not perform well, either because the constraints are not satisfied or the PDE is not solved.

I also wonder if I can find a solution numerically since all of the other components of this PDE are essentially known and it is just a simple linear PDE.

What are your thoughts on this?

Edit

I was reluctant to post the whole problem here because I wanted to keep the question simple. However, it seems that additional details may help straighten things up. The PDE above comes from this paper, which
is a method for stabilizing a class of underactuated mechanical systems. In the development of the controller, the following nonlinear PDEs need to be solved,
$$G^\perp\left(\nabla_q H - M_d M^{-1}\nabla_q H_d + J_2 M_d^{-1} p\right) = 0, \quad (2.9)$$
where $G^\perp$ is a full rank left annihilator of the control input matrix $G$. In (2.9), $M_d(q)$ and $V_d(q)$ are unknown, and $J_2(q, p)$ is a free parameter. The solution to (2.9) may seem complicated but it can be separated into terms that depend on $p$ and terms that are independent of $p$ as follows,
$$G^\perp\left[\nabla_q(p^\top M^{-1} p) - M_d M^{-1}\nabla_q(p^\top M_d^{-1} p) + 2 J_2 M_d^{-1} p\right] = 0, \quad (2.11)$$
$$G^\perp\left(\nabla_q V - M_d M^{-1}\nabla_q V_d\right) = 0. \quad (2.12)$$
My understanding is that though the solutions of (2.11) and (2.12) also satisfies (2.9), the reverse is not true. Note that (2.11) is still a nonlinear PDE and (2.12) is a linear PDE. So supposedly (2.11) is harder to solve, which is also noted by the authors. Moreover, the authors managed to solve (2.11) and (2.12) with constraints on $V_d$ in the examples presented by the papers. So at least in these simple situations, (2.11) and (2.12) can indeed be solved and I believe the solutions are not unique (due to the fact that $J_2$ is a free parameter.

My problem is the following. I wanted to solve these PDEs using Neural Networks. From my experiments, (2.9) can be easily solved if there's no constraint on $V_d$. But we need the constraint on $V_d$ in order for this whole method to work. And by introducing some new cost to the loss, I couldn't get a satisfactory solution for $V_d$.

Now, I also tried to solve (2.11) using NN and I can always find a solution since there's not extra constraint. But substituting this obtained $M_d$ into (2.12) and try to solve for $V_d$ with constraint does not work. In fact, even without adding constraints to $V_d$, (2.12) still cannot be solved with an obtained "correct" $M_d$ from solving (2.11). I think it makes sense because the co-solution of (2.11) and (2.12) is only a subset to the solution of (2.9).

Title: Linear PDE solution with constraints


问:
Faster convergence for minimizing least squares of forward modelling problems
This specific question was raised from optimizing parameters of column experiments in the hydrogeological context. I want to optimize a parameter of interest (in this case $D$), based on experimental data (which is a breakthrough curve, an array of $c$ at the efflux-end of the column with corresponding $t$) and some forward modelling results.

The governing equation within the domain of interest (the column) is this advection-dispersion equation for solute transport (more source and sink terms like reaction/degradation/sorption are omitted for simplicity):
$$\frac{\partial c}{\partial t} = D\frac{\partial^2 c}{\partial x^2} - u\frac{\partial c}{\partial x}$$

Forward modelling was performed with different values of $D$, and the corresponding outcome (as time series of $c$ at the last cell of the modelled domain (Neumann boundary condition at the efflux-end)) was compared with the experimental data. In the end the sum of squared residuals (SSR) between the numerical and the empirical data was calculated.

For optimizing $D$, my approach is:

1). set initial boundaries for $D$, use a couple of values within the ($D$) boundaries and perform the forward modelling respectively.

2). narrow the boundaries of $D$ by selecting those two with the 3rd and 2nd lowest SSR, and reiterate the step 1).

This approach works (like shown in the figures below), but I suppose not efficient enough. How can I reach the convergence of SSR faster in this forward modelling problem (I have no analytical solutions)? Thanks.


问:
How to exploit QR factorization implicitly
I meet a problem when I try to develop an iterative method for discrete inverse problem
$$Ax+e=b$$
where $A\in\mathbb{R}^{m\times n}$ and $e$ is a noise. I want to approximate the true solution $x_{true}$ by using a nonlinear regularization term $R(x)$. Since $m,n$ is very large, the optimal regularization parameter is very hard to choose, and an iterative regularization method is needed.

My idea is using a iterative reweighted process to approximate $R(x)$ at each step $k$. Suppose at step $k-1$, we have computed $x_{k-1}$, then at step $k$, we construct a regularization matrix $L_k$ to approximate $R(x)$, i.e. let $\|L_kx_{k-1}\|_2=R(x_{k-1})$ (obviously $L_k$ depends on $x_{k-1}$ as well as the particular expression of $R(x)$).

Now the problem is to solve $Ax+e=b$ with changed regularization matrix $L_k$ (may be not square). Following the paper "Kilmer M E, Hansen P C and Español M I 2007: A projection-based approach to general-form Tikhonov regularization, SIAM J. Sci. Comput. 29, 315–30", I want to modify the original joint bidiagonalization algorithm by using changed matrix pairs $\{A,L_k\}$ at each step, which can exploit the properties of $L_k$ into $A$ by using changed preconditioner $R_{k}^{-1}$ of QR factorization
$$\begin{pmatrix}
    A \\L_k
\end{pmatrix}=Q_kR_k.$$

This leads to the problem that $R_{k}^{-1}v$ and $AR_{k}^{-1}v$ need to be computed at each iteration. However, the explicit computation of QR and $R_{k}^{-1}$ is unaffordable since $A$ and $L_k$ is very large.

A simple consideration is to solve the problem
$$\begin{pmatrix}
    A \\L_k
\end{pmatrix}z_k = Q_kv$$
by an iterative linear solver to approximate $z_k=R_{k}^{-1}v$, but $Q_k$ is still unknown without an explicit QR factorization.

So is there any iterative method to approximate $R_{k}^{-1}v$ implicitly?


问:
Given a set of 1d-points, find the most probable periodicity that models the points (with possible omissions) as equidistant occurences
I try to detect interference fringes in a bunch of pictures. I projected on one axis, and I was able to detekt the peaks that indicate one of the fringes.
So now I'm having a list with points (e.g. $(100, 165, 220, 275, 341, 461, 519)$.
In the choosen example, the sought for periodicity is 60, the points scatter around $100 + 60*n$, and the 6th points was omitted).

Now the problem is how to detect this most probable periodicity of $60$? If I knew the missing points, I could simply fit with a model $m*n + b$. This fit would yield approximately $m \approx 60$ and $b \approx 100$. But I don't know these missing points, and that's why I'm looking for a better way/ better cost-function to minimize.

How would one detect this periodicity?


问:
adjoint method for degenerate problems
For the sake of the rest of the question, I'm interested in the porous medium equation

$$S\frac{\partial\phi}{\partial t} = \nabla\cdot K\phi\,\nabla\phi$$

where $S$ and $K$ are spatially-variable coefficients and $\psi$ (the hydraulic head) is the field to be solved for.
This is a parabolic PDE, but it's degenerate since the "diffusivity" $K\phi$ of the medium can go to zero whenever $\phi \to 0$.
As a consequence, the $\phi = 0$ contour propagates with finite speed in contrast to the usual infinite propagation speed of disturbances for non-degenerate parabolic PDE.

There are well-posedness results for the porous medium equation and convergent discretization schemes that don't regularize away the degeneracy.
So in practice the weirdness of a diffusion coefficient going to zero can be overcome.

Now suppose we want to estimate the hydraulic conductivity $K$ from observational data.
In order to do that, we'll probably use the adjoint method.
The adjoint of the linearization of the PDE is

$$S\frac{\partial\psi}{\partial t} = -\nabla\cdot K\phi\,\nabla\psi + K\nabla\phi\cdot\nabla\psi + \ldots$$

where the ellipses denote some source terms that don't particularly matter.
This is solved for the adjoint state $\psi$ backwards on some interval $[0, T]$ with a final condition imposed on $\psi|_{t = T}$.
The adjoint linearization is an advection-diffusion equation, again with a diffusion coefficient ($K\phi$) that can and probably does go to zero somewhere.

Does the degeneracy of the adjoint equation pose a problem in practice?
I've seen papers on overcoming the degeneracy of the forward problem, but I haven't seen anyone discuss the impact (if any) on the inverse problem.
I'd add that degenerate problems show up in more fields than just hydrology -- you get similar equations for ice flow and wildfire spread.

Title: adjoint method for degenerate problems


问:
Householder Vector algorithm in Golub and Van Loan
(This is repost of a question first asked on Mathematics. Hopefully there are more people here who have a copy of Golub and Van Loan to hand)

In the 4th edition of "Matrix Computations", Golub and Van Loan present "Algorithm 5.1.1 (Householder Vector)". The first couple of lines (translated into MATLAB-syntax) read:

```<code>m = length(x); sigma = x(2:m)'*x(2:m); v = [1; x(2:m)];
if sigma == 0 && x(1) >= 0
  beta = 0;
elseif sigma == 0 && x(1) < 0
  beta = -2;
else
  ...
</code>```

The else clause handles the case where sigma is nonzero and no code after the provided snippet modifies v if sigma is zero. The matrix form of the resulting Householder transformation is $I-\beta v v^T$.

The elseif clause is a little strange. It doesn't appear in the algorithm listing in the
3rd edition, so it was added for the 4th, presumably for numerical stability. However, it seems to me to generate a $(\beta, v)$ pair that doesn't map to an orthogonal matrix. For example, if x = [-1; 0; 0] then sigma == 0 and x(1) < 0 so we get beta = -2 and v = [1; 0; 0] and a Householder transformation of [3, 0, 0; 0, 1, 0; 0, 0, 1] which is not an orthogonal matrix.

So my questions are:

What benefit is there to handling that case separately, rather than just setting beta = 0 if sigma = 0?

Does the resulting Householder transformation need to be applied differently?


问:
How to numerically integrate Kepler Problem?
I tried to solve a simple Kepler Problem numerically.

I have discrete time steps, a starting position $(x_0,y_0)$ and starting velocity $(u_0, v_0)$.

I used this iteration by calculating the forces and apply Newton's law:

$$u_{i+1} = u_{i}+ F_x/m \cdot \Delta t$$
$$v_{i+1} = v_{i}+ F_y/m \cdot \Delta t$$

$$x_{i+1} = x_{i}+ (u_{i+1}+u_i)/2 \cdot \Delta t$$
$$y_{i+1} = y_{i}+ (v_{i+1}+v_i)/2 \cdot \Delta t$$

I used mean velocity between last and actual time to calculate new (x, y).

Naively, I thought, that by making the timestep small enough this gives good results.
But I observed, that this method is completely bad, because the total energy tends to increase ant what I get is not what I expect:

What method for integrating such problem is the standard? My idea is obviously complete garbage because it is  numerically most unstable. When I calculate total energy it is not constant at all, even when I make the timestep incredible small, I get totally dissatisfying results.

OwnerDisplayName: MichaelW


问:
Model the water vapor fraction of a humid air inside a tower
I am going to model the moisture air inside an energy tower, therefore. I formulated the following equation to calculate the fraction of water vapor, Y in the moist air, where $\rho_w$ is the water vapor density,$\rho$ is the total density $T$ is the temperature, $R_a, R_w$ are the gas constant in dry air aid water vapor respectively. And $p_w$ and $p^s_w$ are the water vapor pressure and saturated vapor pressure
\begin{align}
    Y &= \frac{\tilde{\rho}_w}{\tilde{\rho}} = \frac{\tilde{p}_w}{R_w\tilde{\rho}\tilde{T}} = \frac{\tilde{p}_w}{R_w\tilde{T}(\tilde{\rho}_w+\tilde{\rho})} = \frac{\tilde{p}_w}{R_w\big(\frac{\tilde{p}_w}{R_w}+\frac{\tilde{p}_a}{R_a}\big)}= \frac{\tilde{p}_w}{p_w+\big(\tilde{p}-\tilde{p}_w\big)\frac{R_a}{R_w}} \nonumber \\
    &= \frac{R_a}{R_w}\frac{\tilde{p}_w}{\big(\tilde{p}-\tilde{p}_w\big(1-\frac{R_a}{R_w}\big)\big)} = \frac{R_a}{R_w}\frac{\phi\cdot\tilde{p}^s_w}{\big(\tilde{p}-\phi\cdot\tilde{p}^s_w\big(1-\frac{R_a}{R_w}\big)\big)}
\end{align}

with the relative humidity, $\phi$
\begin{equation*}
    \phi = \frac{\tilde{p}_w}{\tilde{p}^s_w}
\end{equation*}
Obtaining saturation function with relative humidity being unity. i.e. $\phi=1$ and after the pressure split for the small Mach number correction $p=p_0=1$.(not necessary in this question). I obtain the saturation function, $Y^s$ as follows
\begin{equation}\label{1}
    Y^s(T) = = \frac{R_a}{R_w}\frac{ p^s_w(T)}{\big(1-\big(1-\frac{R_a}{R_w}\big)\cdot p^s_w(T)\big)}
\end{equation}

So I expect the $Y^s$ to be 1 when the air is saturated with vapor. So I tried to model this with several saturated vapor pressure functions $p^s_w$. as follows

Antonie Function\begin{equation}
    p^s_w(T) = ptr\cdot exp\Big(17.28-\frac{4102.99}{T-35.719}\Big)
\end{equation}
where $ptr=611.657 Pa$ is the triple point of the water.

Emporical Vapor pressure equations, Swietoslawski and Smith (1938)\begin{align}
    log P_w^s(T) &= A - \frac{B}{T}+\frac{CX}{T}(10^{DX^2}-1)-E(10^{-FY^{5/4}})
\end{align}
where $X=T^2 - 2.937x10^5; A= 8.3074650; B=2.0051x10^3; C=1.3869x10^{-4}; D = 1.1965x10^{-11}; E=4.4x10^{-3}; F = 5.7148x10^{-3}; Y= 647.27 - T$

Saturation vapor pressure of water and Ice, J. Huang (2018)
\begin{align}
    P^s_w = \frac{exp(34.494-\frac{4294.99}{T+237.1})}{(T+105)^{1.57}} \text{ for $T>0$}
\end{align}

Following is the Matlab code which I used to calculate this $Y^s$ using the above saturation functions. But my problem is I did not receive the expected saturation curve here usually got negative values or a constant. The tests are done with the temperature data with $\tilde{T}\in[273.16,333.16]$ and reference temp $T_r=T(0,t)=296.5K$ as the temperature at the bottom of the tower.
What I expect here is that $Y^s$ grows with the temperature $T$ and after certain temperature, $Y^s$ should be constant(ideally 1). But I don't get the required results. I don't know actually my code is wrong or I'm using inappropriate function.

```<code>clc

Rw = 461; Ra = 287.5;  % gas constant of water wapor and dry air
Tr = 296.5;

Pws = Ra/Rw;
Ptr = 611.66; % triple point pressure of water

%-------------------------------------------
% Testing Saturation function

% Temperature range

 T2 = sort(randi([273 400],1,30));% generate a random temperature distribution
 T1 = T2/Tr;

n = size(T1,2);
x = linspace(0,1,n);

Ys1 = zeros(1,n);
Ys2 = zeros(1,n);
Ys3 = zeros(1,n);

% saturation pressure calculation
    for j=1:n
      Ys1(j,n) = antonie(T1(j),Tr,Pws);
      Ys2(j,n) = ss(T1(j),Tr,Pws,Ptr);

      Ys3(j,n) = ijh(T1(j),Tr,Ptr,Pws);
    end

Ys

% Construct the desired plot
subplot(2,2,1)
plot(T1*Tr,Ys1);
xlabel('Temperature, T');
ylabel('Sat. Func. Value, Y^s');
title('Antonie Function');

subplot(2,2,2)
plot(T1*Tr,Ys2);
label('Temperature, T');
ylabel('Sat. Func. Value, Y^s');
title('Swietoslawski and Smith');

subplot(2,2,3)
plot(T1*Tr,Ys3);
xlabel('Temperature, T');
ylabel('Sat. Func. Value, Y^s');
title('Ice, J. Huang');
</code>```

% saturation pressure functions

```<code>function g = antonie(T,Tr,pws)
    T= T*Tr;
    ptr = 611.657;
    g = ptr*exp(17.2799-(4102.99/(T-35.719)));
    Ys = pws*g/(1-(1-pws)*g);
end

% According to the paper: Derivation and use of Antonie function for the
% temp. range 0-374(C°)
function Ys = ss(T,Tr,pws,ptr)
    T = T*Tr;
    A = 8.3074650; B = 2005.1; C= 1.3869*10^(-4); D= 1.1965*10^(-11);
    E = 4.4*10^(-3); F = 5.7148*10^(-3);
    Y = 374.11-T; X = T^2-2.937*10^5;
    g = ptr*exp(A-(B/T)+(C*X/T)*(10^(D*(X^2))-1)-E*(10^(-F*(Y^(5/4)))));
    Ys = pws*g/(1-(1-pws)*g);
end

% Refer : https://journals.ametsoc.org/view/journals/apme/57/6/jamc-d-17-0334.1.xml?tab_body=pdf
function ys = ijh(T,Tr,ptr,pws)
    T = T*Tr;
    t = T+237.1;
    g = ptr*exp(34.494-(4924.99/t))/((T+105)^1.57);
    ys = pws*g/(1+(1-pws)*g);
end
</code>```


问:
Once a method is implemented in a computer program, then does learning the theory become less useful use of time?
Once a method is implemented in a computer program, then does learning the theory become less useful use of time?

This has confused me. I used to view that all things should be studied from first principles, but lately I've thought that this might well be uneconomical.

For example, consider FreeFEM or something. Once it's usable there, then how much theory does one really need to study as an end-user?

But is there some way to gauge what one should do?

My idea was to start with the software and consult theory when needed.

Rather than learn theory before learning applications (the common way in academia).

If the same applies as with APIs in general, then the end-user need not know full implementational details.


问:
Parallel Plate Flow simulation using SIMPLER Alogortithm
I am new to the community and CFD itself and am working on my first CFD program in Scilab based on the SIMPLER algorithm in Scilab i.e. modeling a flow between parallel plates with entrance region. I want to get a Poiseuilli kind of flow profile but somehow I'm not able to reach the desired solution.

It's a 2D-staggered grid arrangement (Central Differencing) like the image below (image just for illustration), the black points being the points for pressure.

And here's my code for reference

https://drive.google.com/file/d/1ludwYAwrkn43DaDQg1lCzWKmduQiIefM/view?usp=sharing

I get satisfactory values for the first couple of iterations but then it starts giving random values (some even negative) for velocity. I think the problem could be with the Pressure equation but even after analyzing for days now, I concluded that the code is perfect theoretically i.e. the discretization of momentum equations, continuity equation, Boundary conditions, and implementation of the SIMPLER algorithm is as per the theory and I have debugged minor errors but still can't point out why it's diverging.

Any help will be appreciated.


问:
Nimfa Lsnmf matrix factorization with sparse input: wrong/different result than dense input
I just started using Nimfa and working on a dummy example. I would like to do a matrix factorization (on a matrix for which we do not have all the entries) and I am assuming the following model (which is for unconstrained matrix factorization from Agarwal book. It calculate the error/cost/residual only on observed entries of the matrix. S is the set of observed ratings of customer i for item j):

I have the following example:

```<code>    r_dense is
    array([[2, 0, 2, 0, 0, 0],
           [2, 2, 2, 0, 0, 0],
           [2, 2, 2, 0, 0, 0],
           [2, 2, 2, 2, 2, 2],
           [1, 1, 1, 2, 2, 2]])

    lsnmf = nimfa.Lsnmf(r_dense, seed='random_vcol', rank=2, max_iter=100)
    lsnmf_fit = lsnmf()
    print('Rss: %5.4f' % lsnmf_fit.fit.rss())
    print('Evar: %5.4f' % lsnmf_fit.fit.evar())
    print('K-L divergence: %5.4f' % lsnmf_fit.distance(metric='kl'))
    print('Sparseness, W: %5.4f, H: %5.4f' % lsnmf_fit.fit.sparseness())
    U = lsnmf_fit.basis()
    V = lsnmf_fit.coef()
    rhat=U.dot(V)
    print(rhat)
</code>```

And the results are:

```<code>    nimfa on dense matrix
    Rss: 2.1661
    Evar: 0.9695
    K-L divergence: 1.5560
    Sparseness, W: 0.4460, H: 0.5728
    [[1.56967575 1.162301   1.56967575 0.         0.         0.        ]
     [2.14928335 1.59677013 2.14928335 0.02665226 0.02665226 0.02665226]
     [2.14928335 1.59677013 2.14928335 0.02665226 0.02665226 0.02665226]
     [2.03517875 1.90490806 2.03517875 2.00634408 2.00634408 2.00634408]
     [0.96082184 1.10669734 0.96082184 1.99283382 1.99283382 1.99283382]]
</code>```

As you can see, the RSS = 2.17. If I plug in the basis into the model I mentioned, the value is 0.81. Hence, Lsnmf function is not using the model I mentioned and for the error uses all the entries. I am a bit confused. Is there a way to force Nimfa function to only use known values for error calculation or I need to write the factorization myself from scratch.

I will be very grateful if you can help me.