问:
Numerical integration of a rapidly varying complex exponential
I have a function $f : \mathbb{R}^2 \mapsto \mathbb{R}^+$ and I wish to numerically evaluate the integral below over a finite domain $\Omega \subset \mathbb{R}^2$

$$
I = \int_\Omega e^{i k \cdot f(\mathbf{x})} \mathrm{d} \mathbf{x}
$$

where $k \in \mathbb{R}$ is a large positive real constant, for example typically $k \approx 10^9$. The function $f$ is bounded over the domain of integration, and is continuous almost everywhere but does have some discontinuities. It is pretty smooth (qualitatively speaking) with gentle variation over its domain and no sharp peaks. In my particular case I know that

$$
0 < f(\mathbf{x}) \leq 0.1
$$

I write it this way to show that there is not much variation in the function itself, and the size of the $k$ parameter is the issue.

I am also able to obtain exact values for $f$ for various points in its domain, but there is no closed-form expression for it.

I am struggling to find a viable way to numerically evaluate $I$. The best way I found was to break up the 2D domain of integration into a mesh of small triangles, evaluate the function at each vertex and integrate piecewise using a linear approximation of the function over each triangle using barycentric interpolation. This works really well for moderate $k$ but as $k$ becomes large it starts requiring far too many triangles to get a good result. I suspect this is because the linear approximation needs to be extremely accurate for it to "survive" the complex exponential, and of course the error increases as $k$ grows.

Are there known methods of numerical integration that "work well" with problems of this form, or are they inherently difficult? What are some things I could try?


问:
2D wave equation is numerically unstable using Finite Difference Method
I'm working with simulating both the heat and wave equation in 2D in a Python code. When simulating the heat equation, I learned about the CFL which I used to get a numerical stable solution.

I found CFL for the 2D wave equation as well, but my solution keeps being numerically unstable no matter how I change the time increments $dt$. I'm presenting the code here:

# Libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib as mpl

def f(X,Y): # f(x,y,0)
    return X**2+Y**2

def g(X,Y): # g(x,y,0)
    return X

def laplacian(arr, row,col, dx ):
    return (arr[row + 1, col]
                            + arr[row - 1, col] + arr[row, col + 1]
                            + arr[row, col - 1]- 4*arr[row,col])/(dx**2)

def simulation(rect,hs,BC, c, frames, eq, eps = 1e-10):

    """ Simulation of the heat equation on a
    2D rectangular grid using the Finite Difference Method """

    X, Y = np.meshgrid(np.linspace(rect[0],rect[1],hs[0]),
                       np.linspace(rect[2],rect[3],hs[1])) # 2D meshgrid

    # Initial Conditions
    Z_init = f(X,Y) # u(x,y,0) = f(x,y,0)
    Z_dot_init = g(X,Y) # u_t(x,y,0) = g(x,y,0)

    zmax = max(Z_init.max(), BC[0], BC[1], BC[2], BC[3])
    zmin = min(Z_init.min(), BC[0], BC[1], BC[2], BC[3])
    zs = [Z_init]

    # Boundary Conditions
    if eq == "Heat":
        Z_init[0] = np.ones(len(Z_init[0])) * BC[0]
        Z_init[-1] = np.ones(len(Z_init[-1]))* BC[1]
        Z_init[:,0] =  np.ones(len(Z_init[:,0]))* BC[2]
        Z_init[:, -1] = np.ones(len(Z_init[:, -1])) * BC[3]

    # Figure settings
    fig = plt.figure()
    ax = plt.axes(projection='3d')
    ax.axes.set_xlim3d(rect[0] + eps, rect[1] - eps)
    ax.axes.set_ylim3d(rect[2] + eps, rect[3] - eps)
    ax.axes.set_zlim3d(zmin, zmax)
    plt.rcParams['mathtext.fontset'] = 'stix'
    plt.rcParams['font.family'] = 'STIXGeneral'

    # Plots for lines constituting the rectangle
    ax.plot([rect[0], rect[1]], [rect[2], rect[2]], [BC[0], BC[0]], color='black', linewidth=2)
    ax.plot([rect[1], rect[1]], [rect[2], rect[3]], [BC[3], BC[3]], color='black', linewidth=2)
    ax.plot([rect[0], rect[0]], [rect[2], rect[3]], [BC[2], BC[2]], color='black', linewidth=2)
    ax.plot([rect[0], rect[1]], [rect[3], rect[3]], [BC[1], BC[1]], color='black', linewidth=2)

    ax.plot([rect[0], rect[0]], [rect[2], rect[2]], [BC[0], BC[2]], color='black', linewidth=2)
    ax.plot([rect[1], rect[1]], [rect[2], rect[2]], [BC[0], BC[3]], color='black', linewidth=2)
    ax.plot([rect[0], rect[0]], [rect[3], rect[3]], [BC[1], BC[2]], color='black', linewidth=2)
    ax.plot([rect[1], rect[1]], [rect[3], rect[3]], [BC[1], BC[3]], color='black', linewidth=2)
    ax.set_title('Temperature development in a \n'
                 'rectangular room', fontsize = 18, fontname = 'STIXGeneral')

    # Infitesimals
    dx = (rect[1] - rect[0]) / ((hs[0] - 1))
    if eq == 'Heat':
        dt = dx**2/(10*c**2)
    else:
        dt = dx/(10*c)

    surf = ax.plot_surface(X, Y, Z_init, alpha=0.7, cmap=plt.cm.jet, vmin=zmin, vmax=zmax)
    cbar = fig.colorbar(surf)
    cbar.ax.set_ylabel('Temperature in ' + r'$^\circ\mathrm{C}$', rotation=270,fontsize = 14, labelpad=20)
    ax.set_axis_off()

    # Finite Difference Method
    if eq == 'Wave':

        for row in range(1, hs[0] - 1):
            for col in range(1, hs[1] - 1):
                Z_init[row,col] = Z_init[row,col] - 2 * laplacian(Z_dot_init,row,col, dx) + 1/2 * c ** 2 * dt ** 2 / (
                    dx) ** 2 * (Z_init[row+1, col] - 4 * Z_init[row, col]
                                + Z_init[row-1, col] + Z_init[row, col+1]
                                + Z_init[row, col-1])

        zs.append(Z_init)

        for iteration in range(2,frames):
            surf.remove()
            for row in range(1, hs[0]-1):
                for col in range(1,hs[1]-1):
                    Z_init[row,col] += 2 * zs[iteration - 1][row, col] - zs[iteration - 2][row, col] + c ** 2 * dt ** 2 / (
                            dx**2) * (zs[iteration - 1][row + 1, col] - 4 * zs[iteration - 1][row, col]
                                        + zs[iteration - 1][row - 1, col] + zs[iteration - 1][row, col + 1]
                                        + zs[iteration - 1][row, col - 1])

            zs.append(Z_init)
            surf = ax.plot_surface(X, Y, zs[iteration-2], alpha = 0.7, cmap =plt.cm.jet, vmin = zmin, vmax = zmax)
            plt.pause(dt * 10)

    if eq == 'Heat':

        for iteration in range(frames):

            surf.remove()
            for row in range(1, hs[0]-1):
                for col in range(1,hs[1]-1):
                    Z_init[row,col] = Z_init[row,col] + c**2 * dt / (dx)**2  *(Z_init[row+1,col] - 4*Z_init[row,col]
                                                   + Z_init[row-1,col] + Z_init[row,col+1]
                                                   + Z_init[row,col-1])

            # Plotting surface with slight pause
            surf = ax.plot_surface(X, Y, Z_init, alpha = 0.7, cmap =plt.cm.jet, vmin = zmin, vmax = zmax)
            plt.pause(dt)

simulation(rect = [-1,1,-1,1], hs = [100,100], BC = [0,0,0,0], c = 3, frames = 1000, eq = 'Wave')
```

I'm adding an image of the user instructions too:

Either way, I've tried searching for how to make my code numerically stable, but with no results. I tried changing the time steps to be even lower, I tried adding more elements into the grid, tried changing the maximum velcity, tried changing the initial conditions. I even tried removing the boundary conditions to only apply for the heat equation.

If you have any ideas on how to help me solve this, it'd be greatly appreciated.

ContentLicense: CC BY-SA 4.0


问:
Consering numerical implementation of gradient based method for control system
I'm trying to reproduce the results in Optimal consensus control of the Cucker-Smale model by Bailo et al. The system is the following,

the adjoint variables,

and the algorithm,

I tried to implement the dynamics in the following manner

def f(X, t, u):
    N = int(len(X)/4)
    x = X[:N]
    v = X[N:2*N]
    p = X[2*N:3*N]
    q = X[3*N:]

    dxdt = v
    dvdt = np.zeros(N)

    dpdt = np.zeros(N)
    dqdt = np.zeros(N)

    v_bar = np.mean(v)

    for i in range(N):
        for j in range(N):
            norm_xij = np.linalg.norm(x[i] - x[j])
            if norm_xij > 0:
                dvdt[i] += a(np.linalg.norm(norm_xij))*(v[j]-v[i])

                dpdt[i] += (a_prime(norm_xij)/norm_xij) * np.inner(q[j] - q[i], v[j] - v[i]) * (x[j] - x[i])
                dqdt[i] += a(norm_xij) * (q[j] - q[i])

        dvdt[i] *= 1/N
        dvdt[i] += u[i]

        dpdt[i] *= 1/N
        dqdt[i] *= 1/N

        dqdt[i] += p[i] - (2/N) * (v_bar - v[i])

        dpdt[i] *= -1
        dqdt[i] *= -1

    return np.concatenate([dxdt, dvdt, dpdt, dqdt])

```

I have 2 questions regarding the implementation ?

For the control $u$ , should´t it be an array (for each $u_i$) ? How can I implement this ?

Is there any way to implement steps 1 and 2 from the algorithm 1 in two separate functon ?


问:
Efficiently compute the Fourier series coefficients of a piecewise trigonometric function
I am searching for the most efficient algorithm to compute the Fourier series coefficients of a number of given functions $f(\theta)$
$$a_k = \int_0^{2\pi} f(\theta)\cos(k\theta) \\
b_k = \int_0^{2\pi} f(\theta)\sin(k\theta)$$
for $k \in {1,...,N}$. The functions $f(\theta)$ are continuous piecewise trigonometric functions, i.e. in each interval $l$ they can be described as follows:
$$f(\theta) = c_{1,l}\cos(\omega\theta) + c_{2,l}\sin(\omega\theta)$$
The functions $f(\theta)$ are also known to be continuous and periodic at $2\pi$. Note that the value of $\omega$ is the same for all sub-intervals.

The goal is to compute $a_k,b_k$ with fastest speed having known the values of $c_{1,l},c_{2,l}$ and $\omega$. Note that there are many number of $f(\theta)$ functions which require great computational efficiency. I know that there is closed form solution for this, but it going to be a trigonometric function itself which at high number of evaluations it becomes costly. One approach that came to mind was using FFT, but then it will require evaluation of $f(\theta)$ at a high number of discrete points which itself requires a lot of trigonometric evaluations. For now, I use the closed form formula. For a 100 number of different $f(\theta)$s each having 100 number of intervals for $N=100$, it takes approximately 1 second. I need to reduce the time to an order of milliseconds if possible.


问:
Is it possible to globally solve a general coordinate-wise monotone nonlinear system 3x3?
Consider a system
$$
F_i(x, y, z) = 0, \quad i = 1, 2, 3
$$
with $F_i$ monotonic w.r.t. $x, y$ and $z$.

The system 2x2 can be easily solved with alternating direction method that will find all its roots.

The problem is to find all roots of the 3x3 system with monotonic coordinate functions. Is it possible in a general case? In other words, it is necessary to find all intersections of 3 surfaces with positive coordinates of normal vectors in 3D space.

Illustration 1: linear case, 1 root

Illustration 2: nonlinear case, 2 roots


问:
Translating poisson equation to sfepy code
As the title suggests, I am trying to convert the Poisson equation:

Into sfepy code that can solve the differential euqation problem numerically.

I am trying to modify the tutorial on lienar elasticity which is similar but I am just causing it to crash.


问:
Numerov Method with Time Varying Potential
Is it possible to use the Numerov method to solve the Time Dependent Schrodinger Equation ($\frac{i\partial\Psi(x, y, z, t)}{\partial t} = \nabla^2\Psi(x, y, z, t) + \Psi(x, y, z, t)V(x, y, z, t)$) with a time varying potential $V(x, y, z, t)$? If so, how?

Thanks!


问:
Can the Crank-Nicolson Method Be used to Solve The Schrodinger Equation with a Time Varying Potential?
I have been following an excellent article about how to use the Crank-Nicolson method to solve the Schrodinger equation. In the article, it starts with a $V(x, y, t)$ but the potential seems to become constant later on. Am I able to recompute the potential at a new time value, take, the previous wave function, and plug it into a new set of matrices ($A$, and $M$ in the article) for the next step?

