问:
How group similar data in a single classes and reduce the error
I have vehicles gps information in my real time traffic application. Averaging that information I know the speed of every road at any time. The problem is that is too much data to send back to the devices, so I need a way to group that speed data into classes so I can to simplify the problem.

Right now I have something like this:

```<code> road_id   time   speed
    1      09:00   45
    1      10:00   47
    1      11:00   38
    2      09:00   78
    2      10:00   77
    2      11:00   61
    3      09:00   51
    3      10:00   48
    3      11:00   41
</code>```

As you can see:

when close to 12pm more cars are in the roads so the speed decrease.

road 1 and road 3 has a similar behavior probably both are regular street.

road 2 is faster, maybe a highway.

So instead of save the speed for each road for each hour. I just want define general road behavior for the day and say that road is that type.

So I will have a table for roads and another for road_type

```<code>road_id      name        type_id
  1        A street         1
  2        B Highway        2
  3        C street         1

type_id  spd_09-10  spd_10-11 spd_11-12 ....
  1         50         50        40
  2         80         80        60
</code>```

So I need an algorithm to create those types efficiently and reducing the error with the original data. I understand if I increase the number of types the error will decrease.

My guess is there is something already related to this problem just don't know how is called, and I don't want invent the wheel again.


问:
How to stochastically estimate the trace of a matrix?
Specifically, the diagonal elements (can possibly both positive and negative) of the matrix can be computed efficiently but the total number is large ($\mathcal O(10^{18})$).

My first thought about this was running an MCMC (Markov Chain Monte-Carlo) with the weight set to be the absolute value of each element; however, this just exactly leads to the sign problem. But I'm still wondering if there is any smart way to tackle that.

ContentLicense: CC BY-SA 3.0


问:
How accurate is cumtrapz in MatLab?
Say I have a set of discrete acceleration data and want to integrate it to get a set of velocity data.  How accurate is the cumtrapz (Cumulative trapezoidal numerical integration) integrator in MATLAB?

Does it follow that the more acceleration data I have (for whatever fixed timescale), the better results I'll get from cumptrapz?

Title: How accurate is cumtrapz in MatLab?


问:
Accuracy of finite difference method for heat equation on a disk
To study an approximation for the heat equation $$\frac{\partial^2 u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\frac{\partial^2 u}{\partial\theta^2}=f(r,\theta)$$
on the disk $D=(0,1)\times(0,2\pi)$ with periodic boundary conditions, we used the following finite difference method
$$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta r)^2}+\frac{1}{r_i}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta r}+\frac{1}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta\theta)^2}$$
where $$r_i=(i-\frac{1}{2})\Delta r, \space\space\space\space i=1,2...N+1$$ $$\theta_j=(j-1)\Delta \theta, \space\space\space\space j=1...M+1$$
$$\Delta r=\frac{2}{2N+1}, \space\space\space\space \Delta \theta=\frac{2\pi}{M}$$
$$M=2N$$
By Taylor expansion, we computed the LTE which should be $$\tau=(\Delta r)^2(\frac{1}{12}u^{IV}(r,\theta)+\frac{1}{6}u^{III}(r,\theta))+\frac{(\Delta \theta)^2}{12}u^{IV}(r,\theta)+O((\Delta r)^4,(\Delta \theta)^4)$$
Thus the method should show second order convergence.

We implemented it in MATLAB with the following code (we tried different values of N and different source functions)

```<code>N = Nvect(index);

M=2*N;

h = 2/(2*N+1);

k=2*pi/M;

for i=1:N

gv(i) =(i-0.5)*h;

end

for j=1:M

gv2(j)=(j-1)*k;

end

[meshr,meshteta] = meshgrid(gv,gv2);

meshr = meshr';

meshteta = meshteta';

f = @(r,teta) 9.*r;

rhs = f(meshr,meshteta);

u = @(r,teta) r.^3;

u_mesh = u(meshr,meshteta);

boundary=u(gv(N),gv2);

rhs=rhs*(h^2);

F = reshape(rhs,[M*N,1]);

for i=1:N

  lambda(i,1)=1/(2*(i-0.5));

end

for j=1:M

   F(j*N)=F(j*N)-(1+lambda(N))*boundary;

end

for j=1:N

beta(j,1)=1/((j-0.5)^2*(k^2));

end

I = speye(M);

I2=ones(N,1);

%I2=speye(m);

D=sparse(diag(beta));

d=ones(M);

T=spdiags([1-lambda, -2*I2, 1+lambda],[1,0,-1],N,N)';

Ap= spdiags([d,d,d,d],[1,-1,M-1,-(M-1)],M,M);

A=kron(I,T)-2*kron(I,D)+kron(Ap,D);

U = A\F;

U = reshape(U,N,M);

u_int = u_mesh(1:N,1:M);

e = max(max(abs(U-u_int)));

end
</code>```

We tried $N = [8,\space16,\space32,\space64,\space128,\space256,\space512,\space1024]$ and the error ratio $\frac{e_{N_{i+1}}}{e_{N_{i}}}$ stays always around 0.5 which shows a first order accuracy.
What could cause such inconsistency between theoretical and computational results? What did we get wrong?

EDIT: As a test, we used a source function s.t. it's easy to compute true solution $U$ by hand so we just compute the difference between $U$ and its approximation and take $e_{N_{i}}=max(U_i-U^{appr}_i)$.


问:
What matrix criteria is "large enough" for Eigen to use a BLAS backend?
Eigen's documentation says that any BLAS back-end can be used (e.g. MKL) to perform the actual matrix calculations.

It is, however, very vague when discussing the criteria as to which the calculations are offloaded (with added emphasis):

When doing so, a number of Eigen's algorithms are silently substituted with calls to BLAS or LAPACK routines. These substitutions apply only for Dynamic or large enough objects with one of the following four standard scalar types: float, double, complex, and complex. Operations on other scalar types or mixing reals and complexes will continue to use the built-in algorithms.

Unfortunately, there is no further guidance as to what is considered "large enough".

So my questions are:

What is the criteria to be considered "large enough"

Does it vary by platform?

Is this dynamically determined?

Can it be controlled by the user?

Background: I have small (25-50 square) matrices, but are large-ish on the embedded platform (Intel Atom Silvermont) we are using. Since this consumes a significant amount of CPU, I am looking to experiment with BLAS back-ends without having to substantially re-write code.


问:
2D Test cases for overset meshes
I developed a parallel 2D overset mesh (chimera) assembler and I would like to test a large scale problem. I don't want to make up a test case myself but rather want something standard. Would you refer me to such test cases?</p>

Title: 2D Test cases for overset meshes


问:
Parallelizing molecular simulation with full configuration energy
First, let just me stress that I'm not a an expert in computation chemistry, so now the problem:

We have GCMC molecular simulation, in the Grand Canonical ensemble, using the standard metropolis algorithm.

We would like to port the program to CUDA and run it in parallel.
There are several techniques to parallelize such simulation to system with local interaction. Unfortunately in our simulation we compute a total energy for the system taking into account all the particles in the system.
We use a primitive model, and given the elementary charge $e$ and the absolute permittivity of vacuum $\epsilon_0$, the interaction between two charges species $i$ and $j$ separated by a distance $r_{ij}$ is calculated acconrding to:
$$\begin{cases}
u^{\text{el}}(r_{ij}) = \frac{q_iq_je^2}{4\pi\epsilon_0\epsilon_sr_{ij}}\\
u^{\text{hs}}(r_{ij}) = \infty,& \quad \text{if } r_{ij} \leq R_i + R_j\\
u^{\text{hs}}(r_{ij}) = 0, &\quad \text{if } r_{ij} > R_i + R_j\\
\end{cases}
$$

The hard radii sphere is set to $2\mathring{A}$. A one-body external field was use to used to correct for the long range electrostatic interaction and to account for the cell boundary contraints along the z-direction acting as hard walls.:
$$
\begin{cases}
v^{\text{ext}}(z_i) = \infty,& \quad \text{when } z_i \geq L_z||\;z_i \leq <0\\
v^{\text{ext}}(z_i) = q_ie\varphi^{\text{ext}}(z_i), &\quad \text{when } 0 < z_i < L_z\\
\end{cases}
$$
where $\varphi^{\text{ext}}(z_i)$ is the potential produced by the external charge distribution.

The full configurational energy of the system is then:
$$ U^\text{conf}= \sum _i^{N_0+N_s} v^\text{ext}(z_i) + \sum_{1<j}^{N_0+N_s} u^\text{el}(r_{ij}) + u^\text{hs}(r_{ij})$$

The heaviest part of the simulation the the calculation of $U^\text{conf}$.

The main iteration is like this:

```<code>for i<-ion_0 to last_ion do:
   actual_energy <- compute_u_conf(i)
   propose_a_move(i)
   new_energy <- compute_u_conf(i)
   accept_or_reject_move(i)
</code>```

Here is where domain decomposition does not fit. Even if we can compute in parallel a move, then the computation of the energy would not be consistent because others particles are moving. This would break the Markov chain and make the system not convergent because the acceptance of the proposed move is based also on the total energy calculated.

It is possible to calculate the energy only locally and using some sort of mean for the rest of the space?
Is there a way to parallelize the calculation of the energy being consistent?
Have someone any idea on who can it be done?
Any help is welcome.


问:
Polynomial approximation - Vandermonde matrix creation - precision
I am trying to fit a polynomial through 340 points in a 3D space, i.e.

$$f(x,y,z) = k$$

I asked previously about the theory behind polynomial interpolation here ->
Polynomial interpolation

Few interesting points came out from the discussion, but I would like to understand how to correctly implement them in Python.

I am writing a script in Python to compute the Vandermonde matrix for a polynomial of a generic order $p$. Now the implementation problem arises.

First method

I have 340 points $(x,y,z)$. Supposing I have a polynomial of 6th order I have

(order+1) ^ number_of_dimensions

degrees of freedom for my polynomial. Therefore, it means that I have 343 degrees of freedom. Since I have 340 points I cannot use such order because otherwise the problem is underdetermined.

If I use a 5th order polynomial I have 216 degrees of freedom, and in this case the problem is overdetermined and the polynomial is an approximating polynomial. In a preliminary phase I do not need the exact interpolation, since I can solve the problem with a complementary solution, but do you think it is feasible the solution given?

Second method

As suggested by the user who helped me in the previous question, actually a polynomial of 6th order has 84 degrees of freedom. This is related to the constraint given by
$$a+b+c<p$$
with $a,b,c$ the exponents in each term $x^a y^b z^c$.

The question which arises is if it more feasible to create a Vandermonde matrix which respects this constraint. Does it give an improved approximated solution compared to a full Vandermonde matrix?


问:
Perfect filtering of high frequencies in 2D FFT (Multidimensional 2/3 Rule)
Let $u_n$ be an array containing discrete values of the function $u(x,y)$. Performing a 2D FFT to this array we obtain $\hat{u_n}$ representing the values of $\hat{u}(k_x,k_y)$. I would like to perform something like the 2/3 rule often used in pseudospectral methods for fluid dynamics:

$$\hat{u}(k)=0\qquad if \quad k>(2/3)k_{max}  $$

I know how to do this in 1D, but I am having trouble implementing it in 2D. I have tried setting to zero all $\hat{u}(k_x,k_y)$ for which $k=\sqrt{k_x^2+k_y^2}>(2/3)k_{max}$ but this method is affecting the energy of my lower frequencies, which I would like to be preserved without numerical dissipation. I would like to obtain something like this:

Which I have only been able to obtain setting a filter in terms of the dependent magnitude (Normalized Energy), and not regarding the wavenumber $k=\sqrt{k_x^2+k_y^2}$ as I would want.

I would appreciate if someone pointed me in the right direction to tackle my problem.


问:
Inverse problems with a discrete set of known parameters
What are the techniques on inverse problems to discover the distribution of parameters from a discrete set of values? For instance, I know that my domain where the PDE is defined is made up of materials A, B and C, but I do not know how they are distributed inside the domain. Working with discrete optimization algorithms seems unfeasible. What are the techniques to solve such a problem?

Edit: I will describe the problem for just two materials, as it seems that the strategy is different than for multimaterials, according to Clason's comments (bang-bang vs multibang)

$$
\underset{u \in U}{\text{min}} \frac{1}{2} \lVert F(u) - y^{\delta} \rVert_Y \\
$$
where
$$
U = \{ u \in X : u \in \{ u_1, u_2 \} ~\text{pointwise}\}\\
F: X \rightarrow Y
$$
$X$ and $Y$ are function spaces and $y^{\delta} \in Y$


问:
lightweight implementation of semiempirical quantum chemistry (e.g. MNDO,AM-1,PM3)
I'm searching for semi-empirical quantum chemistry solver which would be easy to integrate into my own software. I found a few implementations which can be in principle used e.g. MOPAC, ORCA, SQM some of which are integrated e.g. into Molecular Design Toolkit. But I don't want to depend on 1GB packages when I need only a few functions. Also I wan't something where I can dig-in and modify it more easily.

I wonder if there is something more convenient with respect to following aspects:

As much open as possible (GPL or better MIT licence)

compact and simple to understand code-base (small number of lines, clear readability - no fortran77, minimum of external dependencies)

C/C++ (optionally with python interface) ... I don't want to bother with frotran-C or fortran-python interface.

