问:
Fortran compiles for legacy Finite Element Fortran program (1980)
The version of Fortran used come from  Montreal Ecole Polytechical in 1980.
I need a compiler for Fortran for Windows 7 or Windows 8.</p>


问:
Minimize MAE for sum of powers
The model is
$$y_i=\beta\big(x_{i1}^\alpha+x_{i2}^\alpha+...+x_{im_i}^\alpha\big)+\epsilon_i\textrm{ for }i=1, 2, ..., n$$

 I want to minimize the MAE, i.e.
$$\Sigma_{i=1}^n{\big|y_i-\beta\big(x_{i1}^\alpha+x_{i2}^\alpha+...+x_{im_i}^\alpha\big)\big|}$$

 I tried gradient descent, but I'm worried that the solution is not a global minimum.

 Could anyone help me?

ContentLicense: CC BY-SA 4.0

Title: Minimize MAE for sum of powers


问:
Simulating a combustion process
I want to try simulation-(and not experimental)-driven approach to design custom fireplace fuel burners.

What software applications, libraries, code and model templates can I use to model and visualise a combustion process on a 3D object?

I want to design a 3D model of a static solid object, assign properties to it (e.g., choose material to be steel or glass), put it into an environment (Earth's gravity and atmosphere), and simulate the combustion process by controlling the supply of some fuel (liquid or gaseous, e.g., propane or ethanol). The model can be as simple as a cup with a fuel intake or it can be something more advanced, with several chambers, intakes, and vents.

I do not need to bother about the structural integrity of the object, its thermal expansion or other too-complicated combustion effects, like carbon deposition, etc. Maybe, thermodynamics of the object itself can be ignored. The conditions are not meant to be extreme, something that ranges between still liquid fuel burning and a blowtorch at a small exhaust rate.

However, I want to get a more-or-less accurate approximation of the combustion of the fuel that's exhausted by the object. After generating a spark at a certain point, the byproducts of the "burning fuel and oxygen" should influence the aerodynamics of the combustion process.
The whole process has to be visualised. The simulation must capture the dynamics of the combustion, e.g. whether the flame is going to be stable, will oscillate, fade away, propagate to a nearby vent, etc.

Ideally, the software should be graphical and intuitive, capable of realtime simulation and visualisation on a consumer desktop. But, (let's be realistic), more demanding options (computationally and setup-wise) are welcome as well.

I have some experience with Octave/Matlab, Mathematica, and a bunch of C/C++/Python scientific libraries. These are probably the main candidates, but other products (with the cost of up to ~1000$) might also be considered.

Title: Simulating a combustion process


问:
Model and solve nutrition optimization problem: how to?
How to solve/assess the following problem?

Given: $N$ ingredients like apples, bread etc. Mass of an ingredient $j$ is in this simplified model a sum of macro nutrients carbohydrates, proteins and fats in proportions $c_j$, $p_j$ and $f_j$:

$(1)\ M_j = (c_j + p_j + f_j) * x_j $ (in grams)

$(2) \ c_j + p_j + f_j = 1$

Also a nutrition model is provided,

$ (3) A_c + B_p + C_p = 1$.

Problem: find optimal quantities $x_j$ of ingredients so that total mass possibly fulfils the provided nutrition model, i.e.

$ (4) \ M_{total} = x1 * M1 + x_j + M_j * x_j .. + M_N * x_N $

where

$$ (5) \sum_{j=1}^n c_j / M_{total} = A_c$$
$$ (6) \sum_{j=1}^n p_j / M_{total} = B_p$$
$$ (7) \sum_{j=1}^n f_j / M_{total} = C_f$$

From my current knowledge state, I assume this is a linear programming problem which can be approached for example using a solver like Python's PuLP. So I will try to get my head around it but if I'm going in completely wrong direction I hope an answer here will point me to another more suitable method.

Example:

Say we want to consume 1500 kcal (being on a diet, typically adults eat about 2000 plus something kcal which is on a smaller detail level an individual value based on sex, age, weight and physical activity level aka PAL)

A possible recommendation is to have these calories as a "budget" distributed between 55% carbohydrates, 15% proteins and 30% fat. All that keto and low-carb stuff is basically how these three parameters are weighted.

Now we know we should eat accrodingly the "calorical vector" [825.0000000000001, 225.0, 450.0] or recomputed to mass in grams using reference values [201, 55, 48]

A random menu more or less close to that could be then for example like: "1489 kcal = 25g butter, 85g Emmenthal, 180g turkey breast, 800 g banana"

Modeling approach 1

Decision variables: mass values for macro-nutrients; c_banana for total carb value in optimal banana quantity.

Constraints: need to tell the program what separate nutrients contain and what's diet goal

Optimization objective: while eating that less calories, let's optimize for max mass of this healthy food.

Encountered problems

How do I actually model that banana contains 21% carbohydrates? The naive guess,

```<code>model += (c_banana /(c_banana + p_banana + f_banana) == 0.21, "banana_c")
</code>```

leads to the following error:

```<code>TypeError: unsupported operand type(s) for /: 'LpVariable' and 'LpAffineExpression'
</code>```


问:
Understanding Illumination Optimisation Problem
I am a newbie to convex optimisation and I am learning with the aid of CVXPY. I am requesting for clarity on the illumination problem as described in Boyd & Vandenberghe lecture 1 slides here. I have adapted the code from here and tried to visualise what the illuminance from within a room would look like. I have inserted my code and I hope it's not too cumbersome to look at.

I am using four lamps and I have two detectors (patches) labeled Rx0 and Rx1 as seen on the attached image. I desire to have an illuminance of 450 lux at both of these detectors but in my many attempts, none of the detectors is at 450 lux. Rx0 is around 380 lux while Rx1 is around 510. I only get somewhat the exact desired value for one of the detectors if I place one of the detectors at the centre of the room.

Is the objective of this problem to have equal illuminance at all patches or a particular illuminance at just one of the patches? I am not fully understanding here, my thinking is that we are aiming to have equal (desired) illuminance at all patches using the minimum possible power.

I am curious to know if the solutions from CVXPY are approximate or my code is wrong.

Can I formulate the problem so that I end up having two varying illuminance levels at each patch? As in having like two beams from one of the transmitters to two patches? Something more like beamforming? If so, what would be the best approach?

With the definition of illuminance on surface $\mathsf{s}$ as

\begin{equation}
\mathcal{E}_{\textsf{s}}^{\textsf{tx}} = \underbrace{\frac{\mathcal{I}_{N_{l,u}}\cos^{m_{1}}\left(\phi\right)}{d_{l,u}^{2}\cos\left(\psi\right)}}_{\mathcal{E}_{\textsf{TX}}} +  \underbrace{\sum_{j\in \mathcal{L}\backslash\lbrace l\rbrace}\frac{\mathcal{I}_{N_{j,u}}\cos^{m_{1}}\left(\phi\right)}{d_{j,u}^{2}\cos\left(\psi\right)}}_{\mathcal{E}_{\textsf{nTX}}},
\end{equation}

where $\mathcal{I}_{N_{l,u}}$ is the luminous intensity (normal to surface) of one of the four LED arrays, $\phi$ is the irradiance angle, $\psi$ is the incidence angle, $d$ is the distance between the lamp and the surface/patch, $\mathcal{I}_{N_{j,u}}$ is the luminous intensity (normal to surface) of the other three LED arrays, $\mathcal{E}_\textsf{s}^{\textsf{tx}}$ is the illuminance due to all four LED arrays and $\mathcal{E}_{\textsf{r}_u}$ is the required illuminance at patch $u$.

I have formulated my problem as follows;

\begin{align}
\underset{\mathcal{I}_{N_{l,u}}> 0,~\mathcal{I}_{N_{j,u}}\geq 0}{\mathtt{minimise}} \quad& \mathtt{max}_{\textsf{s}}f\left(\mathcal{E}_\textsf{s}^{\textsf{tx}}/\mathcal{E}_{\textsf{r}_{u}}\right), &\label{equ:objectiv2}\\
% \mathtt{subject~to} \quad& \mathcal{I}_{N_{l,u}} \leq \mathcal{I}_{N_{l,u}}^{\max}, &\\
\mathtt{subject~to} \quad& \mathcal{I}_{N_{l,u}} \leq \mathcal{I}_{N_{l,u}}^{\textsf{csk}}, &\\
& \mathcal{I}_{N_{j,u}} \leq \mathcal{I}_{N_{j,u}}^{\max}, &\\
& \mathcal{E}_\textsf{s}^{\textsf{tx}} \geq \mathcal{E}_{\textsf{r}_u}, &
\end{align}
where $f\left(x\right) = \max\left\lbrace x, 1/x\right\rbrace$.

What I am trying to achieve is to have a particular illuminance level at patch $u$. However, the code I have written only gives me the desired illuminance when the patch is at the centre.

Basically, the luminous intensity of one of the LED arrays is being modulated through the use of colour shift keying in visible light communication to transmit some data. In the end I would like to end up having a constraint makes the intensity of the transmitting LED array dominant amongst all four but for now, I need to be able to get the desired illuminance at the two patches labelled $Rx_{0}$ and $Rx_{1}$.

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

np.set_printoptions(precision=3, suppress=True)
plt.rcParams.update({'font.size': 16})
plt.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']  # for \text command
plt.rc('font', **{'family': 'serif', 'serif': ['Palatino']})
plt.rc('text', usetex=True)

tx_array_max_int = 1465  # max_intensity of transmitting array

nttx_array_max_int = 2743  # max_intensity of non-transmitting array

desired_illuminance = 450  # In lux

txrx_geo = np.array(([0.163265306122449, 0.06201550387596899, 0.163265306122449, 0.06201550387596899],
                    [0.0761904761904762, 0.17777777777777778, 0.0761904761904762, 0.17777777777777778]))

n, m = txrx_geo.shape

lamp_intensity = cp.Variable(m)

objective = cp.Minimize(cp.max(cp.maximum((txrx_geo * lamp_intensity) / desired_illuminance,
                                          cp.inv_pos((txrx_geo * lamp_intensity) / desired_illuminance))))

constraints = [lamp_intensity[0] <= tx_array_max_int,
               lamp_intensity >= 0,
               lamp_intensity <= nttx_array_max_int,
               desired_illuminance <= txrx_geo * lamp_intensity]

problem = cp.Problem(objective, constraints)

print(problem.solve())

print("Our solution is: ", problem.status)
print("--------------------------------------------------------------------------------------------")

optimal_intensities = np.ravel(lamp_intensity.value)
print("The optimal intensities (per transmitter array) are: ", optimal_intensities)
</code>```


问:
Shape recognition of tangientially connected polygons using openCV
Lately, I have encountered the following problem, namely how to recognize the shapes of particular polygons. I have realized that quite straightforward to do so when the polygons are separated by using OpenCV library as shown in the following topics
example 1
example 2

My question, however, is quite more complex and concerns how to detect polygons when they share a side,
like for instance in the picture attached to this question

Anyone has any suggestions? Maybe some other tool might be able to perform such an analysis?


问:
Black box optimization
I have a simulation which gives a scalar result depending on the choice of some continuous design variables. I am trying to minimize the output of the simulation. As a first step, I want to study the convexity of the problem. Is there a methodology that allows us to determine whether a black box function  (no analytical form) is convex or not ?

I have already found a pair $\{\vec{x}_1, \vec{x}_2\}$ that does not satisfy the midpoint convexity:

\begin{equation}
f(\vec{x}_1)+f(\vec{x}_2) \geq 2f\left(\frac{\vec{x}_1+\vec{x}_2}{2}\right)
\end{equation}

Is this a sufficient condition to say that the function is nonconvex ? Furthermore, if it is nonconvex how does that affect the choice of the optimizer ?

I have read that evolutionary algorithms could be a good choice for black box optimization. However, I do not have much experience in the field of optimization and I would appreciate it if someone could point me  in the right direction.

Title: Black box optimization


问:
Python routine to calculate shape resonances of H2
I am currently doing a project in which my aim is to write a program that can be used to calculate single and multi-channel shape resonances. So I'm looking at bound states and quasi-bound states.

Specifically, I'm looking at the hydrogen molecule in spherically symmetric potentials.

So the method I follow (renormalized Numerov method from Johnson1977; DOI: https://doi.org/10.1063/1.435384) is I write a function that adds the centrifugal barrier term to the Schroedinger equation since I'm working in radial coordinates. Once that is done, I use the renormalized Numerov method to solve the Schroedinger equation and return my normalized wave function, where I integrate outwards from small r to some r_max that lies in the continuum region. I show the code here since I've already completed this part of the code and it works. Basically if I run it for some start energy for integration (energy that is above the dissociation threshold of the molecule), I get a plot that shows the wave fn oscillating in the bound region and dissociating in the continuum region).

What I want to do next is to match the wave function solution to the linear combination of regular and irregular spherical Bessel and Neumann functions so I can extract the phase shift so I can plot the Lorentzian lineshape so I can fit it and extract the energy and FWHM of my resonances. This process is similar when calculating cross sections where the phase shift is extracted to the compute scattering amplitudes (described in Atkins & Friedman, Chapter 14, Scattering theory).

I don't know how to implement this matching process and was hoping to get some tips as how to use the scipy special package and then Bessel and Neumann functions to do this. Any tips are greatly appreciated.

My understanding is that I have to run the outwards integration for a given N (rotational quantum number) on a grid of energies and use Bessel and Neumann functions to extract the phase shift. I just don't know how to start implementing this.

Thanks in advance for all help and responses. Please let me know if you need additional information and clarification!

Here is also an extract of one of the documents my supervisor gave me for the mathematical process:

The Numerov code:

```<code>import numpy as np
import scipy.constants as const
import scipy.io as scio
from scipy import interpolate
from scipy.integrate import simps
from matplotlib import pyplot as plt

# Define and initialize variables
mu = const.value("proton-electron mass ratio")/2  # reduced mass of H2
EH = const.value("hartree-inverse meter relationship")/100  # hartree energy in cm-1
#hbar = 1

# Load potential energy surface
dat = scio.loadmat('Xp_pot_0_150_0.01.mat')
Xp_pot = np.array(dat['Xp_pot'])

# Interpolate potential on grid
RR = np.arange(0.5, 20, 0.005)
pot_intp = interpolate.PchipInterpolator(Xp_pot[1:, 0], Xp_pot[1:, 1])
yy = pot_intp(RR)
plt.plot(RR, yy, label="H2 Potential")
plt.ylim(-0.62, -0.49)
plt.xlabel('R (a.u.)')
plt.ylabel('Pot (a.u.)')
plt.legend()
plt.show()

pot = np.c_[RR, yy]

# Define effective potential and add centrifugal barrier term
def veff(pot, mu, N):

    vf = pot[:, 1] + N*(N+1)/(2*mu*pot[:, 0]**2)
    effPot = np.c_[RR, vf]

    return effPot

# Plot effective potential
effpot = veff(pot, mu, 7)
plt.plot(effpot[:, 0], effpot[:, 1], label="Effective Potential")
plt.ylim(-0.5-10/EH, -0.5+40/EH)
#plt.ylim(-0.62, -0.5+40/EH)
plt.xlabel('R (a.u.)')
plt.ylabel('Pot (a.u.)')
plt.legend()
plt.show()

# Define Numerov method
def renorm_numerov_outwards(ES, pot, mu, N):
    # Renormalized Numerov integration outwards
    # ES is start energy somewhere in continuum
    # pot is given potential energy surface

    # Determine step size h from potential
    h = pot[1, 0] - pot[0, 0]
    h2 = h**2
    h12 = h2/12

    max_ind = pot.shape[0]

    # Compute Q (energy minus eff. potential)
    Q = 2*mu*(ES-veff(pot, mu, N)[:, 1])  # Choose energy above dissociation threshold

    # Calculate Tn and Un
    Tn = -h12*Q
    Un = (2 + 10*Tn)/(1-Tn)

    ### FUNCTIONS to calculate Rn, Fn, PsiN and PsiN
    def out_int(Un):  # outwards integration of recurrence relation
        Rn = np.zeros(max_ind)
        Rn[0] = 1e300  # Set initial value biggest possible value (inf) so Psi BCs are fulfilled
        for dr in np.arange(1, max_ind):
            Rn[dr] = Un[dr] - 1/Rn[dr-1]

        return Rn

    def out_Fn(Rn):  # computation of recurrence term used to compute wave fn
        Fn = np.zeros(max_ind)
        Fn[max_ind - 1] = 1
        for dr in np.arange(max_ind-2, -1, -1):
            Fn[dr] = Fn[dr+1] / Rn[dr]

        return Fn

    def trans_Fn(Fn):  # functions that transforms Fn to wave fn
        Psi = Fn / (1-Tn)
        return Psi

    def norm_Psi(Psi):  # function that normalized wave fn using Simpson's
        Psi_square_int = simps(Psi**2, pot[:, 0])
        Psi = Psi/np.sqrt(Psi_square_int)
        return Psi

    # todo: write function to determine sign changes and count nodes

    # Integration
    Rn = out_int(Un)

    # Calculate Fn and wavefn
    Fn = out_Fn(Rn)
    Psi = trans_Fn(Fn)

    # Normalize wavefn
    Psi = norm_Psi(Psi)

    return pot[:, 0], Psi
</code>```


问:
Output solution vector in Petsc
I am using petsc to solve a linear elasticity problem discretized by finite elements.The initial mesh is read by a mesh file and the distribution in each processor is done using METIS.I am using only the basics routines from petsc in order to create the submatrices and vectors.After solving I want every processor to write solution in a file.My question is,since petsc re-order the dofs internally(for example processor 0 goes from 0-N etc.) how can I track which dof from initial mesh corresponds to which solution value obtained from petsc solution vector?For example processor 1 has {4,6,8,5,13,9} nodes from initial mesh and the global range of petsc solution vector goes from 5 to 11(using VecGetOwnershipRange).How can I know which value is assigned for node 5?</p>

Title: Output solution vector in Petsc


问:
Mesh partitioning with METIS
I am trying to use METIS-5.1.0 edition in order to partion a FE mesh. For demostration purposes I created  2x2 rectangle mesh and tried to partition it. However, I notice a weird behaviour in my code. When i am using 2 parts to partion the mesh everything works file, but if i increase the number of partions at 3 or 4 I am not getting what I expect. For example for partioning parts=4 I get {0,0,2,2} as partitioing output, while I expect something like {0,1,2,3}.Any idea why is this happening?

Here is my example code:

Title: Mesh partitioning with METIS


问:
How to make a time-parametrization slower around a point but not too slow?
For animation purposes (see below) I need to use a parametrization which 'slows down' near a specific value;

More precisely, I am looking for a 'nice' monotonically increasing  function $r:[0,1] \to [0,1]$ satisfying $
r(0)=0, r(1)=1, $ whose derivative is relatively small around the point
$r^{-1}(\frac{1}{2})$.

After some playing around with polynomials, I came up with the following example:
$$
r(t)=\frac{1}{2} \big((2 t - 1)^3 + 1\big).
$$
This choice of $r$ satisfies $r'(t)=0 \iff t=\frac{1}{2}=r^{-1}(\frac{1}{2})$, so $r'(r^{-1}(\frac{1}{2}))=0.$

However, this is "too slow" for me - it turns out that having derivative zero makes the animation looks like it "stopped" (depending of course on the uniform discrete time-step size that I am using).

Question: How can I make an adaptation which makes the animation still slow but not "too slow" around $r^{-1}(\frac{1}{2})$?

Here is my animation: It is a shrinking disk - which starts at radius $1$, and shrinks to a point. After it reaches to $r=\frac{1}{2}$ it starts 'twisting' as well as shrinking (the twisting and shrinking occurs simultaneously beyond that stage). I want to focus attention on the part around where the twisting starts.

(If it's relevant I am working with Mathematica, but I guess there should be an analytic platform-independent solution to this problem.)


问:
Neural networks trained with subspace iteration algorithms for solving SCF eigenvalue problems
I am toying with the idea of using current libraries available with SCF (Self-Consistent Field) subspace iteration codes to train an ANN (Artificial Neural Network) to solve for SCF eigenvalue problems. Are there any obvious deficiencies with this idea?

Firstly, I would like to make sure that it is possible to train an ANN to approximate an algorithm, in this case, a numerical algorithm. I would hope that the ANN would provide some speedups over the original algorithm.

Since I don't have the requisite knowledge to form the machine learning problem in a suitable and original mathematical framework, I will rely on other people's numerical codes to train an ANN to do what I'm looking to accomplish. This would be the starting point for further exploration.


问:
What is this QR-factorization-based preconditioning called?
I have recently started to delve into someone else's code, and there is a part in there I don't quite understand. The authors of the code use some form of pre-conditioning to speed up the optimization. I wanted to inquire if you are familiar with this pre-conditioning, and could recommend some literature to read up on the method. Unfortunately, no name was provided, so I'll describe the author's approach in the following.

Assume we have a collection of $N$ points in some $D$-dimensional space $\mathbb{R}^D$. This space is defined in some polynomial basis, and we have an $N\times D$ matrix $\boldsymbol{B}$ containing the coordinates of the $N$ points in that basis. We now want to optimize some objective function $f(\theta,\boldsymbol{B})$ for a fixed $\boldsymbol{B}$ and variable $\theta$ (a vector of length $D$). We do this with a quasi-Newton optimizer, specifically BFGS.

For the pre-conditioning, the authors now first

square the matrix $\boldsymbol{B}$ element-wise, then sum it along the rows ($N$-dimensions), and take the square root of the resulting vector. This yields a $D$-vector $\psi_n$.

The authors then create a diagonal square matrix $\boldsymbol{D}^{-1}$ with $\frac{1}{\psi_n}$ as the diagonal entries.

With the matrix $\bf{D}^{-1}$, the authors obtain a matrix $\boldsymbol{\Psi}=\boldsymbol{B} \cdot \boldsymbol{D}^{-1}$. I assume this normalizes the dimensions of $\boldsymbol{B}$.

Then, the authors take the QR factorization of $\bf{\Psi}$ and only keep the $N\times D$ upper triangular matrix $\boldsymbol{R}$

With this, we authors calculate a "preconditioning matrix" $\boldsymbol{A}=\boldsymbol{D}^{-1}\cdot \boldsymbol{R}^{-1}$. This matrix must be $D\times N$.

Apparently, this allows us to do a "regularized optimization" $f(\theta)$=$f(A\cdot\theta^\top)=g(\theta^\top)$. Their comments suggest that this somehow turns the optimization problem more 'spherical', reducing steep gradients in some dimensions. (Apparently, this makes it less difficult for the algorithm to walk up ridges, and then close in on the optimum.) After optimization, the results from the regularized objective $\theta^*_r$ can then be converted back by re-scaling it with $\theta^*=\boldsymbol{A}\cdot \theta^*_r$.

Do you know what method the authors are using here?


问:
Numerical simulation of magnetic dipole in a inhomogeneous magnetic field
Goal:

I want to use Python to illustrate how a magnetic dipole with magnetic moment m2 moves in a non-homogeneous magnetic field in a 2D-Plane. This field is generated by another magnetic dipole with moment m1. However, I get results that are not plausible.

For the calculation of the trajectory of m2 I use the equations from wiki - force between magnets

Theoretical Background:

Given is a magnetic dipole with constant moment m1 and constant position r1!. m1 generates a magnetic flux density

$${\displaystyle \mathbf {B}_1 (\mathbf {r} )={\frac {\mu _{0}}{4\pi }}{\frac {3\mathbf {\hat {r}} (\mathbf {\hat {r}} \cdot \mathbf {m} )-\mathbf {m} }{|\mathbf {r} |^{3}}}\qquad (1)}$$

where r is a vector from the center of the magnetic dipole to the location where the magnetic field is measured. In this non-homogeneous magnetic field a magnetic dipole with moment m2 and position r2 is placed. Its magnetic field can be neglected B2 = 0.

Only the dynamics of m2 shall be considered. According to the linked Wikipedia page, m2 will be subjected to a torque D and a force F. The torque is neglected D=0. The magnetic moment m2 is therefore constant.

Due to the acting force, the position vector r_2(t) depends on time. Its time evolution is obtained from the solution of the equation of motions, which is my objective:

$$ \frac{d^2}{dt^2}\mathbf{r}_2 (t) = \binom{\dot{r}_{2x}}{\dot{r}_{2y}} = \mathbf{F}$$

According to Wikipedia, m1 exerts the following force on m2: (r = r2 - r1, r = ||r||_2):

$$\begin{align} \mathbf{F} &= \nabla |\mathbf{m}_2\cdot \mathbf{B}_1| \qquad(3)
\\\\ &= {\frac {1}{ r^{5}}}\left[(\mathbf {m} _{1}\cdot \mathbf {r} )\mathbf {m} _{2}+(\mathbf {m} _{2}\cdot \mathbf {r} )\mathbf {m} _{1}+(\mathbf {m} _{1}\cdot \mathbf {m} _{2})\mathbf {r} -{\frac {5(\mathbf {m} _{1}\cdot \mathbf {r} )(\mathbf {m} _{2}\cdot \mathbf {r} )}{r^{2}}}\mathbf {r} \right] \qquad (4)
\\ &= \frac{\mathbf{m}_1 \mathbf{r}}{r^5} \mathbf{m}_2 + \frac{\mathbf{m}_2\mathbf{r}}{r^5}+\frac{\mathbf{m}_1\mathbf{m}_2}{r^5}\mathbf{r} - \frac{5(\mathbf{m}_1 \mathbf{r})(\mathbf{m}_1\mathbf{r})}{r^7}\mathbf{r}\\\\
&= \mathbf{c}_1 + \mathbf{c}_2 + c_3 (\mathbf{r}_2 - \mathbf{r}_1) - c_4 (\mathbf{r}_2 - \mathbf{r}_1) \qquad (5)\\\\
&= \mathbf{c}_1 + \mathbf{c}_2 + \mathbf{r}_1 (c_4 - c_3) + \mathbf{r}_2 (c_3 - c_4)\qquad (6)\\\\
&= \mathbf{c}_1 + \mathbf{c}_2 + \mathbf{c} + d\cdot\mathbf{r}_2, \quad \mathbf{c}_1, \mathbf{c}_2, \mathbf{c}\in \mathbb{R}^2, d\in\mathbb{R} \qquad (7)
\end{align} $$

The introduced variables c1,c2, c, d should be self-explanatory. Now we can define 4 ODEs and solve them numerically with the Runge-Kutta method:

$$\begin{align}
\dot{r}_{2x} &= v_{2x} \qquad (8)\\\\
\dot{r}_{2y} &= v_{2y} \qquad (9)\\\\
\dot{v}_{2x} &= c_{1x} + c_{2x} + c_x + d\cdot r_{2x} \qquad (10)\\\\
\dot{v}_{2y} &= c_{1y} + c_{2y} + c_y + d\cdot r_{2y} \qquad (11)\\\\
\end{align}$$

What to expect?

An attractive force should act on the magnet m2, so that

$$ t \gg 1: \mathbf{r}_2 = \mathbf{r}_1 $$

My Results

m2 is repelled instead of being attracted.

Probably I made a mistake in the implementation. But I have already checked it several times and found no error... Does anyone of you know what went wrong

Python 3 Code

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

"""
Important quantites that will influence the dynamics

m1:  Magnetic Moment of the dipole that creates a constant magnetic field
m2:  Magnetic Moment of a dust particle
r1:  Position of m1: r1(t) = const
r2:  Position of m2: r2(t) will depend on time
"""

# Constants
m1 = np.array([0, 1], float)
m2 = np.array([0, -1], float)
r1 = np.array([0, 0], float)

# Initial conditions of the of m2 R = [r2_x, r2_y, v2_x, v2_y]
R = np.array([0, 3, 0, 0], float)

def distance(r1, r2, n):
    """
    computes the distance ||.||_2 between position r1 and r2
    """

    return np.linalg.norm(r1 - r2) ** (n)

def pre_factors(r1, r2, m1, m2):
    """
    Computes the prefactors in the equation of motion
    d^2 r2 / dt^2 = F = c1 + c2 + c + r2 * d   r2,c1,c2,c € R^2, d € R
    """

    r = r2 - r1
    l = distance(r1, r2, 5)

    c1 = ( np.dot(m1, r) * m2 ) / l
    c2 = ( np.dot(m2, r) * m1 ) / l

    c3 = np.dot(m1, m2) / l
    c4 = ( 5 * np.dot(m1, r) * np.dot(m2, r) ) / distance(r1, r2, 7)

    c = (c4 - c3) * r2
    d = c3 - c4

    return c1, c2, c, d

def f(r1, R, m1, m2):
    """
    now we can compute the righthand side of the four 1st order ODEs:
    d/dt r2_x = v2_x
    d/dt r2_y = v2_y
    d/dt v2_x = c1_x + c2_x + c_x + d*r2_x
    d/dt v2_y = c1_y + c2_y + c_y + d*r2_y
    """

    r2x = R[0]
    r2y = R[1]
    v2x = R[2]
    v2y = R[3]

    c1, c2, c, d = pre_factors(r1, R[:2], m1, m2)

    f_r2x = v2x
    f_r2y = v2y
    f_v2x = c1[0] + c2[0] + c[0] + d * r2x
    f_v2y = c1[1] + c2[1] + c[1] + d * r2y

    return np.array([f_r2x, f_r2y, f_v2x, f_v2y], float)

# Constants for the Simulation
a = 0.0  # t_0
b = 20.0  # t_end
N = 1000 # number of iterations
h = (b-a) / N # time step, e.g. h = (100 - 0) / 1000 = 0.1s

# Create lists to store the computed values
tpoints = np.arange(a,b+h,h)

xpoints, ypoints = [], []
vxpoints, vypoints = [], []

for dt in tpoints:
    #print(R)
    xpoints.append(R[0])
    ypoints.append(R[1])
    vxpoints.append(R[2])
    vypoints.append(R[3])

    # Runge-Kutta-4th Order Method
    k1 = dt * f(r1, R, m1, m2)
    k2 = dt * f(r1, R + 0.5 * k1, m1, m2)
    k3 = dt * f(r1, R + 0.5 * k2, m1, m2)
    k4 = dt * f(r1, R + k3, m1, m2)
    R += (k1 + 2*k2 * 2*k3 + k4)

plt.style.use('seaborn')
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 3))
axes[0].plot(tpoints, xpoints, label = '<span class="math-container">$r_{x_2}(t)$</span>')
axes[0].plot(tpoints, ypoints, label = '<span class="math-container">$r_{y_2}(t)$</span>')

axes[0].legend()

axes[1].plot(tpoints, vxpoints,'r', label = '<span class="math-container">$v_{x_2}(t)$</span>', alpha = 0.5)
axes[1].plot(tpoints, vypoints, 'k',label = '<span class="math-container">$v_{y_2}(t)$</span>', alpha = 0.5)
axes[1].legend()

fig.suptitle("Results for: <span class="math-container">$r_1 = (0,0), r_2 = (0,3),  m_1 = (0,1), m_2 = (0,-1)$</span>")
#plt.savefig("2.pdf" , format='pdf',  bbox_inches="tight")
fig.tight_layout()
</code>```


问:
General question related to BFGS
Based on my basic understanding of the BFGS method, the algorithm will iterate until the gradient norm is less than or equal to a set value called "gtol" in the case of Python.

However, when using this method, and checking the output the following is showing:

Iterations: 2
Function evaluations: 76
Gradient evaluations: 13

"Desired error not necessarily achieved due to precision loss."

This got me confused. Shouldn't the algorithm iterate multiple times to reach convergence? How did it reach convergence by only iterating twice? In this case, why can we specify a maximum number of iterations if the algorithm converges this quickly?

Title: General question related to BFGS


问:
Improve code of logarithmic quantizer
I am implementing a logarithmic quantizer which is defined as follows:

$$ q(u) = \begin{cases} u_i , \frac{u_i}{1+\delta} < u < \frac{u_i}{1-\delta} \\
0 , 0 \leq u \leq \frac{u_o}{1+\delta} \\
-q(-u), \ u < 0 \end{cases} $$

where the different parameters are computed based upon these expressions:

$$ u_i = p^{1-i}u_o, \ i = 1,2,\dots $$
$$ \delta = \frac{1-p}{1+p} $$
$$ u_{min} = \frac{u_o}{1+\delta} $$

More information regarding the quantizer can be found at this paper: Adaptive Backstepping Control of Uncertain Nonlinear Systems with Input Quantization.

In order to code a function to implement it in MATLAB, I wrote the following piece of code:

```<code>function q_u = logarithmic_hysteretic_quantizer(u,step,u_min)

    u_o = u_min*(1+step);
    p = (1-step)/(1+step);

    if u < 0

        q_u = -logarithmic_hysteretic_quantizer(-u,step,u_min);

    elseif ( (u >= 0) && (u <= u_o/(1+step)) )

        q_u = 0;

    else

        i = 1;

        while (1)

            u_i = p^(1-i) * u_o;

            if ( (u > u_i/(1+step)) && (u <= u_i/(1-step)) )

                q_u = u_i;
                break;

            end

            i = i + 1;
        end
    end
end
</code>```

Now, my issue is to improve the code as much as I can. For example, the while(1) loop, which codes the different quantization levels, is something that should go away and be replaced. Any thoughts would be really appreciated.

Title: Improve code of logarithmic quantizer


问:
Are FDS scheme and ROE-FDS scheme the same?
In SU2 documentation about the numerical schemes (https://su2code.github.io/docs_v7/Convective-Schemes/) FDS is mentioned as a "standalone" method.

Yet, when looking online for more information about FDS scheme, it seems it is just a preconditioning for the ROE scheme. And this seems coherent with the fact that has reported here it is a necessary treatment for the poor convergence of compressible codes (the main cause is the artificial dissipation) for this type of simulations.

Hence my question: are FDS scheme and ROE-FDS scheme the same with FDS being only a treatment to make ROE scheme work well in incompressible cases?


问:
Sparse linear solver in fortran working with REAL16
I need some (direct) sparse linear solver for fortran, which works with REAL16 data type. Any suggestions? Both Pardiso and MUMPS support only REAL8. (identical question: https://math.stackexchange.com/questions/4016856/sparse-linear-solver-in-fortran-working-with-real16?noredirect=1#comment8293130_4016856)</p>


问:
Best way to compute given functional with accuracy:
I need to plot the following functional with accuracy:

$$
I(x,s) =\int_0^\infty\mathrm dy \frac{F(x + \mathrm iy,s) − F(x −\mathrm iy,s)}{\mathrm e^{2πy}-1},
$$
Where $ F(z,s) = \dfrac{1}{z^s\Gamma(\sin^2[π\Gamma(z)/(2z)])} $.

And let us restrict $s\in[0,1]$

What is the most efficient way of computing this integral?

Nature of functional as $x\rightarrow\infty$ from the computation ?

I computed some relatively small values using Mathematica which suggests the function is oscillatory with damping. But I need big values greater than x=100 and at s=1


问:
Finding the extrema of a transition probability function for a quantum walker on a graph

The goal

Implement some Python code to find the extrema points of a function that is strongly oscillating.

The background

Let $G$ be a connected graph with $n$ points with Laplacian matrix $L(G)$. We may represent each vertex of $G$ as a canonical basis vector $|k\rangle$ for $k\in \{1, ..., n\}$, spanning a position Hilbert space $\mathcal{H}$. A notable element of such a space is the flat superposition
$$|s\rangle=\frac{1}{\sqrt n}(|1\rangle+|2\rangle+...|n\rangle).$$

Within such a model, assume you are given the Hamiltonian
$$ H(\gamma, w, G)=\gamma L(G)-|w\rangle\langle w|$$
where $|w\rangle$ is one of the position basis vectors and $\gamma$ a real parameter, and consider the probability function $$P_{G,n}(t,\gamma)=|\langle w| U_H(t,\gamma)| s\rangle|^2$$ for a system prepared in the flat state to evolve under $H$ into $|w\rangle$.

For a general $G$, $P$ is a strongly oscillating function. Calculating the right maximum for a pair $(\gamma, t)$ given $G, n$ and $|w\rangle$ as fixed parameters is challenging, because most algorithms will rather find something local and then halt.

My implementation

The flat state $|s\rangle$ and the target state $|w\rangle$ can be implemented in Python as

```<code>def flat(n):
    s = []
    for j in range(n):
        s.append(1/m.sqrt(n))
    return np.array(s)

def target(n,p):
    w=[]
    for j in range(n):
        if j==p:
            w.append(1)
        else: w.append(0)
    return np.array(w)
</code>```

The probability function is then implemented as

```<code>def probability(x,G,n,p):
    t,gamma=x
    w=target(n,p)
    U=expm(-1j*t*(gamma*nx.laplacian_matrix(G).todense()-np.outer(w,w)))
    psi=np.dot(U,flat(n))

    return -abs(np.dot(w,psi))*abs(np.dot(w,psi))
</code>```

where I have used an intermediate state $|\psi(t)\rangle=U_H(t)|s\rangle$.

The Scipy library has a function minimize in the scipy.optimize module (note: the reason why the probability function returns a negative quantity is that there is no direct maximizing function, so we find the maxima of $f$ by finding the minima of $-f$). The documentation can be found here. I have used it to implement the minimizer function:

```<code>def minimizer(G,n,p,x0):
    res=minimize(probability,x0,args=(G,n,p))
    return [-res.fun, res.x[0]]
</code>```

returning an array with $-P_{\text{max}}$ and the corresponding $(t_{\text{max}}, \gamma_{\text{max}})$. The argument x0 is an array containing the initial guess on $t$ and $\gamma$ respectively.

The issue

scipy.optimize.minimize works only after setting an initial guess x0, and the result of the optimization strongly depends on such a choice. As a consequence, obtaining a graph of the parameters that maximize probability as a function of the number of vertices $n$ becomes impossible, as any x0 is only good in a range $(n_1, n_2)$ with endpoints depending on the chosen graph.

So here I ask for help: is there any way to make sure minimizer finds the global minima of probability (and in particular the first of them, if there are many) consistently in a dimension range $n \in (0, N\ge 100)$, regardless of the underlying graph?

I hope everything is as clear as possible, and thanks in advance to whoever might answer.

An example: the complete graph

For $n=20$, the probability function for the complete graph looks like this:

The complete graph is a good benchmark for testing the code as we know that $P_{\text{max}}=1$ for $t_{\text{max}}=\pi\sqrt{n}/2$ and $\gamma=1/n$ regardless of the chosen vertex $|w\rangle$.

I have written a snippet of code that plots either of these quantities as functions of $n$, taking $|w\rangle=|0\rangle$:

```<code>#i=0 plots probability_max,  #i=1 plots t_max, #i=2 plots gamma_max
#plots in the interval [n=2, n=maxi)
def plotComplete(x0,maxi,p,i):
    zPoints = []
    yPoints = []
    refX=[]
    refY=[]

    for z in range(2,maxi):
        y=round(minimizer(nx.complete_graph(z),z,p,x0)[i],5)
        zPoints.append(z)
        yPoints.append(y)

    pyplot.plot(zPoints, yPoints, 'b')

    if i==1:
        for z in range(2,maxi):
            y=m.pi*m.sqrt(z)/2
            refX.append(z)
            refY.append(y)

        pyplot.plot(refX, refY, 'r')

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

Calling plotComplete([0.1,1],100,0,0) (notice that the initial guess on $t$ is $0.1$; changing it to, for example, $0$ or $1$ gives wildly different results) gives as expected the flat line

for the first $100$ complete graphs. Calling plotComplete([0.1,1],100,0,1) plots $t_{\text{max}}$ with the same initial guess:

Here, I have made two overlapping plots: the red one is the theoretical value $\pi\sqrt{n}/2$, the blue one is the implemented minimizer. For the most part up to around $n=60$ the two plots are in almost perfect agreement, then after that the initial guess $0.1$ is no more good enough to produce the correct result. The blue spikes appear any time the wrong minimum is found.


问:
Cauchy problem ill-posed?
Find the solution to the Cauchy problem consisting of the wave equation :
$$u_{xx}-u_{yy}=0$$ together with initial conditions: $$ u(x,0)=0,$$ $$u_{y}(x,0)=g(x)$$ for some known initial datum $g$.

Is the problem well-posed or ill-posed?

Why?

My thoughts:

By changing variables we have :
$$ u(x,y) = \xi(x-y)+\eta(x+y) $$
First , we note that the initial condition $u(x,0)=0$ gives us:

$$\begin{align*}
0 =& u(x,0) \\
=&\xi(x-0)+\eta(x+0) \\
=&\xi(x)+\eta(x)
\end{align*}$$

We then derived the first equation $$\xi(x)+\eta(x)=0.$$
Next,we determine the derivative $u(x,y)$ with respect to $y$:

$$\begin{align*}
u_{y}(x,y) =& \frac{\partial}{\partial y}(\xi(x-y)+\eta(x+y)) \\
=& -\xi'(x-y)+\eta'(x+y)
\end{align*}$$

Next the initial condition $$u_{y}(x,0)=g(x)$$ which implies:
$$\begin{align*}
g(x) =& u_{t}(x,0) \\
=& -\xi'(x-0)+\eta'(x+0) \\
=&-\xi'(x)+\eta'(x)
\end{align*}$$

We have then derived the second equation $-\xi'(x)-\eta'(x)=g(x)$
as well.

Solving the equation $$\xi(x)+\eta(x)=0$$ to $\eta(x)$ in the first equation we obtain: $$\xi(x)=-\eta(x).$$

Differentiate both sides of the previous equation with respect to $x$:
$$\xi'(x)=-\eta'(x)$$

Next , we replace $\eta'(x)$ by $-\xi'(x)$ in the equation $$-\xi'(x)-\eta'(x)=g(x)$$ we take:

$$\begin{align*}
g(x) =& -\xi'(x)+ \xi'(x)   \\
     =& -\xi'(x)- \xi'(x)
\end{align*}$$

Combining terms :

$$-2\xi'(x) = g(x)$$

and dividing each side by $-2$

$$ \xi'(x) = -\frac{1}{2}g(x)$$

Integrating each side wrt $x$:

$$ \int_{x_{0}}^{x} \xi'(x)dx = \int_{x_{0}}^{x} -\frac{1}{2}g(x) dx $$
$$ \xi(x) -\xi(x_{0})  = -\frac{1}{2} \int_{x_{0}}^{x} g(\xi) d\xi  = $$

Add $\xi(x_{0})$ to each side:

$$\xi(x) =-\frac{1}{2} \int_{x_{0}}^{x} g(\xi) d\xi + \xi(x_{0})$$

We know that the general solutuon is $u(x,y) = \xi(x-y)+ \eta(x+y)$ with $\eta(x) = - \xi(x)$ and $\xi(x)-\frac{1}{2} \int_{x_{0}}^{x} g(\xi) d\xi + \xi(x_{0})$

\begin{align*}
u(x,y) =& \xi(x-y)+ \eta(x+y) \\
     =& \xi(x-y) -  \eta(x+y) \\
     =& -\frac{1}{2} \int_{x_{0}}^{x-y} g(\xi) d\xi + \xi(x_{0})+\frac{1}{2} \int_{x_{0}}^{x+y} g(\xi) d\xi - \xi(x_{0}) \\
     =& -\frac{1}{2} \int_{x_{0}}^{x-y} g(\xi) d\xi +\frac{1}{2} \int_{x_{0}}^{x} g(\xi) d\xi  \\
     =& \frac{1}{2} \left( \int_{x_{0}}^{x+y} g(\xi) d\xi -\frac{1}{2} \int_{x_{0}}^{x-y} g(\xi) d\xi  \right)\\
     =& \frac{1}{2} \int_{x_{0}-y}^{x+y} g(\xi) d\xi
\end{align*}

For the difference $$u(x,y)- \tilde{u}(x,y)$$ we will have :

\begin{align*}
u(x,y)- \tilde{u}(x,u) =& \frac{1}{2} \int_{x_{0}-y}^{x+y} g(\xi) d\xi -  \frac{1}{2} \int_{x_{0}-y}^{x+y} \tilde{g}(\xi)) d\xi \\
 =&  \frac{1}{2} G(\xi) + c -  \left( \frac{1}{2} \tilde{G}(\xi)+ k \right)
\end{align*}

Knowning that $c+k=0$
then we have that the problem is ill-posed.

Am I correct?

Title: Cauchy problem ill-posed?


问:
Fourier transform in finite element
I have a finite element solver where I am using tetrahedral elements. I am solving for electric potentials and then calculate the current densities in each element, which are constant in each element. Now, assume that I want to calculate Fourier transforms for current density in the $x$ direction ($J_x$).
$$F(J_x(x,y,z))(k_x,k_y,k_z)=\sum_e\int_e{J_{x}}_eexp(-i2\pi(k_xx(\xi,\eta,\tau)+k_yy(\xi,\eta,\tau)+k_zz(\xi,\eta,\tau)))d\xi d\eta d\tau$$
Where I utilized the Gaussian quadrature rule to calculate the integral. I am taking Fourier space points to be on a rectangular grid. Therefore, I was expecting the inverse of the transform (I am using inverse FFT) to be a set of slice interpolations of the current density. However, the inverse transform is meaningless. I wanted to ask is there anyone who tried such a method (FT from a FE model) and can you tell me what is wrong overall because clearly, something is. Thanks in advance and have a nice day.

Title: Fourier transform in finite element


问:
How to solve this boundary value problem which has more unknown than equation on MATLAB
I need your helps about solving the problem below with MATLAB. I am trying to solve 2D Stress Wave Propagation problem by using FDTD(Finite difference time domain) method on the cylindrical coord. I have discretize the domain and now ı am trying to apply boundary conditions. I have the boundary condition equation below for i=1 and j=2,...,ns-1(AB boundary).

This is my discretization;

The equation (AB Boundary Condition);

The right-hand side of the equation is known. But the left side is unknown.

Linear equation system for the equation; for i=1 and j=2,...,4(ı choose the j = 4 for the last term to make simpler explanation)

u_(i,j)^(k+1) is the radial displacement and v_(i,j)^(k+1) is theangular displacement for the (k+1). time step.

A1,B1,C1,A2,B2,C2,A3,B3,C3 are coefficients.

D1,D2,D3 are constants.

I want to find all u_(i,j)^(k+1) and  v_(i,j)^(k+1). But there are just 3 equations and 8 unknowns. I dont know the way how to solve this with matrices like Ax=B. It is not make sense.

Unknowns :  u_1,2^(k+1), u_1,3^(k+1),u_1,4^(k+1) , v_1,1^(k+1), v_1,2^(k+1), v_1,3^(k+1), v_1,4^(k+1) and v_1,5^(k+1)

Matlab code for the equations:

```<code>AB1 = (lambda(i,j)+2*nu(i,j))*(-3*u(i,j,k+1)/(2*d_r))+(lambda(i,j)/R(i,j))*(((v(i,j+1,k+1)-v(i,j-1,k+1))/(2*d_s))+u(i,j,k+1));
AB2 = -Pi-(lambda(i,j)+2*nu(i,j))*((4*u(i+1,j,k+1)-u(i+2,j,k+1))/(2*d_r));
AB1 = AB2; % AB1 IS UNKNOWN,  AB2 IS KNOWN
</code>```


问:
Why is the problem infeasible?
Given $\mathbf V_t=\mathbf v_t\mathbf v_t^H$ where $\mathbf v_t=\left(e^{j\theta_{1}},e^{j\theta_{2}}\right)^H$:
\begin{equation*}
\begin{array}{ll}
\underset{\mathbf V}{\operatorname{minimize}} & 1 \\
\text { subject to } & \operatorname{diag}(\mathbf V)=1\\
&\|\mathbf V\|_*-\operatorname{trace}\left[\boldsymbol\lambda\boldsymbol\lambda^H\left(\mathbf V-\mathbf V_t\right)\right]-\|\mathbf V_t\|_2\leq 0
\end{array}
\end{equation*}
where $\|\bullet\|_*$ and $\|\bullet\|_2$ is the nuclear norm and 2-norm. $\boldsymbol\lambda$ is the leading eigenvector of $\mathbf V_t$.

Since $\mathbf V=\mathbf V_t$ is a feasible solution, the optimization problem is feasible. Why does CVX show the problem is infeasible?

clear;clc;close all;

Vt = rand(2,1);
Vt = (exp(1i*Vt)*exp(1i*Vt)');
[lambda,~] = eigs(Vt,1,'largestabs');
cvx_begin sdp
    variable V(2,2) hermitian semidefinite
    minimize(1)
    subject to
        diag(V) == 1
        norm_nuc(V) - real(trace(lambda*lambda'*(V-Vt))) - norm(Vt,2) <= 0
cvx_end
```

Title: Why is the problem infeasible?


问:
Interpolate location based on triangular sensor array
There are 3 sensors (A, B, C) on a plane, located in the corners of a (known) equilateral triangle.

I want to calculate the (2D) location of an object (X) inside that triangle.

One sensor returns one scalar value per reading, but in order to compensate for environmental effects, I always use the values relative to each other in order to get stable values.

This means I get one value between 1.0 and 0 for each pair of sensors depending on whether the object is on the median line between two sensors (ratio=1.0), directly in one corner of the triangle (ratio= ~0.3) or even further away.

Here you can see some plots that show the idea.

The red dots are the sensor locations. I have plotted selected ratios of the measurements depending on the position of the object X.

So if X moves on a plotted curve, the ratio of sensor readings stays constant.

My goal in the end is to find X's location depending on the three sensor ratios.

Ratio B/A :

All Ratios (Ratios A/B and B/A highlighted) :

While I can make calibration measurements, these take a lot of time and are only practical for roughly 1/10th of the targeted spatial resolution. (i.e.: I can make a calibration measurement each meter, but I am targeting a 0.1m resolution.) I am, however not limited on how I select the location for the calibration measurements.

I have three values per measurement and have to somehow do a reverse interpolation.

Unfortunately, the physical model behind it is very complex (and proprietary) so I don't have a way to use it as a starting point. Also there are a lot of nonlinearities and other unknowns in the system.

What I've tried:

My first attempt was to fit circles to my calibration measurements, and reverse-interpolate on the positions of center and radius. So I would calculate a circle per sensor ratio and look for the point where the three circles would (roughly) intersect. Unfortunately this did not really work out - especially as the circle's radius approaches infitity if my sensor-ratio approaches 1.

So in the end I've ditched the geometrical approach and have looked more into numerical solutions, but I haven't found a solution that could work in my case.

I hope my explanation was clear and precise. I'm grateful for any pointers in the right direction!

(I'm currently working with Python/numpy)

More information as answer to the question in the comments

I have spent some time thinking about the reverse interpolation approach.
I figured that I would need a way to reverse interpolate each 3D surface and obtain a 2D function (per sensor pair).
I would then have to find the intersections of the three 2D functions.
Ideally, I would then find the single point where the error is minimal.

I am not sure how to do this.

OwnerDisplayName: Prima


问:
Nonsymmetric permutations for LU factorisation of symmetric matrix
Let $A$ be a symmetric matrix. It is then well known that computing the LU factorisation of $PAP^T$ instead of $A$ for a suitably chosen permutation matrix $P$ can greatly reduce fill-in. My question is: Can we show that considering nonsymmetric permutations, i.e. $PAQ^T$, cannot reduce fill-in even further?


问:
Poisson equation, stiffness matrix positive definiteness, Dirichlet boundary conditions
I have a question regarding the positive definiteness of the stiffness matrix. Specifically, I believe that it should be positive definite only when at least one Dirichlet point is given, so I would like a clarification on which is the precise moment that the stiffness matrix becomes positive definite. Let us consider a very simple model problem:

$$
-\Delta u(\pmb{x}) = f(\pmb{x}), x \in \Omega \\
u(\pmb{x}) = g(\pmb{x}), \pmb{x} \in \Gamma_D \\
\partial_n u(\pmb{x}) = 0, \pmb{x} \in \Gamma_N.
$$

Using the divergence theorem this can be rewritten as:

$$
\int_{\Omega} \nabla u \cdot \nabla v = \int_{\Omega}fv \\
u(\pmb{x}) = g(\pmb{x}), \, \pmb{x} \in \Gamma_D.
$$

A space $V_h \subset \mathcal{H}^1(\Omega)$ is given, with basis functions $\phi_1, \ldots, \phi_d, \phi_{d+1},\ldots,\phi_{d+m}$. Let $Dir = \{1,\ldots,d\}$ and $Ind = \{d+1,\ldots,d+m\}$. The solution would then be of the form:

$$u_h(\pmb{x}) = \sum_{k=1}^{d+m}u_k\phi_{k}(\pmb{x}), \pmb{x} \in \Omega,$$

and the linear system of equations that must be solved is:

$$
\sum_{k=1}^{d+m}u_k\int_{\Omega}\nabla\phi_i \cdot \nabla\phi_j = \int_{\Omega}f\phi_i, \, i \in Ind \\
u_i = g_i, \, i \in Dir.
$$

The full stiffness matrix and right-hand side are:

$$
W_{ij} = \int_{\Omega}\nabla\phi_i \cdot \nabla\phi_j, \, i,j \in Dir\cup Ind \\
f_i = \int_{\Omega}f\phi, \, i \in Dir \cup Ind.
$$

Then the system can be written more concisely as:

$$\pmb{W}|_{Ind \times Dir\cup Ind}\pmb{u} = \pmb{f}|_{Ind},$$

or if the Dirichlet columns are removed:

$$
\pmb{W}|_{Ind \times Ind}\pmb{u}|_{Ind} = \pmb{f}|_{Ind}-\pmb{W}|_{Ind \times Dir}\pmb{g}.
$$

I know, that if even a single Dirichlet node is provided, then the system has a unique solution, which to me implies that $\pmb{W}|_{Ind \times Ind}$ is positive definite. However, if no Dirichlet node is provided, one is left with pure Neumann boundary conditions, and an additional constraint must be provided, for example (https://fenicsproject.org/olddocs/dolfin/2016.1.0/python/demo/documented/neumann-poisson/python/documentation.html):

$$\int_{\Omega}u = 0 \implies \sum_{k=1}^{m}u_k\int_{\Omega}\phi_k = 0$$.

To me, requiring such an additional constraint implies that the matrix $\pmb{W}|_{Ind \times Ind}$ is not positive definite and that it has a zero eigenvalue. Yet, the answer here: In FEM, why is the stiffness matrix positive definite? claims that it is positive definite without referring to the number of Dirichlet nodes. Clearly I am missing something, and I would like to understand what.

Edit:

At first glance $\pmb{W}|_{Ind \times Ind}$  seems oblivious to Dirichlet nodes, since it includes only terms involving basis functions that correspond to non-Dirichlet nodes. However, introducing Dirichlet nodes modifies the non-Dirichlet basis functions (even non-Dirichlet basis functions have to vanish at a Dirichlet node, thus implicitly affecting the terms in $\pmb{W}|_{Ind \times Ind}$). Is anyone aware of a reference that makes this argument rigorous? More precisely a proof that the matrix $\pmb{W}|_{Ind \times Ind}$ becomes positive-definite upon the introduction of such a Dirichlet point?


问:
Linearising Nonlinear Coupled Partial Differential Equations - Alfvénic Diffusion
I am trying to solve the following coupled partial differential equations with a finite difference scheme:
$$\partial_tf+v\partial_zf+\partial_z\frac{1}{W}\partial_zf=0$$
$$\partial_tW+v\partial_zW-\partial_zf=0,$$

where both f and W depend on t and z.
In order to solve both equations simultaneously I am trying to use vector notation.
However the diffusion term in the first equation is non-linear.
How can I linearise this?

Background:

I have previously solved them in a semi-implicit scheme one after another, using the result of one equation in the other one. This however is explicit and requires very small timesteps, making the program to slow to be useful.

The basic idea is now to write them in a vector/matrix form and solve them simultaneously - again using finite differences.
This requires the equation to be of the form:
$$\begin{array}{c}
f \\
W
\end{array}=\vec{f}$$

$$\partial_t\vec{f}=\textbf{C}_0\vec{f}+\textbf{C}_z\partial_z\vec{f}+\textbf{C}_{zz}\partial_{z}^2\vec{f},$$

where all matrices $\textbf{C}$ are independ of $\vec{f}$.

The full problem consists of an additional dimension treated in LOD and a discontinuity in $v$ and $W$. Therefore, I do not expect standard ODE solvers to work.

Is there any way (substitution, or other) to make this possible?


问:
How are finite volume method boundary conditions implemented without using ghost-cells?
I'm currently trying to implement my own FVM code in cpp, but when I try to calculate the laplacian of a test function, given by
\begin{align}\phi_0=\sin(2\pi x)\sin(2\pi y),\end{align}
I get unphysical boundary values as seen in the figure. I used a Dirichlet boundary condition with 0 value on every boundary.

My code follows the methods given in "The Finite Volume Method in Computational Fluid Dynamics" by Moukalled, F. et al. according to which the laplacian (diffusion in this case) of a field is discretized for a cell using
\begin{align}\int_V\nabla\cdot(\Gamma\nabla\phi)dV\approx\sum_{f_i=1}^{n_{f,i}}\Gamma_{f_i}\frac{\mid\textbf{S}_f\mid}{\mid\textbf{d}_n\mid}(\phi_N-\phi_C)+\sum_{f_b=1}^{n_{f,b}}\Gamma_{f_b}\frac{\mid\textbf{S}_f\mid}{\mid\textbf{d}_n\mid}(\phi_b-\phi_C),\end{align}
where $\phi_C$ is the field value at the cell center, $\phi_N$ is the field value at the neighbour face(s), $\phi_b$, is the field value at the boundary, $\mid\textbf{S}_f\mid$ is the magnitude of the surface normal of the face, $\mid\textbf{d}_n\mid $ is the distance between the cell centers of the owner cell and the neighbour cells.

The above grid represents a test case with equal distances between cell centers ($\delta x$ and $\delta y$) and equal grid spacing ($\Delta x$ and $\Delta y$). So, as an example, the discretization of element 0 in the above grid is given by the algebraic relation \begin{align}a_C\phi_C+\sum a_N\phi_N=b_C , \end{align}
where \begin{align}\sum a_N\phi_N=a_{Nn}\phi_n+a_{Ne}\phi_e+a_{Ns}\phi_s+a_{Nw}\phi_w=\Gamma\Bigg(\frac{\Delta x}{\delta y}\phi_n+\frac{\Delta y}{\delta x}\phi_e+0\phi_s+0\phi_w\Bigg)=\Gamma\phi_3+\Gamma\phi_1 ,\end{align}
and

\begin{align}a_C = a_{Cn}+a_{Ce}+a_{Cs}+a_{Cw}=\Gamma \Bigg[-\frac{\Delta x}{\delta y}+ \Bigg(-\frac{\Delta y}{\delta x}\Bigg)+\Bigg(-\frac{\Delta x}{\frac{1}{2}\delta y}\Bigg)+\Bigg(-\frac{\Delta y}{\frac{1}{2}\delta x}\Bigg)\Bigg]=-6\Gamma\end{align}

$b_c$ is a zero vector due to the boundary conditions. Repeating this for all cells we get the following matrix of coefficients

\begin{align}\textbf{A}=\Gamma\begin{bmatrix}-6&1&0&1&0&0&0&0&0\\1&-5&1&0&1&0&0&0\\0&1&-6&0&0&1&0&0&0\\1&0&0&-5&1&0&1&0&0\\0&1&0&1&-4&1&0&1&0\\0&0&1&0&1&-5&0&0&1\\0&0&0&1&0&0&-6&1&0\\0&0&0&0&1&0&1&-5&1\\0&0&0&0&0&1&0&1&-6\end{bmatrix}\end{align}
The coefficients at the boundaries are larger since the cells centers are closer to the faces than the interior point(s). This is the exact same matrix I get from my code but the resulting field is incorrect as seen in the figure below

I don't think there is anything wrong with my code because the resulting coefficient matrix for the test case with grid size 3x3 agrees with my hand calculations. I don't want to use ghost cells either, since I want to have the possibility of generalizing the grid. Maybe some of you can shed some light on my problem.

I guess my main question is that is this a valid way of implementing the FVM or have I missed something? I wasn't sure what to include in this post so it might be lacking, but please do ask for more information in that case.

Also, I know that I haven't been thorough in explaining all the nomenclature but as most of it is pretty standard notation I think you'll manage.


问:
Boundary layer and high order finite element methods
When modeling boundary layers, it is typically required to have a highly refined mesh along the boundary. My question is if there has been success with high order finite element methods to overcome the need for highly refined meshes. Can a high order basis function of a finite element successfully capture the steep gradient in a boundary layer? Can anyone point me to examples? Thanks.</p>


问:
What is an integrated modelling tool?
In this document they describe their code as an integrated modelling tool. I am trying to understand how this is different from a regular modelling tool? Most codes include different modules that do different tasks, but I have never heard anyone refer to this as an integrated code. What does a code have to do to become an integrated code as opposed to a regular code?</p>

Title: What is an integrated modelling tool?


问:
Calculating magnetic flux density
I have a geometry where current density distribution is constant. I can calculate the $z$ component of magnetic flux density according to Biot-Savart law as following:

$$B_z(x,y,z) = \int\limits_x\int\limits_y\int\limits_z
\frac{(y-y')J_x(x',y',z')-(x-x')J_y(x',y',z')}{((x-x')^2+(y-y')^2+(z-z')^2)^{\frac{3}{2}}}dx'dy'dz'\, .$$

I can calculate this integral numerically using MATLAB on various points $(x,y,z)$. My question is, how do we handle infinities which occur when $(x,y,z)=(x',y',z')$.

Title: Calculating magnetic flux density


问:
Does non-dimensionalizing ODEs affect the stiffness of the equations?
Does non-dimensionalizing ODEs affect the stiffness of the equations? Can it improve the stability of numerical methods like ode45,ode113 in MATLAB?

I am trying to solve 2 eqns. which might involve very tiny step size to get stability like 1e-11. I want to improve the speed of the numerical method. Can non dimensionalizing help in this regard?

The ODEs (I.V: Rbar, N) are shown in link 1. The output should be close to that shown in 2

My working code for a test problem:

```<code>clc
clear
close all
tic

T  = 273 + 160 ;
Xo = 0.06;

[SS,delGv,rc,Z,rstar,delGstar,I_nu] = nucleate2(Xo, T, 0);

t  = [1 1e7];
y0 = [rstar,1];
options = odeset('InitialStep',1e-3,'MaxStep',100);
[t1,y1] = ode23(@(t,y)NOCmodel2(t,y,T),t,y0,options);

figure(1)
set(gcf, 'Position',  [100, 100, 900, 800])
loglog(t1,y1(:,1),'-b');
xlabel ('time(s)')
ylabel ('<span class="math-container">$R_{mean}$</span>','Interpreter','latex')
set(findall(gcf,'-property','FontSize'),'FontSize',20,'fontweight','bold')

figure(2)
set(gcf, 'Position',  [100, 100, 900, 800])
loglog(t1,y1(:,2),'-r');

xlabel ('time(s)')
ylabel ('<span class="math-container">$logN (/mm^3)$</span>','Interpreter','latex')
set(findall(gcf,'-property','FontSize'),'FontSize',20,'fontweight','bold')

figure(3)
set(gcf, 'Position',  [100, 100, 900, 800])
semilogx(t1,y1(:,2),'-r');
xlabel ('time(s)')
ylabel ('<span class="math-container">$N (/mm^3)$</span>','Interpreter','latex')
set(findall(gcf,'-property','FontSize'),'FontSize',20,'fontweight','bold')

toc

function [yprime] = NOCmodel2(t,y,T)
Rm = y(1);
N  = y(2);

%A. Known variables
Xp         = 1;
Xo         = 0.06;
R          = 8.314 ;                             %Gas constant
kB         = 1.38064852 * 1e-23 ;                %R/avogadro no.
Na         = 6.02214 * 1e23;                     %avogadro no.
IncXe      = 0.01;                               %Equilibrium mol frac.
sigma      = 0.13;                               %interfacial free energy J/m2
Di         = 5e-20;                              %aDiffusivity
Vinc       = 1.6*1e-29*Na;                       %m3/mol
VFe        = Vinc;                               %Vinc;

%B. Calculated properties

alpha      = VFe/Vinc;
Rzero      = 2*sigma*Vinc/R/T;
IncXieq(1) = IncXe ;
IncXi(1)   = IncXieq(1)*exp(2*sigma*Vinc/(R*Rm*T));

if Rm <0
   error('negative Rm')
elseif N < 0
   error('negative N')
end

% Solute balance
X_O = (Xo - alpha*(4/3)*pi*Rm^3*N*Xp)/(1 - (alpha*(4/3)*pi*Rm^3*N*Xp));

%C. number density
[~,~,Rstar,~,Rcrit,~,C1] = nucleate2(X_O, T, t) ;

f_coars = 0 ;
if Rm < 1.01*Rstar && Rm > 0.99*Rstar
    f_coars = 1-(1000*((Rm/Rstar) - 1)^2) ;
end

D1 = (4/27) * ((IncXieq(1))/(alpha*Xp - IncXieq(1))) * (Rzero*Di/Rm^3) * ((((Rzero*IncXieq(1))/(Rm*(Xp - IncXieq(1)))) * ((3/(4*pi*Rm^3))) - N) - 3*N);

dNdt = abs(min(0,(C1+D1)/abs(C1+D1)))*f_coars*D1 + max(0,(C1+D1)/abs(C1+D1))*C1;

%D. mean radius
A1 = (Di/Rm)*((X_O - IncXi(1))/(alpha*Xp - IncXi(1))); A2 = (1/N)*dNdt*(Rcrit - Rm);

A = A1 + A2 ;
B = (4/27) * ((IncXieq(1))/(alpha*Xp - IncXieq(1))) * ((Rzero)*Di/Rm^2) ;

%fprintf('t = %1.1e, A1 =  %2.2e, A2 =  %2.2e, B =  %2.2e, C =  %2.2e, D =  %2.2e, fc = %2.2f, Rm = %2.2e, dNdt = %2.2e, xo = %2.2e, N = %2.2e, Rcrit = %2.2e\n',t,A1,A2,B,C1,D1,f_coars,Rm, dNdt, X_O,N,Rcrit )

yprime = [f_coars*B + (1-f_coars)*A ; abs(min(0,(C1+D1)/abs(C1+D1)))*f_coars*D1 + max(0,(C1+D1)/abs(C1+D1))* C1];

end

function [SS,delGv,rc,Z,rstar,delGstar,I_nucl] = nucleate2(Xo, T, t)
R         = 8.314 ;                             %Gas constant
kB        = 1.38064852 * 1e-23 ;                %R/avogadro no.
Na        = 6.02214 * 1e23;                     %avogadro no.
IncXe     = 0.01;                               %Equilibrium mol frac.
sigma     = 0.13;                               %interfacial free energy J/m2
Di        = 5e-20;                              %aDiffusivity
Vinc      = 1.6*1e-29*Na;
VFe       = Vinc;                               %Vinc;
a         = 4.04e-10 ;                          %lattice parameter FCC iron

Vm1       =  Vinc;                              %m3/mol
Vm        =  Vm1/Na;                            %m3
SS        =  Xo/IncXe;                          %reaction quotient of bulk steel
delGv     = -R*T*log(SS)/Vm1;                   %chemical driving force (J/m3)
rc        = -2*sigma/delGv ;                    %critical radius (m)
Z         =  Vm/(2*pi*rc^2)*sqrt(sigma/kB/T);   %Zeldovich factor
rstar     =  rc + 0.5*sqrt(kB*T/(pi*sigma)) ;   %(1/(Z*sqrt(pi))); %m
delGstar  =  16*pi*(sigma)^3/(3*delGv^2) ;      % J

betastar  = 4*pi*rc^2*Di*Xo/a^4;
tow       = 0.5/(betastar*Z);

PreEx     = betastar * Z * Na/VFe;
I_nucl    = PreEx * exp(-delGstar/kB/T) * exp(-tow/t);

end
</code>```

}


问:
CVXOPT intermediate step valuation stepping out of function domain of defintion
I am using CVXOPT, particularly to solve a nonlinear convex optimization problem. Either the objective function or the constraints involve some functions that are only defined in a strict subset of $\mathbf R^n$. For example, the function is the logarithm $\ln$. It seems that the package may well use symmetric primal dual method which allows infeasible points during the iteration. However, that would generate domain error from the functions. Does extending function to the whole of the Euclidean space help? For the example of the logarithm, I extend it so that $\ln(x)=-\infty$ or $-M$, and $\frac1x=\infty$ or $M$ $\forall x\le0$ for some very large positive $M$. Is this the practice to get around the problem? What are the more efficient ways to handle the problem?

I present below the problem I am reduced from the actual problem I am dealing with. The size of the latter is much bigger than the former.
$$z:=\begin{bmatrix}x \\ y\end{bmatrix}.$$
\begin{align}
\min_{z\in \mathbf R^{+2}} &\frac12 (z-z_0)^TA(z-z_0) \\
& g(x,y):=xa^{\frac y x b}\le c  \\
&\Longleftrightarrow f(x,y):= x\ln\frac xc+y\, b\ln a<0,  \\
&z >0
\end{align}
for some given $2\times 1$ matrix $z_0$ and positive $a, b, c$. Both $g$ and $f$ are convex with respect to $z$. So this is a convex optimization problem.

I am using the Python package CVXOPT to solve this problem. I currently adopt the logarithmic formulation for the constraint, i.e. $f$. CVXOPT requires me to supply the gradient and Hessian of the constraint. I list the nonzero partial derivatives up to the second order below.
\begin{align*}
    \frac{\partial f}{\partial x} &= \ln\frac xc+1,   \\
    \frac{\partial f}{\partial y} &= b\ln a,  \\
    \frac{\partial^2 f}{\partial x^2} &= \frac1x.
\end{align*}


问:
C.S.R method in finite element matrix assembly
I have solved the 2D Poisson equation using finite element method with simplex triangular element in MATLAB.

First, I generated the triangular mesh using pdetool in Matlab and it gives me three matrices: p (position of each node), t (element with three vertices) and e. Then, I derived the local stiffness matrix like this,

for I=1:i_elems

  conec(i_elem,1)=t(1,i_elem);
  conec(i_elem,2)=t(2,i_elem);
  conec(i_elem,3)=t(3,i_elem);

  node_1=conec(i_elem,1);
  node_2=conec(i_elem,2);
  node_3=conec(i_elem,3);

    NN_e(1,1)=(1/6) * elem_area(i_elem);
    NN_e(1,2)=(1/12)* elem_area(i_elem);
    NN_e(1,3)=(1/12)* elem_area(i_elem);
    NN_e(2,1)=(1/12)* elem_area(i_elem);
    NN_e(2,2)=(1/6) * elem_area(i_elem);
    NN_e(2,3)=(1/12)* elem_area(i_elem);
    NN_e(3,1)=(1/12)* elem_area(i_elem);
    NN_e(3,2)=(1/12)* elem_area(i_elem);
    NN_e(3,3)=(1/6) * elem_area(i_elem);

% and then I assembled that like this,

    NN(node_1,node_1)=NN(node_1,node_1)+NN_e(1,1);
    NN(node_1,node_2)=NN(node_1,node_2)+NN_e(1,2);
    NN(node_1,node_3)=NN(node_1,node_3)+NN_e(1,3);
    NN(node_2,node_1)=NN(node_2,node_1)+NN_e(2,1);
    NN(node_2,node_2)=NN(node_2,node_2)+NN_e(2,2);
    NN(node_2,node_3)=NN(node_2,node_3)+NN_e(2,3);
    NN(node_3,node_1)=NN(node_3,node_1)+NN_e(3,1);
    NN(node_3,node_2)=NN(node_3,node_2)+NN_e(3,2);
    NN(node_3,node_3)=NN(node_3,node_3)+NN_e(3,3);

end
```

NN matrix will be a sparse matrix. I just want to use the c.s.r technique and make the assembling process more efficient but I don't know how to assemble the global stiffness matrix in this way. I would be grateful if someone could help me.


问:
Providing Jacobian as LinearOperator in scipy.optimize.root
I asked this question a few days ago on stackoverflow, but I figure scicomp.stackexchange is probably a better place. Sorry for the double post.

I want to solve a system of nonlinear equations using scipy.root. For performance reason, I want to provide the jacobian of the system using a LinearOperator. However, I cannot get it to work. Here is a minimal example using the gradient of the Rosenbrock function, where I first define the Jacobian (i.e. the Hessian of the Rosenbrock function) as a LinearOperator.

import numpy as np
import scipy.optimize as opt
import scipy.sparse as sp

ndim = 10

def rosen_hess_LO(x):
    return sp.linalg.LinearOperator((ndim,ndim), matvec = (lambda dx,xl=x : opt.rosen_hess_prod(xl,dx)))

opt_result = opt.root(fun=opt.rosen_der,x0=np.zeros((ndim),float),jac=rosen_hess_LO)
```
<p>Upon execution, I get the following error :</p>
```TypeError: fsolve: there is a mismatch between the input and output shape of the 'fprime' argument 'rosen_hess_LO'.Shape should be (10, 10) but it is (1,).
```

What am I missing here ?

EDIT : A partial answer to this question can be found here.


问:
Numerov method for solving Schrödinger equation
I have just begun learning computer science to apply it to Physics and I am trying to write a code for solving Schrödinger's equation of the harmonic oscillator (setting $V=\frac{x^2}{2}$) in one dimension.

I know the basics of how to do this, but I would like to know whether some resources are available, or if someone could help me how to write my program (I feel more confident with Python, but C, C++, or any other language works fine too).

I know I should:

Choose a value for the maximum angular momentum $\ell$ and the quantum number $n$.

Seek an energy interval $E_i, E_{i-1}$ where I can look more accurately for the solution.

Find the best approximate value of energy through an algorithm (ex: secant formula).

Compare the interval in energies with a certain accuracy $\epsilon$ established at the beginning.

Thank you in advance for your help!


问:
Spline interpolation for vector-valued data in 3D space
I have output from a 3D linear elasticity finite element simulation which uses linear tetrahedral elements, such that the displacement is continuous over the nodes but the gradient is not ($C_0$ continuous). I wish to "upgrade" the displacement data so that it can be reported as a $C_1$ continuous quantity which I can directly differentiate to produce the gradient, which would yield a consistency that I find appealing. Instead of directly interpolating gradient data over the nodes (or using a higher order basis function, which of course would make this exercise moot), I wish to interpolate the primary variable, displacement, as a quadratic polynomial with splines.

For 1D quadratic splines, each element or segment takes the form:
$$f_e(x)=a_e+b_ex+c_ex^2$$
$$\nabla f_e(x)=b_e+2 c_ex$$
To form the system of equations, the function value and gradient are matched at end points. And since this is 1D, each end point connects to at most one other segment.
$$f_e(x_i)=f_{e+1}(x_i)$$ and $$\nabla f_e(x_i)=\nabla f_{e+1}(x_i)$$
Such is the predictable nature of this arrangement that you are guaranteed to have $3(n-1)$ variables and at least $n+2(n-1)=3n-2$ equations. This can be supplemented with an additional assumption or solved in a linear least squares fashion. Awesome.

For the 3D case, it looks similar:
$$f_e(x)=a_e+x^Tb_e+x^Tc_ex$$ and $$\nabla f_e(x)=b_e+2c_ex$$where $x$ is a $3$ x $1$ vector, $a_e$ is a scalar, $b_e$ is a $3$ x $1$ vector, and $c_e$ is a $3$ x $3$ matrix. Because $x^Tc_ex$ is a quadratic form, only the symmetric part of $c_e$ is of any consequence, so it can be safely assumed $c_e$ is a symmetric matrix with 6 independent variables:

$$f_e(x)=a_e+x^Tb_e+\begin{pmatrix}x_0x_0 & x_0x_1 & x_0x_2 & x_1x_1 & x_1x_2 & x_2x_2\end{pmatrix}^T\begin{pmatrix}c_{e00} & c_{e01} & c_{e02} & c_{e11} & c_{e12} & c_{e22}\end{pmatrix}$$
$$\nabla f_e(x)=b_e+\begin{pmatrix}2x_0 & x_1 & x_2 & 0 & 0 & 0\\\ 0 & x0 & 0 & 2x_1 & x_2 & 0\\\ 0 & 0 & x_0 & 0 & x_1 & 2x_2\end{pmatrix}^T\begin{pmatrix}c_{e00} & c_{e01} & c_{e02} & c_{e11} & c_{e12} & c_{e22}\end{pmatrix}$$

The hard part is that each node is not limited to 2 neighboring elements, in some cases it can be upward of 20! In one very modest case (I usually work with much finer grids), I have approximately 2000 tetrahedra and 500 nodes, with 51,000 (!) combinations where I must enforce gradient continuity. For example: if two tetrahedra (T1, T2) share 3 nodes (A, B, C), I must enforce gradient continuity for at A, B and C using the spline variables for T1 and T2.

This yields a sparse matrix ($10^{-4}$ sparsity) with $4t+3*51,000 = 161,000$ equations and $10t=20,000$ variables, and a $161,000$ x $3$ right-hand side (each component of displacement is independently interpolated). This is obviously way more complicated than calculating the original finite element system, which seems counterintuitive to me. Nevertheless, it does work...slowly. In the case above, I solved the original finite element system in 500 milliseconds and the interpolation took 30 seconds. I didn't even have the memory needed (it blew up to 120 gigs or something silly like that) to interpolate my other grids which takes a whopping 4 seconds to solve per time step. I am using SuiteSparse's SPQR solver, which produces a least squares solution using a QR factorization of the giant sparse system.

Questions:

What do you call this interpolation procedure? Splines? Super-splines? B-Splines? I can't make heads or tails of what I read in literature on anything beyond 1D splines.

Am I "double counting" some conditions and making the system overly complicated? I thought I could ignore a few conditions, such as pairs of tetrahedra that touch at one point only, but ignoring even that destroys the interpolation. I am dumbfounded by how much more complicated this interpolation procedure is than the original finite element calculations.

Any ideas on how to simplify this system? I am confident that the algorithms I am using are top notch -- I'm just giving them too much work.

Examples plotting pressure (trace of the gradient) without interpolation:

and with the procedure above:


问:
Differences in using Clausius-Duhem inequality vs Principle of Virtual Work/Power in derriving constitutive equations?
I am a novice getting my toes wet in continuum mechanics and nonlinear elasticity. I have seen papers that use both approaches in developing constitutive connections to compliment balance equations of various sorts. My question is the energy inequality more general than the virtual work approach? I am tackling a visco-elastic problem where I will need to account for the time dependent relaxation of the material. Can this not be handled by virtual work since energy will be lost rather than conserved?

I ask because virtual work seems a bit more straight forward but I wonder if it is at the expense of applicability. I am looking for the pros and cons of both approaches.

Thanks in advance!


问:
Difference between wave vector and density matrix in numerical calculation of Schrödinger equation
I solved Schrödinger equation for a following tow-level time-dependent Hamiltonian numerically in two ways:

```<code>import numpy as np

def H(t):
    return np.array([[t,0.5],[0.5,-t]],dtype="complex64")
</code>```

First, the state is treated as wave vector:

```<code>from scipy.integrate import solve_ivp
def schrodinger(t, X):
    dXdt = -1j * H(t).dot(X)
    return dXdt

T=80.
t_list = np.arange(-T, T, 2*T / 1000.0)
init_state= np.array([1.,0.],dtype="complex64")
solution = solve_ivp(fun=lambda t, X: schrodinger(t, X), t_span=[-T,T], y0=init_state, t_eval=t_list, method="RK45", vectorized=True)
</code>```

Second, the state is treated as (vectorized) density matrix:

```<code>def schrodinger_rho(t, X):
    unit = np.eye(2, 2, dtype="complex64")
    dXdt = -1j * (np.kron(unit, H(t)) - np.kron(H(t).T, unit)).dot(X)
    return dXdt
T=80.
t_list = np.arange(-T, T, 2*T / 1000.0)
init_state= np.array([1.,0.,0.,0.],dtype="complex64")
solution_rho = solve_ivp(fun=lambda t, X: schrodinger_rho(t, X), t_span=[-T,T], y0=init_state, t_eval=t_list, method="RK45", vectorized=True)
</code>```

The numerical results are as follows:

Here is a question : why cannot the state preserve its unitarity, say $P_0+P_1=1$, when the state is treated as wave vector? The larger the value of T, the stronger this tendency of violation becomes. Moreover, the numerical results become worse when I use method="BDF".

I want to treat the state as wave vector because of the numerical cost. Is there any way to improve this phenomenon?


问:
1D Poisson equation and quadratic basis functions assembly
I'm solving the simple Poisson problem $$-u''(x)=1$$ in the interval $[0,1]$ with $u(0)=u(1)=0$.

I discretised my domain as done here, i.e. with

nodes = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]
and
elements = [[0, 1, 2], [2, 3, 4], [4, 5, 6], [6, 7, 8]]

I'm trying to assemble the matrix and the r.h.s by looping over the elements, and I found the same results of here (see pag. 22 for the matrix and load vector)

$$A_K=\frac{1}{3h}
\begin{bmatrix}
7 & -8 & 1 \\
-8 & 16 & -8 \\
1 & -8 & 7
\end{bmatrix}
$$

and $$F_K=\frac{h}{6}
\begin{bmatrix}
1 \\
4 \\
1
\end{bmatrix}$$

This simple codes implements the idea, but my numerical solution is scaled by a factor of $(2h)$: if you set u/(2*h) instead of u, you'll get the right solution exact at grid points).

I think the source of error should be in how I define the step $h$, because those midpoints are confusing me a bit, but I can't figure out how to fix this

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

a = 0.0
b = 1.0
N = 8
h = (b-a)/(N)

nodes = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]
elements = [[0, 1, 2], [2, 3, 4], [4, 5, 6], [6, 7, 8]]

def assembly_stiff_quadratic():
    A = np.array([[7,-8,1],[-8,16,-8],[1,-8,7]])
    return (A)/(3.0*(h))

def assembly_load_quadratic():
    return ((h)/6.0)*np.array([1,4,1])

numK = len(elements)
numV = len(nodes)
A = np.zeros([numV,numV])
f = np.zeros(numV)

for K in range (0,numK):
    loc2glb = elements[K]
    AK = assembly_stiff_quadratic()
    print(K,loc2glb)
    bK = assembly_load_quadratic()
    A[np.ix_(loc2glb,loc2glb)] += AK
    f[np.ix_(loc2glb)]+=bK

f = f[1:-1]
A = A[1:-1,1:-1]
u = np.linalg.solve(A,f)

def sol(x):
    return -(0.5)*x*(x-1)

x = np.linspace(0,1,N+1)
plt.plot(x[1:-1],u,'o')
plt.plot(x[1:-1],sol(x[1:-1]),'r')
</code>```

My question is not about debugging the code, rather about how to do the assembly in this case. I posted the code because I think it is straightforward to read and easy to test.

EDIT:

After the comments, I wrote the following snippet which produces the following order of convergence plot

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

a = 0.0
b = 1.0
def solvePoisson(N):
    h = 2*(b-a)/(N)

    elements=[]
    for i in range(0,int(N-1)):

        if i%2==0: elements.append([j for j in range(i,i+3)])
    nodes= np.linspace(a,b,N+1)

    def phi0(x):
        return (1-x)*(1-2*x)

    def phi1(x):
        return 4*x*(1-x)

    def phi2(x):
        return x*(2*x-1)

    def rhs(x):
        return np.sin(x)

    def assembly_stiff_quadratic():
        A = np.array([[7,-8,1],[-8,16,-8],[1,-8,7]])
        return (A)/(3.0*(h))

    numK = len(elements)
    numV = len(nodes)
    A = np.zeros([numV,numV])
    f = np.zeros(numV)

    for K in range (0,numK):
        loc2glb = elements[K]
        AA = nodes[loc2glb[0]]
#        BB = nodes[loc2glb[2]]
        AK = assembly_stiff_quadratic()
        bK = h*np.array([integrate.quadrature(lambda x: rhs(AA+h*x)*phi0(x),0,1)[0],integrate.quadrature(lambda x: rhs(AA+h*x)*phi1(x),0,1)[0],integrate.quadrature(lambda x: rhs(AA+h*x)*phi2(x),0,1)[0]])

        A[np.ix_(loc2glb,loc2glb)] += AK
        f[np.ix_(loc2glb)]+=bK

    f = f[1:-1]
    A = A[1:-1,1:-1]
    u = np.linalg.solve(A,f)
    return u

def sol(x):
    return np.sin(x) - x *np.sin(1)

N = 8
x = np.linspace(0,1,N+1)
plt.plot(x[1:-1],solvePoisson(N),'o')
plt.plot(x[1:-1],sol(x[1:-1]),'r')
</code>```


问:
Results blow up when number of intervals is increases (Yee algorithm FDTD, dielectric sphere)
I have been trying to write a program that analyses EM wave scattering by a dielectric sphere for a project.

The reference is Sadiku's book Numerical Methods in electromagnetics Edition 3.

Now the book has code in Matlab but I have to use Python.

I rewrote the code in Python and ran it, but the values Ey/Einc and Ez/Einc blow up, much more so when I increase number of intervals.
Even at n = 50 they do not match the results.

Here is my code

```<code>'''

import numpy as np
from matplotlib import pyplot as py
import warnings

'''
The followings constants are defined globally, as they are used in the code.
All are in SI units.
'''

pi = np.pi
c = 3e8
e0 = 8.85e-12
mu0 = 4*pi*1e-7

'''
DOCSTRING for Taflove(er, mu, sigma, nu, num)
This funciton is used to calculate some parameters and Taflove's constants.
The inputs are relative permittivity and permeability, conductivity,
frequency of incoming wave and grid number.
Three tuples are returned.
First tuple contains lambd, delta and delta_t parameters.
Second tuple contains R constants.
Third tuple contains C constants.
'''

def Taflove(er, mur, sigma, nu, num):

    lambd = c/nu
    delta = lambd/num
    delta_t = delta/(2*c)

    R = delta_t/e0
    R_a = (c*delta_t/delta)**2
    R_b = delta_t/(mu0*delta)

    C_a = (1 - R*(sigma/er))
    C_b = R_a/er

    return np.array((lambd, delta, delta_t)), np.array((R, R_a, R_b)), np.array((C_a, C_b))

'''
DOCSTRING for
'''

er = (1, 4)          # Relative Permittivity
mur = (1, 1)         # Relative Permeability (1 for non-magnetic)
sigma = (0.1, 0)       # Conductivity (0 for perfect dielectric/lossless medium)
   # Conductivity at boundary is taken to be non zero to prevent abrupt results
nu = 2.5e9      # Frequency of EM wave = 2.5 GHz (for Microwave)
num = 40        # Due to computing power
radius = 4.5e-2 # Radius of sphere is 4.5 cm

P_list, R_list, C_list1 = Taflove(er[0], mur[0], sigma[0], nu, num)
print(P_list)
print(R_list)
print(C_list1)
input('check1')
_, _, C_list2 = Taflove(er[1], mur[1], sigma[1], nu, num)
C = [C_list1, C_list2]
R, R_a, R_b = R_list
lambd, delta, delta_t = P_list
#1 = air
#2 = sphere
#E = np.array([0, 0, 1])                 # Incoming EM wave is linearly polarized in z direction
#B = np.array([1, 0, 0])/c               # B = k cap cross E/c (in x direction)
#k = 2*pi*np.array([0, 1, 0])/P_list[0]  # k = 2pi/lambda in direction of wave (in y direction)
Size = (19, 39, 19)                     # Grid/Lattice size. More in y because wave is propagating in y driection
Center = (19.5, 20, 19)                 # Center of sphere
WMax = 2                                # Steps in one program looop
TMax =20*5                            # Number of Timesteps
Media = 2                               # Number of Media
JPW = 3                                 # y coordinate of wave

I, J, K = Size
Ex = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Hx = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Ey = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Hy = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Ez = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
Hz = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2, WMax+1))
x_ = np.linspace(0, I+1, I+2)
y_ = np.linspace(0, J+1, J+2)
z_ = np.linspace(0, K+1, K+2)
print(x_)
print(y_)
print(z_)
radius_unit = radius/P_list[1]
print(radius_unit)
input('check')

pos_func = lambda R : np.sqrt(sum(np.array(R) - np.array(Center))**2)

XMed = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2))
YMed = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2))
ZMed = np.zeros((Size[0]+2, Size[1]+2,Size[2]+2))
for i in range(Size[0]+2):
    for j in range(Size[1] + 2):
        for k in range(Size[2]+2):
            if pos_func((i + 0.5, j, k))<radius_unit:
                XMed[i, j, k] = 2
            else:
                XMed[i, j, k] = 1

            if pos_func((i, j+ 0.5, k))<radius_unit:
                YMed[i, j, k] = 2
            else:
                YMed[i, j, k] = 1

            if pos_func((i, j, k+ 0.5))<radius_unit:
                ZMed[i, j, k] = 2
            else:
                ZMed[i, j, k] = 1

print(XMed)
Ey1 = np.zeros(J+2)
Ez1 = np.zeros(J+2)
print(len(Ey1))

print("Setting up completed...")
input('a')
###########################
t3 = 2
t2 = 1
t1 = 0
for n in range(TMax):
    if n%10 == 0:
        print(str(100*n/TMax) + " % complete.")

    #Applying soft truncation conditons
    for k in range(K+1):
        for j in range(J+1):
            for i in range(I+1):
                if i == 0:
                    if k != K and k!=0:
                        Hy[0, j, k, t3] = (Hy[1, j, k-1, t1] + Hy[1, j, k, t1] + Hy[1, j, k+1, t1])/3
                        Hz[0, j, k, t3] = (Hz[1, j, k-1, t1] + Hz[1, j, k, t1] + Hz[1, j, k+1, t1])/3
                    elif k == K:
                        Hy[0, j, k, t3] = (Hy[1, j, k-1, t1] + Hy[1, j, k, t1])/2
                        Hz[0, j, k, t3] = (Hz[1, j, k-1, t1] + Hz[1, j, k, t1])/2
                    else:
                        Hy[0, j, k, t3] =(Hy[1, j, k, t1] + Hy[1, j , k+1, t1])/2

                        Hz[0, j, k, t3] = (Hz[1, j, k, t1] + Hz[1, j, k+1, t1])/2
                if j == 0:
                    Ex[i, 0, k, t3] = Ex[i, 1, k, t1]
                    Ez[i, 0, k, t3] = Ez[i, 1, k, t1]
                elif j == J:
                    Ex[i, j, k, t3] = Ex[i, j-1, k, t1]
                    Ez[i, j, k, t3] = Ez[i, j-1, k, t1]
                if k == 0:
                    if i != I and i!= 0:
                        Ex[i, j, 0, t3] = (Ex[i-1, j, 1, t1] + Ex[i, j, 1, t1] + Ex[i+1, j, 1, t1])/3
                        Ey[i, j, 0, t3] = (Ey[i-1, j, 1, t2] + Ey[i, j, 1, t2] + Ey[i+1, j, 1, t1])/3
                    elif i == 0:
                        Ex[0, j, 0, t3] = (Ex[i, j, 1, t1] + Ex[i+1, j, 1, t1])/2
                        Ey[0, j, 0, t3] = (Ey[i, j, 1, t1] + Ey[i+1, j, 1, t1])/2
                    else:
                        Ex[i, j, 0, t3] = (Ex[i, j, 1, t1] + Ex[i-1, j, 1, t1])/2
                        Ey[i, j, 0, t3] = (Ey[i, j, 1, t1] + Ey[i-1, j, 1, t1])/2
                #Applying Yee's algorithm
                Hx[i, j, k, t3] = Hx[i, j, k, t2] + R_b*(Ey[i, j, k+1, t2] - Ey[i, j, k, t2]\
                                                         + Ez[i, j, k, t2] - Ez[i, j+1, k, t2])

                Hy[i, j, k, t3] = Hy[i, j, k, t2] + R_b*(Ex[i, j, k, t2] - Ex[i, j, k+1, t2]\
                                                         + Ez[i+1, j, k, t2] - Ez[i, j, k, t2])

                Hz[i, j, k, t3] = Hz[i, j, k, t2] + R_b*(Ex[i, j+1, k, t2] - Ex[i, j, k, t2]\
                                                         + Ey[i, j, k, t2] - Ey[i+1, j, k, t2])
                if k == K:
                    Hx[i, j, k, t3] = Hx[i, j, k-1, t3]
                    Hy[i, j, k, t3] = Hy[i, j, k-1, t3]
                if j !=0 and j !=J and k !=0:
                    M = int(XMed[i, j, k] - 1)
                    Ex[i, j, k, t3] = C[M][0]*Ex[i, j, k, t2] + (C[M][1]/R_list[2]) *(Hz[i, j, k, t3]\
                                      -Hz[i, j-1, k, t3] + Hy[i, j, k-1, t3] - Hy[i, j, k, t3])

                if k!=0:
                    M = int(YMed[i, j, k] - 1)
                    if i!=0:
                        Ey[i, j, k, t3] = C[M][0]*Ey[i, j, k, t2] + (C[M][1]/R_list[2]) *(Hz[i-1, j, k, t3]\
                                      -Hz[i, j, k, t3] + Hx[i, j, k, t3] - Hx[i, j, k-1, t3])
                        if Ey[i, j, k, t3]!=0:
                            #print(Ey[i, j, k, t3])
                            #input(str(i) +','+ str(j)+',' + str(k))
                            pass
                    else:
                        Ey[i, j, k, t3] = C[M][0]*Ey[i, j, k, t2] + (C[M][1]/R_list[2]) *(0\
                                      -Hz[i, j, k, t3] + Hx[i, j, k, t3] - Hx[i, j, k-1, t3])
                if j !=0 and j!=J:
                    M = int(ZMed [i, j, k])-1
                    if M == 0:
                        Cam = 1
                    else:
                        Cam = C[M][0]
                    if i !=0:
                        Ez [i, j, k,t3] = Cam*Ez[i, j, k, t2]\
                                          + (C[M][1]/R_list[2]) *(Hy[i, j, k, t3]\
                                      -Hy[i-1, j, k, t3] + Hx[i, j-1, k, t3] - Hx[i, j, k, t3])

                    else:
                        Ez [i, j, k,t3] = Cam*Ez[i, j, k, t2]\
                                          + (C[M][1]/R_list[2]) \
                                          *(Hy[i, j, k, t3]\
                                      -0 + Hx[i, j-1, k, t3]- Hx[i, j, k, t3])

                    if j == JPW:
                        Ez[i,j,k, t3] = Ez[i,j,k,t3] + np.sin(2*pi*nu*delta_t*(n+1))
                if k==K:
                    Ex[i, j, k + 1, t3] =  Ex[i, j, k , t3]
                    Ey[i, j, k + 1, t3] =  Ey[i, j, k , t3]

            if k == K and n>TMax-num-1:
                temp = abs(Ey[I, j, K-1, t3])
                if temp>Ey1[j]:
                    Ey1[j] = temp#/(Ey1[j]+0.0000001)
                temp = abs(Ez[I, j, K, t3])
                if temp > Ez1[j]:
                    Ez1[j] =temp#/(Ez1[j]+1e-10)
    t1 = int(t2 + 0.)
    t2 = int(t3 + 0.)
    t3 = int(t3%3)

Integr = np.linspace(1, 100, 100)
print(Ey1)
for i in range(len(Ey1)):
    if Ey1[i]<1:
        py.plot(Integr[i], Ey1[i], 'r*')
py.xlabel('j')
py.ylabel('Computer Ey/Einc ')
py.grid()
py.show()

print(Ez1)
for i in range(len(Ez1)):
    if Ez1[i]<1:
        py.plot(Integr[i], Ez1[i], 'r*')
py.xlabel('j')
py.ylabel('Computer Ez/Einc ')
py.grid()
py.show()

</code>```

Does anyone know what could be the bug?
I have been checking this only for the last week.
Thanks.


问:
Solving Laplace equation with constraint on boundary
I have found the following PDE problem in a paper:

Essentially, we have a rectangular domain where there is a unknown interface $z=\xi(x)$ (liquid-air interface) separating the domain into two medium (air on top for $i=a$ and liquid at the bottom for $i=b$). We need to solve Laplace equation (3) subject to Dirichlet B.C. at $z=0$ and $z=a+1$ and Neumann B.C. at $x=0$ and $x=\pi/k$ (Eq.4-7). At the interface $z=\xi(x)$ we have the usual matching conditions for electric field and potentials. This is just a standard problem of Laplace equation if $z=\xi(x)$ is known. However, $z=\xi(x)$ is unknown, and it has to satisfy an addition equation (10). $\tau_e$ is the electric stress which depends on $\phi_i$. $K$ is a scalar unknown which needs to be determined as well. $z=\xi(x)$ is subject to Neumann B.C. (eq. 13) and an integral constraint (eq.12). So in total we need to solve for $\phi_i(x,z)$, $\xi(x)$ and $K$. The paper just said it is solved by FEM, but I have no idea how to set it up. Are there any general strategy for solving these kind of Laplace equation where the interface is also an unknown?

Update:
added a diagram for better illustration of the problem


问:
Divergence on wave equation simulation
I'm currenly working on my own PDE solver for non-linear simulations in python. I've done succesfully simulations for KdV and Fisher's equation, but now I'm playing with second order derivatives in time.

The equation I'm using is 1D wave equation, with $\Delta x = 0.5$ and $\Delta t = 0.01$. For simplicity I'm using $c = 1$ and solving by defining an auxiliar field $v(x,t)$, so the equation looks like

$$\partial_t u = v(x,t); \qquad \partial_t v = c^2 \partial_{xx} u$$

Currently I'm using my own RK4 integrator for time and finite differences up to $O(h^2)$ for space. My problem is that as time passes by, with $u(x,0) = \sin(\alpha x)$ and $\partial_t u(x,0) = 0$ the solution seems to grow until it numerically diverges. I tried with gaussian initial conditions, such as $u(x,0) = \exp(-\beta x^2)$, $\partial_t u(x,0) = 0$ and $u(x,0) = 0$, $\partial_t u(x,0) = \exp(-\gamma x^2)$and the solutions seem to irradiate waves which are not part of the analytical solution for these cases.

The boundary conditions are not fixed for the gaussian cases, I just tried these initial condions in order to find what's going on.

Does someone faced the same issues? Could it be for the order on finite differences? I'm pretty new at simulations, so anything would be helpful :)

Here are the plots for the situations mentionatres above, respectively.

Title: Divergence on wave equation simulation


问:
What is temporal order of accuracy of the PISO algorithm?
A few Computational Fluid Dynamics (CFD) codes implement the so called PISO (Pressure-Implicit with Splitting of Operators) algorithm for pressure-velocity coupling.

My concern is what is actual temporal order of accuracy of the algorithm?

It probably isn't enough to discretize time-dependent term in transport equations by second order algorithm (like second order Backward-Euler - BDF2), since we still use mass fluxes from the last time-step in the momentum equations?

Any suggestion of the papers that deal with this topic is welcome, or something relevant to this topic from projection/operator-splitting methods that would help to find an answer.


问:
Energy cannot decent during optimization despite non-zero gradient
Assume we have an (at least) 2nd-order differentiable energy $f(x), x\in R^n.$ And $n$ is very big. Mathematically, I think it is impossible to find a point $\bar{x}$ where the energy cannot be decreased along its negative gradient when the gradient is not equal to zero. However, such kind of case can be occasionally found when an energy is optimized numerically, such as using the interior point method. In the optimization software, there is always a return code that describes that "the step size is too small and cannot make progress". Is it because of the numerical issues? Or is it that my mathematic conclusion is wrong? If it is from numerical issues, what kind of issue could lead to this problem?


问:
Large scale nonsymmetric eigenvalue problems in practical applications
I am trying to determine need for solving large scale nonsymmetric eigenvalue problems in both industry and academia.

I am interested in any kind of problem where we cannot assume that matrices are symmetric. This includes the standard eigenvalue problem, the generalized eigenvalue problem as well as the quadratic eigenvalue problem. I am interested in sparse and dense problems.

MINOR EDIT: I am not interested in problems that only require the computation of a few eigenpairs. This excludes cases where, say, a shifted inverse power iteration will suffice.

The ideal answer to this question includes a reference to at least one application where it is absolute necessary to solve a large scale nonsymmetric eigenvalue problem.

The vast majority of existing libraries only support symmetric standard and generalized eigenvalue problems. I think that there are two good reasons for this:

A large number of applications produce symmetric eigenvalue problems.

The nonsymmetric eigenvalue problem is more difficult than the symmetric eigenvalue problem.

LAPACK contains a complete stack for solving dense nonsymmetric and generalized eigenvalue problems, but the routines for computing eigenvectors are based on BLAS level 1 and level 2 functions. ScaLAPACK never completed the support for nonsymmetric eigenvalue problems and the existing routines for computing eigenvectors are either vulnerable to overflow or do not scale at all. The new StarNEig library can compute all the eigenvalues/eigenvectors for nonsymmetric standard and generalized eigenvalue problems. It applies to matrices/matrix pencils that are real and have real/complex eigenvalue/eigenvectors. It can be substantially faster than LAPACK with parallel BLAS as well as ScaLAPACK.

I have isolated some applications where symmetry is lost.

Example 1 In structural engineering we find the quadratic eigenvalue problem $$(\lambda^2 M + \lambda C + K)x = 0.$$
This problem is equivalent to the generalized eigenvalue problems $$Ax = \lambda B x$$ where $$A = \begin{bmatrix} 0 & K \\ K & C \end{bmatrix}, \quad B = \begin{bmatrix} K & 0 \\ 0 & - M \end{bmatrix}, \quad y = \begin{bmatrix} x \\ \lambda x \end{bmatrix}. $$ We see that $A$ and $B$ are symmetric when the mass matrix $M$, the stiffness matrix $K$ and the dampning matrix $C$ are all symmetric. However, the eigenvalue problem for a damped gyroscopic system takes the form $$(\lambda^2 M + \lambda (C + G) + K)x = 0$$ where $G$ is a skew-symmetric matrix. We see that the symmetric matrix $C$ has been replaced with a general nonsymmetric matrix $C + G$. We see that the equivalent generalized eigenvalue problem is no longer symmetric.

I imagine that rotating equipment such as jet engines, gas turbines and propellers can be represented using damped gyroscopic systems, but I do not have the knowledge to say this for a fact.

Example 2 There electrical engineering we find many examples of eigenvalue problems that not symmetric. Two generalized eigenvalue problems can be found in this paper which surveys the available software. The examples are power networks from Europe. A circuit for which the nodal conductance matrix is nonsymmetric due to the presence of a voltage controlled current source is presented in this paper. In this paper, the authors need all the left/right eigenvectors of a standard eigenvalue problem.

My fundamental problem is that I lack the engineering knowledge and the key words necessary to track down papers that deal with large scale nonsymmetric eigenvalue problems. I can build small toy examples that are nonsymmetric, but I cannot say with authority that their larger siblings are remotely relevant in practice. This is where I hope that this community can help me.


问:
Lumped mass matrices for higher-order finite elements for CFD
Given that some of the mass lumping techniques, for example, row-sum lumping does not produce practically viable lumped mass matrices for all the element shapes, what are the techniques used for mass lumping for higher-order elements? I am particularly interested in applications in CFD.

I know that deal.ii has a tutorial on the projection method where the diffusion term is treated implicity. So, mass lumping is not necessary here.
FreeFem uses explicit step for the diffusion term but demonstrates the implementation with P1 elements only.

Are there any libraries that solve the diffusion term explicitly using lumped-mass matrices for the velocity field when LBB-stable combinations using higher-order elements are employed?

I would appreciate any useful references on the topic. Thanks in advance.