Additional context: I see in the post there are terms
$$V(x, y, z)\psi(x, y, z) = \frac{1}{2}\left(V_{ij}^{n+1}\psi_{ij}^{n+1} + V_{ij}^{n}\psi_{ij}^{n}\right)$$
$$a_{ij} = \left(1 + 2r_x + 2r_y + i\frac{\Delta t}{2}V_{ij}^{n+1}\right)$$
$$b_{ij} = \left(1 + 2r_x + 2r_y + i\frac{\Delta t}{2}V_{ij}^{n}\right)$$
Where $b_{ij}$ and $a_{ij}$ are used each step to solve $\psi(x, y, t)$, this is for a 2D scenario in the post so I think $z$ in the first equation may have been intended to be $t$ as $n$ is the time index and the term $\psi(x, y, t)V(x, y, t)$ appears earlier in the post. I think $a$ and $b$ may not have a time index because the only term in them that depends on time is the potential, but if you don't have a time-varying potential, they would remain the same values throughout each step.  So here, the author I think, excludes it because they use a constant potential ($V_{ij}^{n} = V_{ij}^{n+1}$). But if you did not have a constant potential, could you simply put $n$ it back onto $a$ and $b$ and recompute $a_{ij}^{n+1}$ and $b_{ij}^{n}$ each step? And if so, what are the limits, are there limits to how the potential can change from step to step?


问:
discrete definition of curl $ \nabla \times \vec{A}$ on a 3D grid?
I have the data for 3D vector field $\vec{A}$ (with components $\vec{A_1}$, $\vec{A_2}$ and $\vec{A_3}$) sampled on a 3D grid with integer indices i, j and k.

Assuming that only the third component $\vec{A_3}$ is non-zero, I have the expression for the components of curl: $\vec{B} = \nabla \times \vec{A}$ as:

```<code>B1 = -(A3[i][j][k] - A3[i][j + 1][k] + A3[i + 1][j][k] - A3[i + 1][j + 1][k])/ (2 * dx[2])

B2 = (A3[i][j][k] + A3[i][j + 1][k] - A3[i + 1][j][k] - A3[i + 1][j +1][k]) / (2 * dx[1])

B3 = 0
</code>```

How can the above definition be expanded for a vector $\vec{A}$ with all three components non-zero on a full 3D grid?

I have seen the answer (discrete definitions of curl $\nabla \times F$?) for a similar question on a 2D grid, but I am not sure how to expand it for a 3D grid.


问:
Is it possible to realign sliding-window FFT-filtered output on time-line with the original signal when it is used for noise-filtering in real-time?
In real-time systems, an output data can be taken into sliding-window Fourier Transformation directly without waiting for new data. Then Fourier Transformed data can be altered (remove high frequencies, etc) before Inverse Fourier Transformation, to get a low-pass filter that works for all data except for the first N data points.

One of the problems of this filtering is that the filtered signal follows the real signal with a phase that depends on window size. So, would it be possible to adjust phases of each frequency domain element to make the IFFT produce some signal closer to real signal (but with less noise) on time-line?

What if we generate N more elements from the truncated (for filtering) frequency array (before IFFT) to use it as prediction to do a 2x sized FFT+truncate+IFFT where middle point is the last real-time-signal? How would it work? Would it have disadvantages for some signal / noise types?


问:
Crank Nicolson Simulation Not Preserving Probability?
I have written a Crank-Nicolson simulation based on this post and the code it links too.

def placeDiagonals(
                center : np.array,
                innerShifted : np.array,
                inner : np.array,
                outer : np.array,
                rowLength : int,
                extent : int,
                math,
                sparse
            ):
        return sparse.spdiags(
                math.array([
                        outer,
                        inner,
                        center,
                        innerShifted,
                        outer
                    ]),
                math.array([-extent, -1, 0, 1, extent]),
                rowLength,
                rowLength
            )

def createStepMatrix(
            simulator,
            currentPotential,
            scalar,
            shiftedInnerDiagonal,
            innerDiagonal,
            outerDiagonal
        ):
    math = simulator.math
    sparse = simulator.sparse
    unknownFactorCount = simulator.unknownFactorCount
    centerDiagonal = (scalar * 1j * simulator.timeStep / 2) * currentPotential[1:-1, 1:-1].ravel()
    centerDiagonal = 1 + math.sum(scalar * simulator.stepConstants * 2.0) + centerDiagonal
    stepMatrix = placeDiagonals(
            centerDiagonal,
            shiftedInnerDiagonal,
            innerDiagonal,
            outerDiagonal,
            unknownFactorCount,
            len(currentPotential) - 2,
            math,
            sparse
        )
    return stepMatrix

def createCurrentStepMatrix(simulator, currentPotential):
    return createStepMatrix(
            simulator,
            currentPotential,
            1,
            simulator.knownInnerDiagonalShifted,
            simulator.knownInnerDiagonal,
            simulator.knownOuterDiagonal
        )

def createNextStepMatrix(simulator, nextPotential):
    return createStepMatrix(
            simulator,
            nextPotential,
            -1,
            simulator.unknownInnerDiagonalShifted,
            simulator.unknownInnerDiagonal,
            simulator.unknownOuterDiagonal
        )

class Simulator:
    def __init__(self, profile : SimulationProfile):
        self.profile = profile
        self.grid = self.profile.grid
        self.dimensions = self.profile.dimensions
        self.timeStep = self.profile.timeStep
        self.spaceStep = self.profile.spaceStep
        self.math = self.profile.math
        self.sparse = self.profile.sparse
        self.linalg = self.profile.linalg
        self.stepConstants = self.math.ones(self.dimensions) \
                * (-self.timeStep / (2j * self.profile.spaceStep ** 2))
        self.waveFunctions = [
                applyEdge(self.profile.initialWaveFunctionGenerator(self.grid))
            ]
        self.potentials = [
                self.profile.potentialGenerator(self.grid, 0.0),
                self.profile.potentialGenerator(self.grid, self.timeStep)
            ]
        self.unknownFactorCount = (self.grid.pointCount - 2) ** self.dimensions
        self.diagonalLength = self.unknownFactorCount
        self.knownInnerDiagonal = self.math.ones(self.diagonalLength) \
                * self.stepConstants[0]
        if self.profile.edgeBound == True:
            extent = int(self.math.sqrt(self.diagonalLength))
            self.knownInnerDiagonal.reshape((extent, extent)).T[-1] = 0
            self.knownInnerDiagonal = self.knownInnerDiagonal.reshape(self.diagonalLength)
        self.knownInnerDiagonalShifted = self.math.roll(self.knownInnerDiagonal, 1)
        self.knownOuterDiagonal = self.math.ones(self.diagonalLength) \
                * self.stepConstants[1]
        self.unknownInnerDiagonal = -1 * self.knownInnerDiagonal # These two just to save a bit of compute
        self.unknownOuterDiagonal = -1 * self.knownOuterDiagonal # time multiplying these large arrays by -1
        self.unknownInnerDiagonalShifted = -1 * self.knownInnerDiagonalShifted

    def compute(self, time, unknownStepMatrix, knownStepMatrix):
        math = self.math
        sparse = self.sparse
        waveFunctionVector = self.waveFunctions[-1][1:-1, 1:-1].reshape(
                (self.diagonalLength, 1)
            )
        independantTerms = knownStepMatrix @ waveFunctionVector
        #independantTerms = math.matmul(knownStepMatrix, waveFunctionVector)
        nextWaveFunction = self.linalg.spsolve(
                    sparse.csr_matrix(unknownStepMatrix),
                    independantTerms
            ).reshape(tuple([self.grid.pointCount - 2] * self.dimensions),)
        self.waveFunctions.append(math.pad(nextWaveFunction, 1))
        return self.waveFunctions[-1], self.potentials[-1]

    def simulate(self, timePoints : int, printProgress : bool = False):
        math = self.math
        timeStep = self.timeStep
        #timePoints = round(maxTime / timeStep)
        time = timeStep
        for ii in range(1, timePoints):
            self.potentials.append(
                    self.profile.potentialGenerator(self.grid, (ii + 1) * self.timeStep)
                )
            knownStepMatrix = createCurrentStepMatrix(self, self.potentials[-2])
            unknownStepMatrix = createNextStepMatrix(self, self.potentials[-1])
            self.compute(time, unknownStepMatrix, knownStepMatrix)
            time = (ii + 1) * timeStep
            if printProgress == True:
                print("(" + str(ii) + "/" + str(timePoints) + ")")

    def processProbabilities(self):
        math = self.math
        self.probabilities = math.array(list(map(
                lambda waveFunction :  math.sqrt( \
                        math.real(waveFunction) ** 2 \
                        + math.imag(waveFunction) ** 2 \
                    ).astype(math.float64),
                self.waveFunctions
            )))
        self.probabilityDecibles = 0
        #self.probabilityDecibles = np.array(list(map(
        #        lambda probabilities : -10 * math.log10(probabilities),
        #        self.probabilities
        #    )))
        return self.probabilities, self.probabilityDecibles

def makeWavePacket(grid, startX, startY, sigma=0.5, k = 15 * np.pi, math = np):
    unnormalized = math.exp((-1 / (2 * sigma ** 2)) * ((grid.x - startX) ** 2 + (grid.y - startY) ** 2)) \
            * math.exp(-1j * k * (grid.x - startX))
    totalProbability = math.sum(math.sqrt(math.real(unnormalized) ** 2 + math.imag(unnormalized) ** 2))
    return unnormalized / totalProbability
```
<p>^^ Most of that generates a diagonalized Hamiltonian of the form (with no "perfect" boundaries):</p>
```[[A00 r_x ___ r_y]
 [r_x A01 r_x ___]
 [___ r_x A02 r_x]
 [r_y ___ r_x A03]]
```
<p>I call it like this</p>
<pre class="lang-py prettyprint-override"><code>pointCount : int = 50
spatialStep : float = 1e-2
temporalStep : float = (spatialStep ** 2) / 4
courantNumber : float = 1.0
assert (temporalStep / spatialStep) <= courantNumber
profile = SimulationProfile(
        makeLinspaceGrid(pointCount, 1, 2, False, float, np),
        lambda position : makeWavePacket(position, .8, 0.5, .05, 15 * np.pi, np),
        lambda position, time : position.x * 0.0, #doubleSlit(position, .4, .1, .05, .0125, 20000),
        temporalStep,
        spatialStep,
        False,
        True
    )
simulator = Simulator(profile)
simulator.simulate(1500, False)
probabilities, probabilityDecibles = simulator.processProbabilities()
</code>```

Here is the question, I noticed that the sum of all probabilities is not 1, and it seems to increase then oscillate around a certain term.

Total Probability over all spatial coordinates Vs. Time:

The only value close to 1 (within 4 digits of precision) is the first one, which I normalized (this graph is scaled up 10x to 100x if I dont, but retains the shape).

I was reading this previous post an answer by one user states that the "most obvious method" method is not unitary in time and therefore does not conserve probability, seems to basically state C-N eliminates this problem.

So I am confused, is there something about Crank-Nicolson that I am missing that would make probability not be conserved? Or is it my code, or my boundary conditions? Do I have to normalize each frame? How can I conserve probability?

To summarize: Is it because of C-N method that probability is not being conserved in my simulation, or is it something to do with something I did, my code or my conditions, etc.?

Thank you!


问:
Toggling Constraints in Mixed Integer Programming
Are there MIP solvers that allow certain constraints to be toggled based on the value of a binary variable? My current situation is that I'm approximating the desired behavior by using constraints of the form $$A_{i}x\leq b_i + C\delta_{i}$$ where $x$ and $\delta$ are decision variables, $\delta_{i}\in{0,1}$, and $C$ is a large constant. Thus, constraint $i$ is relaxed when binary variable $\delta_i=1.$ This solution isn't ideal because there is no guarantee that the constant $C$ is large enough to relax the constraints. Moreover, there appear to be numerical issues with the Gurobi solver when choosing too large a value for $C$.

I am not well-versed in how MIPs are solved. Is the solution method amenable to an approach where the constraints are entirely removed from the problem whenever $\delta_i=1$?

As I mentioned, I'm currently using Gurobi as a solver. I'm mainly interested in keeping the constraints affine because this leads to my problem being a mixed integer linear program.


问:
System of ODEs with stochastic process
Consider system of ordinary differential equations

$$
\begin{align}
\frac{dx}{dt}&=\sigma(y-x)\\
\frac{dy}{dt}&=\rho x-y-xz\\
\frac{dz}{dt}&=xy-\beta z
\end{align}
$$

with $\sigma=10$, $\rho=28$ and $\beta=\frac{8}{3}$ and with initial condition $(x_0,y_0,z_0)=(1.508870, -1.531271, 25.246091)$

I am trying to embed the system in a stochastic environment by adding an exponentially colored noise process to the parameter $\sigma$: $$d\tilde{\sigma}_t=-\alpha\tilde{\sigma}_tdt+\gamma dW_t,\sigma_t=10+\tilde{\sigma}_t$$

I got following system of stochastic differential equations