This is in principle continuation of my previous question
Simple open-source Quantum chemistry or DFT code in C/C++. While Libint is very nice software, full Hartree-Fock is too slow for me. I really would like something cheaper, like semiempirical methods and DFTB.


问:
Optimal Control using Dynamic Programming - Optimizing for Furthest Distance
So I have been investigating a problem to get a glider with control of its elevator to fly as far as possible from any given initial state. To keep this simple, we will view this in 2D space with the following differential equation:

\begin{align}
\dot{\boldsymbol{q}} = \dot{\begin{bmatrix}
x \\
y \\
\theta \\
\phi \\
\dot{x}\\
\dot{y}\\
\dot{\theta}
\end{bmatrix}} &= \begin{bmatrix}
\dot{x} \\
\dot{y} \\
\dot{\theta} \\
u \\
-\left(f_w \sin\theta + f_e \sin\left(\theta + \phi\right)\right)m^{-1} \\
\left(f_w \cos\theta + f_e \cos\left(\theta + \phi\right)\right)m^{-1} - g\\
\left(f_e (l \cos\phi + l_e ) - f_w l_w\right) I^{-1}
\end{bmatrix} \\
\end{align}

and the following relationships:

\begin{align}
f_w &= \rho S_w |\dot{\boldsymbol{x}}_w|^2 \sin\alpha_w \\
f_e &= \rho S_e |\dot{\boldsymbol{x}}_e|^2 \sin\alpha_e \\
\alpha_w &= \theta - \tan^{-1}\frac{\dot{y}_w}{\dot{x}_w} \\
\alpha_e &= \theta + \phi - \tan^{-1}\frac{\dot{y}_e}{\dot{x}_e} \\
\dot{x}_w &= \dot{x} + l_w \dot{\theta} \sin\theta \\
\dot{y}_w &= \dot{y} - l_w \dot{\theta} \cos\theta \\
\dot{x}_e &= \dot{x} + l\dot{\theta}\sin\theta + l_e\left(\dot{\theta} + \dot{\phi}\right)\sin\left(\theta + \phi\right)\\
\dot{y}_e &= \dot{y} - l\dot{\theta}\cos\theta - l_e\left(\dot{\theta} + \dot{\phi}\right)\cos\left(\theta + \phi\right) \\
\dot{\boldsymbol{x}}_w &= \dot{x}_w \hat{e}_x + \dot{y}_w \hat{e}_y \\
\dot{\boldsymbol{x}}_e &= \dot{x}_e \hat{e}_x + \dot{y}_e \hat{e}_y \\
\end{align}

where $u$ is the control, essentially a choice for $\dot{\phi}$, $\phi$ is the relative elevator angle with respect to the pitch angle $\theta$, $x$ and $y$ are the horizontal and vertical positions, $\dot{x}$ and $\dot{y}$ are the horizontal and vertical speeds, $\dot{\theta}$ is the angular velocity of the glider, $f_w$ is the net aerodynamic force magnitude from the main wing, $f_e$ is the net aerodynamic force magnitude from the elevator wing, $g$ is gravitational acceleration constant, and the other constants are tied to glider physical traits.

The values being used for the various constants are the following:

\begin{align}
m &= 0.05 \\
g &= 9.81 \\
\rho &= 1.292 \\
S_w &= 0.1 \\
S_e &= 0.025 \\
I &= 6 \cdot 10^{-3} \\
l &= 0.35 \\
l_w &= -0.03 \\
l_e &= 0.05
\end{align}

and the initial condition I primarily use to test is the following:

\begin{align}
\boldsymbol{q}_0 = \begin{bmatrix}
x_0 \\
y_0 \\
\theta_0 \\
\phi_0 \\
\dot{x}_0\\
\dot{y}_0\\
\dot{\theta}_0
\end{bmatrix} &= \begin{bmatrix}
0 \\
2 \\
0 \\
0 \\
6 \\
0\\
0
\end{bmatrix} \\
\end{align}

I am experimenting with using Dynamic Programming to tackle this problem when $u$ is constrained such that $-1 \leq u \leq 1$. Since Dynamic Programming is memory intensive for large state spaces, I recognized that for an optimal distance controller, I don't actually need the first two states, $x$ and $y$. With this change, I defined $\boldsymbol{q} = \lbrack \theta, \phi, \dot{x}, \dot{y}, \dot{\theta} \rbrack^T$ along with the associated differential equations truncation. I also define the discrete dynamical system using the following:

\begin{align}
\boldsymbol{q}_{k+1} &= \boldsymbol{q}_k + \Delta t \dot{\boldsymbol{q}}\left(\boldsymbol{q}_k, u_k\right)\\
&= f\left(\boldsymbol{q}_k,u_k\right)
\end{align}

To go along with this change, I made the overall optimization problem to maximize the following:

\begin{align}
V &= \sum_{i=1}^N \Delta t \dot{x}_i - \gamma \dot{y}^2_i\\
\text{subject to}& \begin{matrix}
-1 \leq u_k \leq 1 \\
\boldsymbol{q}_{k+1} = f\left(\boldsymbol{q}_k,u_k\right)
\end{matrix}
\end{align}

because the cost function should approximate the value for $x$ at the end of a flight, which is obviously what I would have optimized if using the full system of equations.

After doing some experiments, it seems the cost function chosen doesn't really work well for reasons I am unsure of. However, if I change the cost function to the following, it performs much better:

\begin{align}
V &= \sum_{i=1}^N \theta_i^2 + \gamma \dot{\theta}_i^2
\end{align}

for some $0 \leq \gamma \lt 1$. I chose this second cost function thinking one thing that might help a long flight is the glider remaining level instead of diving too soon and losing a lot of energy. It works decent, but I am still wondering why the first isn't doing too well.

With all this said, is there any problems with the first optimization formulation that stand out?


问:
Numerical integral of oscillating function with known zeros
I have a function that I need to numerically integrate from $0$ to $+\infty$, given by:

$$I = \int_0^{+\infty} \mathrm{d}x\,x\,T^2(x)f(x)$$

where $T^2$ is an interpolated function that goes to $1$ at small $x$ and decays as $x^{-4}$ at large $x$.

$f$ is an oscillating function with known zeroes:

$$f(x) = j^2_1(\alpha x)\times\cos(\beta\log x)$$

where $j_n$ is the spherical Bessel.

Since I know the zeros $z_n$ of the integrand (thus having control on its oscillations), is it a good idea to split the integral in sub-integrals (from $z_n$ to $z_{n+1}$) and sum all the contributions?

My idea was to change variables ($x = e^u$) and do the sub-integrals with SciPy's scipy.integrate.quad.


问:
Strange solutions using Finite Element Analysis
I've implemented the Finite Element Method to model the heat transfer between two different materials where one material is surrounded by the other.

When I run the model I'm getting some strange results. Whenever I run the model, I get what appears to be a sinc function through space and a linear change through time. My gut has been telling me that a part of my implementation isn't right but I haven't been able to figure out what is going on. The positional axis is the position of each individual element, not the radial distance.

The model that I'm currently using is a one dimensional model over the radius of a cylinder with heat being generated within the submerged material and the solution is trying to solve how heat moves through the whole system. Due to symmetry in the system and ignoring how heat moves through the length of the cylinder, I'm modelling the elements as a series of concentric circles with a known thickness. Thinking maybe it was a problem with how I implemented the model, I implemented a two dimensional model looking at a small square to get a proof of concept and I found a similar result in the spatial dimensions.

The 1 dimensional problem I'm solving is
$$\rho c_p \frac{\partial T}{\partial t} = k \frac{\partial}{\partial r}\left(r \frac{\partial T}{\partial r}\right)$$

The domain can be seen in the picture below. The inner radius has the set of parameters defined as $k_1, \rho_1, c_{p1}$ and the outer radius has the set of parameters defined as $k_2, \rho_2, c_{p2}$. The inner material is much smaller than the outer material, such that the image is not to scale.

The only boundary condition that I care about is the boundary between the two materials as I'm assuming the two materials are at the same temperature and the outer radius is insulated
$$k_1 \frac{\partial T}{\partial r}|_{r=R_1} = k_2 \frac{\partial T}{\partial r}|_{r=R_1}$$

The source term I'm using is a constant $Q$ and only being applied to the elements in the smaller radius once I get the model working properly I'll refine the source term some more. The large spike in the middle of the first image is the elements caused by the heat source - currently using two elements for now once I'll add more once things work.

I have no formal training with the Finite Element Model and the only understanding I have is from reading numerous books and papers I've scrounged through trying to understand the model, so the problem I'm experiencing could be a simple oversight on my part or a pitfall that I've known nothing about. If anyone could point me in the right direction on how to solve this problem I would be grateful.


问:
What can be done with Finite Element Method and not with the Finite Volume Method, and vice versa?
What are some applications where you would absolutely go for either FEM, but not FVM, or vice versa? What are some applications where both methods are equally suited?

I worked with the FEM so far and have written various codes, but I sometimes wonder whether I should also use FVM in future.


问:
Numerically solving generalized eigenproblem with Neumann conditions
I am interested in finding the eigenvalues/eigenfunctions of problems such as

$$ \partial_{xx} u = \lambda \partial_{yy} u, $$

which can be solved as the generalised eigenvalue problem

$$ \mathbf{A} \vec{u} = \lambda \mathbf{B} \vec{u}, $$

using the finite difference method. My question is about how to incorporate Neumann conditions into the $\mathbf{A}$ and $\mathbf{B}$ matrices. Examples I have found involved changing the right hand side of the equation in equations such as $ \mathbf{A} \vec{u} = \vec{b} $, i.e. not eigenvalue problems, but this is not possible in this case. I am not necessarily specifically interested in solving the PDE above, I am just looking for a general method.


问:
Integrating the 2d vorticity equation on periodic boundaries
This question is a follow-up of https://stackoverflow.com/questions/44718160/solve-a-linear-system-for-fft-coefficients

At some time (kt), the FT of the vorticity (omega) satisfies:

```<code>(1-nu*L*dt/2)fftomega(kt+1)=(1+nu*L*dt/2)fftomega(kt)-fft_J*dt
</code>```

where fftomega is the 2d fft of omega, J(psi,omega) is the Jacobian (see below) and the streamfunction del^2(psi)=omega. L is an array of coefficients that when multiplied by say fftomega return the fft of the laplacian of omega. For a 6x6 L is

0   -1   -4   -9   -4   -1

-1   -2   -5  -10   -5   -2

-4   -5   -8  -13   -8   -5

-9  -10  -13  -18  -13  -10

-4   -5   -8  -13   -8   -5

-1   -2   -5  -10   -5   -2

The fft of the Jacobian is computed as follows (using matlab/octave notation for brevity):

```<code>Dx_psi=real(ifft2(Dx.*fftpsi));                    %v
Dy_psi=real(ifft2(Dy.*fftpsi));                    %-u
Dx_omega=real(ifft2(Dx.*fftomega));
Dy_omega=real(ifft2(Dy.*fftomega));
Jacobian = - Dy_psi.*Dx_omega + Dx_psi.*Dy_omega;

fft_J = fft2(Jacobian);
</code>```

In practice, the Jacobian is integrated with a predictor corrector scheme due to Gadzag (the -fft_j*dt in the equation above is schematic). If arrays such as fftomega have the same dimension as omega, then the equation above can be solved as a system of linear equations, with LAPACK for instance.  However since the physical variables are real, procedures such as FFTW allow for output fft's to be non-square (r2c,c2r).  In this case it seems like the only way to integrate the equation is to expand these fft's back out to square, giving up the advantage afforded by the r2c,c2r scheme. Alternatively, one could make the physical variables complex with null imaginary part, and use the c2c scheme.

Has anyone found an alternative that makes it possible to use r2c,c2c?


问:
Best algorithm for inversion of matrix spanning many orders of magnitude
I have a very similar problem to the one described in Calculate inverse of dense matrix with entries of very different magnitude. The reason why I am opening a new question is because as far as I understood the answer given to the problem in that question was too specific for the problem considered.

My issue is that I have a matrix C which is fairly large (say typically 500x500 entries) and the entries are spanning some orders of magnitude, typically from $10^3$ to $10^{-16}$ and I need to INVERT this matrix. I really need to do this, as in fact in a following stage I need to use the inverse matrix in some multiplications with other matrices. This is also the reason why I am more interested in the PRECISION of the inversion rather than the performance in terms of speed.

Therefore my obvious question would be: what is the best algorithm to use to get the inverse matrix in this case?

I tried Singular Value Decomposition, but I noticed that it is producing unstable results.


问:
Scaling for a nonsymmetric eigenvalue problem
I have an eigenvalue problem emerging from the internal vibro-acoustic coupling.

The eigenvalue problem is nonsymmetric but it was proven in literature that it results in  real eigenvalues and eigenvectors due to the relations between the off-diagonal coupling matrices. The eigenvalue problem reads as,

$$A \phi = \lambda B \phi$$

where

$$ A = \begin{bmatrix} K_1 &K_2\\ 0 &K_3\end{bmatrix}$$

and

$$ B = \begin{bmatrix} M_1 &0\\ -K_2^T &M_2\end{bmatrix}$$

I am using the Arnoldi technique to solve this eigenvalue problem. However, I always have to apply a scaling in order to find the right eigenvalues and eigenvectors. What I am applying is the simple jacobi scaling such as

