问:
How does Tannehill impose boundary conditions when coding the Parabolized Navier Stokes on an Implicit Finite Differences Scheme?
I'm trying to implement the scheme he describes on his book "Computational Fluid Mechanics and Heat Transfer" on Chap.9 and I'm having trouble imposing BC.

I don’t get how he imposes them. I get up to the point where, following the Beam and Warming scheme and after linearizing each flux vector with the flux vector U, he gets to solve a linear system of a tridiagonal matrix with block matrixes inside (A, B,and C he calls them).

I understand the BC must be imposed on the
first and last 5 rows of that matrix (1 for each variable of the flux vector). The problem is, if that linear system is written with the deltaU vector, the BC conditions must be imposed on the deltaU variables, and I don't know how.

I've already tried some approaches on my
code, but none of them seem to work, as the Boundary Conditions I'm imposing always seem to oscillate.


问:
I would like to fit a data to integral type model but cannot figure it how?
I have cumulative data for concentration versus time. I would like to fit the following model to the data:

$$\frac{1}{\Gamma(n) \bar{t}}\int_{0}^{\frac{n t}{\bar{t}}}z^{n-1}  e^{-z}dz$$

$n$ and $\bar{t}$ are the variables that I would like to fit. I tried to make a function with integral included but I wasn't successful would you please give me a hint on how to proceed?

ContentLicense: CC BY-SA 4.0


问:
Estimating $n$ column norms of $A$ by using less than $n$ calls to $Ax$?
Suppose I have an $m \times n$ matrix $A$ that is only accessible through matrix-vector products $Ax$. These are expensive, is there a way to approximate vector $x$ of $n$ column norms of $A$ by using less than $n$ calls to the matrix-vector product?

$$
\begin{aligned}a_{1}=\left[ A\right] \begin{pmatrix}
1 \\
0 \\
0
\end{pmatrix}\\a_{2}=\left[ A\right] \begin{pmatrix}
0 \\
1 \\
0
\end{pmatrix}\\a_{3}=\left[ A\right] \begin{pmatrix}
0 \\
0 \\
1
\end{pmatrix}\\\\
x=\langle\left\| a_{1}\right\|,\left\| a_{2}\right\|,\left\| a_{3}\right\| \rangle\end{aligned}
$$

Note that I can estimate a vector of row norms by sampling a few random vectors $y$ from standard normal, vector of row norms squared is proportional to entries of $E[(Ay) \odot (Ay)]$

Edit random idea, if we let $y_{1,2}^T=\langle y_1, y_2\rangle$ with $E[y_{1,2} y_{1,2}^T]=I$ and $\hat{y}_{1,2}^T=\langle y_1, y_2, 0, 0, \ldots\rangle$ then

$$ E[\|A\hat{y}\|^2]=\DeclareMathOperator{\Tr}{Tr} E \Tr [(A\hat{y})(A\hat{y})^T]=\|a_1\|^2+\|a_2\|^2=t_{1,2}$$

Perhaps there's a strategy of choosing a good collection of overlapping subsets $\{s_1,s_2,\ldots\}$ such that $\{\|a_1\|^2,\|a_2\|^2,\ldots\}$ could be recovered from corresponding totals $\{t_{s_1},t_{s_2},\ldots\}$


问:
Problem with Time Evolving Block Decimation for small finite chains
Preliminary definitions:

