问:
How to find two points within defined region in this constrained optimization problem?
I am doing a project related to robotics where I am using fmincon function from matlab to minimize the distance between the points x_0 and x_f. Basically task is to make manipulator move from x_0 to x_f such that it never goes out of torus region. So I think it is a constrained optimization problem.
I am having difficulty in figuring out what is the strategy to find the two intermediate points (as shown in figure below) so that given x_0 and x_f, it gives us two points [x1,y1,z1] and [x2,y2,z2] to use as via points.

Can you please figure out how should I write the objective function.
Code so far is below:

```<code>    function [x, fval, history] = solver( x0, xn, R, r )

% x0 = initial position, xn = final position, R = distance from the center of the tube to the center of the torus
% r is the radius of the tube.
% plotting torus
th=linspace(0,2*pi,60); % e.g. 60 partitions along perimeter of the tube
phi=linspace(0,pi,60); % e.g. 60 partitions along azimuth of torus
[Phi,Th]=meshgrid(phi,th);
x=0+(R+r.*cos(Th)).*cos(Phi);
y=0+(R+r.*cos(Th)).*sin(Phi);
z=0+r.*sin(Th);
s=surface(x,y,z);
s.EdgeColor='none'
axis equal;
daspect([1 1 1 ])
hold off
alpha(s,0.2);
% defining constraints and objective function
history = [];
options = optimset('OutputFcn', @myoutput);
[x, fval] = fmincon(@(x)objfun(x,xn),x0,[],[],[],[],[],[],@noncol,options);

    function stop = myoutput(x,optimvalues,state)
        stop = false;
        if isequal(state,'iter')
            history = [history; x];
        end
    end

    function z = objfun(x,p)
        z = sqrt((p(1)-x(1))^2+(p(2)-x(2))^2+(p(3)-x(3))^2);
    end
    function [y, yeq] = noncol( x ) % x ==== [ x(1)   x(2)    x(3) ]
        y = (sqrt(x(1)^2+x(2)^2)-R)^2 + x(3)^2-r^2;
        yeq =[];
    end
end
</code>```

<img src="https://i.stack.imgur.com/yJnX7.jpg" alt="Aim is to go from $x_0$ to $x_f$ within this toroid region">


问:
How to implement register blocking for 3D finite-difference stencil computations
I am in need to optimize the performance of the 3D solver that integrates a system of equations for seismic wave propagation.
Unsurprisingly, the function that implements the finite-difference stencil computations, is the bottleneck of the computations.

I read this article, where in Section 3.2 the authors say that the register blocking allowed them to obtain substantial speedup of the 3D stencil computations.

I do not understand how to implement this technique.
Is it possible to find somewhere the description of register blocking in clear form?

ContentLicense: CC BY-SA 4.0


问:
Given a list of intervals, find region that is contained by the largest number of those intervals
Start with 1d case. Say I have lots of 1d intervals $[s_i, e_i]$ and I want to find an interval $[s^*, e^*]$ to maximise the count of interval $i$ such that $[s_i, e_i]\supseteq [s^*, e^*]$.

1d case is easy. All I need to do is to sort the $s_i$ and $e_i$ together, eg $[s_1, s_2, e_1, s_3, e_3, e_2]$. Then go from left to right, $+1$ if it's $s$, $-1$ if it's $e$. And locate the maximum.

How can I deal with the case in 3d? Here intervals are 3d balls with L1 distance from a point $\{(x,y,z): |x-x_i|+|y-y_i|+|z-z_i|\leq d_i\}$.

I can cast it as binary linear programming, but I don't think it's solvable with 1000 binary variables plus 6000 auxiliary variables to deal with absolute value. I wonder if there is a polynomial time solution.

ps: I don't need the entire region, any point inside the max region is good.


问:
Find points that form curve segments in a point cloud
I have an ultrasonic image represented by a point cloud as shown below:

I need to extract features from this image. These features can be easily recognised by a human brain since they form continuous lines/curves. I am currently using a modified density based clustering algorithm to extract the features. A result is shown below:

where points inside those black borders are extracted.

However, the result is not perfect. If I zoom in at one of the clusters, you can see the points actually form two curves/lines in the cluster, one inside the red rectangle and the other inside the blue rectangle.

For this image, I can apply some different parameters for db-clustering algorithm locally to sub-divide the cluster but if I were to process hundreds of images, this became impossible.

Hence, I am asking if there is any algorithm that can recognise the lines/curves in a point cloud like such.

Update:

Link to sample data


问:
Library for solving multidimensional (n > 3) hyperbolic PDE systems
Does there exist a library (in any programming language) for solving (numerically) systems of multidimensional first-order linear PDEs in the form

$$\mathbf{u}_{t}+\hat{A}(\mathbf{x})\mathbf{u}_{\mathbf{x}}=\mathbf{f}(\mathbf{u}),$$
$$\mathbf{u}(\mathbf{x},0)=\mathbf{g}(\mathbf{x}),$$

where $\mathbf{u}=(u_1(\mathbf{x}),\ldots,u_m(\mathbf{x})),\mathbf{x}=(x_1,\ldots,x_n),\mathbf{f}(\mathbf{u})=(f_{1}(\mathbf{u}),\ldots,f_{m}(\mathbf{u})),$ and $\hat{A}(\mathbf{x})$ is a matrix of variable coefficients?


问:
Implementing boundary condition
I'm studying the transport of species A in the blood vessels,

$$\frac{\partial C}{\partial t} =  D\frac{\partial ^2 C}{\partial x^2}-v\frac{\partial C}{\partial x}$$

At x=0, I want to use the Neumann boundary condition

$\frac{\partial C}{\partial x} =$ constant

From the literature, I have the average concentration of the species A in blood.
But, I am not sure what would be a reasonable value of the constant.

I would like to ask for suggestions on how to define the concentration gradient at the inlet.

Title: Implementing boundary condition