```<code>P1 = diag(1./sqrt(diag(K1)))
P2 = diag(1./sqrt(diag(K2)))

P = diag(P1,P2)
</code>```

eventually, what I am solving is the generalized eigenvalue problem with this scaling matrix, namely,

$$ P A P = \begin{bmatrix}
P_1 K_1 P_1 &P_1 K_2 P_2\\
0 &P_2 K_3 P_2\end{bmatrix}$$

and

$$ P B P = \begin{bmatrix}
P_1 M_1 P_1 &0\\
-P_2 K_2^T P_1 &P_2 M_2 P_2\end{bmatrix}$$

then I go back to the original space before scaling. It is fine up to this point but there is one more follow up step I would like to accomplish. Later on, I need to solve a linear system which involves $A$, such as,

$$A x = b$$

and I would like to use the factor I generated during the above eigenvalue solution in order to perform this solution. But if I perform the scaling, my original $A$ transforms to $PAP$ so I can not use that for a forward-backward solve and I need a new factorization for $A$.

And as mentioned above, I can only find the right eigenvalues and eigenvectors with this scaling although my original matrices are not that ill-conditioned. I am looking for a way to solve my eigenvalue problem without scaling and use the available factorization at a later stage, but so far, I have no luck without scaling. I am using the built in solvers (or the Factorize package from Tim Davies) in my MATLAB implementations for the Arnoldi iterations.

Any help/idea is greatly appreciated on what I can try to solve this problem without scaling.


问:
How to determine which parameter I should linearize?
Assumed I have the following equations with the parameters $\vec{A}(z, t)$, $\vec{B}(z, t)$ and the constants $a,\,b,\,c$:
$$
\begin{split}
\partial_tA&=a\partial_zA\\
\partial_tB&=b\partial_zB+A\cdot B\cdot c
\end{split}
$$
Now the usual approach is to split the equations for an implicit scheme like
$$
\left(\frac{1}{dt}-\frac{\partial_z}{2}\right)A_1=\left(\frac{1}{dt}+a\frac{\partial_z}{2}\right)A_0
$$
which then can be solved.
The split is done using
$$\begin{split}
\frac{A_1-A_0}{dt}&=\partial_tA\\
\partial_zA&=\frac{\partial_zA_1+\partial_zA_0}{2}\\
\frac{A_1-A_0}{dt}&=a\frac{\partial_zA_1+\partial_zA_0}{2}\\
\left(\frac{A_1}{dt}-a\frac{\partial_zA_1}{2}\right)&=\left(\frac{A_0}{dt}+\frac{a\partial_zA_0}{2}\right)\\
\left(\frac{1}{dt}-a\frac{\partial_z}{2}\right)A_1&=\left(\frac{1}{dt}+\frac{a\partial_z}{2}\right)A_0
\end{split}$$

For the equation system above I can then write (without the second term in the second equation):
$$\begin{split}
C&=\begin{pmatrix}
A\\B\end{pmatrix}\\
e&=\begin{pmatrix}
a\\b\end{pmatrix}\\
d&=\begin{pmatrix}
0\\0\end{pmatrix}\\
\left(\frac{1}{dt}-e\partial_z-\frac{d}{2}\right)C_1&=\left(\frac{1}{dt}+e\partial_z+\frac{d}{2}\right)C_0\end{split}$$
But when including the second term, I can choose between
$$d=\begin{pmatrix}B_0\cdot c\\0
\end{pmatrix}$$
or
$$d=\begin{pmatrix}0\\A_0\cdot c
\end{pmatrix}$$
depending on how I calculate
$$A\cdot B\cdot c$$
Is there a method to determine which $d$ I should use?


问:
Can x-ray back-projection be converted to hard-field magnetic induction tomography?
This is a question about hard-field back-projection as used in x-ray tomography, applied magnetic induction tomography. Al-Zeibak and Saunders have shown that x-ray filtered backprojection can be applied to an exciter and receiver coil, with a straight line of magnetic flux between them along the central axes of the 2 coils, taking the place of a straight line beam of an x-rays, to scan/image and an object (plastic container of saline solution).

In place of the summation of the x-ray attenuation coefficients $\mu(x,y)$,

$$p_\theta(r) = ln \lgroup \frac{I}{I_0} \rgroup = - \int \mu(x,y) ds$$

https://en.wikipedia.org/wiki/Tomographic_reconstruction

they appear to detect the summation of voltages at receiver coil:

$$p_\theta(r) = \int V(x,y) ds$$

Where $V(x,y)$ is the voltage in the object at $(x,y)$

Would it be possible to convert this (Al-Zeibak and Saunders method) to a square Helmholtz-coil assembly, providing a uniform AC magnetic field (i.e. straight lines of flux) and a square array of air-cored inductors as the receiving/sensor coils. So that a metallic object being scanned would generate a summation of voltages for each approximate straight line of flux of the Helmholtz coils?


问:
How can I numerically solve a saddle point problem with repeated constraints?
I am interested in numerically solving the following constrained minimization problem;  Find the value of $x\in \mathbb{R}^n$ that minimizes $f$ where
 $f\colon \mathbb{R}^n\to \mathbb{R}$ is defined as
$$
f(x) = \frac{1}{2}(x,A x)-(x,b),
$$
for the matrix $A$ being large sparse, symmetric and positive definite, subject to
$$
   B^\mathrm{T} x = g.
$$

Here the matrix $B\in \mathrm{R}^{n\times m}$ is large sparse, and $m$ is large but small compared to $n$.  So far so good, this is just the normal saddle point problem.  The wrinkle here is that $B$ may have linearly dependent columns, and that I know $g\in \mathrm{Range}(B^\mathrm{T})$, at least up to machine precision.