I have a code that uses TEBD (https://en.wikipedia.org/wiki/Time-evolving_block_decimation) to perform time evolution of a given Matrix Product State (MPS), lets call it $|\Psi\rangle$. I am starting from a product state (Neél state) $|\uparrow,\downarrow,\uparrow,..,\downarrow\rangle$ for a chain of $L=10$ spins with spin $S=1/2$. I work with Open Boundary Conditions (OBC). The state is represented in the standard $\Gamma-\Lambda$ notation:
\begin{eqnarray}
|\Psi\rangle=\sum_{\sigma_{1},...,\sigma_{L}}\Lambda_{0}\Gamma_{1}\Lambda_{1}\Gamma_{2}\Lambda_{2}...\Lambda_{L-1}\Gamma_{L}\Lambda_{L}|\sigma_{1},...,\sigma_{L}\rangle,
\end{eqnarray}
where $\Lambda_0=\Lambda_L=1$, and since the Neél state is a low-entangled state, initially all $\Lambda_{n=1,...,L-1}\equiv\text{diag}(1,0)$, i.e. I am fixing the initial bond dimension of ALL bonds to $\chi=2$. Therefore, initially, all $\Gamma,\Lambda$ are composed by $2\times 2$ matrices, and no loss of information on that initial state occurs.

The state is then time evolved (quenched) using the Heisenberg Hamiltonian:
\begin{eqnarray}
H=\sum_{i=1}^{L-1}\left(S_{i}^{x}S_{i+1}^{x}+S_{i}^{y}S_{i+1}^{y}\right) + \Delta\sum_{i=1}^{L-1}S_{i}^{z}S_{i+1}^{z}=\sum_{i=1}^{L-1}h_{i,i+1}.
\end{eqnarray}
The unitary evolution operator is decomposed using a 2nd order Trotter decomposition, acting on even and odd bonds to update the state after a single time step $dt$.

Problem:

As stated above, we start from a very low entangled state. After application of the unitary operator, and several $dt$ steps, I observe a growth of entanglement as expected. In particular, the $\Lambda_{n}$ matrices increase their dimension from $\chi=2$ to a dimension dependent on the bond position, i.e.:
\begin{eqnarray}
\text{dim}(\Lambda_{n})=d^{n},\hspace{10pt}\forall n\in[1,L-1].
\end{eqnarray}
However, after all dimensions reach these values, I observe no more growth in the $\Lambda_{n}$ dimensions, with a maximum dimension of $\chi=32$ at the middle bond ($\Lambda_{5}$), but fixed $\text{dim}(\Lambda_{1})=\text{dim}(\Lambda_{L-1})=2$, etc...

The problem is that way after this happens, I have no way to measure the discarded weight after a single SVD decomposition, as this always gives me 0, i.e., no eigenvalues of the $\Lambda$ matrices are discarded anymore. My feeling is that this is a wrong behaviour, but on the other hand, for a system of $L=10$ spin 1/2, the total amount of coefficients $c_{\sigma_{1},...,\sigma_{L}}$ is bounded to a total of $2^{L}=1024$ for $L=10$.

Is this normal in small systems approached by TEBD? Or on the contrary, should the bond dimensions of all $\Lambda_{n}$ keep growing and growing, in accordance with the entanglement growth featured by the quench? I am confused about this; my observables seem to reproduce some benchmarking results and I have observed that when $L$ is increased, I can reproduce the thermodynamic limit results. But why is this failing for the small chain case? How can we record the discarded weight after each dt step?


问:
Choice between using Moore-Penrose inverse and G2 inverse
Moore-Penrose inverse for an arbitrary matrix $X\in \mathbb{R}^{n \times p}$ is defined by a matrix $X^\dagger$ satisfying all of the Moore-Penrose conditions, namely

\begin{align}
(1) \;\;\;& XX^\dagger X=X \\
(2) \;\;\;& X^\dagger X X^\dagger = X^\dagger \\
(3) \;\;\;& XX^\dagger = (XX^\dagger)^T \\
(4) \;\;\;& X^\dagger X=(X^\dagger X)^T
\end{align}

which can be found by singular value decomposition (SVD). If a matrix $G$ satisfies only the first condition $(1)$, then it is called a generalized inverse (G-inverse). In addition, if it satisfies both two conditions $(1)$ and $(2)$, then it is called a G2-inverse.

My question is rather general. Since there is always a unique Moore-Penrose inverse for $X$ by SVD, why and when should we also consider G or G2-inverse? Is there any computational advantages of using G2-inverse instead of Moore-Penrose inverse?

There is an example of using G2 inverse in ANOVA e.g. https://www.stat.purdue.edu/~fmliang/ANOVA.pdf

Further information, G2SWEEP algorithm is suggested in SAS Technical Report R106 (https://support.sas.com/documentation/onlinedoc/v82/techreport_r106.pdf). However, I don't see how G2-inverse is better than Moore-Penrose inverse.


问:
Calculating the species mass consumption from implicit reaction-term in diffusion-reaction equation
The 1D diffusion equation with a chemical source term has the following form:
$$\frac{\partial Y}{\partial t} = D \frac{\partial^2 Y}{\partial x^2} - k Y,$$
where $Y$ is the molar concentration of the diffusing species (mol/m³), $D$ is the diffusion coefficient (m²/s), and $k$ is the reaction rate coefficient (1/s).

The boundary conditions are:
$$\frac{\partial Y}{\partial x}|_{x=0} = 0, $$
$$-D \frac{\partial Y}{\partial x}|_{x=xN} = mtc (Y_c - Y), $$
where $Y_c$ is the species concentration outside the simulation domain (mol/m³), and $mtc$ is the mass transfer coefficient (m/s).

I can solve this equation numerically stable using an implicit approach:
$$Y(t=t_0) = A\cdot Y(t=t_0+dt),$$
where the matrix $A$ has a term that depends on $k$ in the diagonal only.

My Problem is the following:
I need to know how much species mass was converted by the chemical reactions at each grid-point, so I can calculate the heat generated/consumed at each grid-point as well as know how much of an additional stationary solid phase was converted in the process.
As I solve the equation implicitly I do not know how to split the diffusive contribution and the reactive contribution to the final $Y$ at each grid-point. It is further complicated by mass not being conserved in the system due to the flux at the right BC.

Without reactions I can always calculate the difference of $Y$ before and after a time step and I know how much mass was transported in total to or from this grid-point in the given time.
If I add the reaction term then this difference is a combination of mass being converted due to the reaction and mass being transported by diffusion.

I have spent some time researching reaction-diffusion equations and their solutions, however, I have not been able to find any information about how to calculate how much of my species was converted chemically at each grid-point.

The only workaround I have found so far was to explicitly calculate a chemical source term using the initial value of $Y(t=t_0)$ at a given time-step and then including it as an explicit source in the above equation. The problem with this approach was that this leads to me having to severely limit the size of the time-step such that only a fraction of the available mass at a grid-point is converted in each time step. However, this leads to unsustainably small time-steps for some examples I encountered.

I would be grateful for any help regarding this issue. I doubt there is an exact way to split the contributions of diffusion and reaction and am happy for workarounds or approximations of all sort.


问:
A Finite Element Method for a first order PDE?
I want to develop a finite element method to solve for $u(x,t)$ the PDE:
$$u_t+c u_x= \frac{-c}{x}u$$
where $c$ is a constant.

so I am trying the following ( as Rothe's method? ) :

Letting $k= t_n- t_{n-1}$ ( we begin by time discretization? ) we can write
$$\frac{u_i^n-u_i^{n-1}}{k}+ c\left(\frac{u_i^n-u_{i-1}^{n}}{\Delta x}\right)=\frac{-c}{x}u_i^n$$
Which gives the time stepping formula $$\left(1+ \frac{kc}{x}+ \frac{kc}{\Delta x}\right)u_i^n= u_i^{n-1}+ \frac{kc}{\Delta x}u_{i-1}^{n}$$
I can't avoid to notice that this is a direct finite difference scheme, which can be implemented without using finite elements.

How should I try the weak formulation of the PDE? I need help with that
$$\int_\Omega u_t\phi(x) dx +c \int_\Omega u_x\phi(x) dx= \int_\Omega \frac{-c}{x}u\phi(x) dx$$
$$\int_\Omega u_t\phi(x) dx +cu\phi(x)- c\int_\Omega u\phi'(x) dx= \int_\Omega \frac{-c}{x}u\phi(x) dx$$
How to handle the time derivative of $u$ ?


问:
Splitting system of equations into linear and nonlinear part and solving separately
I was working on a problem recently (calculating all flows in a network given input and output flows, basically what Hardy-Cross tries to achieve) which can be formulated as a well-determined system of $N$ equations (with $N$ variables), of which $M$ equations are linear and $K$ are non-linear (so $M + K = N$, and generally $K < M$).

The easiest way to solve this was to just use a nonlinear solver and put in all $N$ equations. However, taking inspiration from the Hardy-Cross method I decided to instead split the problem into a linear part and a non-linear part, then first solve the linear part directly and proceed to solve the $K$-dimensional non-linear part iteratively. More concretely, the steps I took were as follows:

Suppose the linear part is $A x=b$. Find a solution $x_{lin}$ that satisfies $A x_{lin} = b$ using least squares ($A$ is underdetermined so this is possible).

The space of solutions to the linear part can be written as $\{v \mid A (v - x_{lin}) = 0 \}$. This space will be $K$-dimensional. Build the space of solutions by constructing a basis of the nullspace of $A$. Store them as the columns of a matrix $A_{null}$.

Suppose the nonlinear part is given by a function $f: \mathbb{R}^N \to \mathbb{R}^K$ and the objective $f(x) = 0$. To solve this, we'll look for a solution in the $K$-dimensional space we constructed before. We can transform a given 'reduced' solution $z \in \mathbb{R}^K$ into a 'full' solution $x \in \mathbb{R}^N$ by doing $x = A_{null} z + x_{lin}$, and the reverse is $z = A_{null}^* (x - x_{lin})$ where $^*$ indicates the (left) psuedo-inverse. We can therefore transform our function $f$ into a 'reduced' variant $g: \mathbb{R}^K \to \mathbb{R}^K$ where $g(z) = f(A_{null} z + x_{lin})$.

Solve $g(z) = 0$ using any nonlinear solver. The solution $z_{sol}$ can be transformed into a solution for the full system of equations as $x_{sol} = A_{null} z_{sol} + x_{lin}$.

This method seems to be quite a bit faster than just solving the whole system of equations using the same nonlinear solver, and to give the same results (in my situation). Now, I don't think I'm the first one ever to have thought of this, so my question is: is this (similar to) an existing method or principle, and what is it called? It would be nice to base my application of this on something a bit more rigorous. I haven't been able to find anything comparable, but I might simply not have the right context or search terms ;)


问:
Numerical methods for $u_t+c u_x= \frac{-c}{x}u$?
I am looking for possible numerical methods to solve the PDE
$$u_t+c u_x= \frac{-c}{x}u$$

I am particularly interested in a Finite elements method, although I am also curious if you can expose some other methods.


问:
From 3D to 2D with a STL file
I would like to do a 2D projection from a 3D geometry saved in a stl file and know the distance between the two projected planes. In order to explain better the concept I will start with an almost ideal model or simple example:

Let's say I have a minimal cube model such as the vertex are comprissed between x=-+1, y=-+1 and z=-+1. Therefore, I have 8 vertex, and I have 12 faces (triangles). If I want to get the projection of the cube at the middle point of z. I will end up with 2 triangles and 4 vertex, being x=-+1, y=-+1 and z=0. I also will now that the distance between the z planes of the projection are 2 in all the vertexs (whatever the units).

Now the real problem:

I have a 3D fracture of a surface (in .stl format). The z planes are in this case curve surfaces (more or less parallel walls) and therefore they are not in an exact position. Something similar as what I am trying to convey can be found here https://www.mdpi.com/1996-1073/11/6/1409/htm (figure 6 b, or 6d or 6f).

I am wondering if somebody has an idea of how I can do a 2D plane from the 3D plane (in stl format) and know the distance between the two irregular curve surface walls. Any software (if free better), or library for coding such thing could be helpful.

Thanks

Best regards,

Dani

Title: From 3D to 2D with a STL file


问:
Dense least-squares with millions of variables
Suppose $X$ is a dense $m\times n$ data matrix and we seek to find $w$  by approximately iterating least-squares filter equation:

$$w = w - \mu (X'X)^{-1}X'(Xw-b)$$

What are known approaches for $10^9<m<\infty$ and $n\approx 10^6$?

This setting comes up in Machine Learning and people tend to invent their own techniques for approximate pseudo-inverse of $X$, which in practice has low stable-rank, without referencing scientific computing literature.

For instance KFAC method addresses large $n$ by writing $x\approx a\otimes b$ and using this factorization to approximate $(X'X)^{-1}\approx (A'A)^{-1}\otimes (B'B)^{-1}$. Making this update numerically stable requires considerable tweaking.

Any literature pointers appreciated!

Edit Jul 19 crossposted on mathoverflow


问:
Numerical integration of a 2D hemisphere discrete dataset where limits are unknown
I am trying to compute the integral of a 2D hemisphere dataset $f_r \, (\theta, \phi)$ where $\theta \in [0, \pi / 2[$ and $\phi \in [0, 2\pi]$. I am making the measurements myself, so I can choose the number of values and their coordinates.

My first idea was to use the different points to triangulate my definition set, and sum the mean value of the vertices of each triangle divided by their surface.

But my problem will appear when $\theta \rightarrow \pi/2$, because I cannot have any measurement for these values. Is there a good strategy to solve this problem ?


问:
Derivation of the second Piola-Kirchhoff tensor
I tried many formulas to find the components of the second Piola Kirchhoff.  I need help to derive equations 27-29 on reference 1.

[![enter image description here][1]][1]

we have
$$
\mathbf{C}=C_{11} \mathbf{T} \otimes \mathbf{T}+C_{12}(\mathbf{T} \otimes \mathbf{N}+\mathbf{N} \otimes \mathbf{T})+C_{22} \mathbf{N} \otimes \mathbf{N}+C_{33} \mathbf{E}_{3} \otimes \mathbf{E}_{3} \tag{16}
$$
Here, the components $C_{11}$ and $C_{12}$ are expressed in terms of the kinematic quantities of deformation, namely $\left\{u_{1}, u_{2}, \psi, \theta, \kappa_{0}, \mathbf{T}, \mathbf{N}\right\}$, in the following form:
$$
C_{11} \approx \mathbf{t} \cdot \mathbf{t}+2 z\left(\mathbf{t} \cdot \mathbf{b}+\kappa_{0} \mathbf{t} \cdot \mathbf{t}\right), C_{12} \approx \Gamma=\mathbf{t} \cdot \mathbf{n} \tag{17}
$$
we have
.\begin{array}{r}
I_{1}=\operatorname{tr} \mathrm{C}=\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2} \\
I_{2}=\frac{1}{2}\left(I_{1}^{2}-\operatorname{tr} \mathbf{C}^{2}\right)=\lambda_{1}^{2} \lambda_{2}^{2}+\lambda_{2}^{2} \lambda_{3}^{2}+\lambda_{3}^{2} \lambda_{1}^{2} \\
I_{3}=\operatorname{det} \mathbf{C}=\left(\lambda_{1} \lambda_{2} \lambda_{3}\right)^{2} \tag{20}
\end{array}

consider the neo-Hookean material model with the strain energy density function
$$
W=\frac{\mu}{2}\left(I_{1}-3\right),
$$
where $\mu$ is the shear modulus. From Eq. (23), it follows that
$$
\mathbf{S}=\mu \mathbf{I}-P \mathbf{C}^{-1}
$$
With the aid of Eqs. (16), $(20)_{3}$ and (26), the plane stress assumptions $S_{22}=S_{33}=0$ and the incompressibility constraint $I_{3}=1$, one will have
$$
\left.\begin{array}{r}
S_{22}=\mu-\left(C_{11} C_{22}-\Gamma^{2}\right)^{-1} C_{11} P=0 \\
S_{33}=\mu-C_{33}^{-1} P=0 \\
I_{3}=\left(C_{41} C_{22}-\Gamma^{2}\right) C_{33}=1
\end{array}\right\} \tag{27}
$$
This provides a system of three nonlinear algebraic equations in $\left\{C_{22}, C_{33}, P\right\}$. The exact solution of Eqs. (27) is given by
$$
C_{22}=\frac{\Gamma^{2}+\sqrt{C_{11}}}{C_{11}}, C_{33}=\frac{1}{\sqrt{C_{11}}}, P=\frac{\mu}{\sqrt{C_{11}}} \tag{28}
$$
It is worthwhile to note that the solution obtained for $C_{33}$ satisfies the condition $C_{33}=\lambda_{3}^{2}>0$. Now, the nonzero components of the second Piola-Kirchhoff stress $S$ in the Frenet frame become
$$
S_{11}=\mu \frac{C_{11}^{2}-\Gamma^{2}-\sqrt{C_{11}}}{C_{11}^{2}}, \quad S_{12}=\mu \frac{\Gamma}{C_{11}} \tag{29}
$$
[1]: https://i.stack.imgur.com/CwsMm.png

References

Dadgar-Rad, F., & Sahraee, S. (2021). Large deformation analysis of fully incompressible hyperelastic curved beams. Applied Mathematical Modelling, 93, 89-100.


问:
Notation for a time series with discrete time points
I want to mathematically describe a vector containing a series of time points. The time points must lie on an equidistant time grid of 15 minutes of a day. Thus, elements of the vector can contain values that occur at times $[t_1 = $ 00:00$, t_2 = $ 00:15$, ..., t_n = $23:45$]$.

Example: vector $\pmb{t} = [$00:00, 01:15, 02:00, 04:15, 15:45, 18:00, 20:30$]^T$

Now I want to describe the vector with $\pmb{t} = [t_1, \dots, t_m] \in \square$. However, in my case $\mathbb{R}$ does not fit as a set. How would you describe the set of $\pmb{t}$?


问:
Numerical solution to the Tolman-Oppenheimer-Volkoff equations for any equation of state (numerical or analytical)
I've been working on a code to solve the Tolman-Oppenheimer-Volkoff (TOV) equations for a while and recently I've got it right but only for one specific equation of state, the bag model, which is not what I need since it should run for any hadronic equation of state. The code has given only wrong results for any EoS other than the bag model and I've tried everything I could think of but nothing has made it so that the final diagrams for other EoS are correct. I'll put the code and some of the wrong results for some of the EoS I've been studying below.

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

pi = math.pi
km2_gevfm3 = (1/(1.3234*10**(-6)))*0.001
gevfm3_km = 1/(km2_gevfm3)
B = 7.424*10**(-5)
cm = 1/1.4760

eos = np.genfromtxt("C:\\Users\\felip\\Desktop\\TOV\\Equações de estado\\PtableTM1.dat")
eden, p = (gevfm3_km*eos[:,0])/1000, (gevfm3_km*eos[:,1])/1000
press = interp1d(eden, p, kind = 'cubic', bounds_error = False, fill_value = 'extrapolate')
# press = UnivariateSpline(eden, p, k = 5, s = 5)
dpde = press._spline.derivative(nu = 1)
# dpde = press.derivative(n = 1)

def TOV(y, r):
     eden, M = float(y[0]), float(y[1])
     dedr = -(1*(dpde(eden))**(-1))*((M*eden)/(r**2))*(1+press(eden)/eden)*(1+(4*pi*press(eden)*r**(3))/(M))*((1-((2*M)/r))**(-1))
     dMdr = 4*pi*eden*(r**2)
     return [dedr, dMdr]

r0 = 10**(-10)
rm = 5000  #km
N = 1000
rn = np.linspace(r0, rm, N)

emin, emax =   7.7*10e-15, 0.0003
h_e = 1000
e_range = np.linspace(emin, emax, h_e)
R,M,E =[],[],[]

for e0 in e_range:
    M0 = (4*pi*e0*r0**(3))/3
    contorno = [e0, M0]
    soltov = odeint(TOV, contorno, rn, atol = 10e-10, rtol = 10e-14, mxstep=5000)
    rmax = r0
    mmax = M0
    esol = soltov[:,0]
    psol = press(esol)
    k=0
    for i in psol[:-1]:
        h=rn[k+1]-rn[k] + r0
        rmax = rmax + h
        mmax = soltov[k,1]
        k=k+1
        if i <= 0:
                break
    R.append(rmax)
    M.append(mmax*cm)
    E.append(e0*km2_gevfm3)

immax = np.argmax(M)
Mmax = M[immax]
Rmax = R[immax]

f1 = plt.figure()
ax1 = f1.add_subplot(111)
ax1.title.set_text('M x E')
ax1.set_xlabel('<span class="math-container">$(\epsilon)_c [GeV/fm^3]$</span>')
ax1.set_ylabel('<span class="math-container">$M [M_\odot]$</span>')
ax1.plot(E,M,'-',color='blue')
ax1.legend()

f2 = plt.figure()
ax2 = f2.add_subplot(111)
ax2.title.set_text('M x R')
ax2.set_xlabel('<span class="math-container">$R [Km]$</span>')
ax2.set_ylabel('<span class="math-container">$M [M_\odot]$</span>')
ax2.plot(R,M,'-',color='red')
ax2.legend()
plt.show()
</code>```

I couldn't get the ident right but it is correct on my code. I have no idea what could be changed here and I've gone through every line in this code and can't figure out why it shouldn't work for other EoS.

For the BPS EoS, the results are : .


问:
How to define a stretched coordinate perfectly matched layer (PML) parameters for maximum absorbtion?
I have written a MATLAB code for solving Maxwell's equations (electromagnetic wave propagation) in 3D with perfectly matched layer (PML) boundaries. I am using a stretched coordinate PML, but I see reflections from the boundary and it can be tuned by varying the parameters involved like sigma, kappa, and alpha. Different papers use different values for these parameters. Is there a standard way by which the values can be determined for a wave scattering problem with maximum absorption?</p>


问:
Deconvolution of sinc function in spectrum calculation in FTS
In Fourier transform spectroscopy (FTS) I am calculating a broadband interferogram (e.m. frequency 190-300 GHz top-hat), then back-retrieving the spectrum by FT.
Here in the figure, you can see the result: I expect the retrieved spectrum to be the same as the "injected" one but there is the effect of the top-hat bandpass.

How do I correct the top-hat bandpass in python?
I suspect I need to deconvolve the bandpass function with the spectrum (or the interferogram). As shown in the following figure example, the FT of the sinc function (the result of the bandpass gives the same effect).


问:
Solving differential equations with fast oscillations using odeint
I have wrote this code to solve an equation , I know the behavior of this function has very rapid oscillations, when I RUN it gives bogus values for some "m[x]" and some "t"'s, with this error:

C:\Users\dani\anaconda3\lib\site-packages\scipy\integrate\odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning)

```<code>import scipy as sio
import numpy as np
import mpmath as mp
import scipy.integrate as spi
import matplotlib.pyplot as plt
import time
initial_value=np.logspace(24,27,100)
t=np.logspace(0,6,100)
m=np.logspace(0,6,100)
start_time=time.perf_counter()
phi_m={}
phi_m_prime={}
phi=[]
phi_prime=[]
j=0
i=np.pi*2.435*initial_value[0]
while i<(np.pi*(2.435*10**(27))):
    i=np.pi*2.435*initial_value[j]
    phi=[]
    phi_prime=[]
    for x in range (len(m)):
        def dzdt(z,T):
            return [z[1], -3*1.4441*(10**(-6))*m[x]*np.sqrt(0.69)*(mp.coth(1.5*np.sqrt(0.69)*        (10**(-6))*1.4441*m[x]*T))*z[1] - z[0]]
        z0 = [i,0]
        ts = tp/m[x]
        zs = spi.odeint(dzdt, z0, ts)
        phi.append(zs[99,0])
        phi_prime.append(zs[99,1])
    phi_m[j]=phi
    phi_m_prime[j]=phi_prime
    j+=1
end_time=time.perf_counter()
print(end_time-start_time,"seconds")
</code>```

I don't know what is the problem. how can I get correct results? or at least as accurate as possible? or maybe I should rewrite the code in another form? thank you.

UPDATE_1:
I increased the steps and decreased the step sizes in this way:

```<code>tp=[]
t1=np.logspace(0,1,100)
t2=np.logspace(1,3,500)
t3=np.logspace(3,4,700)
t4=np.logspace(4,5,800)
t5=np.logspace(5,6,1000)
</code>```

and used tp.append() for every single elements in t1,t2,...,t5 ; but there is still some false results for some "m[x]"'s like for x=8 to 13, surprisingly results for x<8 and x>13 are not too bad!

UPDATE_2:

once again I increased the steps number:

```<code>tp=[]
t1=np.logspace(0,1,100)
t2=np.logspace(1,3,5000)
t3=np.logspace(3,4,7000)
t4=np.logspace(4,5,8000)
t5=np.logspace(5,6,10000)
</code>```

and used tp.append() for the elements; so now I have "tp" with 20100 steps that has more steps and smaller step sizes in the parts that function has his rapid oscillations; after several hours it's still under running!
I do not know if this method will help or not.


问:
Dense factorization specialized for RBF-FD method
In RBF-FD methods (see Fornberg & Flyer. A Primer on Radial Basis Functions with Application to the Geosciences. SIAM, 2015. Chapter 5.), the finite-difference stencil coefficients for a set of nodes $\{\vec{x}_1,\vec{x}_2,\dots,\vec{x}_n\}$ are found by solving a linear system which involves the block-structured and symmetric matrix

$$
\begin{bmatrix}
A & P \\
P^T & 0
\end{bmatrix},
$$

where $A$ is a matrix involving radial basis functions, e.g. $\phi(r) = r^3$,

$$
A=
\begin{bmatrix}
\phi(|\vec{x}_1 - \vec{x}_1|) & \phi(|\vec{x}_2 - \vec{x}_1|) & \dots & \phi(|\vec{x}_n - \vec{x}_1|) \\
\phi(|\vec{x}_1 - \vec{x}_2|) & \phi(|\vec{x}_2 - \vec{x}_2|) & \dots & \phi(|\vec{x}_n - \vec{x}_2|) \\
\vdots & \vdots &  & \vdots \\
\phi(|\vec{x}_1 - \vec{x}_n|) & \phi(|\vec{x}_2 - \vec{x}_n|) & \dots & \phi(|\vec{x}_n - \vec{x}_n|)
\end{bmatrix}
$$

and $P$ is a Vandermonde matrix for a selected (multi-variate) polynomial basis,

$$
P=
\begin{bmatrix}
1 & x_1 & y_1 & x_1^2 & x_1 y_1 & y_1^2 \\
1 & x_2 & y_2 & x_2^2 & x_2 y_2 & y_2^2 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
1 & x_n & y_n & x_n^2 & x_n y_n & y_n^2 \\
\end{bmatrix}.
$$

An example code for forming such matrices and RBF-FD discretizations is given in
the work of Flyer, Barnett, & Wicker (2016). J. Comput. Phys., 316, 39-62. It is not difficult to write the same thing in Fortran, C, or C++ with calls to LAPACK or BLIS
for the dense matrix factorization. Typical values for $n$ are between 13 and 100, depending also on the desired polynomial order. In practice, the polynomial part is also evaluated using relative coordinates, with the stencil center $\vec{x}_1 =(x_1,y_1)$ taken as the origin.

In a time-dependent problem with stencils that are changing in time, we might need to
recompute the stencils coefficients at each time iteration. Hence it is of interest
to have a fast and efficient decomposition method. This brings me to the following list of questions:

Would a $LDL^T$ factorization with symmetric packed-storage offer any potential benefits?

What about an explicit block-decomposition?

Given the known structure of the matrix with the RBF part involving $\phi$ and the polynomial part $P$ are there any other performance tricks one could exploit?

Is pivoting needed, or could it be avoided with a renumbering of the stencil nodes (for example in reverse order with respect to distance from the stencil center)?


问:
Finite difference on matrices
This question relates closely to another question I asked earlier this week.

Let's say that I want to find a set of matrices $S_i$, with $0\leq i\leq N$, that minimizes some objective function with numerous terms. I can easily compute the gradient of the objective function, with regard to each matrix. This provides me a system of equations where each matrix is linked to its immediate neighbours.

For example, in my case, we can write $S_i=F_i(S_{i-1}, S_{i+1})$, where $F_i$ is some function. We also have some conditions at the boundaries so $S_0$ and $S_N$ are fully determined.

In some sense, this is very similar to a 1D finite difference problem, except that the $S_i$ are matrices, not humble scalars.

For scalars, the problem is normally simple enough to write, for example as a linear problem with one big (and eventually sparse) matrix. There are many ways to solve this type of problem. However, for matrices, it can get quite challenging as matrices do not commute in general so the function $F$ can be more complex to linearize and write as a large matrix, especially if it also involves the inverse of some matrices $S_i$.

Is there any resource on how to solve these kind of problems? Or at least to reformulate it in a way that is more consistent with how people look at this type of problems? For now the best I can do is to solve each matrix from its neighbours over and over, so it takes a long time to eventually converge.

Title: Finite difference on matrices


问:
How does one proceed to solve (big) underdetermined or overdetermined systems of linear equations "nowadays"?
In my numerical linear algebra class we mentioned this problem briefly and according to some other lectures on the internet especially in data driven environments one mostly has to deal with such over- or underdetermined systems. Nevertheless in our lecture we mostly deal with direct methods (mostly matrix decompositions) and iterative methods (CG, GMRES, MINRES, ...) for different types of square! matrices.

Are there similar methods for under- or overdetermined systems?


问:
Computing the overlapping area of curves
I have a set of datapoints that correspond to 5 curves.
Sample dataset and MATLAB inputs:

```<code>scale =  1.5;
x1 = [0,4,6,10,15,20]*scale;
y1 =  [18,17.5,13,12,8,10];
x2 = [0,10.5,28]*scale;
y2= [18.2,10.6,10.3];
x3 = [0,4,6,10,15,20]*scale;
y3 = [18,13,15,12,11,9.6];
x4 = [9,17,28]*scale;
y4 = [5,5.5,7];
x5 = [1,10,20]*scale;
y5 = [3,0.8,2];

plot(x1,y1, '*-', x2, y2,  '*-', x3, y3,  '*-', x4, y4,  '*-', x5, y5,  '*-')
</code>```

I want to find the overlapping area under these curves.

Suggestions on the functions that could be used to find the intersecting area and the x-limits of the intersecting area will be really helpful.


问:
I need MATLAB or Mathematica code to solve this integral over limit 0- λ
I want the code to integrate equation(1) or (2) over the limits using mathematical or mathlab to get equation (3) as the answer of $Z$ vibrational partition function, giving the following additional information $$
\begin{aligned}
Z_{v i b}(\beta) &=\sum_{n=0}^{\lambda} e^{-\beta E_{n}}, \quad \beta=\frac{1}{k T}, \\
\lambda &=-c+\sqrt{A} \pm \sqrt{A-B}
\end{aligned}
$$, $$
Z_{v i b}(\beta)=\sum_{n=0}^{\lambda} e^{\frac{\beta a^{2} \hbar^{2}}{4 m}(B-2 A)+\frac{\beta a^{2} \hbar^{2} B^{2}}{8 m(n+c)^{2}}+\frac{\beta a^{2} \hbar^{2}}{8 m}(n+c)^{2}}
$$ In the classical limit, the summation turns to integral
$$
\begin{array}{l}
Z_{v i b}(\beta)=\int_{0}^{\lambda} e^{\frac{\beta a^{2} \hbar^{2}}{4 m}(B-2 A)+\frac{\beta a^{2} \hbar^{2} B^{2}}{8 m(n+c)^{2}}+\frac{\beta a^{2} \hbar^{2}}{8 m}(n+c)^{2}} d n,............(1) \\
=\int_{c}^{\lambda+c} e^{\frac{\beta a^{2} \hbar^{2}}{4 m }(B-2 A)+\frac{\beta a^{2} \hbar^{2} B^{2}}{8 m \rho^{2}}+\frac{\beta a^{2} \hbar^{2}}{8 m \prime} \rho^{2}} d \rho, \rho=c+n ...........(2)\\
=e^{\beta \Lambda_{3}} \sqrt{\pi}\left[\frac{\left[1-\operatorname{Erf}\left(\Lambda_{1} \sqrt{\beta}\right)\right]-e^{\Lambda_{4} \beta}\left[1-\operatorname{Erf}\left(\Lambda_{2} \sqrt{\beta}\right)\right]}{\Lambda_{5} \sqrt{\beta}}\right]...........(3)
\end{array}
$$

where we have also introduced the following parameters for mathematical simplicity,
$$
\begin{array}{l}
\Lambda_{1}=\frac{1}{\rho} \sqrt{-\frac{a^{2} \hbar^{2} B^{2}}{8 m}}-\rho \sqrt{-\frac{a^{2} \hbar^{2}}{8 m}}, \quad \Lambda_{2}=\frac{1}{\rho} \sqrt{-\frac{a^{2} \hbar^{2} B^{2}}{8 m}}+\rho \sqrt{-\frac{a^{2} \hbar^{2}}{8 m}} \\
\Lambda_{3}=-\frac{a^{2} \hbar^{2} A}{2 m}, \Lambda_{4}=\frac{a^{2} \hbar^{2} B}{2 m}, \Lambda_{5}=\sqrt{-2 a^{2} \hbar^{2}}, \quad c \leq \rho \leq c+\lambda
\end{array}
$$

OwnerDisplayName: tolulope ojuola


问:
Upper bound for conditional relative entropy
I am wondering if the following inequality could hold up to some constants:
$$\mathbb E_X[D_{KL}(P_{Y|X}, Q_{Y|X})]\leq c D_{KL}(P_{X, Y}, Q_{X,Y}),$$
where $c>0$ is some constant, $P$ and $Q$ are two probability measures.

The difficulty now is that KL does not have triangle inequality. If anyone knows whether this is true, I appreciate it a lot!


问:
Is there a permutation used in sparse QR factorizations that better locates small elements on the diagonal of R?
Is there a permutation used in sparse QR factorizations that better locates small elements on the diagonal of R to the end of the diagonal?  As an example, consider the following snippet of MATLAB code:

```<code>A = [1  1  1; ...
    -1 -1  1; ...
     0  0 -1; ...
     0  0  1; ...
     -1 0  0; ...
     0 -1  0];
s = [1e-8;1e-8;1e-8;1;1e-8;1];
B = sparse([A diag(s)]);
[Q,R,p] = qr(sparse(B)');
pp = [5,6,1,2,4,3];
[QQ,RR] = qr(sparse(B(pp,:)'));
</code>```

This matrix is similar to what may occur in an optimization problem with small slack variables.  If we look at the absolute value of the diagonal elements of R, we see a small element in position 5:

I would prefer this element to be at the end of the diagonal of R.  Using the permutation pp, this can be achieved and the the diagonal of RR is seen in this image:

In this case, the sparsity between R and RR is the same.  Ideally, I'd like a diagonal sorted by absolute value similar to what is achieved in a dense QR factorization.  In that case, it's easy to get a rough idea of the effective rank of the original matrix and truncate the matrix R effectively for the applications that I care about.  At the same time, I understand that this would dramatically impact the sparsity of R and for large matrices this is not tractable.  In truth, I don't really need a fully sorted diagonal, but some kind of permutation that mostly puts the small elements at the end.  At the moment, I have some very large matrices and may get some very small elements in the first few dozen diagonal elements, which is undesirable.

As some additional information, if it helps, I understand that MATLAB uses spqr as its underlying QR factorization algorithm and that the parameter spqrtol can be set in spparms to adjust the drop tolerance for small elements.  Alternatively, spqr could just be used directly.  In this case, I'd contend the issue is how to effectively set this tolerance if the issue is small diagonal elements compared to other diagonal elements.

Thanks for any insights.


问:
Using absolute error as the cost function
This is related to my previous post Minimize distance between curves.

I have a dataset with values of multiple curves. An example plot is shown below.
I want to scale the curves (move up/down) so that all curves overlap.

The following is a sample dataset which includes that data points corresponding to 5 curves and coordinate inputs below

```<code>scale =  1.5;
x1 = [0,4,6,10,15,20]*scale;
y1 =  [18,17.5,13,12,8,10];
x2 = [0,10.5,28]*scale;
y2= [18.2,10.6,10.3];
x3 = [0,4,6,10,15,20]*scale;
y3 = [18,13,15,12,11,9.6];
x4 = [9,17,28]*scale;
y4 = [5,5.5,7];
x5 = [1,10,20]*scale;
y5 = [3,0.8,2];

plot(x1,y1, '*-', x2, y2,  '*-', x3, y3,  '*-', x4, y4,  '*-', x5, y5,  '*-')
</code>```

To scale the curves, I need to find the scale factor for each curve.

In the answer posted in the previous post (https://scicomp.stackexchange.com/a/41555/29087) the following has been suggested.

$$
\min_{a} \left(\sum_i (  a \cdot f_i -g_i)^2\right)
$$
which leads to a single equation
$$
0 = \partial_a \left(\sum_i (  a \cdot f_i -g_i)^2\right) = 2 \sum_i f_i (a \cdot f_i -g_i)
$$
and thus
$$
a = \frac{\sum f_i g_i}{\sum f_i^2}
$$

In the above answer, the cost function is the sum of absolute error, and differentiating this gives a.

Excuse me for the naive question,

I would like to understand how the cost function has to be defined if we want to minimize absolute error instead of squared error and how to estimate the scale factor.

Suggestions will be really appreciated.

EDIT:
$\sum_i abs(a*f_i - g_i)$

min: $\sum_i u_i$
such that,

$f_i * a - u_i <= g_i$,

$-f_i * a - u_i <= -g_i$

I started setting up the problem like the below

```<code>% Loss Function (absolute error loss)
% input curves
scale =  1.5;
x1 = [0,4,6,10,15,20]*scale;
y1 =  [18,17.5,13,12,8,10];
x2 = [0,10.5,28]*scale;
y2= [18.2,10.6,10.3];

% y2 is the target function and y1 is the function to be shifted
f = y1;
g = interp1(x2,y2,x1);

% linprog inputs
% x = linprog(f,A,b)
% i = 1:length(x1)
A = [fi - 1; -fi -1];
x = [a ui];
b = [gi -gi];

% obj = sum_i u_i (function to minimize)
% solve system
</code>```

But I am not sure how to proceed after this and solve this system in MATLAB's/scipy's linprog since the function has unknown parameters.

It will be helpful if suggestions are given for solvig the above using the linprog function in scipy library.


问:
Centered finite volume scheme for an advective term on an unstructured/irregular/non-uniform grid
Consider the continuity equation

$$\frac{\partial u}{\partial t} + \frac{\partial \Phi}{\partial x} = 0$$
$$\Phi = au + b\frac{\partial u}{\partial x}$$

Suppose I want to solve the above using centered finite volumes. Below, $U_j \approx \overline{u}_j$ where $\overline{u}_j$ is the spatial average value of $u$ in a finite volume. $\Phi_{i+1/2}$ and $\Phi_{i-1/2}$ are the fluxes at the edges of the finite volume, and $\Delta x_i$ is the thickness of the the volume.

$$\frac{\partial U_j}{\partial t} = -\frac{\Phi_{i+1/2} - \Phi_{i-1/2}}{\Delta x_i}$$

For a structured grid, I can use the following centered approximation for the fluxes

$$\Phi_{i+1/2} = a_{j+1/2}\frac{U_{j+1} + U_{j}}{2} + b_{j+1/2}\frac{U_{j+1} - U_{j}}{\frac{1}{2}\Delta x_{j+1} + \frac{1}{2}\Delta x_{j}}$$

For the advective term above, I have approximated $U_{j+1/2}$ as

$$U_{j+1/2} = \frac{U_{j+1} + U_{j}}{2}$$

My question is what is a more general centered approximation of $U_{j+1/2}$ for an unstructured grid? One possibility is to linearly interpolate to $U_{j+1/2}$ between the center of the finite volumes containing $U_{j+1}$ and $U_{j}$. The result is

$$U_{j+1/2} = \frac{U_{j+1} \Delta x_j + U_{j} \Delta x_{j+1}}{\Delta x_{j+1} + \Delta x_{j}} \tag{1}$$

A reference equation in a textbook would be appreciated.

Thanks!


问:
Technical term for getting stuck at local optimum when building a mechanistic model
I'm looking for a technical term or literature that describes a problem when creating mechanistic computer models.

By mechanistic computer model, I mean a computer model of a physical system where the different components of the physical system are modelled by components of the computer model, each model component represents a different mechanism in the system.

The error in a mechanistic model will be caused by different mechanisms being poorly represented or missing entirely.

The term I'm looking for is to describe the following situation (or literature that addresses this problem):

You improve one of the components of the model, say by adding detail that you know better represents the underlying physical process. However, the model performance as a whole gets worse, because the error in the component you improved was cancelling out the error in another component of the model.


问:
deal.II and curved faces: how can I get the curved description
I'm not a deal.II expert, and while studying step-6 I was reading the documentation of the MappingQ1 class in the deal.II documentation. At some point in the description (https://www.dealii.org/current/doxygen/deal.II/classMappingQ1.html) , there's written:

Note, however, that in 3D the faces of a general, trilinearly mapped cell may be curved, even if the edges are not

I have two main issues with this:

I don't see why this is mathematically true. I don't see how the image of a trilinear map $\phi(x,y,z)$ gives some curved boundary.

Considering the case of an hyper_ball, i.e. a ball in 3D. I see indeed that boundary faces are not straight quadrilaterals (see the attached image). How can I retrieve the "curved" part of the face, from a software standpoint? I mean, I can of course ask to the Mapping object the vertices of the cell or face, but this are not enough to describe the curvature... I don't see how one can query the "curved part" in terms of Points!

I understand somehow this has to do with the concept of Manifold, but I'd like to see how in the code one can actually got this "curved" information.

Any help is highly appreciated.


问:
Numerical solution to integro-differential equation
The time dynamics of an atom interacting with a reservoir of spectral density $J(\omega)$ are obtained by solving the following integro-differential equation:

$$ \frac{\mathrm{d}c(t)}{\mathrm{d}t} = - \int\limits_0^t \mathrm{d}\tau~ c(\tau)f(t-\tau),$$

where

$$f(t-\tau) = \int\limits_0^{\infty}\mathrm{d}\omega ~ J(\omega) \exp(i(\omega_E-\omega)(t-\tau)) $$

I have an analytical solution when $J(\omega)$ is a Lorentzian distribution, however, is there a way to solve the above numerically for $c(t)$, given an arbitrary spectral density?

In particular I don't have an analytical expression for $J$, only the numerical data.


问:
Solution to minimization problem when variables factor in 2 analytical problems
I am asking a follow up question to this question, but I could probably have written it as an answer instead. However, I don't know if what I am doing here makes sense or is too complicated for my problem so I am not sure it's the correct answer!

As I explained in the other question, I have basically a very standard unconstrained function minimization problem. Let's say that the function to minimize $F$ depends on 2 matrices ($A$ and $B$) and that I can compute the gradients $F_A$ and $F_B$. I could directly combine the gradients and use a standard solver but I know that I can solve the equations $F_A=0$ and $F_B=0$ analytically when only one variable is allowed to change. So I am wondering if there is a better approach to my problem.

More specificaly, let's assume without loss of generality that my iteration starts at $(0,0)$ and I am looking to find a better estimation of my minimum. I find $A_0$ and $B_0$ as the matrices that solve $F_A(A_0,0)=0$ and $F_B(0,B_0)=0$. Then I define an auxiliary function $f(x,y)=F(A_0x,B_0y)$ where $(x,y)\in\mathbb R^2$. The idea is that, as the Hessian is positive definite, an approximate minimum should happen for a certain value $0<x<1$ and $0<y<1$.

The reason why I use this auxiliary function is that it computes in a very straightforward way the value of my function on a much smaller manifold, with most dimensions aligned to the gradient at two points.

I can then compute:
$$\frac{\partial f}{\partial x}=\mathrm{tr}(A_0^TF_A(A_0x,B_0y)$$
$$\frac{\partial f}{\partial y}=\mathrm{tr}(B_0^TF_B(A_0x,B_0y)$$
These functions are easy to compute, for example at the four corners of the square $(x,y)\in[0,1]^2$. This gives me $8$ values, with $2$ of them being zeros (they match the condition that $f_x(1,0)=f_y(0,1)=0$).

If I suppose that $f(x,y)$ is a conic section, I can write:
$$f(x,y)=\frac12ax^2+bxy+\frac12cy^2+dx+ey+g$$
And we can get the partial derivatives:
$$f_x=ax+by+d$$
$$f_y=bx+cy+e$$

There are some conditions on the values, for instance as $f_x(1,0)=a+d=0$ and $f_y(0,1)=c+e=0$, I can reduce to only $3$ variables. There are several ways to estimate these but in the end, it's just a simple $2\times2$ matrix linear problem.

Is it a standard approach, to assume that the solution to the bigger problem is simply the combination of two orthogonal partial solutions? I have honestly not seen this approach to solve minimization problems but it seems to make sense in my situation, as my earlier approach is really zigzagging between two minimum values. I still have to try it on my real data though.


问:
Sparse generalized symmetric eigensystem solver
Can anyone recommend a good software for solving generalized symmetric eigenvalue problems of the form,
$$
A x = \lambda B x
$$
where $A,B$ are symmetric and sparse, and $B$ is positive definite?

I have only 1 machine so I would prefer something based on OpenMP rather than MPI if possible.


问:
stable solutions for Large-scale ODEs under boundary value problem
I'm doing FEM and have a problem about getting numerically stable solution for ODEs problems like:
$$ \frac{\mathrm{d}}{\mathrm{d}x}\mathbf{Y} = \mathbf{AY}, x\in[x_1,x_2]$$
in which $\mathbf{Y}$ stands for a large-scale vector (greater than $1000000\times1$), and $\mathbf{A}$ is its corresponding coefficient matrix. The elements in matrix $\mathbf{A}$ vary greatly, from $10^{-50}$ to $10^{50}$.

For small matrix, analytical solutions for such ODEs are simple. However, those does not work for large scale matrix, especially when the transfer distance $(x_2-x_1)$ is large. I used DQM to get a relatively stable solution, but still is not stable enough when the dimension of $\mathbf{A}$ gets bigger.

Is there any stable solutions/ideas for this problem? Thanks!


问:
What are the known or used numerical methods for integration over the sphere $S^2$ ? and what about over $S^3$?
What are the numerical methods available to compute integrals over $S^2$, for the particular integration :
$$\int_{S^2} f(\omega)\,d\sigma(\omega)\ , \quad \text{ where $d\sigma$ is the usual measure on the sphere }S^2 $$
And $\omega= (\theta, \phi)$ is a direction (or a unit vector) given by a point on this sphere.

I am guessing an algorithm can proceed by computing Riemann sums as the usual single variable integration? Then maybe a trapezoidal approximation is used?

Any idea on how to implement this on Matlab or Mathematica ?


问:
Finding spectrum of a Kronecker factored + block-partitioned matrix
I have dense $d\times d$ matrices $A$, $B$, $C$ with $d\approx 1000$ and want to find the top $10^5$ eigenvalues of the following positive definite matrix:

$$
\Sigma=
\left(\begin{matrix}
A\otimes A & C^T\otimes C\\
C\otimes C^T & B\otimes B
\end{matrix}\right)
$$

Cheap: assume $C=0$ and compute $\lambda(A)$, $\lambda(B)$, since $\lambda$ commutes with $\otimes$.

Expensive: power method. Storing $10^5$ eigenvectors is 40 GB.

Is there something in between the two?

Motivation

$\Sigma$ is a

covariance of $x\otimes x$, where $x$ is a partitioned Gaussian vector with mean 0 and
$
E[xx']=
\left(\begin{matrix}
A&C^T\\
C&B
\end{matrix}\right)
$

Gauss-Newton block for neural network linear layer if we assume forward/backward values are jointly Gaussian


问:
Automatic differentiation necessary for large optimal control problems?
I am investigating ways to solve an optimal control problem in an embedded way, preferably in Java. The system is modeled with triple integrator dynamics $u=\dddot{x}$ and solved with multiple shooting, with continuity constraints put on $x$, $\dot{x}$ and $\ddot{x}$ which are solved with Runge-Kutta-4 integration and optimized with an NLP solver.  For $N$ shooting points I have $4N$ variables and about $3N$ equality constraints with the Runge-Kutta integration expressed within the constraints. Safe to say it quickly becomes computationally heavy, as I usually do $N=100$.

I've found it is quite difficult to express numerical integration within constraints to "normal" optimization libraries in a way that works (imsl and knitro). And when I do manage it becomes incredibly inefficient. I've been able to do it efficiently in python with CasADi, which handles integration well, but also performs automatic integration. The only other way I've been able to solve the problem is through AMPL which also supplies the solver (ipopt in this case) with automatic gradients. Other than AMPL the set of available auto-diff libraries for java is very small and none of them as far as I can tell would be viable to use with this kind of problem (here is a list https://www.autodiff.org/?module=Tools&language=Java ).

For instance, ipopt requires one to supply:

objective function,

gradient of the objective,

constraint function values,

Jacobian of the constraints,

Hessian of the Lagrangian function,

which many of I think would not be easy enough to compute with those libraries listed to gain any speed.
So having something that interfaces natively with such a solver would of course be the best.

Now could it be that this problem can be solved in a reasonable time without auto-diff, just that I've implemented it in a bad way in those "traditional" libraries? To my understanding, auto-diff should be faster for problems with a lot of variables and constraints. But how much is a lot of variables and constraints in reality? Would it be feasible to allow a solver to do finite differences instead, even if this problem has to be solved quickly? I don't have much experience working with optimal control so I don't really have any reference points here. Anyone with experience working with similar problems knows if you used auto-diff or not and what realistic difference in speed can one expect here?


问:
Advantages of using BEM for contact problem over FEM
Are there any advantages of using BEM (Boundary Element Method) for contact problem over FEM (Finite Element Method) apart from lesser computational cost?</p>


问:
SLSQP solver scipy with linear subset constraints
I have been trying to solve a least squares problem of the following form:

$$
\begin{equation}
\min_{\vec{x}} \frac{1}{2} \lVert f(\vec{x}) - f_{\text{target}} \rVert_{2}^2 + \alpha\Big( \frac{1-\rho}{2} \lVert \vec{x} - \vec{x}_{\text{orig}} \rVert_{2}^{2} + \rho \lVert \vec{x} - \vec{x}_{\text{orig}} \rVert_{1} \Big),
\end{equation}
$$

where $f$ is a continuous function $ f: \mathbb{R}^{N} \to \mathbb{R}^{K}$ and $\vec{x} \in \mathbb{R}^{N}$. The first term in the objective function above is hence minimizing the distance between the vector $f(\vec{x})$ and a given vector $f_{\text{target}}$. The next two terms are regularization terms ($\lvert \cdot \rvert_{1,2}$ are the $L_{1}$ and $L_{2}$ norms respectively) . $\vec{x}_{\text{orig}}$ is the starting vector.

I have been able to get a solution in a fairly quick time by using the scipy.optimize function with SLSQP solver. However, the addition of linear constraints seems to not work. The constraints are of the form

$$
\begin{align}
l_{1} &\leq \sum_{i \in I_{1}} x_{i} \leq u_{1} \\
l_{2} &\leq \sum_{i \in I_{2}} x_{i} \leq u_{2} \\
&\vdots \\
l_{n} &\leq \sum_{i \in I_{n}} x_{i} \leq u_{n} \\
\end{align}
$$

where $I_{m}$ are partitions of $\{1 \cdots n \}$. So the constraints are on the sums of the elements of the vector $\vec{x}$.

I have attempted to integrate these constraints within the function call optimize through the constraints list.

```<code>def fun_upper_constraint(x, indices, ub):
    return ub - np.sum(x[indices])

def fun_lower_constraint(x, indices, lb):
    return np.sum(x[indices]) - lb

def generate_constraint_dict(indices, ub, lb):
    return {"type": "ineq", "fun": fun_upper_constraint, "args": (indices, ub)}, {
        "type": "ineq",
        "fun": fun_lower_constraint,
        "args": (indices, lb),
    }

</code>```

Adding constraints this way causes the optimizer to not find the optimum with the status message Positive directional derivative for linesearch. Hence, I wonder whether there is a way to encode these constraints in a way that is compatible with what the SLSQP solver expects.


问:
Magnetic field simulation for "quasi-solenoid" - Where to start?
I would kindly ask you which (preferably free) programs / codes do you suggest for the numerical simulation of the problem described below. I am not asking for the full solution, but just want to avoid “reinventing the wheel” as much as possible. I have access to Matlab and Python (Anaconda), which means that I can install programs like Netgen/NGSolve, Salome etc.

For the case of “quasi-solenoid”, through which I put a strong, but short current (discharged from capacitor), I want to calculate the spatial distribution of magnetic and (using the Maxwell–Faraday equation) electric field at the moment of strongest current.
I used the prefix quasi as I want to specify that the distance between the wires and solenoid axis isn’t always constant, but varies quite substantially, though I still believe that I can describe the spatial distribution of wires using some parametric function or something else. Furthermore, though for now this isn’t so relevant, the cross-section of the wire isn’t circular, but rectangular.

Searching through the internet I already found and used some codes, such as https://www.mathworks.com/matlabcentral/fileexchange/50434-biot-savart-magnetic-toolbox, but they do not cover everything that is relevant for me. For example, I doubt that we can apply the Biot-Savart law, which this code uses, for my case.

Any suggestion would be appreciated.


问:
A direct numerical method for determining a relaxation function from a known creep function
Both the creep function and the relaxation function are connected by the convolution integral. The usual method for calculating the relaxation function from a creep function or vice versa is to transfer the creep function into the Laplace domain. In the Laplace domain, the creep function is equal to 1/relaxation.
In my case, my material model is ill-conditioned in the Laplace domain, so I am not able to transfer it back from the Laplace domain to the time domain. Does anyone know of a direct method to determine the relaxation function from a known creep function?
I could not find an answer in the literature, maybe there is a paper or a way to transfer it directly.
The link of my main question.
Numerical instability in the inverse Laplace transform

As a friend here sed im going to write the math problem im facing now.
This is my problem again
Does anyone know a suitable way to transfer this function from the Laplace domain to the time domain?
Working with python and Mathematica did not help me to find the correct answer.

```<code>import mpmath as mp
import numpy as np
import matplotlib.pyplot as plt
def Laplace_Max(s):
    #mp.dps = 100
    A= s**(5/4)/ (
        0.00480931 +
        0.0244077*s**(1/4) +
        0.0129056*mp.exp(-35*s)*s**(1/4) +
        0.00329167*mp.exp(0.707997*s)*s**(1/4) * mp.gammainc(0.0, 35.708*s,mp.inf, regularized=True) -
        0.00530593*mp.gammainc(1.25, 35*s,mp.inf,   regularized=True)
    )
    return A
t_R_inv = np.linspace(1,1200,12000)
Max_R=np.zeros(len(t_R_inv))
for i in range(len(t_R_inv)):
    Max_R[i]=mp.invertlaplace(Laplace_Max, t_R_inv[i], method = 'stehfest', dps = 100, degree = 1000)
plt.plot(Max_R)
</code>```

Update
The inverse laplace of this function!

$\frac{1}{s}^{\frac{5}{4}}\cdot (0.00480931 + 0.02440766 \cdot s^{\frac{1}{4}} +
   0.0129056 \cdot e^{-35 \cdot s} \cdot s^{\frac{1}{4}}+
   0.003291670 \cdot e^{0.7079967\cdot s} \cdot s^{\frac{1}{4}}
    \cdot \Gamma[0, 35.7079967 \cdot s] -
   0.00530593 \cdot\Gamma[1.25, 35\cdot s])$


问:
Discontinuous Galerkin failing to converge Euler equations under p-refinement
I am solving the steady state compressible Euler equations in conservation form in 2D in a rectangular domain with a discontinuous Galerkin (DG) code I have developed. As a test, the boundary conditions are Dirichlet conditions set to the same far-field values. The Ma number is much less than 0.3.

I expect the solution to simply be the far-field values at every point in the domain. Indeed, for constant (q = 0) polynomials defined on the domain, I obtain this solution across both coarse and fine meshes. However, when the polynomial order is increased (q >= 1), the preconditioned Krylov solver (GMRES) can not satisfy the Newton system which describes the DG implementation of the Euler equations. Instead, the solution for all conservation variables is highly oscillatory. The solution on the boundary is correct, however.

Finally, my question is, assuming I have implemented the Newton system correctly (which I have painstakingly checked many, many times), what might be cause of this behavior? Allowing the Krylov solver more iterations decreases the residual of the solution to the Newton system slightly, but doesn't provide indication that the solver would converge given enough iterations. Considering conservation equations for other flux formulations, e.g. scalar advection in a constant velocity field, will converge under some p-refinement.

I suppose I should mention that I have implemented a hybridized DG scheme, where the Riemann solver is implicitly defined by a stabilization term on the numerical flux and I haven't seen flux limiters mentioned. The literature reports that this stabilization should have less effect on the final solution as the polynomial order increases. I have found little change in the behavior of the solution for various stabilization values beyond a small threshold value.

Edit:

In HDG, we introduce additional variables defined on the interfaces of the elements. We also define a numerical flux function defined on the interfaces written in terms of the trace and element variables which is given by

\begin{equation}
    f_h \left( u, \hat{u}, n \right) = F(\hat{u}) \cdot n + \sigma \left( u, \hat{u}, n\right) \left( u - \hat{u} \right)
\end{equation}

where $u$ and $\hat{u}$ are the conservation variables defined on the element and trace, respectively, $F(q)$ is the Euler flux, and $\sigma$ is the stabilization matrix which implicitly defines the Riemann solver. I found the most success setting $\sigma = \tau I$ for large values of $\tau$, though the literature presents other formulations which include characteristics of $\frac{\partial F(q)}{\partial q} \cdot n$. However, these have precluded convergence in my case. The boundary conditions are imposed by solving

\begin{equation}
   \int_{\partial \Omega} b_h \left(u, \hat{u}, n \right) \cdot \mu = 0
\end{equation}

on the domain boundaries where $\mu$ is a test function in the same space to which $\hat{u}$ belongs and $b_h \left(u, \hat{u}, n \right) = f_h \left( u, \hat{u}, n \right) + f_h \left( u_b, \hat{u}, -n \right)$ with $u_b$ denoting the boundary state. After non-dimensionalization, $u_b$ is set to $\left[ \rho: 1, \rho u: 0.1, \rho v: 0, E: 1.78 \right]^T$.


问:
Convergence rate for not smooth solution with classical $P^1$ Lagrangian FEM
I'm using classical $P^1$ finite elements to solve $- \Delta u = f$ with Dirichlet BC in a 2D domain $\Omega$.

I know from theory that the solution is not in $H^1$ for my particular choice of $f$, so all classical convergence estimates are useless. If I compute the $L^2$ error against the analytical solution indeed I don't see the correct slope, as expected. However, after I reach a certain number of DoFs the errors start increasing.

The rates are the following:

```<code>dofs         L2
    81 3.214e-02
   289 1.735e-02
  1089 4.777e-03
  4225 2.633e-03
 16641 1.330e-03
 66049 9.378e-04
263169 1.256e-03
1050625 1.315e-03
</code>```

As you see, the $L^2$ error at some point starts growing. Can this be justified by the low regularity of the solution?

If I take a smooth solution, I have correct rates as I refine the grid, and I don't observe this phenomenon.


问:
Ratio of error norms or norm of error ratio in adaptive step size control?
Step size controllers for ODE solvers with adaptive step size usually track an error estimate $y_{\mathrm{err}}$ and compare it to the current state $y_\mathrm{current}$ to decide if a step can be accepted or the step size needs to be decreased. We consider the case of multi-dimensional ODEs, so all $y$ are vectors in $\mathbb{R}^n$ for some $n$. There seem to be two ways to do this comparison. The first checks if the norm of the ratios between error and some tolerable error scale is less than one, i.e.
$$
\left\|\frac{y_{\mathrm{err}}}{\mathrm{atol} + \mathrm{rtol} \cdot y_\mathrm{scale}}\right\| < 1
$$
where all operations are elementwise and $\left\|\cdot\right\|$ is some norm. $\mathrm{atol}$ and $\mathrm{rtol}$ are the absolute and relative tolerance, respectively, and $y_\mathrm{scale}$ is something like $\max(\operatorname{abs}(y_\mathrm{current}), \operatorname{abs}(y_\mathrm{next}))$. The advantage of this approach is some built-in scale invariance. If some components of $y$ live on a much smaller scale than others, these will cancel out and the components will be considered more evenly in the step acceptance decision. However, from a mathematical perspective, the formula above is quite nonsensical, because elementwise division of vectors is not a meaningful operation in normal vector spaces. Under the consideration of mathematical sensibility, it would make more sense to check for step acceptance using
$$
\frac{\left\|y_{\mathrm{err}}\right\|}{\mathrm{atol} + \mathrm{rtol} \cdot \left\|y_\mathrm{scale}\right\|} < 1.
$$
Now we are taking norms of vectors living in the same vector space that we are solving the equation in and building ratios of real numbers, all of which are well-defined. Nonetheless, I almost exclusively see the former approach in real code and not the latter, and I have to agree that the scale-invariance argument makes a lot of sense in practice.

Is there some theoretical argument to guide the decision on when to prefer one over the other?

For example, diffrax, a python ODE solver for the JAX framework, chooses norm of ratio as described in their code.


问:
Does cblas_dgemm mutate my input matrices?
I have written a matrix class Matrix<T> for which I have implemented a wrapper function for cblas_dgemm.
cblas_dgemm requires to pass my input matrices as double*. However, I want my wrapper to work with const Matrices as well (for which I can only provide a const double*). Within my call to cblas_dgemm I cast the const away with const_cast<double*>. My question is: Can I rely that cblas_dgemm actually only reads my input matrices, so that the const_cast is safe? I know, of course, that eventually the input matrices are the same after the call to dgemm, but maybe BLAS does some temporary changes to my matrix?

```<code>Matrix<double> wrapper_to_cblas(const Matrix<double>&a, const Matrix<double>& b) {
...
cblas_dgemm(/*here I cast const double* to double* */)
...
}
</code>```

Edit: The BLAS implementation I am using actually allows for the use of const double*, so it really seems to be implementation-dependent.


问:
Error in implementation of Crank-Nicolson method applied to 1D TDSE?
Some context, I've posted this question on physics SE and stack overflow. The former had nothing to offer, the latter had a great commenter that agreed with the phase looking off being one of the potential issues and some difference in the solution's dynamics (and also alerted me to the existence of this part of SE, which I had no idea existed!). Hopefully the question has a bit more success here.

I'm attempting to use Crank-Nicolson scheme to solve the TDSE in 1D. The initial wave is a real Gaussian that has been normalised with respect to its probability density. As the solution evolves, a depression grows in the central peak of the real part of the wave, and the imaginary part's central trough is perhaps a bit higher than I expect (image below).

Does this behaviour seem reasonable? I have searched around and not seen questions/figures that are similar. I've tested another person's code from Github and it exhibits the same behaviour, which makes me feel a bit better. But I still think the center peak should just decrease in height and increase in width. The likelihood of me getting a physics-based explanation is relatively low here I'd assume, but a computational-based explanation on errors I may have made is more likely.

I'm happy to give more information. Thanks in advance!

Here's a link to GIF of time evolution with positive dt:

https://imgur.com/a/9W5MjMC

GIF with negative dt:

https://imgur.com/H4kf4HA

This changing of sign of dt came about from seeing this: https://scicomp.stackexchange.com/q/6725
and testing how it changes the solution (which it does a fair bit!).

An update from hardmath's comment: the equation I'm attempting to solve is the time-dependent Schroedinger equation, with no potential.
$$i \hbar  \frac{ \partial }{ \partial t}  u(x, t) = -  \frac{ \hbar^{2} }{2m}  \frac{\partial^2}{\partial x^2} u(x, t)$$
Applying the Crank-Nicolson method results in (for me):
$$u_{j}^{n+1}+\frac{1}{2}\frac{i\Delta t}{\hbar}\frac{-\hbar^{2}}{2m}
\frac{u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{\Delta x^2}=
u_{j}^{n}-\frac{1}{2}\frac{i\Delta t}{\hbar}\frac{-\hbar^{2}}{2m}
\frac{u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{\Delta x^2}$$
Superscript is for time stepping, subscript for space stepping.
$$A_{l} u^{n+1} =  A_{r} u^{n}$$
where:
$$A_{l} = \begin{pmatrix}
1-\alpha & \alpha & \cdots & \cdots & 0 \\
\alpha & 1-2\alpha & \alpha & \cdots  & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & \alpha & 1-2\alpha & \alpha \\
0 & \cdots & \cdots & \alpha & 1-\alpha \\
\end{pmatrix}$$
$$A_{r} = \begin{pmatrix}
1-\beta & \beta & \cdots & \cdots & 0 \\
\beta & 1-2\beta & \beta & \cdots  & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & \beta & 1-2\beta & \beta \\
0 & \cdots & \cdots & \beta & 1-\beta \\
\end{pmatrix}$$
The coefficients $\alpha$ and $\beta$ collect terms arising from applying the Crank-Nicolson method to the form of the Schroedinger equation I'm attempting to solve here (and working with natural units $\hbar=m=1$):
$$\alpha = \frac{1}{4}  \frac{i  \Delta t}{ \Delta  x^{2} }$$
$$\beta = -\frac{1}{4}  \frac{i  \Delta t}{ \Delta  x^{2} }$$
The matrices are diagonally symmetric, except the corners, where the boundary conditions imposed (Dirichlet: $u(x_{min},t)=u(x_{max},t)=0$) result in elements 'outside' the matrices being known as $0$ and changing the corner elements appropriately.

The wave function is $u$, and its superscript indicates its time: so if $n$ is the current time, I am using the Crank-Nicolson method to estimate the wave function at the next time after taking timestep $dt$, $n+1$.

And the part of my code relevant to solving the 1D TDSE:
(pretty much the entire thing except the plotting)

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

# Define function for norm.
def normf(dxc, uc, ic):
    return sum(dxc * np.square(np.abs(uc[ic, :])))

# Define function for expectation value of position.
def xexpf(dxc, xc, uc, ic):
    return sum(dxc * xc * np.square(np.abs(uc[ic, :])))

# Define function for expectation value of squared position.
def xexpsf(dxc, xc, uc, ic):
    return sum(dxc * np.square(xc) * np.square(np.abs(uc[ic, :])))

# Define function for standard deviation.
def sdaf(xexpc, xexpsc, ic):
    return np.sqrt(xexpsc[ic] - np.square(xexpc[ic]))

# Time t: t0 =< t =< tf. Have N steps at which to evaluate the CN scheme. The
# time interval is dt. decp: variable for plotting to certain number of decimal
# places.
t0 = 0
tf = 20
N = 200
dt = tf / N
t = np.linspace(t0, tf, num = N + 1, endpoint = True)
decp = str(dt)[::-1].find('.')

# Initialise array for filling with norm values at each time step.
norm = np.zeros(len(t))

# Initialise array for expectation value of position.
xexp = np.zeros(len(t))

# Initialise array for expectation value of squared position.
xexps = np.zeros(len(t))

# Initialise array for alternate standard deviation.
sda = np.zeros(len(t))

# Position x: -a =< x =< a. M is an even number. There are M + 1 total discrete
# positions, for the points to be symmetric and centred at x = 0.
a = 100
M = 1200
dx = (2 * a) / M
x = np.linspace(-a, a, num = M + 1, endpoint = True)

# The gaussian function u diffuses over time. sd sets the width of gaussian. u0
# is the initial gaussian at t0.
sd = 1
var = np.power(sd, 2)
mu = 0
u0 = np.sqrt(1 / np.sqrt(np.pi * var)) * np.exp(-np.power(x - mu, 2) / (2 * \
                                                                        var))
u = np.zeros([len(t), len(x)], dtype = 'complex_')
u[0, :] = u0

# Normalise u.
u[0, :] = u[0, :] / np.sqrt(normf(dx, u, 0))

# Set coefficients of CN scheme.
alpha = dt * -1j / (4 * np.power(dx, 2))
beta = dt * 1j / (4 * np.power(dx, 2))

# Tridiagonal matrices Al and AR. Al to be solved using Thomas algorithm.
Al = np.zeros([len(x), len(x)], dtype = 'complex_')

for i in range (0, M):
       Al[i + 1, i] = alpha
       Al[i, i] = 1 - (2 * alpha)
       Al[i, i + 1] = alpha

# Corner elements for BC's.
Al[M, M], Al[0, 0] = 1 - alpha, 1 - alpha

Ar = np.zeros([len(x), len(x)], dtype = 'complex_')

for i in range (0, M):
       Ar[i + 1, i] = beta
       Ar[i, i] = 1 - (2 * beta)
       Ar[i, i + 1] = beta

# Corner elements for BC's.
Ar[M, M], Ar[0, 0] = 1 - beta, 1 - beta

# Thomas algorithm variables. Following similar naming as in Wikipedia page.
a = np.diag(Al, -1)
b = np.diag(Al)
c = np.diag(Al, 1)

NT = len(b)

cp = np.zeros(NT - 1, dtype = 'complex_')

for n in range(0, NT - 1):
    if n == 0:
        cp[n] = c[n] / b[n]
    else:
        cp[n] = c[n] / (b[n] - (a[n - 1] * cp[n - 1]))

d = np.zeros(NT, dtype = 'complex_')

dp = np.zeros(NT, dtype = 'complex_')

# Iterate over each time step to solve CN method. Maintain boundary
# conditions. Keep track of multiple values.
for i in range(0, N):
       # BC's.
       u[i, 0], u[i, M] = 0, 0

       # Find RHS.
       d = np.dot(Ar, u[i, :])

       for n in range(0, NT):
           if n == 0:
               dp[n] = d[n] / b[n]
           else:
               dp[n] = (d[n] - (a[n - 1] * dp[n - 1])) / (b[n] - (a[n - 1] * \
                                                                  cp[n - 1]))

       nc = NT - 1
       while nc > -1:
           if nc == NT - 1:
               u[i + 1, nc] = dp[nc]
               nc -= 1
           else:
               u[i + 1, nc] = dp[nc] - (cp[nc] * u[i + 1, nc + 1])
               nc -= 1

       norm[i] = normf(dx, u, i)
       xexp[i] = xexpf(dx, x, u, i)
       xexps[i] = xexpsf(dx, x, u, i)
       sda[i] = sdaf(xexp, xexps, i)

# Fill in final norm value.
norm[N] = normf(dx, u, N)

# Fill in final position expectation value.
xexp[N] = xexpf(dx, x, u, N)

# Fill in final squared position expectation value.
xexps[N] = xexpsf(dx, x, u, N)

# Fill in final standard deviation value.
sda[N] = sdaf(xexp, xexps, N)
</code>```


问:
Methods for binary linear programming
I have an LP problem (linear objective with eq and ineq constraints) in binary variables.

Except for the objective, all the coefficients are integer, mostly in {-1,0,1}. Maybe the objective coeff could be discretized.

I usually use an MI solver (gurobi). Since it solves a moderate problem in my case under a second, I suspect that it may have a special structure that permits employing a different specialized method that will be faster. Can you suggest other methods to solve my problem?

More details in a duplicate question:

https://math.stackexchange.com/questions/4504572/methods-for-binary-linear-programming/4505035

Title: Methods for binary linear programming


问:
Method of alternating projections for three sets when intersection is empty
The method of alternating projections finds a point in the intersection of two (or more) convex sets. If the intersection is empty, and the number of sets equals two, it converges to a pair of points, one in each set, that are as close as possible. What happens when there are three sets, and only one pair of sets has no intersection?

This paper from 2018 reviews the state of the art and I think the conclusion is that there is no simple answer.


问:
Solving simplified 1D plasma fluid equations with finite difference
The following two equations represent a simple model of a plasma where ions are immobile.

$$n\frac{\partial u}{\partial t}+nu\frac{\partial u}{\partial x}=n\frac{d\phi}{dx}-\theta\frac{\partial n}{\partial x}$$

$$\frac{\partial n}{\partial t}+\frac{\partial}{\partial x}(nu)=0$$

Here, $\theta$ is a constant denoting the temperature of the electron fluid and $\phi$ is the electric potential function obtained from an elliptic PDE (Poisson equation).

I'm looking for a (preferably) simple finite difference scheme to solve this model numerically. In particular, I'd like to see plasma oscillations.

I recently became familiar with Lax–Friedrichs method and other explicit methods for solving hyperbolic PDEs. As far as I know they can be used to solve a system of the form

$$\frac{\partial v}{\partial t}+\frac{\partial f(v)}{\partial x}=0$$

where in general $v$ and $f$ are both vector-valued functions. I'm not sure if the equations in question can be written in this form.

Let's consider two special cases. If $\theta=0$ and $n\neq 0$, the first equation can be written as

$$\frac{\partial u}{\partial t}+\frac{\partial}{\partial x}\left(\frac{1}{2}u^2-\phi\right)=0$$

And if we neglect the electric field term, the first equation can be written as

$$\frac{\partial}{\partial t}(nu)+\frac{\partial}{\partial x}\left(nuu+\theta n\right)=0$$

Both of these equations seem to be in the form mentioned above.

What is the standard or most common finite difference scheme to solve the set of the two PDEs at the beginning of this question?


问:
How can Android be utilized in computational science pipelines?
How can Android be utilized in computational science pipelines?

For example, if one wanted to produce a visualization for an Android mobile device. Similarly to how people are using browsers and JS to do visualizations and interactive apps even.

I'm thinking that support for numerical codes in Java is poor compared to C/C++, which would make me think that a way to utilize C/C++ code would be to run these in a server and just poll results for displaying on a Kotlin app or something. But is this all there is to this?


问:
bingham plastic model implementation
I'm looking for a book/paper/code which describes the implementation of a Bingham plastic fluid model. I've tried deriving the equations for implementation using the style of the 1d visco-plastic model for solids described in the first chapter of Simo-Hughes but it seems to me that the analogy between solids and fluids is imperfect, because the Bingham plastic model lacks an elastic spring.

Edit: I'm looking for an exact implementation, not regularized, analogous to the return mapping algorithms of plasticity.

Title: bingham plastic model implementation


问:
How to find an alignment between two overlapping time series?

Context

I'm working with absorbance spectral data. I have two spectra generated for a single analyte. Each spectrum was collected over two overlapping wavelength regions e.g. 8-10 micron and 9-11 micron. The spectrums both contain a few thousand data points where each point is a real number.

Question

Is there a method for the partial alignment of two time-series? One that would align a suffix of one spectrum to the prefix of another? Additionally, is there a way this could be paired with something like DTW? Due to the nature of the instrument, there's usually some misalignment of movement between any two spectrums.

Any general suggestions would be appreciated, however I'm currently working in python so any specific libraries would be additionally helpful.

Edit w/ more context:

There is not a correlating wavelength for each data point. I only have the known range and the data sequence itself. The instrument has a laser which scans over a specified wavelength range. Data collection is passive and begins and stops with the laser within some negligible time difference. However, the sampling frequency can vary slightly throughout the scan which stops a linear interpolation of the correlating wavelengths from being consistent. Consequently, the resulting data then has some degree of variability in its length, usually on the order of 10s or 100s of points. This overall makes it difficult to truncate at a given point since it cannot be guaranteed that the data is being split at the correct start/end of the overlapping region - that is, without some hands on experimentation at least.


问:
some confusions about the implementation of spectral element method for incompressible flow
I am writing a parallel spectral element code of my own based on the three-step splitting scheme.

I have benchmark-tested my Poisson solver(with the first and second Neumann boundary conditions) and Helmholtz solver and both the solver can reach the machine error of 1E-15.
But when I implement the three-step splitting scheme for the Kovasznay flow, my converged error level can only reach an order of -5, I have tried to increase the N and NE, and decrease the time step,  but it doesn't work. It seems the pressure gradient obtained with my code is incorrect.
I have some confusion about the first advection step implementation now, could you please help me out?

In the first advection step, I simply calculate the u-star in a spectral collocation way,i.e. calculate its values on the mesh points,   but I see in some paper (in the following figure ),

that the first step needs to multiply the Weighting matrix W and the inverse matrix Mg
( I don't know what the matrix is ?). I don't understand what exactly does it mean? since the W can be canceled out from the left and right-hand sides. The following figure show my deduction about the implementation :

Could you please show me the detail about how to implement the first step.


问:
how should i find an appropriate time step for FSI modeling of a blood vessel?
i am working on a transient FSI solution with ansys fluent but i don't know how to set time step to get good results what's your suggestion about this guide me please</p>


问:
How to apply Dirichlet and Neumann boundary condition at the corner in Finite-Element-Method?
I have a 2D rectangular domain. The governing equation on this domain is Laplace equation:

∇2f=0

In the left and right edge  there is Neumann boundary conditon : ∂f∂n=a

n is the normal vector to the domain's boundary(here on the left edge it's equal to the negative direction of x axis) and 'a' is a given date and it's a constant.

There is a Dirichlet boundary condition at the bottom and upper  edge.

My problem is how to apply that bc at the corner nodes since there are two different bcs located at the same node.


问:
How to calculate a product of two real functions of large, but opposite magnitude?
I need to evaluate a product of two real functions, namely $F(x)\cdot G(x)$. The function $F(x)$ is a diverging and obscure member of scipy.special, while $G(x)$ is a gaussian. While the product is a well behaved function and of order ~1, its factors by themselves, for x range that I am interested in, are of order $10^{\pm 500}$. My problem is that each of them overflow pythons float. My go-to move is to try numpy.float128 which would cover range of interest quite nicely, if not for the error:

TypeError: ufunc 'pbdv' not supported for the input types, and the
inputs could not be safely coerced to any supported types according to
the casting rule ''safe''

So, function $F$ doesn't support float128. Is there a well known path to deal with this? Only thing I can think of right now is to express offending function as power series and calculate Cauchy product of sums - would that work?


问:
How to do a parametric study of arbitrary 2D surface?
I need to do a parametric study of the performance of a room heater on different rooms by simulating the temperature distribution in there. The problem here is that the rooms are not simple rectangles like shown in a) (in blue are the walls of the room and in red the heater), where I could just simply sweep over length and wide, but they may have more complex shapes like shown in b) to d).

I thought about making a raster and parametrize with 1 or 0 if the given rectangle in the raster is inside the walls, but the number of input variables would be too high and I would to make sure that the parametrization is consistent somehow.

Is there a systematic way to parametrize an arbitrary (limiting the number of nodes could be a possibility to limit complexity) 2D surface?


问:
Guaranteeing that the numerical approximation to $B = K^{-T}K^{-1}$ is positive definite in Zhang's algorithm
In Zhang's Algorithm to determine the intrinsic parameter matrix $K$ (see here for slides talking about this), we instead compute the matrix $B = K^{-T}K^{-1}$ since we obtain a linear system in the variables of that matrix. The general idea is that once we obtain $B$, we use the Cholesky decomposition to decompose $B$ into the form $LL^T$ and the uniqueness of this decomposition upto signs of the columns and the fact that $K$ is positive on the diagonal allows us to retrieve $K$.

However in practice, we cannot solve the system of equations exactly, and the suggested method is to take the least squares optimal solution of the linear system by picking the vector corresponding to the smallest eigenvalue.

My issue is how this guarantees that the $B$ we obtain is positive definite. Should we perform some sort of projection to ensure this is the case or will this just hold in general (or does it not matter)? The Cholesky decomposition requires that the matrix be at least positive semi-definite (and I belive uniqueness requires positive definite).

EDIT: I think there is some confusion in my wording, so I would like to be a bit more precise. The matrix $K$ is unknown and we wish to identify it. We instead look for the matrix $B = K^{-T}K^{-1}$ by solving a linear system and them applying the Cholesky decomposition to find $K$. However we don't solve the linear system exactly for $B$ (due to noise and real world errors) and we only obtain a (least-squares) approximation $\widehat{B}$ from which we calculate $K$. My question is why this matrix $\widehat{B}$ has to be positive definite.


问:
Looking for non-trivial examples of solutions to 3D wave equations?
We have developed a (new) numerical scheme to solve the classical wave equation in 3 dimensions and we aim to publish the results.

We can read in the aim and scope of the journal of computational and applied mathematics :

" ...The computational efficiency (e.g. the convergence, stability, accuracy, ...) should be proved and illustrated by nontrivial numerical examples"

Immediately we can raise the question : What are known (or classical) nontrivial examples of exact solutions of the 3D wave equation? by that I mean the Cauchy problem :
$$u_{tt}=c^2\Delta u,\quad u(X,0)= f(X),\quad u_t(X,0)= g(X)$$
Nontrivial examples against which we can compare the results of the numerical scheme, that is, we need given $f$ and $g$, along with a 'computable' formula for $u(X,t)= u(x,y,z,t)$


问:
Why does the non-orthogonal term in FVM need to be evaluated explicitly?
In FVM, the face normal gradient is split into orthogonal and non-orthogonal part. The orthogonal part is evaluated implicitly, whereas the non-orthogonal part is evaluated explicitly as a source term using known velocity field. My question is why can’t the non-orthogonal part be evaluated implicitly?</p>


问:
Rectangular problem to use for benchmarking a large dense linear solver?
Can anyone suggest a rectangular problem I can use for evaluating a large dense linear solver? ($m,n>1000$) IE, ideally something that's easy to download and has been evaluated for some solvers already.


问:
Anisotropic lines identification algorithm
I am implementing an algorithm to identify the anisotropic lines in an unstructured mesh. This is done in the framework of an Agglomeration algorithm for Non-Structured multigrid applied to FV discretization of N-S equations.

An Anisotropic line can be defined as a line of cells with the same characteristic: to have an aspect ratio that in case of a 2D mesh is more than 4 ($AR>4$) and in 3D more than 8 ($AR>8$). The lines are generated by connecting together the cells along the face of maximum area (or length in 2D). This is quite useful in the geometric multigrid algorithm to preserve the Boundary Layer characteristic structure (made of prism in 2D or of wedge/Hexa in the most classical non-structured meshes). Here after an example from the paper by Mavriplis on the subject.

In the literature I did not find many algoithms to accomplish the task. Right now I managed to implement the algorithm in this paper by Mavriplis (a reference for the community) that is a sort of greedy-front advancing algorithm. I report here the explanation of line identification algorithm given in the article:

The line construction
algorithm begins by pre-computing
the ratio of maximum
to average
adjacent
edge
weight for each vertex.The vertices are then sorted according
to this ratio.The first vertex in this ordered
list is then picked as the starting point for a line. The line is built by addingto the original vertex the
neighboring
vertex which is most strongly connected
to the current vertex,provided this vertex does not
already belong to a line,and provided
the ratio of maximum to minimum edge weights for the current vertex
is greater than 4 ...
The line terminates
when no additional vertex can be found. If
the originating
vertex is not a boundary point,then the procedure
must be repeated
beginning
at the original
vertex,and proceeding
with the second
strongest
connection
to this point.

The algorithm (even if it has been optimized) ends up to be pretty complex (in terms of branching) and does not guarantee in 3D efficient results, even if it is a good compromise among the final result quality (measured in the preservation of the Boundary layer of the NS mesh) and computational efficiency.

I wanted to ask if there are more efficient algorithms in literature (even from other communities, e.g computer vision community or graph community) that are more efficient and able to work globally (avoiding the front advancing local feature).


问:
Finding optimal values from multiple parameter estimation runs
I've performed a parameter estimation repeat (i.e. 1000 parallel runs with the same
initial values of parameters). I am trying to estimate ~20 parameters using measurements from experiments.

After running the estimation task (using Evolution strategy algorithm) for 10 hours (using --maxTime argument), the results were obtained.

When I look at the distribution of the objective values, the values fall in the range  (900, 2000).  There are several combinations of parameter values that I obtain for objective values that lie in the same order of magnitude. I understand increasing the duration of the estimation run may decrease the objective value further, but I think I would still end up with several sets of parameter values.

I would like to ask for suggestions on how to process these results.
If there are some references that I could read to understand how to cluster these results and find an optimal parameter set, that would be of great help.


问:
Solving a system of PDEs with an ODE
I want to solve the following system of equations which consists two PDEs and one ODE:
\begin{align}
    \rho_t+v\rho_x &= 0; \newline
    Y_t+vY_x &= 0 ;\newline
    v_t &= -\frac{1}{(\int_0^1\rho dx)}\Big[-\frac{1}{3}\frac{1}{Fr^2}\int_0^1(\rho-1) dx +\xi v|v|\Big]
\end{align}
where $\rho, Y, T$ are function of $(x,t)$ and $v$ is a function of $t$. And $\xi, Fr,v_0, a ,b,c,d$ are real constants. I am going to solve the above system with the following initial conditions. for $f$ given function
\begin{align}
        \rho(x,0) = f(x),\quad Y(x,0)=0,\quad v(0)=v_0
\end{align}
and following boundary conditions
\begin{align}
    \rho(0,t) &= 1 \quad \text{ if } v(0)>0\newline
    \rho(1,t) &= 1 \quad \text{ if } v(0)<0\newline
    Y(0,t) &= c \quad \text{ if } v(0)>0\newline
    Y(1,t) &= d \quad \text{ if } v(0)<0
\end{align}

It's clear that the first two advection equations and have analytical solutions but since the coefficient, $v$ depends on the time, I first try to solve the ODE for $v_t$ with Euler or Runge Kutta method and then solve the two PDE with the explicit upwind method. But still, I do not know if I'm solving it right. If anyone has a different idea of how to solve this system or has some comments on my approach I would really appreciate it.

Title: Solving a system of PDEs with an ODE


问:
How to obtain the unitary operator to get specific partial trace using searching methods?
Is there a unitary $U_{AB}$ such that, for any density operator $\rho$, we have
$${\rm {Tr}}_A \left[U_{AB} \left(\frac{I_A}{2} \otimes \rho_B\right)U_{AB}^{\dagger}\right]= \frac{\rho_B}{2}+\frac{I_B}{4}
\\
{\rm {Tr}}_B \left[U_{AB} \left(\frac{I_A}{2} \otimes \rho_B\right)U_{AB}^{\dagger}\right]= \frac{\rho_A}{2}+\frac{I_A}{4},$$
where $I_A=I_B=I$ is the identity matrix, $\rho_A=\rho_B=\rho$, A and B are both single qubit systems.

I have thought about decomposing $U_{AB}$ to express the partial trace but failed. I also considered about searching for such unitary, but $U_{AB}$ and $\rho$ are both unknown, which holds me back. Any ideas or comments, both in the analytic way or a computation way, would be appreciated.

Cross-posted and answered on qc.SE


问:
Integral from function approximations
I have some data which I cannot manage to model and fit with a known function, so let’s say that they are a sample from the unknow function $f(x)$, which look a sort of skewed bell-shaped distribution.

Given $a,b \in  \mathbb{R}$ and $\vec{w}$ some parameters of the unknown function:

$$\int_{a}^{b} f(x,\vec{w}=\vec{w_1}) = A \,,$$

I would like to know if there is an approximate method that allow me to sampling or to have a function approximator of $f(x,\vec{w}=\vec{w_2})$ s.t. $$\int_{a}^{b} f(x,\vec{w}=\vec{w_2}) = A^c$$ where $c \in     \mathbb{R}$.
What I'm interested in are the tails of distribution $f(x,\vec{w}=\vec{w_2})$.

Title: Integral from function approximations


问:
Modal representations of nodal tensor product Galerkin elements
Nodal discontinuous Galerkin methods on simplices, like those described in Hesthaven and Warburton, have the nice property that the number of nodes is equal to the minimum number needed to represent a unisolvent order $N$ polynomial in $d$ dimensions,

$$N_p = \frac{\prod_{i=1}^{d}(N+i)}{d!}$$

And we can construct $N_p$ unisolvent, hierarchical, orthonormal basis functions, and define an $N_p \times N_p$ Vandermonde matrix to uniquely map between nodal values and modal coefficients,

$$V \hat{u} = u$$

This is useful in many contexts:

Interpolation from solution nodes to other arbitrary nodes happens via a matrix multiplication like $u^\star = V^\star V^{-1} u$ where $()^\star$ is defined on the other nodes but with the same basis functions

P-multigrid restriction and prolongation matrices are just a matter of removing the highest modes or adding higher zero modes.

Spectral filtering operations are accomplished by mapping to modal space, truncating higher order coefficients, and mapping back to nodal space.

On tensor-product elements, we now have more nodes/coefficient than are required to be unisolvent, and we no longer have hierarchical basis functions (correct me if I'm wrong). How does one perform the above types of operations for tensor product elements? Do we form a Vandermonde matrix with the multidimensional monomials, and just not worry about the ill-conditioning and lack of orthogonality?


问:
Spurious oscillations in solving diffusion problems using finite elements
I've been struggling with this problem for a while so I hope someone can help me here.

I'm trying to solve the McNabb-Foster equations for hydrogen diffusion in metals using a simple 1D finite element simulation in Python. This works fine under isothermal conditions, but I can't make it work for anisothermic conditions (i.e. increasing temperature).
The equations are essentially a non-linear Poisson equation:
$$
\frac{dC}{dt}=D ∇^2 C-N\frac{dθ}{dt}
$$
where
$$
\frac{dθ}{dt}=κC(1-θ)-λθ
$$

Where C is the hydrogen concentration, D, κ and λ are temperature-dependent terms, N is constant and θ is the so-called trap occupancy (traps are sites with high binding energy). Using the backward difference method the first equation can be expressed in matrix form as:
$$
(\frac{1}{∆t}[M]+D[K])\{C\}_{t+∆t}={F}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]N\Bigl\{\frac{dθ}{dt}\Bigl\}_{t+∆t}
$$
Where [M] and [K] are the system matrices. So far, so good, but the tricky part is how to incorporate θ. For the isothermal case I have done this again using the backward difference method, but taking θ from the previous iterative step:
$$
\Bigl\{\frac{dθ}{dt}\Bigl\}_{t+∆t}=\frac{\{θ\}_{t+∆t}-\{θ\}_t}{Δt}=κ\{C\}_{t+∆t}\{1-θ\}_t-λ\{θ\}_t
$$
Combining the above two equations gives:
$$
(\frac{1}{∆t}[M]+D[K]+[M]Nκ\{1-θ\}_t )\{C\}_{t+∆t}=\{F\}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]Nλ\{θ\}_t
$$
Below is the iterative part of my code where I have tried to implement this. Basically, I'm solving the equations twice; once to find the new concentrations C and a second time to find the flux F, which is actually what I need.

The trouble is that I get (spurious?) oscillations at the boundary for reasonable input parameters (see image). These disappear for very small time increments, but it makes the code too slow (I'm running it in a least-squares routine so it needs to be fast). The mesh size has no influence.

Does anyone have any suggestions how I can improve the stability?

```<code>for i in range(1, Iter+1):

     Clold = Cl # save Cl from last iteration
     T = Tstart + i * dt * dTdt # calculate temperature

     # calculate temperature dependent terms
     Dl = D0 * np.exp(-El / R / T)
     kappa = kappa0 * np.exp(-Et / R / T)
     lambd = lambda0 * np.exp(-Ed / R / T)

     # establish reverse difference equation
     LHS = sysm + dt * Dl * sysk + dt * Nt * kappa * sysm * (1-theta)
     RHS = np.dot(sysm, Clold) + dt * np.dot(sysm, (Nt * lambd * theta))

     # apply constraints & boundary conditions
     LHS[0, :] = 0.0
     LHS[0, 0] = 1.0
     RHS[0, 0] = 0.0

     # calculate new Cl
     Cl = np.linalg.solve(LHS, RHS)

     # calculate flux matrix
     LHS = sysm + dt * Dl * sysk + dt * Nt * kappa * sysm * (1-theta)
     RHS = np.dot(sysm, Clold) + dt * np.dot(sysm, (Nt * lambd * theta))
     F = np.dot(LHS, Cl) - RHS

     # calculate new trap occupancies
     theta = theta + dt * (kappa * Cl * (1 - theta) - lambd * theta)

     # store flux for plotting
     Fstore = np.append(Fstore, np.array(-F[0] / dt))
</code>```

EDIT: using an implicit scheme the last equation would look like this:

$$
(\frac{1}{∆t}[M]+D[K]+[M]Nκ\{1-θ\}_{t+∆t} )\{C\}_{t+∆t}=\{F\}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]Nλ\{θ\}_{t+∆t}
$$

So now you have a second unknown term in the equation. The new θ can be given in terms of the new C, i.e.:

$$
θ_{t+∆t}=\frac{∆tκC_{t+∆t}+θ_t}{1+∆tκC_{t+∆t}+∆tλ}
$$

but this makes the first equation non-linear. Any suggestion how to solve this?


问:
Unstable FDTD for axisymmetrical wave propagation
I am writing a 2.5D axisymmetric FDTD solver for Maxwell equations. It operates on three fields on staggered grids ($E_r$, $E_z$, $B_\varphi$).

As a test problem, I am considering a cylindrical waveguide of finite length (open at each side), excited by an axial current density $J_z$ (Gaussian pulse). A wave starts propagating towards the axis, and after it reaches it, the algorithm becomes unstable.

It is probably related to the way I remove the singularity of Maxwell equations in cylindrical coordinates (the approach taken from [1] and [2]).

The $E_z(r=0)$ is calculated from $\oint_C H_\varphi dl = \int_S \partial/\partial t (\varepsilon E_z) dS$, where $H_\varphi = B_\varphi/ \mu_0$, and contour $C$ is a circle with radius $\Delta r/2$ and area $S = \pi (\Delta r/2)^2$.
This leads to the following update equation on axis:
E_z[1,k] = E_z[1,k] + c^2 * DT * 4.0 / DR * B_φ[1,k]

However, this leads to the instability when the wave reaches the axis. If I use 2.0 / DR instead of 4.0 / DR, the instability is gone.

I've checked the source code of other FDTD solvers and the all use 4.0 / DR (for example [3]).

I've also prepared a minimal code reproducing the instability (it's written in Julia): https://gist.github.com/bchaber/a354f8fdedd58dbca9a1712bc69fe739

[1] Hanbing Jin, "Particle-in-Cell Simulation of Electromagnetic Pulse Generated by High-power Laser-target Interaction", p. 22, URL: https://kth.diva-portal.org/smash/get/diva2:1224073/FULLTEXT01.pdf

[2] Yu M Tolochkevych et al 2014 J. Phys., p. 3, URL: https://iopscience.iop.org/article/10.1088/1742-6596/511/1/012001/pdf

[3] https://github.com/ECP-WarpX/WarpX/blob/5761b4bf998eab84ea0f7e4b132026593f3ddf9f/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp#L355


问:
Simplest way to "monotonize" MacCromack method
The MacCormack finite-difference "predictor-corrector" method is well known to generate spurious oscillations near solution discontinuities such as shock waves in gas dynamics equations. Or even in simpler case of scalar advection equation, when advected function contains jumps.

For teaching purposes, I'd like to find a simplest way to make it work with such discontinuities, so that it wouldn't be too complex for students to implement. I'm aware of such methods as adding Davis artificial viscosity or applying Boris-and-Book-like "monotonizators". But are there simpler ways to dampen oscillations to some acceptable values?