问:
function over conditional probability
I need to create a scoring model out of estimated conditional probability functions for two events, A and B. Let 0.5 be the threshold value. Ideally, the probability is in the interval $[0,0.5)$ for A and in the interval $[0.5,1]$ for B. However, my probabilities are estimated with errors. Therefore, for a certain number of events A the probabilities fall into the interval $[0.5,1]$, and for B into $[0,0.5)$.
Let the minimum score min_score $=0$ ($P = 0$) and maximum score `max_score' $= 10$ ($P = 1$).

I need to map those probabilities to scores with a function that would emphasize and enlarge difference in probabilities of A and B in terms of scores (for example, a sigmoidal function). That is the higher the probability (likely event B) the higher the score, the lower the probability (likely event A), the smaller the score. Though for the intersection interval, which I can identify as $[a,b]$ the scores should be around (max_score + min_score)/2. In the interval $[t_1,t_2]$, which is symmetric about $1/2$, scores should be smoothing, while in $(0, t_1)$ and $(t_2, 1)$, the scores should be "more extreme" with a higher confidence.

The problem is $[a,b]$ can be spanned over $[0,1]$, and that effect should be smoothed. Can I apply some kind of regularization here? What function would I use and how would I establish the coefficients?

ContentLicense: CC BY-SA 3.0

Title: function over conditional probability


问:
Does quantum espresso and VASP use same same self-consistent-field procedure?
The codes Quantum Espresso (QE) and Vienna Ab initio Simulation Package (VASP) both use plane wave basis sets and psuedopotentials. Most of the codes in both implementation of DFT uses Fortran code. Do they use similar self-consistent-field procedure to find the solutions of Kohn-Sham equations? </p>


问:
Is this a form of stochastic gradient descent?
I want to minimize the following with respect to parameters $B$.

$$\sum_{k = 1}^{K} f(A_{k}, B)$$

where $A_k$ are $K$ different data-sets and $B$ is a matrix of parameters.

Can I do this by a gradient descent where:

randomly pick one matrix $A_{i}$ among $K$ matrices.

perform a gradient step that minimizes $f(A_{i}, B)$

repeat steps 1 and 2 until convergence.

The function $f$ is the logarithm of a linear function.

After searching online, this looks like a stochastic gradient descent.
Is this really an SGD?


问:
Quadrature methods for peaky integrands
Consider
$$
I = \int_{-L}^L f(x)dx,
$$
where $f(x)$ is real-valued and analytic on $[-L,L]$, but it has a pole in the complex plane whose real part lies in $[-L,L]$. Call it $z_0$, and assume it is a simple pole. If the imaginary part of $z_0$ is small compared to the real part, then $f(x)$ is sharply peaked in the vicinity of $x_0=\mathrm{Re}\{z_0\}$. If this peak is sufficiently sharp and $f(x)$ is slowly varying elsewhere, then basic adaptive quadrature methods (e.g., recursive Simpson or Gauss-Legendre) will spend considerable effort attempting to resolve the peak.

An obvious thing to try is to approximate the value of $x_0$ and break the integral up
$$
I = I_-+I_+ = \int_{-L}^{x_0}f(x)dx+\int_{x_0}^Lf(x)dx
$$
But doing so only pushes the problem to the endpoints of the integration domain instead of the interior. I wonder if there is a standard change of variables to alleviate this behavior, though?

Another option could be to try to evaluate $I$ by means of contour integration in the complex plane and use of the residue theorem. However, this will not fly if the software implementation of $f$ is only defined for real arguments.

Are there other specialized quadrature methods for these kinds of "peaky" integrands? Ideally I'd like to make no further assumptions on $f$ execpt that it's real, analytic, and slowly varying away from $x_0$. But I am also open to answers that require more assumptions on $f$ if there is no general answer.

Title: Quadrature methods for peaky integrands


问:
Block matrix and DSYRK
I want to compute the matrix
$$
A = \sum_{i=1}^N v_i v_i^T
$$
where each $v_i$ is a given vector of length $2500$, so that $A$ is $2500 \times 2500$, and my $N$ is about 2 million. Rather than call DSYR 2 million times, I found I can speed up the calculation significantly, by storing the $v_i$ in block matrices, and then calling DSYRK when the block matrix is full. I define a $\textrm{block_size} \times 2500$ matrix $W_k$ as:
$$
W_k = \pmatrix{ v_1^T \\ \vdots \\ v_b^T }
$$
so that
$$
A = \sum_{k=1}^{\textrm{nblocks}} W_k^T W_k
$$
The reason I set the rows of $W_k$ equal to the vectors $v_i$ and not the columns is because I use C and this way the incX of the rows will be 1 which I think helps to speed up the BLAS operations.

I am finding this procedure is still slower than I would like, and my question is how do I choose the $\textrm{block_size}$ (i.e. number of rows of $W_k$) to get the best performance?

Currently, I am choosing $\textrm{block_size} = 50000$, which I only chose so that the total size of the $W_k$ storage matrix is 1GB ($50000 \times 2500 \times 8 = 1 GB$). Is there a better choice of $\textrm{block_size}$ (number of rows of $W_k$) to get optimal performance?

In case it matters I am using the ATLAS CBLAS library on a machine with 24 cores.

Title: Block matrix and DSYRK


问:
Best way to numerically compute elliptic integrals of the third kind with complex arguments?
I need to compute elliptic integrals of the third kind with complex arguments, preferably in C++. Is there code out there to do this? I have discovered the Arb library, but that does much more than I would like to do, and would be much slower for my particular application.</p>


问:
Computing a Flux Integral in Paraview
I am currently looking into post-processing of simulation data using Paraview. I would like to compute certain integrals of field quantities. As an example, consider the following surface integral of an electric field in a simple situation:

$$\iint_{S}^{} \vec{E}\cdot\vec{\omega}=\Phi $$

The surface $S$ is shown in the following picture:

I can define a surface in ParaView and I can compute the surface normals $\vec{w}$ (white arrows in he picture). In order to compute the integral in question, I guess I would have to use a calculator to compute $\vec{E}\cdot\vec{\omega}$ and then use the "Integrate Variables" Filter.

My problem is, that I can only add the calculator to the "case0001.vtu tree" or the "plane tree" in the pipeline browser such that it's impossible to select both the vectors "electric fields" and "normals" in the calculator to compute the integrand $\vec{E}\cdot\vec{\omega}$.

Does anyone know how I could use a calculator to compute the integrand in Paraview?

Title: Computing a Flux Integral in Paraview


问:
How to simulate water, falling under gravity, and impinging on a curved surface, which is kept/present in a domain, containing air?
TL;DR: How do I simulate a hole, at the bottom of a (full) water tank?

I am attempting to simulate water, flowing out of a hole/slit, at the bottom of a tank (Water Domain) (under the influence of gravity); entering into a domain (vertically downward), that contains air (Air Domain), and impinging on a curved surface. I am using COMSOL Multiphysics 5.4/5.3a. The idea is to simulate (in 2D) and obtain data on:

1. The distribution of water, as it flows along the surface, on two sides of the surface (for a 2D Simulation), for different initial points of impingement.

2. The (distribution of) stress produced, due to the impinging water, on the (different points of the) curved surface.

I have looked through a plethora of papers, as well as tutorials and none of them include, “liquid/gas, flowing out, under gravity, into another domain, containing a different liquid/gas”. I am assuming Laminar Flow and No Slip Boundary Condition, at the hole & walls and the curved surface.

Here, as far as I know, I cannot define an Inlet at the interface of the two domains. I have tried Multiphase (Phase Field), with no inlet and 1 outlet in the Air Domain (situated vertically downwards), but it did not converge. A Surface Plot of Fluid Velocity showed water, entering the Air Domain, but after a few more time steps, COMSOL threw a Did not converge error, thereby preventing the water from reaching the curved surface. I have tried the following, hoping to fix the issue:

Turning the “hole”, into a faucet or nozzle (of sorts), with a smooth downward curve: No change in outcome observed - Still not converging.

Defining a Pressure Point Constraint (1 atm) at the top of the Fluid Domain - No change in outcome observed - Still not converging.

Any help or hint is much appreciated. Also, if this problem (First paragraph - in Bold text) can be done in any other way (preferably in COMSOL), I’d like to hear that, as well.


问:
How to compute large condition number of a matrix in Python?
I have a matrix that is extremely singular, but I am still interested in computing the exact condition number, which is the ratio between the largest and smallest singular values.

Is it possible to do it in Python? I am curious because if we stick with machine precision 10^-16, then any condition number beyond 10^16 could be inaccurate.

OwnerDisplayName: ArtificiallyIntelligence


问:
Hydrogen-like wavefunction as starting guess for atomic solver?
I've been looking into radial solvers for quantum wave equations (Schroedinger and Dirac). In both cases, the suggestion seems to be to go with the "shooting method", with integration schemes of various precision such as the Numerov one.

Now, for a shooting method, one needs to impose starting values coming from boundary conditions. Many texts mention asymptotic behaviours such as Bessel functions near $r=0$ etc. Since my interest is in solving perturbed versions of the basic Coulomb potential, and the solution for this problem is already known,  however, couldn't I just start with the known hydrogen-like solution and work from there? Or is there some reason why this is not indicated?


问:
Can I take advantage of a nearly banded A in AX=b?
I am working on a 1D drift-driffusion problem in a finite-difference (FD) approach.

I hade 3 equations per node ($3N$ in total): electron continuity $E_i$, Poisson $P_i$, hole continuity $H_i$. With corresponding variables: electron concentration $n_j$, potential $V_j$, hole concentration $p_j$.

By rearranging the problem has a banded matrix, the Jacobian $A$ takes the following form, with shape $(3N\times 3N)$

$$A = \begin{pmatrix}
A_{1,1} & A_{1,2} &         &         &         &         &       \\
A_{2,1} & A_{2,2} & A_{2,3} &         &         &         &       \\
        & A_{3,2} & A_{3,3} & A_{3,4} &         &         &       \\
        &         & \ddots  & \ddots  & \ddots  &         &       \\
        &         &         & A_{N-2,N-3} & A_{N-2,N-2} & A_{N-2,N-1} &       \\
        &         &         &         & A_{N-1,N-2} & A_{N-1,N-1} & A_{N-1,N} \\
        &         &         &         &         & A_{N,N-1} & A_{N,N} \\
\end{pmatrix}$$

where

$$
A_{i, j} = \begin{pmatrix}
{\partial E_i}/{\partial n_j} & {\partial E_i}/{\partial V_j} & {\partial E_i}/{\partial p_j} \\
{\partial P_i}/{\partial n_j} & {\partial P_i}/{\partial V_j} & {\partial P_i}/{\partial p_j} \\
{\partial H_i}/{\partial n_j} & {\partial H_i}/{\partial V_j} & {\partial H_i}/{\partial p_j}
\end{pmatrix}
$$

Owing to that data structure, I have been using scipy.linalg.solve_banded to solve for $X$ in

$$AX = b$$

where $b$ is a $(3N \times 1)$ shaped vector.

My problem comes from that I want to add an equation that controls the potential at one end (say on the left) based on the value of total fluxes of $n$ and $p$ at both ends.

As a result the additional equation depends on $(n_1, V_1, p_1, n_2, V_2, p_2)$ but also $(n_{N-1}, V_{N-1}, p_{N-1}, n_N, V_N, p_N)$.

Therefore, if I was to include this new equation in $A$ to form a $(3N+1 \times 3N+1)$ Jacobian $A'$, the bandwith would suddenly become much larger. Though I can optimize this by inserting the new equation in the middle row of $A$, the banded representation would remain very sparse.

My question is then: is there a way I can solve for $X'$ in
$$A' X' = b'$$

while keeping the advantages of the banded structure of $A$?

EDIT addressing comments

Here is a visual example of the Jacobian before (left) and after reordering (right). The last column/line, highlighted in yellow, is the additional equation that fails to fit the banded structure.
Elements in dark gray background are non-zero, elements in white are zeros, elements in light gray background are zero, but a future version of my physical could in principle make them non-zero so an ideal solution would not exclude that possibility.

(click to enlarge)

Disclaimer: I do not pretend to be the author of this ordering method. I have read it a few years ago in a software manual I believe ...


问:
Speeding up the solution of a large set of nonlinear algebraic equations in `sympy`
I have a quite large algebraic equation system to solve, the system is so large, I can't post the example here, so I am posting it to pastebin.

The sympy.solve is taking ages: I started it this morning and it is still running on a single core. I have tried to use the optional arguments from the sympy documentation, like quick=True and warn=True, to speed things up and get warnings if the solver fails, but there was no warning given, and the solver is still dying.

Is there a way multiprocessing can be used or maybe mpi4py with the solve function in sympy. I have access to HPC resources, but I'm not sure how solve can be separated into different MPI processes.

Is there a way in sympy to get an idea if the system is solvable or if the solver is diverging? Something like runtime solver information, but without running in debug mode?


问:
WENO5 scheme in a staggered grid
I'm trying to use the finite-difference WENO scheme to solve the 2D density conservation law with axial symmetry (coordinates $r,z$):

$\frac{\partial \rho}{\partial t}+\nabla \cdot (\rho \vec{v}) = \frac{\partial \rho}{\partial t}+ \frac{1}{r}\partial_r (r \rho v_r)+\partial_z (\rho v_z) = 0$.

$\rho$ is defined at grid points $i,j$, $v_r$ at $i+1/2,j$ and $v_z$ at $i,j+1/2$.

Ideally, I would like use a WENO reconstruction procedure to estimate $f_r = r \rho v_r$ at $i\pm 1/2,j$, and $f_z = \rho v_z$ $i,j \pm 1/2$, from their values at $i,j$ (flux-splitted), and use those estimates to update the density at $i,j$:

$\frac{\partial \rho^{i,j}}{\partial t} = -\frac{1}{r}\frac{1}{dr} \left[f_r^{i+1/2,j}-f_r^{i-1/2,j} \right]-\frac{1}{dz} \left[f_z^{i,j+1/2}-f_z^{i,j-1/2} \right]$

Now, since the velocities aren't defined at $i,j$, I use a simple interpolation (i.e, $v_r^{i,j} = (v_r^{i+1/2,j}+v_r^{i-1/2,j})/2$) to estimate their values at those points.

Doing this, however, does not yield good results. I've tried to do a WENO interpolation of the densities, but that also isn't working as well as I need it to.

Any suggestions on how to handle this issue?

Additionaly, any recomendations on how to properly flux-split in this case? Most flux-splitting procedures seem to reduce to a trivial case when $f(\rho) \propto \rho$.

Title: WENO5 scheme in a staggered grid


问:
Unusual boundary conditions on Matlab
I'm trying to solve the following PDE by Matlab,

$$
u_t-\Delta u = 0, \quad \text{in}\quad \Omega\times (0,T)
\tag{1}
$$
$$
u_t-\Delta_\Gamma u + \partial_\nu u=0,\quad \text{on}\quad \Gamma\times(0,T)
\tag{2}
$$

where $\Omega$ is a bounded domain of $\mathbb{R}^n$, $\Gamma =\partial \Omega$ is the boundary of $\Omega$, $\partial_\nu$ is the normal derivative, and $\nu$ is the outer unit vector.

For simplicity, I considered the one dimensional case when $\Omega$ is the interval $(0,1)$, and discretized the equation following the Finite Difference method (we can neglect the Laplace-Beltrami operator $\Delta_\Gamma$ in this case), but I didn't find how to solve such PDE with the dynamic boundary conditions (2) in a simple way using pdepe or PDE toolbox. Usual Matlab functions use the usual boundary conditions (Dirichlet or Neumann), but the dynamic boundary conditions are not included. So I'm asking how to solve this PDE with Matlab in a simple way.
For the 2D heat equation can we define $\Delta_\Gamma$ by Matlab tools.

Title: Unusual boundary conditions on Matlab


问:
Correct order of convergence

I have a sequence of points which was obtained from an iterative algorithm, and I computed the order of convergence $p$ of the method using the formula

$$
p \approx \frac{\log({\rm err}(k+2))-\log({\rm err}(k+1))}{\log({\rm err}(k+1)) - \log({\rm err}(k))}.
$$

The first computed ratio is much larger compared to the other ones. Without removing this value from the complete list of values, I obtained $p = 1.2615$ as the order of convergence of my algorithm. However, when I computed the order of convergence after removing the outlier (the first calculated ratio), I got $p = 1.0495$ (see figure).

My question is, which of the two is the correct order of convergence of my algorithm?

Title: Correct order of convergence


问:
SDE solver in python: manual determination of integrator step size (dt)
Aim: I am trying to solve a system of SDEs, while using the SDEint package in python 3.x. It is a system of SDEs adapted from and inspired by the Zombie Apocalypse model, published by Munz, 2009. I tried to add an explanation of the terms within the code.

Problem: The problem I would like to solve occurs in the example below (overflow encoutered in double_scalar). Note that the noice term is currently 0, so it is basically a system of ODEs and not yet a system of SDEs.

```<code>import matplotlib.pyplot as plt
import numpy as np
import sdeint
from scipy.integrate import odeint

p, f, e, k = 0.7, 0.0002, 0.05, 0.02    # human birth rate, confilt occurence rate, conflict end rate, zombie killing rate

tspan = np.linspace(0, 20., 100)
y0 = np.array([100000., 0., 0., p])

def ff(y, t):
    Hi = y[0]
    Ci = y[1]
    Zi = y[2]

    # Human are created (exponential)                       --> y[3] * Hi
    # Humans and zombies engage in conflict at rate (f)     --> - f * Hi * Zi
    # the conficts end at rate (e) without a victor         --> + e * Ci
    f0 = y[3] * Hi - f * Hi * Zi + e * Ci

    # Humans and zombies engage in conflict at rate (f)     --> f * Hi * Zi
    # the conficts end at rate (e) without a victor         --> - e * Ci
    # and zombies kill humans in conflict at rate (k)       --> - k * Ci
    f1 = f * Hi * Zi - e * Ci - k * Ci

    # Zombies arise from the earth at a fixed rate          --> 10**7
    # they arise as victor from a conflict                  --> + k * Ci
    # enter conflicts                                       --> - f * Hi * Zi
    # and exit conflicts without killing                    --> + e * Ci
    f2 = 10**7 + k * Ci - f * Hi * Zi + e * Ci

    f3 = 0
    print(f0, f1 ,f2, f3)
    return np.array([f0, f1, f2, f3])

def GG(y, t):
    return np.diag([0, 0, 0, 0])

result = sdeint.itoint(ff, GG, y0, tspan)

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_yscale('log')
ax.set_ylim(bottom=1, top=10**13)
plt.plot(result)
plt.show()
</code>```

I have tried to solve the same system of ODEs (without noise) using Scipy's odeint and this seems to work well due to the possibility of adding several arguments to odeint.

```<code>result = odeint(ff, y0, tspan, mxstep=10 ** 9, rtol=10 ** (-3),
                atol=10 ** (-3), hmax=1)