$$
\begin{align}
dx_t&=(10+\tilde{\sigma}_t)(y_t-x_t)dt\\
dy_t&=(\rho x_t-y_t-x_tz_t)dt\\
dz_t&=(x_ty_t-\beta z_t)dt\\
d\tilde{\sigma}_t&=-\alpha\tilde{\sigma}_tdt+\gamma dW_t
\end{align}
$$

and I want to analyze the strong and weak order of convergence of the Euler-Maruyama scheme. Does the stochastic system have exact solution? If not, what is the best way of acquiring accurate solution that I can use to calculate errors agains numerical approximations in Python?

I was thinking about generating a path of $\tilde{\sigma}_t$ first and using scipy.integrate.solve_ivp to solve the rest of ODEs. However, it seems like time step used in solve_ivp is adjusted and t passed to fun is not integer so I am not sure how to combine pregenerated path with the rest of ODEs.

Another approach I had in mind was to solve the system of SDEs as is but I am not aware of Python libraries that can solve the system of SDEs and I believe that I need to use the same path for both the reference solution and the numerical approximation when calculating errors.

Title: System of ODEs with stochastic process


问:
Why GPU scaling (speedup) is usually worse than CPU scaling?
Let's define scaling, as how linear the speedup of using more than one GPU or CPU is. For example, having 2 GPUs gives you 2x faster execution time.

I have noticed that in many software (e.g. Molecular Dynamics, Gromacs, AMBER, OPENMM), it scales very badly for more than 2 GPUs, but it is not a problem for CPUs (you can use hundreds of CPUs, even located on different nodes). What is stopping GPUs from achieving the same scaling as CPUs? Why CPUs even located physically far away with huge latency, can achieve better scaling than GPUs?


问:
Integration problem
I want to numerically solve integrals of the form,
$$
I = \int_a^b x^k f(x) dx
$$
where $k$ is a given integer, and $f(x)$ is a cubic polynomial, expressed as,
$$
f(x) = c_0 + c_1 (x - a) + \frac{c_2}{2} (x-a)^2 + \frac{c_3}{6} (x-a)^3
$$
When $k \geq 0$, the integrand is a polynomial of degree $k+3$ and a standard Gauss-Legendre approach could be used, however this doesn't seem appropriate for the $k < 0$ case. I have tried two methods so far. I need to calculate these integrals for many different limit pairs $(a,b)$ and for all integers $k \in [-100,100]$.

Method 1

Directly integrating term by term. Essentially I used Mathematica to find the analytic expressions for the integrals of each of the 4 terms. There are multiple cases to be considered (e.g. when $k=-1,-2,-3,-4$ and otherwise). This leads to a large function where each case involves lots of adding,subtracting and multiplying. I get an answer which is close to correct, but it seems I am losing signficant digits somewhere in the complicated expressions for each term, perhaps due to cancellation problems.

Method 2