The problem comes from a finite element code where I generate the constraint matrices by against too many test functions.  I can explain more if anyone is interested. (I Should be testing against all elements in trace of some space,  S, but I don't have that space available to me easily so I test against all elements in S.

How can I efficiently solve this problem?  I am mainly interested in some sort of Krylov / Relaxation method.   Below I outline a method using the SVD of $B$, and a $QR$ method would work too.

This is what I would do for small matrices: (Spoiler Alert: this is very tedious and somewhat boring, I wrote it out in this much detail for my own reference, and to show that I did some work before asking here.)

Let $(U\in \mathbb{R}^{n\times n},\Sigma\in \mathbb{R}^{n\times m}, V\in \mathbb{R}^{m\times m})$ constitute the singular value decomposition of $B$, so that
$$
U\Sigma V^\mathrm{T} = B.
$$

If $r$ is the rank of $B$, then the first $r$ columns of $U$  (call then $U_1$) span the range of $B$ and the last $n-r$ columns of $U$ (call them $U_2$) span the nullspace of $B^\mathrm{T}$.  Similarly, the first $r$ columns of $V$ (call then $V_1$) span the range of $B^\mathrm{T}$ and the last $m-r$ columns (call them $V_2$) span the null space of $B$.  We also define $\Sigma_1$ as $r\times r$ matrix that sits at tthe top-left hand corner of $\Sigma$.

We now decompose $x$ in $x = \hat{x} + x_\perp$ where $\hat{x} \in \mathrm{Null}(B^\mathrm{T})$ and $x_\perp\in (\mathrm{Null}(B^\mathrm{T}))^{\perp}= \mathrm{Range}(B)$.  This means that $x_\perp = U_1\alpha$ for some $\alpha\in \mathrm{R}^{n-r}$ and that $\hat{x} = U_2\beta$ for some $\beta\in \mathbb{R}^{n-t}$.

Subbing this second identity into the constraint we get
$$
  V\Sigma^\mathrm{T} U^\mathrm{T} U_1\alpha = g,
$$

We can now multiply on the left by $V_1^\mathrm{T}$ which gives
$V_1^\mathrm{T}V\Sigma^\mathrm{T} U^\mathrm{T} U_1\alpha = V_1^\mathrm{T}g$, and noting that $V_1^\mathrm{T}V\Sigma^\mathrm{T} U^\mathrm{T} U_1 = \Sigma_1^\mathrm{T}=\Sigma_1$ we get $\alpha = \Sigma_1^{-1} V_1^{\mathrm{T}} g$ and $x_\perp$ is known.

Note that if we had multiplied on the left by $V_2^\mathrm{T}$ we have obtained $0 = 0$ since $g\in \mathrm{Range}(B^\mathrm{T})$ so all the "information" in the constraints was retained when we multiplied on the left by $V_1^\mathrm{T}$.

Then we define $F:\mathrm{R}^{n-r}\to \mathrm{R}$ as
$$
F(\beta) = \frac{1}{2}(x_\perp + U_2\beta,A(x_\perp+U_2\beta) )
- (b, x_\perp + U_2\beta)
$$

Now we find the $\beta$ that minimizes this functional (symmetric positive definite) and we can find $\hat{x}$.


问:
Solve ODE with two unknown functions
I want to solve a diffusion-convection problem numerically. I start from a PDE
of the form $$\partial_t f(x,t) = \partial_x(K(x)\;f(x,t))+\partial_{xx}\;f(x,t)^{m}\;.$$
It is possible to calculate a stationary solution $f_{st}(x)$, which looks something like this
$$
f_{st}(x) = (a \cosh(x-b))^c\;.
$$
Now I try to simplify the original PDE be assuming that the solution is similar to the stationary one, the only difference being that the coefficients $a$ and $b$ are now time dependent, i.e. $a(t)$ and $b(t)$.
This allows me to calculate the RHS of the PDE, reducing it formally to a first order ODE in time. My problem is that I'm not sure how to feed it to the computer, because due to the chain rule the spatial derivatives 'extract' additional instances of $a(t)$ and $b(t)$ out of $f$. In other words, I end up with a ODE like
$$
\partial_t(g(t)\;h(u(t)) = f(g(t),u(t))
$$
Intuitively I'd expect there to be a problem because I have two unknowns and just one equation but I'm trying to retrace a paper that seems to be doing just that. Anyways, thanks for reading, hope you can understand what I mean. </p>

Title: Solve ODE with two unknown functions


问:
Is there a FEM package for R
I am looking for a R package that performs finite element analysis. Is there a fenics or equivalent FEM package for R?</p>

Title: Is there a FEM package for R


问:
Minimizing the products of variables
My problem

Maximize

$$\min_{i} \{\ c_i \cdot \prod_{j \in A(i)} {x_{j}} \prod_{j \in B(i)} {y_{j}} \} $$

Subject  to

\begin{align}
&\sum_{j \in C(k)} x_{j} = 1,\ \forall k \\
&l \leq x_{j} \leq u,\ \forall j\\
&\sum_{j \in D(k)} y_{j} = 1+p,\ \forall k\\
&y_{j} \in \{1,\ p\},\ \forall j
\end{align}

where $c_i, l, u, p$ are positive constants less than 1, and sets $A, B, C, D$ are also given.

What I have tried

I have tried to use ncpol2sdpa to relax the polynomial programming into a semi-definite programming and call the sdpa solver to solve it.

The objective is replaced with a new variable $ \ f$, and the following constraints are added

$$\ f \leq c_i \cdot \prod x_{j} \prod y_{j}, \forall i\, .$$

Each discrete variable $y_i$ is replaced by $y_i = p \cdot (1-z_i) + z_i$, where $z_i = \{0,\ 1\}$ or equivalently, $z_i^2-z_i=0$.

However even for a small problem (5 $x$'s and 8 $y$'s), it took hours for ncpol2sdpa to relax and sdpa to solve the relaxed problem. The level of relaxation is set to be 3.

I also tried to relax the constraint on $y$ to be continues, but still quite slow.

I wonder if I did anything wrong with both tools? Or is there a better method/solver for this problem?

Title: Minimizing the products of variables


问:
Minimize interesting objective function with knowledge of gradient nonlinearity?
I plan on using a Quasi-Newton method (L-BFGS) to minimize a non-linear objective function.

$$ f: \mathbb{R}^n \rightarrow \mathbb{R}$$

The gradient is kind of interesting: as the values of the first set of parameters approach infinity, the gradient of the objective function with respect to the remaining parameters will approach zero. In other words,

$$ \lim_{x_i\to\infty}  \frac{\partial f}{\partial x_j} = 0$$
$$ i < j$$

The important concept: this relationship has no bearing on whether the objective function is being minimized or not. I fear the L-BFGS algorithm will follow a path where the gradient is minimized by this relationship instead of being minimized appropriately.

Can I use this knowledge to aid in my optimization strategy?


问:
Random Number Generation-full period
How can i demonstrate that a LCG has full period? Is there a way to demonstrate it or do I just need to check the parameters(a, c, m)?</p>

Title: Random Number Generation-full period


问:
Solving two coupled PDEs using two methods, how to calculate the error between them in theory?
Assumed I have the following coupled PDEs with the variables $A(z,t)$ and $B(z,t)$ and the constants $c_1,\,c_2$:
$$\begin{split}
\partial_tA&=\partial_zA+c_1\left(A+B\right)\\
\partial_tB&=\partial_zB+c_2\left(A+B\right)
\end{split}
$$
Now I have seen this algorithm for solving those equations:

First, calculate the first equation using
$$\left(1-\frac{dt}{2}\partial_z\right)A_{n+1}=\left(1+\frac{dt}{2}\partial_z\right)A_n+dt\cdot c_1(A_n+B_n)$$
and then use the solution in the second equation as:
$$\left(1-\frac{dt}{2}\partial_z\right)B_{n+1}=\left(1+\frac{dt}{2}\partial_z\right)B_n+dt\cdot c_2(A_{n+1}+B_n)$$
I would have approached it differently by combining $A$ and $B$ into a vector $C$:
$$\begin{split}
C_{n+1}&=\begin{pmatrix}A_{n+1}\\B_{n+1}\end{pmatrix}\\
C_{n}&=\begin{pmatrix}A_{n}\\B_{n}\end{pmatrix}
\end{split}$$
and then expanding the corresponding matrices:
$$
\begin{pmatrix}
1-\frac{dt}{2}\partial_z-dt\frac{c_1}{2}&-dt\frac{c_1}{2}\\
-dt\frac{c_2}{2}&1-\frac{dt}{2}\partial_z-dt\frac{c_2}{2}
\end{pmatrix}C_{n+1}=\begin{pmatrix}
1+\frac{dt}{2}\partial_z+dt\frac{c_1}{2}&dt\frac{c_1}{2}\\
dt\frac{c_2}{2}&1+\frac{dt}{2}\partial_z+dt\frac{c_2}{2}
\end{pmatrix}C_n
$$
Is my approach "more" correct than the initial one? If yes, how can I estimate my error and compare it to the original approach?


问:
Fit constant term in LM algorithm
I'm using the Levenberg-Marquardt algorithm to fit my data with a Gaussian function:
$$
f(x)=a\cdot e^{-\frac{(x-c)^2}{2\sigma^2}}+f_0
$$
$a$, $c$, $\sigma$ and $f_0$ are the fitting parameters. The algorithm correctly computes $a$, $c$ and $\sigma$. Given that
$$
\frac{\partial f}{\partial f_0} = 1
$$
and a $\delta \vec{x}$ solution of the least-squared problem:
$$
\delta \vec{x} = (\vec{J}\vec{J}^T)^{-1} \cdot\vec{J}f
$$
I get that $\delta \vec{x}$ is not sensible to any change in $f_0$; so it says nothing on what to do to minimize the $\chi^2$.

What am I missing? Do I have to switch to minimizing the $\chi^2$ function with the Brent algorithm (though it kinda spells trouble to me) to obtain the right known term or is there a more effective way to get there?

Or am I making it a lot harder than it is?

Title: Fit constant term in LM algorithm


问:
Double Integrating acceleration data to obtain position: 2 Problems
I have a data sample from an accelerometer from my phone (pretty bad accelerometer though). I'm trying to double integrate it in order to obtain the position as a function of time. I'm using a program in Python3 to integrate it, here is the code:

nm = input("Filename: ")

fil1 = open(nm,"r")
i = 0
s = 0.0
fil2 = open("int1.dat","w")

for line in fil1:
    xi,yi = line.split()
    if i == 0:
        x = float(xi) ; y = float(yi)
    if i > 0:
        s = s  + (float(xi)-x)*(float(yi)+y)/2
        x = float(xi)
        y = float(yi)
    i = i + 1
    fil2.write("{} {} \n".format(xi,s))
fil2.close()
#I just integrated the file. Now, i will integrate it again:
fil3 = open("int1.dat","r")
s = 0.0
i = 0
fil4 = open("int2.dat","w")
for line in fil3:
    xi,yi = line.split()
    if i == 0:
        x = float(xi) ; y = float(yi)
    if i > 0:
        s = s  + (float(xi)-x)*(float(yi)+y)/2
        x = float(xi)
        y = float(yi)
    i = i + 1
    fil4.write("{} {} \n".format(xi,s))
fil3.close()
fil4.close()
```

As you can see i'm approximating the area using trapezes.

The acceleration data was obtained in a movement as follows:

Started at a position that i defined as $x=0$.

Moved in straight line towards the end of an table, and stopped.

Waited a few seconds.

Returned to the original position

So, the graph of position x time should me something similar to a "hill": it goes up, than it goes down. It's value $x(t)$ at the peak would represent the total length of the table.

However, things are not working quite well. I took out the noise by removing the mean value of all $|a(t)| < 0.15$. However, it still did not sum up well. The graphs get pretty bad. So, here is my first problem:
The noise stays too big when i integrate it.
And so, my first question is: How can i reduce the noise?

Anoter problem that is happening and that i'm finding odd is that my position doesn't "go and come back" to where it started. It starts in $0$, goes to $-7$ (which is approximately the length of the table) and then, instead of coming back to $0$, it goes all the way down to $-14$. Why is this? I can't seem to find any particular reason for the integral not reckognizing the fact that i am coming backwards.

I tried to upload the graphics, but i couldn't right now. I will upload them as soon as i can. Much appreciate.


问:
Why can I not solve the negative advection equation (backwards in time)?
Suppose we have the negative, inhomogeneous advection equation:
$$\left(\frac{\partial}{\partial x}-\frac{1}{c}\frac{\partial}{\partial t}\right)v(t,x)=u(t,x)\qquad(t\in\mathbb{R}_{+},x\in\mathbb{R})$$
I want to solve it numerically, backward in time and forward in space.

Now, I have tried a few schemes without avail. For example, the Crank-Nicolson scheme going backwards in time:

In particular, we may discretise as

$$\frac{1}{c}\frac{\mathrm{d}}{\mathrm{d}t}v_\ell(t)=\underbrace{\frac{v_{\ell+1}(t)-v_\ell(t)}{\Delta x}-u_\ell(t)}_{G_\ell(t)}+\mathcal{O}(\Delta x)\qquad\Delta x\to 0$$
where we have taken the forward difference of the spatial derivative.

Then applying the Crank-Nicolson scheme yields

$$\frac{1}{c}\frac{v^{n+1}_\ell-v^n_\ell}{\Delta t}=\frac{1}{2}[G_\ell(t^{n+1})+G_\ell(t^n)]=\frac{1}{2}\left[\frac{v^{n+1}_{\ell+1}-v^{n+1}_\ell+v^n_{\ell+1}-v^n_\ell}{\Delta x}-u^{n+1}_\ell-u^n_\ell\right]$$
Hence
$$v^{n+1}_\ell-\frac{c\Delta t}{2\Delta x}(v^{n+1}_{\ell+1}-v^{n+1}_\ell)=v^n_\ell+\frac{c\Delta t}{2}\left[\frac{v^n_{\ell+1}-v^n_\ell}{\Delta x}-u^{n+1}_\ell+u^n_\ell\right]$$
i.e.
$$\left(\mathrm{I}-\frac{c\Delta t}{2\Delta x}\mathrm{A}\right)v^{n+1}_\ell=\left(\mathrm{I}+\frac{c\Delta t}{2\Delta x}\mathrm{A}\right)v^n_\ell-\frac{c\Delta t}{2\Delta x}(u^{n+1}_\ell+u^n_\ell)$$
where
$$\mathrm{A}=\begin{pmatrix}
-1 & 1 & & 0
\\
&\ddots & \ddots
\\
& & -1 & 1
\\
0 & & & -1
\end{pmatrix}$$
Now, since we want to go backwards in time, we replace $\Delta t$ with $-\Delta t$ and we get
$$v^{n-1}_\ell=\left(\mathrm{I}+\frac{c\Delta t}{2\Delta x}\mathrm{A}\right)^{-1}\left[\left(\mathrm{I}-\frac{c\Delta t}{2\Delta x}\mathrm{A}\right)v^n_\ell+\frac{c\Delta t}{2\Delta x}(u^{n-1}_\ell+u^n_\ell)\right]$$
Now, I try to implement this on Matlab via

```<code>%% Parameters
L = 5; % size of domain
T = 5; % measurement time
dx = 1e-2; % spatial step
dt = 1e-3; % time step
x0 = 0; % point of measurement
c = 1; % speed of advection
%%
t = 0:dt:T; % time vector
x = [0:dx:L]'; % position vector
nt = length(t); % number of time steps
nx = length(x); % number of position steps
mu = dt/dx;
I = eye(nx,nx+1); % identity matrix
A = spdiags(ones(nx,1)*[-1 1],0:1,nx,nx);
%% Solve Backward Advection Equation
% preallocate the memory
v = zeros(nx,nt);
u = zeros(nx,nt);
% final condition, sinc function centered around x = .5
v(:,nt) = sinc((x-x0)/dx);
for k = nt-1:1
    v(:,k-1) = (I+(c/2)*mu*A)\((I-(c/2)*mu*A)*v(:,k)+(c/2)*mu*(u(:,k-1)+u(:,k)));
end
</code>```

My solution simply "blows up" on the boundary of the domain and is zero everywhere else. I have similar results using backwards Euler instead of Crank-Nicolson. Where am I going wrong?

I provide two plots:


问:
Simulating periodic boundary condition with long_range interaction
Consider 2 particles in a 2D plane. There is a long-range interaction between two particles:
$$f (r_i-r_j)= \frac{\vec{r_i} - \vec{r_j}}{|r_i-r_j|^3}$$
$$\vec{r_i} =(x_i,y_i)$$

$r_i$ determines the position of particle number $i$ in the plane. Knowing position of particles in the 2D plane and also knowing that distances of particles is more than a threshold $r_0$ (so $f$ do not get infinite), how can one find $f_{tot}$ in the position of all particles in a periodic boundary condition simulation?

$f_{tot} $ is sum of effect of images of $j$ particle in f:

$$f_{tot}(\vec{r_i})= \sum_{n=-\infty}^{\infty} \bigg ( \frac{(x_i - x_j + n L,0)}{|y_i - y_j + n L|^3} \bigg) +  \sum_{n=-\infty}^{\infty} \bigg ( \frac{(0,y_i - y_j + n L)}{|y_i - y_j + n L|^3} \bigg)$$

in above relation, r_i and r_j are known and n is an integer number which changes from -infinity to infinity


问:
Comparison of different methods for solving $\partial_tA=\nabla^2_xA$
For testing I would like to solve the equation
$$
\partial_tA=\nabla^2_xA,\,A_0=\sin(x)
$$
In order to solve it, I use an implicit solver
$$\left(1-\frac{\nabla^2_x}{2}dt\right)A_{n+1} = \left(1+\frac{\nabla^2_x}{2}dt\right)A_n
$$
or the Runge-Kutta-method in two different orders. For solving the whole equation I use the following matlab script:

```<code>clear all;
warning('off', 'all');
x_max = 1;
x_num = 1000;
a0_val = 1;
b0_val = 1;
dx = x_max/20000;
steps = 1/dx;

x_vec = linspace(0, x_max, x_num).';
A0_vec = sin(x_vec);%linspace(1, x_num, x_num).';
ones_vec = ones(x_num, 1);
dh = (abs(x_vec(2)-x_vec(1)));
nabla_matrix = spdiags([-1*ones_vec 16*ones_vec -30*ones_vec 16*ones_vec -1*ones_vec], -2:2, x_num, x_num);
nabla_matrix(1, 1) = 35;
nabla_matrix(1, 2) = -104;
nabla_matrix(1, 3) = 114;
nabla_matrix(1, 4) = -56;
nabla_matrix(1, 5) = 11;
nabla_matrix(2, 1) = 11;
nabla_matrix(2, 2) = -20;
nabla_matrix(2, 3) = 6;
nabla_matrix(2, 4) = 4;
nabla_matrix(2, 5) = -1;
nabla_matrix(x_num, x_num) = 35;
nabla_matrix(x_num, x_num-1) = -104;
nabla_matrix(x_num, x_num-2) = 114;
nabla_matrix(x_num, x_num-3) = -56;
nabla_matrix(x_num, x_num-4) = 11;
nabla_matrix(x_num-1, x_num) = 11;
nabla_matrix(x_num-1, x_num-1) = -20;
nabla_matrix(x_num-1, x_num-2) = 6;
nabla_matrix(x_num-1, x_num-3) = 4;
nabla_matrix(x_num-1, x_num-4) = -1;
nabla_matrix = nabla_matrix / (12*(dh*dh));

A_equation = @(A_vec) ((nabla_matrix*A_vec.'));

A_matrix_single = speye(x_num)-dx/2.*nabla_matrix;
B_matrix_single = speye(x_num)+dx/2.*nabla_matrix;

A_matrix = zeros(steps, x_num);
A_matrix_NL_rk4 = zeros(steps, x_num);
A_matrix_euler = zeros(steps, x_num);
A_matrix_rk4 = zeros(steps, x_num);
A_matrix_rk42 = zeros(steps, x_num);

A_matrix(1,:) = A0_vec;
A_matrix_NL_rk4(1,:) = A0_vec;
A_matrix_euler(1,:) = A0_vec;
A_matrix_rk4(1,:) = A0_vec;
A_matrix_rk42(1,:) = A0_vec;

%Do steps
for i = 2:steps
    A_k1 = A_equation(A_matrix_rk4(i-1,:));
    A_k2 = A_equation(A_matrix_rk4(i-1,:)+A_k1.'*dx/2);
    A_k3 = A_equation(A_matrix_rk4(i-1,:)-A_k1.'*dx+2*dx*A_k2.');
    A_matrix_rk4(i,:) = A_matrix_rk4(i-1,:)+dx*(A_k1/6+2/3*A_k2+A_k3/6).';
end

for i = 2:steps
    A_k1 = A_equation(A_matrix_rk42(i-1,:));
    A_k2 = A_equation(A_matrix_rk42(i-1,:)+A_k1.'*dx/2);
    A_k3 = A_equation(A_matrix_rk42(i-1,:)+0.5*dx*A_k2.');
    A_k4 = A_equation(A_matrix_rk42(i-1,:)+dx*A_k3.');
    A_matrix_rk42(i,:) = A_matrix_rk42(i-1,:)+dx*(A_k1/6+1/3*A_k2+A_k3/3+A_k4/6).';
end

for i = 2:steps
    A_k1 = A_equation(A_matrix_NL_rk4(i-1,:));
    A_k2 = A_equation(A_matrix_NL_rk4(i-1,:)+A_k1.'*dx/2);
    A_k3 = A_equation(A_matrix_NL_rk4(i-1,:)+0.5*dx*A_k2.');
    A_k4 = A_equation(A_matrix_NL_rk4(i-1,:)+dx*A_k3.');
    A_matrix_NL_rk4(i,:) = speye(x_num)\(speye(x_num)*A_matrix_NL_rk4(i-1,:).'+speye(x_num)*(dx*(A_k1/6+1/3*A_k2+A_k3/3+A_k4/6)));
end

for i = 2:steps
    A_matrix(i,:) = A_matrix_single\(B_matrix_single*A_matrix(i-1,:).');
end

figure;
imagesc(A_matrix);
title('Forward step')
figure;
imagesc(A_matrix_rk4);
title('RK4 step')
figure;
imagesc(A_matrix_rk42);
title('RK42 step')
figure;
imagesc(A_matrix_NL_rk4);
title('RK4_NL step')
</code>```

but the results are way of my expectations. I get a more or less correct result when using the implicit approach, but the latter three all explode.

Why is that? Is my general approach flawed?


问:
Time-dependent Schrodinger equation(time-dependent Hamiltonian)
This question was already asked here and this is a suitable form of the equation for numerically solving.
\begin{align*}
i\frac{\partial}{\partial t}u_{\ell}(r,t) = \Bigg(-\frac{1}{2} \frac{\partial ^2}{\partial r^2} + \frac{\ell(\ell+1)}{2r^2} +  V(r)\Bigg)u_{\ell}(r,t)\quad \quad \quad \quad \quad \quad \quad \nonumber\\  + \frac{2}{3}k(t)r^2 \sum_{\ell' = {max}(\ell-2,0)}^{{min}(\ell+2,L_{{max}})} \Bigg(\delta_{\ell,\ell'} - \sqrt{ \frac{4\pi}{5}} \alpha(\ell, \ell') \Bigg) u_{\ell'}(r,t)
\end{align*}
where the coefficient $\alpha(\ell, \ell')$ is given by
\begin{align*}
 \alpha(\ell, \ell') = \int Y^*_{\ell0}(\Omega)Y_{20}(\Omega)Y_{\ell'0}(\Omega) d\Omega
\end{align*}

which the integral can be solved by Wigner-3j coefficients (see  here ).

Can anyone suggest an efficient numerically method for solving the equation?


问:
How to implement boundary conditions in arbitrary geometries, in MD\Lagrangian FD simulations
Let's say I have a random arbitrary geometry, like this one https://ibb.co/dG8fUa (something I just made up). I make one of these in 3D in SolidWorks. I have molecules/particles coming in through the inlet and leaving out of the two outlets. In a standard geometry, like a pipe or square channel, a simple way to implement a wall BC is something like

```<code>if(y_particle>0.3)
  v_particle=-v_particle
</code>```

assuming a square channel height of 0.15 and the origin along the lengthwise axis. This obviously wont work for a complicated geometry with wildly varying X and Y. How might one implement codes which can work on any generated geometry? I am trying for MD simulations through such geometries.


问:
Alternating Schwarz does not converge without Dirichlet conditions on physical boundaries
Note: Thanks to comments, I realized that I have two problems, each which can be described more clearly on its own. This revised question covers the first.

I would like to solve
$$
-\Delta u - 1=0,\quad u(0,y)= 0,\quad u(1,y)=0
$$
over a rectangle using alternating Schwarz. The following is my mesh, where the small disks indicate the internal boundary nodes of the subdomains (the nodes that have Dirichlet conditions that I update on them):

$x = 0$ and $x = 1$ are the physical boundaries of the full mesh. The local solutions do not converge to the global solution. The problem is mostly in the overlapping part of the mesh, here is the difference between the actual solution and the solution given by the local solutions:

If I impose the Dirichlet conditions $u(x,0) = 0$ and $u(x,1/\phi) = 0$, where $1/\phi$ is the height of the mesh, then the local solutions converge to the global solutions.

Should it not also work with only the boundary conditions that I mentioned about earlier?

EDIT:
I solved it by selecting the artifical boundary in this manner:


问:
Roots of transcendental equation involving bessel functions
I need an efficient way to numerically find the first $n$ positive roots $\lambda_n$ of the transcendental equation

$$ \dfrac{J_0 (\lambda_n r) Y_1 (\lambda_n) - J_1 (\lambda_n) Y_0 (\lambda_n r)}{J_0 (\lambda_n r) Y_0 (\lambda_n) - J_0 (\lambda_n) Y_0 (\lambda_n r)} = \dfrac{\delta}{\lambda_n}$$

where $\delta$ is a constant.

I have tried newtzero and chebfun in MATLAB. They are both taking longer than I'd like. The roots are needed for the solution of a PDE (transient heat transfer in one spatial dimension) with time varying boundary conditions using eigenfunction expansions. If anyone knows of a technique to determine the roots even a little faster, it would be much appreciated. Thank you.


问:
Approximating derivative of numerical solution (using ode45) within event function
The issue I have is having to compute the derivative (in real time) of the solution produced by ode45 within the events function.

Some pseudo-code to explain what I'm mean is,

function dx = myfunc(t,x,xevent)
persistent xevent
% xevent is the solution at the event
dx(1) = x(2);
dx(2) = complicated function of x(1),x(2), and xevent;
end

function [value,isterminal,direction] = myeventfcn(t,x)
position = function of x(1), x(2), and dx(2);
isterminal = 1;
direction = 0;
end
```

I know that if I didn't need to use the solution at the event within myfunc I could just compute dx=myfunc(t,x) within my event function and use dx(2), yet since xevent is used within myfunc I can't input xevent.

I know there is a way of inputting constant parameters within the event function, but since it's the solution at the event location, which also changes, I'm not sure how to go about doing that.

My work around to this is to approximate dx(2) using the solution x. What I wanted to know is if it will be satisfactory to just use a finite difference approximation here, using a small fixed step size relative to the step size od45 takes, which is a variable step size.

Thanks for any help!


问:
Finite difference scheme for unconfined aquifer equation
For an unconfined aquifer we have this PDE for the water table position( of course after somehow making the original Boussinesq equation linearized ):

$$ \frac{\partial^2(h^2)}{\partial x^2} + \frac{\partial^2(h^2)}{\partial y^2} = \frac{2s}{k} \frac{\partial h}{\partial t} - \frac{2N}{k}$$

where $N$, $k$ and $s$ are constants .

Would you tell me what type mathematically it is? In respect to $h^2$ it seems elliptic, but physically it is not. Can we treat it and solve it like 2D heat equation? I need to solve it numerically via finite difference. Which FD scheme would you recommend?

I tried it with the simple explicit  and it works OK and I'm aware of instability restrictions. but with methods like , ADE  (alternating-direction explicit), because the second derivatives are in power 2  , it will lead to quadratic equation for each node in each time.


问:
How to formulate Poisson's equation into flux eqution
I have a small 2D system I'm trying to model using a non-linear extension of Darcy's law for fluid flow in porous media. I'm primarily interested in the local flow velocity, not necessarily the pressure, that's later used in the real model I'm working on:

The original equation known as Darcy's Law used in ground water modeling is essentially Poisson's equation, where $\vec{q}$ is the fluid velocity, $\mu$ and $\kappa$ are constants that together describe the hydraulic conductivity of the porous medium, $p$ is local hydrostatic head and f is a source:

$$ -\frac{\mu}{\kappa} \nabla p = \vec{q} $$

$$ \nabla \cdot ( \vec{q} ) = f = -\frac{\mu}{\kappa} \nabla^2 p $$

The non-linear equation I want to solve instead is the following, known as the Darcy–Forchheimer law. As you can see it's essentially a polynomial instead of a linear model like Darcy's law:

$$ \frac {\partial p}{\partial x} = \frac {\mu }{\kappa} \vec{q} - \frac {\rho }{\kappa _{1}}\vec{q}^{2} $$

I've picked up on some authors [1] that go through the calculus and have figured out a way to represent $\vec{q}$ as a function of $p$ by calculating the inverse of the flow/pressure equation, and I see somewhat how they arrive at their methodology, but what I'm interested in is actually calculating the flux ($\vec{q}$) itself, which I use elsewhere in the model I'm developing - I don't have much of a use for $p$ besides in the initial problem setup with some simple Dirichlet boundary conditions.

Somewhat related, I've also read that for some numerical methods such as DG FEM, diffusion dominated problems can very unstable due to the non-directional flux. In this paper [2], the authors investigate some 'equation splitting' methods where they solve one equation for the state variable, and another for the flux, essentially the first two equations I showed above. This apparently relieves some of the instability issues caused by diffusion

Looking at these two, it looks like I have two options:

Use the methodology of [1] to calculate the pressure throughout the system, and then calculate the gradient of the hydrostatic head to generate the local velocity, Or...

Split my equation into two equations and solve them together. I end up having the solve two equations at once, but I need to have the velocity regardless so that doesn't seem so bad. Other than that I know nothing about this method.

I don't know how to implement the equation splitting method as I've never done it before, so I have two questions:

Is this even a good idea? Is splitting the equations up into double the number of equations disadvantageous for any other reason than doubling the number of state equations to solve?

Is there a source that goes into the actual implementation of equation splitting for (preferably) the finite element or finite difference methods. I am thoroughly intimidated by any DG theory I come across, and I'm hoping there is an easier way to solve this problem

[1] https://arxiv.org/pdf/1508.00294.pdf

[2] https://www.cs.utah.edu/~kirby/Publications/Kirby-11.pdf


问:
Optimizing for loops by using Einstein summations instead
I have some Python code containing a couple of for loops, which I would like to optimize by using low-level functions. My first approach is to use numpy.einsum(), but the results don't match the for-loop implementation:

## first term
T11 = 0.125 * numpy.einsum('abcd,efab,cdef,abcd,efcd',
    g[:Nelec,:Nelec,Nelec:,Nelec:], # <ab||rs> -> abcd
    g[:Nelec,:Nelec,:Nelec,:Nelec], # <cd||ab> -> efab
    g[Nelec:,Nelec:,:Nelec,:Nelec], # <rs||cd> -> cdef
    D,                              # e_a + e_b - e_r - e_s -> abcd
    D,                              # e_c + e_d - e_r - e_s -> efcd
)

T12 = 0
for a in range(Nelec):
    for b in range(Nelec):

        for c in range(Nelec):
            for d in range(Nelec):

                for r in range(Nvirt):
                    for s in range(Nvirt):

                        # nominator
                        nom  = g[a,b,Nelec+r,Nelec+s]
                        nom *= g[c,d,a,b]
                        nom *= g[Nelec+r,Nelec+s,c,d]

                        # denominator
                        denom = D[a,b,r,s] * D[c,d,r,s]

                        T12 += nom*denom

T12 *= 0.125

## second term
T21 = 0.125 * numpy.einsum('abcd,cdef,efab,abcd,abef',
    g[:Nelec,:Nelec,Nelec:,Nelec:], # <ab||rs> -> abcd
    g[Nelec:,Nelec:,Nelec:,Nelec:], # <rs||tu> -> cdef
    g[Nelec:,Nelec:,:Nelec,:Nelec], # <tu||ab> -> efab
    D,                              # e_a + e_b - e_r - e_s -> abcd
    D,                              # e_a + e_b - e_t - e_u -> abef
)

T22 = 0
for a in range(Nelec):
    for b in range(Nelec):

        for r in range(Nvirt):
            for s in range(Nvirt):

                for t in range(Nvirt):
                    for u in range(Nvirt):

                        # nominator
                        nom  = g[a,b,Nelec+r,Nelec+s]
                        nom *= g[Nelec+r,Nelec+s,Nelec+t,Nelec+u]
                        nom *= g[Nelec+t,Nelec+u,a,b]

                        # denominator
                        denom = D[a,b,r,s] * D[a,b,t,u]

                        T22 += nom*denom

T22 *= 0.125

## third term
T31 = numpy.einsum('abcd,edfb,cfae,abcd,aecf',
    g[:Nelec,:Nelec,Nelec:,Nelec:], # <ab||rs> -> abcd
    g[:Nelec,Nelec:,Nelec:,:Nelec], # <cs||tb> -> edfb
    g[Nelec:,Nelec:,:Nelec,:Nelec], # <rt||ac> -> cfae
    D,                              # e_a + e_b - e_r - e_s -> abcd
    D,                              # e_a + e_c - e_r - e_t -> aecf
)

T32 = 0
for a in range(Nelec):
    for b in range(Nelec):
        for c in range(Nelec):
            for r in range(Nvirt):
                for s in range(Nvirt):
                    for t in range(Nvirt):

                        # nominator
                        nom  = g[a,b,Nelec+r,Nelec+s]
                        nom *= g[c,Nelec+s,Nelec+t,b]
                        nom *= g[Nelec+r,Nelec+t,a,c]

                        # denominator
                        denom = D[a,b,r,s] * D[a,c,r,t]

                        T32 += nom*denom

print(T11, T12)
print(T21, T22)
print(T31, T32)
print(T11+T21+T31, T12+T22+T32)
```

<p>For my test case I get the following output:</p>

``` 0.016038311274  0.016038311274
 0.014256694984  0.014256694984
-0.044483006105 -0.044483006105
-0.014187999845 -0.014187999845
```

The first column contains the einsum() results, the second column is for the correctly working for loop implementation.

So my question is: Why does einsum() not give the same results (except for the second term $T_2$)? Is there some error in my einsum() implementation or is there a conceptual problem here?

Edit: After fixing one error in T11 and another error in T31 both methods yield same results.

Not necessarily  important, but the equations are for the Møller-Plesset perturbation theory third-order energy correction:

\begin{equation}
 E_0^{(3)} = T_1 + T_2 + T_3 \\
 T_1 = \frac{1}{8} \sum_{abcdrs} \frac{\langle ab||rs\rangle\langle cd||ab\rangle\langle rs||cd\rangle}{(\epsilon_a+\epsilon_b-\epsilon_r-\epsilon_s)(\epsilon_c+\epsilon_d-\epsilon_r-\epsilon_s)} \\
 T_2 = \frac{1}{8} \sum_{abrstu} \frac{\langle ab||rs\rangle\langle rs||tu\rangle\langle tu||ab\rangle}{(\epsilon_a+\epsilon_b-\epsilon_r-\epsilon_s)(\epsilon_c+\epsilon_d-\epsilon_t-\epsilon_u)} \\
 T_3 = \sum_{abcrst} \frac{\langle ab||rs\rangle\langle cs||tb\rangle\langle rt||ac\rangle}{(\epsilon_a+\epsilon_b-\epsilon_r-\epsilon_s)(\epsilon_a+\epsilon_c-\epsilon_r-\epsilon_t)}
\end{equation}

In the code I use:

g[i,j,k,l] = $\langle ij||kl\rangle$ being the anti-symmetrized two electron repulsion integrals

D[a,b,r,s] = $\frac{1}{\epsilon_i+\epsilon_j-\epsilon_r-\epsilon_s}$ with $\epsilon_i$ being the energy of orbital $i$

Nelec for the number of electrons, which equals the number of occupied spin orbitals

Norb for the total number of spin orbitals

Nvirt = Norb - Nelec for the number of virtual (unoccupied) orbitals

orbitals are ordered by energy, thus occupied first, then virtuals

And orbital indices follow the convention:

a,b,c,d running over occupied orbitals only

r,s,t,u running over virtual orbitals only

i,j,k,l running over all orbitals


问:
Iterative linear solvers compatible with automatic differentiation?
I'm using automatic differentiation on a function that contains a sparse nonsymmetric linear system to be solved. I was using BiCGStab to solve this part of the function, but noticed the derivatives produced by the automatic differentiation that were carried through BiCGStab were diverging in some cases. Is there an iterative linear solver with proven convergence using automatic differentiation?</p>


问:
Second Quantization in Matlab
This question may be more suited for physics.stackexchange, but I saw this post was recommended for StackOverflow or Computational Science, so I'm asking my question here.

I am trying to write a program in Matlab to evaluate an expression with fermionic annihilation and creation operators. Define
$$\{f, g \} \equiv fg + gf$$
then the fermionic second quantization operators satisfy
$$\{a_i,a_j \} = \{ a_i^\dagger, a_j^\dagger\} = 0 \qquad \{a_i^\dagger, a_j\} = \delta_{ij}$$
$\forall i, j \in \{1, 2, ..., N \}$ for some given $N$.

$W(t)$ is a $2N \times 2N$ matrix that I know that gives me the time evolution of the operators in the Heisenberg picture;
\begin{align*}
a_i(t) &= \sum_{j=1}^N \left(W_{2i-1,2j-1}(t)a_j + W_{2i-1,2j}(t)a_j^\dagger\right)\\
a_i^\dagger(t) &= \sum_{j=1}^N \left(W_{2i,2j-1}(t)a_j + W_{2i,2j}(t)a_j^\dagger\right)
\end{align*}
i.e.
$$\left(\begin{matrix}a_1(t)\\a_1^\dagger(t)\\\vdots\\a_N(t)\\a_N^\dagger(t) \end{matrix}\right) = W(t)\left(\begin{matrix}a_1\\a_1^\dagger\\\vdots\\a_N\\a_N^\dagger \end{matrix}\right)$$

I want to write a function that takes in $W(t)$ (evaluated at a given $t$) and $j$ and returns
$$\langle \psi | \left(\prod_{i=1}^{j-1}\left(1-2a_i^\dagger(t)a_i(t) \right) \right)\left(a_j(t)^\dagger + a_j(t) \right) | \psi \rangle$$
where
$$| \psi \rangle = \left(\prod_{i=1}^N\frac{1+a_i^\dagger}{\sqrt{2}} \right)|0\rangle = 2^{-N/2}\left(|0...0\rangle + |0...1 \rangle + |0...1 \ 0 \rangle + |0 ...1 \ 1 \rangle + ... + |1...1 \rangle \right)$$
$|\psi\rangle$ is just a normalized, equally weighted linear combination of all states. Note that it does not change with time. The time evolution is encoded in the operators.

My thought was to use Matlab's symbolic algebra, and try and code for what each operator does to a state. But the symbolic algebra does not always preserve order, which is crucial because the operators don't always commute. Also, trying to just naively go through is so computationally expensive and grows exponentially with $N$.

Any ideas?

Title: Second Quantization in Matlab


问:
Do DG methods for the Helmholtz equation always return positive quantities?
Helmholtz Diffusion equation with reaction term:
$$
k\Delta u + u = f ~ \text{in} ~\Omega
$$
$$
\nabla u \cdot \mathbf{n} = 0 ~ \text{in} ~\partial \Omega
$$

For sufficiently small $k$ (relative to the mesh), there can be oscillations (with Lagrange elements) that return a negative $u$ in some regions, specially when $f$ is a step function. I tried to apply a DG method with the Symmetric Interior Penalty formulation, and I am still seeing these negative values. Is this possible to happen? I was told that DG methods preserve positivity and cannot return negative values. Is there any way to preserve positivity in DG methods for elliptic equations?

Preserve positivity: If $f(x)>0 ~\forall x\in\Omega$, then $u(x)>0 ~\forall x\in\Omega$

Edit again: What if the boundary conditions are of zero flux? Edited the problem.


问:
Pseudo-Spectral cosine transform
I'm trying to solve the following equation
$$u_t = u_{xx} + u(1-u^2), u_x(\pm 1) = 0,$$
using the Fourier cosine transform. The nonlinear term gives a convolution which I would rather avoid, which is why I wanted to try a pseudo-spectral method.

In MATLAB I have the following script (simplified as much as possible)

```<code>function TestCosineSeries()

Nx = 32;
x = linspace(-1,1,Nx)';
f = tanh(1000*x);

dt = 1e-3;
Nt = 100;

for k = 1 : Nt

    t = k*dt
    [t,y] = ode45(@(t,y)RHS(t,y,L),[0 dt],f);
    f = y(end,:);

end

end

function duout = RHS(t,u,L)

N = length(u);

uhat = fft(u);
kappa = fftshift((-N/2:N/2-1)');
DDuhat = -pi^2*kappa.^2.*uhat;

DDu = ifft(DDuhat);

duout = DDu + (1-u.^2).*u;

end
</code>```

I imagine I'm doing something wrong in using the cosine series (maybe I can't use FFT). Any help would be appreciated!

Title: Pseudo-Spectral cosine transform


问:
Vectorizing Matrix Multiplication
I would like to do the following operation: I have a "4D" matrix A and a "3D" matrix B. Both A and B are actually 2D matrices, where for A, each element is a 2D matrix, and for B, each element is a 1D matrix or vector. I would like to, for each element in the 2D floor plan, multiply the matrix on A and the matrix in B. Graphically,

Currently, my code is:

```<code>% size(A) = Nx, Ny, n, n
% size(B) = Nx, Ny, n
% size(C) = Nx, Ny, n

for k1 = 1 : Nx
for k2 = 1 : Ny
    C(k1,k2,:) = squeeze(A(k1,k2,:,:)) * squeeze(B(k1,k2,:));
end
end
</code>```

I am looking for a way to "vectorize" the code and avoid loops completely. I am hoping such a vectorization will speed up my program.

Title: Vectorizing Matrix Multiplication


问:
Spherical Gradient along $z$-axis
Let us consider a function $f(r,\theta,\phi)$ that can only be efficiently computed in spherical coordinates. We need its Cartesian gradient. We can implement its analytical spherical derivatives, compute the spherical gradient, $\left(\frac{\partial f}{\partial r},\frac{1}{r}\frac{\partial f}{\partial \theta},\frac{1}{r\sin\theta}\frac{\partial f}{\partial \phi}\right)$, and convert the gradient to Cartesian component.

What would you do with the values along the $z$-axis, where $\theta$ is undefined and usually set to 0, and in the origin, where $r$ is 0 as well? More specifically, what would you do with the terms in the denominators $r$ and $r\sin\theta$?

Title: Spherical Gradient along $z$-axis


问:
Find B-spline coefficients from values on collocation points
The context of my question is how to compute high order derivates on direct numerical simulation of turbulent channel flow. It is of particular interest for fluid dynamics and turbulence research.

I want to find the coefficients of a B-spline curve given the value of the resulting function on a set of collocatioin points. Is there any way to do it without resorting to solve a linear system each time I need to find the coefficients.

My specific case is to interpolate using a 7-th order B-spline scheme with knots given, but the values are set on 'Greville Abscisae'.

The knots are given in this address:
    http://turbulence.pha.jhu.edu/docs/channel/y-knots.txt

And the collocation ponits in this:
    http://turbulence.pha.jhu.edu/docs/channel/y.txt

Though not necessary, you can query any example of such coefficients on this website, just inputing the collocation points on this website:

http://turbulence.pha.jhu.edu/webquery/query.aspx

The quantity I am seeking to interpolate is the velocity field for the channel dataset.

I have never worked with b-splines before, so I might be missing something, but I haven't found anywere a way to compute the coefficients from the collocation points, and this is bothering me.

Solve a linear system for each line (in this case, a 512 x 506) would be bothersome, since I would need to do this for each vertical line on each snapshot of the simulation, and I would like to avoid this as much as possible, but I am not finding any other alternative.


问:
1-D local interpolation of sampled function that is polynomial in 1/x
I have a number of evenly-spaced samples of functions that are polynomial in $1/x$ where $x$ is a continuous variable on $x\in\left[a,b\right]$, i.e. $$f(x_i)=\sum_{n=0}^{N} f_n x_i^{-n}.$$ I want to use some sort of local interpolation function like a hat function (linear interpolation), but it is pretty inaccurate (about 10% relative error in the $L_2$ norm) for values of $x$ near $a$. Near $b$, the interpolation is much better. I haven't yet implemented higher-order local Lagrange interpolation polynomials. I was wondering if anyone knows of a local interpolation method better suited to this problem. For reference, $N$ may be anywhere from about 5 up to 1000, and the number of samples reflects this. I'm not completely wedded to evenly-spaced samples, either, if a good interpolation method requires a different spacing. My thinking is that perhaps concentrating sampling points closer to $a$ and interpolating on a uniformly-spaced transform grid might work.</p>


问:
How would you specify mixed boundary conditions for a 2D PDE in the matrix used for finite differences
I have the following PDE in 2D:

$U_{x} + U_{xx} + U_y + U_{yy} + U_{xy} = f$

where $f$ is a constant.

And I'm trying to create a matrix $A$ to solve the PDE through finite differences: $AU = f$. I will invert $A$ such that $U = A^{-1}f$.

I have mixed boundary conditions in this form:

$0 = a_0U + a_1U_{x} + a_2U_{xx} + a_3U_y + a_4U_{yy} + a_5U_{xy}$

I'm just not sure how I could specify or code this BC equation into my matrix. I did a lot of searching but most of the tutorials are for 1d equations only. I know about the ghost point method but I have only done it for $0 = a_0U + a_1U_{x}$. However, I'm not sure what I should do when it comes to the point at $f(x_1, y_1)$, since it would involve multiple ghost points ($f(x_0,y_1), f(x_1,y_0)$, etc.).

Separately, should $a_{0, ..., 5}$ be a function of $x$ or $y$?


问:
Object-oriented non-linear solving in python
I'd like to build a system in Python, consisting of (broadly speaking) objects which are internally described with (not necessarily just linear) equations, that I can connect with each other - similar to what can be done in OpenModelica. Is there anything like this already available or do I have to use OpenModelica for this?

For now, a steady-state solution is just fine. That's also why I think that OpenModelica would be overkill.


问:
Iterative solver, which balances the equations with the biggest errors
I am trying to understand the mathematical side of an algorithm, which approximates the solution of coulomb's law. Their approach is to balance the two elements with the biggest positive/negative error.

Trying to understand the stability conditions of the algorithm, I am searching if there is an iterative solver with an similiar approach, as the problem can be stated as a linear system (but has a lot more computational complexity when solving the system with usual solvers).

The algorithm does not use any derivatives, but reduces the error induced from two two elements with the biggest residuum on each other to zero and recalculates the residuum of all other elements after the balancing.

The potential they are balancing is defined by $U_i = \sum_j I_{ij} q_j$ at points $x_i$ at the center of triangles $\Delta S_j$ with charge $q_j$, where
$$I_{ij} = \frac{1}{\Delta S_j} \int_{\Delta S_j}k|\vec{x}_i-\vec{x}'|^{-1}dS'.$$

In each step the points m, n with the biggest positive/negative deviation to an equipotential are balanced by transferting a charge $$q=\frac{U_m - U_n}{I_{mm}+I_{nn}-I_{mn}-I_{nm}}$$ giving the new potentials at the points $x_i, x_j$:

$$
U'_m = U_m - I_{mm}q + I_{mn}q \\
U'_n = U_n + I_{mm}q - I_{mn}q
$$

After this transfer $U'_m = U'_n$ holds and all other potentials $U_i$ are updated:

$$U'_i = U_i + q_m(U'_m - U_m) + q_m(U'_n - U_n)$$

After each step, the potentials $U'_n, U'_m$ are equal and in the setting of electrostatics all other potentials come closer to an equipotential on the surface as well.

I suspect, that this way of optimization is already described somewhere else with more numerical analysis about errors, convergence criteria and stability conditions, but was not able to find it in literature, yet.


问:
MPI_Finalize() only some processors
I want to develop a distributed algorithm with MPI which have load balancing capability. But there is a situation that I'm not very clear about, which is, when the number of tasks is less than the number of processors available, I need to detect that and terminate those processors whose ranks (actually rank + 1) are greater than the total number of tasks. In my implementation, I just add a condition to determine that, if the condition is true, do something, otherwise, call MPI_Finalize() to shut down these processors. For your reference, here is the psuedocode:

```<code>MPI_Init();
MPI_Comm_size(MPI_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_WORLD, &proc_rank);
if (proc_rank + 1 < task_counts):
    do something (including some collective communications);
MPI_Finalize();
return 0;
</code>```

This algorithm seems to be quiter obvious. However, it doesn't work. I added some fprintf code after MPI_Finalize(), but this code never got there for those processors whose ranks are out of range.

So what's the problem here?

Title: MPI_Finalize() only some processors


问:
Crank-Nicolson scheme in space for advection equation
Consider the equation

$$\frac{\partial}{\partial t}v(t,x)=\frac{\partial}{\partial x}v(t,x)$$,

for $t,x\in\mathbb{R}$.

I'd like to solve this equation forward in space and backward in time, updating in space given the initial condition in space: $v(t,0)=f(t)$.

Now, we can discretise in time and take the backward difference of the time derivative, which yields

$$\frac{\mathrm{d}}{\mathrm{d}x}v^n(x)=\frac{v^n(x)-v^{n-1}(x)}{\Delta t}\qquad(1)$$

Then the Crank-Nicolson scheme gives

$$\frac{v^n_{i+1}-v^n_i}{\Delta x}=\frac{1}{2}\left[G^n(x_{i+1})+G^n(x_{i})\right],\qquad(2)$$
where $G^n(x)$ is the RHS of equation $(1)$.

Then $(2)$ becomes

$$\frac{v^n_{i+1}-v^n_i}{\Delta x}=\frac{1}{2}\frac{v^n_{i+1}-v^{n-1}_{i+1}}{\Delta t}+\frac{1}{2}\frac{v^n_i-v^{n-1}_i}{\Delta t}$$
i.e.
$$\left(\frac{1}{\Delta x}-\frac{1}{2\Delta t}\right)v^n_{i+1}+\frac{1}{2\Delta t}v^{n-1}_{i+1}=\left(\frac{1}{\Delta x}-\frac{1}{2\Delta t}\right)v^n_{i}+\frac{1}{2\Delta t}v^{n-1}_{i}$$
so that
$$\mathrm{A}\vec{v}_{i+1}+\frac{1}{2\Delta t}\mathrm{B}\vec{v}_{i+1}=\mathrm{A}\vec{v}_{i}+\frac{1}{2\Delta t}\mathrm{B}\vec{v}_{i}$$
for $\mathrm{A},\mathrm{B}$ matrices yet to be specified.

However, this implies that

$$\vec{v}_{i+1}=\vec{v}_i,$$

which would appear to be nonsense. Have I done something wrong?


问:
CFL condition of source term
Consider one dimensional hyperbolic pde
$$u_t+f'(u)u_x=0$$

For the above problem ,CFL condition is $\Delta t\leq \dfrac{\Delta x}{|f'(u)|}.$

But if we include  the source term $,S(u,t),$ which depend upon $u$  and $t$ in above equation i.e.
$$u_t+f'(u)u_x=S(u,t)$$
 What is the influence of $S(u, t)$ in the CFL condition ?.

Title: CFL condition of source term


问:
How to make LAPACK eigenvectors orthogonal like Matlab?
I'm using LAPACK zgeev to calculate eigenvectors of a symmetric complex matrix of high dimensions ($n \approx 2000$). I need these eigenvectors to satisfy

$$\sum_{m=1}^{n} \frac{\lvert m \rangle \langle \overline m \lvert}{\langle \overline m\lvert m \rangle} = I_n$$

Where $\lvert m\rangle$ is an eigenvector and $\langle \overline m\lvert$ is the same eigenvector transposed (Attention: Not conjugated) and $I_n$ is the identity matrix of size $n$.

I don't know if all symmetric complex matrices satisfy this but this matrix has specific structure and it always should:

$n=6$

(Algorithm to generate these matrices are way too complicated to share, But I think one can get their pattern from these two examples)

$n=9$

But the eigenvectors I get from zgeev do not follow the equation I mentioned and they are not orthogonal.

I have tested these matrices in Matlab and eigenvectors generated by Matlab satisfy this condition AND are orthogonal. (i.e $\langle m_i\lvert \overline m_j \rangle=\delta_{ij}$)

How can I get the same eigenvectors from LAPACK zgeev?


问:
BTCS-like method for heat conduction in unstructured triangular grid
I want to write a simple simulation for heat conduction in a unstructured triangular mesh.

I already made it work for a structured rectangular grid with the ADI method, but now I need more complex geometries.

With $d_x=\frac{\kappa \Delta t}{(\Delta x)^2}$ and $d_y=\frac{\kappa \Delta t}{(\Delta y)^2}$ you can say:

$$d_x T_{i+1,j}^{n+1}+d_x T_{i-1,j}^{n+1}-(2d_x+2d_y+1)T_{i,j}^{n+1}+d_y T_{i,j-1}^{n+1}+d_y T_{i,j+1}^{n+1}=-T_{i,j}^n\ ,$$

and you basically have a penta-diagonal equations system.

For the triangular grid, the only thing I came up with till now is to calculate the heat-flux-density over every edge, sum them up for every cell and add them to the temperature of the cell. For me this seems to be a explicit method, which (if it follows the same behaviour than the FTCS approximation) should be pretty instable.

Thats why I want to come up with an implicit method, but till now, I'm not able to.

Has anybody some advice or some estimation on the stability of my explained method?


问:
Color legend location in eps files, exported from paraview v5.4
I am getting this strange error with paraview 5.4.
For example, if I load one of the example visualisations:

and if I export it as eps or pdf, then the color legend is place in the bottom left corner:

This does not happen for jpg or png files.
Anybody know how to fix this?
Thanks


问:
Understanding Davis artificial viscosity
I'm solving 2D Euler's equations in Cartesian coordinates:
$$
\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = 0, \qquad
U = \left( \begin{array}{c} \rho \\ \rho u \\ \rho v \\ e \end{array} \right), \quad
F = \left( \begin{array}{c} \rho u \\ p+\rho u^2 \\ \rho uv \\ (e+p)u \end{array} \right), \quad
G = \left( \begin{array}{c} \rho v \\ \rho uv \\ p+\rho v^2 \\ (e+p)v \end{array} \right)
$$

using finite difference MacCormack method with Davis artificial viscosity which dampens non-physical oscillations. The MacCormack method is pretty simple and not shown here. Davis viscosity works good in my simulations, but is quite difficult; I'm trying to fully understand it in order to apply it cylindrical coordinates later. I hope that some of you may help me with unclear parts of it (shown in bold font below).

Viscous terms are added after 2nd MacCormack step as follows (formulated in 1D for simplicity):
$$
U_j^{n+1}=U_j^{(m)} + \left( D_{j+\frac{1}{2}}-D_{j-\frac{1}{2}} \right),
$$
where $U_j^{(m)}$ denote field values (at $j$-th coordinate point and $(n+1)$-th time layer) calculated by MacCormack scheme as usual, and $D_{j+\frac{1}{2}}, D_{j-\frac{1}{2}}$, are forward and backward differences of U with non-linear coefficients:
$$
D_{j+\frac{1}{2}} := \frac{1}{2} C(\nu_j) \left( 2 - \phi(r_j^+) - \phi(r_{j+1}^-) \right) \Delta U_{j+\frac{1}{2}}, \quad \Delta U_{j+\frac{1}{2}} = U_{j+1} - U_j \\
D_{j-\frac{1}{2}} := \frac{1}{2} C(\nu_{j-1}) \left( 2 - \phi(r_{j-1}^+) - \phi(r_{j}^-) \right) \Delta U_{j-\frac{1}{2}}, \quad \Delta U_{j-\frac{1}{2}} = U_{j} - U_{j-1}
$$

Where $C(\nu)$ is a limited coefficient based on Courant number:
$$
C(\nu) := \begin{cases} \nu(1-\nu) & \;\; \text{if} \;\; \nu \le 0.5 \\ 0.25 & \;\; \text{otherwise} \end{cases}, \qquad \nu = \rho(A(U_j)) \frac{\Delta t}{\Delta x},
$$

$A(U)$ is the Jacobian matrix of vector $F$, $\rho(A(U))$ is it's spectral radius (for Euler's equations it equals $\text{max}(|u|+a, |u|-a)$, where $a$ is the speed of sound).

$\phi(r)$ is simple limiting function: $\begin{cases} 0, & r < 0, \\ 2r, & 0 \le r \le 0.5, \\ 1, & r >0.5\\ \end{cases}$

And finally $r_j^+$,$r_j^-$ are slope ratios:
$$
r_j^+ := \frac{\left( \Delta U_{j-\frac{1}{2}}, \Delta U_{j+\frac{1}{2}}\right)}{\left( \Delta U_{j+\frac{1}{2}}, \Delta U_{j+\frac{1}{2}}\right)}, \qquad
r_j^- := \frac{\left( \Delta U_{j-\frac{1}{2}}, \Delta U_{j+\frac{1}{2}}\right)}{\left( \Delta U_{j-\frac{1}{2}}, \Delta U_{j-\frac{1}{2}}\right)}
$$
Where $\left( \cdot , \cdot \right)$ is scalar product. These ratios are positive in monotone areas and negative in non-monotone areas, so that in general coefficients like $\left( 2 - \phi(r_j^+) - \phi(r_{j+1}^-) \right)$ are equal 2 in non-monotone areas; about 1 in monotone, but "curvy" areas; and equal 0 in monotone (and rectilinear) areas.

If both these ($2-\phi-\phi$) coefficients are equal 2, and both $C(\nu)$ are equal 0.25 we have an effective second derivative of the field value:
$$
D_{j+\frac{1}{2}}-D_{j-\frac{1}{2}} = 0.25 \left( U_{j+1} - 2U_j + U_{j-1} \right).
$$

But where is the division by $\Delta x^2$? If the viscosity would be simply the 2nd derivative, the scheme would look like
$$
U_j^{n+1}=U_j^{(m)} + \frac{\Delta t}{\Delta x^2} \left( U_{j+1} - 2U_j + U_{j-1} \right).
$$
Is $C(\nu)$ effectively equals $\frac{\Delta t}{\Delta x^2}$?

And, what if $D_{j+\frac{1}{2}}$ is non-zero and $D_{j-\frac{1}{2}}$ is zero? Does it still yield an effective 2nd derivative?


问:
adaptive / smart grid choice for dynamic programming
I've a read a paper about investing where they used a dynamic programming approach to solve a finite horizon problem, i.e.

$$\max_{x_t} E[u(W_T)] $$

where $u$ is a utility function and $W_T$ denotes the terminal wealth. The Bellman equation can then be written as

$$ V_t(W_t) = \max_{x_t}E[V_{t+1}(W_{t+1})|W_t]$$
with $V_T=u$. I have a continuous and concave utilitiy function which penalizes not hitting a certain target amount

$$  u(w) =
  \begin{cases}
   K_1 & \text{if } w \geq K_1 \\
   w - c\max{(0,(K_2-w))}^2       & \text{if } w < K_1
  \end{cases}$$
where $K_1>K_2$ are constant and $c>0$ a scaling factor. Above $K_1$, $u$ is constant, below $K_2$ you add a penalty. $c$ is a control parameter how much we should penalize $u$. For $K_1 = 12$, $K_2=10$ and $c= 100$ we get the following picture.

and to see the linearity of $u$ between $K_2$ and $K_1$ another picture:

Since I'm using a backward induction algorithm and discretize the space I've naively tried a uniform space grid $\{x_i\}$ in a fairly large interval $[w_1,w_2]$. Afterwards, I've solved for all the grid points $\{x_i\}$ the problem leading to solutions $y_i$. The pair $(x_i, y_i)$ was then used to approximate $V_t$. For this I've used a shape preserving method. What I've seen from experimenting with the problem that the approximation is very good far out on the left and right side, away from high curvature region. Moreover, the solution becomes quite stationary in these "wing regions", too. What I then thought is to generate a non uniform grid which is denser where curvature is high. The problem is I don't know a priori the shape of $V_t$ (expect for $t=T$) so that it's hard to come up with a rule. Is there a generic way how to construct such non uniform grids? If someone has a concrete example or could point me into the right direction to do my own research, that would be really appreciated. I will potentially add another dimension to the state space so that the same question would arise in a higher dimension. That's why I put the mesh-generation tag as well


问:
How to choose a stable PML for pseudo-spectral method with strongly varying velocity
My friend was working on this, and he asked me about the stability of PML while applying on pseudo-spectral method, I believe his concern was how to choose the difference(if the difference should be chosen with a directional preference, forward/backward, or inflow/outflow) such that the numerical solution is stable over time. I am nowhere an expert on this topic and do not even know how to choose the keywords to google, so here I am requesting some references. Thanks a lot.</p>


问:
Computing dilogarithm
I'm measuring the integral of a quantity which, mathematically, requires the computation of a dilogarithm function.

$$\operatorname{Li}_2(be^{ax})$$

where $b$ and $a$ (are real and) can be positive or negative. Dilogarithm has very nice properties that enable us to compute its value efficiently, For $|z|<1$, the following series converges absolutely:

$$\operatorname{Li}_2(z)=\sum_{n=1}^{\infty}\frac{z^n}{n^2}$$

I'm wondering is there any explicit way to compute or numerically approximate dilogarithm -- especially in COMSOL?

Title: Computing dilogarithm


问:
Discrepancy in estimating boundary stencil for finite difference method
I am trying to estimate the FD stencil for boundary as mentioned in this paper (section 4.1)  using MATLAB. The stencil order (6th) is higher than the one mentioned in paper (4th).
$$ f_1' +\alpha f_2' = \frac{1}{h} (c_1f_1 + c_2f_2 + c_3f_3 +c_4f_4+ c_5f_5 + c_6 f_6 + c_7 f_7 )$$
Each term of both sides is expanded using Taylor series expansion and then derivatives are matched to estimate the coefficients [$c1,c2,c3,c4,c5,c6,c7$ and $\alpha$].

I solved the above system as following:

Above system of equation can be written as M*x=y

Note:

1) The alpha was taken to the left to form a part of unknown vector x

2) the factorial on both side were simplified and adjusted to right side.

I wrote following code for solving the above system of equation.

```<code>LM = zeros(8);
LM(1,1:7)=1;                LM(1,8)=0;
LM(2,2:7)=1:6;              LM(2,8)=-1;
LM(3,2:7)=(1:6).^2;         LM(3,8)=-2;
LM(4,2:7)=(1:6).^3;         LM(4,8)=-3;
LM(5,2:7)=(1:6).^4;         LM(5,8)=-4;
LM(6,2:7)=(1:6).^5;         LM(6,8)=-5;
LM(7,2:7)=(1:6).^6;         LM(7,8)=-6;
LM(8,2:7)=(1:6).^7;         LM(8,8)=-7;
y=[0;1;0;0;0;0;0;0];
x=LM\y;
rats(x)
</code>```

which gives completely different answer than reported one. Both are given below:

Solution i.e. $x =[c1,c2,c3,c4,c5,c6,c7, \alpha$] by

MATLAB estimated solution:     [-69/20, -17/10, 15/2, -10/3, 5/4, -3/10, 1/30, 6]

Solution reported in literature: [-197/60, -5/12, 5, -5/3, 5/12, -1/20, 0, 5]

Any help is appreciated.


问:
Extrapolating to non-fluid cells (for Shallow Water Equations), for a shore/beach?

The water height $h$ and 2d velocity field $(u,w)$ are "extrapolated
  to non-fluid-cells, i.e., setting $h$ equal to the value in the
  nearest fluid cell." [Bridson]

I'm using finite differences.

Looking at the course notes, this extrapolation seems involved: tracking the zero  level set surface $\phi$ (using signed distance from the surface), then taking its gradient $\nabla \phi$.
It seems to be done for all non-fluid cells, not just those adjacent to fluid.

I don't understand why this done (BTW for the 2d SWE case, if the non-fluid cell is lower e.g. the lip of a dam wall, then extrapolating the water height will make large quantities of water appear from nowhere (because the height is so far above the ground).

Is there a simple explanation, and a simpler way to do this?

EDIT Two different quantities are being extrapolated, velocity and height. The velocity extrapolation seems to be needed for semi-Lagranian advection. This is where you advect a quantity $q$ by projecting velocities backwards over a timestep, to find what value of $q$ was at that location (interpolating between the grid locations where $q$ is stored). The issue is that this back-projection can land in a region where there is no $q$... so, you extrapolate the value from the nearest cell that does have $q$ (BTW which is along the normal to the surface).

EDIT The other extrapolation, of height $h$, seems to only be used for pressure acceleration. (Its effect seems to be a Neumann BC, zero derivative normal to surface, by making the other side the same as this one...?)  Height $h$ is derived from depth $d$ and bottom $b$, $h=b+d$, and $h$ isn't stored, and doesn't seem to be used again. Although, logically, when $d$ is advected semi-Lagrangianly, it should need extrapolation, the same for any other $q$...

Any good references or articles with a nice precis/lit review, or survey papers?

[Bridson] I'm implementing the SWE from chapter 12 of Bridson's book Fluid Simulation for Computer Animation.


问:
Solve a PDE BVP with Spectral methods in time and space?
I have a PDE (coming from a Bellman equation with $z$ under brownian motion).  Let $z \in [0,\infty)$ and $t \in [0,T]$.   To sketch the equation:
$$
(r - g(t)) v(t,z) = \pi(t,z) + (\gamma - g(t))v_z(t,z) + \frac{\sigma^2}{2}v_{zz}(t,z) + v_t(t,z)
$$
with a terminal condition $v(T,z)$ and ugly nonlinear boundary conditions involving integrals for $v(t,0)$ and $v_z(t,0)$.  There is a very good chance that the PDE will be stiff and $v(t,z)$ fairly close to $v(T,z)$ for all $T$

I am considering two alternatives:

Use a 2-dimensional Chebyshev basis
$$v(t,z)\approx \sum_{i=0}^{N_z}\sum_{j=0}^{N_t}a_{ij}T_i(z)T_j(t)$$
for some $N_z$ and $N_t$ and $T(\cdot)$ are the chebyshev polynomials adapted to the domain.  Then solve the BVP using spectral collocation methods for the coefficients (i.e., the tensor product of the polynomial rather than any sparse grid method)

Use finite-differences in time, and a 1-dimensional chebyshev in space.  Hence, for each $t_j$ in the time grid,
$$v(t_j,z) \approx \sum_{i=0}^N a^{j}_i T_i(z) $$
and solve for the set of $a^j_i$ for each $j$ in the finite difference scheme, and the set of $N$ nodes at each period.

My question is:
If I don't care all that much for extreme precision, but care a great deal about convergence and stability, is either case reasonable or is there a good reason not to use a 2-dimensional Chebyshev basis?

Without going into details, the $\pi(t,z)$ and $g(t)$ come from an outside set of equilibrium conditions that need to be solved, and it would be very convenient to be able to approximate them as chebyshev polynomials themselves (and solve for all of the approximating coefficients for the $g(t)$ and the above $a_{ij}$ at the same time using auto-differentiation and a good nonlinear solver).

(I found this, Spectral Methods in time but the exact set of problems didn't seem clear from the answer.  The typical advantages of a finite-difference time stepping are less clear in my case).


问:
Optimizing an error function involving rotation vectors in Python
I'm prototyping a system that finds the 3D pose of a object in a video sequence. For this I minimize a error function involving the rotation and translation of the object as parameters and two sets of points (polluted with gaussian noise but no outliers) as data.

I have tried to parameterize the rotation matrix both using 3 euler angles, and a rotation vector. Then I test my system with a synthetic data sequence and I see I obtain substantially less error with the 3 euler angles parameterization that with the rotation vector. I was expecting the opposite (notice the sequence I use doesn't result in any gimbal locks for the euler angles).  This makes me wonder if I'm using the rotation vector parameterization correctly.

The steps I follow are, first I convert the rotation matrix to rotation vector (using the Rodrigues formula available in OpenCV) the resulting 3-vector is part of the parameters of my error function. Inside the error function I convert the rotation vector to rotation matrix and compute the residuals. I'm using python and scipy.optimize.leastsq (which internally uses MINPACK.lmdif) to minimize this error function. leastsq  will compute the jacobian needed for the minimization numerically from the error function.

After reading a bit of theory about the lie group SO(3) and the lie algebra so(3), I understand that the so(3) (for me the rotation vector) corresponds to a tangent plane at the identity of SO(3). So my rotation vector is a local approximation to the corresponding rotation matrix.

What I thought that maybe I'm missing in my minimization is to convert the rotation vector back to a rotation matrix after each minimization step, and then back to a rotation vector (to take a new local approximation to the rotation matrix) for the next step. I can't test this possibility as I'm using a standard minimizer, which doesn't allow to do this.

Now I think that the Jacobian calculated at each minimization step should be enough to direct the rotation vector parameters in the correct direction without the conversions I mention in the above paragraph. However, I'm not fully sure and as I get worse performance with the rotation vector parameterization wrt. the euler angle parameterization, I assume I'm doing something wrong.

Would someone clarify/confirm the correct way of using the rotation vector in a minimization problem?


问:
3D visualization of the results of a 2D simulation in COMSOL
Considering that COMSOL forum is dead and CFD-online is not responsive, I'm gonna try asking my question here:

I have done a 2D simulation in COMSOL multiphysics. I have two independent variables x and t and two dependent variable u1 and u2. Comsol gives me 2D projections of the 3D graphs u1(x,t) and u2(x,t) as below:

What I need is a 3D plot like this:

I have tried exporting the data. There are two options in COMSOL. Text which generates a file with several columns like this

The other options is Unstructured VTKFile (*.vtk). I tried importing both into ParaView but I could not manage to get any plots. Here and here I uploaded the exported data and here is my COMSOL file.

I would appreciate if you could help me know how to do this in COMSOL, ParaView or any other possible software.


问:
Implementing Neumann boundary condition in nonlinear integro-differential equation
Problem

I would like to solve a nonlinear integro-differential equations on a 2D square domain $\Omega$ subject to Neumann boundary conditions using finite differences:
$$\frac{\partial u}{\partial t} = f(u) + \int_\Omega G(x'-x)\left[u(x')-u(x)\right] dx'.$$

The variable $u(x,t)$ depends on space and time, $f$ is a nonlinear function and $G$ is the integral kernel. Neumann (also called no-flux or reflective) boundary conditions are given as
$$ \frac {\partial u}{\partial x}\Bigr|_{\partial\Omega} = 0$$
at the boundary $\partial\Omega$.

Solution Attempt

Usually to implement Neumann boundary conditions, it is either possible to use ghost cells or modify the stencil. Let's try the latter. Starting from the heat (or diffusion) equation in 1d, we know the $2^{nd}$ derivative can be discretized as follows:
$$ \Delta u \rightarrow \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}=\frac{1}{h^2}\left((u_{i+1}-u_i)+(u_{i-1}-u_i)\right),$$

where $u_i$ is shorthand notation for $u(x_i)$ and $h$ is the uniform distance between neighboring points.

With $G=1/h^2$, this formula can be rephrased as a discretized version of the integral equation from the top:
$$ \sum_{j=i-1}^{i+1} G(x_j-x_i)\left[u_j - u_i\right].$$

Differences that involve points $x_j$ outside the domain are zero, due to the Neumann boundary condition:
$$ \frac {\partial u}{\partial x}\Bigr|_{\partial\Omega} \rightarrow \frac{u_j-u_i}{h} \overset{!}{=} 0.$$

Pseudo C code might look like this:

```<code>for(j=max(0,i-n); j<min(n-1,i+n); j++){
   uIntegro[i] += G[j-i]*u[j];
}
</code>```

Question

Coming back to the original problem, given that the integral kernel ranges beyond the first neighbor, let's say it takes into account the first 5 neighbors, is it correct to just set all differences in the integral involving outside points to zero in order to realize Neumann boundary conditions (in 2d as well)? Bonus: If this is not Neumann, what would it be called instead?


问:
How to construct a eigensolver targeting a specific type of matrix
I need to diagonalize such kind of matrix during research: it's n-by-n, with it's upper-left (n-1)-by-(n-1) corner be diagonal while the nth row & column are dense. It's observed that the self adjoint eigensolver is not quite efficient as a solution when n is greater than roughly 2000 and there is no existing eigensolver to take advantage of the sparseness. A possible way I just came up with is to do the LU via the corresponding sparse package first and then attack the two triangular matrices. However I'd love to hear more suggestions & insight. </p>