</code>```

So I thought maybe SDEint uses too large integration steps, however there is no possibility to manually add additional arguments to sdeint.

Potential Solution(?): So I was wondering, if a solution might be to manually set the set size of the integration. I tried to do so by increasing the amount of steps in np.linspace e.g.

```<code>tspan = np.linspace(0, 20., 100000) # 100 --> 100000
</code>```

However, this does not seem to work properly.

Question: Since I am not really familiar with ODEs/SDEs, I am not sure if I am addressing this problem properly. Is this a common way to address the problem of too large integration steps? And might the solution be to increase the step size or am I completely mistaken?

Any help would be greatly appreciated!


问:
Numerical analysis, pivoting and incomplete LU decomposition
When doing LU decomposition, the algorithm will break down if any of the diagonal element $x_{ii}$ is zero.  Therefore, we can use pivoting on the matrix such that $x_{ii}$ is no longer zero. That is instead of looking at the $x_{ii}$ we look at another element $x$ in the matrix.

One problem with LU decomposition is that when our matrix $A$ is sparse, we would like to keep the sparsity pattern of $A$, something that will not happen when we do LU decomposition. Therefore, we can use incomplete LU decomposition, where if $x_{ij}$ of $A$ is zero, we also skip $x_{ij}$ of $L$ or $U$ depending on where $x_{ij}$ is in $A$.

My question is: will we do pivoting for incomplete LU decomposition? If we do, how do we do it? If we don’t, how do we make sure that our algorithm doesn't break?


问:
Piecewise-linear Continuations vs Marching Squares/Cubes
It seems that both piecewise-linear continuation and marching squares are methods to produce iso-contours of a scalar function given the function's values on a grid. It seems that piecwise-linear continuation works for any triangular mesh, not just a grid.  Marching squares also has a variation which works on triangles.

My questions are:

What is the connection between these methods?

They both seem to use linear interpolation to position the contour at the edges of the grid. Are the two methods equivalent under some conditions?


问:
Number of $S_n$-orbits in $P^k(\{1,\dots,n\})$
This is a particular case of a question I asked on Mathematics Stackexchange, question which got no answer so far.

Let $n$ and $k$ be integers with $n\ge1$, $k\ge0$, and let $a(n,k)$ be the number of orbits of the symmetric group $S_n$ on the $k$-th iterated power set
$$
P^k(\{1,\dots,n\})
$$
of the set $\{1,\dots,n\}$.

We have $a(n,0)=1$, $a(n,1)=n+1$, $a(1,k+1)=2^{a(1,k)}$ and
$$
\frac{b(n,k)}{n!}\le a(n,k)\le b(n,k),
$$
where $b(n,k)$ is the cardinality of $P^k(\{1,\dots,n\})$.

We also have $a(2,2)=12$. Indeed, the twelve orbits of the two-element group $S_2$ in $P^2(\{1,2\})$ can be described as follows.

The set $P^2(\{1,2\})$ having sixteen elements, it suffices to give the eight fixed points. The complement of a fixed point being a fixed point, it suffices to give four fixed points $F_1,F_2,F_3,F_4$ such that, for all $i,j$, the subset $F_i$ is not the complement of $F_j$. A possible choice of the $F_i$ is
$$
\varnothing,\ \{\varnothing\},\ \{\{1,2\}\},\ \{\varnothing,\{1,2\}\}.
$$

Can anybody compute, say, $a(2,3),a(3,2)$ and $a(3,3)$?

If you don't feel like computing these three integers, I'd be most grateful if you could nevertheless explain to me a simple (perhaps not efficient) way of computing them.

Sadly, I know nothing about programming languages.

Edit 1. Here are more details about the equality $a(2,2)=12$. We have
$$
P(\{1,2\})=\{\varnothing,\{1\},\{2\},\{1,2\}\}.
$$
So
$$
P^2(\{1,2\})=P(\{\varnothing,\{1\},\{2\},\{1,2\}\})
$$
has $16$ elements, one of cardinality $0$ (that is, with zero elements), which is the empty set
$$
\varnothing,
$$
four of cardinality $1$, which are
$$
\{\varnothing\},\{\{1\}\},\{\{2\}\},\{\{1,2\}\},
$$
six of cardinality $2$, which are
$$
\{\varnothing,\{1\}\}, \{\varnothing,\{2\}\}, \{\varnothing,\{1,2\}\},
$$
$$
\{\{1\},\{2\}\}, \{\{1\},\{\{1,2\}\}, \{2\},\{\{1,2\}\}
$$
four of cardinality $3$, which are
$$
\{\varnothing\}^c,\{\{1\}\}^c,\{\{2\}\}^c,\{\{1,2\}\}^c,
$$
where $A^c$ denotes the complement of $A$ in $P(\{1,2\})$, and one of cardinality $4$, which is
$$
P(\{1,2\}).
$$
Now let $X$ be one of the sixteen above expressions.

If you swap the symbols $1$ and $2$ in $X$, you get either $X$ itself, or a different expression $Y$, and if you perform the same operation $Y$, you get $X$ back.

So the non-invariant expressions (the "non-fixed points") come in pairs.

By inspection one sees that there are $8$ invariant expressions and $4$ pairs of non-invariant expressions, that is $12$ "orbits".

Edit 2. An inductive formula for $a(2,k)$ is given in an edit to this Mathematics Stackexchange post, but I'm still unable to compute $a(3,2)$.


问:
Open Source Packages Implementing Continuous Wavelet and Scaling Functions
I'm looking for an open source software package that provides a fast evaluation of continuous Daubechies/Symmlet wavelet/scaling functions.

GSL only has the discrete wavelets, and PyWavelets comes close to what I want, in that it gives a gives the values of the functions on an equispaced grid of whatever refinement you desire.


问:
Simplification of an optimization objective
Let $G(V,E)$ is a weighted simple graph, where $V$ and $E$ are the set of vertices and Edges. The graph is undirected. Let $A \in \{0,1\}^{n\times n}$ and $W \in R_+^{n\times n}$ be the adjacency matrix and the weight matrix of the graph. Note that both both $A$ and $W$ are symmetric. My target is to output a sparse representation of this graph. I am focusing on edge sparsification that means deleting edges from the graph. In particular I am concentrating on degree constraint subgraph of the original graph as a representation of sparsity. Let the new graph is $G_s(V,E_s)$, where the $E_s$ is the edge set of the new graph while the vertex set remains unchanged. My goal in sparsification is to find a degree constraint subgraph $G_s$ of the original graph $G$, optimizing some distance criteria between them.

A vertex $v \in V$ in $G$ can be thought of a vector, say $W_v$ in $R_+^{|V|}$. Here $W_v$ is the $v$'th row of the matrix $W$ which represents the weights of the adjacent edges of $v$. We define a binary matrix $X \in \{0,1\}^{n\times n}$. $X$ would be the adjacency matrix of the new graph $G_s$. We want to solve the following optimization problem.

\begin{align*}
 \text{min} \quad
    &\sum_{v \in V}||W_v-W_v \circ X_v||_2  \\
    \text{s.t.,}\quad
    &\sum_{i \in V} X_{vi} \leq b_v  \quad \forall v \in V \\
    &X_{vi} \leq A_{vi}  \quad \forall v \in V, \forall i \in V \\
    &X_{vi} = X_{iv}  \quad \forall v \in V, \forall i \in V \\
    &X_{vv} = 0  \quad \forall v \in V\\
    &X_{vi} =\{0,1\}, \quad \forall v \in V, \forall i \in V
\end{align*}

The $(W_v-W_v\circ X_v)$ is the difference vector between the original incident vector and the new incident vector for vertex $v$. Note that $W_v\circ X_v$ is the pairwise multiplication between these two vector also known as hadamard operation. The objective function is the summation of the norm of this difference vector over all the vertices in the graph. We wish to minimize this quantity subject to several constraints. The first constraint bounds the degree of the new graph by a user input positive bound defined for each vertex. The second constraint ensure that the new graph is a subgraph of the original graph that is no new edges are created. The third and fourth constraint is necessary to output a loop free symmetric graph.

Now what I need is the simplification of the objective function. Ideally I want to show equivalence of the objective function to something like this
$$
{\rm max} \sum_{v \in V} ||W_v \circ X_v||_2
$$

Any idea of how to simplify the objective function?


问:
Finite Element Model of Euler-Bernoulli Beam Theory with Isoparametric Element
In the formulation of Euler-Bernoulli Beam Theory, there are two degrees of freedom at a point, $w$ and $\frac{dw}{dx}$. Typically, the finite element model of this theory uses cubic polynomial for interpolation of $w$ using a two noded element as given in Chapter 5 of this book. This element is a subparametric because the geometry is represented by low-order elements than those used to approximate the dependent variable. I need a reference where I can find the FE model for this problem using isoparametric elements.


问:
Open source multiphysics software that can model the interaction of two insulators with a user-specified electrical force law?
Given two solid bodies of masses $m_1$ and $m_2$ with invariable charge distributions $\rho_1(\mathbf{r_1})$ and $\rho_2(\mathbf{r_2})$ in them,* respectively, I would like to model the electrical interaction of the bodies using a user-specified electrical force law.

* $\mathbf{r_1}$ and $\mathbf{r_2}$ are measured with respect to a coordinate system fixed to body 1 and body 2, respectively.

What open source multiphysics software packages could I use to model the motion of these bodies?


问:
PDE discretization on triangular domain
Given the 2D Poisson equation

$$\Delta u = f\\ u(x,0) = g_1(x), 0<x<1\\u(0,y) = g_2(y), 0<y<1\\
\partial_n u (x, 1-x) =0, 0<x<1$$

defined on the domain $\Omega := \{(x,y) \in \mathbb{R}^2_+: x+y<1\}$, discretize the equation using Shortley-Weller approximation with 5-point stencils.

On this domain, however, it appears that all near-boundary points use standard 5-point stencils. Is my understanding correct? I.e., with a uniform mesh non-boundary points' descritization uses the same step $h$ in all directions.

I'm, however, puzzled about the Neumann condition. Is this where Shortley-Weller approximation actually applied? I would appreciate some help with pointing me in the right direction towards discretizing the Neumann condition, and this equation as a whole.

Update: I figure the Neumann condition for this domain can be discretized as follows: since $\partial_n u = \nabla u \cdot n$, we have $\partial_{n} u_{i,j} = \frac{u_{i, j+1}-u_{i,j-1}+u_{i+1,j}-u_{i-1,j}}{h}=0$ (because $n=(1, 1)^T$). Now it looks like we encounter two ghost shells in each case. Am I on the right track?

Title: PDE discretization on triangular domain


问:
Minimizing the ratio of two specific non negative quadratic convex functions
$F$ is $m\times m$ diagonal, with real non negative elements

$D$ is $n \times m$ complex

$P$ is $n \times 1$ complex

$A$ is $m \times 1$ complex.

Minimize $\Gamma(A)$, with respect to $A$.

$$\Gamma(A) = \frac{m^2(DA-P)^H (DA-P) + (FA)^H(FA)}{A^HA}$$

It is known that both numerator and denominator of $\Gamma(A)$ are convex and non negative. Also both the terms of the numerator are individually convex and non negative.

Question:  An numerical optimization algorithm to find global minimum.

Apart from a reagular solution, I am also interested in a gradient descent based method if possible as  the matrices are large. Also $m >> n$.

PS : This question is a specific version of this question.

EDIT : more known information

No constraints on problem but

$\sum P = 0$, I mean sum of elements of matrix $P$ is zero.

Diagonal elements of $F$ are not all zeros.

Also $P^HP \ne 0$.

Rows of $D$ are orthogonal to each other. Also they are linearly independent.


问:
Determine set of dividend and divisor for quotient
I have the following Problem:

I know the range within the results of a division must lie.
[ quotientrange ] Additionally, the quotient should not exceed a certain number of fractions.

I know, that both dividend and divisor must consist of a certain set of digits. Dividend and divisor are further restricted by a range.
Edit: Dividend and divisor are only integers.

My question is: How can I determine a set of dividends and divisors, which satisfis this assumption?

Example:

Dividend should be within the rage [100.00; 999.99] and only consist of the digits 1, 2, 3, 4, 5

Divisor should be within the range [10; 99] and only consist of the digits 1, 2, 5.

Thus, the maximum range of the result is clearly determined by [100.00 / 99; 999.99 / 10]

Question: How can I determine a set [dividend1,dvisor1;....;dividend_n,divisor_n] which satisfies my assumptions?

Brute force is an alternative, certainly, by calculating from an assumed result (the range is known), constructing a valid divisor, performing a multiplication and chcecking, wheater the result satisfies the constraint for the dividend.

However, I wonder if there is a more intelligent approach?


问:
Cardinal B-Splines with derivative information
Have Schoenberg's cardinal B-splines been extended to accept derivative information at each knot, similar to how Lagrange interpolation can be improved by Hermite interpolation?</p>


问:
Structural Analysis Library
Can anyone recommend a structural analysis library that satisfies the following requirements:

C++ API

Simulate both beam elements and shell (slab) elements

Both static and dynamic analysis

Free and/or open-source

We're essentially looking to find a library that can do analysis similar to e.g. SAP2000.

However we are looking to incorporate this analysis in our larger software product, so we'd like the library to have an API that we can include and link against (i.e. has header(s) and dll(s)).

I am aware of, at least, libMesh. Are there other popular libraries I should be considering?

Title: Structural Analysis Library


问:
Cost functions to judge time/memory/accuracy tradeoffs
I am working on an interesting algorithm: Its absolute error is exponential in a parameter $j \in \mathbb{N}$, and for a given $j$, I have complete freedom to choose between an $\mathcal{O}(1)$ time-complexity algorithm with memory consumption ${\sim}2^{j}$, or a $\mathcal{O}(2^j)$ time-complexity algorithm with $\mathcal{O}(1)$ memory consumption. I can also "meet in the middle" and consume any intermediate amount of memory and compute that I wish.

I'm having trouble deciding on how make the compute/memory/accuracy tradeoffs. Have there been cost functions developed that could help answer this question?

Obviously, I could write something like
$$
\mathrm{cost} = c_1 t + c_2 m + c_3 e
$$
for evaluation time $t$, memory consumption $m$ and error $e$, and just pick three positive constants $c_1, c_2, c_3$. But I would prefer to use a cost function that someone has already put some thought into.


问:
Finite element method for Surface integrals using polar coordinates
I am trying to solve a 2D elliptic PDE (see complete electrode model for electrical impedance tomography) using the finite element method (FEM) over a circular region $\Omega$. I have discretized the model, and got the following surface integral

$$\int_e \phi_i(x,y) \phi_j(x,y) \mathrm{d}s, \mbox{ for } i,j = 1,...,n\, .$$

Here $n$ is the number of mesh points, $\phi_i$ are the basis functions and $e$ is an electrode placed on the boundary connecting two adjacent boundary points.

Can I use polar coordinates combined with the standard basis functions in FEM to evaluate these integrals? If so, how can I do this?


问:
Detecting blocks in non-linear system of equations
When solving systems of non-linear equations using Newton's method, it is often observed that the system has an independent sub-system, e.g. :

$$ f(x,y) = 0 $$
$$ g(x,y) = 0 $$
$$ h(x,y,z) = 0 $$

If I am not mistaken, you can (and should) solve first for $x$ and $y$ and then subsequently solve for $z$.
Now when the dimension of the system increases, it is not so easy to spot such independent sub-systems. Are there algorithms adapted to this issue ?

In my specific case, the Jacobian of the system has a fixed structure so even a slow algorithm for detection would be largely beneficial when solving the non-linear system using Newton's method. I feel that this is a common topic in solving systems of equations, so I am open to be redirected to additional ressources that are relevant.


问:
Richardson's Iteration, Gradient Method and Spectral Radius
Richardson's iteration introduce a scalar $\alpha$ to the update formula:

$$ \textbf{x}^{(k+1)} = \textbf{x}^{(k)} + \alpha \textbf{r}^{(k)} $$

And compute $\alpha$ by minimizing the spectral radius:

$$ \min_{\omega}\,\rho(B) = \min_{\omega}\, \rho(I-\omega A)\, . $$

As $I$, and $A$ do not change over the iterations, it might seem it exists only one (i.e. equal for all the iterations) $\omega$ making the max eigenvalue minimal, hence the spectral radius minimal. Since I know gradient method and conjugate gradient perform better by setting $\alpha$ dynamically, I am wondering what am I missing here. Is there any way to see through the spectral radius expression that gradient and conjugate gradient iterations converge faster?


问:
Does adaptive Gauss-Kronrod reuse function evaluations?
I'm curious to know how QUADPACK's QAG routine works.  My understanding is that it begins by calculating on each subinterval the numerical quadrature with a Gaussian-Legendre rule and a nested Kronrod rule to provide estimates of the integral and the local error.  Then, globally the interval with the highest error is bisected and the process is repeated. However, since Gauss-Legendre abscissae do not nest, can any of the previous function evaluations be reused at the next iteration? Or are all previous function evaluations discarded when an interval is bisected and reevaluated? If the latter is true this seems incredibly inefficient, so why would anyone choose such an adaptive quadrature scheme like that in QUADPACK over one where the function evaluations may be reused, such as Newton-Cotes or Clenshaw-Curtis?</p>


问:
Grid Data Interpolation
What are the most sophisticated methods for interpolating a scalar field say Electric or Magnetic Field on a 3-D grid?

I have scalar data on a meshgrid with equal spacing. I would like to use an interpolating method to interpolate those values onto random points I select on the same grid. I think my best bet is using tricubic splines. Am I right? Some of the scalars might not be linear with others being linear though. Any guidance would be appreciated.

Title: Grid Data Interpolation


问:
MATLAB: Compute the Schwarz-Christoffel transformation symbolically
Suppose we have a conformal mapping from the unit disk in the $\omega$ plane onto the exterior of a polygon in the $z$ plane.

The  Schwarz-Christoffel mapping in this case is defined as:
$$f(u) = A - C \int^u w^{-2}\prod_{k=1}^{n} (1-\frac{\omega}{u_k})^{1-\alpha_k} d\omega.$$

Meaning that we determine the vertex $z_k$ with the mapping $f(u_k)$ from the respective prevertex.

I have provided an image from Discroll and Trefethens' book on the subject (PDF sample with the image is freely available).

I would like to compute a factor from an integral over function analytic inside a polygon in the complex plane.

For this I will use Gauss-Kronrod quadrature, specifically MATLAB's quadgk command with 'waypoints'

I ll just pick my function arbitrarily as $g = \exp(z)$.

What I want to compute is:

$$\int{\frac{\exp(f(u))}{u}}du$$

To use quadgk I have to define a function handle for $g$.
But I am confused on how to do that with the SC-Mapping f(u).
I know I can extract the prevertices and interior angles $\alpha_k$ with the SC-toolbox for any polygon like so:

```<code>PSI = extermap(drawpoly)
p = polygon(PSI);
a = parameters(PSI);
</code>```

How do I compute the integral $\int^u w^{-2}\prod_{k=1}^{n} (1-\frac{\omega}{u_k})^{1-\alpha_k} d\omega$ ? It's a polygonal integral in the complex plane so should I use quadgk as well ?

How do I set the complex constants $A$ and $C$ ?

b = a.constant

seems to refer to the logarithmic capacity of the polygonal domain.


问:
Constraining the total volume in Finite Element Methods
I have a diffusion problem which can be broken down to be:

$-\Delta u = f(u) $    on $\Omega ~/~ \Omega_{int}$

$u = 1$  on $\Omega_{int}$

Note that this is an internal Dirichlet constraint to the value 1.0 in the center circle, unlike the standard textbook examples with homogeneous boundary conditions. (see picture) Think of it like a source in the middle, which diffuses into the domain with f(u) = u(1-u) being a logistic growthterm. I want that reaction to stop at a certain maximum total mass.

I want to impose an integral condition which will stabilize the solution:

$P(u) = \int_{\Omega}  u  ~ dx =V_{tot}$

Now, I tried to read into the problem, and I find that it is possible to add penalty terms to the bilinear form of the finite element discretization, but I am deeply puzzled as to how to properly define the penalty term so that my Newton solver will behave correctly.

Multiplying with a test function v, integrating over the domain we may formulate the problem without the volume constraint like this:

find u, v so that

$ r(u,v) = a(u,v)+ l(v) = 0 $

where $r(u,v)$ is a global residuum (To be minimized by the Newton solver), $a(u,v)$ is the canonical discretisation of the diffusive term and $l(v)$ the disc. of the reaction term.

Now in the language of Lagrange multipliers, we have a functional $r(u,v)$ and a constraint $P(u)$ on the volume so that the new Functional would look like:

$L(u,v,\lambda)= a(u,v)+ l(v) -\lambda(P(u)-V_{tot}) = 0 $

How do I apply this to my Finite Element Method in a clean way?

How would I determine $\lambda$?

Is there any readable literature on how to construct volume-constraints via penalty terms?

Is there a more readable/hands-on literature/Paper on Lagrange multipliers in the context of FE-Methods than the work of Ivo Babuska?

How do I make sure that the penalty term, which will probably contain a domain-wide integral will not mess up my Newtonian (since the necessary calculations for the Jacobian are no longer local)?


问:
Existence and uniquness of solution of FVM for Poisson equation
I'm discretizing the following Poisson equation using FVM where the domain $\Omega$ of the solution is a regular hexagon of side $1$ centered about the origin.
$$\Delta u =k,\text{ $k$ constant}\\
\partial_n u = 0\\
\Omega \subset \mathbb{R^2}$$

I'm not really quite sure if my discretization and especially the argument for existence and uniqueness are correct. So I would appreciate some feedback.

Let $V_i\subset \Omega$, $i=1,\dots,6$, be one of the six equilateral triangles comprising the hexagon; so $\bigcup\limits_{i=1}^6 V_i=\Omega$.

Then we obtain a system of six equations, that is a linear system after the following steps:
$$\int_{V_i} \Delta udx = \int_{\partial V_i} \nabla u(x)\cdot n ds =\int_{V_i}kdx$$
$$\int_{partial V_i} \nabla u(x)\cdot nds=\int_{\partial V_i}\partial_n u ds$$
$$\approx \sum\limits_{j\in n_i}\frac{u_j-u_i}{h_{i,j}}l_{i,j}+\int_{\Gamma_i}0ds=\sum\limits_{j\in n_i} \frac{u_j-u_i}{h_{i,j}}l_{i,j}+\text{constant $C$}$$
$$=|V_i|k$$

Where $|V_i|$ means the area of $V_i$ and $n_i$ is the set of cells immediately neighboring $V_i$ (two such cells for a hexagon).

Now, we can use some geometry to deduce that

$$|V_i| = \frac{\sqrt{3}}{4}; h_{i,j} = \frac{\sqrt{3}}{2}; l_{i,j}=1$$

So the system becomes:

$$\frac{2}{\sqrt3}\sum_{j\in n_i} (u_j-u_i) + C = \frac{\sqrt3}{4}k$$

Which can be represented in matrix-vector form as

$$B\vec{u}=\vec{c}$$

where $B$ is a square $6\times 6$ matrix with the following components

$$b_{ij}=\frac{2}{\sqrt3}, i\ne j; b_{ii}=-\frac{2\sqrt2}{3}$$
$$c_i = \frac{\sqrt3}{4}k-C$$

Now, the "most interesting" part: is it enough to argue that there is a unique solution and it exists because $B$ is non-singular?


问:
Paring Large Numbers Question
I am writing some simple code to raise a base value to a power
then perform some operations on the output.  I have included the code I am
using now, below ( see end of this message ), also see immediately below
what this looks like when I run my script from the console ;

```<code>C:\PYTHON34>python tst.py
INPUT Value: 431.1
INPUT Power Value; 1.9907
RESULT 1739.554505641426658257842063903808593750000000000000000000000000
</code>```

What I need to be able to do is:

Store whole number value to left of the decimal point without decimals to a variable = x

Grab entire decimal value to the right of the decimal point and store in a variable = y

Parse (y) according to some simple rules. Here it is where it gets trick for me.

What I want to do is examine (y) to see if there are any leading leading and trailing zeros

In the example 1739.554505....., if this value was instead, something like any of the following

```<code>1739.0554505
1739.00554505
1739.000554505
</code>```

Then, I want to extract the whole number less any continuous leading zeroes.

Next, I want to cut the number so that when there are four continuous trailing zeros, all zeros starting from the first one that started off the first four continuous zeros and all following zeroes are truncated.

So in our example

```<code>1739.554505641426658257842063903808593750000000000000000000000000
</code>```

(y) becomes 55450564142665825784206390380859375

Next, I want to take the mod of above value by %1999. This returns 1407.

Next, I want to join this value AS the decimal value to the right as

```<code>1739.1407
</code>```

Ok, so now why am I stuck at this point?

I'm still fairly new to learning Python so I don't actually know what the best ways are to do this.

I'm thinking I could maybe just do a trim on the value to the right of the decimal, take that, convert the whole thing to a string, then chop up each digit and store each individually to an array then code a rule set for stepping through strings in array cells, but this also seems like it could take more CPU  time than I'd like.  I'm hoping to perform this sort of operation fast, as fast as possible.

Also, I do not even know how to best approach this even if I had to store this as a string into an array, and if I understand correctly string processing tends to takes more CPU time than number calculations.

Any insight, direction, code suggestions greatly appreciated!

This is the code I am currently using

```<code>import time
e0=time.time()
c0=time.clock()

