问:
Stochastic differential equation system (SDE) : overflow encountered in double-scalars
I'm trying to integrate the following SDE system from Dekker et al. [1]

$$\begin{cases}
\frac{dx}{dt}=a_1x^3+a_2x+\phi+\zeta_x\\
\frac{dy}{dt}=b_1z+b_2(\kappa(x)-(y^2+z^2))y+\zeta_y\\
\frac{dz}{dt}=c_1y+b_2(\kappa(x)-(y^2+z^2))z+\zeta_z
\end{cases}
$$

To do that I'm using the Euler-Maruyama method which is basically a forward-Euler scheme plus a Gaussian noise term where the mean and the standard deviation parameters are picked up from the Table 2 in Dekker et al. article.

We have

Noise mean : 0

Noise variance : 0.1

Integration time : 500

Time step : 0.5

The parameters $a_1, a_2,...$ as well as the coupling term $\kappa(x)$, the integration domain, the step time and the initial conditions are all picked from Table 2.

However, running that simulations I get an RuntimeWarning: overflow encountered in double-scalars. Seems like step time would be too big but that's what is recommended in the article.

I'm wondering why I get this error and how I could fix it?

As you can see, after a few iterations, the decimal of x, y and z become exponentially big. It doesn't come from the stochastic term as the print (i.e the array [x, y, z]) in the below capture is the print of what comes out from the forward-euler alone.

Code

# 1D EDO for Fold bifurcation paramaters
a1 = -1
a2 = 1

# 2D EDO for Hopf bifurcation parameters
b1 = b2 = 1
c1 = -1
c2 = 1

# coupling parameters
gamma1 = -0.1
gamma2 = 0.12

# initial conditions
x0 = -0.5
y0 = 1.0
z0 = -1.0

r0 = 5.0

# time range
t_init = 0
t_fin = 500
time_step = 0.01

# gaussian noise parameters
mean = 0.0
variance = 0.1

# parameters for stochastic plot
a1_stoch = -1
a2_stoch = 1
b1_stoch = 0.1
b2_stoch = 1
c1_stoch = -0.5
c2_stoch = 1
gamma1_stoch = -0.2
gamma2_stoch = 0.3
time_step_stoch = 0.5
```
<pre class="lang-py prettyprint-override"><code>import numpy as np
from consts import (a1_stoch, a2_stoch, b1_stoch, b2_stoch, c1_stoch, c2_stoch, gamma1_stoch, gamma2_stoch)

# linear coupling parameter
# proposed by Dekker et al. article
def gamma(x):
  return gamma1_stoch + (gamma2_stoch * x)

def fold_hopf_stoch(v, phi):
  print(v, phi)
  return np.array([
    a1_stoch * (v[0] ** 3) + a2_stoch * v[0] + phi,
    b1_stoch * v[2] + b2_stoch * (gamma(v[0]) - (v[1]**2 + v[2]**2))*v[1],
    c1_stoch * v[1] + c2_stoch * (gamma(v[0]) - (v[1]**2 + v[2]**2))*v[2]
  ])

</code>```

import numpy as np
from consts import (mean, variance)

"""
forward euler method with a stochastical term
x_{i+1} = x_i * dt + zeta + sqrt(dt)
(forward-euler + gaussian noise * sqrt(dt))
"""
def forward_euler_maruyama(edo, v, dt, *args):
  zeta = np.random.normal(loc=mean, scale=np.sqrt(variance))
  print(zeta, edo(v, *args) * dt)
  return edo(v, *args) * dt + zeta * np.sqrt(dt)

```
<pre class="lang-py prettyprint-override"><code>import numpy as np
import matplotlib.pyplot as plt

from tqdm import tqdm

from .edo import (fold_hopf)
from .edo_stoch import (fold_hopf_stoch)
from .rk4 import rk4
from .euler import forward_euler_maruyama

from consts import (x0, y0, z0, t_init, t_fin, time_step, time_step_stoch)

plt.style.use("seaborn-whitegrid")
np.set_printoptions(precision=3, suppress=True)

class time_series():
  def __init__(self):
    # initial conditions
    self.initial_conditions = np.array([[x0, y0, z0]])
    # forcing parameter
    self.phi = -2
    # time
    self.t0 = t_init
    self.tN = t_fin
    # stochastic number of iterations
    self.niter = 10
    # legend
    self.legends = ["<span class="math-container">$x$</span> (leading)", "<span class="math-container">$y$</span> (following)", "<span class="math-container">$z$</span> (following)"]

  def solve(self, solver, edo, dt, nt):
    # [x0, y0, z0]
    #v = self.initial_conditions[0]

    # vector mesh -- will receive the result
    v_mesh = np.ones((nt, 3)) #
    # set inital conditions
    v_mesh[0] = self.initial_conditions[0]
    print(edo, dt)

    for t in tqdm(range(0, nt - 1)):
      v_mesh[t + 1] = v_mesh[t] + solver(edo, v_mesh[t], dt, self.phi)

      # increase forcing parameter
      self.phi += 0.001

    return v_mesh

  def basic(self):
    dt = time_step
    nt = int((self.tN - self.t0) / dt)
    time_mesh_basic = np.linspace(start=self.t0, stop=self.tN, num=nt)

    # [ [x0, y0, z0], [x1, y1, z1], ..., [xN, yN, zN] ]
    results = self.solve(rk4, fold_hopf, dt, nt)
    return time_mesh_basic, results

  def stochastic(self):
    dt = time_step_stoch
    nt = int((self.tN - self.t0) / dt)
    time_mesh_stoch = np.linspace(start=self.t0, stop=self.tN, num=nt)

    stochastic_results = np.ones((self.niter, nt, 3))
    # compute a lot of simulations
    for i in tqdm(range(0, self.niter)):
      results = self.solve(forward_euler_maruyama, fold_hopf_stoch, dt, nt)
      stochastic_results[i] = results
      return
    return time_mesh_stoch, stochastic_results

  def plot(self):
    #time_mesh_basic, basic_results = self.basic()
    time_mesh_stoch, stochastic_results = self.stochastic()
    return

    # 3 lines, 1 column
    fig, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(15, 7))
    fig.suptitle("Série temporelle")

    ax1.plot(time_mesh_basic, basic_results[:,0])
    ax1.plot(time_mesh_basic, basic_results[:,1])
    ax1.plot(time_mesh_basic, basic_results[:,2])

    ax1.set_xlabel("<span class="math-container">$t$</span>")
    ax1.set_ylabel("<span class="math-container">$x$</span>, <span class="math-container">$y$</span>, <span class="math-container">$z$</span>")
    ax1.set_xlim(0,500)
    ax1.set_ylim(-1.5,1.5)
    ax1.legend(self.legends, loc="upper right")
    ax1.set_title("Basique")

    print(stochastic_results)

    ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,0])
    ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,1])
    ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,2])

    ax2.set_xlabel("<span class="math-container">$t$</span>")
    ax2.set_ylabel("<span class="math-container">$x$</span>, <span class="math-container">$y$</span>, <span class="math-container">$z$</span>")
    ax2.set_xlim(0,500)
    ax2.set_ylim(-1.5,1.5)
    ax2.legend(self.legends, loc="center left", bbox_to_anchor=(1,0.5))
    ax2.set_title("Stochastique")

    plt.tight_layout()
    plt.show()

</code>```

References:

Dekker, M. M.; von Der Heydt, A. S.; Dijkstra, H. A. Cascading transitions in the climate system. Earth Syst. Dyn. 2018, 9 (4), 1243–1260. DOI: 10.5194/esd-9-1243-2018.


问:
Initial condition sensitivity issue in MINLP with Gekko
I'm having a problem with sensitivity to initial conditions in an MINLP problem with Gekko. My test case is finding $a$ and $b$ such that $\pi^a/b = \pi^2/6$. So the solution is $a=2$ and $b=6$. I set the initial values of $a$ and $b$ to $a0$ and $b0$ in the code. But i get different solutions depending on what these initial values are. Here are some examples:

```<code>Result 1: a0 = 5, b0 = 5, a = [6.0], b = [585.0], objective = 0.0015337359503
Result 2: a0 = 5, b0 = 6, a = [5.0], b = [185.0], objective = 0.0092263914506
Result 3: a0 = 2, b0 = 6, a = None, b = None, objective = 0.0
Result 4: a0 = 0, b0 = 0, a = [0.0], b = [0.0], objective = 1e+99
Result 5: a0 = 2, b0 = 5, a = [2.0], b = [6.0], objective = 0.0
</code>```

Result 5 gives the correct solution, but there seems to be many different issues here. Why does result 1 and 2 miss the correct solution? Why does the solver give up when given the correct values? (Result 3). What is going on in result 4? How do i set up Gekko to find the correct solution here, independent of initial values, and keeping lb = 0?

```<code>import math
from gekko import GEKKO

m = GEKKO(remote = False)
m.options.SOLVER = 1 # MINLP

pi = m.Const(math.pi)
c = m.Const(math.pi**2 / 6)

a0 = 5
b0 = 5

a = m.Var(value = a0, lb = 0, ub = 5000, integer = True)
b = m.Var(value = b0, lb = 0, ub = 5000, integer = True)

m.Obj(abs((pi**a)/(b)-c))

try:
    m.solve(disp = False)
except:
    a.value = None
    b.value = None

print(f'Results: a0 = {a0}, b0 = {b0}, a = {a.value}, \
b = {b.value}, objective = {m.options.objfcnval}')
</code>```

ContentLicense: CC BY-SA 4.0


问:
Solving coupled PDEs with self-consistency condition
I am figuring out how to attack a problem (the Usadel equations of superconductivity) in which I need to solve a set of nonlinear PDEs for the fields $\{G_i (r)\}$
$$ U(G_i(r), \nabla G_i(r), \Delta(r) ) = 0\,, $$
together with a (nonlinear, algebraic) self-consistency condition
$$V({G_i (r)}, \Delta(r) ) = 0\,.$$

The usual way to solve this problem in the literature is to discretize the PDEs using finite differences and with a segregated approach: solving the PDEs given $\Delta(r)$ for the fields $G_i(r)$, then computing the self-consistent field $\Delta$ solving the non-linear equation, then use the updated $\Delta(r)$ for a new iteration of PDEs and so on until convergence.

I was wondering about the feasibility of using a finite element approach and maybe even solving the system of equations as a fully-coupled problem. In particular, I would like to use FEniCS. Is it doable? Do you know any example of a similar computational problem I can study?

I guess coding the solver for the PDEs would be fairly simple and it would probably be easy to find some examples. What worries me is the self-consistency condition as it's a non-linear algebraic equation and I have not a clear idea on how to write it down in a finite-element environment.

When I tried something similar (for a different problem) in Comsol I found some convergence issues as the lack of gradient terms would cause strong spatial variations of the fields.

Edit:

