问:
Strain propagation from surface to bulk in COMSOL
I am trying to simulate strain propagation from the surface into the bulk. I have a rectangular semiconductor block (~2 μm thick) on top of which metal gates (~25 nm thick) are deposited as seen in the picture below of xz-plane crosssection (zoomed near the gates). The x-axis is from left to right, the z-axis is from bottom to top and the y-axis is going into the picture.

I am fixing the bottom boundary with a constraint, other surfaces are free and all initial values for all domains are zero. Using Linear Elastic Material > Initial Stress and Strain (applied in the Metal Gate domains) I am putting a strain of $\epsilon_{xx} = 10^{-2}$ and $\epsilon_{yy} = 10^{-2}$. After a stationary study, I have plotted strain $\epsilon_{xx}$ in the same picture below (red is positive and blue is negative).

According to strain propagation, the positive surface strain under the contacts (and negative surface strain under the gaps) are supposed to propagate all the way until the bottom (with an exponential decay). However, in my simulation, the strain reverses at a certain depth. Even if I change my initial strain from $10^{-2}$ to $10^{-3}$, this inversion still occurs at that depth. This is not supposed to happen.


问:
Governing equations vs Transport equations
This is a basic question. But I did not find any explanations for this. How are governing equations, like mass, momentum, energy conservations equations, different from 'Transport equation'?. Is a conservation equation, like energy equation (heat equation), a transport equation for the temperature?  </p>

ContentLicense: CC BY-SA 4.0


问:
Why does GMRES converge much slower for large Dirichlet boundary conditions?
I'm trying to numerically solve a simple Laplace equation in 2D, with a nonlinear source term:

$\nabla^2 u = u^2$

with boundary conditions as $u=0$ everywhere except for $y=1$ where $u=u_0$. I'm using scipy's newton_krylov solver (method="lgmres") to minimized the discretized equation using the finite-difference method.

Here's the problem: for relatively small values of $u_0$ such as $u_0=1$, LGMRES converges relatively fast with fewer than 7 iterations. However, when I increase $u_0$ to $10$ or $100$, the the solver converges much slower and needs on the order of O(1000) iterations to converge.

Is this an expected behaviour of Newton-Krylov solvers? If so, what can I do to alleviate the issue?


问:
least squared optimization
I want to decompose a list of 3D vectors  $X_j$ as linear combination of five 3D verctors $C_k$

$$X_j= \sum_{i=1}^{5}{w_{ji}C_i}$$ both $X_j$ and $C_i$ are 3 components vectors

$$C= \begin{bmatrix} C_{11} & C_{12}  & C_{13}  \\ C_{21} & C_{22}  & C_{23} \\ ... \\  \\ C_{51} & C_{52}  & C_{53} \end{bmatrix}$$

$$W= \begin{bmatrix} W_{11} & W_{12} &...   & W_{15}  \\ W_{21} & W_{22} &...   & W_{25}  \\ ... \\  \\ W_{m1} & W_{m2} &...   & W_{m5}  \end{bmatrix}$$

$$X= \begin{bmatrix} X_{11} & X_{12}  & X_{13}  \\ X_{21} & X_{22}  & X_{23} \\ ... \\  \\ X_{m1} & X_{m2}  & X_{n3} \end{bmatrix}$$
i think this is a least squared optimization problem, but i am unable to correctly formulate it in the standard form $\left \|b-Ax \right \|^2$ and implement it in cvxpy.

$$\left \| X_j- \sum_{i=1}^{5}{w_{ji}C_i}  \right \|^2$$

Title: least squared optimization


问:
Efficient algorithm to determine the intersection volume of simple convex polyhedra
TLDR: Is there an efficient algorithm to compute the intersection of polyhedra with 8 or fewer vertices?

I have two sets of FEM meshes for one geometry (one exhibiting a skin effect). I have to transfer a field output of one to the other by a python script. The field output and desired input are both element dependent variables, so builtin nodal transferral would interpolate and extrapolate, losing me some precision. Therefore I need to check the meshes against intersection (which can be done by the Separating Axes Theorem), compute interesting volumes, times them by the volumetric output and sum them up for each target element.

The 2D case has been touched on, but I've only been able to find papers on the 3d problem concerned with arbitrarily complex polyhedra, requiring equally complex data structures to even reference the geometry.

I'm dealing only in hexahedra, tetrahedra, and wedge elements, so I'm hoping that someone's pinned down an algorithm for these 'simple' cases.


问:
WENO scheme on curvilinear coordinates
I've been developing a curvilinear FVM code. So far I've implemented the PPM scheme and am looking into adding WENO schemes. So far I've been discretizing the grid metrics using a second-order central.

For the WENO scheme will I need to use a higher-order scheme (5th order?) for the grid metrics. When I use the WENO scheme in curvilinear coordinates I get numerical oscillations, which makes me believe it has to do with the grid metrics. I've implemented the WENO scheme on a purely cartesian grid to verify the implementation was done right.

Title: WENO scheme on curvilinear coordinates


问:
Integrators for Nonlinear/Stiff PDE
It was suggested I ask this question in this section. Anyway:

I have a particular nonlinear PDE of the form

$$
u_t(x,t)=iu_{xx}(x,t)+f(x,u(x,t))
\tag{1}
$$

Where f is some nonlinear function. With boundary/initial conditions

$$ u(0,t)=u(L,t)=0, u(x,0)\approx0
\tag{2}
$$

Where $L$ is the length of the region considered. The equation is homogeneous so the initial conditions cannot be zero everywhere. Thus the interior points need only be close to zero, just not exactly zero.

My first step to use finite difference on the derivative

$$
u_{xx}(x_l)=\frac{1}{\Delta x^2}(u_{l+1}-2u_l+u_{l-1})
\tag{3}
$$

Giving now a system of stiff nonlinear ODE

$$
\textbf u_t=iD\textbf u+\textbf f(x_l,\textbf u)
\tag{4}
$$

Where $D$ is is a matrix with $-2$ on the diagonal, $1$ on the first off diagonals, zeros everywhere else, and all divided by the spatial step size squared.

My most successful attempt is to use partial backwards Euler where the derivative is taken at the forward time step and the nonlinear part is taken at the "old" time step, that is:

$$
\frac{1}{\Delta t} (\textbf u^{n+1}-\textbf u^{n})=iD\textbf u^{n+1}+\textbf f(x_l,\textbf u^n)
\tag{5}
$$

Which yields the solution

$$
\textbf u^{n+1}=(\frac{1}{\Delta t}I-iD)^{-1}(\frac{1}{\Delta t}\textbf u^n+\textbf f(x_l,\textbf u^n))
\tag{6}
$$

Where $I$ is the identity matrix. The problem with this method is the solutions distort if I change the time step. I have no idea why.

I tried solving the fully implicit problem via Newton-Raphson but the solutions always converge to zero. I do not know why. It's possible I implemented it incorrectly.

I decided to move on to the method of lines where I would use some kind of established time integrator.

Because of my issues with iterative methods I've looked into non-iterative integrators. The first one I attempted was a Rosebrock method. This, however, gives the incorrect solutions (the Rosebrock solutions are symmetric, the real solution should be asymmetric). And like in my semi-implicit method, changing the time step also changes the solutions. Next, I attempted exponential integrators. I used both a one-stage and a seven-stage integrator. Both gave me the same hideous solution. It practically looks like noise.

I'm interested in using the HHT-$\alpha$ integrator, but I don't understand how to apply it to my system. Specifically, I was interested in using this paper. Could anyone explain it to me? How do I find the $M$ and $Q$ matrices? They have several distinct Jacobians, which one do I find if I take the Jacobian of the right-hand side of equation 4? How do I find the other Jacobian(s)?

When they perform the stability analysis on the simple equation $u_t=\lambda u$ I don't see how they even obtain the resultant equations. Considering this integrator was designed for rigid body mechanics I'm having difficulty translating the method to a more generic system of ODE.

And I know two of the papers are behind paywalls. I have the pdfs for both. Am I allowed to share them?

Does anyone know of some other integrators that might be useful?

Some properties that might be useful to know: The real and imaginary solutions constantly oscillate from positive to negative in time. It can be shown there is no steady-state solution for $u$. However, the term $|u|^2$ does reach a steady state.

MAJOR EDIT:

The full set of equations can be found in the paper "Effects of excitation inhomogeneity in semiconductor lasers pumped by an electron beam." As far as I can tell, this paper simply does not exist in PDF form on the internet. I had to scan a copy from my school's library. The equations to study are