import sys
import math
from sympy import mpmath
from mpmath import *
mpmath.mp.dps = 10000

inputx = float(input('Enter inputx Value: '))
powerx =float(input('Enter Power Value; '))

inputx_squared = float((mpmath.power(inputx, powerx)))%1999

print('\r')
print('RESULT ', '%.60f' % inputx_squared)

elapsed_time=time.time() -e0
cpu_time=time.clock() -c0

print('\r')
print("CPU TIME ;", cpu_time)
</code>```

Title: Paring Large Numbers Question


问:
How can I maximise orthonormality between degenerate eigenvectors using ARPACK?
I am using ARPACK's zndrv1 to diagonalise a matrix (the context is quantum chemistry). While all vectors have a norm 1, as expected, vectors corresponding to degenerate eigenstates aren't always orthogonal (often their inner product will be roughly 0.01).

Are there any parameters I can adjust to improve orthonormality (even if it means longer compute time)?

I understand there are alternative libraries and routines out there. I am specifically curious about ARPACK.


问:
Determine truncation error of PDE discretization
The equation is $$\frac{\partial}{\partial x}\left(u\frac{\partial u}{\partial x}\right)=f(x)\\ 0<x<1, u(0)=u(1)=0$$

I'm discretizing this PDE using FVM as follows:
$0=x_0=x_{1/2}<x_1<x_{3/2}<\dots<x_{N+1/2}=x_{N+1}=1$

$h_i:=x_{i+1/2}-x_{i-1/2}$

$h_{i+1/2}:=x_{i+1}-x_i$

$$\int\limits_{x_{i-1/2}}^{x_{i+1/2}} \frac{\partial}{\partial x}\left(u\frac{\partial u}{\partial x}\right)dx=(uu_x)(x_{i+1/2})-(uu_x)(x_{i-1/2})=\int\limits_{x_{i-1/2}}^{x_{i+1/2}}f(x)dx$$

Using finite differences, obtain

$$u_{i+1/2}\frac{u_{i+1}-u_i}{h_{i+1/2}}-u_{i-1/2}\frac{u_{i}-u_{i-1}}{h_{i-1/2}}=h_if_i$$
where $u_0=u_{N+1}=0$, $f_i:=\frac1h_i\int\limits_{x_{i-1/2}}^{x_{i+1/2}}f(x)dx$

So now comes the part where I got really stuck. I want to determine the truncation error of this discretization. For this I substitute the Taylor expansion of the exact solution into the equation:

$$r_i:=\frac1h_i\left(u(x_{i+1/2})\frac{u(x_{i+1})-u(x_i)}{h_{i+1/2}}-u(x_{i-1/2})\frac{u(x_{i})-u(x_{i-1})}{h_{i-1/2}}\right)-f_i$$

$$=\frac{1}{h_i}\left(u(x_{i+1/2})\left(u_x(x_i)+\frac{h_{i+1/2}}{2}u_{xx}(x_i)+O(h^2_{i+1/2})\right)-u(x_{i-1/2})\left(u_x(x_i)+\frac{h_{i-1/2}}{2}u_{xx}(x_i)+O(h^2_{i-1/2})\right)\right)-f_i$$

I'm not sure what to do with the terms $u(x_{i+1/2})$ and $u(x_{i-1/2})$ to arrive at the truncation error, which I think will be $O(h^2)$. I'd appreciate some help with this.

I've already tried to expand these terms about $x_{i-1/2}$ but that wasn't quite useful.


问:
Discretizing a parabolic PDE with finite volume method
I want to discretize the following parabolic PDE:

$$u_t = \nabla\cdot(\alpha(x)\nabla u)- \beta u\\
x\in\Omega \subset \mathbb{R}^2\\
\partial_n u = 0\\
u(t,0) = u_0(x)\ge 0, \alpha(x)>0$$

Given the Neumann boundary condition above, one can integrate this PDE:

\begin{alignat}{2}
\frac{\mathrm d}{\mathrm dt}\int_\Omega \! u\,\mathrm dx&=\phantom{-}\int_\Omega \! \nabla\cdot (a(x)\nabla u)\,\mathrm dx \:\:&-\int_\Omega \! \beta(x) u \,\mathrm dx \\
&=\phantom{-}\int_{\partial\Omega}\! a(x)\nabla u\cdot n \,\mathrm ds&-\int_\Omega \! \beta(x) u \, \mathrm dx \\
&=\phantom{-}\int_{\partial\Omega} a(x)\partial_n u \,\mathrm ds&-\int_\Omega \! \beta(x) u \, \mathrm dx \\
&=-\int_\Omega \! \beta(x) u \, \mathrm dx
\end{alignat}

The form of the discretized linear system I want to get is

$$\frac{d\vec{u}}{dt}=A\vec{u}-\vec{c}$$

Unfortunately, I don't think this is exactly what I'm obtaining, despite I do think I use a correct discretization technique. To this end, please consider my steps below.

I discretize the two integrals using the finite volume method:

$$\frac{d}{dt}\sum\limits_{j\in n_i} u_j l_{ij}=-\sum\limits_{j\in n_i} \beta_{ij} u_j l_{ij},$$

where $n_i$ is the set of neighboring cells to the cell $V_j\subset \Omega$ and where $l_{ij}$ is the length of the boundary of $V_j$ between $u_i$ and $u_j$.

Now one can discretize it in time to obtain:

$$\sum\limits_{j\in n_i} \frac{u_j^{n+1}-u_j^n}{\Delta t}l_{ij}=-\sum\limits_{j\in n_i} \beta_{ij} u_j^n l_{ij}$$

So we obtain the linear system
$$L\frac{d\vec{u}}{dt} = -A\vec{u},$$

where $L$ is the matrix containing $l_{ij}$ and each entry of $A$ contains $-\beta_{ij}l_{ij}$.

Am I correct? I think I'm quite confused as I don't see where each term of the matrix should go, and where the vector $\vec{c}$ is supposed to be. Moreover, what does one do with the $l_{ij}$ on the LHS?


问:
Translating grid with extrusion speed
I am putting into MATLAB code the equations that describe a plastic extrusion process. From a paper, I found I should use a spatial grid that translates with the extrusion speed, being the reference coordinate system the centerline of the extrudate.

In particular, $s = 0$ is the extruder die and the length of the domain is time dependent as the extrudate gets longer and longer. In the paper, new grid points are continually added to the solution domain at $s = 0$ (the extruder die) to account for the extrusion of fresh fluid. How would I do this?

Title: Translating grid with extrusion speed


问:
3D Tollmien-Schlichting Waves Imposed in a Channel Flow (Are Physics correct?, etc)
So I am trying to do some further tests on a 2nd-order code Incompressible Navier Stokes equations, by studying transition to turbulence in a Poiseuille flow. Specifically, I'm interested to see whether the code is energy-conserving by modelling a 3D TS wave in a Poiseuille flow. Note that I have tested the code for a 2d inviscid Taylor vortex and the code holds. Now, I have calculated the perturbation velocities via the Orr-Sommerfeld/Squire equations, and I get the Eigenspectrum and corresponding Eigenfunctions (velocities/wall-normal vorticity) for the following Poiseuille flow:

$U^0_1 = 1-y^2,\,\, -1 \le y \le 1$

With a Reynolds number based on the centreline velocity and the channel's half-height

$Re_{cl}=\frac{U_{max}(h/2)}{\nu}=5000$

From here, the driving pressure gradient is just:

$\Pi_1 = -\frac{2}{Re}$

For the linear perturbation I choose a streamwise and spanwise wavenumbers of $\alpha = 1.12, \beta= 2.1$, respectively. For the eigenvectors, the most unstable mode has been chosen (that where Real(Lambda) is the closest to zero). All the results are the following:

<img src="https://i.stack.imgur.com/kgtGT.jpg" alt="Eigenvectors for $\lambda \approx -0.07 - 0.36i$">

I have imposed the velocity/vorticity perturbations into my code, without success (meaning not getting the transition), so clearly there's something I'm not quite getting right. Now, my questions are:

Orr-Sommerfeld Modelling: Are my results obtained (meaning perturbation velocities/vorticity) reasonable? Are these correct?

Equations to solve for temporal transition modelling: In my opinion, the equations to be solved when studying spatio-temporal transitions are the incompressible Navier-Stokes equations (INSE) for the perturbation velocities, that is:

$ u_{i,t} +u_j u_{i,j} + u_j U^0_{i,j} + U^0_j u_{i,j} = -p_{,i}+\frac{1}{Re}u_{i,kk}+\Pi_i$

Which results from the splitting of the instantaneous fields into base  and perturbation fields, $\widetilde{\phi_i} = \phi^0_i + \phi_i$ in the INS, and then substracting the INSE of the base flow. Now my questions are: Is the previous equation correct? I'm having doubts about the pressure gradient $\Pi_i$. Second: Should be these the equations I should be solving for simulating temporal transitions? Or, Should I go just with the classic INSE for the poiseuille flow and superposed fluctuations?


问:
Biconvex problem whose objective function depends on only one variable
I am solving the following biconvex problem:
$$\min_{x,y} f(y)$$
$$s.t. ~~ g(x) \leq 0$$
$$~~~~~h(x,y) = 0$$
$$x \in X, y \in Y$$
where $X$ and $Y$ are compact convex sets, $g(x)$ and $f(y)$ are convex and $h(x,y)$ is affine for fixed $x$ or $y$.
I wonder if there is any efficient algorithm to solve this type of optimization problem, not restricted to algorithms for general biconvex optimization problems.


问:
Connection between piecewise linear basis functions and RELU activation function
ReLU activation is defined as follows
$$\sigma(x)=\max(0, x).$$
Let's assume that I have deep network of 1 hidden layer, than output from my layer has form
$$ f(x)= \sigma(Wx +b), $$
where matrix W represents the weights, b is bias and x stands for input data (output of previous layer).

What is connection  between $f(x)$ and piecewise linear function, defined as
$$g(x_j)=\sum_{i=N}^N \alpha_i \psi_i(x_j),$$
where $\alpha_i$ stands for the coefficients and $\psi_i$ for the basis function.

More precisely, I would like to identify basis functions and coefficients in the definition of $f(x)$.


问:
Can I use the Schur basis returned by ARPACK in a restart capacity?
Reading ARPACK documentation, I see that ARPACK will return an "orthogonal basis for the invariant subspace corresponding to the eigenvalues in D" if eigenvectors are not requested. Can this subspace be exploited in some kind of restart capacity? Documentation discusses restarts but I'm not sure how these have been implemented in the routines.

I know you can set the residue vector to an initial guess, but I'd like to know if a pre-generated subspace can be used.


问:
Numerical solution to N-dimensional diffusion on simplex?
Assume I have a system of at least (but generally only) $N+1$ points in an $N$-dimensional space ($N > 3$ is possible). At each of these points $x_i, i=1,...,N+1$ I know an initial potential/concentration/quantity $P_i$. I can evaluate the possible connections between all points. I also have an isotropic diffusion coefficient $D$. I am now interested in calculating diffusion between these three points (with mass conservation). I have sketched an example case in 2-D.

This problem sounds rather easy, but I quickly realized that it isn't. Classic finite difference approaches require (to my knowledge) a structured gid, which I don't have. Finite volume could handle unstructured grids, but I cannot afford tessellation in high dimensions (say, $N > 100$). For both approaches I would have calculated the first potential derivative along each line (e.g., $x_1$-->$x_2$), and assumed that the same derivative half a distance outside (dashed lines) is zero. Unfortunately, I don't think this assumption works if I have more than $N+1$ points. Finite elements might be an option, but I have yet to find a tutorial which would cover such a case.

Do you have any idea how to solve this system?

Edit: Some additional clarification concerning my goal: I want to simulate the equilibration between the $N+1$ vertices. I am only interested in solutions at the vertices themselves, and diffusion would occur along the edges. The specific equation I want to solve is the heat equation:

$$\frac{dP_i}{dt}=K\Delta P_i$$


问:
What will PDE discretization matrix look like for time and space?
Please note: this question is not a duplicate of this question since, while the PDE is the same, the nature of this question is different, i.e. the other question treats a different aspect of this PDE.

Here is the PDE system I am working with:

$$\partial_tu = \nabla \cdot (a(x)\nabla u)-\beta(x)u\\
\partial_nu=0, x \in \Omega \subset \mathbb{R}^2\\
\beta(x)>0$$

Integrate the PDE:

$$\int_\Omega \partial_t u=\int_\Omega \nabla\cdot (a(x)\nabla u)dx-\int_\Omega \beta(x)udx$$

Apply the Divergence theorem:

$$\int_\Omega \partial_t u=\int_{\partial\Omega} a(x)\nabla u\cdot nds-\int_\Omega \beta(x)udx$$

Apply the Neuman boundary condition:

$$\int_\Omega \partial_t u=\int_{\partial\Omega} a(x)\partial_n uds-\int_\Omega \beta(x)udx=-\int_\Omega \beta(x)udx$$

Obtain a new PDE:

$$\partial_t u+\beta(x)u=0$$

This PDE is in essence an ODE, so we can get the following solution for this ODE:

$$u(x,t)=e^{-\beta(x)t}D(x),$$
where $D(x)$ is an integration "constant".

Discretize the ODE:

$$\frac{u_{ij}^{n+1}-u_{ij}^n}{\Delta t}=-\beta_{ij} u_{ij}^n$$
or
$$u^{n+1}_{ij} = u^n_{ij}(1-\Delta t \beta_{ij})$$

Now we can assign lexicographical ordering to each $i,j$ and denote it by a function $p(i,j)$, thus rewriting the discretized PDE about as follows:

$$u^{n+1}_{p(i,j)} = u^n_{p(i,j)}(1-\Delta t \beta_{p(i,j)})$$

Am I correct to deduce that for each $n$ we have a linear system of the form

$$U^{n+1}=AU^n$$

where $A$ is constant and $A = diag(1-\beta_{ij},\dots , 1-\beta_{ij})$? If this is true then is it also correct that there are thus as many such linear systems as there are $n$'s (say, $M$) and to solve this discretized system we need to solve $M$ such linear systems?

I think, alternatively, we can probably extend the system to include not only $i,j$ in the lexicographical ordering, but also to include $n$ in it, so as to get only one linear system to solve.


问:
Analytic vs discrete understanding of PDE
The PDE I am working with:

$$\partial_tu = \nabla \cdot (a(x)\nabla u)-\beta(x)u\\
\partial_nu=0, x \in \Omega \subset \mathbb{R}^2\\
\beta(x)>0$$

Integrate the PDE:

$$\int_\Omega \partial_t u=\int_\Omega \nabla\cdot (a(x)\nabla u)dx-\int_\Omega \beta(x)udx$$

Apply the Divergence theorem:

$$\int_\Omega \partial_t u=\int_{\partial\Omega} a(x)\nabla u\cdot nds-\int_\Omega \beta(x)udx$$

Apply the Neuman boundary condition:

$$\int_\Omega \partial_t u=\int_{\partial\Omega} a(x)\partial_n uds-\int_\Omega \beta(x)udx=-\int_\Omega \beta(x)udx$$

Obtain a new PDE:

$$\partial_t u+\beta(x)u=0$$

This PDE is in essence an ODE, so we can get the following solution for this ODE:

$$u(x,t)=e^{-\beta(x)t}D(x),$$
where $D(x)$ is an integration "constant".

Discretize the ODE. After discretization I should obtain a linear system of the following form:

$$\frac{du}{dt}=Au - c$$

But when I discretize this system I obtain the following:
$$\frac{du_{ij}}{dt}=-\beta_{ij} u_{ij}$$

So one could say that $A$ in this case has the form diag$(-\beta_{ij}, \dots, -\beta_{ij})$, but then where is the vector $c$ coming from?

Also, one could observe that, since $\beta(x)>0$,
$$\lim\limits_{t\to \infty} u(x,t) = 0$$.
But how can we deduce the same from the discretized system?


问:
Oscillations when solving parabolic heat equation with FTCS
I'm wondering if someone could help me out, or point me in a direction of how I can understand the following oscillations that occur when I solve the Porous Medium Equation $$u_t = u_{xx}^{m+1}$$ using FTCS method. The following plots occur for a the following initial condition
$$
f(x) =
\begin{cases}
1 & \text{if}\ |x| \geq 1 \\
1 - e^{-1/(1-|x|^2)} & \text{if}\ |x| < 1
\end{cases}
$$
with Dirichlet boundary conditions as $u(x_0, t)= u(x_M,t) = 1$.

By adjusting the step size for space and time, I end up with the given plots below. The surprising result for me, who just learned about finite difference methods, is that the numerical solution does not converge for values where $$\frac{k\ \text{max}\{u^{m}\}}{h^2} < \frac{1}{2}.$$
How come? It seems to converge for values $< \frac{1}{4}$ instead.


问:
How to numerically calculate the transition dipole integral in periodic systems?
Now I have wave functions $\psi_a$ and $\psi_b$ of two states in Gaussian CUBE format. I'd like to evaluate the transition dipole moment integral $\pmb\mu$ between these two states. As my simulation cell is cubic periodic with length of $L$, I will need the Berry phase formula to calculate $\pmb\mu$:
$$\pmb\mu= \rm{IMAG}\; \rm{LOG}\; \langle \psi_a|\exp(\it i\frac{2\pi}{\it L}\pmb r)|\psi_b\rangle\; $$
My implementation is as follows: $n$ is the number of grid points of CUBE file in each direction. $b$ is the step length. So the lattice constant $L$ is $nb$. $\pmb r$ is the coordinate. psi_a and psi_b are the volumetric data which are changed to 1D arrays. $dd=\exp(i2\pi \pmb r/L)$.

However, the above implementation didn't give me correct numbers as I expected. So please can anyone help me point out the problems of the implementation?

OwnerDisplayName: xmW


问:
High-accuracy numerical differentiation
I have a $200 \times 200$ matrix representing the values taken by a function over an equally spaced grid. I would like to perform derivatives on it.

I am interested in its gradient (i.e. its derivative in direction $x$ and its derivative in direction $y$) and in its Laplacian.

I work on the Matlab platform, where I use the built-in functions gradient and del2 which work very well but their accuracy is limited by the fact that they make use of a small number of points.

Reading this Wikipedia page about Finite difference coefficients, one can understand that there are finite-difference schemes to perform numerical differentiation in a more accurate way. The price to pay to have an increased accuracy is of course to use more complex formulas, which include more points.

I would like to know if there is a library where there are functions capable of doing this job which -I repeat- is: computing the first-order and the second-order derivatives of a 2D matrix (representing the values of a not-explicitly-known function) with a user-defined accuracy.

It would be great if someone could suggest me Matlab libraries, but also C/C++ libraries could work, I guess.

Title: High-accuracy numerical differentiation


问:
Use of non-typical values of $\theta$ in theta-methods
The theta-method is a popular solution for solving time-transient PDEs (or ODEs), which consists of solving the general equation for each time step:

$$
\frac{u^{n+1} - u^{n}}{\Delta t} + (\theta f(u^{n+1}, t^{n+1}) + (1-\theta)f(u^{n}, t^{n})) = 0
$$

Some values of $\theta$ are often use and their pros and cons are well documented, notably:

$\theta = 0$ is purely explicit, and as such has a low computational cost

$\theta = 0.5$ (AKA the Crank-Nicolson method) has a precision in $O(\Delta t^2)$, vs $O(\Delta t)$ for other values

$\theta = 1$ is unconditionally stable for oscillations

I have also sometimes encountered the value $\theta = \frac{2}{3}$, referred as the Galerkin method, for example in [1] (p.1814), without any mention of its advantages over other values for $\theta$. The fact that the Galerkin method also refers to the much more popular method of transforming a PDE problem through a weak formulation doesn't help getting information.

Also, on a specific problem I'm working on, I find that the best results are obtained with a value of $\theta = 0.9$.

Do you know of any resources or have more info about the use non-typical values for $\theta$ in a theta-scheme? What are the possible advantages and caveats of such values?

[1] Dalhuijsen, A. J. and Segal, A. (1986), Comparison of finite element techniques for solidification problems. Int. J. Numer. Meth. Engng., 23: 1807-1829. doi:10.1002/nme.1620231003


问:
Instability in Lattice Boltzmann Solver
I wrote a Lattice Boltzmann Solver in Rust a little while back using both the BGK approximation and the TRT (two relaxation time) method on the D2Q9 lattice. In both cases I run into major stability problems even in simple examples. Here is a video of the output from the simulation of a left-to-right flow with a cylindrical obstruction (using the bounce-back method): link. This video shows momentum density using lightness as magnitude and color as direction. I can't remember if it was generated with BGK or TRT but I experienced similar failures with both collision operators. I'm using convolution with a [[1, 1, 1], [1, 0, 1], [1, 1, 1]] stencil to implement the streaming step. I tried following the steps described in The Lattice Boltzmann Method: Principles and Practice for getting stability in BGK and TRT simulations, but haven't had any luck. Am I doing something obviously wrong?

Title: Instability in Lattice Boltzmann Solver


问:
Need help applying Implicit Eulers Method together with Newtons Method on Burgers' Equation
From the inviscid Burgers' equation: $\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x} = 0$, I get the discretization $\frac{u_i^{n+1}-u_i^n}{\Delta t}+\frac{(u_i^{n+1})^2-(u_{i-1}^n)^2}{\Delta x}=0$. Having the starting conditions $u_0(x)=sin(\pi x)$, $x\in[0,1]$, I wish to plot $u$ after $n$ timesteps. However, I cannot seem to get it to work and would appreciate any help. Python code:

```<code>import numpy as np
import matplotlib.pyplot as plt
import scipy as py
from scipy.sparse.linalg import gmres