For this method, I attempted to derive a simple quadrature rule (simlar to Gauss-Legendre). We proceed by assuming,
$$
\int_a^b x^k f(x) dx = w_1 f(x_1) + w_2 f(x_2)
$$
with the weights $w_1,w_2$ and nodes $x_1,x_2$ to be determined. Writing out the equations for the cases $f(x) = 1,x,x^2,x^3$ gives the following nonlinear system of equations.
$$
\begin{align}
w_1 + w_2 &= \left( \begin{array}{cc} \frac{1}{k+1} (b^{k+1} - a^{k+1}), & k \neq -1 \\\\ \log{\frac{b}{a}}, & k = -1 \end{array} \right. \\\\
w_1 x_1 + w_2 x_2 &= \left( \begin{array}{cc} \frac{1}{k+2} (b^{k+2} - a^{k+2}), & k \neq -2 \\\\ \log{\frac{b}{a}}, & k = -2 \end{array} \right. \\\\
w_1 x_1^2 + w_2 x_2^2 &= \left( \begin{array}{cc} \frac{1}{k+3} (b^{k+3} - a^{k+3}), & k \neq -3 \\\\ \log{\frac{b}{a}}, & k = -3 \end{array} \right. \\\\
w_1 x_1^3 + w_2 x_2^3 &= \left( \begin{array}{cc} \frac{1}{k+4} (b^{k+4} - a^{k+4}), & k \neq -4 \\\\ \log{\frac{b}{a}}, & k = -4 \end{array} \right.
\end{align}
$$
This system of equations can be solved using a root-finding method. In my experiments, this method provides a more accurate solution than Method 1. However the drawback here is that the equations need to be solved for every pair $(a,b)$ (unlike Gauss-Legendre which can be solved on $[-1,1]$ and then scaled to any $[a,b]$). Also, I had to use a Monte-Carlo approach to find appropriate initial conditions for the weights and nodes in order to get convergence for certain parameters $k,a,b$. So this method, while more accurate, is not very efficient.

Can anyone think of a better approach which is both accurate and efficient?

===EDIT===

Based on the comments below, I propose the following test case:

```<code>k = 3
a = 1.017265193370166e+00
b = 1.018359131242105e+00
c0 = 1
c1 = -2.742386086956186e+03
c2 = 5.013787633287243e+06
c3 = -4.583247149493306e+09
</code>```

According to Mathematica's NIntegrate with 60 decimal digits of Working Precision, the answer is

```<code>truth = 0.00028808155205599365876989193328612177562110
</code>```

With my code for Method 1 (given below), I get the following results in double and quad precision:

```<code>double_result = 2.880813553929329e-04 relerr = 6.826645419623634e-07
quad_result   = 2.880815520711622e-04 relerr = 5.265359890594515e-11
</code>```

The quad precision result is acceptable, but the double precision result is a bit disappointing, and it seems I should be able to get more accuracy. Here is my double precision code, which just integrates term by term:

```<code>      result = c0 / (1.0 + k) * powdiff(b, a, 1 + k);
      result += c1 / (2.0 + k) * powdiff(b, a, 2 + k);
      result -= a * c1 / (1.0 + k) * powdiff(b, a, 1 + k);
      result += 0.5 * c2 / (3.0 + k) * powdiff(b, a, 3 + k);
      result -= a * c2 / (2.0 + k) * powdiff(b, a, 2 + k);
      result += 0.5 * a * a * c2 / (1.0 + k) * powdiff(b, a, 1 + k);
      result += c3 / 6.0 / (4.0 + k) * powdiff(b, a, 4 + k);
      result += 0.5 * a * a * c3 / (2.0 + k) * powdiff(b, a, 2 + k);
      result -= 0.5 * a * c3 / (3.0 + k) * powdiff(b, a, 3 + k);
      result -= a * a * a * c3 / 6.0 / (1.0 + k) * powdiff(b, a, 1 + k);
</code>```

The powdiff function computes $a^n - b^n$ is such a way as to minimize catastrophic cancellation:

```<code>static double
powdiff(const double a, const double b, const int n)
{
  double a_to_n = pow(a, (double) n);
  double result = -a_to_n * expm1(n * log1p((b - a) / a));
  return result;
}
</code>```

Any suggestions on how to improve the result are most welcome.

Title: Integration problem


问:
Accuracy of the Crank-Nicolson method for non-linear, inhomogeneous heat equation
I am currently coding a solution to the following PDE:

$\frac{\partial T }{\partial t} =\frac{\partial}{\partial \theta}(A(\theta ,\phi )\frac{\partial T }{\partial \theta}) +\frac{\partial }{\partial \phi}(B(\theta ,\phi )\frac{\partial T}{\partial \phi}) + F(\theta ,\phi)T$ with periodic boundary conditions in polar coordinates.

I have been able to get some code which solves the PDE via the Crank-Nicolson method and matrix inversion, the only problem that I have is with the accuracy of the method, according to Wikipedia, the Crank-Nicolson method is unconditionally stable but ''the approximate solutions can still contain (decaying) spurious oscillations if the ratio of time step $\Delta t$ times the thermal diffusivity to the square of space step, $\Delta x^2 $, is large (typically, larger than 1/2 per Von Neumann stability analysis).''.

How can I check the accuracy of my method in a similar way that Wikipedia has done? Obviously taking a small timestep would improve the method, but would also increase runtime, so I'd rather only increase it if my method's accuracy is worse than say 1% from the actual answer. I am also unsure whether the periodic boundary conditions aid or detriment to the accuracy of my solution as a whole.

I hope someone can point me in the right direction.

edit: I have reposted this from math.stackexchange as I was told there was a higher chance of getting an answer here. I will cross post any useful answers between the 2.

edit: I have gone with the Backward Euler method, thanks to @David for the help though.


问:
Grid search for bi-level optimization
Apologies if this isn't the best place to ask this question, and further apologies for such a basic question (I am a secondary school graduate and have not learned very much yet). Please direct me to a better place to ask this if this is not the best.

Suppose we wanted to solve the following optimization problem:
$$\inf_{x \geq 0}\sup_{y \in [0, 1];\ z > 0} f(x, y, z),$$
where $f(x, y, z)$ is some objective function with a closed form that can be specified in terms of parameters $x$, $y$, $z$. How would we implement this optimization, say using the Python programming language?

I am only vaguely familiar with implementing optimization in the Bayesian setting (eg using variational inference) or by a grid search. But I do not think Bayesian optimization works here as there are any priors here (I could be wrong though) and I'm not quite sure how to discretize these parameters in a reasonable way. Any guidance about these, or anything at all, will be quite helpful.

(The motivation for the optimization problem above comes from my personal study of robust optimization models from the paper https://arxiv.org/pdf/1908.05659.pdf.)

Title: Grid search for bi-level optimization


问:
How to solve evolution equation numerically?
How do I solve the following evolution equation numerically:

$$
\dot{\boldsymbol{\mathcal A}} = -\lambda A_1|\boldsymbol{s}-\boldsymbol{\alpha}|^2\left[(\boldsymbol{n}_r:\boldsymbol{\alpha})\boldsymbol{n}_r\otimes\boldsymbol{n}_r+\frac{3}{2}A_2\boldsymbol{\mathcal A}\right]
$$

Here $\boldsymbol{\mathcal A}$ is a fourth-order tensor, $\boldsymbol{s}$ is the deviatoric stress, $\boldsymbol{\alpha}$ is the back stress, and $\boldsymbol{n}_r$ is $\frac{\boldsymbol{s}-\boldsymbol{\alpha}}{|\boldsymbol{s}-\boldsymbol{\alpha}|}$. $A_1, A_2$ are material constants, and $\lambda$ is a parameter that is dependent on $\boldsymbol{s},\boldsymbol{\alpha},\boldsymbol{\mathcal A}$.

This is the evolution equation from the paper "Directional distortional hardening in metal plasticity within thermodynamics", where they got this equation from Clausius-Dirhem inequality.


问:
How to get an "optimal point" for refinement in FEM adaptive mesh refinement?
Consider the following 1D problem

\begin{align*}
\begin{cases}
\displaystyle
-\frac{d^2u}{dx^2} = f(x), \hspace{0.5cm}  x\in (a,b) \\[4mm]
u(a) = u_{a}, \ \ u(b) = u_{b}
\end{cases}
\end{align*}

I know how to compute a finite element approximation $u_{h}$ of this problem. I'm working in an h-adaptive method, it works fine but when I refine an element, I divide it in the element midpoint.

I want to refine the element at an "optimal point". Suppose I want to refine the element $I_{k} = [x_{k},x_{k+1}]$. In order to find the "optimal point" $x^{*}\in (x_{k},x_{k+1})$ for the refinement of $I_{k}$, I know that this could be done with the Dual Weighted Residual, that is, we consider the finite element residual

\begin{align*}
R = \frac{d^2u}{dx^2}+ f
\end{align*}

and the value $\eta_{k} = r_{k}w_{k}$, that relates the local residual $r_{k}$ and the weight $w_{k}$ (which comes from the dual problem).  But I don't understand how this is related to the optimal point $x^{*}$.

Do you know how to use these quantities in order to find the point $x^{*}$?

Do you know books or articles where I can study this problem?


问:
Is there an obvious n-dimensional elasticity complex?
Finite element exterior calculus helps to furnish nice stable elements for many common elliptic PDE by looking at the de Rham complex

$$\Lambda^0(\Omega) \stackrel{d_1}{\to} \Lambda^1(\Omega) \stackrel{d_2}{\to} \ldots \stackrel{d_n}{\to} \Lambda^n(\Omega)$$

where $\Lambda^k(\Omega)$ are the $k$-forms on the domain $\Omega$ and $d_k$ is the exterior derivative.
For example the mixed form of the Poisson equation and various other vector Laplace equations are nicely described by other parts of the sequence.
Moreover, the de Rham complex is relatively easy to describe in arbitrary dimensions, not just $n = 2$ or 3.

There are also chain complexes that describe the elasticity equations in dimensions 2 or 3.
For example, in dimension 2, we can use the sequence

$$\phi \mapsto \partial_k\partial_k\phi\cdot\delta_{ij} - \partial_i\partial_j\phi,$$
$$\sigma_{ij} \mapsto \partial_j\sigma_{ij}.$$

The composition of these two maps is zero so it's a perfectly nice complex.
There's a similar generalization to 3D.
What, if any, is the natural generalization of the elasticity complex to arbitrary dimensions?
I've only ever seen it described specifically for dimensions 2 or 3.


问:
Under what specific circumstances can a matrix lose symmetry?
This is a general question: In this blog post, the author points out that under a unitarily invariant norm, the closest symmetric matrix to the square matrix $M$ is $\frac 1 2 (M + M^T)$. The general question I have is: When might a matrix which ought to be symmetric lose symmetry? I understand generally that this may result from rounding errors or sensor errors. The question is: In what concrete situations might this happen?

The author suggests some examples but I'd like to see more: For example, a finite difference approximation to a Hessian matrix can be nonsymmetric even though the Hessian is symmetric. In some cases, lack of symmetry is caused by rounding errors.

These examples are rather vague though.


问:
Scheme Solving CCZ4 Formalism Numerically
For my thesis, I have to implement the CCZ4 formalism into an existing code (I tried to type the evolution equations here but does not seem to work).

However, I can't seem to come up with a scheme to solve these equations. For instance, here, the authors give an explicit scheme on how to solve the ADM equations (see page eight). I am looking for something similar. I tried something similar but that did not get me far. Furthermore, I am aware of the codes such as  GRChombo (see here) that are able to solve these equations but I can't find anywhere in what order they solve these equations. I tried looking at the source code of GRChombo but I could not anything useful that would solve my problem (I am not familiar with C++)

Hence, I am looking for references or ideas on how I could solve my problem.

Thanks in advance!


问:
Eigenvalues of same operator expressed in two different orthonormal basis are coming out different
I have an operator $H$. I express $H$ as a matrix in the orthonormalized $\{ |e > \}$ basis. Then I diagonalize it to obtain eigenvalues, let's say for example $H$ is $6 \times 6$ and the eigenvalues are $1, 5, 7, 9, 3, 8$. Now I express $H$ as a matrix in a different orthonormalized $\{ |f > \}$ basis obtained through unitary transformation of the $\{ |e > \}$ basis. Now I diagonalize this matrix. In principle, I should get same result $1, 5, 7, 9, 3, 8$. But when I run the code, only four eigenvalues are matching - let's say I get $1, 5, 7, 9, -4, 15$. Notice that the last two eigenvalues don't match, but peculiarly their sum is same as before, i.e., $3 + 8 = -4 + 15 = 11$. I am using ${numpy.linalg.eig}$ . Please enlighten me.


问:
Advice for a beginner to a computational physics problem
I am a beginner in computational physics but have a background in physics, math, and programming (from online courses, videos, and textbooks) and I would like to know how should I:

Approach and think of any general computational physics problem and how I might find a solution and apply it to a real-world physical experiment that implements the model with a physical device to test

I do not have any specific problem in mind right now (as I invalidated my previous physics issue and will not continue working on it, as seen in the comments).

Which textbooks and/or online courses or resources I could learn from that will help me with general computational physics problems and solutions. I would like to learn universal principles and rules of thumb and any seasoned tips from anyone here.

I am looking to use COMSOL for the simulations in case anyone was wondering which software I will use.


问:
Is there any book about fundamental FEM theory using similar terminology as Comsol MultiPhysics?
I am considering  to solve complexed PDE systems, like in this post, using Comsol MultiPhysics. The PDEs are different from the General Form provided by Comsol. The Weak Form module may be worth a shot.

But I am no export in FEM theory. Even if I become one after learning, there might be a distance from FEM theory to Comsol weak form implementation, say different terminology.

So, I was wondering if there are some books or materials that details the fundamentals of FEM empolying terminology similiar to Comsol, maybe like this website given by Comsol but with much more details as a textbook on FEM theory.
Any suggestions would be appreciated.


问:
How is the transfer function of a state-space representation numerically computed?
This question is a duplicate of this question that I asked on dsp stack exchange. However, nobody had the answer there but the question seems more appropriate on this forum (If not feel free to tell me, it's my first post here).

The question is the following :

I've got my state space representation
$$ \dot{X} = AX+BU $$  $$Y = CX$$

with Laplace we get :
$$Y(s) = C(A-sI)^{-1}BU(s) = H(s)U(s)$$  And I want to evaluate the coefficients of the rational function matrix : $$H(s) = C(A-sI)^{-1}B$$

my first idea was to use the cofactor formula :
$$(sI-A)^{-1} = \frac{1}{det(sI-A)}Com(sI-A)^T$$ with $Com(sI-A)^T$ the matrix of the cofactors of $sI-A$.

my second idea was
$$A-sI = -s(I-\frac{1}{s}A)$$

and then if you have $(\frac{1}{s}A)^{n}  \xrightarrow[n \to \infty]{} 0$ for $|s|$ above a certain threshold:
$$(I-\frac{1}{s}A)^{-1} = \sum_{k=0}^{\infty}(\frac{1}{s}A)^{k} \approx \sum_{k=0}^{N}(\frac{1}{s}A)^{k}$$

Although I feel like it's a bad idea somehow.

I also checked the ss2tf scipy function code but I really don't get what they do (especially in the loop l 278/280).

Can anyone enlighten me?


问:
Strategies to solve an equation with a polynomial and a numeric function
I have to solve numerically an equation of the following form:

$$
\sum_{n=0}^m c_n x^n = f(x) x^k
$$

Where the $c_n$ are real values, $k$ is an integer and $f$ can only be evaluated numerically.

The current implementation in our code is using a numerical approach (specifically using scipy.optimize.brentq) but I wondered:

If I could find all the roots somehow

If there was known strategies that exploit the polynomial parts.

EDIT:
I realized I have a few properties on $f$:

$f$ is in $C^\infty(\Bbb{R}, \Bbb{R})$

$f$ is monotonous (non increasing)

$\lim_{x\rightarrow+\infty} f(x) = -\infty$

$\lim_{x\rightarrow-\infty} f(x) = +\infty$


问:
Are there better methods than Levenberg-Marquardt for large-scale nonlinear least squares problems?
I am using Levenberg-Marquardt to solve a nonlinear least squares problem with 16 million residuals and about a million variables. I have a fast matrix vector product for both the Jacobian and its transpose, and I can also compute the norms of the columns of the Jacobian. To compute the step, I am using LSQR. This scheme is probably best when the Jacobian is formed explicitly, and forming the Jacobian is quite expensive. So once we have it, it pays to solve for the step quite accurately. In my case, it seems wasteful to evaluate the matrix vector product many times at the same point, as opposed to making a small step every time a matrix vector product is evaluated. This would be the case in the conjugate-gradient method, for example, but its convergence is much slower. I’ve tried limiting LSQR to just a few iterations but that doesn’t seem to help.

Is there a better method for when we have only a fast matrix vector product?


问:
Pyton get fem solving simple 2D differential equation?
I need to solve the differential equation:

$$u_t = \Delta u + \sin t$$

On a 2D domain with homogeneous Dirichlet boundary conditions. Tod o this I am trying to use the python package getfem.

I am struggling quite a bit trying to set things up especially the parts of $u_t$ and $\sin t$ as I am not sure how to define them based on the docs.

If I understand properly this:

import getfem as gf
import numpy as np
import math
import pyvista as pv
from pyvirtualdisplay import Display

center = [0.0, 0.0]
dir = [0.0, 1.0]
radius = 1.0
angle = 0.2 * math.pi

mo = gf.MesherObject("cone", center, dir, radius, angle)

h = 0.1
K = 2
mesh = gf.Mesh("generate", mo, h, K)

# Select boundary.
outer_faces = mesh.outer_faces()
OUTER_BOUND = 1
mesh.set_region(OUTER_BOUND, outer_faces)

sl = gf.Slice(("none",), mesh, 1) #Turn mesh into triangles

mfu = gf.MeshFem(mesh, 2)
elements_degree = 2
mfu.set_classical_fem(elements_degree)
```
<p>Will give me a cone mesh for the domain (a triangle). Which is what I want. And the line:</p>
<p><code>mfu = gf.MeshFem(mesh, 2)</code></p>
<p>Is creating an mesh with a 2D vector field on top of it (I think?).</p>
<p>Then I can do this:</p>
```md = gf.Model("real")
md.add_fem_variable("u", mfu)
```

To introduce the variable I am solving for.

But from here on out I don't know how to set up the equation and tell getfem to solve it.


问:
Approximating eigenvalues of DPR1 matrix with special properties
In my application, I have a sum of diagonal $A$ and rank-1 $B$

$$T=\underbrace{\text{diag}(1-2\alpha h+2\alpha^2 h^2)}_A + \underbrace{\alpha^2 hh^T}_B$$

Where $h$ is a vector $\in \mathbb{R}^d$ with positive entries decaying faster than $1/i$, $\alpha>0$ and $\|T\|<1$.
Can eigenvalues of $T$ be approximated faster than $O(d^2)$?

Earlier answer pointed to $O(d^2)$ algorithms for finding full eigendecomposition for general DPR1 matrix, but I'm wondering if the specific structure is useful here.

Motivation: this matrix describes evolution of errors of SGD with step-size $\alpha$ when used to solve Ordinary-Least Squares problem with Hessian $\text{diag}(h)$ and normally distributed observations. In contrast to standard gradient descent where error evolves as $\text{diag}(1-2\alpha h+\alpha^2 h^2)$


问:
boundary condition at rotation axis of a spherically-symmetric system
The quantity I am interested in is not the rotation rate $\Omega$, but I will use this quantity nonetheless to make the problem clearer. I am interested in a spherically-symmetric system and in the boundary condition of the rotation rate $\Omega=u/r\sin\theta$. Now, at $\theta=0$, we have an undetermined form $0/0$ and we must use de L'Hopital rule:
$$
\Omega (r,0)=\text{lim}_{\theta\rightarrow 0}\frac{du/d\theta}{r\cos\theta}\neq 0
$$

which is not zero as expected. Using central differences (I assume this is still valid):

$$
\frac{du}{d\theta}|_{\theta=0}=\frac{u_1-u-1}{2\Delta\theta},
$$

where the index for radius was ignored (the minus 1 should be as an index but it somehow doesn't work within the fraction). Now, $u_{-1}=-u_1$ and:

$$
\frac{du}{d\theta}|_{\theta=0}=\frac{u_1}{\Delta\theta}=\frac{\Omega_1 r\sin\theta_1}{\Delta\theta}\simeq\Omega_1 r\quad\quad\text{[for $\theta_1=\Delta\theta$ small]}.
$$

Finally:
$$
\Omega_0=\Omega_1.
$$

Now, if I decide to discretize the equations regularly in $\mu=\cos\theta$, so that:

$$
\sin\theta=\sqrt{1-\mu^2}\quad\quad\frac{1}{d\theta}=-\frac{\sin\theta}{d\cos\theta}=-\frac{\sqrt{1-\mu^2}}{d\mu},
$$

then:

$$
\Omega(r,1)=\text{lim}_{\mu\rightarrow 1}\frac{du/d\mu}{rd\sqrt{1-\mu^2}/d\mu},
$$

and

$$
\Omega(r,1)=\text{lim}_{\mu\rightarrow 1}-\frac{\sqrt{1-\mu^2}}{r\mu}\frac{du}{d\mu}=0
$$

What am I missing here? This should definitely not be zero at the pole (or anywhere).


问:
Interface condition for 1D Helmholtz equation using finite element method
I want to implement a 1D Helmholtz equation with jump condition. The domain is $x=[0,1]$ and both ends have Dirichlet boundaries($p$=0). The 1D strong formulation is;

$$c^2\nabla^2p + w^2p=0 \qquad \mathrm{in} \quad \Omega,$$

and

$$\nabla p(0)=\nabla p(1)=0 \qquad \mathrm{at} \quad \partial \Omega,$$
$$\nabla p^-(0.8)=K[[p]]=K(p^+(0.8)-p^-(0.8)) \qquad \mathrm{at} \quad x=0.8, $$

where $c$ is the speed of sound, $\omega$ is the eigenvalue, the $p$ is the eigenfunction. $^-$ and $^+$ superscripts denote upstream and downstream directions and $K$ denotes a complex coefficient.

The weak(finite element) formulation without jump condition term is;

$$-\int_\Omega c^2\nabla u \nabla v dx + \int_{\partial\Omega} c^2\nabla u v ds + w^2\int_\Omega u v dx = 0$$

where $u$ and $v$ are test and trial functions. I just wonder how can I include the jump condition in the weak formulation? And how it will affect the implementation of stiffness and mass matrices?


问:
2D FEM element with 6 DOF per node (including rotation)
I'm writing finite element code in 3D, and I wish to connect vertical (Euler-Bernoulli) beams to a horizontal plate. The deformations of a beam are well described, but for the plate it is more difficult: for in-plane deformations of the plate I use quadrilateral elements, and for out-of-plane deformations I use Kirchhoff-Love plate theory. However, there is one degree of freedom that is missing in the plate: the one that matches the torsional degree of freedom of the beam. This would require in-plane rotation of the plate. Is there any theory about in plane rotation of a plate?

So, I would like to implement 2D elements in a 3D environment, which includes rotation in all degrees of freedom (so 6 DOFs per node). Which element type is suitable for that purpose?


问:
A way to solve nonsmooth stiff ODEs
Let us considered the following ODEs
\begin{align*}
    \dfrac{dX}{dt} = F(X), \tag{1.1}
\end{align*}
where the unknown $X\in \mathcal{D} \subset \mathbb{R}^l$ and it can be stiff. These ODEs can be solved numerically by classical Newton's method if $F$ is differentiable almost everywhere on the considered domain
\begin{align*}
    \dfrac{X^{n+1}-X^n}{\Delta t} = F\left(X^{n+1}\right).
\end{align*}
In cases of $F$ is not differentiable everywhere, I have read in a lecture that there is a way to solve numerically by replacing the system (1.1) with the following system
\begin{align*}
    X^{n+1} - X^n - (1-\nu)\Delta t F(X^{n+1}) & = 0 \in \mathbb{R}^l\tag{1.2a}\\
    \det\left[I - (1-\nu)\Delta t \nabla F(X^{n+1}) \right] - w &= 0 \in \mathbb{R} \tag{1.2b} \\
    \min(\nu,w) &= 0 \in \mathbb{R} \tag{1.2c}
\end{align*}
where $\nu$ and $w$ are two new unknowns to be sought for simultaneously to $X^{n+1}$. I wonder how they can do like this and what is technique they use?

Title: A way to solve nonsmooth stiff ODEs


问:
It is possible to solve integro-differential equations using in Fenics?
I am interested in solve the following integro-differential equation:

\begin{align}
    \frac{\partial{\rho(\theta, t)}}{\partial{t}} = D \frac{\partial{\rho(\theta, t)}}{\partial{\theta^2}} - \beta \frac{\partial{}}{\partial{\theta}}\left(\rho(\theta, t) \left[ \int_{\theta}^{\theta + \pi} \rho(\theta', t) d\theta' - \int_{\theta - \pi}^{\theta} \rho(\theta', t), d\theta' \right] \right)
\end{align}

As far as I know, it is possible to solve the diffusion-convection equation using FEM in Fenics, but I am not sure if it is possible to implement the integrals in Fenics. Do you think it is possible? Any idea or advice would be very welcome.

I will be happy to share more details if necessary.


问:
How to calculate the reynolds number for airborne wind energy system?
I am asking for some help on how to calculate the reynolds number for an AWES. I know the formula of Re=rhovl/viscosity. I have selected the air foil as Clark y with chord length 3.72m (which will be the characteristic length). I specifically need the values of rho, v, viscosity with changing altitude. And also if it is based on geographical location India it would be helpful.


问:
A staggered grid for an eigenvalue problem (linear stability analysis)
I'm interested in extending the concept of a staggered grid (commonly used to solve the incompressible Navier-Stokes equations) to a linear stability analysis context. For example, we can consider flow within a two-dimensional rectangular duct, as in Theofilis, Duck, and Owen in the Journal of Fluid Mechanics, 2004. The base flow here is $\bar{u}_z(x,y)$ only, just to simplify our discussion. The staggered grid (below) can be called a marker-and-cell (MAC) or a particle-in-cell (PIC) mesh.

The inviscid temporal linear stability governing equations (solved as an eigenvalue problem) would be:

$-i\beta\hat{u}-\partial\hat{p}/\partial x=-i\omega\hat{u}$

$-i\beta\hat{v}-\partial\hat{p}/\partial y=-i\omega\hat{v}$

$-\frac{\partial \bar{u}_z}{\partial x}\hat{u} - \frac{\partial\bar{u}_z}{\partial y}\hat{v}-i\beta\hat{w}-i\beta\hat{p}=-i\omega\hat{w}$

$\frac{\partial\hat{u}}{\partial x }+ \frac{\partial\hat{v}}{\partial y} - i\beta\hat{w}=0$

Wouldn't combining differentiation with interpolation into the same matrix cause errors? Particularly around the boundaries!

A scheme could be:

Differentiation performed by finite differences

Discretize the continuity and z-momentum equations at the cell centers

Discretize the x-momentum equation at the x-edges

Discretize the y-momentum equation at the y-edges

Interpolate the $u_x$ values onto cells in the x direction and interpolate the $u_y$ values onto cells in the y direction for the z-momentum equation. Then multiply by the matrix for $\partial\bar{u}_z/\partial x$ or $\partial\bar{u}_z/\partial y$, respectively.

All other terms that need to be transferred from one grid to the other involve derivation and interpolation. For example, $dp/dx$ in the x-momentum equation could be treated as: create an interpolation matrix I from the pressure cells to the x-edges and then apply a differentiation matrix D to the result i.e., DI.

Steps 4 and 5 seem likely to introduce errors. Also, this scheme seems to be the only way to obtain a square eigenproblem.

Am I missing something?

Also, would the scheme need to be changed if a spatial linear stability analysis were performed instead?

I have seen Interpolation of velocities on staggered grid (in PIC), which explores interpolation, and is where I obtained the photo of the staggered grid.


问:
meshing of parametric 3d models
I would like to do CFD analysis on a parametric model and would like to obtain mesh deformation in the parametric space.

As a concrete 2D example let us consider a house shape with a bottom rectangle and a triangle sitting on top. This shape has 3 parameters: (A) the width of the rectangle, (B) the height of the rectangle and (C) the height of the triangle. Let us choose A=10, B=8 and C=4, and then I can calculate the coordinates of this 5-point polygon from this. Moreover, I can also calculate the partial derivatives of these coordinates with respect to A, B and C.

So we have a boundary polygon where the partial derivatives of the coordinates with respect to a list of parameters are known. I would like to mesh this polygon in a way that (1) it is "nice" for CFD analysis, which many tools can achieve, (2) I would like to know the partial derivatives of the added mesh points with respect to my original parameter space.

I have not found any way to do (2). Basically I have arbitrary number of parameters assigned to initial points and would like to do meshing so that these parameters are linearly propagated (when adding a midpoint of a line segment then the parameters should be averaged). Is there any way to do this?

Title: meshing of parametric 3d models


问:
How to include zero flux boundary conditions?
I am trying to solve the following differential equation in the domain of $\theta \in [0, 2 \pi]$ using finite differences scheme:

For $0< \theta \leq \pi$
\begin{align}
\rho_i^{n+1}=\rho_i^{n}+D\frac{\Delta t}{\Delta\theta^2}[\rho^n_{i+1}-2\rho^n_i+\rho^n_{i-1}]-\beta \frac{\Delta t}{\Delta\theta}[\rho^n_i-\rho^n_{i-1}]I_{1}(\theta, t)
+\beta\Delta t\rho_{i}^n[2\rho_{i}^n(\theta+\pi)-2\rho_{i}^n(\theta)]
\end{align}

\begin{align}
I_{1}(\theta, t) = \int_{\theta}^{\theta+\pi} \rho(\theta', t)d\theta' - \int_{0}^{\theta} \rho(\theta', t) - \int_{\theta+\pi}^{2\pi} \rho(\theta', t) d \theta'
\end{align}

For $\pi<\theta < 2\pi$:
\begin{align}
\rho_i^{n+1}=\rho_i^{n}+D\frac{\Delta t}{\Delta \theta^2}[\rho^n_{i+1}-2\rho^n_i+\rho^n_{i-1}]-\beta\frac{\Delta t}{\Delta t\theta}[\rho^n_i-\rho^n_{i-1}]I_{2}(\theta, t)
+\beta\Delta t\rho_{i}^n[2\rho_{i}^n(\theta-\pi)-2\rho_{i}^n(\theta)]
\end{align}

\begin{align}
I_{2}(\theta, t) = \int_{0}^{\theta-\pi}\rho(\theta', t)d\theta' + \int_{0}^{2\pi} \rho(\theta', t) d\theta' - \int_{\theta - \pi}^{\theta} \rho(\theta') d\theta'
\end{align}

I have added zero flux boundary conditions using:
\begin{align}
    \frac{\partial{\rho}}{\partial\theta} &\approx \frac{\rho_{i}^{n} - \rho_{i-1}^{n}}{\Delta \theta} = 0
\end{align}

Then, $\rho_{i}^{n} = \rho_{i-1}^{n}$ should assure zero flux boundary conditions.

Here is my code:

```<code>
#1. Import libraries
import numpy as np
import matplotlib.pyplot as plt
import scipy as sci

from scipy import integrate
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

#2. Setup
ntheta = 41
dtheta = 2 / (ntheta - 1)
nt = 500     # the number of timesteps
D = 1.     # Diffusion parameter
beta = 5. #0.  # Advection parameter CHANGE 0 to 5 and viceversa
sigma = .2 # For Counrat number
dt = sigma * dtheta**2 / D # dt is defined using sigma

rho = np.zeros(ntheta)                   # Array with ntheta elements all equal to 1.

rhon = np.zeros(ntheta) # rhon, to advance the solution in time
theta = np.linspace(0, 2*np.pi, ntheta)

for i in range(len(theta)):
  if theta[i] <= (np.pi - 0.25):
    low = i

for j in range(len(theta)):
  if theta[j] <= (np.pi + 0.25):
    up = j

rho[low:up] = 2 #the area under the curve should be 1.
angle = []

for i in range(len(theta)):
    if theta[i]==np.pi:
        angle=i

#3. Solve the EDO (iterate finite difference scheme)
for n in range(nt):  #iterate through time
    rhon = rho.copy() ## copy of the existing values of rho into rhon
    for i in range(1,ntheta-1): # iterate through angle (space)
        if theta[i] <= np.pi:
          #Using second-order central, first-order backward difference techniques, generalization of the Second Fundamental Theorem of Calculus.
            rho[i] = rhon[i] + D * dt / dtheta**2 * (rhon[i+1] - 2 * rhon[i] + rhon[i-1]) - beta * dt / dtheta * (rhon[i] - rhon[i-1]) * ((sci.integrate.simpson(rhon[i:i+angle], theta[i:i+angle])) - sci.integrate.simpson(rhon[:i], theta[:i]) - sci.integrate.simpson(rhon[i+angle:], theta[i+angle:])) - 2 * beta * dt * rhon[i] * (rhon[i+angle] - rhon[i])
        else:
            rho[i] = rhon[i] + D * dt / dtheta**2 * (rhon[i+1] - 2 * rhon[i] + rhon[i-1]) - beta * dt / dtheta * (rhon[i] - rhon[i-1]) * ((sci.integrate.simpson(rhon[:i-angle], theta[:i-angle])) + sci.integrate.simpson(rhon[i:], theta[i:]) - sci.integrate.simpson(rhon[i-angle:i], theta[i-angle:i])) - 2 * beta * dt * rhon[i] * (rhon[i] - rhon[i-angle])

    rho[0] = rho[1]
    rho[len(rhon)-1] = rho[len(rhon)-2]

plt.plot(theta, rho)
plt.ylim (0, 2)
plt.axhline(y=0.16, linestyle='--', color='k') # Area = b*h i.e 1 = 2*pi*h => h = 1/2pi
plt.xlabel("Theta")
plt.ylabel("Rho")
plt.show()
</code>```

The initial area under the curve is equal to 1, then, in zero flux this should be constant $Area = 1 = b h = 2 \pi h$, $h \approx 0.16$. However, this seems to works only when $\beta=0$ (i.e the area is always 1.), otherwise, the area under the curve is less than 1 (e.g. when $\beta = 5.$)

Please, could you help me to find my mistakes or misunderstandings?


问:
Looking for an algorithm to find points where lines intersect most frequently
I have a file that contains the coordinates of two points that make up a line. When I plot these lines, I notice that there are areas where the lines intersect more often. These points are identified as the areas where the yellow patterns are the thinnest on the plot. To make it easier to locate these points accurately, I plot them as a 2D histogram.

However, even with the 2D histogram, it's still a bit challenging to identify these points precisely. Therefore, I am wondering if anyone knows of an algorithm that can be used to find these points where the lines intersect most frequently. Any suggestions or ideas would be greatly appreciated.

Thank you!


问:
Which dense matrices are hard to invert?
Suppose I'm solving $Ax=b$ for dense $m\times d$ matrix $A$. For which $A$ is this hard to do?

More concretely, is there any work on estimating the error after $k$ steps of iterative solver, $k\le d$, for some specific solvers, based on properties of $A$? For random isotropic vector $x$ and large $d$, a useful property proved here is that

$$T^k x \approx c \operatorname{Tr}(T^k)$$

Hence this comes down to understanding how fast $\operatorname{Tr}(T^k)$ decays where $T$ is the iteration matrix of corresponding solver. Which properties of $A$ result in $\operatorname{Tr}(T^d)$ being large?

Classical approach formulates difficulty in terms of condition number $\kappa(A)$, but this comes into play after many steps, ie $k>d^2$. When $d$ is huge, we can only afford $k=O(d)$ steps. For instance, modern language model might have $d\approx$10B and be trained for 1B steps. A linear analogue of this problem is to solve $Ax=b$ where $A$ is a dense $10B\times200B$ matrix, in less time than it takes to read the matrix $A$.

For "single step gradient descent", you can show that an expression involving second/third moments of the matrix determine the difficulty of approximate inversion for gradient descent using $k=1$, but I feel like I'm reinventing the wheel here

$$\frac{\text{loss}_1}{\text{loss}_0}\approx1-\frac{\operatorname{Tr}(H^2)^2}{\operatorname{Tr}(H^3)\operatorname{Tr}H}$$

where $H=A^TA$


问:
Analytical Equation of the gaussian 1D wave equation with periodic Boundary condition
I am trying to validate the 1D analytical wave equation with a numerical solution with periodic boundary conditions. I have implemented the periodic boundary condition for the numerically calculated wave but have not been able to implement it for the analytical solution. The analytical solution for the the 1D gaussian wave is
 U(x,t) = e^((x-x0-c*t)**2)* cos(k0*(x-x0-c*t))
where x0 is the starting point of the wave (the midpoint of the wave would be at this location), c being the speed, and k0 is the wave number. The value of x varies (i.e, the domain) varies from -L to L.

I want if the moving wave crosses the domain L the wave packet again starts with -L (staring of the domain).
Any suggestions would be of great help.

What changes are required to be made to the equation to get an analytical solution with periodic boundary conditions implemented?


问:
Convolution/weighted average of two arrays in Python
I have an equation that I need to calculate numerically, but I am having doubts about my approach. I am cross-posting this question from Stack Exchange, because I am not getting any responses.

This is the equation I need to program:

where $\chi$ is a volume mixing ratio in the stratosphere, $\chi_0$ is a volume mixing ratio in the troposphere and G is an age spectrum (a distribution of times it takes an air parcel to travel from troposphere to stratosphere).  G is an inverse Gaussian:

where t' is defined as the transit time $\Delta$t divided by $\Gamma$.
The tropospheric volume mixing ratios as a function of time form an almost linear curve. I'll put the curves below (ignore the dashed lines). The time lags are with respect to 2009, i.e. a time lag of 1 corresponds to 2008 etc. Now I need to apply the first equation to these two curves numerically. I have tried numpy.convolve(chi_0[::-1], G, mode='same') and variations thereof, but I'm having doubts about the results, so any feedback is greatly appreciated. Do I have to multiply by the time step? Even then, however, the result seems wrong to me (too small). I have to "weigh" the tropospheric curve with the age spectrum.
So to state it clearly, my question is: How do you convolve the tropospheric curve with the spectrum G, as prescribed by the first equation?

For those interested, the method is described in Stiller et al. 2012 under "non-linearity correction".

The result I get when I do a convolution of the age spectrum and the tropospheric reference curve looks something like this:

This is just for one Γ. You can do it for any Γ you like. According to the method, the modeled SF6 (i.e. the convolution) at measurement time (2009) should be larger than the original SF6 (red dashed line), but in my case it is always smaller, unless I drop the *dt in the code, then it is always way too large. If you calculated the new SF6 for several different values of Γ, you apply the lookup method and pick the Γ that causes a value of SF6 that is closest to your measurement, to get the actual mean age (and not just the lag time). So if I understand the method correctly, the blue curve should come above the red dashed line to get an SF6' > SF6.

Here is the code:

```<code>import numpy as np
import xarray as xr
import gdf_tools as GDF
import matplotlib.pyplot as plt
import pandas as pd
import sys
from datetime import datetime as dt
import time
from scipy import optimize
import scipy as sp
from scipy import interpolate
import datetime
from scipy.signal import savgol_filter
from scipy import stats
from scipy import signal

#---------------------------------------------------------------------------
#                             FUNCTIONS
#---------------------------------------------------------------------------

def smooth(y):
    return signal.savgol_filter(y, 21, 2)

def func(x, a, b):
    return a*x + b

def spec(t, *args):
    gamma, d = args
    tp = t / gamma
    return ( 1 / ( 2 * d * np.sqrt( np.pi * tp**3 ) ) ) * \
        np.exp( ( -(gamma**2) * ((tp - 1)**2 )) /( 4 * (d**2) * tp ) )

def getYfromX(y, xval, x):
    return np.interp(xval, x, y)

def getXfromY(y, yval, x):
    return np.interp(yval, y, x)

#---------------------------------------------------------------------------
#                               MAIN
#---------------------------------------------------------------------------
# Read tropospheric reference curve
ref=pd.read_csv('~/pyscripts1/aoa/refcurve/refcurve.csv')
ref.index=ref["time"].astype("datetime64[D]")

t = ref['decimal_time'].values # Decimal time
t_trunc = t[t<2009] # Truncate at measurement time
z = ref["GML_Global_SF6"].values
ref_trunc =  z[:len(t_trunc)] # Truncate at measurement time

# Linear approximation of the reference curve
a,b = np.polyfit(t_trunc, ref_trunc, 1)
fit = func(t_trunc, a, b)

# Smoothed reference curve
y=smooth(ref_trunc)

# Age of air (lag time) sequence to loop over
ages = np.linspace(3.5, 8.5, 200)

# Modeled (i.e. "propagatec") SF6 vmr
concentrations = []

# The measured vmr in 2009 is approx 5.25e-12
# Calculate a first guess age of air by projecting this measurement on
# the tropospheric reference and calculating the lag time
proj=getXfromY(y, 5.25, t_trunc)
print("First guess: ", 2009-proj)

# Loop over lag times/ages
for i in ages:
    gamma = i   # Mean age
    delta = np.sqrt(0.7*gamma)  # Spread of the age spectrum
    deltat = 2009 - t_trunc   # Lag times
    dt = t[1] - t[0]    # Time increment

    # Calculate age spectrum
    kernel = spec(deltat, gamma, delta)
    kernel[np.isnan(kernel)] = 0

    # Normalise age spectrum
    I = np.trapz(kernel, t_trunc)
    kernel = kernel / I

    # Calculate "weighted average" of the tropo. ref.
    model=np.convolve(ref_trunc, kernel, mode="same")*dt

    concentrations.append(model[-1])

# Plot kernel
plt.figure()
plt.plot(t_trunc, kernel)

# Plot tropospheric reference curve and the model
plt.figure()
plt.plot(t_trunc, model, label="convolution")
#plt.plot(t_trunc, ref_trunc, label="trop")
#plt.plot(t_trunc, fit, label="linear fit")
plt.plot(t_trunc, y, label="smooth trop")
plt.axhline(y=5.25, ls='--', c='red')
plt.axvline(x=proj, ls='--', c='black')
plt.legend()
plt.grid()
plt.xlabel("Time")
plt.ylabel("SF6 [vmr]")
plt.tight_layout()

# Plot modeled SF6 as a function of age
# This is the "look-up" method described in Fritsch et al. 2020
plt.figure()
plt.plot(ages, concentrations)
plt.axhline(y=5.25, ls='--', c='red')
plt.xlabel("Mean age of air")
plt.ylabel("SF6 [vmr]")
plt.grid()
plt.tight_layout()

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


问:
Are my boundary conditions in my stiffness matrix correct?
I am trying to find a P1 Lagrange finite element solution to the following ODE:

$\begin{cases}-u''+u'+u=f~~~~~~~~\text{in}~~(0,1)\\ u(0)=1, u'(1)=0\end{cases}$

Where $f(x)=-2e^{x}+2\left(1-x\right)e^{x}+\left(1-x\right)^{2}e^{x}$.

My procedure:

First I obtained the variational formula $\sum_{i=1}^{n+1}U_ia(\phi_i,\phi_j)=<f,\phi_j>$ where each $\phi_i$ is a hat function and $a(u,v)=\int_0^1u'v'+u'v+uv$, then solve our ODE is equivalent to solve $AU=b$ where $A=[a(\phi_i,\phi_j)]$, $U=[U_i]$ and $b=[<f,\phi_i>]$.

My next step is take $A(1,:)=[1,0,0,...,0]$, $A(n,:)=[0,0,...,1,1]$ and $b(1)= 1, b(n)=0$ to ensure that the solution meets the initial conditions (Is it correct?).

If I had my initial conditions correct I would not understand what my problem is since I checked in my code my matrix $A$ and vector $b$ and they are correct and yet I get an incorrect solution.


问:
Scipy.root not converging even when provided with initial guesses very close to solution
I've made a previous question here and also in SO wondering why only the fsolve solver converges for the simple one dimensional unsteady conduction problem

$$ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2T}{\partial x^2}$$

$\alpha = 0.001$

$T(x=0,t) = 60°C$

$T(x=L,t) = 25°C$

$T(x,t=0) = 25°C$

One of the answers alluded to initial values, from which i started to test and, to make things easier, i changed the value of alpha to

$\alpha = 0.000001$

Which means that the result profile shall be very close to simply

$T(x=0,t) = 60°C$

$T(x=L,t) = 25°C$

$T(x,t) = 25°C$

Which is what i used as initial condition for this new code

```<code>import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from scipy.optimize import fsolve
from scipy import optimize
import pandas as pd

t_max = 25
x_max = 1
N_points_1 = 100
N_points_2 = 11
dt = t_max/(N_points_1-1)
dx = x_max/(N_points_2-1)
N_variables = N_points_1 * N_points_2
alpha = 0.000001
T_x_0 = 60
T_x_L = 25
T_t_0 = 25

def EDP (v):
    T_N =  np.ones((N_points_1, N_points_2))

    #Variables
    k=0
    for i in range(N_points_1):
        for j in range(N_points_2):
            T_N[i,j]=v[k]
            k=k+1

    #Model
    sys = []
    ##Boundary Condition x = 0
    for i in range(N_points_1):
        BC1 = T_N[i,0] - T_x_0
        sys.append(BC1)

    ##Boundary Condition x = L
    for i in range(N_points_1):
        BC2 = T_N[i,-1] - T_x_L
        sys.append(BC2)

    ##Initial Condition
    for j in range(1,N_points_2-1):
        IC1 = T_N[0,j] - T_t_0
        sys.append(IC1)

    ## Energy Balance
    for i in range(N_points_1-1):
        for j in range (1,N_points_2 -1):
            EB = alpha*(T_N[i,j+1]-2*T_N[i,j]+T_N[i,j-1])/dx**2 - (T_N[i+1,j]-T_N[i,j])/dt  #partial(T_N,t_1,1)[i,j] - alpha*partial_2(T_N,x,2)[i,j]
            sys.append(EB)

    sys=np.array(sys)

    return sys

initial_guess = 25*np.ones((N_points_1,N_points_2))
initial_guess[:,0] = T_x_0
initial_guess[:,-1] = T_x_L
initial_guess = initial_guess.flatten()
print(np.sum(EDP(initial_guess)))
0.3464999999999999
</code>```

As you can see the value of the sum of the values of each equation is already close to 0

fsolve does converge, but even so, it takes a long number of iterations. Others either do not converge or straigth up diverge

```<code>Solution_T = fsolve(EDP,initial_guess,maxfev = 1200)

Solution_T2 = optimize.anderson(EDP,initial_guess, verbose = True,maxiter = 1000)
104:  |F(x)| = 2.36365e+93; step 1
OverflowError: (34, 'Result too large')

Solution_T3 = optimize.newton_krylov(EDP,initial_guess,verbose=True,maxiter = 1000)
999:  |F(x)| = 0.00372199; step 1
NoConvergence: [60.         24.99999991 24.99999972 ... 24.99990242 24.99977307
 24.99968681]

Solution_T4 = optimize.broyden1(EDP,initial_guess,verbose=True,maxiter = 1000)
659:  |F(x)| = 4.37351e+151; step 1
OverflowError: (34, 'Result too large')

Solution_T5 = optimize.broyden2(EDP,initial_guess,verbose=True,maxiter = 1000)
999:  |F(x)| = 2.81888e+10; step 1.74723e-09
NoConvergence: [ 6.00000000e+01 -2.63924280e+06 -2.98161396e+04 ...  3.39690111e+08
  2.14121070e+07 -2.54528649e+08]
</code>```

So, what actually is happening here? I can only converge these algorithms when i feed the solution of fsolve into them, from which the objective function is already below the tolerance for convergence. I know that fsolve did converge, but i am just running tests for much larger system of equations, from which the large scale solvers, those above besides fsolve, are required.


问:
step size cycling in semismooth newton for convex problem
I am using backtracking linesearch to globalize a (semismooth) newton solver to minimize a (strongly semismooth) strongly convex function , and I am observing something strange (which may be a bug). For some problems, I observe that the accepted step size alternates between 0 and 1 in successive iterations. On other problems, it seems to work fine.

My question is: Is it possible for Newton method to cycle (on a convex problem) when using backtracking linesearch?


问:
Form of nonlinear diffusion equation
Consider the following nonlinear diffusion problem,

$$
\frac{\partial u}{\partial t} = x^{-2}\frac{\partial}{\partial x}\left(x^2 u^4 \frac{\partial u}{\partial x}\right), \quad 0 < x < 1
$$

We can notice that if $v = u^5/5$, we can simplify the diffusion term to

$$
\frac{\partial u}{\partial t} = x^{-2}\frac{\partial}{\partial x}\left(x^2 \frac{\partial v}{\partial x}\right), \quad 0 < x < 1
$$

The time derivative can also be replaced recognizing that

$$
\frac{\partial u}{\partial t} = u^{-4} \frac{\partial v}{\partial t} = \frac{u}{5v}\frac{\partial v}{\partial t}
$$

While the transformation eliminates the non-linearity in the diffusion term, we potentially face trouble evaluating the nonlinear mass matrix, say if the initial condition where $u = 0$.

Q1: Is there a way out of this situation? Perhaps by solving for both $u_t$ and $v_t$ (or $\dot{v}$) in time?

Q2: Does the idea of performing the semi-discretization in terms of $v = u^5 / 5$ have merit?

If possible I'd like to use the pdepe solver in MATLAB. The original form is compatible with pdepe, however I thought the solution may be more accurate if semi-discretized for $v$ instead.

Edit (additional information)

For a bit more context, this problem appears as Example A on pg. 15 in the report:

Bakker, M. (1977). Software for semi-discretization of time-dependent partial differential equations in one space variable:(preprint). Stichting Mathematisch Centrum. Numerieke Wiskunde, (52/77).

A PDF (3MB) can be downloaded here.

The full problem is to solve:

$$
\frac{\partial u}{\partial t} = x^{-2}\frac{\partial}{\partial x}\left(x^2 u^4 \frac{\partial u}{\partial x}\right), \quad 0 < x < 1
$$

given the boundary conditions,

$$
\frac{\partial u}{\partial x} = 0, \; x = 0\\
5 \int_0^1 x^2 u(x,t) dx + [u(1,t)]^5 = 1
$$

and initial condition

$$
u(x,0) = \begin{cases}0, \; 0 \leq x < 1\\
0, \; x = 1 \end{cases}
$$

In the original report the author chose to semi-discretize the right-hand side (or flux term) as a function of $v = u^5/5$. With pdepe it is not possible to do so, as the gradient $\frac{\partial v }{ \partial x}$ would only be available if $v(x,t)$ was solved for as a second PDE. (Implementing the integral is a second problem...)

Here is what the solution looks like, obtained with linear finite elements on an equidistant mesh with 80 elements:

The sharp penetrating front in the solution presents quite a challenge for numerical schemes.

The steady-state solution can be obtained from the second boundary condition, by substituting $u(x,t) = u_\infty$, which results in the non-linear equation
$$
\frac{5}{3} u_\infty + u_\infty^5 = 1
$$
with solution $u_\infty \approx 0.565$.

Title: Form of nonlinear diffusion equation


问:
Solving Laplace for Velocity Potential in Constricted Channel
I am trying to solve the 2D Laplace equation numerically to give the velocity potential of a fluid flowing in a channel with a constriction:

$u_{xx}+u_{yy}=0$

There is a constriction in the channel at $0\leq x\leq 1$, described on the lower boundary ($y=0$) by $y_l(x)=2tx(1-x)$, where $t << 1$. Similarly, at the upper boundary ($y=1$), the constriction is mirrored: $y_u(x)=1-y_l(x)$.

For $x<0$ and $x>1$, the condition at the boundary is given by $\frac{\partial u}{\partial y}=0$.

Along the lower constriction, the boundary is conditioned by $\frac{\partial u}{\partial y}-(1+\frac{\partial u}{\partial x})\frac{dy_l}{dx}=0$ (and equivalently for the upper boundary with $\frac{dy_u}{dx}$). Note, since $t<<1$, I am approximating the boundary to apply to $y=0$ and $y=1$, not on the actual $y_l$ and $y_u$ boundaries.

To start with, I have discretized the Laplace equation, setting the x and y step size to be equal ($\Delta x=\Delta y=h$) to get a stencil for each interior point of the problem:

$u_{i,j}=\frac{1}{4}(u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1})$.

Then, I discretize the $x<0$ and $x>0$ boundary conditions along $y=0$ using central finite differences:

$\frac{\partial u}{\partial y}\approx \frac{u_{i,1}-u_{i,-1}}{2h}, \frac{\partial u}{\partial y}=0 \implies u_{i,1}-u_{i,-1}=0$

Then subbing in from the stencil to eliminate $u_{i,-1}$ as it's outside the domain: $u_{i,0}=\frac{1}{4}\left(2u_{i,1}+u_{i-1,0}+u_{i+1,0}\right)$ (similarly for $y=1$ to get an expression for $u_{i,N}$).

For $0\leq x\leq 1$, I discretize the given boundary condition using central finite differences:

$\frac{u_{i,1}-u_{i,-1}}{2h}-2t(1-2x)-2t(1-2x)\frac{u_{i+1,0}-u_{i-1,0}}{2h}=0$

Making the same substitution for $u_{i,-1}$, I get:

$u_{i,0}=\frac{1}{4}\left(2u_{i,1}+u_{i+1,0}+u_{i-1,0}-2t(1-2x)\left(u_{i+1,0}-u_{i-1,0}\right)-4ht(1-2x)\right)$

and similarly for $u_{i,N}$ for the $y=1$ boundary. At each end on the x axis, I have the condition that $u_x=0$ as $x\to \pm\infty$.

Are these discretizations correct? When implementing them in Python, I get a strange plot:

Darker red means higher potential $u$, white means $u=0$. The two thick red bars indicate $x=0$ and $x=1$ respectively. I would expect my solution plot to somewhat look like the first image in this post but there seems to be a semicircle of 0 potential and of max potential where I would expect just a semicircle of 0 potential (indicating the constriction).

Note, the x and y axis units are grid point index (i,j) rather than the actual x and y values.

What is going wrong?


问:
SPH example codes for solving differential equations
I'm quite new to the smoothed particle hydrodynamics method, so I'm trying to find Matlab or Python code to help me understand it. Specifically, I want to know how to solve the one-dimensional diffusion equation like in this paper,

$$\frac{\partial u}{\partial t} = \lambda \frac{\partial^2 u}{\partial x^2}$$

SPH is reviewed in many good introductory texts, eg. Monaghan and Price papers. I understand the math/theory, but I'm trying to find a detailed simple code to see how things work numerically. I know there are software packages like DualSPHysics, but I prefer something basic for solving the above differential equation. Can you please help?


问:
Problem with the simulation for pulse propagation solving schrodinger equation using split step method
Based on this code, I was trying to work on a project for my coursework, but I am not able to understand the result. I am not sure, if i am doing anything wrong with the code but the data doesn't make sense.

```<code>#!/usr/bin/env python
# coding: utf-8

#
#
# ```
# # This is formatted as code
# ```
#
# Split-Step Method for Nonlinear Schrodinger Equation
#
#
# <span class="math-container">$i\frac{∂U}{∂z}=-\frac{sign(k^"_0)}{2L_D}\cdot\frac{∂^2U}{∂\tau^2} +  \frac{1}{L_N}\lvert U(z,\tau)\rvert^2\cdot U(z,\tau)$</span>

# In[ ]:

import numpy as np
import sympy as sp
from scipy.fftpack import fft, ifft, fftshift, ifftshift, fftfreq
import matplotlib.pyplot as plt
import math
get_ipython().run_line_magic('matplotlib', 'qt')

# Simulation Parameters

# In[228]:

c = 3e8  #m/s
lmbda=0.8                    # Wavelength in micrometer
pul_t = 10e-15  # Pulse Duration in seconds
I_0 = 1e12                   # 1TW/cm**2
nAr= 19.4e-19                #
nSi = 2.19e-16               # cm**2/W
w0 = 2*np.pi*c/(lmbda)  # Central frequency of the Pulse
chirp=0                      # Initial Chirp
L=1e3                       # in units of micrometer                    # Propagation Length (1cm)

# In[229]:

N  = 2**15       # Number of points
dt = 0.1e-15     # Time resolution [s]

# Define function for getting FFT frequency range from time

# In[230]:

def getFreqRangeFromTime(time):
    return fftshift(fftfreq(len(time), d=time[1]-time[0]))
t=np.linspace(0,N*dt,N)      # Time step array
t=t-np.mean(t)               # Center so middle entry is t=0
f=getFreqRangeFromTime(t)
w=2*np.pi*f

# In[231]:

def nderivative(x,n):
  n_val = n.subs(x,lmbda).evalf()      # Returns Value of refractive index
  d2n = sp.diff(n,x,2)                 # Evaluates second order derivative
  d2n_val = d2n.subs(x,lmbda).evalf()  # Computes value of d^2n/dlambda^2
  return d2n_val, n_val

x = sp.Symbol('x')
n_silica = (1 + 0.6961663*x**(2)/(x**(2) - 0.0684043**(2)) + 0.4079426*x**(2)/(x**2-0.1162414**2) + 0.8974794*x**2/(x**2-9.896161**2))**0.5
n_argon = 1 + 6.7867e-5 + 3.0182943e-2/(144-x**(-2))
#n_testfunc =  x**3 + 2*x**2
d2n_val1,n_val1 = nderivative(x,n_argon)
d2n_val2,n_val2 = nderivative(x,n_silica)
#dtestfunc,testfunc = nderivative(x,n_testfunc)
print(d2n_val1)
d2n_val2

# In[275]:

# Calculate Dispersion Parameter
def disp(lmbda,pul_t,d2n_val):
   k_0 = ((lmbda*10**(-6))**3)/(2*np.pi*c**2)*(d2n_val*10**12)
   l_d = pul_t**2/(np.abs(k_0))          # Calculate Dispersion Length
   return float(l_d),complex(k_0)

l_d1,k_0ar= disp(lmbda,pul_t,d2n_val1)           # l_d in meters
l_d2,k_0si= disp(lmbda,pul_t,d2n_val2)
print(l_d1,l_d2, k_0si)

# In[239]:

# Calculate Nonlinear Length in Micrometer
def Nlength(n2):
  return lmbda*10**(-6)/(2*np.pi*n2*I_0)
ln_ar=Nlength(nAr)
ln_Si=Nlength(nSi)
print(ln_ar,'micrometer',ln_Si)

# In[272]:

#Gaussian Pulse input Characteristics

def GaussianPulse(t,chirp,pul_t):
    A=1
    phi = w0*t + 0.5*chirp*t**2                #Define Phase
    E_t = A*np.exp(-2*np.log(2)*(t/pul_t)**2 )*(1+0j)
    E_t = E_t/np.max(np.abs(E_t))             #Normalisation
    return E_t,phi

E_t,phi = GaussianPulse(t,chirp,pul_t)
# #Plot in time
# plt.plot(t,np.abs(E_t)**2)
# plt.xlabel('time')
# # plt.xlim([-0.5e-10, 0.5e-10])
# plt.show()

# Pulse Spectrum
E_f=fftshift(fft(E_t))*dt
# E_f=E_f/np.max(np.abs(E_f))
plt.plot(f,np.abs(E_f)**2)
# # plt.xlim([0,5])
# plt.xlabel('freq')
# plt.show()

# SSFM

# In[274]:

#Main SSFM Loop
# gamma = 10e-3
# beta2 = 100e3*1e-30

def SSFM(E_t,E_f, l_n,k_0,l_d):
    Nz = 2**8
    dz = L/Nz
    E_t2 = np.zeros((Nz+1,N ) )*(1+0j)
    E_S = np.zeros((Nz+1,N ) )*(1+0j)
    E_S[0,:] = E_f
    E_t2[0,:]=E_t
    gamma = 1/l_n
    beta2 = np.sign(k_0/l_d)
    nonlinearity=1j*gamma*dz
    disp_and_loss=np.exp((1j*beta2/2*(w)**2)*dz )
    for z in range(Nz):
#     #     E_nl = np.exp((-1j)*1/(l_n)*np.abs(E_t2)**2*dz)*E_t2
#     #     Efft = fftshift(fft(E_nl))
#     #     argm = (-1j*np.sign(k_0)/(2*l_d)*(w-w0)**2)*dz
#     #     E_l = np.exp(argm)*Efft
#     #     E_t2 = ifft(ifftshift(E_l))

        E_t*=np.exp(nonlinearity*np.abs(E_t)**2) #Apply nonlinearity
        spectrum = fftshift(fft(E_t))*dt*disp_and_loss #Go to spectral domain and apply disp and loss
        E_t = ifft(ifftshift(spectrum))/dt #Return to time domain

        E_t2[z+1,:]= E_t
        E_S[z+1,:] = spectrum
    return E_t2, E_S

E_Ar_t, ES_Ar=SSFM(E_t, E_f, ln_ar,k_0ar,l_d1)
E_Si_t, ES_Si=SSFM(E_t, E_f, ln_Si,k_0si,l_d2)
# plt.plot(t, np.abs(E_Si_t[0,:])**2)
# plt.plot(t, np.abs(E_Si_t[-1,:])**2)
# # plt.plot(t, np.abs(E_Ar_t)**2)
# plt.xlim([-0.5e-13, 0.5e-13])
# plt.show()
# def Spectrum(E_t):
#   E_Spectrum =fftshift(fft(E_t))
# #   E_Spectrum = E_Spectrum/(np.max(np.abs(E_Spectrum)))
#   return E_Spectrum

# plt.plot(f, np.abs(ES_Ar[0,:])**2)
plt.plot(f, np.abs(ES_Ar[-1,:])**2)

plt.plot(f, np.abs(ES_Si[-1,:])**2)
</code>```

The output result in terms of the spectrum is:

I'd appreciate it if anyone could help me understand this result. The expected result is broadened spectrum with a decrease in amplitude due to dispersion. However, I don't understand the noise. Is something wrong with the code? Although the formulas and algorithm seem to be correct, the result is unexpected.


问:
Symplectic (or alike) integrator for system with Coulomb singularity and time-dependent potentials
I am trying to calculate classical trajectories for a single positive ion and a single electron inside an RF trap. Therefore, I am dealing with a two-body system that possesses:

Coulomb potential with a singularity

time-dependent potential that is periodic in time and quadratic in space (smooth, relatively weak - compared to the Coulomb). However, asymptotically it becomes dominant.

I would like to calculate the classical trajectories that start from an arbitrary point in the phase space. The time interval I am interested in, is up to 10 periods of the trap field.

I've had some limited success with Tao's symplectic integrator for charged particle in electromagnetic fields (J. Comp. Phys. 327, 245 (2016)). However, the number of time steps necessary to get converged trajectories, is quite large and I would need to perform some systematic scan of the parameter space.

There are several ideas I have in mind but I am afraid that I am going to reinvent a wheel with them. Some of them:

When the particles approach each other, the time-dependent field can be perhaps neglected or treated perturbatively. So using different integrators in different regions of the phase space could help. However, it seems from the classical references that it would take away the symplecticity.

I could regularize the COulomb singularity, add an additional degree of freedom and make the Hamiltonian non-separable. Then I could try to solve it in the extended phase space.

Maybe it would be worth trying some existing simulation software (LAMMPS?). However, how is the singularity treated there?

Is there "canonically" appropriate integrator for this kind of situations please? Is any of the strategy mentioned above worth the shot? Thanks for any hint.


问:
Numerical integration library interfacing with eigen
I am looking for a numerical integration library like this one. The examples look very appealing but I see that all test functions use very barebones C arrays.

Do you have any recommendations of libraries for numerical integration over functions of multi-variate arguments that bind easily to Eigen?

I suppose in a pinch, one could always use Eigen's .data() methods and element accessors to emulate the C arrays in their examples.

I just think it would be nice to interface with Eigen in a slightly more streamlined way.

EDIT: To make things more precise consider the following problem:

We have a known function $f(\lambda,a)$ and want to numerically integrate over $\lambda$, where $\lambda \in R^{m}$, so we are looking at a high-dimensional integral and possibly need to use Monte Carlo integration or high-dimensional quadrature  methods.

Any leads would be highly appreciated!


问:
Solving 2D Poisson Eq with Dirichlet BC's in Python
Question:

I am trying to solve the following PDE:

$$
    \begin{align*}
    u_{xx} + u_{yy}
    =
    \begin{cases}
    - \cos(x) \quad -\pi/2 \leq x \leq \pi/2, \\
    0 \quad \text{otherwise}
    \end{cases}
    \end{align*}
$$
The domain is the rectangle with vertices $(-π, 0), (-π,2), (π,0)$ and $(π,2)$. The boundary conditions are Dirichlet: $u$ must be zero on the boundary.

Code:

```<code>import numpy as np
from scipy.sparse import diags, csr_matrix
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import math

x_min = -np.pi #Left endpoint of x interval
x_max = np.pi  #Right endpoint of x interval
y_min = 0      #Left endpoint of y interval
y_max = 2      #Right endpoint of y interval

#mesh sizes and grid points
dx = 0.05
dy = 0.01
nx = int((x_max - x_min)/dx)-1
ny = int((y_max - y_min)/dy)-1
n = nx*ny
x = np.linspace(x_min+dx, x_max-dx, nx)  # Grid points in x-direction
y = np.linspace(y_min+dy, y_max-dy, ny)  # Grid points in y-direction
xg, yg = np.meshgrid(x, y)

# Create the diagonals
ones = np.full(n, 1)

a=np.ones(ny,dtype='float')
b=np.array([0.0],dtype='float')
c=np.concatenate((a,b),axis=0)
lower_diagonal = np.tile(c,nx)
upper_diagonal = np.tile(c,nx)

# Create the offsets
offsets = [0, -1, 1,ny,-ny]

# Create the sparse matrix
A = diags([-2*ones/dx**2 + -2*ones/dy**2,lower_diagonal/dy**2,upper_diagonal/dy**2,ones/dx**2,ones/dx**2], offsets, shape=(n, n), dtype=float)

# Construct the RHS
b = np.zeros((ny,nx))

for i in range(nx):

    if abs(x[i]) <= np.pi/2:

        b[:, i] = -np.cos( x[i] )

b1 = b.reshape(n,1)

#Solve the linear system using a sparse matrix solver
As = csr_matrix(A)
bs = csr_matrix(b1)
u = spsolve(As, bs).reshape(ny, nx)

#Plot the solution
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection='3d')
surf = ax.plot_surface(xg, yg, u)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u');
</code>```

Here is what my solution looks like:

https://imgur.com/a/kznhrZ5

Here's the correct solution:

https://imgur.com/a/oAWq38K

However, I can't figure out where my code is wrong. The only non-trivial thing I have done is generated the matrix A using scipy. I think my matrix is correct, though. :/ With dx = 1.2 (nx = 4), dy = 0.4 (ny = 4), my matrix is of the following form:

https://imgur.com/a/5elDvhG

I think this is correct, no?


问:
Schur complement formulation of linear system
Consider a system of the following form:

$$(A+K)x=b$$

where $A$ is symmetric, positive definite and block diagonal (in fact, a block diagonal matrix made of stiffness matrices arising from FEM discretizations at multiple timesteps, grouped together). Multiplying our linear system by $P^t=K^tA^{-1}+I$, we get $$(K+K^t + K^tA^{-1}K+A)x = P^tb$$

We assume that $K$ is non-symmetric, has full rank and that $K+K^t$ is positive definite, so that this is a symmetric and positive definite system of equations. The only inconvenient here is the presence of $A^{-1}$, and the fact that one needs to solve a (block diagonal) system to obtain $P^tb$.

Yet another way is to further rewrite everything recognizing a Schur complement:

$$
\left\lbrack
  \begin{array}{cc}
  A& -K^t\\
  -K^t& -K-K^t-A
\end{array}
\right\rbrack \left\lbrack \begin{array}{c} p\\ x\end{array} \right\rbrack=
\left\lbrack \begin{array}{c}-b\\ -b\end{array}\right\rbrack
$$

Which one of the three systems would you rather work with, and with which solver would you work? I am interested in wall-clock time and parallelizability

We assume our blocks to be large, e.g. of size $1e8-1e10$.

To connect the first two formulations, an idea might be to solve formulation one with pre-conditioned GMRES, and use $P^t$ as preconditioner. Applying $P^t$ can be done with the (preconditioned) conjugate gradient, and this is highly parallelizable since $A$ has a diagonal block structure.

But in this way, since we are solving essentially formulation two, we'd still need to come up with a good preconditioner for it.

Edit. Also because formulation two seems much worse conditioned than formulation one, as numerical experiments show.


问:
Solving 2D Poisson Eq with mixed BC's in Python
I am trying to numerically solve the Poisson's equation

$$
u_{xx} + u_{yy} = - \cos(x) \quad \text{if} - \pi/2 \leq x \leq \pi/2 \quad  \text{0 otherwise}
$$

The domain is the rectangle with vertices $(-π, 0), (-π,2), (π,0)$ and $(π,2)$. The boundary conditions are mixed:

Dirichlet: $u(x,0) = 0$
Neumann: $u_y(x,2) = 0$
Periodic: $u(-\pi,0) = u(\pi,0) = 0$.

I have the code written, but I can't figure out why my code isn't plotting the correct solution.

My solution:

https://imgur.com/a/fnDEvvl

Correct solution:

https://capture.dropbox.com/4hzsLwatNsc05Fhl

Here is my code:

```<code>import numpy as np
from scipy.sparse import diags, csc_matrix
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import seaborn as sns
import math

x_min = -np.pi #Left endpoint of x interval
x_max = np.pi  #Right endpoint of x interval
y_min = 0      #Left endpoint of y interval
y_max = 2      #Right endpoint of y interval

nx = 50                                     # Number of grid points in x
ny = 50                                     # Number of grid points in y
dx = (x_max - x_min)/nx                    # Spacing in x direction
dy = (y_max - y_min)/ny                    # Spacing in y direction
n = nx*ny                                  # Dimension of system
x = np.linspace(x_min, x_max-dx, nx)       # Grid points in x-direction
y = np.linspace(y_min+dy, y_max, ny)       # Grid points in y-direction
xg, yg = np.meshgrid(x, y)

# Create the diagonals
ones = np.full(n, 1)

a=np.ones(ny-1,dtype='float')
b=np.array([0.0],dtype='float')
c=np.concatenate((a,b),axis=0)
upper_diagonal = np.tile(c,nx)

d=np.ones(ny-2,dtype='float')
e=np.array([2.0,0.0],dtype='float')
f=np.concatenate((d,e),axis=0)
lower_diagonal = np.tile(f,nx)

# Create the offsets
offsets = [0, -1, 1, ny, -ny, (nx-1)*ny, -(nx-1)*ny]

# Create the sparse matrix
A = diags([-2*ones/dx**2 -2*ones/dy**2,lower_diagonal/dx**2,upper_diagonal/dx**2,ones/dy**2,ones/dy**2,ones/dy**2,ones/dy**2], offsets, shape=(n, n), dtype=float)

# Construct the RHS
b = np.zeros((ny,nx),dtype=float)

for j in range(nx):

    #for i in range(ny):

    if abs(x[j]) <= np.pi/2:

        b[:,j] = -np.cos( x[j] )

b1 = b.reshape(n,1)

#Solve the linear system using a sparse matrix solver
As = csc_matrix(A)
bs = csc_matrix(b1)
u = spsolve(As, bs).reshape(ny, nx)

#Plot the solution
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection='3d')
surf = ax.plot_surface(xg, yg, u)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('u');
plt.show()
</code>```

The only mistake I can think of is in the construction of the matrix $A$. Because I have been able to write a similar code where I can solve a number of BVPs with all Dirichlet conditions. So I must not be taking into account the Neumann and periodic boundary conditions correctly into account. However, I think the matrices I generate are correct. For instance, here's the heatmap of the matrix when $n_y = 7, n_x = 4$:

https://imgur.com/a/ep4Qun8

For simplicity, I have shown the matrix prior to having divided each entry by $\delta x, \delta y$. The distribution of $-4, 1, 2$'s seems right to me, so I don't know what's going wrong. Suggestions?


问:
why does the average temperature vary a lot from my setting temperature in nve?
I set my t_start and t_stop as 0.52 (lj unit), and I ran the nvt first. after it reached equilibration I ran the nve but found the temperature of the last 1000steps went up to 0.56 while my coworker still had the normal result around 0.52. we even used the same code but our results are quite different. why does the average temperature vary so much from the one I set before?


问:
Poisson equation solution in a semiconductor structure
I am trying to solve the $\textbf{1-D}$ Poisson equation for a semiconductor structure at equilibrium (There is no external bias applied).

$\textbf{Background}$

\begin{equation}
\frac{d^2V}{dx^2} = -\frac{\rho(V)}{\epsilon}\\
\rho(V) = q(N_D-N_A+p(V)-n(V))
\end{equation}

When there are two different semiconductors placed in contact, there will be a redistribution of charge, causing a dependence of charge on potential (the bands will bend).

$\textbf{Solving with the Newton-Rhapson method}$

I am trying to solve the above non-linear equation using the newton-rhapson method. Taking the central difference, the second derivative becomes:
\begin{equation}
\frac{d^2V}{dx^2}|_i = \frac{(V_{i+1}-2V_i+V_{i-1})}{dx^2}
\end{equation}

The issue with solving this now is that this derivative cannot be computed for the two end points of the structure. What are the boundary conditions to be applied on the edges?

Another issue I'm facing is when computing the charge densities by including a non-parabolic coefficient $\alpha$ (for non-parabolic bands), I have to compute these integrals in every iteration -
\begin{equation}
n = \int^{\infty}_{Ec}F_{1/2}(E)g_n(E)dE \hspace{10pt}p= \int^{E_v}_{-\infty}F_{1/2}(E)g_p(E)dE
\end{equation}
where $F_{1/2}$  is the 1/2nd order Fermi-distribution and $g(E)$ is the density of states.
I have to compute the integrals with $E_c\text{ and }E_v$, but those are unknown. How do I initiate and update their values? How are each of these energies ($E_c,E_v,E_f$) related to each other and the potential of the device?