$
\frac{\partial E}{\partial t}=\frac{\omega}{2}[\frac{i}{k^2}\frac{\partial ^2 E}{\partial x^2}+\frac{(1-i\chi_1)\epsilon ^"-\epsilon ^"_0}{\epsilon_r}E
\tag{7}]
$
$
\frac{\partial \epsilon^"}{\partial t}=\frac{\sigma}{k}g_0g(x)-(\frac{1}{\tau_1}+2\sigma I)\epsilon^"
\tag{8}
$
$
I=\chi_0 |E|^2
\tag{9}
$

making the change of variables to nondimensionalize the problem, I let $x=\xi/k$, $t=s\tau_1$, and $E=u_0U$ ($u_0=1 V/m$). This changes equations 7 and 8 to

$
\frac{\partial U}{\partial s}=\gamma(i\frac{\partial^2 U}{\partial \xi^2}+(\alpha \epsilon^"-\beta)U)
\tag{10}
$
$
\frac{\partial \epsilon^"}{\partial s}=G_0g(\xi/k)-(1+\delta|U|^2)\epsilon^"
\tag{11}
$

Where

$
\gamma=\frac{\tau_1\omega}{2}=1.07757*10^7
$

$
\alpha=(1-i\chi_1)/\epsilon_r=0.0775194 - 2.86484*10^{-6}*i
$

$
\beta=\epsilon^"_0/\epsilon_r=1.73648*10^{-7}
$

$
G_0=\tau_1\frac{\sigma}{k}*g_0=3.56747*10^6
$

$
\delta=2\tau_1\sigma\chi_0=0.0000111495
$

$
k=258195.
$

And the function g(x) is

$
g(x)=a*\exp[-b^2(x-\mu)^2]Sech(cx)
$

The coefficients are $a=465.076$, $b=824.326$, $\mu=5.629393×10-3$, $c=-5633.573$.

The domain goes from $x=0$ to $x=L=6*10^{-3}$.

The first thing to notice, is that if you use backwards euler on equation (11) you can find an exact expression (I'm dropping the double primes).

$
\epsilon^{n+1}_l=(\epsilon^n_l+\Delta sPG_0g(\xi_l/k))/(1+\Delta s(1+\delta |U^{n+1}_l|^2))
\tag{12}
$

Where I've added a new term P. This term is varied. At some values the solutions converge to zero, above a threshold, they (should) converge to a proper solution. For me, threshold seems to be about 2*10^-7.

If you apply the method I used in equation (6) you'll find it will converge with any initial condition you throw at it. For instance, setting the internal points to $U(\xi,0)=\exp(i\xi)$

Title: Integrators for Nonlinear/Stiff PDE


问:
Shape functions in Euler Bernoulli Beam Equation
Does anyone have a intuitive explanation of why Hermite polynomials have to be utilized as the shape functions in the FEM solution of the Euler Bernoulli Beam 4th order ODE?

I have been learning FEM on my own and can't figure out why any other 2nd order polynomial can't be utilized in the place of the cubic Hermite ones, especially because the weak form transfers two derivatives to the test functions...


问:
Partial differential equation with convolution integral
The GP equation  is

$$i\frac{\partial u}{\partial z}+\nabla^2u+|u|^2u+\int e^{-[(x-x')^2+(y-y')^2]}|u(x',y')|^2 dx'dy'u(x,y)=0$$

with Neumann boundary condition. The initial condition is given by a Gaussian function, and $x, y\in [-10,10], z>0$.

Question: How do we deal with partial differential equations with a convolution integral?


问:
Tracking channel states using Machine Learning
I am new in AI and would like to apply machine learning to estimate the channel states.
I have a set of data. It is a matrix of 10000*8. Each row of this matrix is regarding a time step, i.e., 1st row = current time step (t), 2nd row = next time step (t+1), and so on. Each column is related to one transmitter, and I have 8 different transmitters. Each transmitter at each time step estimates the channel status and assigns one value form the set {0,1,2,-1}. for example the n-th row has a form of [ -1 0 0 0 1 0 0 2].

Knowing the status of the channel in the time step t I would like to know what is the channel states in time step (t+1).

I used the MLP and got the mse = 0.04 but the activate function that I used is either 'tansig' or 'logsig'. the results then is in [0, 1] or [-1, 1]  and I don't know how to convert them to {-1, 0, 1, 2}. I have also applied LSTM but the mse error is 0.31. So I don't know why I get such a large mse. Any suggestions would be greatly appreciated.


问:
Can Taylor methods be used effectively on stiff ODEs?
Cleve Moler has stated that "all numerical methods for stiff odes are implicit." However, I don't know whether this statement is a mathematical fact, or an simply an observation. Moreover, many numerical methods for ODEs assume that the only thing we can do with the right-hand side is point evaluation. If, for example, we can use automatic differentiation on the right-hand side of $x' = f(x,t)$, we are led into Taylor methods for ODEs.

But can Taylor methods be used effectively on stiff equations? If not, are there implicit Taylor methods which patch this problem up?

Note: There is always a problem about the definition of stiff. Will Cleve's definition work in this context?

A problem is stiff if the solution being sought varies slowly, but there are nearby solutions that vary rapidly, so the numerical method must take small steps to obtain satisfactory results.


问:
How can I solve the matrix optimization problem where denominator and numerator are different?
I want to solve the following maximization problem in $\mathbf{X}\in {\mathbb{R}}^{{m} \times {n}}$

\begin{eqnarray}
\begin{split}
\quad\max_\mathbf{X} \frac{\mathbf{Tr}(\mathbf{X}^\top \mathbf{Q} \mathbf{X})}{\mathbf{Tr}((\mathbf{X-A})^\top \mathbf{P} (\mathbf{X-A}))}
\end{split}
\end{eqnarray}

where $\mathbf{A}\in {\mathbb{R}}^{{m} \times {n}}$, $\mathbf{P}\in {\mathbb{R}}^{{m} \times {m}}$ and $\mathbf{Q}\in {\mathbb{R}}^{{m} \times {m}}$ and $\mathbf{P}$, $\mathbf{Q}$  are positive definite matrix.

It is well known that if $\mathbf{A}=\mathbf{O}$, it is easily solved by eigenvalue problem. However, I am not sure how to solve such a problem.


问:
A non linear ode with boundary conditions at infinity
I want to solve the non-linear ODE

$$\frac{d^2}{dx^2}y=a(y+y^3)$$

With the boundary conditions that
$$\lim_{x\to \pm \infty} y(x) =0$$

I am not aware of any analytical method for solving this kind of problem. So, I want to solve it numerically. Which method would be suitable for this ODE? I am aware of different methods for finite boundary conditions.

Also, does anybody know if it has an analytical solution to which I could compare my numerical solution?


问:
Find all recurring subgraphs/patterns of maximal size in a single undirected, labeled, connected graph
I would like to identify all subgraphs of maximal size (maximum number of nodes) that are recurrent in a single undirected, labeled, connected graph. I provide exemples of input and expected output below.

I'm a freshman to graph theory so the terminology may not be the most appropriate. Please correct me if you read between the lines.

By labeled graph I mean the nodes are labeled and labels can be redundant.

By connected graph I mean all nodes are connected.

By undirected graph I mean edges are not oriented/directed.

By all subgraphs of maximal size I am not sure if I mean all possible non-overlapping isomorphisms.

So I have a single graph that I know will contain some symmetries (it contains subgraphs that will be repeated N times in the graph). And I want to find the repeated patterns.

I understand this can be a (very) complex problem but the graphs I have to use will never contain more than 150-200 nodes, so that's already something.

I have to do this using python 3.x and I can use any libraries. I have tried some brute-force detection but it gets complicated quickly and I feel like I shall use some clustering, maybe with trees representations. Possibly it can be solved with the right functions from packages likes networkX, but I can't figure it out atm.

Also, in case there is a central node, I don't really care if it is included in multiple subgroups.

I can't tell you how grateful I already am for any insights you may provide :)

Exemple 1

Exemple 2

Exemple 3

But this last criterion is really bonus at this point.


问:
What's the more efficient way to solve this matrix equation?
This is intended to be a more generic question not about a specific system. Given a hermitian matrix $H(x_1,\dots,x_n)$ depending non-linearly on some real parameters $x_1,\dots,x_n$. We want these to be determined numerically such that
$$H=0$$
and assume we already know that a set of parameters always exists such that this holds.

Would it be more efficient to directly solve the above equation or to optimise its Frobenius norm $||H||_F^2$ (which would just be a scalar)? Also, which method would you apply?


问:
Does Boost provide a template implementation of the wedge product?
Does the boost C++ library implement the computation of the wedge product?

The wedge product is mentioned here, but it is not very clear (to me at least) whether there is a template implementation of it available.


问:
FDM discretization of equation on the boundary
In order to simulate the following equation using FDM
$$u_t(t,x)-u_{xx}(t,x)=0, \quad (t,x) \in (0,1)\times (0,1)$$
$$(u_t(t,x)-u_{x}(t,x))\rvert_{x=0}=0, \quad t \in (0,1)$$
$$(u_t(t,x)+u_{x}(t,x))\rvert_{x=1}=0, \quad t \in (0,1)$$
$$u(0,x)=u_0, \quad x \in (0,1)$$
$$u(0,0) =\mu_0, \quad u(0,1)=\mu_1,$$
where $f,g,u_0$ are given functions and $u_1,u_2$ are constant. I'm struggling in the discretization to obtain the corresponding linear system.

For the first equation taking $u_j^n=u(t_n,x_j)$ and using explicit Euler scheme, I found
$$u_j^{n+1}=ru_{j-1}^n+(1-2r)u_j^n +ru_{j+1}^n,$$
with $r=\dfrac{\Delta t}{(\Delta x)^2}$, but I don't know how to discretize the second and third equations. For the $x$-derivatives I think I should use centered formula by adding two ghost points as in Neumann b.c, and the same for time derivatives but this is a bit complicated and yields much of equations.

Maybe I should treat the equations like a coupled system through the boundary but I don't know how to do that. What is the right FDM discretization. Thanks for any help.

PS: My background to FDM is very simple (just theoretical) and this is my first time to simulate PDE. So, what is your advice to this in a right way? (methods, languages, programs...etc).


Title: Check if two points are symmetrics/asymmetrics

Body:
I am working on an Android app which lets a group of children draw whatever they want in a specific area. I need to check if the lines and figure that they draw are symmetric. The problem is that since the app is designed to let them draw freehand, I don't know how to compare the lines by considering certain flexibility. For example, I would consider symmetric these two "triangles" in the image:

even if they have not the same inclination and the same length. Do you have any suggestion? These are the points of the two shapes, do you notice any constant or something similar I can use to find what I need?

```<code>1303.9062, 245.88171
1306.9531, 239.48767
1309.055, 234.89392
1309.9219, 231.82782
1312.2241, 231.00659
1312.8906, 228.64291
1312.8906, 226.54477
1315.9375, 227.98535
1315.9375, 231.00659
1319.2817, 236.99414
1319.9219, 243.00916
1322.1042, 244.02539
1322.8906, 248.47485
1322.8906, 251.52356
1322.8906, 255.54462
1324.9263, 256.98926
1325.9375, 260.74615
1325.9375, 265.4226
1325.9375, 267.45148
1325.9375, 269.0193
1330.1698, 269.0193
1328.9062, 275.0343
1328.9062, 275.9956
1328.9062, 279.01685
1328.9062, 281.98315
1331.9531, 285.0044
1331.9531, 285.0044
</code>```

Triangle in the bottom:

```<code>1133.9062, 950.9956
1133.9062, 955.37256
1136.1552, 959.9429
1136.9531, 965.0901
1138.6583, 968.71606
1139.9219, 972.3617
1142.586, 976.9783
1143.9062, 982.42285
1146.608, 985.6449
1146.9531, 990.5189
1149.9219, 999.83105
1149.9219, 1002.55286
1153.7056, 1005.82275
1152.8906, 1010.1306
1154.862, 1011.9697
1155.9375, 1016.5016
1155.9375, 1018.0122
1155.9375, 1020.9785
1155.9375, 1024.9885
1159.9219, 1024.9885
1163.43, 1017.47314
1162.8906, 1013.48035
1165.9375, 1008.792
1167.0302, 1000.8804
1168.9062, 995.6416
1171.9531, 990.8357
1175.2189, 986.532
1175.9375, 981.4552
1178.9062, 975.46765
1181.9531, 970.21094
1181.9531, 966.84827
1186.4062, 966.9807
1184.9219, 962.5906
1184.9219, 960.99316
1187.8906, 960.99316
1187.8906, 956.98315
1187.8906, 956.98315
</code>```


问:
How to set up and solve acceleration-limited trajectory optimization problems?
I've been trying to learn how to solve simple acceleration-limited trajectory planning problems. I'm working in C++ and I've been using the Eigen library to do linear systems solving. I'm doing the nonlinear-programming part myself because that's part of what I want to learn.

Here's a blog post with some animated GIFs that show the sort of problem I'm working with. It's a one-dimensional trajectory made of two cubic segments, so the start, midpoint and endpoint positions are fixed, the start and end velocities are fixed, and there's a limit on maximum acceleration. The free variables are the durations of the two segments, and the midpoint velocity. I'm trying to minimize the total time of the maneuver.

For movement in a line like this the acceleration constraints can be broken down into separate constraints on minimum and maximum acceleration at each end of the two cubics, for a total of eight constraints. This works pretty well; I wrote an interior-point solver (using the Boyd and Vandenberghe book) and it seems to converge reliably.

I want to extend the problem to two dimensions, which would make the acceleration limit be a disc. For a first step I stayed in the one-dimensional problem, but replaced my eight constraints with four constraints based on squared acceleration at the segment endpoints.

However, with the set of constraints based on squared acceleration it's pretty easy to come up with problem states that the solver won't be able to make progress from. I think this is because the constraint gradient info is not helpful if you're on the side of the constraint function where the objective moves away from the constraint edge? It's also possible I'm just getting my math wrong.

Is there a better way to pose and solve this type of problem? Perhaps I should be using an active-set method, so the problematic constraints aren't in effect until the solution gets over to the side where they should apply?


问:
What are some good debugging habits for numerical simulation?
I'm currently writing a lid drive cavity CFD code on python. Currently, my code has some issues (values jumping bear b.c). I was wondering what are some good habits in debugging numerical codes. Hopefully, I would be able to learn some which I may apply to my future projects.</p>


问:
Simulate friction arround steady-state point
When simulating a friction dynamic system, I have come to a problem when approaching steady-state.

In the example below, a box moves along a surface with a certain friction coefficient:

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

import numpy as np
import matplotlib.pyplot as plt

G = 9.8
MU = 0.3

TIME_STEP = [0.1, 0.01, 0.001]
SIM_TIME = 10 # s

fig, ax = plt.subplots()

for time_step in TIME_STEP:
    t = np.linspace(0, SIM_TIME, int(SIM_TIME/time_step))
    speed = 10 # m/s
    s = []

    for i in t:
        speed += -np.sign(speed)*G*MU*time_step
        s.append(speed)

    ax.plot(t, s, label = "Time step {:1.3f} s".format(time_step))

ax.set_xlabel('Time (s)')
ax.set_ylabel('Speed (m/s)')
ax.legend(loc='upper right')
</code>```

Since the friction is only dependent on the speed unit vector, the system oscillates around the steady-state point instead of converging to 0. The oscillation vary depending on the time_step simulated, but it is there always.

The difficulty is that in the real system, the steady-state point is not known, so approaches like making the speed equal to zero when crossing the steady-state point are not applicable.

Is there a way for this to be coded in order to avoid those steady-state oscillations, and without loosing physical integrity?

EDIT

The equations for the given system are:

$v_{0}=10$

$\dot{v}=-\frac{v}{\left | v \right |}\cdot g\cdot \mu $


问:
What's a good numerical/optimization software package for solving the 2-D optimal stopping problem?
I am looking for a numerical software package to help me solve the 2-dimensional "free boundary" PDEs that arise in optimal stopping problems.  In one dimension a standard optimal stopping problem in Economics is when to exercise an expansion option.  In the end, solving this problem comes down to solving the following ODE

$$
\mu x u' + \frac{1}{2}\sigma^2 x^2 u'' + x k_0= r u
$$
such that
\begin{eqnarray}
u(0) & =& 0 \\
u(\overline{x}) & =&  \frac{\overline{x}k_1}{r-\mu} - p\\
u'(\overline{x}) & = &  \frac{k_1}{r-\mu}.
\end{eqnarray}
A solution to this problem is a function $u$ together with the location of the boundary $\overline{x}$ - this is the so called "free boundary". This problem has a closed form solution.  In two dimensions, the problem no longer has a closed form solution.  The problem I want to solve is
$$
 x k_0 + \mu x \frac{\partial u}{\partial x}+ rw\frac{\partial u}{\partial w} + \frac{1}{2} x^2 \left(\sigma_x^2\frac{\partial^2 u}{\partial x^2}  + \sigma_w^2 \frac{\partial^2 u}{\partial w^2} + 2\sigma_x\sigma_w \frac{\partial^2 u}{\partial x\partial w}\right)  = r u
$$
such that
\begin{eqnarray}
u(0,w) & = & 0\\
u(\overline{x}(w),w) & = & v(\overline{x}(w),w) - p \\
\frac{\partial}{\partial x}u(\overline{x}(w),w) & = & \frac{\partial}{\partial x}v(\overline{x}(w),w) \\
u(x,0) &= & \ell x \\
\frac{\partial}{\partial w}u(x,\overline{w}) & = & -1 \\
\frac{\partial^2}{\partial w^2}u(x,\overline{w}) & = & 0 \\
\end{eqnarray}

where $v$ is some known function of $x$ and $w$.  The solution of this problem will then be a function $u$ together with functions $\overline{x}(w)$ and $\overline{w}(x)$ that give the locations of the boundaries.

Do you know of any software packages that can solve these free boundary problems?  If so, can someone point me towards a reference or tutorial that might help?

ContentLicense: CC BY-SA 3.0


问:
Is there any catch on using `zgemm3m` vs regular `zgemm`?
I've just (to my embarrassment) encountered a BLAS-like extension of a matrix-matrix product subroutine gemm in Intel MKL: gemm3m. This subroutine (particular versions: cgemm3m and zgemm3m) allows performing matrix-matrix multiplication for complex-valued matrices using fewer arithmetic operations.

The gemm3m documentation claims that it

...reduces the time spent in matrix operations by 25%, resulting in significant savings in compute time for large matrices.

Looking at the provided error analysis in the Application Notes, I don't see anything "criminal":
$$
\hat{C}=\text{fl}(C_1+iC_2)=\text{fl}\big((A_1+iA_2)(B_1+iB_2)\big)=\hat{C}_1+i\hat{C}_2
$$
$$
||\hat{C}_1-C_1||\leq 2(n+1)u||A||_\infty||B||_\infty+\mathcal O(u^2)\\
||\hat{C}_2-C_2||\leq 4(n+4)u||A||_\infty||B||_\infty+\mathcal O(u^2)
$$
where $A,B,C\in\mathbb C^{n\times n}$ are complex matrices, $A_{1,2},B_{1,2},C_{1,2}\in\mathbb R^{n\times n}$ are their real and imaginary parts, respectively, $i=\sqrt{-1}$; $\hat{C}\in\mathbb C^{n\times n}$ and $\hat{C}_{1,2}\in\mathbb R^{n\times n}$ are the result of floating-point operations on $A$ and $B$ accoring to the gemm3m matrix-matrix multiplication algorithm. $|u|<\epsilon_\text{mach}$ if the floating-point arithmetic is IEEE-754 and no underflow\overflow happens.

So, is there any catch on using zgemm3m vs regular zgemm? Is there a situation where I should avoid using zgemm3m?


问:
Minimizing a polynomial with millions of monomials
I need to minimize a single polynomial $P(x_1,x_2,...,x_n)$ with the constraint that for each $i$, $0\leq x_i \leq 1$. The number of variables in my practical problem is at most $50$. The degree is at most $10$. The number of monomials is about 2 million.

Are there any tools that can handle such minimization problems? An approximate solution which is guaranteed to be off by some multiplicative factor would also be welcome.


问:
Is semi-Lagrangian 1D advection identical to upwind Euler, when $|u|<\Delta x/\Delta t$?
This looked true to me, I worked through the algebra for 1D and confirmed it. But no one seems to mention it, so maybe I'm missing something...

I'm using Bridson's SIGGRAPH 2007 course notes. [5MB pdf], p. 21 (page 33 in the pdf).

Semi-Lagrangian: back-trapolate the velocity field, to find the "particle" that would end up here ($x_i$), then linearly interpolate to get the value of that particle (pressure or whatever), which will become the value here ($q_i$).

Here, I assume velocity $u>0$ (to the right) and $\Delta t u < \Delta x$ (so the particle's previous position $x_P$ lies within the interval $[x_{i-1},x_i]$).

back-trapolate:
$$x_P = x_i - \Delta t u_i$$

interpolate:
$$\alpha = \frac{x_P - x_{i-1}}{\Delta x}$$

$$q_i^{t+1} = (1-\alpha)q_{i-1} + (\alpha)q_i$$

Upwind Euler
Choose forward or backward Euler, depending on the direction information is coming from, i.e. "upwind". This is stable, with my velocity restriction.

Reusing the above velocity restriction (to the right; moves less than a cell-width in a timestep).

$$q_i^{t+1} = q_i - \frac{\Delta t u_i(q_i - q_{i-1})}{\Delta x}$$

This is identical to the semi-Lagrangian above. (I won't go through all the algebra here, but I began by expanding the semi-Lagrangian into one line and used $x_i - x_{i-1} = \Delta x$).

They're also identical for $u<0$, just using the cell to the right i.e. $x_i$ and $x_{i+1}$ instead of $x_{i-1}$ and $x_i$.

Semi-Lagrangian advection enables jumping multiple cells (which I exclude with my velocity restriction), the only difference seems to be handling the velocity direction by determining which cell interval the particle lies in.

What am I missing? Perhaps higher-order "semi-Lagrangian" schemes don't correspond to upwind schemes so neatly?

Or maybe, jumping multiple cells is very useful sometimes? (In my case, I need the velocity restriction anyway).

I prefer upwind Euler, because it's simpler to calculate, less complex and easier to understand. But I've got into trouble before by "improving" schemes: my confidence exceeded my competence. So I'm checking this time.


问:
Long AMPL model preparation time
We deal with a large-scale linear optimization problem (~50000 variables and ~4000000 constraints). We use AMPL Studio modeling environment for problem modeling and then calling linear solver (CPLEX, Gurobi, etc.). However, before calling a solver, the preparation process (generation of LP or MPS file I guess?) of such model takes ~950 seconds.

Profiling shows that most of the time is consumed by AMPL's genmod function. Is it normal that preparation of such model (for the solver) takes so long?

If so, could you suggest an alternative way to generate model faster? As far as I am aware, AMPL is among the leader's for such job compared to other AMLs.

Title: Long AMPL model preparation time


问:
Regarding Amdahl's balanced system law
One of the paper titled "Rules of Thumb in Data Engineering" (Jim Gray et. el.) mentions some calculations based on Amdahl's balanced system law.

Link to paper: https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/ms_tr_99_100_rules_of_thumb_in_data_engineering.pdf

Initially, the paper provides some rules of thumb.

"10. Amdahl’s revised balanced system law: A system needs 8 MIPS/MBpsIO, but the instruction rate and IO rate must be measured on the relevant workload. (Sequential workloads tend to have low CPI (clocks per instruction), while random workloads tend to have higher CPI.)

Alpha (the MB/MIPS ratio) is rising from 1 to 4. This trend will likely continue.

Random IO’s happen about once each 50,000 instructions. Based on rule 10, sequential IOs are much larger and so the instructions per IO are much higher for sequential workloads."

Later, the paper provides following example.

"Amdahl’s balanced system law becomes more complex to interpret in the new world of quad-issue pipelined processors. Table 2 summarizes the following analysis. In theory, the current 550 MHz Intel processors are able to execute 2 billion instructions per second, so Amdahl’s IO law suggests that each 550 MHz processor needs 160 MBps of disk bandwidth (all numbers rounded). However, on real benchmarks, these processors demonstrate 1.2 clocks per instruction (CPI) on sequential workloads (TPC-D,H,R) and 2.2 clocks per instruction on random IO workloads (TPC-C, W) [7,8]. These larger CPIs translate to 450 MIPS on sequential and 260 MIPS on random workloads. In turn, Amdahl’s law says these processors need 60 MBps sequential IO bandwidth (~450/ 8) and 30 MBps random of IO bandwidth (~260/8) per cpu respectively (for tpcH and tpcC). A recent tpcH benchmark by HP [5] used eight 550 MHz processors with 176 disks. This translates to 22 disks per cpu, or about 70 MBps of raw disk bandwidth per cpu and 120 MBps of controller bandwidth per cpu (consistent with Amdahl’s prediction of 60MBps). Amdahl’s law predicts that system needs 30MBps of IO bandwidth. Using 8KB pages and 100 IO/s per disk implies 38 disks per processor – a number comparable to the 50 disks Dell actually used [4]."

I could not understand two parts of it (marked in bold). Need help to understand it.

1). How Amdahl's IO law suggests 160 MBps here (or is it 60MBps as indicated later (in bold))?

a. According to Amdahl's IO law, 50,000 instructions will need 1 IO (but of what size??)

b. Hence, 2 billion instructions (= 2 * 10^9 instrs = (50,000) * (40,000)) will need 40,000 IOs (but of what size?)

c. Should we multiply this 40,000 by 4 (because of quad-issue pipelined processor)? This will give 160,000 IOs.

d. Should we assume 1000 Bytes for each IO? If this is the case, it will give 160MBps!!

e. After this I am lost. Are my assumptions in (c) and (d) correct?Please help

2) How tpcH benchmark by HP translates to 70MBps of raw disk bandwidth per cpu and 120 MBps of controller bandwidth per cpu?

a. I could understand that if 176 disks are used with 8 processors, it is about 22 disks per cpu.

b. But, how this translates to 70MBps of raw disk bandwidth per cpu and 120MBps of controller bandwidth per cpu?

Title: Regarding Amdahl's balanced system law


问:
Solving a simple Shallow Water model
Hi. I have a question at Mathematics and they suggested post here, once it's not common. I transcript as following. Many thanks

I need to solve with basic methods this simple Shallow Water Model:

$$\begin{bmatrix}h\\ hv\end{bmatrix}_t+\begin{bmatrix}hv\\ hv^2+\frac{1}{2}gh^2\end{bmatrix}_x=\begin{bmatrix}0\\ 0\end{bmatrix}$$

where $h$ is the height of the water, $v$ is the horizontal velocity of the fluid and $g$ is the gravitational constant.

Writing in quasilinear form, I've got

$$\begin{bmatrix}h\\ hv\end{bmatrix}_t+\begin{bmatrix}0&1\\ gh-v^2&2v\end{bmatrix}\begin{bmatrix}h\\ hv\end{bmatrix}_x=\begin{bmatrix}0\\ 0\end{bmatrix}$$

The Jacobian can be diagonalized as $J=RDR^{-1}$ where

$$R=\begin{bmatrix}\dfrac{1}{v-\sqrt{gh}}&\dfrac{1}{v+\sqrt{gh}}\\ 1&1\end{bmatrix}$$

$$R^{-1}=\begin{bmatrix}\dfrac{v^2-gh}{2\sqrt{gh}}&\dfrac{1}{2}-\dfrac{v}{2\sqrt{gh}}\\ \dfrac{gh-v^2}{2\sqrt{gh}}&\dfrac{1}{2}\bigg(\dfrac{v}{\sqrt{gh}}+1\bigg)\end{bmatrix}$$

$$D=\begin{bmatrix} v-\sqrt{gh}&0 \\ 0 & v+\sqrt{gh}\end{bmatrix} $$

I've done the change of variables

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}=R^{-1}\begin{bmatrix}h\\ hv\end{bmatrix}$$

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}=\begin{bmatrix}\frac{h}{2}(v+\sqrt{gh})\\ \frac{h}{2}(v-\sqrt{gh})\end{bmatrix}$$

Geting the system

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}_t+\begin{bmatrix}v-\sqrt{gh}&0\\0&v+\sqrt{gh}\end{bmatrix}\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}_x=\begin{bmatrix}0\\ 0\end{bmatrix}$$

ie

$$\overline{h}_t+(v-\sqrt{gh})\overline{h}_x=0\\
\overline{k}_t+(v+\sqrt{gh})\overline{k}_x=0$$

Trying to solve this conservation laws by the characteristic method, I've got to the first

$$x(t)=(v_0(x_0)-\sqrt{gh_0(x_0)})t+x_0\qquad (i)$$

and to the second

$$x(t)=(v_0(\overline{x}_0)+\sqrt{gh_0(\overline{x}_0)})t+\overline{x}_0\qquad (ii)$$

where $v_0,h_0$ are the initial conditions.

To find the solution at time $(x,t)$, I need to solve $(i)$ for $x_0$ and take $\overline{h}(x,t)=\overline{h}_0(x_0)$ and solve $(ii)$ for $\overline{x}_0$, taking $\overline{k}(x,t)=\overline{k}_0(\overline{x}_0)$.

The solution in time $(x,t)$ may be the solution of the system

$$\begin{bmatrix}\overline{h}_0(x_0)\\ \overline{k}_0(\overline{x}_0)\end{bmatrix}=\begin{bmatrix}\frac{h}{2}(v+\sqrt{gh})\\ \frac{h}{2}(v-\sqrt{gh})\end{bmatrix} $$

That is

$$hv(x,t)=\overline{h}_0(x_0)+\overline{k}_0(\overline{x}_0)$$

$$h(x,t)=\bigg(\dfrac{\overline{h}_0(x_0)-\overline{k}_0(\overline{x}_0)}{\sqrt{g}}\bigg)^{2/3}$$

Or better

$$hv(x,t)=\frac{h_0(x_0)(v_0(x_0)+\sqrt{gh_0(x_0)})+h_0(\overline{x}_0)(v_0(\overline{x}_0)-\sqrt{gh_0(\overline{x}_0)})}{2}$$

$$h(x,t)=\bigg(\dfrac{h_0(x_0)(v_0(x_0)+\sqrt{gh_0(x_0)})-h_0(\overline{x}_0)(v_0(\overline{x}_0)-\sqrt{gh_0(\overline{x}_0)})}{2\sqrt{g}}\bigg)^{2/3}$$

I would like to know:

1) If my idea is correct (mainly about characteristics: is it valid? I think there is troubles with this part)

2) If these calculus are sufficient to "solve" the system. In other words, what mean the "solution" of the system? The question that I am dealing is only the following:

Find the solution to the model.

I think this is something vague once I don't have any data. What do you understand as the "solution of the system"?

3) And if the characteristics crosses or occur some problem with this method? Or if (i) and (ii) don't have solutions? With some simulations that I've tried the equations (i) and (ii) don't had solutions for a great number of points $(x,t)$. Without the initial conditions, how can I know about possible troubles?

Many thanks.

==========================

Edit

If we wish to study waves with very small amplitude, then we can linearize these equations around some constant state $(\hat{h},\hat{h}\hat{v})$ to obtain a linear system.

$$\begin{bmatrix}h\\ hv\end{bmatrix}_t+\begin{bmatrix}0&1\\ g\hat{h}-\hat{v}^2&2\hat{v}\end{bmatrix}\begin{bmatrix}h\\ hv\end{bmatrix}_x=\begin{bmatrix}0\\ 0\end{bmatrix}$$

Getting the system

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}_t+\begin{bmatrix}\hat{v}-\sqrt{g\hat{h}}&0\\0&\hat{v}+\sqrt{g\hat{h}}\end{bmatrix}\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}_x=\begin{bmatrix}0\\ 0\end{bmatrix}$$

where

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}=R^{-1}\begin{bmatrix}h\\ hv\end{bmatrix}=\begin{bmatrix}\dfrac{\hat{v}^2-g\hat{h}}{2\sqrt{g\hat{h}}}&\dfrac{1}{2}-\dfrac{\hat{v}}{2\sqrt{g\hat{h}}}\\ \dfrac{g\hat{h}-\hat{v}^2}{2\sqrt{g\hat{h}}}&\dfrac{1}{2}\bigg(\dfrac{\hat{v}}{\sqrt{g\hat{h}}}+1\bigg)\end{bmatrix}\begin{bmatrix}h\\ hv\end{bmatrix}$$

The solution is

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}(x,t)=\begin{bmatrix}\overline{h}_0(x-(\hat{v}-\sqrt{g\hat{h}})t)\\ \overline{k}_0(x-(\hat{v}+\sqrt{g\hat{h}})t)\end{bmatrix}$$

Or better

$$\begin{bmatrix}h\\ hv\end{bmatrix}(x,t)=R\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}(x,t)$$

$$\begin{bmatrix}h\\ hv\end{bmatrix}(x,t)=\begin{bmatrix}\dfrac{1}{\hat{v}-\sqrt{g\hat{h}}}&\dfrac{1}{\hat{v}+\sqrt{g\hat{h}}}\\ 1&1\end{bmatrix}\begin{bmatrix}\overline{h}_0(x-(\hat{v}-\sqrt{g\hat{h}})t)\\ \overline{k}_0(x-(\hat{v}+\sqrt{g\hat{h}})t)\end{bmatrix}$$

For instance, let take $\hat{h}=1$ and $\hat{v}=0$, considering $g=1$ (in some dimension).

So,

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}=\begin{bmatrix}\dfrac{-1}{2}&\dfrac{1}{2}\\ \dfrac{1}{2}&\dfrac{1}{2}\end{bmatrix}\begin{bmatrix}h\\ hv\end{bmatrix}=\begin{bmatrix}\dfrac{(hv-h)}{2}\\
\dfrac{(h+hv)}{2}\end{bmatrix}$$

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}_t+\begin{bmatrix}-1&0\\0&1\end{bmatrix}\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}_x=\begin{bmatrix}0\\ 0\end{bmatrix}$$

The solution is

$$\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}(x,t)=\begin{bmatrix}\overline{h}_0(x+t)\\ \overline{k}_0(x-t)\end{bmatrix}=\begin{bmatrix}\dfrac{(h_0v_0-h_0)(x+t)}{2}\\ \dfrac{(h_0v_0+h_0)(x-t)}{2}\end{bmatrix}$$

Or better

$$\begin{bmatrix}h\\ hv\end{bmatrix}(x,t)=R\begin{bmatrix}\overline{h}\\ \overline{k}\end{bmatrix}(x,t)=\begin{bmatrix}-1&1\\1&1 \end{bmatrix}\begin{bmatrix}\dfrac{(h_0v_0-h_0)(x+t)}{2}\\ \dfrac{(h_0v_0+h_0)(x-t)}{2}\end{bmatrix}$$

$$\begin{bmatrix}h\\ hv\end{bmatrix}(x,t)=\begin{bmatrix}\dfrac{(-h_0v_0+h_0)(x+t)+(h_0v_0+h_0)(x-t)}{2}\\ \dfrac{(h_0v_0-h_0)(x+t)+(h_0v_0+h_0)(x-t)}{2}\end{bmatrix}$$

Taking $v_0\equiv0$ and $h_0=1+\frac{2}{5}e^{-5x^2}$ as in 1, we have

$$\begin{bmatrix}h\\ hv\end{bmatrix}(x,t)=\begin{bmatrix}\dfrac{h_0(x+t)+h_0(x-t)}{2}\\ \dfrac{-h_0(x+t)+h_0(x-t)}{2}\end{bmatrix}$$

$$\begin{bmatrix}h\\ hv\end{bmatrix}(x,t)=\begin{bmatrix}\dfrac{1+\frac{2}{5}e^{-5(x+t)^2}+1+\frac{2}{5}e^{-5(x-t)^2}}{2}\\ \dfrac{-1-\frac{2}{5}e^{-5(x+t)^2}+1+\frac{2}{5}e^{-5(x-t)^2}}{2}\end{bmatrix}$$

$$\begin{bmatrix}h\\ hv\end{bmatrix}(x,t)=\begin{bmatrix}1+\frac{1}{5}(e^{-5(x+t)^2}+e^{-5(x-t)^2})\\ \frac{1}{5}(-e^{-5(x+t)^2}+e^{-5(x-t)^2})\end{bmatrix}$$

1 CROWHURST, Peter. Numerical Solutions of One-Dimensional Shallow Water Equations. Australian Mathematical Science Institute, School Of Computing And Mathematics, Charles Sturt University, Feb

Title: Solving a simple Shallow Water model


问:
Minimize number of math operation of a specific matrix vector multiplication?
Let's say we have a Matrix M and a column vector v like below

multiply

equals

Assume we can only perform multiplication, addition and substraction operation.
With normal approach we need 3 multipliction and 3 addition to calculate row1
row1 = v_1 + 2v_2 + 2v_3 + 4v_4.
However we can do it smartly by

```<code>t1 = v_1 + 2*v_3
t2 = v_2 + 2*v_4
row1 = t1 + 2 * t2
row2 = t2 - t1
row3 = ...
</code>```

Totally we only need 4 multiplication and 5 addition/substraction to do the calculation.

Is there any algorithm or approach to figure out the minimum number of operations?

P.S. I know sub-expression elimination algorithm can be used for this kind of problem.
However I am trying to figure out whether it can be solve by linear algebra method or explaination.


问:
Dividing Point Cloud into voxels
I am reading a paper on implementing convolution neural network for a 3D point cloud. In this paper, they are dividing the point cloud into voxels. Is there any easy way to do it using point cloud library?</p>

Title: Dividing Point Cloud into voxels


问:
How to do a Hoshen–Kopelman algorithm in Matlab with periodic boundary conditions?
Say I have an array s with elements of 0 and 1. From my understanding, a Hoshen–Kopelman algorithm for s will be (in 3 dimensions d=3):

```<code>label=labelmatrix(bwconncomp(s,2*d));
</code>```

If I want to do the same thing but now s will have periodic boundary conditions (specifically in 3 dimensions if it matters), how do I do that?

Edit: After searching some more I stumbled upon this:

https://uk.mathworks.com/matlabcentral/fileexchange/66079-cc2periodic?s_tid=mwa_osa_a

Edit: If this is off-topic then I would like to move this to stack overflow if possible.


问:
Find solution to Polynomial Sequences without going through variety in Sage
I'm computing the Groebner basis of an ideal defined over the QQ ring. Once I have this Groebner basis, I would like to obtain a set of values that satisfy the equations in the Groebner basis. I know that the full set is going to be the variety of the ideal, but since this object is huge, I might not be interested in finding those values.
In Maple, after computing the Groebner basis, I'm able to call the solve() method on it and even set a maximum number of solutions I want to obtain. As a small example:

```<code>P.<x,y,z,t>=PolynomialRing(QQ,4)
I = P.ideal(x*(x-1), y*(y-1), z*(z-1), t - x*y*z)
gb = I.groebner_basis()
</code>```

Here I could have called I.variety() or gb.variety() and obtained the same set of solutions:

```<code>sage: gbI.variety()
[{y: 0, z: 0, t: 0, x: 0},
 {y: 0, z: 0, t: 0, x: 1},
 {y: 1, z: 0, t: 0, x: 0},
 {y: 1, z: 0, t: 0, x: 1},
 {y: 0, z: 1, t: 0, x: 0},
 {y: 0, z: 1, t: 0, x: 1},
 {y: 1, z: 1, t: 0, x: 0},
 {y: 1, z: 1, t: 1, x: 1}]
</code>```

But I would like to know if it is possible to call a method like (I can do this in Maple):

```<code>solve(gb,[max_sol=2])
</code>```

Such that I can obtain a subset of the variety instead of the whole set. My motivation is that the size of the initial system of polynomials that I have is considerably larger than this example, and finding the feasible solutions on the reduced Groebner basis is more manageable. I might also not be interested in all the elements in the variety. Finally, if I transform the Groebner basis in an ideal itself and try to compute the variety on that object

```<code>gbI = ideal(gb)
gbI.variety()
</code>```

I find the following error

```<code>RuntimeError: error in Singular function call 'groebner':
int overflow in hilb 3
error occurred in or before standard.lib::stdhilb line 299: `    intvec hi = hilb( Id(1),1,W );`
expected intvec-expression. type 'help intvec;'
leaving standard.lib::stdhilb
leaving standard.lib::groebner
</code>```


问:
How to set the errorbars in a different color from the plot in Matlab?
When making a plot in Matlab using errorbar(...) the color of the error bars is the same as the plot. How do I set them to be in a different color?

I tried looking for a way to do it in here:

https://www.mathworks.com/help/matlab/ref/matlab.graphics.chart.primitive.errorbar-properties.html

https://www.mathworks.com/help/matlab/ref/errorbar.html

But I couldn't find it.


问:
Equi-order in pressure correction schme of Navier-Stokes equation
I am wondering if there is an stabilized equi-order scheme in pressure correction scheme in solving Navier-Stokes equation? Usually P2-P1 element combination is used to solve NS equation, and a PSPG stabilizing term could be added if NS equation is solved by P1-P1 element in a coupled way. My question is, is there a conterpart in pressure correction equation?</p>


问:
Can I use matplotlib to plot the surface of a 3D body?
If Matplotlib could volume render, I would not ask this question. But it can't. Can I however instead use Matplotlib to plot the surface of a 3D body? I.e. is there a way to (i) triangulate the outside or 'skin' of this body, and (ii) can that surface be plotted?

Note that the body is regularly sampled in the 3 dimensions (a simple 3D Numpy array) but can be quite complex, with indentations and bits sticking out, etc...

PS: I know there are other packages that can volume render, but I really would like to stick to Matplotlib.


问:
Verification on pressure predictor method for CFD code
I have developed a python code for a lid-drive cavity model. However, my results are not converging. The algorithm of my code looks like this:

Euler Momentum Equation looks like this:

$$\frac{u^{n+1}-u^{n}}{\Delta t}=-\frac{1}{\rho}\vec{P}^{n+1}-\vec{U}^{n}\cdot\nabla\vec{U}^{n}+\nu\nabla^{2}\vec{U}^{n}$$

1) The first step is to update time (current time = $\Delta t$ + current time) and enforce boundary conditions on $U$ and $V$.

2) Next step is to discretize the above formula to obtain predicted velocities:

$$\frac{u^{*}-u^{n}}{\Delta t}=-\vec{U}^{n}\cdot\nabla\vec{U}^{n}+\nu\nabla^{2}\vec{U}^{n}\quad\Rightarrow\quad u^{*}= u^{n}-\Delta t\left(\vec{U}^{n}\cdot\nabla\vec{U}^{n}-\nu\nabla^{2}\vec{U}^{n}\right)$$

3) The calculated  $U^*$ is used to find the pressure values $P^{n+1}$ using the Poisson's equation below:

$$\Delta P^{n+1}=-\frac{\rho}{\Delta t}\nabla \cdot U^{*}$$

4) Use $P^{n+1}$ values to find $U^{n+1}$ :

$$u^{n+1}=u^{n}+\Delta t\left(-\frac{1}{\rho}\vec{P}^{n+1}-\vec{U}^{n}\cdot\nabla\vec{U}^{n}+\nu\nabla^{2}\vec{U}^{n}\right)$$

5) Enforce boundary conditions, plot velocities and pressure and repeat 1-5 again until the end time.

I think I missed something that it does not show correct results. The source code is shown below:

```<code>import numpy as np
import math
import matplotlib.pyplot as plot
#--------------------------------------------------------------------------------------------------------------------------
#1.0 INPUTING MESH VALUES AND CREATING MESH GRID
Lx=1    # length at x direction
dx=0.05     # x step size
Ly=1    # length at y direction
dy=0.05     # y step size
i=math.ceil(Lx/dx) #1.3 Indexing ith cell faces
j=math.ceil(Ly/dy) #1.3 Indexing jth cell faces
#1.4 initializing velocity matrices
v=np.zeros((i+1,j+2))
u=np.zeros((i+2,j+1))
#--------------------------------------------------------------------------------------------------------------------------
#2.0 INPUT INITIAL VALUES
current_time=0
end_time=200
dt=0.025
#--------------------------------------------------------------------------------------------------------------------------
#3.0 INPUT CONSTANTS
mu=1/400 #Input kinematic viscosity
rho=1 # Input density
#--------------------------------------------------------------------------------------------------------------------------
# 4.0 FUNCTIONS
def b_conditions(uu,vv): # TO ESTABLISH BOUNDARY CONDITIONS
    vv[0]=0 #top
    vv[(len(vv)-1)]=0 #bot
    uu[:,0]=0 #left
    uu[:,(len(uu[0])-1)]=0 #right
    #3.2 Force boundary conditions at face
    u_top=1
    uu[0]=-uu[1]+2*u_top
    u_bot=0
    uu[(len(uu)-1)]=-uu[i]+2*u_bot
    v_left=0
    vv[:,0]=-vv[:,1]+2*v_left
    v_right=0
    vv[:,(len(vv[0])-1)]=-vv[:,j]+2*v_right
    return uu,vv
def velo_tem(u,v,j,i,dx,dy,dt,mu): # TEMPORAL DISCRETIZER FOR VELOCITIES
    vs=np.zeros((i+1,j+2))
    us=np.zeros((i+2,j+1))
    for y in range(1,j):
        for x in range(1,i+1):
            dudx2=(u[x][y-1]-2*u[x][y]+u[x][y+1])/(dx*dx)
            dudy2=(u[x+1][y]-2*u[x][y]+u[x-1][y])/(dy*dy)
            ududx=u[x][y]*(u[x][y+1]-u[x][y-1])/(2*dx)
            vdudy=0.25*(v[x][y]+v[x][y+1]+v[x-1][y]+v[x-1][y+1])*(u[x-1][y]-u[x+1][y])/(2*dy)
            us[x][y]=u[x][y]+dt*(mu*(dudx2+dudy2)-ududx-vdudy)
    for y in range(1,j+1):
        for x in range(1,i):
            dvdx2=(v[x][y-1]-2*v[x][y]+v[x][y+1])/(dx*dx)
            dvdy2=(v[x+1][y]-2*v[x][y]+v[x-1][y])/(dy*dy)
            udvdx=0.25*(u[x+1][y-1]+u[x][y-1]+u[x+1][y]+u[x][y])*(v[x][y+1]-v[x][y-1])/(2*dx)
            vdvdy=v[x][y]*(v[x-1][y]-v[x+1][y])/(2*dy)
            vs[x][y]=v[x][y]+dt*(mu*(dvdx2+dvdy2)-udvdx-vdvdy)
    b_conditions(us,vs)
    return us,vs
#--------------------------------------------------------------------------------------------------------------------------
#5.0 INITIATE LHS FOR LINEAR EQUATION (PRESSURE POISSON'S EQUATION)
#5.1 x direction
lhs_x=-2*np.eye(j,j)+np.eye(j,j,k=1)+np.eye(j,j,k=-1)
lhs_x[0,0]=lhs_x[j-1,j-1]=-1
lhs_x=np.kron(np.eye(i,i),lhs_x)/(dx*dx)
#5.2 y direction
lhs_y=-2*np.eye(i,i)+np.eye(i,i,k=1)+np.eye(i,i,k=-1)
lhs_y[0,0]=lhs_y[i-1,i-1]=-1
lhs_y=np.kron(np.eye(j,j),lhs_y)/(dy*dy)
#5.3 Initiate LHS
ln=j*i
lhs=np.zeros((ln,ln))
xx=np.arange(ln).reshape((i,j))
xx=np.transpose(xx).reshape(ln)
for x in range(ln):
    for y in range(ln):
        lhs[x,y]=lhs_x[x,y]+lhs_y[xx[x],xx[y]]
#5.4 i dont know the significance of this but was told to put value 1 at P(0,0)
lhs[ln-1]=0
lhs[ln-1,ln-1]=1
#--------------------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------------------
# CALCULATION STARTS HERE
#--------------------------------------------------------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------------------------------
# Initializing boundary conditionS
b_conditions(u,v)
# Conditions to stop when time reaches end_time
while current_time<=end_time:
#--------------------------------------------------------------------------------------------------------------------------
    #CFL control of dt
    print("current time: ",current_time)
    cfl=dt*(np.max(abs(u))/dx+np.max(abs(v))/dy)
    print("cfl: ",cfl)
    if cfl>=1:
        dt=cfl/(3*(np.max(abs(u))/dx+np.max(abs(v))/dy))
    #--------------------------------------------------------------------------------------------------------------------------
    #Time update
    current_time=current_time+dt
    #--------------------------------------------------------------------------------------------------------------------------
    #Initializing predictor velocities cells - step 2
    us,vs=velo_tem(u,v,j,i,dx,dy,dt,mu)
    #--------------------------------------------------------------------------------------------------------------------------
    #5.5 Creating Tridiagonal Matrix (RHS) for Poisson Equation
    n=0
    rhs=np.zeros(ln)
    for x in range(i):
        for y in range(j):
            rhs[n]=((us[x+1,y+1]-us[x+1,y])/dx)+((vs[x,y+1]-vs[x+1,y+1])/dy)
    rhs=-rho*rhs/dt
    #--------------------------------------------------------------------------------------------------------------------------
    #5.6 Solving pressure field
    p1=np.linalg.solve(lhs,rhs)
    p=np.reshape(p1,(i,j))
    #--------------------------------------------------------------------------------------------------------------------------
    #6.0 UPDATE VELOCITY FIELD
    vn=np.copy(v)
    un=np.copy(u)
    for y in range(1,j):
        for x in range(1,i+1):
            dudx2=(u[x][y-1]-2*u[x][y]+u[x][y+1])/(dx*dx)
            dudy2=(u[x+1][y]-2*u[x][y]+u[x-1][y])/(dy*dy)
            ududx=u[x][y]*(u[x][y+1]-u[x][y-1])/(2*dx)
            vdudy=0.25*(v[x][y]+v[x][y+1]+v[x-1][y]+v[x-1][y+1])*(u[x-1][y]-u[x+1][y])/(2*dy)
            un[x,y]=u[x,y]-dt*((p[x-1,y]-p[x-1,y-1])/(dx*rho)+ududx+vdudy-mu*(dudx2+dudy2))
    for y in range(1,j+1):
        for x in range(1,i):
            dvdx2=(v[x][y-1]-2*v[x][y]+v[x][y+1])/(dx*dx)
            dvdy2=(v[x+1][y]-2*v[x][y]+v[x-1][y])/(dy*dy)
            udvdx=0.25*(u[x+1][y-1]+u[x][y-1]+u[x+1][y]+u[x][y])*(v[x][y+1]-v[x][y-1])/(2*dx)
            vdvdy=v[x][y]*(v[x-1][y]-v[x+1][y])/(2*dy)
            vn[x,y]=v[x,y]-dt*((p[x-1,y-1]-p[x,y-1])/(dy*rho)+udvdx+vdvdy-mu*(dvdx2+dvdy2))
    u=un
    v=vn
    b_conditions(u,v)
    #--------------------------------------------------------------------------------------------------------------------------
    # VISUALIZATION & PLOTTING
    # Interpolation Function
    def inter(q11,q21):
        return (q11+q21)/2
    #// Generation of new loop
    #Interpolation and Coordinates of U
       #U Visualization
    u_i=np.zeros((i,j))
    for x in range(i):
        for y in range(j):
            u_i[x,y]=inter(u[x+1,y],u[x+1,y+1])
    #V Visualization
    v_i=np.zeros((i,j))
    for x in range(i):
        for y in range(j):
            v_i[x,y]=inter(v[x,y+1],v[x+1,y+1])
    #MeshGrid
    xx,yy=np.meshgrid(np.linspace(0.5,j-0.5,num=j),np.linspace(0.5,i-0.5,num=i))
    yy=np.flip(yy,0)
    fig, ax = plot.subplots()
    ax.contourf(xx,yy,p)
    ax.quiver(xx,yy,u_i,v_i)
    plot.show()
</code>```

Editted:

Boundary conditions are all zero except for u=1 at the top (v left=v right=v top=v bottom=0, u left=u right=u bottom=0). I use neuman boundary conditions for pressure so pressure derivative normal to grid walls = 0.

I've screenshot a plot. It does not show any influence from V, thus there's no swirling of the fluid in the cavity.


问:
Coupling 1D-3D problem in FENICS
I was wondering whether any of you knew if FENICS has the capability to solve coupled 1D-3D problems that are linear non-iteratively? As an example, pipes embedded in a porous domain will provide pressure sources (at their terminal ends) to the Darcy equation. This can be solved with a linear solver if I can concatenate the finite element matrix with the Poiseuille flow in the pipes. I can't quite see how to do this in FENICS, if it is at all indeed possible. Any suggestions, or directions towards more appropriate forums, would be massively appreciated.</p>

Title: Coupling 1D-3D problem in FENICS


问:
Why not use the preconditioned residual as termination criterion for preconditioned CG?
I have a Poisson equation with wildly varying material parameters (1 .. 1000), wildly varying element sizes (5 nm .. 100 um) and some quite anisotropic (tetrahedral) elements (100 nm x 100 um). I use (a C++ port of) Dan Spielman's approxChol as a preconditioner for a preconditioned CG method. It converges pretty fast (around 20 iterations), but the convergence test of the preconditioned CG method sometimes stops the iteration much too early, especially if I start the iteration with a "good" initial guess of the solution. The termination criterion is that the L2 norm of the residual is smaller than 1e-6 times the L2 norm of the right hand side. I looked around at some implementations of CG and some documents describing CG, and they all seem to use this as a termination criterion.

I wonder whether my trouble with the convergence test is caused by the bad condition of my matrix, and whether it would be a good idea to use the preconditioned residual instead as termination criterion. But why does nobody else does this? Perhaps because one would have to additionally compute $M^{-1}b$ in order to evaluate $b^TM^{-1}b$? Could one use $x_0^Tx_0 + r_0^TM^{-1}r_0$ instead ($r_0^TM^{-1}r_0$ gets computed anyway) to avoid that additional cost? (But the additional cost seems small anyway, so ...)


问:
How to translate Python infinitely big integer into Cython?
I want to factorize very huge numbers (e.g. 100-bits, 200-bits numbers) with Cython.

Hi everyone, I implemented the Elliptic Curve Method for factorization into Python 3.6. Now, I want to speed up my code using Cython (version 29.13). I'm a beginner in Cython world but I know that Cython works better if I define type for variables. So, I converted all my Python classes into Cython classes and now I would like to typize variables. I read that Cython automatically convert the "cdef int" declaration into classical Python integer with infinite length but this is not the case. When I try to factorize numbers like this "5192296858543544183479685583896053", I get an OverflowError because of "int is to big to convert into C long". There exists any ways to declare huge integer to speed up my code? The only variables without type declaration are the variables that could be very huge integer.

PS: I have already tried to use the cpython type uPY_LONG_LONG (unsigned long long) but it was useless because I always got the same error.

Thank you for paying attention to me and thanks in advance for your answers.

[UPDATE]

If I declare something like this:

```<code>cdef int function():
    cdef int a
    a = 2**100
    return a
</code>```

I get an OverflowError because of 2**100 is too huge to cast it into an integer.

If I import the long type from cpython, I get the same error:

```<code>from cpython import long as Long
cdef Long function():
      cdef Long a
      a = 2**100
      return a
</code>```

If I import the int type from cpython, I get no error but I have no speed up:

```<code>from cpython import int as Integer
cdef Integer function():
      cdef Integer a
      a = 2**100
      return a
</code>```

If I analyse the C++ code created as translation, I notice that the variable a has been declared as a pointer to a PyObject.
This is exactly the same translation that I get if I don't declare the variable. So, maybe in this context there are no difference.
I cannot improve all the for loops I use because of I have something like this:

```<code>for x in range(p):
     .....
</code>```

But if p is a huge integer and Cython declare p and x as pointers to PyObject, Cython can translate this loop into a C loop to speed up it.


问:
A maximization problem, with motivation in machine learning
Consider the minimization problem described this paper. Let $f_{\lambda}$ be the minimizer. As a part of extending my work, I am able to show the following facts

$$\lim_\limits{\lambda \to 0}\|f_{\lambda}\|_{L^2} = 0$$ and $$\lim_\limits{\lambda \to \infty}\|f_{\lambda}\|_{L^2} = 0$$

My problem now is (as I would like to extend my work), find $\lambda \in (0,\infty)$ for which $\|f_{\lambda}\|_{L^2}$ is maximum. Appreciate your suggestions to solve this problem.

The minimization problem from the linked paper is given below for the self containment of the post. If given that $k>\frac{m}{2}$, the paper proves that there is a unique minimizer for the functional $C(f)$ in the set $S$.

It is given that $k>\frac{m}{2}$

My progress

Partial Progress :

Progress : I am able to derive the corresponding PDE equations for the problem.

Let $f(\lambda,x) = f_{\lambda}(x)$. Then

The first equation corresponds to maximizing $\|f_{\lambda}\|$, while the second PDE is for the minimization problem associated with the parameter $\lambda$.

The second equation (minimization problem), given any $\lambda$, I can solve for $f(\lambda,.)$ either using linear algebra or steepest descent algorithm, which I have described in my article. Now I need to use this solution and the first equation to obtain $\lambda$, which is a problem I am facing.

Trying to solve using linear algebra, by formulating the discrete version of the problem using Fourier series coefficients and Plancheral theorem, I get stuck at the matrix problem described here.

More Partial Progress

An Iterative algorithm which is a modified steepest descent.

Initialize $f$.

Assuming some $\lambda$ and assuming gradient of $C_{\lambda}(f)$ wrt $f$ be $\nabla_f C_{\lambda}(f)$, and if we were to update $f$
with this gradient as in we do in steepest descent, it would be
$f^u_\lambda = f - \delta \nabla_f C_{\lambda}(f)$, where $\delta$ is
a constant learning rate. Now set
$\frac{\partial\|f^u_\lambda\|}{\partial \lambda} = 0$ and solve for
$\lambda$. Let the root be $\lambda_0$.

Update $f = f^u_{\lambda_0}$. (update $f$ as in steepest descent, but using $\lambda$ value as $\lambda_0$ which was computed in step
2.)

check some convergence criterion and if not met, go to step 2.

I have implemented this numerically and it converges as desired. Need to work on the proof.

PS : This was first posted on MO by me, 3 months back. Link and cross posted on math.SE here.


问:
Implementation of boundary conditions for 1D Euler equations
I'm trying to solve 1D Euler equations with gravity in spherical coordinates using a finite-difference TVD MacCormack method on a non-uniform grid of $N$ components, following the method provided in Appendix A in this paper. My aim is to evaluate stationary solutions.

The discretization for the source term in the second equation should guarantee that hydrostatic equilibrium, if reached, is maintained. I am using fixed BC at the lower boundary for all phisical variables and free-flow BC at the upper boundary; since the TVD correction requires 2 cells on both side of each cell, BC have to be imposed using two ghost cells:

\begin{equation}
\begin{cases}
\rho_1=\rho_2=1 \\\
u_1=u_2=0 \\\
p_1=p_2=1
\end{cases}
,\qquad
\begin{cases}
 \rho_{N}=\rho_{N-1}=\rho_{N-2} \\\
 u_{N}=u_{N-1}=u_{N-2} \\\
 p_{N}=p_{N-1}=p_{N-2}
\end{cases}
\end{equation}

I was testing the code by using as IC a particular equilibrium solution, i.e.:

\begin{equation}
\begin{cases}
\rho_{st}(r)=e^{-r}\\\
v_{st}(r)=0\\\
p_{st}(r)=(1+r)e^{-r}\\\
\end{cases},
\end{equation}

with $g=-\frac{\partial \phi}{\partial r}$, $\phi=\frac{r^2}{2}$ and $Q_\text{net}=0$. The code works well for all middle-grid points except near the boundary; in fact, if we consider the predictor step for the $(N-2)$-th component of the second equation at the first time step we have

\begin{equation}
\begin{aligned}
\big(\rho u & r^2\big)_{N-2}^{1/2} =\big(\rho u r^2\big)_{N-2}-\lambda \Big[\big(\rho u^2 r^2\big)_{N-1}-\big(\rho u^2 r^2\big)_{N-2}\Big]\\\
&-\frac{\lambda}{2}\Big[ \Big(\big(\rho r^2\big)_{N-1}+\big(\rho r^2\big)_{N-2}\Big)\Big(\phi _{N-1}-\phi _{N-2}\Big)+\big(r^2_{N-1}+r^2_{N-2}\big)\big(p_{N-1}-p_{N-2}\big)\Big]
\end{aligned}
\end{equation}

where $\lambda=\frac{\Delta t}{r_{N-1}-r_{N-2}}$.

In the first time step, since $v(r)=0$ and given the BC for pressure at the outer boundary, the equation reduces to

$$
\big(\rho u r^2\big)_{N-2}^{1/2}=\Big(\big(\rho r^2\big)_{N-1}+\big(\rho r^2\big)_{N-2}\Big)\Big(\phi _{N-1}-\phi _{N-2}\Big)
$$

This clearly shows that the equilibrium at this grid point is not preserved. The same thing happens near the lower boundary in the corrector step for index $3$. If I run the code for more time steps, the error propagates throughout the interior and the solution deviates from equilibrium overall points.

Is this due to my choice of BC? Or it's a more general problem of the reconstruction of MacCormack method? Do you think that using WENO schemes should be better to achieve better stationary solution (I'll have to modify the equations in the future like in the previously mentioned paper)?

If you have any reference for useful books or papers, it would be appreciated!


问:
Reducing run time of a numerical calculation using a mex file in Matlab
I wrote a Matlab code that involves doing a numeric calculation (relaxation), but it is quite slow. I learned of the possibility of using a mex file to run a C code and integrate it into Matlab, so I was thinking of doing the numerical calculation (which is relatively simple but involves loops and takes time) in C, and the rest (before and after) in Matlab.

The part of my Matlab code where the calculation is done:

```<code>    % evolution of the potentials %
    % note : for the index directions with periodic boundary conditions: index=mod(index-1,L)+1 . for index=index+1 it is mod(index,L)+1 , and for index=index-1 it is mod(index-2,L)+1 %
    for i_t=1:max_relaxation_iterations
        for q=1:length(i_eff_V_bounded) % this is set instead of running i=2:(L-1), j=1:L , k=1:L and ending up going over sites that are 0 in our effective system %
            i=i_eff_V_bounded(q);
            j=j_eff_V_bounded(q);
            k=k_eff_V_bounded(q);
            V0=V(i,j,k);
            V1=( V(i+1,j,k)+V(i-1,j,k)+V(i,mod(j,L)+1,k)+V(i,mod(j-2,L)+1,k)+V(i,j,mod(k,L)+1)+V(i,j,mod(k-2,L)+1) )/( system(i+1,j,k)+system(i-1,j,k)+system(i,mod(j,L)+1,k)+system(i,mod(j-2,L)+1,k)+system(i,j,mod(k,L)+1)+system(i,j,mod(k-2,L)+1) ); % evolving the potential as the average of its occupied neighbors %
            V(i,j,k)=V0+(V1-V0)*over_relaxation_factor; % evolving the potentials in time with the over relaxation factor %
            delta_V_rms(i_t)=delta_V_rms(i_t)+(V1-V0)^2; % for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg %
            delta_V_abs(i_t)=delta_V_abs(i_t)+abs(V1-V0); % for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg %
            delta_V_max(i_t)=max(abs(V1-V0),delta_V_max(i_t)); % for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg %
        end
    end
</code>```

So in C it should be something like:

```<code>#include <stdio.h>

int mod(int x,int N) /* a function for the modulo operator (instead of the remainder operator which is the % operator) assuming N is positive (x can be negative) */
{
return (x%N+N)%N;
}

double d_abs(double x) /* a function for the absolute value operator */
{
if x<0
    {
    return -x;
    }
else
    {
    return x;
    }
}

double max(double x,double y) /* a function for the max operator */
{
if x>y
    {
    return x;
    }
else
    {
    return y;
    }
}

/* evolution of the potentials */
/* note : periodic boundary conditions for the j,k directions */
void potentials_evolution(int max_relax_iters,int N_eff_occ_sites,int i_eff_V_bounded[],int j_eff_V_bounded[],int k_eff_V_bounded[],int system[][][],over_relax_fact,double V[][][],double delta_V_rms[],double delta_V_abs[],double delta_V_max[])
{
int i_t,q,i,j,k;
double V0,V1;
for(i_t=0;i_t<max_relax_iters;i_t++)
    {
    for(q=0;q<N_eff_occ_sites;q++) /* going over only the occupied sites left in our effective system */
        {
        i=i_eff_V_bounded[q];
        j=j_eff_V_bounded[q];
        k=k_eff_V_bounded[q];
        V0=V[i][j][k];
        V1=( V[i+1][j][k]+V[i-1][j][k]+V[i][mod(j+1,L)][k]+V[i][mod(j-1,L)][k]+V[i][j][mod(k+1,L)]+V[i][j][mod(k-1,L)] )/( system[i+1][j][k]+system[i-1][j][k]+system[i][mod(j+1,L)][k]+system[i][mod(j-1,L)][k]+system[i][j][mod(k+1,L)]+system[i][j][mod(k-1,L)] ) /* evolving the potential as the average of its occupied neighbors */
        V[i][j][k]=V0+(V1-V0)*over_relax_fact; /* evolving the potentials in time with the over relaxation factor */
        delta_V_rms[i_t]=delta_V_rms[i_t]+(V1-V0)*(V1-V0); /* for each t at a given p, we sum over (V1-V0)^2 in order to eventually calculate delta_V_rms_avg */
        delta_V_abs[i_t]=delta_V_abs[i_t]+d_abs(V1-V0); /* for each t at a given p, we sum over |V1-V0| in order to eventually calculate delta_V_abs_avg */
        delta_V_max[i_t]=max(d_abs(V1-V0),delta_V_max[i_t]); /* for each t at a given p, we take the max of |V1-V0| from all the sites in order to eventually calculate delta_V_max_avg */
        }
    }
}
</code>```

And so in Matlab I will replace the part of my Matlab code shown above with something like:

```<code>potentials_evolution(max_relax_iters,N_eff_occ_sites,i_eff_V_bounded,j_eff_V_bounded,k_eff_V_bounded,system,over_relax_fact,V,delta_V_rms,delta_V_abs,delta_V_max);
</code>```

How do I implement this? I tried looking for a simple way to do it but I couldn't figure out how to properly do it.

Note 1: This numeric calculation is done not just once but many times for different systems that are generated randomly (there is a for loop going over the different systems).

Note 2: My C is quite rusty.


问:
What is the best methodology for physics simulators of large floating base rigid body systems?
I want to implement a physics simulator for large floating base rigid body systems from scratch. The Rigid Body Dynamics Systems (RBD) should typically have the following characteristics:

About ~50 degrees of freedom (DoF) plus the 6 DoF of the floating base (e.g., a legged robot, a humanoid robot, etc.).

It will always have tree structure (no parallel robots)

The simulation needs to be fast.

There will be contact forces (e.g., during walking, grasping or manipulation)

What would be the best physics formulation in terms of numerical accuracy and computational complexity? Lagrange appears too slow, but Newton Euler is often hard to comprehend.


问:
How to solve advection equation using semi-lagrangian method?
I am working on something that involves solving an advection equation $\partial{x}/\partial{t}+\vec{u}\cdot\nabla{x}=0$ in 3D. I discretized the space into 3d cartesian grid and used the Semi-Lagrangian method to solve the equation.

By 'Semi-Lagrangian method' I mean the method that Jos Stam proposed in his paper 'Stable Fluids'(People in the computer graphics community might be familiar with him). The method states that to get the new value of $x$ at current position, we trace back in time and use the value of $x$ at the 'traced-back' position. It is based on 'the method of characteristic'.

Imagine I am advecting some objects with a velocity field and $x$ is some attribute of the objects. I update $x$ according to the advection equation.

I have one question: What if I get the value of $x$ outside the boundary of the object while tracing back? I can't use this value to update the value of $x$ inside the object since it's undefined. Is there any method to solve this problem?


问:
Pseudospectral method for Rayleigh-Benard with constant temperature gradient
$$
\nabla\cdot \mathbf{u} = 0 \\
\frac{\partial \mathbf{u}}{\partial t}+\left(\mathbf{u}\cdot \nabla\right)\mathbf{u} = -\nabla p+\nu\nabla^2\mathbf{u}+\alpha g\theta\mathbf{e}_z\\
\frac{\partial\theta}{\partial t}+ \left(\mathbf{u}\cdot \nabla\right)\theta = \kappa\nabla^2\theta+\mathbf{u}_z\frac{\Delta T}{H}
$$

Can I apply Fourier-spectral method for 3D-Rayleigh-Benard convection with a constant temperature gradient in the vertical direction?

I want to take all the directions as periodic.

There are papers which have done this with free-slip boundary condition in the vertical direction using the Fourier method.

But there is a paper that has done it use Lattice-Boltzmann method with triply-periodic boundaries.

Edit:
$\theta$ here is the temperature fluctuation and not the actual temperature.
Knowing that the equations can be solved using periodic boundary conditions.

I'm not sure if the question would be helpful to anyone (since it leaves out some crucial detail).


问:
Heisenberg Model python : Specific heat capacity for spin 2
I have the correct plot for specific heat capacity when I am using the formula which is $C_V$ = differentiation of entropy with respect to temperature. However, When I try to calculate $C_V$, by using formula

$$C_V =  \frac{\langle E^2\rangle - \langle E\rangle ^2}{T^2}\, ,$$

I am getting a strange result.

I am expecting the second figure just like the first one. Specific heat doesn't need to go negative however it has many negative values in the second figure.

Following is the code, any help would be immensely welcome.

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

np.random.seed(789645)
# CONSTANTS INSTANTIATIONS          ----------------------------------------------------
J =1.0  ;Jeff = 2*J;D = 0.*Jeff;h = 0;Spin = 2
n = int(2*Spin + 1)  #size of the matrix
k = Spin*(Spin+1)   # total value spin can take
T1 = 0.1; T2 = 5;S = np.arange(-Spin,Spin+1e-5)   #possible values of S
mm=30;T = np.linspace(0.1,6.,mm);beta = [1./a for a in T]
#---------------------------------------------------
Su = np.zeros((n,n));Sd = np.zeros((n,n));Sz = np.zeros((n,n))

np.fill_diagonal(Sz,S)   # this gives Sz which is diagonal matrix. we
#Assume   Sz to be the same for both the lattice A and B
for i in range(n-1):
    Su[i,i+1] = np.sqrt(k - (S[i]+1)*S[i])
    Sd[i+1,i] = np.sqrt(k - (S[i]+1)*S[i])

Sx=0.5*(Su+Sd);Sy=(-0.5j)*(Su-Sd)

def ham(mA,mB,site):
    if site == "A":
            Ham =-D*Sz*Sz-h*Sz+Jeff*(mB['x']*Sx+mB['z']*Sz)- \
              Jeff/2.0*(mA['x']*mB['x'] + mA['z']*mB['z'])*np.eye(n)
    else:
            Ham =-D*Sz*Sz-h*Sz+Jeff*(mA['x']*Sx+mA['z']*Sz)-\
              Jeff/2.0* (mA['x']*mB['x'] + mA['z']*mB['z'])*np.eye(n)
    return np.linalg.eigh(Ham)

def expectation(b,en,ev,spin):
    Z=np.sum(np.exp(-b*en))
    ave = 0.
    for i in range(n):
        ave += np.exp(-b*en[i])*(np.matmul(np.asmatrix(np.reshape \
               (ev[:,i],(n,1))).getH(),np.matmul(np.asmatrix(spin), \
               np.asmatrix(np.reshape(ev[:,i],(n,1))))))
    return ave[0,0]/Z

def FreeEnergy(b,en):
    Z=np.sum(np.exp(-b*en))
    FE=-1/b*np.log(Z)
    return (FE,Z)

nn = 1   # samplin number
ff = []  # Free energy save here
mA_list = []  # expectation of <S>A
mB_list = []  # expectation of <S>B
MxA = [];MzA = [];MxB = [];MzB = [];Mx = [];Mz = [];ent = [];Cv = []
mA = {'x':np.random.rand(),'z':np.random.rand()}
mB = {'x':np.random.rand(),'z':np.random.rand()}

for bt in beta:
    change=1
    itr=0
    while change > 1e-5:

            itr += 1
            mA_new={}
            mB_new={}

            enA,evA = ham(mA,mB,"A")
            mA_new['x']=expectation(bt,enA,evA,Sx)
            mA_new['z']=expectation(bt,enA,evA,Sz)

            enB,evB = ham(mA,mB,"B")
            mB_new['x']=expectation(bt,enB,evB,Sx)
            mB_new['z']=expectation(bt,enB,evB,Sz)
            change=abs(mA_new['x']-mA['x'])+\
                   abs(mA_new['z']-mA['z'])+\
                   abs(mB_new['x']-mB['x'])+\
                    abs(mB_new['z']-mB['z'])
            p=0.2
            mA['x']=p*mA['x']+(1-p)*mA_new['x']
            mA['z']=p*mA['z']+(1-p)*mA_new['z']
            mB['x']=p*mB['x']+(1-p)*mB_new['x']
            mB['z']=p*mB['z']+(1-p)*mB_new['z']

    mA_list.append((mA['x'],mA['z']))
    mB_list.append((mB['x'],mB['z']))
    (FreeA,ZA)=FreeEnergy(bt,enA)
    (FreeB,ZB)=FreeEnergy(bt,enB)

    # These are the average energy  <E>
    UA = np.sum(np.dot(enA,np.exp(-bt*enA)))/ZA
    UB = np.sum(np.dot(enB,np.exp(-bt*enB)))/ZB
    UT = UA + UB

    ent.append((UT-FreeA-FreeB)*bt)   # This gives the entropy, F = U - TS

    #  These are <E^2>, I think the problem lies here.
    UA2 = np.sum(np.dot(enA**2,np.exp(-bt*enA)))/ZA
    UB2 = np.sum(np.dot(enB**2,np.exp(-bt*enB)))/ZB
    UT2 = UA2 + UB2
    Cv.append(((UT2 - UT**2))*bt**2)

sp_heat = np. diff(ent, n =1)
deltaT = (T2-T1)/mm
sp_heat[:] = T[:mm-1] * [a/deltaT for a in sp_heat]

# The first plot is correct. Which is obtained from the derivative of
# Entropy, where the second plot is incorrect
# and should appear as first is obtained from the fluctuation of energy
#<E^2> -<E>^2)*bt**2

#plt.plot(T[:mm-1], sp_heat)

plt.plot(T, Cv)
</code>```


问:
Classification of multiobjective optimization algorithms
I am looking for a good (canonical?) overview paper(s)/book(s) on the classification of multiobjective optimization algorithms. I am focused on obtaining a representative set of Pareto optimal solutions generally for real-world optimization problems in electromagnetics; however, the reference should not be field-restrictive.

The reference ideally should not focus only on evolutionary methods (this is the case of overview papers that I encountered), but probably compare the method families and give a broader overview.

Particular goals:

have it as a to-go reference for the topic

get the overview of major families and subfamilies of multiobjective optimization algorithms. In particular, Wikipedia has

Most a posteriori methods fall into either one of the following two classes: mathematical programming-based a posteriori methods, where an algorithm is repeated and each run of the algorithm produces one Pareto optimal solution, and evolutionary algorithms where one run of the algorithm produces a set of Pareto optimal solutions.

I wonder, what would be the methods (hybrid?) that does not fall in either of those classes.

get it as a reference to use as a "justification" for choosing or not choosing evolutionary-based methods


问:
How to deal with a huge system of ODEs in Boost ODEINT?
I am using the C++ library ODEint, which is part of Boost, to solve an extremely large system of coupled ODEs - in particular 1975 equations with large rational functions in the coefficients. In the provided examples by the authors the ODE system is always explicitly written in the source code. Unfortunately my source code file is 37MB large, if I do it in the same manner. The compilation process is now running for 14 hours and has still not finished. </p>

<p>Is there another approach to cope with this huge system? Or is it possible to load the ODEs during runtime? </p>


问:
Need suggestions about technical difficulties - Drying + Pyrolysis of coal particle
In OpenFOAM by default, the FireFOAM is well supported for solid pyrolysis modeling. With that in mind, I managed to built my solver for a modified version of pyrolysis (for dry coal - without moisture) based on FireFOAM. Now in the case of wet coal pyrolysis, the pyrolysis model should be developed [absence of air - no combustion] involving the influence of moisture drying.

--> Drying and Pyrolysis occur simultaneously in different regions of the same particle (however both phenomenon should be modeled altogether or modeled separately) in OpenFOAM?

(I) Basically, Pyrolysis is modeled based on solidChemistry ==> Chemistry model for solid thermodynamics, solved using Arrhenius type equation.
~/OpenFOAM/OpenFOAM-2.1.1/cokeovenGasFOAM/regionModels/pyrolysisModels/reactingOneDim

PYROLYSIS

solve (fvm::ddt(rho_) == - solidChemistry_->RRg() ); //continuity equ

fvScalarMatrix YiEqn ( fvm::ddt(rho_, Yi) ==> solidChemistry_->RRs(i)
); //species equ.

(II) In my case, I need to model the continuity equation of drying as:

DRYING

solve ( fvm::div(rhoU_, Yi) == - XXXX_->RRg() ); //continuity equ

Here, this XXXX cant be used as solidChemistry technically, because the drying is not solved using Arrhenius type equation.

How to proceed with this problem further,

(i) Should I need to couple drying model externally ==> such that it should satisfy the real flow physics (Pyrolysis and Drying occur simultaneously) ?

(ii) Should I need to create sub-models based on DRYING?

I would really appreciate if anyone could give me a hint where to look at or point me about the implementation in OpenFOAM.

Thank you!!


问:
3D log density plot in ParaView
I have a .csv file of thousands of (x,y,z) point particles that I would like to visualize in ParaView. I am able to plot a 3D scatterplot using the TableToPoints, and by decreasing point opacity I can see some substructure in the particles' arrangement.

I'm looking for a better option to easily see different regions where particles are clustered or sparse. Is there a 3D density plotting option in ParaView (preferably using a log scale for density). i.e. a 3D histogram where I can see regions of high and low particle density?

Title: 3D log density plot in ParaView


问:
How to compute the following Crank-Nicolson scheme for the viscous Burgers' Equation?
I am attempting to replicate results from this article.

I'm not sure why but my results are completely different and wrong. For example, the exact solution with parameters ($x=0.1$, $T=0.01$, $Re=0.1$ and $\Delta t=0.0001$) is $0.11461$ but my approximated solution gives $-175.4007$.

Is anyone able to spot my mistake? I've included my code (in R) below.

The Thomas algorithm was implemented.

I've asked a very similar question in the Math SE website recently but I think with the focus on the coding as well, this perhaps might be a better place to post this question.

In case you can't access the article, these are the relevant information:

$0<l<N$ represents spatial point, $0<j<M$ represents time, $k=T/M$, $h=1/N$ and $x_i=ih$.

The initial condition for $0<x<1$ is given as
$$\phi(x,0)=\exp\left({-\frac{Re}{2\pi}\left[1-\cos(\pi x)\right]}\right)$$

and boundary conditions for $0\leq t \leq T$ are
$$\phi_x(0,t)=\phi_x(1,t)=0.$$

Then the Crank-Nicolson finite difference discretization equations are given as:

For $l=0$
$$-r\phi_{l+1,j+1}+(1+r)\phi_{l,j+1}+0=r\phi_{l+1,j}+(1-r)\phi_{l,j}$$

For $l=1:N-1$
$$-\frac{r}{2}\phi_{l+1,j+1}+(1+r)\phi_{l,j+1}-\frac{r}{2}\phi_{l-1,j+1}=\frac{r}{2}\phi_{l+1,j}+(1-r)\phi_{l,j}+\frac{r}{2}\phi_{l-1,j}$$

For $l=N$
$$0+(1+r)\phi_{l,j+1}-r\phi_{l-1,j+1}=r\phi_{l-1,j}+(1-r)\phi_{l,j}$$

where $r=k/Reh^2$.

The solution via Hopf-Cole transformation is then calculated as:

$$u_{l,j}=-\left(\frac{1}{Re}\right){\frac{\phi_{l+1,j}-\phi_{l-1,j}}{h\phi_{l,j}}}$$

```<code># Parameters
Re=0.1
T=0.01
input_dt=0.0001
N=80

M=T/input_dt
h=1/N
k=T/M
r=k/(Re*(h^2))

# Spatial discretization
x_discrete=numeric()
x_discrete[1]=0
for( i in 1:(N+1)){
    x_discrete[i] = (i-1)*h
}
x_discrete[N+1]=1

# Calculate initial conditions
u_init=exp(-(Re/(2*pi))*(1-cos(pi*x_discrete)))
# Adjust u_init with BC
u_init[1]=0
u_init[N+1]=0

# matrix to store solutions
u=matrix(0, (N+1), (M+1))
# Add in initial conditions
u[,1]=u_init

# Crank-Nicolson
# For each time-step
for(j in 1:(M)){
    # A,B,C, D set up from:
    # https://www.sfu.ca/~rjones/bus864/notes/notes2.pdf
    A=numeric()
    B=numeric()
    C=numeric()
    D=numeric()

    # For each spatial node (x values)
    for(l in 1:(N+1)){
        if(l==1){
            # LHS
            A[l]=-r
            B[l]=(1+r)
            C[l]=0
            #RHS
            D[l]=r*u[l+1,j]+(1-r)*u[l,j]
        }else if(l==(N+1)){
            # LHS
            A[l]=0
            B[l]=(1+r)
            C[l]=-r
            #RHS
            D[l]=r*u[l-1,j]+(1-r)*u[l,j]
        }else{
            # LHS
            A[l]=-r/2
            B[l]=(1+r)
            C[l]=-r/2
            #RHS
            D[l]=(r/2)*u[l+1,j]+(1-r)*u[l,j]+(r/2)*u[l-1,j]
        }
    }

    # Set up tridiagonal-matrix
    # Method taken from:
    # https://stackoverflow.com/questions/49536400/create-a-tridiagonal-matrix
    m = diag(B)
    m[row(m) - col(m) == 1]  = A[-1]
    m[row(m) - col(m) == -1] = C[-length(C)]

    # Solve for solutions
    sol=base::solve(m)%*%t(t(D))

    # Update solution matrix
    u[ ,(j+1)]=sol

    u[1, (j+1)]=0       # Adjust according to BC
    u[(N+1), (j+1)]=0       # Adjust according to BC
}

# Exact solution at x=0.1 should be 0.11461 (from Table 1)
# BUT mine gives -175 (from the line of code below)
-(1/Re)*(u[9+1,101] - u[9-1,101])/(h*u[9,101])
</code>```


问:
How one can simulate a system given by differential equation?
I want to simulate a diffusion environment given by the differential equation
$$\frac{\partial u(x,y,t)}{\partial t}=D\left(\frac{\partial^2 u(x,y,t)}{\partial x^2}+\frac{\partial^2 u(x,y,t)}{\partial y^2}\right),$$
in a constrained environment with reflecting and absorbing boundaries and certain initial conditions. I know it may be done with the help of particle simulation method. But I want to know is there a way, in general, to obtain numerically the concentration profile $u(x,y,t)$   by using directly the differential equations?

I am actually trying to verify some of the ambiguities. I have used the particle simulation method to verify obtained PDF (normalized concentration profile for my case) but it gives me different results and I am not sure whether particle simulation method is correct to validate the obtained result?


问:
Evaluate Nth root of a rational to a correctly rounded float
Excuse my lack of vocabulary for I have no formal training in this field, which is also why I ask this question - it may be trivial or it may be impossible.

I want to evaluate an expression in the following form down to a binary float with "correct" rounding and I'm worried that double-rounding will corrupt my result:

$$\sqrt[c]{\frac{a}{b}} \qquad a \in \mathbb N \quad b,c \in \mathbb N_{> 0}$$

By binary float I mostly mean a plain old double, but since I use MPFR it could have any finite precision.

By correct rounding I mean as if the complete expression is evaluated to infinite precision and then rounded correctly, without any chance of "double rounding" corrupting the result.

I use MPFR and GMP to work with these expressions using arbitrary sized integers and floats, and so far I divide the work into two steps, rounding to a larger float in between:

$$t = a/b \qquad t:\mathrm{e.g.\ binary128}$$
$$y = \sqrt[c]{t} \qquad y: \mathrm{e.g.\ binary64}$$

Using MPFR I can calculate both $a/b$ and $\sqrt[c] t$ correctly rounded to any arbitrary precision, but I don't know how to do these two steps in a way that guarantees correct rounding for the complete expression!

Here is a code example of what I'm currently doing. Since it's written in C it is a bit clumsy (no overloads or default arguments) but I hope you can see the steps even if you're not fluent in this language or the libraries I use:

extern unsigned long int a, b, c;     // b and c != 0
MPFR_DECL_INIT( y, DBL_MANT_DIG    ); // Fits a double exactly
MPFR_DECL_INIT( t, DBL_MANT_DIG*2  ); // Twice the precision I want
mpfr_set_ui   ( t, a,    MPFR_RNDN ); // Exact conversion
mpfr_div_ui   ( t, t, b, MPFR_RNDN ); // Correctly rounded
mpfr_rootn_ui ( y, t, c, MPFR_RNDN ); // N'th root correctly rounded
```

I'm expecting one of these to hold true:

Twice the precision for the intermediate is not enough

Twice the precision for the intermediate is more than needed - I only need N extra bits

I need to do it differently but it's still possible in a finite time and space

It's not possible to guarantee correct rounding because of reason


问:
Fast convergence of smoothing of periodic noise
I have essentially periodic data from a simulation (not exactly periodic but is qualitatively fairly periodic), and I'd like to take an average or noise filter of some sort that I can get a well converged average/smoothed value. Ideally I'd like something that converges to a steady value as the square of the number of samples, but am unsure if this is possible. I'd like a filter that for approximately 10000 samples would be almost invariant with any further added points and the filtered value would be converged to approximately machine zero. Is there such a filter? I'd appreciate any help for either specific methods or families of methods I can look at.

Thanks


问:
Error on the fit parameters when several good fits exist
I am using the reduced chi-squared statistic to determine the goodness of fit. I run several simulations and determine that a parameter 'p' has a certain range of values that all give values between 0.9-1.1 for the reduced chi-squared. Visually all the fits are deemed acceptable. Is there are a way to estimate the error on the parameter p and is there any literature that deals with the acceptable chi-squared values?

Suggestions along the lines of MCMC are welcomed. Importantly, any frequentist approaches to this problem are also welcomed


问:
Are there well-known methods for navigating on kd-trees?
When you have a mesh, there are many well-known methods to navigate it, as for example using a half-edge data structure, that allows easy circulation around faces and vertices.
Are there similar structures/algorithms for kd-trees, that allow you to easily find adjacent and overlapping (half)faces?

I am mostly interested in circulating around an edge to find all adajacent (half)faces and cells, what seems to get rather complicated especially on imbalanced kd-trees. In addition, identifying overlapping halffaces which neighbor faces are overlapping the same opposite halfface would be useful.

I wonder if there are some standard algorithms. It is hard to search for them, because every book and every internet search result is full of standard algorithms for searching nearest neighbors in euclidean space using a kd-tree, without mentioning how to navigate in the kd-tree mesh or graph for queries with respect to the tree itself.

Something that works analogous to a halfedge structure would be really nice, but I suppose most algorithms will rather traverse the tree.

If there are no specific algorithms, I would really like suggestions how to work with a collection of cuboids, i.e. the geometric representation of the kd-tree bounding boxes, for querying adjacency, overlaps and similar geometric properties.

Halfedges

For reference short introduction to halfedge structures on meshes. A halfedge structure looks like this:

(Image taken from the CGAL manual)

It makes all important queries really easy:

Circulate around a face: Follow the next pointer until you're back at the start edge.

Find a neighbor face: Follow the opposite pointer to the halfedge of the other face.

Find the edges around a vertex: Follow the next pointer, then the opposite pointer to find the next incident edge.

The structure has many more useful properties, like allowing to follow the boundary of a mesh using the next points of boundary edges and circulating around holes, in the same way as circulating inside a face.

KD-Trees

These operations are only well defined on polygonal meshes. A kd-tree looks like this:

(Picture taken from the Wikipedia article on kd-trees)

The tree is build by subdividing an initial box. You start in one dimension and when you subdivide a box, you always use the next dimension (e.g. x, y, x, y for 2D or x, y, z, x, y, z for 3D). There are several ways to choose the splitting plane, depending on your actual application.

Such a tree is no mesh, as it has T-junctions by design. If you try to build a halfedge structure in this example tree, you have edges, that have multiple opposite edges, as you see in the rightmost column of the tree.

In addition, such a tree has overlapping (half)edges like the blue edge on the left side of the example tree. When you try to split it into half edges, you get four halfedges, each adjacent to two opposite edges.

Adjacency in a KD-Tree

If you want to use such a tree for example to create a mesh of the points contained in the tree, you need adjacency information between the boxes.

Here is an example in which I used a kd-tree to create a mixed quad/triangle mesh using a kd-tree that separates the vertices:

(Note that not all boxes with overlapping edges are connected, see below for more on how this graph is created)

Neighbors sharing an edge are neighbors in the corresponding graph as well (and the construction guarantees that it is the same edge). But boxes that do not share the same parent node, may be adjacent to each other. In the upper left you see a triangle, that corresponds to the subtree:

```<code>|\
A |\
  B C
</code>```

in which the neighboring leaf nodes are connected with the neighbor cell of their parent cell. This is still a rather easy query and the intersections are uncomplicated as well. There are no overlapping edges and an edge has either one or two non-overlapping opposite edges.

But you also see many triangle and quads, that connected boxes that are far apart from each other in the tree with no obvious way to find the needed adjacency information.

Finding neighbors in the tree

I created this graph using the following algorithm:

```<code>start at all* splits:
  traverse the two subtrees in parallel as follows:
    if the dimension is the same as the split:
      follow+ the subtrees in the split direction in both trees (e.g. right, left)
    else:
      follow+ both trees in the first direction (e.g. left, left) AND
      follow+ both trees in the second direction (e.g right, right)
    if both new nodes are leaf nodes, connect them

* all O(n) inner tree nodes.
+ When you arrive at a leaf node, store it and continue traversing the other subtree as before.
</code>```

This way you stick to the splitting plane and traverse it always in the same directions for both subtrees, without taking the actual coordinates into account to find boxes to connect. These edges result in the quads and triangles in the example above.

But this approach does not even find all adjacencies. For example most of the quads should be split into triangles, because there are additional adjacencies, that are not found by the algorithm.

Even the second from left quad  at the bottom can be split into two triangles, with regard to adjacency, its edge intersection is just very short and thus harder to see than the other intersections.

(Finding these triangles is only one application, on the other hand it could be quite useful to identify what connections could create a quad and still produce a useful mesh. The mesh in the example may be good for some applications and bad for others)

Are there techniques for more efficient and more general queries?

The combinatoric approach does not guarantee that the neighbors actually share a non-zero edge intersection and as you see in the example, the edges connecting their centers (or other points inside the cells) may cross other boxes in the tree, what cannot be detected by this algorithm.

In addition, the algorithm starts at $n$ splitting planes and traverses the full subtrees in time $n log n$, what seems rather inefficient.

I am now looking for techniques for faster queries to find neighbors of cells, edges, and vertices. I would for example like to start at an edge of a box in the tree and find all neighbor faces (and thus neighbor boxes). This would not only involve faces that are adjacent to the edge, but also neighbors of the parent cell and more complex neighborhoods that may occur in a kd-tree.

Papers I found so far

Rope trees, that construct additional links inside the tree for more efficient ray intersections:

Havran, V., Bittner, J., & Zára, J. (1998, April). Ray tracing with rope trees. In 14th Spring Conference on Computer Graphics (pp. 130-140). (PDF)

Using nested kd-trees allows for some of the queries:

Greß, A., & Klein, R. (2004). Efficient representation and extraction of 2-manifold isosurfaces using kd-trees. Graphical Models, 66(6), 370-397. (Website)


问:
Implementation of Lanczos method that returns tridiagonal matrix
The Lanczos method can be used to obtain extremal eigenpairs of sparse symmetric or hermitian matrices. I know there are several implementations of the Lanczos method (as well as Arnoldi, Davidson, and others) out there, but I am specifically interested in one that returns the tridiagonal matrix, not just the eigenvalues and eigenvectors. The way I am doing this now is using my own Lanczos implementation in Fortran. Is there a maintained library that does the same efficiently, preferably in python?


问:
Grids for atmosphere simulation with finite volumes on the globe
I am currently in the early construction process of building a simple CFD model of a rotating planetary atmosphere. The planet should be allowed to tilt significantly, so that a time-dependent source function for heating the atmosphere can also go overhead the poles.

Because of this, I have to consider grids other than the popular longitude-latitude grid in atmospheric (earth) science. My first instinctive move was to look for a tesselation of the sphere which is more agnostic towards the planetary rotation.

I found a classic Sadourny paper that also provides a simple algorithm how to generate said tesselation using triangles. I implemented this and it works nicely to generate grids of arbitrary refinement.

However, it turns out (Staniforth, grid review) that triangles are bad and one should use quadrilaterals instead. A quadrilateral tesselation based on a icosahedron as lowest resolution grid is still possible.

Now the latter review I've mentioned did comment on grid problems using finite difference methods, which i do not intend to use. My first question would now be:

Does a triangular tesselation of spheres with finite volume methods have the same problems as finite difference methods with noise and accuracy? Also I imagine, computing a FV method like the WAF (weighted averaged flux, mentioned in the Toro book, chapt. 16 to be easily implementable for complex geometries??) on triangles to be awkward, because of the geometrical enlargment of a wave into the next triangle.

and because of those considerations (on the risk that that's probably a more opinion-based question)

Would numerics on a sphere simplify if I just abandon the triangular tesselation and go for a quadrilateral one, by pairwise joining triangles? Is there maybe a similar review like the Staniforth (see link above) just for finite volume methods?


问:
Second fundamental form - Maple
I would like to know the command/line in Maple 16 or similar to obtain the second fundamental form tensor for a given metric. I've managed to obtain Rienmann and Ricci tensor, even Weyl, but I can't find anything on the internet that would help me to obtain the K_mu_nu of a generic spherically symmetric metric such as:

```<code>gmetric := evalDG(A(r)*`&tensor`(dt, dt)+B(r)*`&tensor`(dr, dr)+r^2*`&tensor`(dtheta, dtheta)+r^2*sin(theta)^2*`&tensor`(dphi, dphi))
</code>```

Any suggestion?

Thank you in advance

OwnerDisplayName: Mitsumako

Title: Second fundamental form - Maple