#x
x0 = 0 #x min
xe = 1 #x max
N = 100 #Nodes
dx = (xe-x0)/N
x = py.linspace(x0,xe,N)

#t
t0 = 0 #time
te = 1 #max time
dt = 0.0001 #timestep
M = int((te-t0)/dt) #timesteps

#u
u0 = py.sin(py.pi*x)

ts = 50

def EulerMethod(ustart):
    U = np.zeros((N,M))
    U[:,0] = ustart
    t=t0
    for i in range(N-1): #one step Euler forward
        U[0,1] = 0
        U[i+1,1] = U[i+1,0] - (dt/dx)*(U[i+1,0]**2-U[i,0]**2)
    for n in range(ts):
         U[:,n+2] = NewtonMethod(U[:,n+1],U[:,n]) #Eulers implicit method with Newtons Method
         t=t+dt
    return U

def NewtonMethod(u,c):
    #c = u_n
    for i in range(3):
        F = np.zeros((N,1))
        J = np.zeros((N,N))
        J[N-1,N-1] = 1/dt + u[N-1]/dx
        F[0] = (u[0] - c[0])/dt + (u[0]**2)/2*dx
        for j in range(N-1):
            J[j,j] = 1/dt + u[j]/dx
            J[j+1,j] = 1/dt + u[j]/dx
            F[j+1] = (u[j+1] - c[j+1])/dt + (u[j+1]**2 - u[j]**2)/2*dx
        du = gmres(J,-F)[0]
        u = u + du
    return u