The specific form of Usadel equations depends on the system studied and approximation used, however an example of the PDEs (taken from the supplementary material of https://journals.aps.org/prb/supplemental/10.1103/PhysRevB.103.L020504) is the following:
$$
\begin{cases}
(D \nabla^2 - 2 \omega_n) f^s_n(r) = - 2 \pi \Delta(r) + 2 i \vec{h}(r)\cdot\vec{f}^t_n(r) \\
(D \nabla^2 - 2 \omega_n) \vec{f}^t_n(r) =  2 i \vec{h}(r) f^s_n(r) \\
\end{cases}
$$

while the self-consistency equation may look like
$$
-\Delta(r)\ln(\tau) = 2 \tau \sum_{n=0}^{n_{cutoff}} \Big(\frac{\Delta(r)}{2 \omega_n}-f^s_n(r)\Big)
$$
In this case, $f^{s,t}_n(r)$ are the unknown fields in the PDEs while $\Delta(r)$ is the self-consistent field (unknown as well), while $\tau$ and $\vec{h}$ are the external parameters. In the numerical implementation, $n$ are integer numbers comprised between 0 and a cutoff which can be placed around a hundred. The domain is a simple rectangle.


问:
Finite element method for an equation requiring switch between spectral and temporal domain
Some equations (such as the non-linear schrödinger equation for pulse propagation) are more easily solved in the spectral form, but still need a representation in the temporal domain to calculate physical parameters such as the energy of the electric field. This results usually in the following steps (with $E$ the electric field in the temporal domain, and $\hat{E}$ the electric field in the spectral domain):
$$E_0\underbrace{\rightarrow}_{\text{FFT}}\hat{E}_0\underbrace{\rightarrow}_{\text{Propagation}}\hat{E}_1\underbrace{\rightarrow}_{\text{IFFT}}E_1$$
which requires two different "grids". Is such a calculation possible using FEM, and if yes, how (after requiring two different grids)?

Background to this question is that I have additional equations coupled to that problem which preferably can be solved using FEM and do not require the transformations from the temporal domain to the spectral domain and back (after they only depend on $E_x$, not on $\hat{E}_x$). For this coupling I'd like to avoid having to have two different grids (one for the non-linear schrödinger equation, one for the FEM-covered equation) which I have to interpolate data in between.


问:
Convergence of the Roothaan-Hall equations
Suppose that we are given a time-independent quantum mechanical system whose wavefunction depends on three space coordinates. Let $F$ be the Fock operator of the system. Suppose also that we have a set of functions $f_k : \mathbb{R}^3 \to \mathbb{C}$, $k \in \mathbb{N}$, so that $$\mathrm{clos} \;\mathrm{span} \{f_k : k \in \mathbb{N}\} = L^2(\mathbb{R}^3) .$$ Let $n \in \mathbb{Z}$ and $B_n = \{ f_k : k = 0,\ldots,n\}$. Let $$F_n C_n  = E_n S_n C_n$$ be the Roothaan-Hall equations of the system in basis set $B_n$. What are the necessary and sufficient conditions for the solutions of the Roothaan-Hall equations to converge to the exact solutions of the Hartree-Fock equations of the system when $n \to \infty$? Does $\{ f_k : k \in \mathbb{N}\}$ have to be a Riesz basis of $L^2(\mathbb{R}^3)$? We may consider only the ground state of the system.


问:
Solving two field Schrodinger-Poisson system numerically
I want to solve the system of Schrodinger-Poisson equations numerically:

\begin{align}
    \chi_1''(r) + \frac{2}{r}\chi_1'(r)&=2U(r)\chi_1(r) \\
    \chi_2''(r) + \frac{2}{r}\chi_2'(r)&=2\left(\frac{m_2}{m_1}\right)^2U(r)\chi_2(r) \\
     U''(r)+ \frac{2}{r}U'(r) &=   \chi_1^2(r) +  \left(\frac{m_2}{m_1}\right)^2  \chi_2^2(r) \label{eq:sys4}
\end{align}
where $\chi_1(r)$ and $\chi_2(r)$ are scalar fields, U(r) is the gravitaional potential (shifted by some constant). This problem is a BVP, where $\chi_1(r)$ and $\chi_2(r)$ approach zero as r tends to infinity, and U(r) approaches some constant, but I can rewrite the problem as an IVP where the initial conditions are (for $\delta \ll 1$):

\begin{align}
    \chi_1(\delta) &=1 + \frac{U_0}{3}\delta^2\\
    \chi_1'(\delta) &=2\frac{U_0}{3}\delta + \frac{1}{15} \left( 1 + 2U_0^2 + \left(\frac{m_2}{m_1}\right)^2 B_0^2 \right) \delta^3\\
    \chi_2(\delta) &= B_0 \left(1 + \frac{1}{3}(U_0 + \Delta \gamma) \left(\frac{m_2}{m_1}\right)^2 \delta^2\right)\\
    \chi_2'(\delta)&= \frac{B_0}{3}\left(\frac{m_2}{m_1}\right)^2 \left( 2(U_0 +\Delta \gamma ) \delta + \frac{1}{5}\left(1 + \left(\frac{m_2}{m_1}\right)^2(2(U_0 + \Delta \gamma)^2 + B_0^2)\right)\delta^3 \right) \\
    U(\delta) &= U_0 + \frac{1}{6}\left(1 + \left(\frac{m_2}{m_1}\right)^2 B_0^2 \right) \delta^2 \\
    U'(\delta) &= \frac{1}{3}\left(1 + \left(\frac{m_2}{m_1}\right)^2 B_0^2 \right)\delta + \frac{2}{15}\left(U_0 + \left(\frac{m_2}{m_1}\right)^4(U_0 + \Delta \gamma)B_0^2\right)\delta^3 \label{eq:sys5}
\end{align}

I want to find the solution without any nodes (the ground state solution). When calculating, I can fix $B_0$ and $\frac{m_2}{m_1}$ to be some constant, and now I have to vary two parameters, $\Delta \gamma$ and $U_0$ in order to obtain the nodeless solution.

Can someone suggest me an algorithm that could solve this system of equations and obtain the values of $\Delta \gamma$ and $U_0$ for which the solution is nodeless?

The version with only one field $\chi$ can be solved using the shooting method where the parameter that is varied is $U_0$ (this is the only parameter in that case). In the case of the two fields, the shooting method doesn't look feasible.


问:
Errors in Integral Estimate of Gaussian using Trapezoidal Rule
I'm trying to estimate the percentage error in computing the integral of a Gaussian via composite trapezoidal rule versus via an exact formula. To do this I've generated a gaussian with mean 0, standard deviation $\sigma = 1$ and amplitude $A = 1$ in Python (plotted below), on a grid of $N$ points in the range $[-2\sigma,2\sigma]$, and used scipy.integrate.trapz to compute the approximate integral $A$ for a range of different $N$.

For the exact integral I used the formula $I = \int_{-x}^{x}e^{\frac{-t^2}{2\sigma^2}} dt = \sqrt{2\pi}\sigma\text{erf}(\frac{x}{\sqrt{2}\sigma})$ to compute the exact integral between $\pm$ two standard deviations of the mean as $I = \int_{-2\sigma}^{2\sigma}e^{\frac{-t^2}{2\sigma^2}} dt = \sqrt{2\pi}\sigma\text{erf}(\sqrt{2})$, or

$I = \sqrt{2\pi}\text{erf}(\sqrt{2})$

in this $\sigma = 1$ case. The % error is then $E = \frac{|I - A|}{I}\cdot 100$. Below is a plot of $\frac{1}{E}$ versus $N$

and for reference a plot of $E$ verus $N$.

Now I would expect $E \propto \frac{1}{N^2}$ but this is not seen here. Would anyone be able to help me figure out why?

I've included my code below. Thank you in advance!

```<code>from scipy.integrate import trapz
from scipy.special import erf
import matplotlib.pyplot as plt
import numpy as np

errs = np.array([]) # Errors in integral estimates
Ns = np.array([]) # Sample sizes tried

def get_percent_err(N, errs, Ns):
    A = 1 # Amplitude of Gaussian
    std = 1 # Standard deviation of Gaussian
    r = 2*std # Extent of domain
    dt = 2*r/N # Step
    t = np.linspace(-r,r,N) # Domain
    gauss = A*np.exp(-((t)**2)/(2*(std**2))) # Gaussian function

    actual_area = np.sqrt(2*np.pi)*A*std*erf(np.sqrt(2)) # Exact gaussian integral from -2*std to +2*std
    est_area = trapz(gauss, dx=dt) # Estimated integral via trapezium rule
    percent_err = (abs(actual_area - est_area)/actual_area)*100 # Percentage error

    errs = np.append(errs, percent_err)
    Ns = np.append(Ns, N)

    return errs, Ns

# Sweep through a range of sample sizes and obtain errors for each
min_samples = 100
max_samples = 100100
sample_step = 100
samples = np.arange(min_samples,max_samples,sample_step)

for N in samples:
    errs, Ns = get_percent_err(N, errs, Ns)

# Plot reciprocal of error in integral against sample size
plt.figure()
plt.plot(1/Ns,errs)
plt.xlabel("1 / Samples (1 / N)")
plt.ylabel("% Error in Integral (E)")
plt.title("Plot of % Error (E) in Integral versus 1 / Samples (1 / N)")
</code>```

UPDATE:

For reasons that are unknown to me (though I believe it has something to do with numerical precision), this code, adapted from this link, gives an expected result.

Code:

```<code>import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
from scipy.integrate import trapz as trapz_scipy
from pandas import DataFrame

def trapz(f,a,b,N=50):
    x = np.linspace(a,b,N+1) # N+1 points make N subintervals
    y = f(x)
    dx = (b - a)/N
    T = trapz_scipy(y, dx=dx)
    return T

def get_err(N, errs):
    A = 1
    std = 1
    f = lambda x : A*np.exp(-x**2/(2*std**2))
    a = -2*std; b = 2*std;

    I = np.sqrt(2*np.pi)*A*std*erf(np.sqrt(2))
    T = trapz(f,a,b,N)
    return np.append(errs, 100*np.abs(I - T)/T)

min_samples = 10
max_samples = 1000
sample_step = 10
sample_sizes = np.arange(min_samples,max_samples+sample_step,sample_step)

errs = np.array([], dtype=np.float64)

for s in sample_sizes:
    errs = get_err(s, errs)

err_vs_sample_size = DataFrame({'% Error' : errs, 'Samples' : sample_sizes})

plt.figure()
plt.plot(1/sample_sizes**2,errs)
plt.xlabel("1 / Samples (1 / N)")
plt.ylabel("% Error in Integral (E)")
plt.title("Plot of % Error (E) in Integral versus 1 / Samples (1 / N)")
</code>```

Error plot:


问:
Approximation in the derivation of the Arc Length method
I am studying the proof of the Arc Length method in section 2.2 of this thesis. In equation (2.2) the author introduces the supplementary conditions

$$
(\Delta {\bf u} + \delta {\bf u})^T \cdot (\Delta {\bf u} + \delta {\bf u}) + \psi^2(\Delta \lambda + \delta \lambda)^2 \cdot ({\bf q}^T \cdot {\bf q}) = \Delta L^2
$$
thereby adding the final equation to the system

$$
[K_T] \cdot \delta {\bf u} - \delta \lambda \, {\bf q} = -R({\bf u}^\prime , \lambda^\prime)
$$

The author gives the resulting system as

$$
\begin{bmatrix}
[K_T] & -{\bf q} \\
2 \Delta {\bf u}^T & 2 \psi^2 \Delta \, ({\bf q}^T \cdot {\bf q})
\end{bmatrix}
\cdot
\begin{bmatrix}
\delta {\bf u} \\
\delta \lambda
\end{bmatrix}
=
- \begin{bmatrix}
{\bf R} \\
\delta \lambda
\end{bmatrix}
$$

Where

$$
A = -[\Delta {\bf u}^T \cdot \Delta {\bf u} + \psi^2 \Delta \lambda^2 \, ({\bf q}^T \cdot {\bf q}) -\Delta u^2]
$$

It seems to me that the terms $\delta {\bf u}^T \cdot \delta {\bf u}$  and $\delta \lambda^2$ have been dropped. Is this because the terms are because they are small and nonlinear? And if so, is this why the Arc Length method is iterative i.e. in order to iteratively eliminate the errors caused by dropping these terms?


问:
Semi-analytical/empirical modelling of wall boundary conditions in advection-diffusion-reaction equation with distributed source
Let's suppose I need to numerically solve a 3D steady-state transport equation of the form
$$
\nabla \cdot (\mathbf{u} c) = \nabla \cdot (D \nabla c) - \lambda c + S
$$
where $c$ is the transported quantity, $\mathbf{u}$ is a velocity field, $D$ is the diffusivity (optionally varying with position), $\lambda$ is a uniform linear decay rate and $S$ a spatially varying distributed source.

Let's also assume that the equation is to be solved in a pre-existent Finite Volume framework, coupled with the "simultaneous" solution of velocity and other quantities, from which $S$ is obtained.

The equation is subject, at specified wall boundaries, to a homogeneous mixed-type BC:
$$
- D \nabla c \cdot \mathbf{n} + \gamma c = 0
$$
where $\mathbf{n}$ is the outward pointing wall-normal direction and $\gamma$ represent a "deposition velocity". The mixed-type BC is used to stay a little general, but indeed $\gamma$ should ideally tend to $+\infty$ to approximate a homogeneous Dirichlet BC.

Preliminary analysis has shown that the combination of the values of physical parameters for the problem at hand (e.g. low diffusivity, slow decay) together with large deposition velocity and a non-zero source at walls produces dramatically steep gradients in very thin wall layers. If my interest is in both grasping the 3D distribution of $c$ and the deposition fluxes at walls, this would require mesh refinement close to walls up to unaffordable levels.

To give an idea, it can be easily verified that the analytical solution of the corresponding problem in a fully-developed laminar flow between parallel plates, with uniform source $S$ and uniform properties, yields a solution (in the limit of infinite $\gamma$) which approaches a uniform value of $S/\lambda$ everywhere except for close to the walls, where it rapidly falls to zero according to an exponential law with a characteristic length of $\sqrt{D/\lambda}$.

My question then is: is there some semi-analytical and/or semi-empirical method to treat the wall boundary conditions in such a situation? Specifically, I'm looking to some way to model the wall gradient instead of resolving it in the wall layer.

I'm thinking of something similar to what is achieved with the Wall Function approach in standard turbulence modelling.

The setup with distributed source term is not very common and I'm struggling a little bit to find something suitable. Maybe I'm just using the wrong keywords.

Anyway, any comment/advice/suggestion/reference would be greatly appreaciated.


问:
discretization of advection diffusion with variable coefficients
I am looking for help to find a somewhat stable FD numerical scheme for the advection diffusion equation posed on a curve $(x(r),y(r))$. The equation becomes

$$u_t=\alpha(r) u_r +\beta(r) u_{rr}+f(r,t),$$

where $\alpha(r)= \phi(r)\phi'(r)+\phi(r) $, $\beta(r)=\phi(r)\phi(r)$, and $\phi(r)=\frac{1}{\sqrt{x'(r)^2+y'(r)^2}}$.

I'm specifically trying to use periodic boundary conditions, but any will do. I have tried to use a centered difference discretization combined with a forward Euler, but it is diffusing very quickly. Any advice would be appreciated!


问:
Software for solving large systems of linear equations over gf(2)
What available solvers are there for linear equation solver over GF(2) (Boolean), capable of dealing with large sparse systems (in the 10k - 100k variables range)?</p>


问:
Newton-Raphson for estimating Weibull distribution does not converge
I've been trying to estimate the two-parameter (a,b) Weibull distribution (loc. param. = 0).

$$f(t;a,b)=\frac{b}{a}\left(\frac{t}{a}\right)^{b-1}\exp(- \left(\frac{t}{a}\right)^b) $$

To find the estimates I am using maximum likelihood and NR method.

```<code>t = np.array(list(range(1, 1000))) #data input
def gradient(a,b):
    grad_a = -10*b/a + b/a*np.sum((t/a)**b)
    grad_b = 10/b - 10*(np.log(a)) + np.sum(np.log(t)) - np.sum(((t/a)**b)*np.log(t/a))
    grad_matrix = np.array([grad_a, grad_b])
    return grad_matrix

def hessian(a,b):
    hess_a = 10*b/a**2 + (b*(b+1)/a**2)*np.sum((t/a)**b)
    hess_b = 10/(b**2) + np.sum((t/a)**b* (np.log(t/a))**2)
    hessians = np.array([hess_a, hess_b])
    return hessians

iters = 0
a0, b0 = 2,6 #initial condition

while iters < 4:
    if hessian(a0,b0).any() == 0.0:
        print('Divide by zero error!')
    else:
        a = a0 - gradient(a0,b0)[0]/hessian(a0,b0)[0]
        b = b0 - gradient(a0,b0)[1]/hessian(a0,b0)[1]
        print('Iteration-%d, a = %0.6f, b= %0.6f, e1 = %0.6f, e2 = %0.6f' % (iters, a,b,a-a0,b-b0))
    if gradient(a,b)[0]==0 and  gradient(a,b)[1]==0:
        break
    else:
        a0,b0 = a,b
        iters = iters +1
print(a,b)
print(iters)
</code>```

Output:

Iteration-0, a = 1.714391, b= 6.656262, e1 = -0.285609, e2 = 0.656262

Iteration-1, a = 1.490482, b= 7.249362, e1 = -0.223909, e2 = 0.593099

Iteration-2, a = 1.309805, b= 7.794966, e1 = -0.180677, e2 = 0.545604

Iteration-3, a = 1.160878, b= 8.303167, e1 = -0.148927, e2 = 0.508201

a goes to zero, b to infinity, neither of which is a correct estimation. It feels like I am doing something wrong in defining gradients and hessians, though these seem to be correct first and second derivatives of the log-likelihood function.


问:
Boundary conditions for an FEM approximation of the Laplace operator
Using FEM, I want to approximate the Laplacian

$$u = \nabla \cdot \nabla h \, ,$$

where $h(x,y)$ is an FEM approximated scalar field on the same mesh, i.e. piecewise differentiable.

I am using MOOSE to solve the following (hopefully correct!) weak formulation of the above equation:

$$\int_\Omega uw_i \, \text{d}\Omega = -\int_\Omega \nabla h \cdot  \nabla w_i \, \text{d}\Omega + \int_\Gamma \nabla h \cdot \vec{n}w_i \, \text{d}\Gamma \,,$$

where $w_i$ is the i-th test function and $\vec{n}$ is the normal vector on the boundary $\Gamma$ of the computational domain $\Omega$. I am not using any boundary conditions on $u$.

The FEM solution that I obtain for this weak formulation is inaccurate at the boundaries (see below). I have tested this by prescribing the field $h=\text{sin}(x)$ on a 2D rectangular domain from $-2\pi \le x \le 2\pi$
and solving with quadratic Lagrange shape functions for both $h$ and $u$. (it is also inaccurate using linear shape functions, but it is easier to spot with quadratic shape functions).

Is it possible to obtain an accurate $u$ at the boundaries with this approach? Am I missing a boundary condition/contribution?

Background:

I am dealing with a thin film approximation of the Navier-Stokes equations in 2D FEM, where the film height is modeled by a dependent variable $h(x,y)$. In order to include surface tension as a contribution to the pressure gradient, I want to find the curvature of the film surface $u(x,y)$ that I approximate by the second derivative

$u \approx \nabla \cdot \nabla h$ .

My idea is to introduce curvature $u$ as a variable, solve the above equation by FEM, and couple $u$ into the momentum equation.


问:
How to find a lot of (if not all) local minima / critical points of a function?
Briefly stated, I would like to find "all" local minima / critical points of a function.

This function comes from the discretization of a continuous problem with infinitely many degrees of freedom (dofs), but after discretization, there remains actually rather "few" dofs. Let's say, n_dof ~ 1000 - 10k.

This function (let's denote it $S$) is a discretization of the action of a physical system. Glossing over some details, saying that a set of parameters $x$ "follow the rules of physics" is equivalent to saying that "$x$ is a critical point of $S$", i.e. " $\nabla S (x) = 0$". Since local minima are critical points, which might be easier to find numerically, I would also be happy with answers specifically addressing minima.

I would qualify my current approach as "Monte-Carlo like". It can be summed up as :

I first designed an efficient way to find a nearby critical point from a given estimate, using Newton-like methods, dealing with conditioning, etc.

I then randomly select starting points for the algorithm, and apply the above "local" method to find a critical point

I then check the newly found critical point against my database of already found points. If it is unknown, I add it to the database, else I start over.

Although this method is very easy to implement, it offers little theoretical guarantees. I don't know how many critical points "remain to be found" (say, within a compact set of parameters). In practice, the local method very often converges towards points that were already found (this morning : 100k initial guesses => 50 unique critical points). I don't know when to stop the algorithm, or what would be the expectation of finding a new critical point in less than a given number of iterations.

Do you have any suggestion on how to find more unique critical points / all critical points ?

Remark: I think there might be links to Morse theory which I know only little about.


问:
Numerical Range of a matrix in Python
In the mathematical field of linear algebra and convex analysis, the numerical range or field of values of a complex

$n\times n$ matrix $A$ is the set

$$W(A)=\left\{{\frac  {{\mathbf  {x}}^{*}A{\mathbf  {x}}}{{\mathbf  {x}}^{*}{\mathbf  {x}}}}\mid {\mathbf  {x}}\in {\mathbb  {C}}^{n},\ x\not =0\right\}$$

where $\mathbf{x}^*$ denotes the conjugate transpose of the vector $\mathbf {x}$ .

In engineering, numerical ranges are used as a rough estimate of eigenvalues of A. Recently, generalizations of the numerical range are used to study quantum computing.

A related concept is the numerical radius, which is the largest absolute value of the numbers in the numerical range, i.e.

$$r(A)=\sup\{|\lambda |:\lambda \in W(A)\}=\sup _{{\|x\|=1}}|\langle Ax,x\rangle |.$$

I want to implement the numerical range of a matrix in Python. I found the following code but contains a function command mlab.find which no longer exists. I also found how to replace it but the output is not the desired one.
Any help?

import numpy as np
import matplotlib.pyplot as plt

def find(condition):
    res,  =  np.nonzero(np.ravel(condition))
    return res

def numerical_range(A, resolution=0.01):
    A = np.asmatrix(A)
    th = np.arange(0, 2*np.pi + resolution, resolution)
    k = 0
    w = []
    for j in th:
        Ath = np.exp(-1j*j)*A
        Hth = (Ath + Ath.H)/2
        e,r = np.linalg.eigh(Hth)
        r = np.matrix(r)
        e = np.real(e)
        m = e.max()
        s = find(e == m)
        if np.size(s) == 1:
            w.append(np.matrix.item(r[:,s].H*A*r[:,s]))
        else:
            Kth = 1j*(Hth - Ath)
            pKp = r[:,s].H*Kth*r[:,s]
            ee,rr = np.linalg.eigh(pKp)
            rr = np.matrix(rr)
            ee = np.real(ee)
            mm = ee.min()
            sm = find(ee == mm)
            temp = rr[:,sm[0]].H*r[:,s].H*A*r[:,s]*rr[:,sm[0]]
            w.append(temp[0,0])
            k += 1
            mM = ee.max()
            sM = find(ee == mM)
            temp = rr[:,sM[0]].H*r[:,s].H*A*r[:,s]*rr[:,sM[0]]
            w.append(temp[0,0])
        k += 1
    plt.plot(w)
    return None

A = np.matrix([[.6,.5], [.6,.5]])
numerical_range(A, resolution=0.01)
plt.show()
```</span>
</code>```

Title: Numerical Range of a matrix in Python


问:
Find time step for Euler method in numerical solving of a system of non linear differential equations
I have a system of non linear differential equations in the form :
$$\frac{dy_i}{dt}=\sum_j a_{ij} y_i y_j $$.

I first tried to solve it with Python suing scipy.integrate.odeint but it is very slow.

I'd like to solve it with an Euler method : $$y_i(t+\delta t)=y_i(t)+\frac{dy_i}{dt}\delta t $$.

How should I chose the value of $\delta t$ to ensure a good convergence ?

What I tried is the following.

Let's define the error $\epsilon$ :
$$\epsilon=\sqrt{\sum_i(\delta t \times dy_i/dt)^2} $$

Starting from a low enough initial guess for $\delta t$, I looked wether the error increase or descrease.

If the error decreases : $\epsilon(t)>\epsilon(t+\delta t)$, I change $\delta t\rightarrow 2\delta t$

If the error increases : $\epsilon(t)<\epsilon(t+\delta t)$, I change $\delta t\rightarrow \delta t/2$

But the calculation is very slow and can even diverge.

Another possibility I tried :

If the error decreases : $\epsilon(t)>\epsilon(t+\delta t)$, I change $\delta t\rightarrow 2^n\delta t$ with $n$ so that $\epsilon(t+2^{n}\delta t)<\epsilon(t)<\epsilon(t+2^{n+1}\delta t)$

If the error increases : $\epsilon(t)<\epsilon(t+\delta t)$, I change $\delta t\rightarrow \delta t/2^n$ until $\epsilon(t+\delta t/2^{n})<\epsilon(t)<\epsilon(t+\delta t/2^{n-1})$

but this scheme is not very efficient too.


问:
Blown-up iterates in Gauss-Newton method
I am working on a non-linear least squares problem with standard form, in which I need to calibrate a parameter vector $\Theta$ to a set of inputs $\mathbf{x}$ and outputs $\mathbf{y}$:

$$\begin{align}
& \min_\Theta ||R(\mathbf{x},\mathbf{y};\Theta)||_2^2
\\[4pt]
&\quad \text{st: } R(\mathbf{x},\mathbf{y};\Theta):=\mathbf{y}-f(\mathbf{x};\Theta)
\end{align}$$

More specifically, I want to fit the 2-dimensional vector $\Theta=(\theta_1,\theta_2)$ to $\mathbf{y}=\{y_1,y_2\}$ based on inputs $\mathbf{x}=\{x_1,x_2\}$. To do so, I have implemented a Gauss-Newton algorithm. The function $f$ is non-linear and requires Monte Carlo methods to evaluate (I am fixing a given simulation seed to avoid noise during the fitting procedure).

I have started by writing a 1-dimensional Newton algorithm for benchmarking purposes. I have then tried to fit $\theta\in\{\theta_1,\theta_2\}$ to one point $y\in\{y_1,y_2\}$ $-$ so yielding four fitting programs in total. All four optimizations above have worked, meaning satisfactory convergence in few iterations.

Next, I have written a Gauss-Newton routine to fit $\Theta$ to  $\mathbf{y}$, and this is where my problems start. My Gauss-Newton update at the first iteration is way too large compared to what I was obtaining when fitting individual parameters. For example:

When fitting $\theta_1$ to $y_1$ using Newton, on the first iteration I was obtaining a Newton update approximately equal to $0.002$. However, under the Gauss-Newton algorithm, I get an update 3 orders of magnitude larger for $\theta_1$, i.e. approximately $-4.268$.

Similarly, when fitting $\theta_2$ to $y_2$ using Newton, on the first iteration of the Newton algorithm I get an update of approximately $2.6$, whereas under Gauss-Newton I get again 3 orders of magnitude more, $1616.5$.

I've done the following checks:

My Gauss-Newton algorithm yields the same results than my Newton one when applied to individual parameter/data-point fitting.

The four entries on the Jacobian matrix from the first iteration of Gauss-Newton, coincide with the derivatives I was obtaining on the first iteration of each one of the four single parameter-target problems above.

Using the solutions from the individual Newton fits as initializers for $\Theta^{(0)}$.

Below is the Jacobian I am getting on the first iteration to 7 decimals, as well as the initial state of the squared residuals $R^2_0$:
$$
J_{\text{iter}_1}=\begin{pmatrix}
4.998\times 10^{-4} & 1.320\times 10^{-6} \\
1.520\times 10^{-3} & 4.021\times 10^{-5}
\end{pmatrix}
\qquad
R_0^2 =\begin{pmatrix}
1.189\times 10^{-6}\\
1.048\times 10^{-6}
\end{pmatrix}
$$
such that the Gauss-Newton update is:
$$\Theta^{(1)} = \Theta^{(0)}-J_{\text{iter}_1}^{-1}R_0^2$$

I am unsure what the problem is here. Starts to look like there is a fundamental issue with applying the Gauss-Newton method, but for the problem I am tackling this algorithm should be applicable.

Any pointers as to what the problem might be or what to look for to continue debugging are appreciated.


问:
Plot a function (constant on each square=CELL) in paraview using vtk file
I'm trying to plot the solution of a numerical method using a vtk file in paraview.

The structure is very simple, and I going to explain with an example. I have a mesh of the unit square [0,1]x[0,1] consisting on 4 squares (CELLS) of equal dimensions, and then this mesh has 9 POINTS. I have a constant solution on each one of these 4 squares.

I'm trying to use the following vtk file, but does not work, What I'm doing wrong?

```<code># vtk DataFile Version 3.0
Little example
ASCII

DATASET UNSTRUCTURED_GRID
POINTS 9 DOUBLE
  0.0  0.0  0.0
  0.5  0.0  0.0
  1.0  0.0  0.0
  0.0  0.5  0.0
  0.5  0.5  0.0
  1.0  0.5  0.0
  0.0  1.0  0.0
  0.5  1.0  0.0
  1.0  1.0  0.0

CELLS 4 20
  4  0  1  3  4
  4  1  2  4  5
  4  3  4  6  7
  4  4  5  7  8

CELL_TYPES 4
  9
  9
  9
  9

CELL_DATA 4
SCALARS pressure DOUBLE
LOOKUP_TABLE DEFAULT
    0.00000
    1.00000
    2.00000
    3.00000
<span class="math-container">```
```


问:
Error $L_{2}$ convergence in Finite Element for Poisson Equation
I have written a Matlab code to solve the equation $-u'' = f$ with conditions $u(0) = u'(1) = 0$ on the domain $x \in [0,1]$. I tested the code with $f(x) = -2, \forall x\in [0,1]$. I check the plot with one of element $n=8$, and the result $u = u^{h}$ agrees:

However, I am testing the convergence error $L_{2}$ for the problem $n=8,16,32,64$. The error plot I get is not giving me the result of $O(\Delta x^{2})$ where $\Delta x$ is the fixed step size. I am trying to find if there is a bug. Why is a big change in my error between $n=16$ and $n=32$? It is probably something small in my code, but for some reason I can't figure it out. Thanks in advance!!!! Furthermore, I use a uniform grid $[0,1]$ so $\Delta x (h)$ is a fixed step size.

This is the code that I solve with Matlab and a driver to test the convergence:

```<code>% script to generate the uniform mesh for Finite element:
% dom = [a b] where a < b is the domain of our problem
function [xgrid, h] = generate1d_uniform(n)
h=1/(n-1);
xgrid=linspace(0,1,n);
end

% implement the solver for Finite Element for Poisson equation
function uh = fem1d(n,f)

% generate the uniform mesh:
[x,h] = generate1d_uniform(n);
% load vector
b = loadVector1D(x,h,f);
A = StiffMat1D(x,h);

% since u0 is 0
Ae=A(2:end,2:end);
fe=b(2:end);
uhe=Ae\fe;

% adjust the boundary condition u(0)=0
uh=[0;uhe];
end

% function to assembly load vector b:
% f is the source function; x is from the domain generate by 1d mesh
function b = loadVector1D(x,h,f)
n = length(x)-1;
b = zeros(n+1,1);

for i = 1:n
   b(i) = b(i) + f(x(i))*(h/2);
   b(i+1) = b(i+1) + f(x(i+1))*(h/2);
end

end

% function to load the stiffness matrix for A:
% the local stiffness matrix is of the form Ak = 1/h [1 -1; -1 1]
% adjust the A(1,1) = 1 and A(n+1,n+1) = 1 for boundary conditions in this
% case u'(1) = 0 and u(0) = 0
function A = StiffMat1D(x,h)
n = length(x)-1;
Ak = spdiags(ones(n+1,1)*[-1 2 -1],-1:1,n+1,n+1);
% Adjust the boundary condition u'(1) = 0
Ak(1,1) = 1;
Ak(end,end) = 1;
A = (1/h)*Ak;
end

%%% DRIVER for convergence error %%%
function driver1
clear all;
clc;
format short;
%number of elements set:
nset=[8,16,32,64];

% source function in part (a)
f = @(x) -2.*(0<=x & x<=1);

% exact solution:
u = @(x) x.^2 - 2*x;
%define the error L2 set to store the values:
errorL2=zeros(size(nset));
hset=zeros(size(nset));
% iteration:
niter=length(nset);

    for i=1:niter
        n=nset(i);
        disp(n)
        errorL2(i)=compute_error(u,n,f);
    end

% write the table for error L2:
T1=table(nset', errorL2');
T1.Properties.VariableNames ={'n elements', 'Error L2'};
writetable(T1,'TableError.csv');
disp(T1);
figure(1)
loglog(errorL2);
grid on;
grid minor;
end

% function to compute the L2 for an input of u, n and f source function
function errL2 = compute_error(u,n,f)

%solve the Poisson problem first.
[xgrid,h]=generate1d_uniform(n);
uh=fem1d(n,f);
% exact solution
uexact=u(xgrid);

% compute the error
errL2 = (h^1/2)*norm(uexact(:)-uh(:),2);
end
</code>```


问:
Help with finding an objective function for optimization
Im working in comsol and needed help figuring out the objective function for the following situation (Comsol does have multiple objective function option as min max or sum of)
I have 7 batteries connected in parallel all at the same voltage discharging. I need the current from each as uniform as possible with the maximum (the battery closest to the tab) also being as small as possible. Currently i used the standard deviation and the max of the 7 batteries as objectives and used the sum. What do you suggest?


问:
Why does this integral converge faster than normal rectangle or trapz integration?
I was looking for the fastest converging method to integrate a family of functions. After some tries, an old-school colleague suggested me a method that he used to use in excel to perform such task.

It relies on a simple procedure:

sample the function at log-spaced sample such that $x_i = x_0 \alpha^i$ for $i=0...N$

compute the values $\hat{y_i} = x_i f(x_i)$

compute the integral as $\int f(x) \approx \log(\alpha)\sum{\hat{y}} $

This "trick" is based on the following equation:

$$
\int_{x_0}^{x_1}f(x)dx=\int_{x_0}^{x_1}xf(x)\frac{dx}{x} = \int_{x_0}^{x_1}xf(x)d(\log(x))
$$
since $x_i$ is logarithmically spaced,

$d(log(x)) = const. = \log(\alpha)$.

Therefore, applying the trivial rectangles rule,

$\int f(x) \approx \sum{x_i f(x_i)\left(\log(x_{i+1})-\log(x_i)\right)} = \log(\alpha)\sum{\hat{y}} $

Now... at a first look this looks like a trick to perform an integration with the rectangle rule with logarithmically spaced point and avoiding the element-wise multiplication imposing $d(log(x)) = const.$

However, when I tested this versus more standard rectangle rules, it's convergence speed was still much faster.

The question therefore is: why is this method so fast and how does it differ from a normal rectangle-rule with log-spaced evaluation points?

To support my claim, here are the results of few test I ran.

With $f(x) = \exp(-x)$

Since $\exp(-x)$ could be a very special case in conjunction with the logarithmically spaced points, I also tried a less "special" function (that is also of the type that I need to integrate). This is: $f(x) = \frac{1}{\left(1-\left(\frac{x}{4}\right)^2\right)^2 + \left(2\cdot0.05\cdot\frac{x}{4}\right)^2}$. (Some of you might recognize this to be the dynamic amplification factor of a resonant system with natural frequency equal to 4 and damping ratio equal to 5%.)

Despite a less clear and monotonic convergence, the "magic" formula still converges faster than the rectangle rule with the same evaluation points.

Why does this happen?

Note: I initially thought this to be related with this (Integral in log-log space), but I'm now starting to think that this is a quite different topic

EDIT 1:

Following the requests of @Maxim Umansky and BlaB, here is the convergence plot in loglog scale for the two examples. It looks to me that the four methods converge with a similar power-law, but the "magic" one starts with a much lower error. (NOTE that the evaluation points for the second to fourth lines are the same)

$f(x) = \exp(-x)$

$f(x) = \frac{1}{\left(1-\left(\frac{x}{4}\right)^2\right)^2 + \left(2\cdot0.05\cdot\frac{x}{4}\right)^2}$

EDIT 2:

Please find below the python 3 code I used to generate these plots

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

fun = 'admittance'

if fun == 'exp':
    fun = lambda x: np.exp(-x)
    true_int = 1
elif fun == 'admittance':
    f0 = 4
    xi = 0.05
    fun = lambda x: 1/np.sqrt((1 - (x/f0)**2)**2 + (2*xi*x/f0)**2)
    true_int = 36.546222700939126

N = np.unique(np.round(np.logspace(1, 8, 1000)))

#%% plot function
plt.figure()
x = np.linspace(0.001, 100, 200)
y = fun(x)
plt.plot(x,y)
plt.xlabel('x')
plt.ylabel('f(x)')

#%% open new figure
plt.figure()

#%% Ractangle w/ linspaced x
err = []
for n in N:
    x = np.linspace(0.001, 100, int(n))
    y = fun(x)
    dx = x[1]-x[0]

    err.append(np.sum(y)*dx - true_int)

plt.loglog(N, np.abs(err), label='Ractangle w/ linspaced x')

#%% Ractangle w/ logspaced x
err_left = []
err_right  = []
for n in N:
    x = np.logspace(-3, 2, int(n))
    y = fun(x)
    dx = np.diff(x)

    err_left.append(np.sum(y[:-1]*dx) - true_int)
    err_right .append(np.sum(y[1: ]*dx) - true_int)

plt.loglog(N, np.abs(err_left), label='Ractangle w/ logspaced x (left)')
plt.loglog(N, np.abs(err_right ), label='Ractangle w/ logspaced x (right)')

#%% ACA
err_ACA = []
for n in N:
    x = np.logspace(-3, 2, int(n))
    y = fun(x)
    xy = x*y
    dlogx = np.log(x[1]) - np.log(x[0])

    err_ACA.append(np.sum(xy) * dlogx - true_int)

plt.loglog(N, np.abs(err_ACA), label=r'<span class="math-container">$\log(\alpha)\sum{xf(x)}$</span>')

#%%
plt.legend()
plt.xlabel('Number of integration points')
plt.ylabel('|Error|')
</code>```


问:
Why does Eigen allocate a temporary to evaluate A.noalias() = B.transpose()*C in parallel?
I wrote a program which iteratively transforms data using matrix multiplications. To minimize the number of large memory allocations, I use two roughly equal-sized std::vector<double>'s data and temp_data to store the current data and the new data respectively. For the matrix operations, I wrap each buffer into an Eigen::Map<Eigen::MatrixXd>, called data_matrix and temp_data_matrix respectively. Then, I need to compute the following operation:

```<code>temp_data_matrix.noalias() = data_matrix.transpose()*factor;
</code>```

where factor is simply an Eigen::MatrixXd which is much smaller than data_matrix or temp_data_matrix. Thanks to .noalias(), this line does not allocate any large temporaries normally. However, when multi-threading the program suddenly allocates a temporary (roughly) equal to the size of temp_data! In practice this can increase memory usage by around 50% so I'm curious about a few things:

Why does this only happen when multi-threading? What algorithmic aspect requires a temporary buffer on top of the output buffer?

Why does this not happen with the line temp_data_matrix.noalias() = factor*data_matrix? After all, my operation seems more cache-efficient: because we transpose data_matrix, we simply need to compute the scalar product of each pair of columns from data_matrix and factor, which should be fast considering both matrices are column-major. However, I'm not considering the effects of blocking here, which I know can be very important for performance.

What would be the fastest way to implement the aforementioned line without allocating any temporaries of the order of data or temp_data? In case the parallel algorithm needs some "scratch space", can't it use temp_data_matrix for that? I don't care about copying factor because it's relatively small in practice. Right now the fastest solution seems to be to switch back to single-threading for this line, but I hope this can be improved. I tried using lazyProduct but this turned out to be much slower than the single-threaded solution.

Appendix: The observations above were made using Massif (valgrind) both in my original program and with the example code below. I compiled the example program using the command g++ -O3 -DNDEBUG -I /usr/include/eigen3 -fopenmp eigen_memory_test.cpp -o eigen_memory_test.

```<code>#include <iostream>
#include <Eigen/Dense>
#include <algorithm>
#include <random>
#include <cstring>

int main(int argc, char** argv) {

    bool parallel = (argc > 1 && strcmp(argv[1], "parallel") == 0);
    std::cout << "parallel: " << (parallel  ? "true" : "false") << std::endl;
    bool lazyproduct = (argc > 2 && strcmp(argv[2], "lazyproduct") == 0);
    std::cout << "lazyproduct: " << (lazyproduct  ? "true" : "false") << std::endl;
    bool test = true;

    // Initialize parallellism
    int threads = 4;
    Eigen::setNbThreads(threads);

    // Initialization and memory allocation
    size_t columns = 100000;
    std::vector<size_t> rows_per_step = {100, 90, 80, 70, 60, 50};
    size_t steps = rows_per_step.size() - 1;
    std::vector<double> data(rows_per_step[0]*columns);
    std::vector<double> temp_data(rows_per_step[0]*columns);

    // Generate random data
    std::random_device rnd_device;
    std::mt19937 mersenne_engine {rnd_device()};
    std::uniform_real_distribution<double> distribution(-1.0, 1.0);
    auto generator = std::bind(distribution, std::ref(mersenne_engine));
    std::generate(data.begin(), data.end(), generator);

    // If testing, print first elements to verify later
    if (test)
        std::cout << "data[:3]: [" << data[0] << ", " << data[1] << ", " << data[2] << "]" << std::endl;

    // Process steps
    for (size_t step = 0; step < steps; step++) {

        size_t initial_rows = rows_per_step[step];
        size_t new_rows = rows_per_step[step + 1];

        // Choose factor
        Eigen::MatrixXd factor = (test ?
            Eigen::MatrixXd(Eigen::MatrixXd::Identity(initial_rows, new_rows)) :
            Eigen::MatrixXd(Eigen::MatrixXd::Random(initial_rows, new_rows))
        );

        // Transform data
        Eigen::Map<Eigen::MatrixXd> data_matrix(data.data(), initial_rows, columns);
        Eigen::Map<Eigen::MatrixXd> temp_data_matrix(temp_data.data(), columns, new_rows);
        if (!parallel)
            Eigen::setNbThreads(1);
        // In my real problem setting, the transpose in the following matrix multiplication is actually necessary
        if (lazyproduct)
            temp_data_matrix.noalias() = data_matrix.transpose().lazyProduct(factor);
        else
            temp_data_matrix.noalias() = data_matrix.transpose()*factor;
        Eigen::setNbThreads(threads);

        // Transpose back to original data
        // This doesn't happen in the real problem, but I put this here to keep this example simpler
        Eigen::Map<Eigen::MatrixXd> transposed_data_matrix(data.data(), new_rows, columns);
        transposed_data_matrix.noalias() = temp_data_matrix.transpose();

    }

    // Print first elements to make sure matrix operations aren't optimized away
    std::cout << "data[:3]: [" << data[0] << ", " << data[1] << ", " << data[2] << "]" << std::endl;

    return 0;

}
<span class="math-container">```
```


问:
The physical meaning of conservative mass in diffusion equation
I am working on 1-D mass transient diffusion in a radial domain (spherical object) using finite volume method. My equation reads
$$\frac{\partial C}{\partial t} = D\left[\frac{\partial^2 C}{\partial r^2} + \frac{2}{r}\frac{\partial C}{\partial r}\right]\, .$$

In a sense of conservative mass within the the radial domain, plus a closed BC ( constant concentration of A species on the surface), does that mean I should end up with concentration gradient among the discretized nodes? In other words, at first time step (forward in time) for internal nodes, Should I get distributed concentrations that the sum of whole internal nodes will give me the concentration on the surface?


问:
Galerkin Least-Squares stabilization for FEM solution advection (hyperbolic) equations
I am playing with Galerkin Least-Squares stabilization to solve advection diffusion problem in the context of the finite element method. This works very well for steady-state advection-diffusion equation, even in the High Péclet limit, when the problems are solved in the steady state.

The equation we solve is :
$$\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = D \nabla^2 \alpha  $$

with $ \mathbf{u}$ a velocity vector field which is known and $\alpha$ a scalar. $D$ is a diffusion coefficient.
The weak form (without boundary terms) is :
$$ \int \phi_i \frac{\partial \alpha}{\partial t} d\Omega + \int \phi_i \mathbf{u} \nabla \alpha d\Omega + \int D \nabla\phi_i \cdot \nabla\alpha d\Omega = 0$$

Adding GLS stabilization leads to the following:
$$ \int \phi_i \frac{\partial \alpha}{\partial t} d\Omega + \int \phi_i \mathbf{u} \nabla \alpha d\Omega + \int D \nabla\phi_i \cdot \nabla\alpha d\Omega + \int (\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha - D \nabla^2 \alpha ) (\tau \mathbf{u} \cdot \nabla \phi_i) d\Omega = 0 $$

Where the additional term that is added effectively consists in the strong form of the equation. We define the stabilization coefficient $\tau$ as:

$$\tau = \left( \left(2 \frac{|\mathbf{u}|}{h} \right)^2+ \left(4 \frac{\mathbf{D}}{h^2} \right)^2 \right)^{-\frac{1}{2}}$$

where $h$ is the size of the cell which is calculated using the diameter of a sphere of equal volume to that of the cell. Consequently $\tau$ has units of time.

Solving the steady-state version of the equation :
$$\mathbf{u} \cdot \nabla \alpha = D \nabla^2 \alpha  $$

works perfectly fine for any Péclet number. Obviously the GLS term introduces artificial diffusion, but it ensures that the value remain strictly positive and bounded. As the mesh is refined, the solution gets better and better. My issue appears in transient simulations. There, I get loss of positivity and oscillations (overshoots) in my results. It is nothing too bad (nothing like regular Galerkin would give me), but for example, the advection of a step-function without diffusion leads to overshoot of alpha and the bounded character is not preserved.

I feel like something is missing for my GLS stabilization when my equations are transient, but I  have not been able to find literature that discuss this problem exactly. What am I missing? Is it because of how I define my $\tau$? For information, I use BDF-type schemes (BDF1 or BDF2) for time integration.

Edit 1:
Following the (Finite Element Method : Theory, Implementation and Application by Larson) I understand that a full GLS formulation would require me to add another term that would be :

$$ \int (\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha - D \nabla^2 \alpha ) (- D \nabla^2 \phi_i)$$

but this term is trivially zero for P1 elements, so I do not think it is a necessary component here. At least this is what I saw in the source below:

Edit 2

Results obtained for steady-state advection-diffusion at Pe=1000 with a coarse mesh

Advection of a sharp tanh function. Initial condition and after first step (note the overshoot)

Advection of a step function. Initial condition and after first step (note the overshoot)


问:
I need help developing a process to map many individual records (from a database) into larger groups
There are many rows of data in a database table, each of which contains many columns.
ID  (numeric - 12345, 12378, etc)
Name 1 -  AAA,BBB,CCC
Name 2 -  dd, ee, ff
Name 3 - Ggg, Hhh, Iii

The objective is to come up with rules based on a combination of these  values to map to a set of groups.
Group ID    1 - 1,000

There is a large mapping table that applies rules to the individual records and maps each one to it's group ID.  The code is spaghetti and has gaps and duplications, and I would like a better approach if you can suggest one.

Here is a small example:

if Name 1 in (...) then  use this process below:

if ID = 856743 THEN Group ID = 2000

if ID in (1,2,3... 200) THEN Group ID = 2001

if ID = 377 OR Name 2 in (X,Y,Z)  THEN Group ID = 2003

if Name 2 in (U, V, W)  THEN Group ID = 2007

Else Group ID = 2011

I am thinking that I need to take all the column permutations and make them dimensions.
The problem is there are many segments of mapping like this:
if a column value is in (X,Y,Z) and the ID is one of (a large set) then the value is ...

From time to time more IDs are added, and who knows if new column values are added.

Should I just have a mapping table listing:
ID, and each dimension, and what the Group ID is?
That would necessitate adding new mapping rows when new IDs are added.  The current process sometimes handled that automatically because it just takes a record that matches Name 1 and Name 2,

ID  N1 N2 Group ID
1   1   1     10
2   1   1     11
2   1   2     12

2   1     13

2   2     14

Thoughts?


问:
Why is bounding a surface with a capsule is better than with a cylinder to detect intersections?
In this article:
https://www.geometrictools.com/Documentation/IntersectionOfCylinders.pdf

the writer says: "If you plan on using cylinders
for bounding volumes in a real-time graphics engine—think twice. A better alternative to a cylinder is a
capsule."

Well, I am planning on bounding surfaces with cylinders (actaully I already have). What are the reasonings behind this statement?

Notice - I am mainly interested if two surfaces do not intersect. Some sort of a quick rejection method by bounding the surfaces with solid cylinders.


问:
Stability plot of upward difference implicit time
I am analyzing the stability of 1D convection (advection) equation as shown in the picture. When I derive the equations as shown I want to get rid of the complex number.

I`m asking if those stability plots make sense or not

I want to get rid of a complex number


问:
Check if LinearOperator is symmetric
I have a scipy.sparse.linalg.LinearOperator object. I'd like to check if its associated matrix is symmetric without actually instantiating the matrix in the most computationally efficient way possible. Does anybody know a good way to do this?

Title: Check if LinearOperator is symmetric


问:
Tool that calculates total drag for a submarine
Can you please suggest a tool (software or online tool) that can analyze a submarine hull shape and give the total drag at a specific velocity.

A simplified tool would work too, where I select a general geometry, then input the dimensions and speed, and get the drag force.

I'm reading Submarine Hydrodynamics by Martin Renilson, where he writes "For a typical submarine the total submerged pressure resistance does not exceed 10–20% of the skin friction resistance. Thus, the skin friction dominates the resistance of a deeply submerged submarine."

I'm assuming that this is the case for slow speed only, and that pressure resistance is the dominating resistance at high speed, or submarine hulls wouldn't be so streamlined. Or am I wrong?

My application is for low speed, and I need a non conventional hull shape that is flat (box shaped), and I'm trying to understand how much penalty I will have to pay for this.

The standard drag equation is too simplified in my case, because it approximates both skin drag and pressure drag into one variable (the drag coefficient).


问:
A priori error estimates - finite element method - mixed boundary conditions
Consider the problem
$$  \left\{\begin{array} {rcl}
-\Delta u & =  0 & \text{ in } \Omega \\
u & =  0 & \text{ on } \Gamma_D \\
\frac{\partial u}{\partial n} &= g &\text{ on } \Gamma_R \\
\frac{\partial u}{\partial n} &= 0 &\text{ on } \Gamma_N
\end{array}\right.  $$
Where $g$ is a trace of a $H^{1+\varepsilon}(\Omega)$ function. An example of such a configuration in is shown in the Figure below.

I am aware that depending on the function $g$ the solution to this problem might be less regular than $H^2(\Omega)$. I am interested in reading about the a priori error estimates that can be obtained when solving this problem with triangular piecewise linear finite elements (P1, Lagrange).

My questions are:

Is it possible to indicate a reference which treats in detail this classical case?

What is the best convergence result that can be expected in this case in terms of the regularity properties of $u$? (i.e. what is the best exponent one can get in the estimate $\|u-u_h\|_{H^1(\Omega)}\leq Ch^\gamma$ when $u$ is not necessarily in $H^2$)


问:
FEM port Boundary definition for electromagnetics and wave guides
We are currently in the process of implementing ports in our EM FEM simulation SW.
We have come across the definition of boundary conditions for the ports, and we do not understand the equation for the boundary. In the book “The Finite Element Method in Electromagnetics” they start from the understandable idea of it being a super position of the incoming and reflected wave. And because it is a Homogeneous waveguide the field follows the analytical TE10 mode.

To find the equation for the field they take the curl of both sides. We do not understand how this makes sense in a physical way. The $\hat{n}$ vector is the normal to the waveguide port. We assume that this is used to isolate the field along the port surface.
Do you have any suggestions on how and why this is done?


问:
How does the error work for the Strang Splitting?
We know in Strang splitting that the splitting error in the steady state solution is proportional to $h^2$. I want ask 2 things:
If this error in the steady state solution is the global error?
If we are in in a previous step before the steady state solution for example $t=0.3$ sec (with $t_{total}=1$ sec) will be the splitting error proportional to $h^2$ or it changes ?
Thanks!


问:
Optimization algorithm to find coinciding root and minimum
Are there any optimization algorithms aimed at finding a coinciding root and a (local) minimum of a multi-variable function f. Say it is known analytically that these should coincide for the function.

I am interested in any such algorithm but more specifically for the case where f is continuous, at least twice differentiable and very expensive to compute (so internal costs for the algorithm can be neglected, only function calls matter). What would you use if the gradient could be computed for free? What if the Hessian can also be computed for free? What if the root and the minimum only coincide approximately, i.e. the minimum $f(x_c)\in $ $[-\epsilon,\epsilon]$ for some small epsilon.


问:
Problem with Numerical DAE Implementation of Nernst-Planck Equations
I write to you all today in order to clarify if I implemented my model correctly into MATLAB.

I will outline the following: analytical governing equations, my discretization for the equations and boundary conditions, and outline the implementation into ode15s. Ultimately, my problem with implementation is that steady-state is not achieved as I reach an infinite accumulation with and without a volumetric rate term. I want to make sure that my implementation is correct before reevaluating the governing equations.

System description - arbitrary cation and anion react with gaseous species in porous electrode. Cell potential is specified. BV kinetics is utilized but not included here. Concentrations are mol/m^3. Fluxes are mol/m^2/s. Volumetric charge is A/m^3.

The gas balances work and are not included.

Please see the images attached.

Governing equations:

Discrete form of the equations:

How to get boundary condition for anionic species:

Numerical solving procedure - method of lines with ode15s. Finite differences with second order central and first order central. Electroneutrality equation applies but is not used in DAE solver. Total volumetric charge conservation is used as algebraic constraint. To deal with the imaginary points , I solve for the cation, anion, and potential value via the boundary values.

5 equations (3 listed + 2 not-listed) / node. n refers to number of nodes. Length of cell is on the range of a 100 microns.

Three big questions:

The boundary condition for the anionic species (i.e. a) is derived utilizing the electroneutrality as applied to the cationic boundary condition. Is this sufficient to specify the boundary value problem?

I have a 5xNode# by 5xNode# mass matrix with 1 for equations 1-4 on the diagonals and 0 for equation 5. Is this mass matrix constructed correctly?

The arbitrary reference potential is set as 0 V. I specify in my mass matrix a value of 1 for the charge balance as I set the potential. I set my yt for the charge balance at node n = 0 as the initial value is set as the boundary condition. Is this causing an instability problem or should the mass matrix be set to 0 for this equation?

Thanks for any help!


问:
Using event function to shift between 2 sets of ODE
I have a bunch of ODEs I am trying to solve using ode45 in MATLAB. I have hidden the details of the equations to keep it simple (so as to build a general algorithm)!

\[\frac{dR}{dt}= F_{1}(R,N) \quad\text{........................(Eq.1)}\]
\[\frac{dN_{1}}{dt}= F_{2}(t) \quad\text{.............................(Eq.2)}\]
\[\frac{dN_{2}}{dt}= F_{3}(R,N) \quad\text{.......................(Eq.3)}\]
\[\frac{dN}{dt}= \frac{dN_{1}}{dt} \quad \text{or} \quad \frac{dN_{2}}{dt} \quad \text{based on a criteria of } max(-\frac{dN_{1}}{dt},\frac{dN_{2}}{dt}) \]

So, basically I want is, to solve for R, N for 2 cases. The ODE solver should stop when the criteria is met and should restart with Case 2 by solving a different ODE set.

\[\text{Case 1}:  \frac{dR}{dt}; \frac{dN_{1}}{dt}\]
\[\text{Case 2}:  \frac{dR}{dt}; \frac{dN_{2}}{dt}\]

I know I need an event function to do this. Unfortunately, I am new to this for which your help is solicited (to check if my code is functioning properly). I also can't think of any test functions to test the algorithm. Let me know if any other info is required.

```<code>clc
clear

y0   = [0.1,0.1,0.1,0.21]; %random initial nos.
t0   = 0;                  %random no.
tEnd = 10;                 %random no.

opt = odeset('RelTol', 1e-3);

y = y0;
t = t0;

State = 1;

while t(end) < tEnd
    yprime = myFunction(t, y, State);

    if yprime(2) > yprime(3)
        State = 1;
    else
        State = 2;
    end

    fcn = @(t,y) myFunction(t, y, State);

    opt = odeset('Events', @(t, y) myEvents(t,y,State));
    [at, ay] = ode45(fcn, [t(end), tEnd], y(end, :),opt);

    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
end

plot(t,y(:,4))

function [yprime] = myFunction(t, y, State)
yprime = [5*y(1); 4*y(4) + 77*t ; 30*y(4) - 7*t]; %random funcs.
if State == 1
    yprime(4) = yprime(2);
    disp('state is 1')
else
    yprime(4) = yprime(3);
    disp('state is 2')
end
end

function [val,isterm,dir] = myEvents(t,y,State)

fun = myFunction(t, y, State);

val      =  fun(2) - fun(3);
isterm   =  0;
dir      =  0;

end
</code>```


问:
Stability of Crank-Nicholson for advection diffusion equation for spatial discretization other than finite differences second-order centered
Crank Nicholson is a time discretization method (see 4th equation here). From what I see around, you can use different space discretization, such as Finite elements. But for the linear advection-diffusion equation, I see on books the von Neumann stability proof only when second-order centered spatial finite differences are used for both first and second-order derivatives (see e.g. here). My questions are:

If I use upwind, first or higher-order, higher-order centered differences, or finite elements, is the scheme always unconditionally stable for the linear advection-diffusion equation? Any reference?

Does the unconditional stability of Crank Nicholson also hold for only advection or only diffusion, and for any physical Peclet number?


问:
Comparison between stability and accuracy of various Finite Difference schemes
I`m Analyzing the 1D convection equation (PDE) for stability, consistency, and accuracy.

I know Both upwind schemes (explicit and implicit) show better stability with a higher number of waves for courant number less than one. The upwind scheme, in general, overcomes the problem of stability in the centered difference scheme. However, the upwind schemes scarify some accuracy to gain stability.

Now I`m confused about the stability and accuracy of both the Centered difference – trapezoidal and Lax Wendroff Schemes.

I have attached photos of the analysis I conducted but not sure how to comment. I`m also confused about the order of schemes of:

a) Backward difference – explicit time b) Backward difference – implicit time c) Centered difference – trapezoidal time d) Lax Wendroff

In the attached graphs Beta is courant number, N-lambda is no. of intervals per wave length, r is the amplification factor.


问:
Unexpected CUDA processing time dependency on thread count
When calling a kernel function, the number of threads per block should ideally be a multiple of the warp size. This yields more efficient use of resources and lower processing times. However, there seems to be another factor that periodically decreases the processing time. Shown below, the processing time gets offset every 32 threads, while an additional speed-up takes place at every multiple of 11 threads per block. What are the reasons that could lie behind this behaviour?

The GPU in question is the GeForce GT 730, running the kernel function attached at the bottom of this post. For timing purposes, it is invocated in a loop using:

kernel_generate_image[(16,16),(1,i+1)](px, 32)

where px = np.zeros([1024,1024])

```<code>@cuda.jit
def kernel_generate_image(image, T):

    # Calculate the thread's absolute position within the grid
    x = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x
    y = cuda.threadIdx.y + cuda.blockIdx.y * cuda.blockDim.y

    # Set stride equal to the number of threads we have available in either direction
    stride_x = cuda.gridDim.x * cuda.blockDim.x
    stride_y = cuda.gridDim.y * cuda.blockDim.y

    for i in range(x, image.shape[0], stride_x):
        for j in range(y, image.shape[1], stride_y):
            image[i, j] = (sin(i*2*pi/T+1)*sin(j*2*pi/T+1)*0.25)
</code>```


问:
tiny typo in Numerical Recipes Eq. 9.4.6
The Numerical Recipes Forum http://numerical.recipes/forum/ is closed, so I will record a tiny typo here for the benefit of others who may wonder about this.  (This typo is not in the software, but only in the text, and is currently in the online version of the textbook at http://numerical.recipes/book):

In Eq. 9.4.6 (describing the convergence rate of the Newton-Raphson algorithm), the "-" sign on the RHS should be removed.  This can be easily verified by substituting 9.4.3 into 9.4.5 and expanding for small epsilon_i to keep the lowest order non-zero term.  One can also see from Fig. 9.4.1 that if f' and f'' are positive then a Newton method converges from the right (positive errors).

A possible source of confusion in looking at other references is that some people define an error eps_j in an estimate x_j as eps_j =  x - x_j, where x is the exact answer, but the more traditional definition (as used in most statistics references) is eps_j = x_j - x.  Numerical Recipes uses the more traditional definition in Eq. 9.4.3, so the "-" sign in 9.4.6 is wrong.


问:
Neighbor lists in overdamped Langevin MD simulations
The following question is about the implementation of molecular dynamics simulations with neighbor lists.  In cases where the paths of individual particles are (i) stochastic and (ii) somewhat rough, without well-defined velocities, it seems reasonable (for performance' sake) to maintain neighbor lists in a way that may allow occasional unphysical overlaps between particles, but only with a low probability, and that otherwise prevents the vast majority of particles from overlapping.  Particles may be capable of moving distances of order $N\sqrt{\delta t}$ in $N$ steps, but owing to randomness are typically confined to a ball whose radius is of order $\sqrt{N\delta t}$.  In this case, would it be possible to reduce overlap error even further with a 'replay' step, where configurations with unphysical overlaps are discarded by restoring the simulation to an earlier snapshot, or could doing so introduce subtle statistical anomalies that would be difficult to address?  The main risk seems to be in systematically undersampling orbits that are poorly represented by the neighbor list; that is, a (possibly small but nontrivial) portion of orbits, in which two particles approach each other closely from a moderate distance within a short time scale, would be omitted.  Is this actually the case?  And if so, are there additional and possibly more subtle fallacies being swept under the rug?  (Consider, for example, the opposite limit where the neighbor list is rarely updated at all, unphysical overlaps are relatively common, and almost all orbits are omitted except those where all particles avoid each other.)


问:
Why is this Advection-convection model with insulating boundary losing mass?
I am trying to model a 1-d advection-convection numerically, using an upwind scheme. I'm using the following equation to calculate the value of internal cells:

$$C_x^{t+1} = C_x^{t} + D\frac{\Delta t}{\Delta x^2} ( C_{x-1}^t-2 C_{x}^t+C_{x+1}^t) + u\frac{\Delta t}{\Delta x}(C_{x+1}^t-C_{x}^t)$$

and I'm using the Robin boundary condition (a.k.a. insulating boundary). In other words, I'm trying to guarantee no mass will leave the model environment. The boundary value is calculated by setting the leftmost and right most cells as:

$ C_1^{t+1} = D \frac{C_2^{t}}{u+D} $ and $ C_n^{t+1} = C_{n-1}^t\frac{(u\Delta x+D)}{D} $

Where $C$ is the concentration, $D$ the diffusion coefficient, $u$ is the velocity and $x$ is the current cell x-coordinate which goes from $1$ to $n$.

I understand this scheme should be numerically stable and accurate when  both the following conditions are met:

$ u\frac{\Delta t}{\Delta x} + 2D\frac{\Delta t}{\Delta x^2} \leq1 $ and $\frac{u\Delta x}{D} < \frac{2}{1- u\frac{\Delta t}{\Delta x}} $

However, if I run this model with parameters $D=0.05$,$u=0.1$,$\Delta t=0.001$,$\Delta x=1$, there is some mass loss after 10,000 time steps, even tough the conditions above are met. It may seem a lot, but since my timesteps are very small it means less $t<10$ when that starts happening.

Any idea, of what could cause this loss of mass?

Below there's code to conduct this simulation in R.

```<code>transfun = function(time,state,parms) {
    D = parms<span class="math-container">$D
u = parms$</span>u
    dt = parms<span class="math-container">$dt
dx = parms$</span>dx
    n  = length(state)

    new = numeric(n)

    #Boundary cells
    state[1]= state[2]* D/(u+D)
    state[n]= state[n-1]* (u*dx+D)/D

    # Internal cells
    for(a  in 2:(n-1)) {
        state[a] = state[a] +D*dt/dx^2 * ( state[a-1] + -2*state[a] + state[a+1]) +
                 u * dt/dx * (-state[a] + state[a+1])
    }

    #Boundary cells
    state[1]= state[2]* D/(u+D)
    state[n]= state[n-1]* (u*dx+D)/D

    return(list(state))
}

## Hand made loop
# Ensure parameters follow the stability rule: 0<= C^2 <= 2s <= 1
# where C = u*dt/dx and s = D*dt/dx^2
D=0.05
u=0.1
dt=0.001
dx=1
C = u*dt/dx
s = D*dt/(dx^2)
n = 15
T = 10000
print(C+2*s < 1)
print(u*dx/D<2/(1-C))

resul = matrix(NA,T,n)
resul[1,] = rep(0,n)
resul[1,ceiling(15/2)] = 1
for(a in seq(2,T,1)) {
    resul[a,] = transfun(1,state=resul[a-1,],parms=list(D=D,u=u,dt=dt,dx=dx))[[1]]
}
</code>```

,


问:
Confused on how Method of Manufactured solutions works?
I am new to computational science and I am trying to wrap my head around how MMS works. I am solving the time independent Helmholtz equation as a simple test of the technique

so my starting equation is following with homogeneous boundary conditions on a line segment

$\Delta u +k^2u = F $

So I am supposed to take something that is sufficiently differentiable. So, following other people I pick the following for my 1D case

$ \bar{u} = \sin(\pi x)$

My new source term created from MMS would be

$G(x) = -\pi^2\sin(\pi x) + k^{2}\sin(\pi x) - F $

Now here is where I am stuck. I picked a u and now I have a new F. Somehow I am supposed to be able to verify my rate of convergence of my actual problem with this fictitious one. Are the points the same?

Any advice would be appreciated.


问:
Continuum solutions for the Dirac equation in Coulomb potential - numerical codes
Following the representation used in [1, pag. 11] the solution of the Dirac equation in polar coordinates for energy $E$ is of the type:

$$ \psi_{E\kappa m}(\bf{r})= \dfrac{1}{r} \Bigg( \begin{matrix}
P_{E\kappa} (r) \chi_{\kappa}^m(\theta,\phi)\\
iQ_{E\kappa} (r) \chi_{-\kappa}^m(\theta,\phi)\\
\end{matrix}\Bigg)\ ,$$

where $P$ and $Q$ represent the large and small radial components, respectively. While $\chi_{\kappa}^m(\theta,\phi)$ is the spherical spinor function.

The book gives the solutions for the radial equation in terms of the solution of the Kummer's confluent hypergeometric equation, that is:

$$\dfrac{d^2Y(\rho)}{d\rho^2}+(b-\rho)\dfrac{dY(\rho)}{d\rho}-aY(\rho)=0$$.
Thus, the large and small components are given as:

$$P_{E\kappa}\propto\rho^{\gamma}e^{-\rho/2} [X(\rho)+Y(\rho)]$$

$$Q_{E\kappa}\propto\rho^{\gamma}e^{-\rho/2} [X(\rho)-Y(\rho)]$$

where:

$$X(\rho)\propto \biggl(aY(\rho)+\rho\dfrac{dY(\rho)}{d\rho}\biggr).$$

and $\rho$ is proportional to the radius.

Depending from the type of solution searched (bound or continuum) $P$ and $Q$ assume different shapes.
The code GRASP2K [2] provides the possibility to calculate the Dirac bound wavefunctions, implementing the multiconfiguration Dirac-Hartree-Fock method.

I would like to know if anyone knows a code to calculate the continuum wavefunctions for a given atom or I should try to implement the analytical solutions.

References

[1] Relativistic quantum theory of atoms and molecules, by I. Grant.

[2] https://www-amdis.iaea.org/GRASP2K/


问:
Does shift-invert method has invertibility issue?
Please note that I have nearly zero background on numerical methods.

I understand the shift invert method as described in SciPy Tutorial

The main argument of the above link is as follows. Suppose we want to solve a generalized eigenvalue problem $A\mathbf x = \lambda M \mathbf x$. Then, the shift invert method transforms the problem into
$$(A-\sigma M)^{-1} M\mathbf x = \frac{1}{\lambda-\sigma} \mathbf x.$$

My question: What if $A - \sigma M$ is not invertible? In particular, I am interested in the case $M=I$ and $\sigma=0$. Then, the above equation reduces to $A^{-1}\mathbf x = \lambda^{-1} \mathbf x$. But, I know that $A$ has zero as an eigenvalue, so the above-mentioned problem happens indeed!

The motivation of question is as follows. Using eigsh function in SciPy, I obtained that the eigenvalue with the smallest magnitude is of order $10^{-4}$, but it is supposed to be zero!

I used $\sigma=0$ in the eigsh function parameter, and my matrix is sparse Hermitian with size $10000\times 10000$.


问:
Solution of an underdetermined system stemming from a PDE with Neumann BC
Consider the Poisson's equation in 1D with homogeneous b.c.'s $\mathrm{d} \phi/\mathrm{d} x=0$ with the seven point Laplacian (1 -54 783 -1460 783 -54 1 / 576 on a uniform grid). The resulting system is underdetermined because the solution is known only up to an additive constant. It only exists if the compatibility condition is fulfilled (sum of the RHS is 0).

When being solved with LAPACK (dgbtrf and dgbtrs) on the uniform grid an accurate solution is recovered even when the system is left underdetermined.

However, I am having problems with the non-uniform grid (the second derivative in difference form is derived by applying a discrete gradient to field of a discrete gradients in intermediate grid points).

When the grid is perturbed, the very same procedure fails. It either returns NaN or the solution is so large in magnitude that there is a large numerical error in each point. E.g. this system:

$\frac{1}{576\overline{\Delta x}^2} \begin{pmatrix}-624.5 &&675 && -50.5 \\
1350 &&-2700 && 1350 \\
-50.5 && 675 && -624.5\end{pmatrix} \begin{pmatrix}x_1 \\ x_2 \\ x_3\end{pmatrix} = \begin{pmatrix}0.216 \\ -0.864 \\ 0.216\end{pmatrix}, \overline{\Delta x} = 1/3$

When $n=3$ and grid is symmetric, like here, it is enough to add another equation that fixes the solution in one point, e.g., $x_1 = 0$. For that I just multiply $A_{1,1}$ by 2. I get the solution

$\frac{1}{576\overline{\Delta x}^2} \begin{pmatrix}-1249 &&675 && -50.5 \\
1350 &&-2700 && 1350 \\
-50.5 && 675 && -624.5\end{pmatrix} \begin{pmatrix}0 \\ 0.02048 \\ 0\end{pmatrix} = \begin{pmatrix}0.216 \\ -0.864 \\ 0.216\end{pmatrix}, \overline{\Delta x} = 1/3$

that is also a solution of the original system.

However, when the grid is asymmetric or just larger, the same procedure yields solution that is accurate, except for the value of $x_1$, where the residue is typically 1e-4 or a similar number a few orders more or less.

The RHS is adjusted to zero mean by subtracting the weighted average

$b_0 = b - \frac{\sum_1^n b(i)\Delta x_i}{\sum_1^n \Delta x_i}$.

What is the proper procedure to solve such a system? Why do I get the non-zero residuum after adding $x_1=0$ to the first row in a non-trivial case?

It is hard to show a complete example due to rounding, but something like

$\frac{1}{576\overline{\Delta x}^2} \begin{pmatrix}-535.58897243  && 582.46753247 &&  -46.87856004 \\
774.43609023 && -1671.42857143  && 896.99248120 \\
-57.92258424 &&  895.23809524 && -837.31551100\end{pmatrix} \begin{pmatrix}x_1 \\ x_2 \\ x_3\end{pmatrix} = \begin{pmatrix}0.198 \\ -0.852 \\ 0.588\end{pmatrix}, \overline{\Delta x} = 1/3$

giving, after the modification,

$x^T = \begin{pmatrix}1.0558979027522729E-004 &&2.0058136568986706E-002 &&-2.3505247184514636E-002\end{pmatrix}^T$
with residuum

$\begin{pmatrix}8.8363636363628095E-004&& 0 &&0\end{pmatrix}^T$.

Maybe the non-zero $x_1$ is a que?


问:
High quality constrained optimization C++ library with matrix free second order solver?
I'm working with large scale constrained optimization problem. Some of my constraints can be non linear. Currently i'm using IPOPT. Quality is good by my Hessian computation too slow. It seems that i can implement hessian-vector problem much more effective. Unfortunately i cannot find how to use matrix free solver with IPOPT.

Can anyone suggest me good library which supports matrix free solver? Preferably in C/C++/fortran. Currently i only found SCIPY trust-constr solver but it's purely python.


问:
Applications where random number generation is a speed bottleneck
Pseudorandom numbers are used massively in Monte Carlo methods, e.g. in physics simulations. The fastest random number generators, such as the xoroshiro or xorshiro families by Vigna (see eg here) are extremely fast, taking approximately three CPU cycles for generating a 64-bit random integer. Using parallel instances and the SIMD registers present in modern CPUs (SSE, AVX) the speed can be further increased by a factor 2-3.

I would like to know if there are any real-world applications for which random number generation is a speed bottleneck (when using a fast generator). It seems to me that modern generators are so fast that almost any processing of the generated random data is likely to take longer than the generation step itself.


问:
HLLC Riemann solver with shock test 2 - extension to low densities
I am currently using the HLLC solver to solve a 1-D system of Euler equations with very satisfactory results.

However, there are cases where my solution produces low-density, high velocity states, which lead to strange heating at the origin. The heating is strange, as I am running with a $\gamma_{adiabatic}$ close to 1, hence the simulation is quais-isothermal.

I've double and triple-checked my boundary conditions, and am fairly certain its not that. My other idea is that this situation is resembling the anomalous shock tube test-2 heating, as seen in HLLC solvers in general (circled in red in the image below).

Hence my question:

As the HLLC solver is otherwise a nearly perfect solver, I was curious whether there is a fix for this near-vacuum overheating. I am mainly asking for a reference.

Toro's references in the HLLC chapter turned out interesting papers, but not a fix for this specific issue. I am also aware of the keyword 'drying and wetting' for near-vacuum problems with the shallow water equations, but searching for this in conjunction with the Euler system and fixing the overheating amounted to nothing. Therefore I am turning to the community.


问:
Insert a boundary condition without removing periodicity assumption
I've to perform a multi-fluid internal flow simulation with a code which intrinsically assumes periodicity on a given direction (say z) on a cartesian grid. Nonetheless, the problem I'm trying to solve does not involve this kind of condition, in fact the fluids are running out from the computed domain through the two z extreme edges (one of them at sonic speed, the other at almost sonic speed) after being generated inside it from a given constant source. The simulation evolves the two fluids profiles up to a quasi-stationary state, when it is stopped.

Due to the code structure, I can't change the coordinate set (hence, I need to keep z as the direction on which the two fluids are getting out of the computed domain). I was then trying to figure out a way to solve this problem without removing the periodicity assumption, which would require a painful amount of coding. I thought of adding a penalized buffer layer at both z direction extremes on which I could impose the boundary conditions I need.

The problem is I must have Dirichlet BCs for velocities, and homogeneous Neumann BCs for fluids density at z extremes, and this would lead to a constant accumulation of the two fluids in this buffer layer, which will at least give some troubles to reach a quasi-stationary state for the simulation to end. I thought of adding a densities sink inside this buffer layer for this reason but I'm not sure about how this could impact the simulation inside the "real" computed domain (since one of the two fluids is not actually flowing at sonic speed).

Hence my questions are:

how do you think the densities sink inside the buffer layer would influence the simulation inside the "real" computed domain?

do you have any other idea to deal with this problem?

Thanks a lot in advance to anyone who will pay attention to this.


问:
Multigrid: Gauss–Seidel eigenvalues and eigenvectors
I am trying to figure out the questions from Multigrid Tutorial by Briggs. However I am stuck on these two questions.

https://www.researchgate.net/publication/220690328_A_Multigrid_Tutorial_2nd_Edition

(a) Show that the eigenvalue problem for the Gauss-Seidel iteration matrix, $R_{G} \mathbf{w}=\lambda \mathbf{w}$, may be expressed in the form $U \mathbf{w}=(D-L) \lambda I \mathbf{w}$, where
$U, L, D$ are defined in the text.

(c) Show that the eigenvector and eigenvalue are associated with $\lambda_{k}$ = $\cos ^{2}\left(\frac{k \pi}{N+1}\right)$ is $w_{k, j}=\left[\cos \left(\frac{k \pi}{N+1}\right)\right]^j \sin \left(\frac{j k \pi}{N+1}\right)$.

My attempt

For part a, I derived $R_{G}=(D-L)^{-1} U$ where D, L, U is the diagonal matrix of A. I am not sure how to get the answer as shown in question (a).

I am lost to show the eigenvector and eigenvalue relationship for part c.


问:
continuous analogues of Newton's method
Suppose we want to minimize some convex functional $J(u)$ where $u$ lives in some Banach space $V$.
The classical Newton method

$$\mathrm d^2J(u_n)(u_{n + 1} - u_n) = -\mathrm dJ(u_n)$$

can be viewed as an explicit Euler discretization of the dynamical system

$$\tau\;\mathrm d^2J(u)\frac{\mathrm du}{\mathrm dt} = -\mathrm dJ(u),$$

that takes a constant timestep $\tau$ at every step.
You can think of a primitive damped Newton method as taking a timestep that is some fraction of $\tau$, and a globally-convergent Newton line search method as using a variable timestep that is tuned to make sure that the objective decreases.
Several papers have remarked on this (see below).

Are there improved, practical solution methods based on viewing Newton's method as a dynamical system?
For example, what happens if, instead of doing a line search at every iteration, you were to use a classical adaptive Runge Kutta method with stepsize control like RK-1/2 or similar?
More generally, is there some respect in which this dynamical systems viewpoint is practically relevant, or is it just a mathematical curiosity?

References:

https://doi.org/10.1080/01630560802484294

https://doi.org/10.1142/9789812792686_0012

In Russian, but arguably the first: http://www.mathnet.ru/links/fec5f826de87c3778ef0ddea8e20b256/ivm2977.pdf

Title: continuous analogues of Newton's method


问:
Comparison of computational complexities of MD versus MC simulations
In my humble understanding MD simulations of systems with short-range(like LJ interactions) and long-range(electrostatic) has a computational complexity $O(N . log(N))$. What will be the computational complexity of a MC(monte-carlo) for the same system? Also what will be the complexities for both cases if electrostatic computations are neglected?


问:
Do computational scientists think about the statistical aspects of using least squares?
Correct me if I'm wrong, but when least squares is used in the computational science community, it's typically not in the context of regression. It could be used to solve for gradients in discretized PDEs and for other things. In these kind contexts, is there any relevance to consider considering statistical aspects such as checking p values of the least squares coefficients, ensuring that the residuals $\hat{y} - y_{observed}$ are normally distributed, etc..?

I don't think I've ever seen the word "regression" brought up when using least squares in the scientific computing community, and that gave me the impression that computational scientists and statisticians use least squares in near orthogonal manners.


问:
hardware agnostic (GPU/CPU) sparse linear algebra C++ solvers framework/technology
I am looking for recommendations on matured C++ solver for Linear Sparse Algebra problems.
The goal is to select between more or less GPU hardware agnostic libraries/frameworks that can be compiled on both Linux and Windows with quite high control on host-device memory.
My current research narrows down to the following candidates:

CUDA - matured, more than BLAS, cross platform, run on NVIDIA only

OpenCL - matured, competitor to CUDA

HIP - less matured, but picks lot from CUDA, runs on both NVIDIA & AMD too, no cost on porting between OpenCL and CUDA

Intel One API (DPC++) - implements SYCL, seems to be only MKL/LAPACK wrapper now, not sure how it works with GPU system, since MKL is more CPU-centric - seems immature

Magma - matured, looks like big competitor to CUDA, smaller community

I am close to start with CUDA, since it has a lot of examples in the toolkit, and then potentially migrate to HIP (with hipify) or to Intel's DPC++ SYCL guide.

Do you have a link/article/opinion that will help me in challenging these considerations?


问:
Jacobian matrix cutoff in ODE solver
I am studying an implementation of a 3rd semi-implicit Runge Kutta method (siRK3) from the book by Villadsen & Michelsen (1978), Solution of differential equation models by polynomial approximation, Prentice-Hall.

The code expects an analytically evaluated Jacobian matrix $J = \partial F / \partial y$. In the calculation of the new function values $y^{n+1}(t + h)$, the matrix $I -h a J$ is assembled by overwriting the storage used by the Jacobian matrix (array df):

```<code>    do i = 1, n
      do j = 1, n
        df(i,j) = -h * a * df(i,j)
        if (abs(df(i,j)) < DF_TOL) df(i,j) = 0.0_dp
      end do
      df(i,i) = df(i,i) + 1.0_dp
    end do
</code>```

In the original code, DF_TOL = 1.0d-12. Afterwards an LU factorization of the matrix df is performed and used to calculate the new function values by back-substitution.

What is the purpose of truncating the Jacobian matrix values, and why choose this particular cutoff value? Is this type of truncation needed at all?

Title: Jacobian matrix cutoff in ODE solver


问:
Weak form of Elliptic problem with mixed Dirichlet & Neumann conditions
Let $\Omega \subset \mathbb{R}$ in a bounded polygon domain and $f:\Omega \to \mathbb{R}$ known function.We split the boundary into two parts $\partial \Omega_{1}$ and $\partial \Omega_{2}$ such that $\partial \Omega =\partial \Omega_{1} \cup \partial \Omega_{2}.$
The elliptic problem of  Dirichlet & Neumann conditions is :

\begin{align*}
-\Delta u +u &=f \quad in \quad \Omega \subset\mathbb{R}^2 \\
u&=0 \quad in \quad \partial \Omega_{1}   \\
\frac{\partial u}{\partial \mathbf{n}}&=1 \quad in \quad \partial \Omega_{2}
\end{align*}

It was asked to write the problem in a weak form in a suitable solution space $\mathcal{V}$ and the bilinear form $a(\cdot,\cdot):\mathcal{V} \times \mathcal{V} \to \mathbb{R}$ as long as a norm $||\cdot||_{\mathcal{V}}$ and show that is indeed a norm.

My beginning thought :

\begin{align*}
-\Delta u+u &=f \\
-\Delta u v+uv &=fv \\
-\int_{\Omega}\Delta u v dV+ \int_{\Omega}uv dV&=\int_{\Omega}fv dV \\
\int_{\Omega}\nabla u \cdot \nabla v dV- \int_{\partial \Omega_{2}}v \nabla u \cdot \vec{n} dS  + \int_{\Omega}uv dV&=\int_{\Omega}fv dV \\
\int_{\Omega}\nabla u \cdot \nabla v dV- 1 + \int_{\Omega}uv dV&=\int_{\Omega}fv dV \\
\end{align*}

So $$a(u,v) = \int_{\Omega} u'(x) \cdot  v'(x) dx- 1 + \int_{\Omega}u (x)\cdot v(x) dx=\int_{\Omega}f(x) \cdot v(x) dx \quad \text{(1)}$$

in $\mathcal{V}_{h}$ subspace of $\mathcal{V}$

$$a(u,v) = \int_{\Omega} u_{h}'(x) \cdot  v_{h}'(x) dx- 1 + \int_{\Omega}u(x) \cdot v_{h(x)} dx=\int_{\Omega}f(x) \cdot v_{h}(x) dx  \quad \text{(2)} $$

Subtracting (1) from (2) we obtain :
$$\int_{\Omega} \left(u'(x)- u'_{h}(x)\right) v'_{h}dx = 0$$

but $$||u ||_{\mathcal{V}} = \left( \int_{\Omega}(v'(x))^2 dx\right)^{1/2}$$ is indeed a norm for the space $\mathcal{V}$

\begin{align*}
||u - u_{h} ||_{\mathcal{V}} &= \inf_{v_{h} \in \mathcal{V}_{h}} \\
||u - u_{h} ||_{\mathcal{V}}^2&=  \int_{\Omega} \left( u'(x) -v'_{h}(x)  \right)\left( u'(x) -v'_{h}(x)  \right) dx \\
&= \int_{\Omega} \left( u'(x) -v'_{h}(x)  \right) u'(x)   dx -
\underbrace{\int_{\Omega} \left( u'(x) -v'_{h}(x)  \right) u'_{h}(x)   dx}_{G.O = 0 }\\
&= \int_{\Omega} \left( u'(x) -v'_{h}(x)  \right) u'(x)   dx - \underbrace{
\int_{\Omega} \left( u'(x) -v'_{h}(x)  \right) v'_{h}(x)   dx}_{G.O = 0 }
\end{align*}

Therefore $$||u - u_{h} ||_{\mathcal{V}}^2 =\int_{\Omega} \left( u'(x) -v'_{h}(x)  \right)\left( u'(x) -v'_{h}(x)  \right) dx \leq ||u-u_{h}||_{\mathcal{V}}||v-v_{h}||_{\mathcal{V}}  $$

and by the definition of inf, we must also have:

$$\inf_{v_{h} \in \mathcal{V}_{h}}   ||u-u_{h}||_{\mathcal{V}} \leq ||v-v_{h}||_{\mathcal{V}} $$

right calculation (post)

\begin{align*}
\int_{\partial \Omega}v \nabla u \cdot \vec{n} dS &=\int_{\partial\Omega_{1}} v \cdot \frac{\partial u}{\partial n} - \int_{\partial \Omega_{2}} v \cdot \frac{\partial u}{\partial n} dV\\
&=0 + \int_{\partial\Omega_{2}} v \cdot 1 dV \\
&= \int_{\partial\Omega_{2}} v  dV
\end{align*}

$$\int_{ \Omega}\nabla u \cdot \nabla v dV +\int_{\partial \Omega_{2}} v  dV+ \int_{ \Omega}uv dV=\int_{\Omega}fv dV$$


问:
Algorithm for Gaussian random field generation via Karhnunen--Loeve expansion
I have difficulty understanding how to generate a Gaussian random field from Karhunen--Loeve expansion.

Precisely, I need to generate a zero-mean field with covariance function $C(x, y)$. We find eigendecomposition of $C(x, y)$ from the equation:

$$
\int_0^1 C(x, y)\, e(y)\, \mathrm dy = \lambda \,e(x)
$$

I discretize the above equation in the following way. For grid $0=x_0 < \dots x_I=1$ with $I+1$ points and the grid size $1/I$, I use midpoint method to obtain

$$
\frac1I \sum_{j=0}^I C(x_{i+1/2}, y_{j+1/2}) \, e(y_{j+1/2}) = \lambda \, e(x_{j+1/2})
$$
which can be compactly represented as matrix eigendecomposition problem for matrix $\frac1I C$:
$$
\frac1I C e = \lambda e
$$
$C \in \mathbb R^{I\times I}$, $e \in \mathbb R^I$, $\lambda \in \mathbb R$. Solving this problem, we obtain $I$ eigenvalues and eigenvectors, and then approximated random field is
$$
\kappa(x, \omega) \approx \sum_{k=1}^N \sqrt(\lambda_k) y_k(\omega) e_k
$$
where $y_k(\omega)$ are standard Normal Gaussian variables.

The problem that I have is that the generated random field is sensitive to the grid size $I$, and I do not understand where the problem is and how to fix it.


问:
Relative interior requirement in Slater's condition
I'm reading Convex Optimization by Boyd and Vandenberghe. This is how they describe Slater's condition:

What I don't understand is why it is necessary to enforce that $x$ be in the relative interior of $\mathcal{D}$, where

$$
\mathcal{D} := \bigcap_{i = 1}^m \text{dom}(f_i),
$$

i.e., it is the intersection of the domains of the $f_i$'s. Doesn't the fact $f_i(x) < 0$ for all $i$ immediately imply that $x \in \text{relint}(\mathcal{D})$? So why state the theorem this way?

As a side note, Wikipedia gives the exact same seemingly extraneous condition:

Please let me know what I'm misunderstanding. Why is this additional condition necessary? Thank you!

EDIT: I probably have some fundamental misunderstanding, so I'm adding this edit to be more clear about what I'm thinking. If I'm not wrong, the fact that the constraints $f_i$ are convex means that the intersection of the spaces spelled out by $f_i(x) \le 0$ for all $i = 1, \dots, m$ is convex. The feasible region is then the intersection of this convex space with an affine space, that being the affine space spelled out by $Ax = b$.

It seems to me that any point $x$ that satisfies $f_i(x) < 0$ will not only be in $\text{relint} \mathcal{D}$, it will even be in the relative interior of the feasible region described above. This is just because the convex functions $f_i(x) = 0$ form the boundary of the feasible region (in a semi-intuitive sense). The fact that $x$ satisfies $f_i(x) < 0$ for all $i$ means that it will be an interior point since it doesn't lie on the boundary.