UM = EulerMethod(u0)
plt.plot(x,UM[:,ts])
</code>```


问:
Properly sizing ParaView plots for LaTeX
I would like to make a visualization with ParaView for a paper. In order to have a high quality look, I would like the image to be, say, 300 dpi (or vector graphics). Moreover, since I know the column width of the paper, I'd like to somehow specify that width in ParaView so I can adjust the size of the visualization, colorbars, etc., to make it fit within the column. I'd also like to specify the text size and font to be the same as in LaTeX. In other words, I do not want to rescale the image afterwards, because that makes the text too small or too large. In other programs I usually figure out how to do this, but so far no luck in ParaView. Even if I set the font to, say, 10 pt, it turns out way too small when I include the figure in LaTeX with scale=1. Any help appreciated.


问:
extrapolation/interpolation in fmincg.m
Can you tell me these equations come from where in MATLAB fmincg.m?

```<code>A = 6*(f2-f3)/z3+3*(d2+d3);  % cubic fit
B = 3*(f3-f2)-z3*(d3+2*d2);
z2 = (sqrt(B*B-A*d2*z3*z3)-B)/A;
</code>```

It is either cubic extrapolation or interpolation.
I have looked for these equations in many websites and books of numerical optimization but but I do not find these equations even in the book "Numerical Optimization by Peter Glynn, Stephen M. Robinson p:57". It is not same in the fmincg.m

Title: extrapolation/interpolation in fmincg.m


问:
Finding second excited state of Schrödinger equation with secant Runge Kutta method
In our assignment, we are required to find the energies of the ground state and the first two excited states of the Schrödinger equation in a harmonic potential:

$$V = \frac{50 x^2}{(10^{-11})^2}\, .$$

I have the Python code for solving Schrödinger equation using secant Runge Kutta method:

```<code>import numpy as np

m = 9.1094e-31
hbar = 1.0546e-34
e = 1.6022e-19
L = 5.2918e-11
N = 1000
h = L/N

def V1(x):
    return 0.0
    #return (50*x**2)/(10**-11)**2

def V2(x):
    return 0.0
    #return (50*x**4)/(10**-11)**4

def f(r, x, E, V):
    psi = r[0]
    phi = r[1]
    fpsi = phi
    fphi = (2*m/hbar**2)*(V(x)-E)*psi
    return np.array([fpsi,fphi], float)

def solve(E, V):
    psi = 0.0
    phi = 1.0
    r = np.array([psi, phi], float)

    for x in np.arange(0,L,h):
        k1 = h*f(r, x, E, V)
        k2 = h*f(r+0.5*k1, x+0.5*h, E, V)
        k3 = h*f(r+0.5*k2, x+0.5*h, E, V)
        k4 = h*f(r+k3, x+h, E, V)
        r += (k1 + 2*k2 + 2*k3 + k4)/6
    return r[0]

print("Start program")

E1 = 0.0
E2 = e
psi2 = solve(E1, V1)

target = e/1000
while abs(E1-E2)>target:
    psi1, psi2 = psi2, solve(E2, V1)
    E1, E2 = E2, E2-psi2*(E2-E1)/(psi2-psi1)

print("For V(x) = V0x^2/a^2:")
print("E=", E2/e, "eV")
</code>```

However, there are two problems that arise that causes me not being able to find the solution to this assignment. And I hope that someone can help me with it.

First of all, when the function V1(x) returns 0.0, the program accurately outputs the energy of ground state for V1(x) = 0, which is E= 134.28637169369105 eV. However, when we change to V1(x) = V = (50x^2)/(10^-11)^2, which is what our assignment requires us to find, the program creates the following error:

```<code>C:/Users/wormw/Desktop/PHYS 3142 - Computational Methods in Physics/Assignment 5/Assignment 5.py:29: RuntimeWarning: overflow encountered in double_scalars
  return np.array([fpsi,fphi], float)
C:/Users/wormw/Desktop/PHYS 3142 - Computational Methods in Physics/Assignment 5/Assignment 5.py:54: RuntimeWarning: invalid value encountered in double_scalars
</code>```

And the output of E is E = nan eV. So, I have no idea why this would happen and how to fix it.

Secondly, the assignment also requires us to find the first two excited state, but I have no idea how to modify the code to obtain the eigenvalue for getting the excited states

Thirdly, the assignment also states that Try using x = -10a to +10a, with the wave function psi = 0 at both boundaries. Should I change the value of L at the start?


问:
Solve ODE with non-negative and maximization constraints
My task is to solve $$\eta_k\frac{d^2C_k}{dz}(z)=-e_k, k = 1,2,3$$

$$C_k\ge0$$

$$C_1(0)=0, C_2(0)=A, C_3(0)=0$$
$$C_1(L)=B, \frac{dC_2}{dz}(L)=0, \frac{dC_3}{dz}(L)=0$$
with $$e_1 = -\beta_1-\beta_3$$
$$e_2=-\beta_1-\beta_2$$
$$e_3=\beta_2-\beta_3$$

satisfying the constraints $$-e_k\le \frac{\eta_kC_k}{K_k+C_k}$$ where $\eta_k, K_k$ are some constants, where $\beta_1,\beta_2,\beta_3$ are such that $$\max 19\beta_1+0.5\beta_2+16\beta_3$$ is obtained.

I try to solve it by iterations, treating the ODE and the optimization as two separated steps. For every iteration,

First solve the ode to get some $C_k$;

then I fix $C_k$ and solve the linear programming to get $e_k$.

I repeat the above two steps to iteratively get next $C_k$, "$e_k$", and so on. My main problem is I do not know how to properly guarantee $C_k\ge0$ when solving the ODE. Any reference to papers or software packages that enforce such constraint is very much appreciated!

Without being able to guarantee that constraint, the above iterative process does not seem to converge and give weird oscillations.


问:
Time sampling changes solution
I'm currently trying to solve a problem using numerical methods. The set-up is rather long, so I apologize in advance...

TL;DR: My solutions change depending on how big my steps are and I don't know why.

Set-up

The first part of this problem can be described by a Lagrangian $L(s,\theta)$, which is given by
$$L = \frac{1}{24}(L^2\dot\theta^2+4b^2\dot\theta^2+24s\dot\theta^2+12\dot s^2+ 12 b\dot s \dot\theta+24gs\sin\theta -12bg\cos\theta).$$
This Lagragian gives, by means of the Euler-Lagrange equations,
$$
\begin{align*}
 \dot\theta^2 +g \sin\theta &=\ddot s + \frac{1}{2}b\ddot\theta\\
\left(\frac{1}{12}L^2+\frac{1}{3}b^2\right)\ddot\theta+\frac{1}{2}b\ddot s &= \frac{1}{2}bg\sin\theta-2\dot s\dot\theta-gs\cos\theta,
\end{align*}
$$
where you can assume $L=0.1,\,b=0.02$ and $g=9.81$.
I also have the initial conditions given by:
$$s(0)=0.1,\quad \dot s(0)= 0.3,\quad u(0)= 0,\quad \dot u(0)= 0.$$
The generalized coordinates $s,\theta$ and the real trajectory are connected over
$$x = \begin{pmatrix}s\cos\theta + \frac{b}{2}\sin\theta\\ -s\sin\theta +\frac{b}{2}\cos\theta\end{pmatrix}.$$

The first part of the problem is to solve this set of differential equation up to the time $t^*$ when the condition $\ddot x(t^*) = 0$ is fulfilled, or in terms of $s$ and $\theta$ until we have
$$-2 \sin\theta(t^*) \dot s(t^*)\dot\theta(t^*) - \cos\theta(t^*) s(t^*)\dot\theta^2(t^*) - \frac{1}{2}b\sin\theta(t^*) \dot\theta^2(t^*) + \cos\theta(t^*) \ddot s(t^*) + \frac{1}{2} b\ddot\theta \cos\theta(t^*) -s(t^*) \ddot\theta(t^*) \sin\theta(t^*) \overset{!}{=} 0.$$
After this time $t^*$ the new trajectory is given by

$$\vec{x}= \begin{pmatrix}x\\ - \left[x\tan(|\alpha|)+\frac{gx^2}{2v_0\cos^2(|\alpha|)}\right]\end{pmatrix} \quad\text{and}\quad \theta(t)= \dot \theta(t^*) \cdot t + \theta(t^*)$$
where $\alpha = \theta(t^*)$ and $v_0 = \dot s(t^*)$.

My attempt to solve this goes as follows: First defined

$$y = \begin{bmatrix}s\\\dot s\\ \theta\\ \dot\theta\end{bmatrix} \Longrightarrow \dot y = \begin{bmatrix}\dot s\\\ddot s\\\dot \theta\\ \ddot\theta\end{bmatrix}$$
I solved the Euler-Lagrange eq. for $\ddot \theta$ and $\ddot s$ and inserted them into $\dot y$ which gives rhs

```<code>def rhs(y, t, g=9.81, L=0.01):
    return np.array([y[1],
                    (1/(b**2 + L**2 + 24.0**y[0]**2))*
                    (6.0*b*g*np.cos(y[2])*y[0] - 7.0 *b**2 *g*np.sin(y[2]) - g* L*2.0 *np.sin(y[2]) - 24.0 *g *y[0]**2*np.sin(y[2])
                        - 24.0*b* y[0]* y[1]* y[3] - 8.0* b**2* y[0] * y[3]**2 - 2.0* L**2* y[0]* y[3]**2 - 48.0*y[0]**3* y[3]**2),
                    y[3],
                    -12.0*(-g *np.cos(y[2]) * y[0] + b* g* np.sin(y[2]) + 4* y[0] * y[1]* y[3] + b *y[0] *y[3]**2)/(b**2 + L**2 + 24 *y[0]**2)
                    ])
</code>```

I then used scipy's odeintalgorithm to solve the problem in the follwoing way

```<code>T = np.linspace(0,10,10000)
b = 0.02
lstart = 0.02
ystart = np.array([np.sqrt(lstart**2+ (b/2.0)**2), v, 0.0, 0.0])
y, data = odeint(rhs, ystart, T, full_output = True, args=(g, L))
</code>```

This worked out and plotting the solutions gives results that I expected.

The next thing to do was to figure out $t^*$. For this I interpolated the solutions and inserted them into the defining condition to find $t^*$. The following function is supposed to return the index of T at which this happens

```<code>def find_tstar_index(T, y):
    # Find point where bread leaves table
    s_spl = UnivariateSpline(T, y[:, 0], s=0, k=4)
    u_spl = UnivariateSpline(T, y[:, 2], s=0, k=4)

    s_spl_d = s_spl.derivative(n=1)
    s_spl_2d = s_spl.derivative(n=2)

    u_spl_d = u_spl.derivative(n=1)
    u_spl_2d = u_spl.derivative(n=2)

    searching_for_zeros = -2.0*np.sin(u_spl(T))*s_spl_d(T)*u_spl_d(T) \
                        - np.cos(u_spl(T))* s_spl(T)* u_spl_d(T)**2 \
                        -1/2.0* b* np.sin(u_spl(T))*u_spl_d(T)**2  \
                        + np.cos(u_spl(T))*s_spl_2d(T)\
                        + 1/2.0* b* np.cos(u_spl(T))* u_spl_2d(T) \
                        - s_spl(T)*np.sin(u_spl(T))*u_spl_2d(T)

    tstar_index = np.where(searching_for_zeros < 1e-25)[0][0]
    return tstar_index

tstar_index = find_tstar_index(T, y)
</code>```

The next step would be to set up $\vec{x}$ and $\theta (t)$ using $t^*$

The actual problem starts here: I now definie the follwoing function

```<code>angle = lambda t: y[tstar_index, 3] * t + y[tstar_index, 2]
</code>```

The thing is now that depending on the T defined in the beginning of my code the function angle changes significantly. For example

```<code>T_A = np.linspace(0,1,10000)
T_A = np.linspace(0,10,10000)
</code>```

give completely different results for the intervall [0,1] when plugged into the respective angle function... There has to be an error in my logic but after hours of trying to figure out what the problem is I'm kind of lost...

Complete Code

I'll post down here my complete code so that you can try it out yourself, but up to line 142 everything seems to work fine (T_fall is the same array regardless of what T is...)

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

g = 9.81
L = 0.1
#T = np.linspace(0,10,10000)
b = 0.02
lstart = 0.02
table_hight = 1.3

vs = [0.7]

def find_tstar_index(T, y):
    # Find point where bread leaves table
    s_spl = UnivariateSpline(T, y[:, 0], s=0, k=4)
    u_spl = UnivariateSpline(T, y[:, 2], s=0, k=4)

    s_spl_d = s_spl.derivative(n=1)
    s_spl_2d = s_spl.derivative(n=2)

    u_spl_d = u_spl.derivative(n=1)
    u_spl_2d = u_spl.derivative(n=2)

    searching_for_zeros = -2.0*np.sin(u_spl(T))*s_spl_d(T)*u_spl_d(T) \
                        - np.cos(u_spl(T))* s_spl(T)* u_spl_d(T)**2 \
                        -1/2.0* b* np.sin(u_spl(T))*u_spl_d(T)**2  \
                        + np.cos(u_spl(T))*s_spl_2d(T)\
                        + 1/2.0* b* np.cos(u_spl(T))* u_spl_2d(T) \
                        - s_spl(T)*np.sin(u_spl(T))*u_spl_2d(T)

    tstar_index = np.where(searching_for_zeros < 1e-25)[0][0]
    return tstar_index

def wurfparabel (x, v, beta):
    b = np.abs(beta)
    return - (x * np.tan(b) + (g * x**2)/(2* v**2* np.cos(b)**2))

def draw_bread(x_p, u, col):

    x_b = x_p + np.array([-b/2.0 * np.sin(u), -b/2.0 * np.cos(u)])
    x_left_down = x_b + np.array([-L/2.0 * np.cos(u), L/2.0 * np.sin(u)])
    x_left_up = x_left_down + np.array([b * np.sin(u), b * np.cos(u)])
    x_right_up = x_left_up + np.array([L * np.cos(u), -L * np.sin(u)])
    x_right_down = x_right_up + np.array([- b * np.sin(u), - b * np.cos(u)])

    plt.plot([x_b[0], x_left_down[0], x_left_up[0], x_right_up[0], x_right_down[0], x_b[0]],
        [x_b[1], x_left_down[1], x_left_up[1], x_right_up[1], x_right_down[1], x_b[1]], col)

    pass

def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(np.abs(array) - value)).argmin()
    return idx

def x_p(y):
    return np.array([y[:, 0] * np.cos(y[:, 2]) + b/2.0 * np.sin(y[:, 2]),
                    - y[:, 0] * np.sin(y[:, 2] + b/2.0 * np.cos(y[:, 2]))
                    ])

plt.figure()
def test(T):

    def bread_falling_trajectory(v):
        """
        Function that calculates the trajectory of a piece of toast falling from a table with the initial velocity v.
        """

        # Define starting values. The  first value is chosen in such a way that the Center of mass is of the same hight as the zero point.
        ystart = np.array([np.sqrt(lstart**2+ (b/2.0)**2), v, 0.0, 0.0])

        x_sol_start = np.array([ystart[0] * np.cos(ystart[2]) + b/2.0 * np.sin(ystart[2]),
                        - ystart[0] * np.sin(ystart[2] + b/2.0 * np.cos(ystart[2]))
                        ])

        # Define the right hand side of the differential equations.
        def rhs(y, t, g, L):
            return np.array([y[1],
                            (1/(b**2 + L**2 + 24.0**y[0]**2))*
                            (6.0*b*g*np.cos(y[2])*y[0] - 7.0 *b**2 *g*np.sin(y[2]) - g* L*2.0 *np.sin(y[2]) - 24.0 *g *y[0]**2*np.sin(y[2])
                                - 24.0*b* y[0]* y[1]* y[3] - 8.0* b**2* y[0] * y[3]**2 - 2.0* L**2* y[0]* y[3]**2 - 48.0*y[0]**3* y[3]**2),
                            y[3],
                            -12.0*(-g *np.cos(y[2]) * y[0] + b* g* np.sin(y[2]) + 4* y[0] * y[1]* y[3] + b *y[0] *y[3]**2)/(b**2 + L**2 + 24 *y[0]**2)
                            ])

        # Solve the ODE
        y, data = odeint(rhs, ystart, T, full_output = True, args=(g, L))

        # Get the trajectory
        x_sol = x_p(y)

        #Find the index of x_sol at which the bread leaves the table.
        tstar_index = find_tstar_index(T, y)

        # Lösung ab Freier Fall
        # EQM is given by \ddot theta = 0, solving this and inserting the boundry conditions leads to winkel
        winkel = lambda zeit: y[tstar_index, 3] * zeit + y[tstar_index, 2]
        print(y[tstar_index, 3], y[tstar_index, 2] )

        # for the spatial part we essentially have a free falling rod, so we solve the eom and construct trajectory
        v0 = 1/2.0 * np.sqrt(4*y[tstar_index,1]**2 +4*b*y[tstar_index,1]* y[tstar_index,3] + (b**2 + 4*y[tstar_index,0]**2)*y[tstar_index,3]**2)
        beta = y[tstar_index, 2]
        trajectory_ = lambda x: wurfparabel(x, v0, beta)
        y_0 = np.abs(x_p(y)[1][tstar_index]-trajectory_(x_p(y)[0][tstar_index]))

        trajectory = lambda x: trajectory_(x)  + y_0

        # Creat a plot of the results
        x_free_fall = np.linspace(x_sol[0][tstar_index], 0.45, N)

        sol_dict = {'x_sol_start': x_sol_start, 'table_traje': x_sol,
                    'x_free_fall': x_free_fall, 'free_fall_traj': trajectory(x_free_fall), 'y': y}

        return sol_dict, tstar_index, winkel

    for k, v in enumerate(vs):
        dicct, index, winkel = bread_falling_trajectory(v)
        # Plotting all the results
        #plt.subplot(3,1,k+1)
        plt.plot(dicct['x_sol_start'][0],dicct['x_sol_start'][1], "b*")
        # Plot the points where the bread leaves the table.
        plt.plot(dicct['table_traje'][0][index], dicct['table_traje'][1][index], "bo")
        # Plot the trajectroy up until this point
        plt.plot(dicct['table_traje'][0][:index+1], dicct['table_traje'][1][:index+1], label=r"<span class="math-container">$v$</span>={}".format(v))
        # Draw the bread
        draw_bread(dicct['table_traje'][:,0:index:100], np.pi/6.5 + dicct['y'][0:index:100,2], "g")

        x_p_free_fall = np.array([dicct['x_free_fall'],dicct['free_fall_traj']])

        spec_index = find_nearest(x_p_free_fall[1], table_hight)

        delta_T = T[-1]/N
        T_fall = np.linspace(0, spec_index * delta_T, spec_index)

        angle_free_fall = np.pi/6.5 + winkel(T_fall)

        draw_bread(x_p_free_fall[:,0:spec_index:100000], angle_free_fall[0:spec_index:100000], "r")

        # Plot the trajectroy after the bread has left the table
        plt.plot(dicct['x_free_fall'], dicct['free_fall_traj'])

        plt.ylim(-1.5, 0.1)
        plt.xlim(-0.1, 1.5)
        plt.plot(0.0, 0.0, "b*")
        plt.xlabel(r"<span class="math-container">$x$</span> [m]")
        plt.ylabel(r"<span class="math-container">$y$</span> [m]")
        plt.legend()

    return T_fall, winkel(T_fall)
N=1000000
ta, winekla = test(np.linspace(0,1,N))
tb, wineklb = test(np.linspace(0,10,N))
plt.show()
</code>```

Note that the huge test function is just to play around with the different T.

EDIT: To illustrate the problem, two plots with different time samplings (N=1000,10000,...) and different time lengths (the first picture for T_end = 7s and the second one for T_end=10s).

Title: Time sampling changes solution


问:
Numerical methods for the continuity equation with Sobolev vector field
Consider the continuity equation
$$
\partial_t \rho(x,t) + \operatorname{div} (b(x,t) \rho(x,t)) = 0, \qquad t \in [0,T], \quad x \in \mathbb R^N,
$$
with $b \in L^1((0,T), W^{1,p}(\mathbb R^N))$.

Have any numerical methods been used to approximate the (unique renormalized) distributional solution of this problem? Can you point out any references on this topic?

What if we replace $W^{1,p}$ with $BV$?

Clearly, if $b$ is smooth, the problem is much more classical.

I've asked a related question on MathOverflow and a specific implementation problem (with a specific $b$ of bounded variation) on Mathematica StackExchange.


问:
Boundary conditions for a Non-linear Schrödinger equation using an extended crank nicolson scheme
I try to solve numerically the following PDE for $E(r, z)$ with a cylindrical symmetrie (i. e. $E(r, z) = E(-r, z)$).

$\frac{\partial E}{\partial z} = \frac{i}{2k} \Delta E + \mathcal{N}(E)$

Where $\Delta$ is the Laplace operator in transversal direction and $k$ a real number. I want to use the Crank-Nicolson scheme to treat the Laplace operator and the Adams-Bashforth scheme to treat the nonlinearity ($\mathcal N(E)$). $r$ is defined by $r_j = 0 + j\Delta r$ for $j = 0 \dots N$. Since the Laplace in cylinder coordinates is given by: $\Delta E = \frac{\partial^2E}{\partial r^2} +  \frac{1}{r}\frac{\partial E}{\partial r} $ one gets the following representation of the Laplace operator:

$\Delta E_j^n = E^n_{j-1} -2 E^n_j + E^n_{j+1} + \frac{1}{2j}(E^n_{j+1} - E^n_{j-1})$

The given PDE equation is therefore given by the following, where the nonlinearity is treated by the Adams-Bashforth scheme.

$E_j^{n+1}  - E_j^n = i\delta (\Delta_j E^{n+1}_j + \Delta_j E^{n}_j) + (3/2 \mathcal N^n_j - 1/2 N^{n-1}_{j})$

where $\delta = \frac{\Delta z}{4 k \Delta r^2}$. From this $E^{n+1}_j $ can be expressed the following way:

$E^{n+1}_j = L_-^{-1} [L_+E^n_j + 3/2\mathcal N^n_j - 1/2\mathcal N^{n-1}_j]$ with the following matrixes $L_-$ and $L_+$.

$ L_{\pm} = \left( \begin{array}{rrrr}
1\mp2i\delta  & \pm i\delta v_0   \\
\pm i\delta u_1 & 1\mp2i\delta & \pm i\delta v_1\\
 &  &  &  & \\
 &  &  &  & \\
 &  &  & \pm i\delta u_N & 1\mp2i\delta\\
\end{array}\right)
$

Where $u_j = 1 - 1/(2j)$ and $v_j = 1 + 1/(2j)$.

For my problem the following boundary conditions are given: $ E(r = r_{max}, z) = 0 $ and $\frac{\partial E}{\partial r} |_{r=0}$. The first one for $r=r_{max}$ I can easily incorporate by changing the last row of $ L_{\pm} $.

The other one ($\frac{\partial E}{\partial r} |_{r=0}$) gives me trouble. I somehow have to change the entries of the first row of $L_{\pm}$, but I don't know how. I got some working boundary conditions for the case the nonlinearity is 0, but they break as soon as I add a nonlinearity.

I am greatful for any help.


问:
Computing the structure factor from positions and radial distribution function
I'm currently analysing some spatial point patterns that come from some fluid dynamics simulations and I'm having some difficulty computing the structure factor, $S(\pmb{k})$, from both the positions of the points and the radial distribution function, $g(\pmb{r})$. I have been following the structure factor Wikipedia page and also this paper on arxiv but I seem to be getting results that don't make sense.

Below is my class object to compute the radial distribution function for said spatial point patterns (I've omitted the help data and assertions from the code for clarity). The class takes in a numpy array or list of numpy arrays, finds the locations of the nonzero values, computes the distances between each nonzero value and then computes the radial distribution. Note that here, we consider an array cell a point if it has a nonzero value

```<code>import numpy as np
import bisect

class pair_correlation_function():

    def __init__(self,data,annulus_width,boundary=None):

        self.data = [data] if type(data) is np.ndarray else data
        self.dr = annulus_width
        self.boundary = boundary

    def get_positions(self):

        positions = list()

        for item in self.data:

            rows,columns = np.nonzero(item)
            positions.append(np.array(list(zip(rows,columns))))

        return positions

    def RDF(self):

        positions = self.get_positions()

        Lx,Ly = np.shape(self.data[0])
        area = Lx*Ly
        radii = list(x for x in range(int(Lx/(2*self.dr))))

        all_radial_distributions = list()

        for item in positions:

            if self.boundary == 'periodic':

                item_new = np.vstack(np.array([np.abs(item[k]-item[k+1:]) for k in range(len(item)-1)]))
                item_new[item_new > Lx/2] = Lx-item_new[item_new > Lx/2]
                norms = [np.sqrt((position**2).sum(axis=0)) for position in item_new]
                norms = sorted(norms)

            else:

                norms = [np.sqrt(((item[k]-item[k+1:])**2).sum(axis=1)) for k in range(len(item)-1)]
                norms = sorted(np.hstack(norms))

            number_particles = len(item)
            item_rdf = list()

            for r in radii:

                i = bisect.bisect_left(norms,r*self.dr)
                j = bisect.bisect_left(norms,(r+1)*self.dr)
                if i != len(norms) and j != len(norms): particle_count = len(norms[i:j])

                normalisation = (2*r+self.dr)*np.pi*self.dr*number_particles**2
                bin_value     = 2*area*particle_count/normalisation

                item_rdf.append(bin_value)

            all_radial_distributions.append(item_rdf)

        radii = np.array(radii)*self.dr

    return radii,all_radial_distributions
</code>```

Below is a sample array (a triangular lattice) with which to play around (we are assuming periodic boundary conditions here because my actual data comes from a simulation with periodic boundaries)

```<code>tri_period = np.array([[0,0,0,0],[0,1,0,1],[0,0,0,0],[1,0,1,0]])
triangular_lattice = np.tile(tri_period,(16,16))
dr = .025

radii, RDF = pair_correlation_function(triangular_lattice,dr,'periodic').RDF()
</code>```

and plotting gives

Here, we see a series of Bragg peaks as expected (though I'm not sure why the peaks segment into smooth, monotone decreasing bands? My actual data doesn't do this though and looks correct.) Now, if I try to compute the structure factor given the positions (see arxiv paper above, equation $(2)$)

```<code>def sf_from_positions(positions,box_dimension):

    sf = list()
    modes = list(x for x in range(1,int(box_dimension)))

    dk = 2*np.pi/box_dimension

    for h in modes:

        k_vec = np.array([1,1])*dk*h
        summation = 0

        for position in positions:

            summation += np.exp(-1j*k_vec.dot(position))

        sf.append(abs(summation)**2/len(positions))

    return sf, modes

rows,columns = np.nonzero(triangular_lattice)
positions = list(zip(rows,columns))
dimension = np.shape(triangular_lattice)[0]
factor, wave_modes = sf_from_positions(positions,dimension)
</code>```

and plotting gives

which looks off, because the Fourier transform of Bragg peaks should also be Bragg peaks (see this arxiv paper, page 12).

So, my first question is, what is going wrong in my computation of the structure factor? Can anyone see what I have done wrong? And secondly, I was wondering how I would integrate the radial distribution function to get the structure factor, through the formula $$S(\pmb{k}) = 1 + 4 \pi \rho \int_{\mathbb{R^{n}}} [g(\pmb{r})-1] e^{-i \pmb{k} \cdot \pmb{r}} d \pmb{r} $$ I tried to use the cumtrapz module in scipy but it gave me a structure factor that was negative (impossible) and that oscillated with an increasing wave packet.

Sorry if the question is overly involved. Thanks for your help.


问:
Parameter estimation using fmincon
This is a follow up to my previous question posted here.
I am solving an optimization problem using fmincon in MATLAB. There are no equality constraints in my model.  I want to pass ODEs as constraint to fmincon using, non-linear constraint, nlcon argument.
The function model  contains odes

```<code>function dz = model(t, z, c)
              dz(1) = ..
              dz(2) = ..
</code>```

I'd like to ask how the ode intergrator has to be called in the place of nlcon in fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon).

```<code>%p0 = Initial Value of Parameters
A = [];
b = [];
Aeq = [];
beq = [];
nlcon  = ode45(@(t,z)model(t,z,x), tSpan, z0); % Is this correct?
p = fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon);
</code>```

Since the ode function is already called in nlcon , is it required to call the function model again inside the objective function?

```<code>function cost = objective(c,expZ,tSpan,z0)
sol = ode45(@(t,z)model(t,z,c), tSpan, z0); % Is this step required?
ModelZ = sol(end,:)%  the solutions obtained from model at end point of tSpan
cost = ModelZ-expZ; %expZ contains steady state values of z1,z2..
</code>```

Any help would be highly appreciated.

Title: Parameter estimation using fmincon