问:
How to apply Neumann boundary conditions in Newton's method
Suppose that I have a very long and tedious set of differential equations. After discretization, I can get a mapping $f:\mathbb{R}^N \to \mathbb{R}^N$ such that solutions $\phi =(\phi_0,...,\phi_{N-1})$ to the discretized diffferential equations are exactly the zeros of $f$, i.e., $f(\phi)=0$.

Now say I want to find solutions $\phi$ which satisfy some Neumann boundary equation, e.g.,
$$
\frac{\phi_1-\phi_0}{dx}=1
$$
Of course, I can always define another function
$$
F(\phi)= \left(\frac{\phi_1-\phi_0}{dx} -1, f(\phi) \right):\mathbb{R}^{N} \to\mathbb{R}^{N+1}
$$
And just use Newton's method on $F(\phi)$. However, since the derivative/Jacobian $DF$ is not a square matrix in this case, things are a little more complicated. Is there an easy way to modify Newton's method to incorporate Neumann conditions?

Indeed, in this post, the OP  deals with Dirichlet boundary conditions, e.g., $\phi_0=0$,, by removing row $0$ of the Jacobian $Df$ and replacing it with the Jacobian of $\phi \mapsto \phi_0$, which is just the Dirac delta at $0$, i.e., $(Df)_{0j} \mapsto \delta_{0j}$. This way the "modified" Jacobian is a still a square matrix. However, this doesn't seem to work if I were to replace row $0$ of $Df$ with the Jacobian of $\phi \mapsto (\phi_1-\phi_0)-dx$, i.e., $(Df)_{0j}\mapsto \delta_{1j}-\delta_{0j}$.

I have also read posts, e.g., here, where they substitute the Neumann boundary condition into the differential equation. However, since my $f$ mapping is quite long and tedious (and not just a simple second order derivative), this would be a little cumbersome for me to do by hand.

Are there any simple alternatives to applying Neumann boundary conditions?

EDIIT. Actually, I realized there was a typo in my code. My initial method works.


问:
Best way to solve system of quadratic forms
I have a system of equations that have the following structure. Let $x\in\mathbb{R}^m$ and let $x_k$ be the $k$-th element of $x$. Let $H_k\in\mathbb{R}^{m\times m}$ for $k=1,\ldots, m$. I need to find the $x$ satisfying,
\begin{align}
x_k = x^\top H_k x
\end{align}
for $k=1,\ldots,m$. I know that I can solve this system via fixed point iteration or via Newton's method. Are these the best ways? Are there non-iterative methods?

ContentLicense: CC BY-SA 4.0


问:
Anyway to escape ODEintWarning (scipy)?
I am trying to fit a differential equation to some data and obtain the parameters of the underlying model. This requires me to try out various parameter values, but this often gives me an ODEintWarning. Having a few of these warnings is OK, but having too many will cause my program to crash quite often (I work using Spyder and if I click somewhere on the GUI while the program is running, Spyder will freeze).

Is there a way to detect an ODEintWarning and escape it? In pseudocode, I essentially want something like:

```<code>for parameter in parameters_list:
    try:
        x = odeint(dxdt, x0,t, *parameters)
    except (at the first sign of an) ODEintWarning:
        pass (and go on to the next set of parameters to be tested)
</code>```

Or, is my best option to detect bad parameter values in advance and try not to feed them into ODEint in the first place?

Title: Anyway to escape ODEintWarning (scipy)?


问:
LSMR: update $\|x\|$ cheaply
In Fong, Saunders: LSMR: An iterative algorithm for sparse least-squares problems, section 3.3 one reads

Since only the last diagonal of $\tilde{R}_k$ and the bottom 2x2 part of $\hat{R}_k$ change each iteration, this estimate of $x_k$ can again be updated cheaply. The pseudocode, omitted here, can be derived as in section 3.1.1.

As a consequence of the omission, all LSMR codes I've found compute $\|x\|$ explicitly (and not cheaply) via a dot product.

Can anyone fill in the details here?

Title: LSMR: update $\|x\|$ cheaply


问:
Split Operator 2D-->3D
I am trying to modify my simulations on population dynamics using the split-operator method from 2D two-level real Hamiltonain $$H_{2D} =T(x,y)\otimes1_2 +\begin{pmatrix} -z & y\\
y & z\\
\end{pmatrix},$$ to a complex 3D two level Hamiltonian
$$H_{3D} =T(x,y,z)\otimes1_2 + \begin{pmatrix} -z & x+iy\\
x-iy & z\\
\end{pmatrix}.$$ To simplify things for myself, I've thrown away the potential contributions and I am just focusing on the kinetic terms ensuring that norm/energy conservation is good.
Loosely speaking, the split operator method follows as
$$e^{-iH\Delta t} \approx e^{-iT\Delta t}e^{-iV\Delta t} +\mathcal{O}(\Delta t)^2$$ and for our case $e^{-iH\Delta t} \approx e^{-iT\Delta t}$, with no error as it can be shown that $\text{Error } =[T,V] = 0$.
In the 2D case:

```<code>import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import scipy.special as scl
import numpy.matlib as mat
import scipy.fftpack as fft
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.widgets import Slider, Button

##Split-Operator ###

# Constants

ω = np.array([1,1]) #Frequency for each coordinate

gs = np.array([0,0]) #Inital Wavepacket shifts

g = 0.0825 # g-vector:(g_j-g_i)/2

h = 0.0430

s_y = 0# s_y vector (g_jy+g_iy)/2

s_z = 0.125 # s_y vector (g_jz+g_iz)/2

#Location of conical intersection: (0,0,0)--> you can translate to coordinates in article if desired

### If you want to switch what surface WP is propagating on, change here###
iS = 0 # Intial starting state

y0c = gs[1]

z0c = gs[0]

##intial Momenta

kIz = 0

kIy = 0

##Set up a grid #Lets see how this goes

My = 128*2

Mz = 128*2

#Number of states

N = 1

#Number of time steps

Tsteps = 600

dt = 0.005

# Grid Lengths

Ly = 10

Lz = 10

LyT = Ly*2

LzT = Lz*2

#Grid of M points

y0 = np.linspace(-Ly,Ly,My)

z0 = np.linspace(-Lz,Lz, Mz)

#Parameters
#k0[1xM] = Grid of M momenta points from 0->L

k0y = np.linspace(-My*np.pi/LyT,My*np.pi/LyT-2*np.pi/LyT,My)

k0z = np.linspace(-Mz*np.pi/LzT,Mz*np.pi/LzT-2*np.pi/LzT,Mz)

##Properties

##Postion and momenta grids##

y0op = (np.tile(y0,(Mz,1))).T

z0op = np.tile(z0,(My,1))

k0yop = (np.tile(k0y,(My,1))).T

k0zop = (np.tile(k0z,(Mz,1)))

##Inital wavefunction: 2D gaussian
ψ_0 = np.zeros((N,My*My),dtype = 'complex')
σ_x = np.sqrt(2/ω[0])
σ_y = np.sqrt(2/ω[1])
σ_z = np.sqrt(2/ω[1])
temp = np.tile(np.exp(-((z0-z0c)/(σ_z))**2)*np.exp(1.j*kIz*z0),(My,1))

temp_y = np.tile(np.exp(-((y0-y0c)/σ_y)**2)*np.exp(1.j*kIy*y0),(Mz,1)).T

temp_1 = temp_y*temp

ψ_0[iS,:] = temp_1.reshape(1,My*Mz)##Gaussian wavepacket

##
ψ_0[iS,:] = ψ_0[iS,:]/np.sqrt(ψ_0[iS,:]@ψ_0[iS,:].T)##Normalized wavefunction

#Now we need to calculate the propagators

##Kinetic T = exp(i*hbark^2.2mdt/hbar)

TP = np.exp(-1.j*(np.tile((k0y**2)/2,(Mz,1)).T+np.tile((k0z**2)/2,(My,1)))*dt)

T =  np.tile((k0y*k0y)/2,(Mz,1)).T + np.tile((k0z*k0z)/2,(My,1))

ψ_0

##Potential Energy propagator V

ψ = ψ_0 #Initialization of wavepacket

tR = np.linspace(0,dt*Tsteps,Tsteps)
tR
ek = np.zeros((Tsteps, N), dtype = 'complex')

e = np.zeros((Tsteps, 1), dtype = 'complex')

ev = np.zeros((Tsteps,1), dtype = 'complex')

norm = np.zeros((Tsteps,1), dtype = 'complex')

meabs = np.zeros((Tsteps,N), dtype ='complex').T
#print(psi) Uncheck to make sure everything is proper
for t in range(Tsteps):
    ψ
##T propgators

    for n in range(N):
        temp2 = ψ[n,:]
        temp3 = temp2.reshape(My,Mz)
        temp4 = fft.fftshift(fft.fft2(temp3))
        ek[t,n] = np.real(np.sum(np.sum(np.conj(temp4)*T*temp4))/np.sum(np.sum(np.conj(temp4)*temp4)))
        temp5 = temp4*TP#Kinetic Propagator
        temp6 = fft.ifft2(fft.fftshift(temp5))
        ψ[n,:] = temp6.reshape(1,My*Mz)

####Now we check the energy conservation####
    ###ek(t)+ev(t) = constant

    meabs[:,t] = np.sum(np.conj(ψ)*ψ,1) ## population on each dibat
    e[t] =  ek[t,0]*meabs[0,t]
    norm[t] = np.real(np.sum(np.sum(np.conj(ψ)*ψ_0)))
    ψ_0 = ψ

fig = plt.figure(figsize = (5,5))
plt.plot(np.real(meabs[0,:]))
plt.xlabel('Time')
plt.ylabel(' <span class="math-container">$S_{0}$</span> Population')
plt.show()

plt.xlabel('Time')
plt.ylabel(' <span class="math-container">$S_{0_{adi}}$</span> Population')
plt.show()

plt.imshow(np.abs(ψ[0,:].reshape(Mz,Mz)))
plt.show()
plt.plot(np.abs(e))
plt.xlabel('Time')
plt.ylabel('Energy')
plt.show()

plt.plot(np.abs(norm))
plt.xlabel('Time')
plt.ylabel('Norm')

#g.s ψ

fig = plt.figure(figsize = (15,15))
X,Y = np.meshgrid(y0,z0)
ax = plt.axes(projection='3d')
surf =ax.plot_surface(X,Y,np.abs(ψ[0,:].reshape(My,Mz)),cmap = cm.inferno)
ax.view_init(25, 50)
plt.xlabel('x')
plt.ylabel('y')
fig.colorbar(surf, shrink=0.5, aspect=5)
</code>```

Everything looks good.

But when I try to extend my model to 3D, I am quite confused on how I should set up my momenta/position surfaces, and as such, how does this effect my Guassian wavepackets and kinetic propagators. Here I added the transpose to the 2D $y-component$  and to the $x-component$ as well in the $3D$ case. but I do not have a basis for doing so.
Below is my code for the 3D:

```<code>import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import scipy.special as scl
import numpy.matlib as mat
import scipy.fftpack as fft
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

##Split-Operator ###
# Constants
ω = np.array([1,1,1]) #Frequency for each coordinate
gs = np.array([0,0,0]) #Inital Wavepacket shifts
iS = 0 # Intial starting state
y0c = gs[1]

z0c = gs[2]

x0c = gs[0]

##intial Momenta

kIz = 0

kIy = 0

kIx = 0

##Set up a grid

My = 128*2

Mz = 128*2

Mx = 128*2

#Number of states

N = 1

#Number of time steps

Tsteps = 200

dt = 0.005

# Grid Lengths

Ly = 10

Lz = 10

Lx = 10

LyT = Ly*2

LzT = Lz*2

LxT = Lz*2

#Grid of M points

y0 = np.linspace(-Ly,Ly,My)

z0 = np.linspace(-Lz,Lz, Mx)

x0 = np.linspace(-Lx,Lx, Mz)

#Parameters
#k0[1xM] = Grid of M momenta points from 0->L

k0y = np.linspace(-My*np.pi/LyT,My*np.pi/LyT-2*np.pi/LyT,My)

k0z = np.linspace(-Mz*np.pi/LzT,Mz*np.pi/LzT-2*np.pi/LzT,Mz)

k0x = np.linspace(-Mx*np.pi/LxT,Mx*np.pi/LxT-2*np.pi/LxT,Mx)

##Properties

##Postion and momenta surfaces##

y0op = (np.tile(y0,(Mz,1))).T

z0op = np.tile(z0,(Mx,1))

x0op = (np.tile(x0,(My,1))).T

k0yop = (np.tile(k0y,(My,1))).T

k0zop = (np.tile(k0z,(Mx,1)))

k0xop = (np.tile(k0x,(My,1))).T

##Inital wavefunction: 2D gaussian
ψ_0 = np.zeros((N,My*My),dtype = 'complex')
σ_y = np.sqrt(2/ω[1])
σ_z = np.sqrt(2/ω[2])
σ_x = np.sqrt(2/ω[0])

temp = np.tile(np.exp(-((z0-z0c)/(σ_z))**2)*np.exp(1.j*kIz*z0),(Mx,1))

temp_y = np.tile(np.exp(-((y0-y0c)/σ_y)**2)*np.exp(1.j*kIy*y0),(Mz,1)).T

temp_1 = temp_y*temp*np.tile(np.exp(-((x0-x0c)/σ_x)**2)*np.exp(1.j*kIx*x0),(My,1)).T

ψ_0[iS,:] = temp_1.reshape(1,My*Mz)##Gaussian wavepacket

##
ψ_0[iS,:] = ψ_0[iS,:]/np.sqrt(ψ_0[iS,:]@ψ_0[iS,:].T)##Normalized wavefunction

#Now we need to calculate the propagators

##Kinetic T = exp(i*hbark^2.2mdt/hbar)

TP = np.exp(-1.j*(np.tile((k0y**2)/2,(Mz,1)).T+np.tile((k0x**2)/2,(My,1)).T+np.tile((k0z**2)/2,(Mx,1)))*dt)

T =  np.tile((k0y*k0y)/2,(Mz,1)).T +np.tile((k0x*k0x)/2,(My,1)).T + np.tile((k0z*k0z)/2,(Mx,1))

ψ = ψ_0 #Initialization of wavepacket

tR = np.linspace(0,dt*Tsteps,Tsteps)

ek = np.zeros((Tsteps, N), dtype = 'complex')

e = np.zeros((Tsteps, 1), dtype = 'complex')

ev = np.zeros((Tsteps,1), dtype = 'complex')

norm = np.zeros((Tsteps,1), dtype = 'complex')

meabs = np.zeros((Tsteps,N), dtype ='complex').T
#print(psi) Uncheck to make sure everything is proper
for t in range(Tsteps):
    ψ
##T propgators

    for n in range(N):
        temp2 = ψ[n,:]
        temp3 = temp2.reshape(My,Mz)
        temp4 = fft.fftshift(fft.fft2(temp3))
        ek[t,n] = np.real(np.sum(np.sum(np.conj(temp4)*T*temp4))/np.sum(np.sum(np.conj(temp4)*temp4)))
        temp5 = temp4*TP    #Kinetic Propagator
        temp6 = fft.ifft2(fft.fftshift(temp5))
        ψ[n,:] = temp6.reshape(1,My*Mz)

####Now we check the energy conservation####
    ###ek(t)+ev(t) = constant

    meabs[:,t] = np.sum(np.conj(ψ)*ψ,1) ## population on each dibat
    e[t] =  ek[t,0]*meabs[0,t]
    norm[t] = np.real(np.sum(np.sum(np.conj(ψ)*ψ_0)))
    ψ_0 = ψ

fig = plt.figure(figsize = (5,5))
plt.plot(np.real(meabs[0,:]))
plt.xlabel('Time')
plt.ylabel(' <span class="math-container">$S_{0}$</span> Population')
plt.show()

###Visualization###
plt.plot(np.abs(e))
plt.xlabel('Time')
plt.ylabel('Energy')
plt.show()

plt.plot(np.abs(norm))
plt.xlabel('Time')
plt.ylabel('Norm')

#g.s ψ

fig = plt.figure(figsize = (15,15))
X,Y = np.meshgrid(y0,z0)
ax = plt.axes(projection='3d')
surf =ax.plot_surface(X,Y,np.abs(ψ[0,:].reshape(My,Mz)),cmap = cm.inferno)
ax.view_init(25, 50)
plt.xlabel('x')
plt.ylabel('y')
fig.colorbar(surf, shrink=0.5, aspect=5)

</code>```

If anyone could help guide me through this problem, or explain to me the correction it would help greatly.

Thanks, and if further information is required please let me know.

Title: Split Operator 2D-->3D


问:
Solving general initial value problem $\mathbf{f}(\mathbf{x}, \mathbf{x'}, t)=\mathbf{0}$
The common form of initial value problem that can be solved using ODE integrator is
$$
\mathbf{x'}=\mathbf{g}(\mathbf{x}, t)
$$
where $\mathbf{x'}=\partial\mathbf{x}/\partial t$. The initial conditions $\mathbf{x}(t=0)$ are given. There has been a lot of prewritten solver that can be used to solve the equation above (e.g. scipy.integrate.solve_ivp) that typically uses Runge-Kutta algorithms.

However, if I change the form above into a more general one, i.e.
$$
\mathbf{f}(\mathbf{x}, \mathbf{x'}, t)=\mathbf{0},
$$
there is no algorithm I can find that can solve it (assuming it contains sufficient information to be solved). Please note that the first form above is just a specific case of the second form, i.e.
$$
\mathbf{f}(\mathbf{x}, \mathbf{x'}, t)=\mathbf{x'}-\mathbf{g}(\mathbf{x}, t)=\mathbf{0}.
$$

Is there any algorithm available out there that can solve the general form of initial value problem above?
Any reference to literatures or codes are appreciated!


问:
Numerov method for 1D potential barrier problem: numerical errors
I wanted to solve a 1D potential barrier problem. It seemed easy enough to do: set a travelling wave coming from the left, part of it gets reflected by the barrier, and part of it gets gets transmitted on the right:
$$\Psi(z) = \left\{
\begin{array}{ccc}
e^{ikz} + r e^{-ikz} & \mathrm{if} & z<=0 \\
te^{ikz} & \mathrm{if} & z>=L
\end{array}
\right.
$$
Then I create an array containing the potential $V(z)$, an array for the wave function, and I start with the boundary condition $\Psi(N\Delta z)$, $\Psi((N-1)\Delta z)$, where $(N-3)\Delta z =L$ and $N$ is the number of grid points.

Then I use Numerov's method to get from $N$ to $0$ and find the wave function.

To check, I calculate $e^{ikz}$, and I compare to the solution for $\Psi$ if $V(z)=0$ and the energy of the particle is $k^2/2m$. They should be identical, but they are slightly different. It seems that the finer grid doesn't help much.

Is there an issue I'm missing? I also tried to make the grid a lot finer than the period of the oscillation along $z$.

OwnerDisplayName: Magicsowon


问:
How to solve the minization problem with $Z^{T}Z=I$
I have the following problem
$$\underset{Z}{\min}\left\|X-AZBH\right\|_F^2+\left\|H-B^{T}Z^{T}A^{T}X\right\|\\
\text{subject to } Z^{T}Z=I$$

How to solve the variable $Z$?


问:
Interpolation question
I have a set of data $(x_i,y_i,y'_i)$, $i=1,\dots,N$ and I want to fit an interpolating curve $f(x)$ which matches both the data $y_i$ and the first derivatives $y'_i$ at the nodes $x_i$,
\begin{align}
f(x_i) &= y_i \\
\frac{d}{dx} f(x_i) &= y'_i
\end{align}
At first I thought that piecewise cubic Hermite interpolation would be well suited for this, however this method imposes no constraints on the second derivatives, and I am finding that the interpolated curve has discontinuous second derivatives at the endpoints, which is undesirable.

So I am think I should impose the additional constraint that the second derivative of $f(x)$ is continuous across the nodes $x_i$. If I consider piecewise polynomials between the nodes, it looks like I would need at least degree 4 polynomials (cubics don't have enough degrees of freedom to match the second derivative constraint). So if define a polynomial of degree 4, $p_k(x)$ on $[x_k,x_{k+1}]$, $k=1,\dots,N-1$, then these piecewise polynomials would have to satisfy these constraints:
\begin{align}
p_k(x_k) &= y_k \\
p'_k(x_k) &= y'_k \\
p_k(x_{k+1}) &= y_{k+1} \\
p'_k(x_{k+1}) &= y'_{k+1} \\
p''_k(x_{k+1}) &= p''_{k+1}(x_{k+1})
\end{align}
The above equations represent a total of $5(N-1) - 1$ constraints. There are a total of $5(N-1)$ unknown parameters if each $p_k(x)$ is a quartic polynomial. This gives me 1 additional degree of freedom, which I could use, for example, to constrain the second derivative either at the beginning or end of the interval (but not both!)

I am a little surprised that such a scheme exhibits this asymmetry, where I can only constrain either $p''_1(x_1)$ or $p''_{N-1}(x_N)$ but not both. Does anyone know of a reference for doing an interpolation like this, and whether there exists a high-quality implementation in some numerical library?

Title: Interpolation question


问:
Feasibility of synthetic, simulated predictions w/o matching against empirical physical measurement?
Feasibility of synthetic, simulated predictions w/o matching against empirical physical measurement?

This has puzzled me lately.

So, while it's sometimes common to approach modelling and model development by making empirical observations and measurements first (e.g. small experiments). And then developing a model that "matches" against these observations.

However, is it possible to also take a "synthetic" approach to prediction in the sense of:

if the phenomenon is complex

then

why not create a model that measures general patterns that do not depend on "observation-level" phenomena? I.e. why not create a model that does not require collecting empirical observations.

However, as it sounds a bit "unorthodox", then are there perhaps resources that would discuss the feasibility of this idea?

An example:

Could one predict the growth of cancer by extrapolating from properties that increase the risk of cancer growth, instead of creating a model that matches against empirically observed cancer growth patterns?


问:
Advection diffusion equation using Crank-Nicolson with total flux and Diriclet BCs
I am trying to model the 1D advection-diffusion equation:

$${\partial c \over \partial t} = D_c{\partial^2 c \over \partial x^2} -u{\partial c \over \partial x}.$$

With Robin boundary conditions that has a total flux BC at the left boundary (x=0):

$$uc-D_c{\partial c \over \partial x}=uC_f,$$

and a zero-flux BC at the right boundary (x=L):

$${\partial c \over \partial x}=0.$$

Using the Crank-Nicolson method the problem can be discretized and an example of this is presented in this paper https://publications.waset.org/10010828/analysis-of-one-dimensional-advection-diffusion-model-using-finite-difference-method. The spatial domain is discretized in points $i$ from $0-N$ while the temporal domain is discretized with $j$ from $0-M$. Crank-Nicolson in equation form reads

$$PC^{j+1}=QC^j+R.$$

in which $C^{j+1}$ is the (N-1)x(1) matrix containing the unknown values of the current timestep, $C^j$ is the (N-1)x(1) matrix containing the known values of the previous timestep. $R$ is of shape (N+1)x(1) (I assume because it includes points i=0 and i=N) and contains boundary conditions.

I understand the discretization of the non-boundary nodes, but how exactly are the BCs done?
Also, for N nodes, why are the $P$ and $Q$ matrices (N-1)x(N-1)?


问:
Conserved current for a 1d oscillator using Maple
I found out the conserved current for harmonic oscillator with angular frequency 1, particle falling in gravity close to ground ($g=1$) using maple. I'm unable to understand the result:

$$J[t](t, x(t), \frac{dx(t)}{dt}) = F1(x(t)^2 + (\frac{dx(t)}{dt})^2, -arctan(x(t)/\frac{dx(t)}{dt}) + t)$$

For particle falling under gravity:

$$J[t](t, x(t), \frac{dx(t)}{dt}) = F1(-t + \frac{dx(t)}{dt}, t^2/2 - t*\frac{dx(t)}{dt}+ x(t))$$

What do these terms mean? Also, how does one calculate conserved current for these systems?
I thought that conserved current is a 4-vector that satisfies continuity equation. I'm unable to make that connection here.


问:
Appropiate Artificial Boundary Conditions for the radial part of the Klein Gordon equation?
I am trying to simulate the following equation using FDTD

$ \left(- \partial^2_t + \partial^2_x + V(x) \right) \psi(x,t) =0 $

subjected to the initial conditions $\psi(x,0) = f(x),~ \partial_t \psi(x,0) = g(x)$, the boundary conditions being that it is completely outgoing on the right boundary, and completely ingoing at the left boundary.

Obviously, I need to implement either absorbing/artificial boundary conditions (ABCs) or something fancy like Perfectly Matching Layers. Focusing on the former, I have tried using Mur's first order ABCs and the standard non-reflecting ABCs (following Everstein) just as it is done for the good old wave equation. However, I am getting a pretty huge reflected wave. Can someone suggest how to approach the problem and minimize reflections at the artificial boundary? Or, how I may construct a better ABC?


问:
Solve simultaneous differential equations with embedded functions and a parameter estimation
The aim is to solve the below equations and plot $m$ with time, i.e. $\frac{dm}{dt}$

$k$ is unknown and needs to be estimated. For the parameter estimation, the below values in the table for m versus t can be used. This is a dummy dataset, as I would like to get an idea if this problem is solvable. If needed, an initial guesstimate for $k$ can be provided $6\times10^{-5} \mathrm{m^2/(bar \cdot hr)}$.

I have added an initial idea for the code, but not sure if my implementation of the $A(H)$ and $P(H)$ functions in the odes is correct. The code as it is doesn't run. Also I need a way to estimate for $k$ whilst fitting to the provided $m$ vs $t$ data.

Surface area in $\mathrm{m}^2$:
$$A(H)=\pi r^2+\frac{2}{r}(1+\frac{\rho_0}{\rho_H}H)\frac{w_0}{\rho_0}$$
where $A(H)_{t=0}=2\pi rL+2\pi r^2$. Note that $H$ is differentiated later in the equation $\frac{dH}{dt}$.

Pressure in bar:
$$P(H)=5004.6H+150.3$$

Volumetric flowrate in $\frac{\mathrm{m^3}}{\mathrm{hr}}$:
$$\frac{dV}{dt}=\frac{k}{h}A(H)P(H)$$

Hydration in $\mathrm{hr}^{-1}$
$$\frac{dH}{dt}=\frac{\rho_H}{w_0}\frac{k}{h}A(H)P(H)$$
where $H=\frac{w_0}{w_h}$

Mass flowrate in $\frac{\mathrm{kg}}{\mathrm{hr}}$
$$\frac{dm}{dt}=C_d\frac{dV}{dt}$$

Where

$L=0.106 \mathrm{m}$

$D=0.00264 \mathrm{m}$, $r=D/2$

$\rho_0=1320\frac{\mathrm{kg}}{\mathrm{m}^3}$ (material 1)

$\rho_h=1000\frac{\mathrm{kg}}{\mathrm{m}^3}$ (water)

$w_0(t=0)=6.48\times10^{-4}$ (decreases with time)

$w_H(t=0)=0$ (increases with time)

$L=0.0015 \mathrm{m}$

$C_d=300 \frac{\mathrm{kg}}{\mathrm{m}^3}$

Time (hr)
Release Rate (kg/hr)

2
.000012

4
.000018

6
.00002

8
.000015

10
.00002

12
.000015

14
.0000075

16
.000001

18
.000002

20
.000001

22
0

24
0

```<code>function dydt = diffunTAR(~,y,k, h, A_H, P_H, rho_h, rho_o, wo)
dydt = zeros (2,1);
D = 2.64/1000 %mm/1000 to convert to m
r=D/2
rho_o = 1320 %kg/m3
rho_h = 1000 % kg/m3
wo=648/1000000 %weight component 1 in mg converted to kg
h= 1.5/1000 %1.5mm expressed in m
        function A_H= pi*(D/2)^2 + (2/r)*(1+(rho_o/rho_h*y(2)))*(wo/rho_o);
        end
            function P_H = 5004.6*y(2) + 150.39;
            end
%dV/dt
dydt(1)= (k/h)* function A_H* function P_H;
%dH/dt
dydt(2)= (rho_h*k/wo*h)*function A_H*function P_H;
end

tspan = linspace(0,21*24); %time in hrs for 21 days
y0 = [0 6.4800e-04];
[t,y]=ode45(@(t,y) diffunTAR(t,y, k, h, A_H0, rho_h, rho_o, wo), tspan, y0);
</code>```

Edit
I realised that my previous syntax had errors, the below might be improved, but still gives me an error, displayed at the end.

```<code>function dydt = diffunTAR(~,y,k,h,rho_h,rho_o,wo,L,r)
dydt = zeros (2,1);
%dV/dt
dydt(1)= k/h* AH * PH;
%dH/dt
dydt(2)= (rho_h*k/wo*h) * AH * PH;
AH
PH
function PH
    (5004.6*y(2))+150.3;
end
function AH
    pi*(r^2)+(2/r)*(1+(rho_o/rho_h*y(2)))*(wo/rho_o);
    setInitialConditions (AH, 2*pi*(r)*L+2*pi*(r)^2);
end
end
</code>```
```<code>D = 2.64/1000 %mm/1000 to convert to m
r=D/2
L = 10.6/100 %cm/100 to convert to m
wo=648/1000000 %weight osmotic tablets in mg converted to kg
k=6e-5 %guess for permeability
h= 1.5/1000 %1.5mm expressed in m
rho_o = 1320 %kg/m3
rho_h = 1000 % kg/m3
tspan = linspace(0,21*24); %time in hrs for 21 days
y0 = [0 6.4800e-04];
[t,y]=ode45(@(t,y) diffunTAR(t,y,k,h,rho_h,rho_o,wo,L,r), tspan, y0);
</code>```

Error

Error using diffunTAR/AH
Too many output arguments.

Error in diffunTAR (line 4)
dydt(1)= k/h* AH * PH;

Error in FitODETAR>@(t,y)diffunTAR(t,y,k,h,rho_h,rho_o,wo,L,r)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);


问:
Contact analysis does not converge due to the projection falls outside valid domain
I implemented Node-To-Surface contact algorithm (Wriggers, Peter, Computational contact mechanics., Berlin: Springer (ISBN 3-540-32608-1/hbk). xii, 518 p. (2006). ZBL1104.74002.). The code is done by first forming the contact element, comprising of a slave node on the ball and a master line of the base. In the example of a ball penetrating into a soft block, the simulation could run until a point that the projection within a contact element falls outside the corresponding master in the Newton-Raphson iteration. Subsequently, the implicit NR loop does not converge anymore. Does anyone have experience dealing with this issue?

This is the snapshot of the system before failing

This is the issue, noted the black circle, the slave node (blue) located close to one of the master nodes (purple), leading to the issue above.


问:
How do we approximate the numerical error a numerical scheme (e.g Runge Kutta, Euler etc) makes without having access to an analytical solution?
So I recently encountered this question in my head while taking my Scientific Computing class, where the lecturer talked about computing numerical error of a scheme.
My guess would be that we take a much higher order scheme and treat its solution close to the analytical one and compute the error for the scheme we want but this obviously doesn't answer the question for very high order numerical schemes.</p>


问:
Solve Rational Equation for Root Music in MATLAB
I'm trying to estimate DOA in the Hybrid architecture using root music so I need to solve the attached equation to find the roots for the Root_Music equation in Matlab. Does anyone have an idea for solving this?

$$f_\theta(\theta)=g^H(\theta)\mathbf{a}^H_{M}(\theta)\mathbf{E}_{N}\mathbf{E}^H_{N}\mathbf{a}_{M}g^H(\theta)=0$$

$$g(\theta)=\frac{1-\exp^{j\frac{2\pi}{\lambda} M d\sin({\theta})}}{1-\exp^{j\frac{2\pi}{\lambda}d\sin({\theta})}}$$

$$\mathbf{a}_{M}=[1,\exp^{j\frac{2\pi}{\lambda} M d\sin({\theta})},.......,\exp^{j\frac{2\pi}{\lambda} (K-1)M d\sin({\theta})}]$$

$$f_\theta(\theta)=\frac{2-\exp^{-j\frac{2\pi}{\lambda} M d\sin({\theta})}-\exp^{j\frac{2\pi}{\lambda} M d\sin({\theta})}}{2-\exp^{-j\frac{2\pi}{\lambda}d\sin({\theta})}-\exp^{j\frac{2\pi}{\lambda}d\sin({\theta})}}\sum_{m=1}^K\sum_{n=1}^K\exp^{-j\frac{2\pi}{\lambda} (m-1)d\sin({\theta})}\mathbf{C}_{mn}\exp^{j\frac{2\pi}{\lambda}(n-1)d\sin({\theta})}=0$$

$K$ is the number of subarrays, $M$ is the number of elements in each subarray, $d$ is the antenna elements separation and we decompose covaiance matrix of the received signal into signal subspace and noise subspace $\mathbf{E}^H_{N} $

$$z=\exp^{j\frac{2\pi}{\lambda} M d\sin({\theta})}$$

$$\mathbf{C}=\mathbf{E}_{N}\mathbf{E}^H_{N}$$

$$f_z(z)=\frac{2-z^{-1}-z}{2-z^{\frac{-1}{M}}-z^{\frac{1}{M}}}\sum_{m=1}^K\sum_{n=1}^Kz^{-(m-1)}\mathbf{C}_{mn}z^{(n-1)}=0$$

so I need to get the $2K-2$ roots of this equation to estimate $\theta$

i solved the second part only using:
$$\sum_{m=1}^K\sum_{n=1}^K\exp^{-j\frac{2\pi}{\lambda} (m-1)d\sin({\theta})}\mathbf{C}_{mn}\exp^{j\frac{2\pi}{\lambda}(n-1)d\sin({\theta})}=\sum_{l=-K+1}^{K-1}\mathbf{C}_{l}z^l=0$$
where $\mathbf{C}_{l}$ is the sum of the elements along the lth diagonal of $\mathbf{C}$

is it possible to consider only the roots of the second part??


问:
Comparing block versus non-block Krylov methods for handling multiple right-hand-sides
Suppose I wish to solve a linear system $AX=B$ iteratively where $A$ is an $m\times m$ matrix and $X,B$ are $m \times s $ matrices (not single vectors). Instead of solving $s$ independent systems I'm exploring ways to solve them simultaneously so that matrix-vector products may better use caching (i.e. calculating $AX$ should require only a single read of $A$ whereas $Ax_1,\ldots,Ax_s$ would require $s$ reads of $A$.)

One approach would be simply to feed this as a linear system of size $ms$ into e.g. GMRES.

I have seen other so-called Block-GMRES algorithms which use instead a block-Arnoldi factorization.

The block-GMRES approaches appear a little more complex to implement and it's not clear if the added complexity is worthwhile.

Is it known whether block-GMRES outperforms non-block-GMRES in terms of:

memory required to achieve solution

number of matrix-matrix products $AX$ required to compute solution

Improved communication profile (maybe block-GMRES requires fewer synchronization points? or can pipeline better?)


问:
Differences between Two-Grid Correction and V-Cycle Scheme
I found in A Multigrid Tutorial: Second Edition by William L. Briggs, Van Emden Henson, Steve F. (here) the following Schemes:

Two-Grid Correction Scheme, p. 37:

And V-Cycle Scheme, p. 40

I want to check that I understood it correctly, Two-Grid is just a multigrid scheme when the algorithm going one level coarser. It's useful for cases where it possible to solve $A^{2h}e^{2h}=r^{2h}$ equation, otherwise it implemented recursively which basically creating the V-Cycle. Did I understand it right? Is there any additional advantages and disadvantages of Two-Grid Vs. V-Cycle?


问:
How to position oneself, if one has little practical engineering background or function, but may have a broader insight into mathematical modeling?
How to position oneself, if one has little practical engineering background or function, but may have a broader insight into mathematical modeling?

Or, lets say:

I could write models to prove properties about the thing perhaps, but I wouldn't necessarily build that
thing.

And this is not the same as "writing some Matlab alongside work to see what I'm doing".

So how does a "software person" position oneself so that it's also viable towards "practical", w/o having the experience of a "practical engineer"?

This is also perhaps why the degree contents for mech. engr. or elec. engr. or software is different, while they're a bit of overlap. But the profiles are different.

So what's it that a "computational methods expert" looks for instead of what a mech. engr. using Matlab or similar looks for?

What I think:

Such computational methods expert should be concerned about implementing software that aids engineering analysis. W/o necessarily being the one doing such analysis. But studying the design of such software systems.


问:
Find tuples of points from multiple sets
Given n sets of points in general position in dimension 2 (n typically small, 2-6), can one find tuples of points, one from each of the sets, which are close in some sense (the closest, mutual distances less than some threshold)?  I'd be interested even in the usual name for this problem.


问:
Decomposing a banded matrix
Suppose we have a linear algebra problem with a banded matrix A which has nonzero entries on the main diagonal, two nearest sub-diagonals, and two other sub-diagonals (such band structure often arises in the context of numerical solution of PDEs in 2D),

$
A x = b
$

and we write A as a sum of a 3-diagonal matrix $T_1$ (nonzero entries on the main diagonal and two nearest sub-diagonals), and a 3-banded matrix $T_2$ (nonzero entries on the main diagonal and two sub-diagonals not nearest to the main one),

$
A = T_1 + T_2
$

to attempt solving the linear system iteratively, using an efficient solver for 3-diagonal matrices,

$
T_1 x^{k+1} = b - T_2 x^k
$

One can symmetrize it by alternating $T_1$ and $T_2$ on the left-hand side.

When would this strategy be successful? If this method works, there must be a standard name for it - what is it? Is this a variation of the ADI method? What is the optimal way to decompose $A$ into $T_1$ and $T_2$?

Title: Decomposing a banded matrix


问:
How to troubleshoot numerical instability using finite difference for steady-state non-linear heat conduction equation
I have a problem which I believe is numerical instability when trying to solve a heat conduction equation using finite difference. The short version is that when the parameter $I=80.3$ I get the blue curve below (positive values and with one local maximum, which is physically reasonable), but when $I=80.4$ I get the orange curve which is not physically correct. I would appreciate any help troubleshooting this problem.

Here is the long description, which I divided into physical modelling (for completeness's sake), finite difference discretization and numerical implementation (where I believe lies the problem):

1: Physical Modelling.

The problem I am trying to solve is of a thin straight metallic wire, whose both ends are held at constant temperatures $T_0$ and $T_L$. Through the wire flows a current $I$ which heats up the wire through Joule heating. The electrical resistivity is temperature-dependent as $\rho(T) = \rho_0 + \rho_1 \cdot T$

2: Finite Difference Discretization.

I divide the wire in $Nx$ nodes, for each inner node I write the following balance equation

$$q_{i-1} + q_{i+1} + q_e = 0 $$

Where $q_{i-1}$ is the heat coming from node $i-1$, $q_{i+1}$ is the heat coming from node $i+1$ and $q_e$ the heat generated by Joule heating.

I apply Fourier law $q = -k \cdot A \cdot \Delta T / \Delta x$, where $k$ is heat conductivity, $A$ is area, $\Delta T$ temperature difference and $\Delta x$ is the length between two nodes.

$$ \frac{k \cdot A}{\Delta x} \cdot (T_{i-1} - T_i) + \frac{k \cdot A}{\Delta x} \cdot (T_{i+1} - T_i) + (\rho_0 + \rho_1 \cdot T_i) \cdot \frac{\Delta x}{A} \cdot I^2 = 0 $$

Reordering I get

$$ \alpha \cdot T_{i-1} + (-2\alpha + \rho_1 \cdot \frac{\Delta x}{A} \cdot I^2) \cdot T_i + \alpha \cdot T_{i+1} = -\rho_0 \cdot \frac{\Delta x}{A} \cdot I^2  $$

3: Numerical Implementation.

I put the left side of the previous equation in a sparse tridiagonal matrix $A$, the right side in a vector $b$ and use scipy to solve the system of equations using float64. Here is the Python code:

```<code># Imports
import numpy as np
import scipy.sparse.linalg as la
import scipy.sparse as sp
import matplotlib.pyplot as plt

#% Geometric specifications

D_L = 0.4E-3 # in m
L_L = 12E-3  # in m
T_0 = 87  # in °C
T_L = 175 # in °C

# Material constants

rho0 = 1.6E-8 # Ohm m
rho1 = 6.6E-11 # Ohm m / K
k  = 3.94E2 # W/m K

Nx = 30 # Number of intervals

x   = np.linspace(0, L_L, Nx+1)
d_x = x[1] - x[0]

for I in np.array( [80.3, 80.4] ):

    u = np.zeros(Nx+1)

    # Representation of sparse matrix and right-hand side
    main  = np.zeros(Nx+1)
    lower = np.zeros(Nx)
    upper = np.zeros(Nx)
    b     = np.ones(Nx+1)

    # Construction of matrix equations

    A_c   = 3.14*(D_L)**2 /4
    alpha = k * A_c / d_x
    R     = rho0 * d_x / A_c
    Pe    = R * I**2
    b     = -Pe * b

    # Precompute sparse matrix
    main[:]  = -2*alpha + rho1 * d_x / A_c * I**2
    lower[:] = alpha
    upper[:] = alpha

    # boundary conditions
    main[0]  = 1; upper[0] = 0; b[0]  = T_0
    main[-1] = 1; lower[-1]= 0; b[-1] = T_L

    A = sp.diags( diagonals=[main, lower, upper], offsets=[0, -1, 1],
              shape=(Nx+1, Nx+1), format='csc' )

    u[:]   = la.spsolve(A, b)

    plt.plot( x, u, label='I='+str(I) )

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

Does somebody know how to troubleshooting this problem? Many thanks in advance!


问:
Slow convergence of Stokes solver used with the Immersed Boundary method
I am using Immersed Boundary Method to simulate elastic particles in 3D Stokes flow. Specifically, one has $\nabla ^2 \mathbf{u}-\nabla p + \mathbf{f}(t) = 0$, $\nabla \cdot \mathbf{u} \; $, where $\mathbf{f} (t)$ is calculated based on instantaneous configuration of some particle assembly that is advected by flow $\mathbf{u}$. More explicitly, $ \dot{\mathbf{r}}_i = \mathbf{u}(\mathbf{r}_i) \;$, $\mathbf{f}(t)=\mathcal{F}(\{\mathbf{r}_i\})$, $\mathcal{F}$ denoting some function(al) that depends on all $\mathbf{r}_i$:s.

Within the Immersed Boundary Method, $\mathbf{u}$ is discretized on a spatial (non-moving) grid. To simulate a time step, one first solves a discretized linear system to obtain $\mathbf{u}$ for a given instantaneous value of $\mathbf{f}$, then moves particles at their positions $\mathbf{r}_i$:s with the velocity of the flow at those positions: $\mathbf{r}_i (t+\delta t) = \mathbf{r}_i (t) + \delta t \mathbf{u}(\mathbf{r}_i,t) $. Finally, one updates $\mathbf{f}$ according to the newly obtained values of $\mathbf{r}_i$:s.

Broadly, I can see two ways of dealing with Stokes dynamics:

(1) The method described above, requiring to solve a large linear system in every step.

(2) Simulating transient flow by adding $\rho \partial _t \mathbf{u}$ term and artificial compressibility to the equation.

In my hands, using method (1) is slow since it requires solving a linear system in every time-step. I use iterative MINRES solver implemented in PETSc, but the number of iterations to reach convergence tends to be large (some 4000).

Method (2) does not require solving a linear system, but introduces a small parameter which degrades accuracy unless the parameter is made really small (requiring small time step).

My question is twofold:

(1) Using method (1), I do recycle my solution from a previous time-step as an initial guess for the iterative solver. Yet, typically, the number of iterations required seems large (~4000). Am I doing something wrong?

(2) If using method (2), is there some "clever" way to correct for non-zero inertia and obtain something closer to the limit of zero Reynolds? Perhaps something based on a combination of both strategies?


问:
How to write a simple finite element solver in python in order to solve Poisson equation in 2D
I would like to write a simple finite element solver in python in order to solve 2D Poisson equation and then visualize it.

$$
-\nabla^{2} u(x,y)=f(x,y), \quad x,y \quad in \quad \Omega\\
u(x,y) = u_D (x,y), \quad x,y \quad in \quad \partial \Omega
$$

Why do I wanna do that? I would like to learn how write a simple finite element solver using python and since the Poisson equation is like the "hello world!" for the computational science I would like to work with it.
I already had look at FEniCS, deal.II and some other solvers but as I mentioned I would like to write my own simple script.

Would you introduce any simple to understand source where I can learn how to solve and code this equation using finite element method?


问:
Numerical evaluation of Duhamel's integration
I am trying to numerically evaluate the following Duhamel's integration:

$$
x = \frac{-1}{\omega_d} \int_0^t \ddot{x}_g (\tau) e^{-\zeta \omega_n(t - \tau)} \sin{\left( \omega_d (t - \tau) \right)} d\tau
$$
where $\omega_d$, $\zeta$ and $\omega_n$ are constants.

From our data, $\ddot{x}_g(t_i)$ are given at each discrete time increment $t_i$, but I am struggle to compute $x(t_i)$, or the Duhamel's integral, correctly?

In my MATLAB script, I have tried to store $e^{-\zeta \omega_n(t_i)} \sin{\left( \omega_d (t_i) \right)} $ to an array H, and then the numerical value of $\ddot{x}(t_i)$ for $i = 1, \cdots N$ is stored inside another array xdd (same dimension as H). I then tried to use the convolution function inside MATLAB, i.e.:

```<code>X = conv(xdd, H)
</code>```

However, it does not look like my attempts are working out. My questions, along with what I have tried, are:

From the documentation, the resulting length for X is length(xdd) + length(H) - 1, but shouldn't the dimension of X be the same with xdd (used to store $\ddot{x}$)? I would expect that at each value of $\ddot{x}$, there would be a corresponding value of $x$. Therefore, X and xdd should have the same dimension. But that doesn't seem to be the case here...

I have also tried the 'same' for the Shape options from the MATLAB documentation above, but that also doesn't seem to give me the right answers. The resulting X that I got is at a much higher frequency that what we would expect analytically.

Is this the correct strategy to arrive at the numerical solution for $x$?


问:
Overflow in Solving numerically the 1D Kuramoto-Sivashinsky Equation using spectral methods
I'm having trouble with the following implementation of the KS model (see below) found on

Solving numerically the 1D Kuramoto-Sivashinsky Equation using spectral methods

When tN > 300 an overflow occurs on the line:

u_hat2[:,j+1] = (1 / nx) * np.fft.fftshift(np.fft.fft(u[:,j+1]**2))

I tried using a workaround with the convolution theorem with no use, any help would be much apreciated.

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

from matplotlib.pyplot import cm

nu = 1
L = 100
nx = 1024

t0 = 0
tN = 200
dt = 0.05
nt = int((tN - t0) / 0.05)

# wave number mesh
k = np.arange(-nx/2, nx/2, 1)

t = np.linspace(start=t0, stop=tN, num=nt)
x = np.linspace(start=0, stop=L, num=nx)

# solution mesh in real space
u = np.ones((nx, nt))
# solution mesh in Fourier space
u_hat = np.ones((nx, nt), dtype=complex)

u_hat2 = np.ones((nx, nt), dtype=complex)

# initial condition
u0 = np.cos((2 * np.pi * x) / L) + 0.1 * np.cos((4 * np.pi * x) / L)

# Fourier transform of initial condition
u0_hat = (1 / nx) * np.fft.fftshift(np.fft.fft(u0))

u0_hat2 = (1 / nx) * np.fft.fftshift(np.fft.fft(u0**2))

# set initial condition in real and Fourier mesh
u[:,0] = u0
u_hat[:,0] = u0_hat

u_hat2[:,0] = u0_hat2

# Fourier Transform of the linear operator
FL = (((2 * np.pi) / L) * k) ** 2 - nu * (((2 * np.pi) / L) * k) ** 4
# Fourier Transform of the non-linear operator
FN = - (1 / 2) * ((1j) * ((2 * np.pi) / L) * k)

# resolve EDP in Fourier space
for j in range(0,nt-1):
  uhat_current = u_hat[:,j]
  uhat_current2 = u_hat2[:,j]
  if j == 0:
    uhat_last = u_hat[:,0]
    uhat_last2 = u_hat2[:,0]
  else:
    uhat_last = u_hat[:,j-1]
    uhat_last2 = u_hat2[:,j-1]

  # compute solution in Fourier space through a finite difference method
  # Cranck-Nicholson + Adam
  u_hat[:,j+1] = (1 / (1 - (dt / 2) * FL)) * ( (1 + (dt / 2) * FL) * uhat_current + ( ((3 / 2) * FN) * (uhat_current2) - ((1 / 2) * FN) * (uhat_last2) ) * dt )
  # go back in real space
  u[:,j+1] = np.real(nx * np.fft.ifft(np.fft.ifftshift(u_hat[:,j+1])))
  u_hat2[:,j+1] = (1 / nx) * np.fft.fftshift(np.fft.fft(u[:,j+1]**2))

# plot the result
fig, ax = plt.subplots(figsize=(10,8))

xx, tt = np.meshgrid(x, t)
levels = np.arange(-3, 3, 0.01)
cs = ax.contourf(xx, tt, u.T, cmap=cm.jet)
fig.colorbar(cs)

ax.set_xlabel("x")
ax.set_ylabel("t")
ax.set_title(f"Kuramoto-Sivashinsky: L = {L}, nu = {nu}")<span class="math-container">```
```


问:
Hessian-free preconditioner for non linear least squares
I am solving a nonlinear least squares problem using Gauss Newton method. Due to the large dimension of the problem, I use the Hessian-free approach. As a linear solver I use either MINRES or CG.  To speed up convergence, I began to figure out which preconditioners can be easily added.
First, I decided to try the Fisher preconditioner that Martens talks about in his work:
$$
M = \left[ diag \left( \sum_{i=1}^{D} \nabla f_{i}(\theta)\odot \nabla f_{i}(\theta)\right) + \lambda I )\right]^a
$$

In the article I use $\alpha$ = 0.75, but in my case it only worked at 0.25. The optimization error began to decrease faster, but the linear solution error increased relative to identity.
Then I came across an interesting work Preconditioning for Hessian-Free Optimization by Robert Seidl.

In it, I was interested in the SPAI Preconditioner, where, using multiplication, we pull out the elements with the greatest gradients.
To begin with, I decided to pull out the real diagonal elements corresponding to the greatest magnitude of the gradients, and use them to calculate the preconditioner. Leave the remaining elements of the diagonal equal to 1.

If I take the entire deagonal, the convergence is good, as soon as I reduce the number of elements by at least 2%, it becomes worse than the identity preconditioner. And this is strange, I assumed that any addition of information to the preconditioner should not make it worse than the identity preconditioner. A similar situation if I make SPAI preconditioner like Robert Seidl. Any advice on this?


问:
Trouble inverting complex matrix with numpy and scipy
I have some matrix-valued, complex data $Z(f)$ with $f\in\{f_0,f_1,\dots\}$ and $Z(f_i)$ being a 3x3 matrix. I require the inverse $Z^{-1}(f)$ in my workflow. After encountering some problems with my results, I cold narrow the issue down to the inversion, which seems to be performed not very accurately by numpy.linalg.inv and to a lesser extend also by scipy.linalg.inv.

As you can see in the plot of the off-diagonal element $Z_{12}$ below, at intermediate values $10^3<f <10^5$ applying the numpy inversion twice does not yield the original data $Z$. In the plot, the scipy inversion seems to have no issues, but using scipy to compute $Z^{-1}Z$ also does not result exactly in the unity matrix at these values of $f$.

While the condition number of $Z(f)$ diverges for $f\rightarrow 0$ as seen in the second plot below, I would like to understand why my issues can occur at intermediate $f$ values, and how to choose an inversion algorithm that minimizes such problems.


问:
switching boundary conditions to simulate a reciprocating flow in COMSOL
This is a re-posted question which was poorly defined earlier

I have been trying to simulate reciprocating flow of a convective fluid through a heated channel in COMSOL-Multiphysics. The schematic of the domain is as follows:

There is reciprocating flow of a convective fluid through the channel through the highlighted boundaries L and R, while it is being heated from the bottom $(q'')$. I will now explain the problem in more detail. A sinusoidal mass flux at some tempertature $T_i$ (only the positive values in the mass flux curve, otherwise zero where the sinusoid turns negative) has been defined at L. Similar mass flux profile with the same frequency, amplitude but with a phase lag of $\pi$ is applied at R . In this scenario, for a cycle of time-period $T$, the flow should reciprocate between L and R.

For the half cycle $0<t<T/2$, the boundary L should act as an inlet, while R should act as the oulet. During the next half of the cycle, i.e., $T/2<t<T$, the fluid enters from R and leaves through L. Thus R must now acts an inlet and L as the outlet.

It must also be noted that the simulation also models heat transfer. So for $0<t<T/2$, L should be acting as if it admits a fluid at temperature $T_f$ and R should act as the outlet (ususally modeled as zero conductive flux). The phenomenon will be opposite for the next half, i.e., $T/2<t<T$.

I want to achieve this switching of boundary condition in each half cycle inside a transient simulation in COMSOL, for some finite number of cycles until the system reaches a quasi-steady state.


问:
Rosenthal equation for multi track
Rosenthal's equation lets one calculate the temperature profile of a moving point heat source analytically for thin and thick plates. For simplicity I use the equation for thick plates defined as:
$$T-T_0 = \frac{Q}{2\pi kR}\exp\bigg(-\frac{v(R+x)}{2a}\bigg)$$
Where $T$ is the temperature one wants to know, $T_0$ is the initial temperature, $Q$ is the heat source heat flux, $R$ is the radius from heat source center, $v$ is the heat source velocity, $k$ the thermal conductivity, $c$ the thermal capacity, $\rho$ the density, $a=\frac{k}{\rho c}$ and $x$ is the current position of the heat source in traveling direction.

The equation was initially developed to calculate the temperature for single track welding applications. I am wondering how it would be possible to use this equation for multi track welding (where one track is located right next to the last track) and how to capture the effect of the previous track in the new track. Can I simply update the material properties in case of temperature dependent properties and use an updated $T_0$ for the new track?

Title: Rosenthal equation for multi track


问:
How amenable is this 2D Frenkel–Kontorova-like energy minimization problem in Python to the use of a modest PC + GPU? (Heavy reliance on indexing)
@Richard's answer to Going to try to move some of my scipy/numpy calculation to a new GPU, how to avoid disappointing results? is quite helpful, and as promised I've added a simple running example below of the kind of project I've started.

While there may be "high level libraries" out there that can do this, I always code/script my projects myself first (to the extent feasible) in order to get insight/appreciation into the problem.

I have a layer of atoms on a surface and am trying to do a 2D Frenkel–Kontorova-like energy minimization using at first simplex then once that works, with simulated annealing. The number of atoms could range from 100 to 10,000.

I'll be looking at deviations from coincident lattices (see also related code golf).

Each evaluation of energy is calculated from the atomic positions with respect to the substrate lattice and nearest-neighbor bond distances (and later, angles as well) within the layer, and will be done with numpy. In the future I may need to call scipy.special functions also.

The script below first builds the problem and stores it in the object bob for convenience (plotting, interrogation, rendering with Blender, etc.) but then extracts a few constants that are stored in args and three arrays:

atom_positions an N x 3 array of positions of atoms, order arbitrary except that the first element is the centermost of the group. These are the values that a minimization routine will adjust.

atoms_1 and atoms_2 are numpy indices for the atom_position array that are used to identify all pairs of atoms for which bond energies should be calculated; they are in this case nearest neighbors and to avoid double counting there are three per atom on the interior and fewer on the edges/corners. The order is again arbitrary.

There is an alternative way to formulate the numpy position array for minimization by making it a masked array. It might be slightly larger to add an extra "dummy atom" all the way around, and bond energies would be calculated by simple shifts, something like x_masked[1:, :] - x_masked[:-1, :] letting the predefined masks deal with atoms on the edge with fewer bonds.

note 1: Later I will be introducing defects in the lattice and including honeycomb structures (e.g. Xenes). In this case indexing has the advantage that you can introduce these quickly by editing the index arrays. Trying to do the defects in a regular array may become impossible.

note 2: The example script calls a generic Scipy minimization routine for convenience. I will ultimately be scripting a simulated annealing algorithm and experimenting to find a good annealing schedule, but that can't really happen until I get this running faster, likely by GPU augmentation of a PC.

Question: How amenable is this 2D Frenkel–Kontorova-like energy minimization problem in Python to the use of a modest PC + GPU? Is heavy reliance on indexing going to bite me later as I scale it up to 10,000 atoms?

Tiny example (runs in seconds on a laptop) with nmax=6.  With nmax=20 for example you get 1261 atoms and a much more challenging problem to calculate.

```<code>bob = Problem('bob', sig=0.2, k_bond=20, a_sub=a_substrate, a=a_lattice, h=height,
          R=R, origin=origin, nmax=6)
</code>```

note 3: per request I have tried to do some profiling using scaline but the results seem trivial, to minimize energy most time is spent in the scipy minimizer and the line of numpy in my function mini_me that calculates the energy to minimize. Results are here https://pastebin.com/GfneKJhB and in nicely formatted html here https://pastebin.com/ZJAL8gqv

note 4: colors of the atoms and the bonds reflect their relative energies.

```<code>import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib import colors as mcolors
from scipy.optimize import minimize
import time

class Problem():
    def __init__(self, name=None, nmax=12, a_sub=1., k_vert=1., a=1.,
                 h=1.0, k_bond=2., R=0., origin=np.zeros(2), sig=0., seed=42):
        self.name = str(name)
        self.nmax = int(nmax)
        self.a_sub = float(a_sub)
        twopi = 2 * np.pi
        r3o2 = np.sqrt(3) / 2
        g = (twopi / (self.a_sub * r3o2)) * np.array([[0, 1], [r3o2, 0.5], [-r3o2, 0.5]]).T
        self.g_sub = g[::-1] # for purposes of calculating hexagonal substrate potential
        self.k_vert = float(k_vert)
        self.a = float(a)
        self.h = float(h)
        self.k_bond = float(k_bond)
        self.R = float(R)
        self.origin = np.array(origin, dtype=float)
        self.sig = float(sig)
        self.seed = int(seed)
        self.initial_hex_lattice(origin=origin)
        np.random.seed(self.seed)
        self.random = self.sig * np.random.random(self.x0_init.shape)
        self.x_initial = self.x0_init + self.random
        self.x = self.x_initial.copy() # this array will be updated
        # self.x_final = None
        self.n_atoms = self.x_initial.shape[0]
        self.define_nearest_neighbor_bonds()
        self.calculate_bond_energies()
        self.calculate_surface_energies()

    def initial_hex_lattice(self, origin):
        """generates an N x 3 array of atom locations, origin atom is first, wight
           Gaussian noise in 3D. The reduction to 1/6 slice by hex-hex symmetry is optional"""
        i, j = np.mgrid[-self.nmax:self.nmax+1, -self.nmax:self.nmax+1]
        keep = (np.abs(i + j) <= self.nmax) * ~((i == 0) * (j == 0)) # add origin back later
        i, j = [thing[keep] for thing in (i, j)]
        r3o2 = np.sqrt(3) / 2.
        x = self.a * (i + 0.5 * j)
        y = self.a * r3o2 * j
        s, c = [f(np.radians(self.R)) for f in (np.sin, np.cos)]
        x, y = c * x - s * y, c * y + s * x
        z = np.zeros_like(x)
        # origin = np.array([0, 0, self.h])
        self.x0_init = origin + np.vstack([np.zeros(3), np.stack([x, y, z], axis=1)])
        self.i = np.hstack(([0], i))
        self.j = np.hstack(([0], j))
        self.ij = np.stack((self.i, self.j), axis=1)

    def define_nearest_neighbor_bonds(self):
        """calculate index pairs for all nearest-neighbor bonds in the hexagonal
           lattice"""
        a_bonds, b_bonds, c_bonds = [], [], []
        n_atom = np.arange(self.n_atoms)
        A, B, C, D = (self.i < self.nmax, self.i + self.j < self.nmax,
                      self.j < self.nmax, self.i > -self.nmax)
        self.A, self.B, self.C, self.D = A, B, C, D
        self.atoms_1a = n_atom[A * B]
        self.atoms_1b = n_atom[B * C]
        self.atoms_1c = n_atom[C * D]
        self.atoms_2aij = self.ij[self.atoms_1a] + np.array([1, 0]) # to the right
        self.atoms_2bij = self.ij[self.atoms_1b] + np.array([0, 1]) # to the above-rith
        self.atoms_2cij = self.ij[self.atoms_1c] + np.array([-1, 1]) # to the above-left
        self.atoms_2a = np.array([np.where(np.all(self.ij == atom , axis=1))[0][0]
                             for atom in self.atoms_2aij])
        self.atoms_2b = np.array([np.where(np.all(self.ij == atom , axis=1))[0][0]
                             for atom in self.atoms_2bij])
        self.atoms_2c = np.array([np.where(np.all(self.ij == atom , axis=1))[0][0]
                             for atom in self.atoms_2cij])
        self.atoms_1 = np.hstack((self.atoms_1a, self.atoms_1b, self.atoms_1c))
        self.atoms_2 = np.hstack((self.atoms_2a, self.atoms_2b, self.atoms_2c))

    def calculate_bond_energies(self, which='x'):
        """harmonic oscilator potential for bond lengths, minimum at 'a'
           with spring contsant 'k_bond'"""
        array = getattr(self, which)
        vectors = array[self.atoms_2] - array[self.atoms_1]
        self.bond_lengths = np.sqrt((vectors**2).sum(axis=1))
        self.bond_energies = self.k_bond * (self.bond_lengths - self.a)**2

    def calculate_surface_energies(self, which='x'):
        """ applies smooth sinusoidal-like hexagonal pattern with lattice constant 'a'
            to a  N x 2  grid of x-y points. Minimum at origin."""
        array = getattr(self, which)
        E = np.cos((array[:, :2, None] * self.g_sub).sum(axis=1))
        result = (1.5 + E.sum(axis=1)) / 4.5
        surface_energy_xy = 1. - (1.5 + E.sum(axis=-1)) / 4.5 # "sinusolidal" in x, y
        surface_energy_z = self.k_vert * (self.x[:, 2] - self.h)**2 # harmonic osc. in z
        self.surface_energies = surface_energy_xy + surface_energy_z

    def get_energy(self, which='x'):
        self.calculate_bond_energies(which)
        self.calculate_surface_energies(which)
        return self.bond_energies.sum() + self.surface_energies.sum()

    def potential_yx(self, yx):
        """ applies smooth sinusoidal-like hexagonal pattern with lattice constant 'a'
            to a 2 x M x N y-x grid of points. Minimum at origin."""
        E = np.cos((yx[:, None, ...] * self.g_sub[::-1][..., None, None]).sum(axis=0)) # sum over x, y then cos()
        result = 1. - (1.5 + E.sum(axis=0)) / 4.5 # then sum over three directions
        return result

    def plot_it(self, hw_pixels=250, cmap_sub='gray', cmap_atoms='inferno',
                cmap_bonds='jet', linewidth=2.5, vmin=0., vmax=1.2, which='x'):
        array = getattr(self, which)
        hw = 1.05 * np.abs(array[:, :2]).max() # largest cartesian distance of an atom from center + 10%
        self.plot_scale = hw / hw_pixels  # 0.05 angstroms per pixel
        self.yx = self.plot_scale * np.stack(np.mgrid[-hw_pixels:hw_pixels+1,
                                                 -hw_pixels:hw_pixels+1], axis=0)
        self.substrate = self.potential_yx(yx=self.yx)
        self.plot_extent = [self.yx[1].min(), self.yx[1].max(),
                            self.yx[0].min(), self.yx[0].max()]
        print('self.substrate.shape: ', self.substrate.shape)
        # energy histograms
        a, b = np.histogram(self.surface_energies, bins=(self.n_atoms >> 2))
        c, d = np.histogram(self.bond_energies, bins=(self.n_atoms >> 2))
        # color coded bond lines
        my_cmap = plt.get_cmap(cmap_bonds)
        colors = my_cmap(self.bond_energies / self.bond_energies.max()) # normalized for cmap
        linez = np.stack((self.x[bob.atoms_1][:, :2], self.x[bob.atoms_2][:, :2]), axis=1)
        ln_coll = LineCollection(linez, colors=colors, linewidths=linewidth)

        fig = plt.figure(constrained_layout=False, figsize=[9, 7.5])
        gs = fig.add_gridspec(nrows=4, ncols=5, left=0.08, right=0.92,
                              bottom=0.07, top=0.96, hspace=0.15)
        ax1 = fig.add_subplot(gs[:3, :])

        # fig, ax = plt.subplots(1, 1, figsize=[9, 7.5])
        ax1.imshow(self.substrate, origin='lower', extent=self.plot_extent,
                  cmap=cmap_sub, vmin=vmin, vmax=vmax) # surface in background
        ax1.add_collection(ln_coll) # then bonds
        x, y = self.x.T[:2]
        ax1.scatter(x, y, c=self.surface_energies, cmap=cmap_atoms, s=60,
                   edgecolors='k', zorder=2) # then atoms on top of the bonds
        ax1.set_ylabel('Angstroms')
        ax2 = fig.add_subplot(gs[3:, 1:4])
        ax2.plot(b[1:], a, label='surface energies')
        ax2.plot(d[1:], c, label='bond energies')
        ax2.legend()
        ax2.set_yscale('log')
        ax2.set_xlabel('energy')
        plt.show()

def ang(i, j, k, l):
    root3 = np.sqrt(3)
    th_ij = np.arctan2(root3 * j, 2*i+j)
    th_kl = np.arctan2(root3 * l, 2*k+l)
    return np.degrees(th_ij - th_kl)

#### NEXT
# print old/new mean origin
# print old/new mean rotation
# print old/new mean bond length

# is there any collective behavior?

# BEGIN

a_lattice = 3.5
a_substrate = 2.9
height = a_lattice
R = 9.

origin = [0, 0, height]

bob = Problem('bob', sig=0.2, k_bond=20, a_sub=a_substrate, a=a_lattice, h=height,
              R=R, origin=origin, nmax=6)

bob.plot_it()

def mini_me(atom_positions_flattened, args):
    g_sub, k_vert, a, h, k_bond, atoms_1, atoms_2 = args
    atom_positions = atom_positions_flattened.reshape(-1, 3)
    # bond energies
    vectors = atom_positions[atoms_2] - atom_positions[atoms_1]
    bond_lengths = np.sqrt((vectors**2).sum(axis=1))
    bond_energy = k_bond * ((bond_lengths - a)**2).sum()
    # surface energies
    E = np.cos((atom_positions[:, :2, None] * g_sub).sum(axis=1))
    result = (1.5 + E.sum(axis=1)) / 4.5
    surface_energies_xy = 1. - (1.5 + E.sum(axis=-1)) / 4.5 # "sinusolidal-like" in x, y
    surface_energies_z = k_vert * (atom_positions[:, 2] - h)**2 # harmonic osc. in z
    surface_energy = (surface_energies_xy + surface_energies_z).sum()
    return bond_energy + surface_energy

atom_positions = bob.x.copy()

parameters = ('g_sub', 'k_vert' , 'a', 'h', 'k_bond', 'atoms_1', 'atoms_2')

args = [getattr(bob, attr) for attr in parameters]

if True:
    print('start minimization!')
    print('mini_me(atom_positions, args): ', mini_me(atom_positions, args))
    maxiter = 1000
    options = {'maxiter': maxiter}
    t_start = time.process_time()
    wow = minimize(mini_me, atom_positions.flatten(), args=args,
                   tol=1e-2, options=options) # method='Nelder-Mead', bounds=bounds
    process_time = time.process_time() - t_start
    print('wow.fun: ', wow.fun)
    print('wow.message: ', wow.message)
    print('process time: ', process_time)
    print('number of iterations: ', wow.nit)
    print('iterartions per second: ', wow.nit / process_time)

bob.x_final = wow.x.reshape(-1, 3) # insert final results back into object

print('bob.get_energy(which="x_initial"): ', round(bob.get_energy('x_initial'), 3))
print('bob.get_energy(which="x_final"): ', round(bob.get_energy(which='x_final'), 3))

bob.plot_it(which='x_final')
</code>```


问:
Numerical calculation of out-of-time order correlators (OTOCs)
I want to numerically calculate OTOCs for different quantum mechanical systems. Consider the following Hamiltonian
$$H=p_x^2+p_y^2+x^2y^2$$

and I want to calculate the following OTOC
$$C_T(t)=-\left<[x(t),p_x(0)]^2\right>_\beta$$
where
$$\left<\mathcal{O}\right>_\beta=Z^{-1}\text{Tr}(\mathcal{O} e^{-\beta H}),\qquad Z=\text{Tr}(e^{-\beta H})$$

You can look into the details for the recipe for calculating the OTOC in this paper: "Out-of-time-order correlators in quantum
mechanics", in section 2.

My question is that I can only manage to calculate OTOC's using Mathematica, but it is very slow with numerical calculations. The essential idea behind the calculation of OTOCs is to find the energy eigenfunctions and energy eigenvalues for the given Hamiltonian. Then using the OTOC recipe, we can easily calculate the OTOC. Therefore I am looking for a programming language that can most efficiently calculate the eigenfunctions and energy eigenvalues of a given Hamiltonian, then calculate the matrix elements of the position operator $x_{nm}=\left<n\right|x\left|m\right>$ using parrallel computation. Parrallel computation is key since the matrix elements of the position operator form a quite large matrix if we want accurate results.

Do you know what programming language the following papers use

Out-of-time-order correlators in quantum mechanics

Lyapunov Exponent and Out-of-Time-Ordered Correlator's Growth Rate in a Chaotic System

Does scrambling equal chaos?

Gauging classical and quantum integrability through out-of-time ordered correlators


问:
How to assign initial velocity field and handle pressure-velocity coupling in FVM?
I am trying to solve the 2D incompressible Navier-Stokes equations for laminar flow over a backward facing step using the finite volume method.

This is the plot that I generated of a generic mesh showing u1, u2, and p nodes to help myself derive the discrete integrals in the conservation equations:

The grid is staggered to avoid pressure checkerboard-ing. I am doing explicit time integration and I am using a uniform grid. Except near the boundaries, I apply central differencing scheme for derivatives, and I perform second order linear interpolations to obtain values at the control volume faces.

I am following the procedure from Ferziger and Peric in this picture:

Currently my solution is blowing up, and I have a few reasons for why I suspect this to be the case, but I have some general questions that may help me narrow down why this is happening.

When I assign an inflow boundary condition of u1, how do I ensure that the initial velocity field is divergence-free? For example, if I assign a u1 inflow condition and the rest of the u1 nodes are 0, it seems like the velocity field will of course not be divergence free.

I am handling the pressure velocity coupling by taking the divergence of the momentum equation which leads to the equation shown in the picture below from Ferziger and Peric. I then apply Gauss's theorem to both sides of the equation to turn this into a surface integral equation that I can apply the finite volume method to.

Is this typically how the velocity-pressure coupling is handled for explicit time integration?

In all my reading and searching, I have not seen it done this way. I have read about pressure-correction methods like SIMPLE, but based on the procedure in the image above, it seems that I do not need to apply a pressure correction method.

Does it mean that if the velocity field at time n is divergence free, I can just solve the pressure Poisson equation in integral form and then use this pressure field to obtain a divergence free velocity field at time n+1? Or is this not an adequate way to handle the pressure-velocity coupling?

Any help would be greatly appreciated.


问:
Symmetric Matrix in Eigen C++
I am aware of a symmetric matrix type in uBLAS as ublas::symmetric_matrix matrix. Is there an equivalent for this in Eigen library that can be used to construct one or do we need to explicitly check if a matrix is symmetric or not which could be a time expensive solution (O(n^2)). Any guidance is appreciated.</p>

Title: Symmetric Matrix in Eigen C++


问:
Efficient way of calculating a cumulative integral with prefactor
I have a grid of points $x_i$ and corresponding function values $y_i=f(x_i)$. I'm interesting in something like the cumulant of $f$, but it has an awkward prefactor. The desired quantity we'll call

$$z(x) = e^{x}\int_{x}^{\infty}e^{-q}f(q)dq$$

And this quantity is needed for each grid point, $z_i=z(x_i)$

Let us assume for numerical sake that $\infty$ is really a sufficiently large $q$ such that $f(q)$ goes smoothly to zero, and that $f(q)$ is a nicely behaved, non-negative function for any $q$ we might care about. In this particular context, its a probability distribution. Note that I don't always have a functional form for $f$, so an analytical solution isn't particularly helpful here.

Were it not for the $e^x$ prefactor, integrating this would be rather simple, numerically. Either use scipy.integrate.cumulative_trapezoid or do the trapezoidal rule yourself by taking

integrand = np.exp(-x) * y
almostZ = dx * (np.cumsum(integrand[::-1])[::-1] - integrand/2 - integrand[-1]/2)
Z = almostZ * np.exp(x)
```
<p>where <code>dx</code> is the grid spacing for the <span class="math-container">$x_i$</span>. Done this way, you get <span class="math-container">$z_i$</span> in roughly <span class="math-container">$O(N)$</span> scaling, which is great.</p>
<p>The problem arises when <span class="math-container">$x$</span> gets sufficiently large that either <span class="math-container">$e^x$</span> or <span class="math-container">$e^{-x}$</span> are up against the limits of numerical precision, even though the <em>product</em> of the two terms is a reasonably sized number. Doing this in a loop is the current state of the "answer":</p>
<pre class="lang-py prettyprint-override"><code>for i in range(len(x)):
    Z[i] = np.trapz(np.exp(x[i] - x[i:]) * y[i:], x[i:])
</code>```

This ensures that there's no really huge numbers multiplied by really small numbers, and is instead a sum of reasonably sized numbers. But it also involves a lot of double counting and $O(N^2)$-ness, which is a bit of a pain when this is something that's being calculated as integrating an ODE (via solve_ivp).

Is there a trick to doing this sort of computation, or am I just doomed to $O(N^2)$?


问:
Quadrature rules for non-linear finite element problems
For solving linear problems stemming from PDEs with the FEM, such as the Poisson equation or the wave equation, it is customary to use the "simplest" numerical quadrature that exactly integrates the quantities that define the mass/stiffness matrices. More precisely, given shape functions $\phi_k$ for nodes $k$ on an element $K$, we use a quadrature rule that exactly integrates the expression

$$
\int_K \phi_i \, \phi_j \, \mathrm{d} x
$$

for each combination of $i$ and $j$ for the mass matrix, and
$$
\int_K \nabla \phi_i \, \nabla \phi \, \mathrm{d} x
$$
for the stiffness matrix. It turns out that for other linear problems, such as linear elasticity, we can use the same quadrature rule because the components of the matrix are made up of constants times products of the basis function gradients.

However, for many non-linear problems, such as non-linear elasticity with a Neo-Hookean material model, or other models that involve transcendental functions, we cannot exactly evaluate the integrals with standard quadrature rules.

My understanding is that the rule-of-thumb in these cases is to use the same quadrature as you would use for linear problems; presumably the quadrature error due to inexact integration is small compared to other sources of error. Another criterion that seems very important is to make sure that we do not lose rank of the matrix. In other words, the rank of our approximate matrix should be the same as that of the matrix obtained through exact integration.

I am looking for literature that discusses the choice of quadrature in a non-linear but smooth standard FEM setting (ignoring discontinuous material parameters or non-standard discretizations). In particular, I have several questions whose answers are not directly obvious to me:

Can we pose "minimum" requirements on a quadrature rule so that we do not lose rank of our matrix compared to exact integration? Is using the same quadrature as for the linear problem sufficient?

Are there prominent examples where we need higher accuracy than the "linear problem quadrature"?

Are there other concerns than the ones I listed here? (*)

Any advice, pointers or insights would be most welcome!

(*) In some special settings, it may be important to capture certain characteristics. For example, with a Neo-Hookean elastic model you may want to have enough quadrature points to reasonably ensure that the determinant of your deformation gradient remains positive throughout the element, in order to prevent local inversion of the element.


问:
How do I identify negative group speeds?
This question is a continuation of one of my other questions.
I've been trying to show that collocated (non-staggered) grids can suffer from negative group speeds in the linearized shallow water system. I do this by solving the equations for a simple, right-traveling wave crest, so that any left-traveling component can be identified as having negative group-speed. For my initial conditions I use a Gaussian wave profile with initial velocities $u_0(x) = h_0(x)\sqrt{g/h_0(x)}$.

My wave indeed moves in the right direction, but I have trouble with the interpretation. To the left of the crest, a trough appears with a small bump in the middle. This trough moves to the left. I see the same thing on a staggered grid. I can't spot a clear difference between the staggered and the non-staggered grid and I can't see negative group speed.

Is this trough normal? How can I identify the negative group speeds?

EDIT:
A comment suggests that I split the wave up into modes and calculate the speed of the modes separately. I'm not really familiar with such decompositions in practice, so I used the discrete fourier transform from this post. I use $2\pi$freqs as the wavenumbers, but I'm not sure if that's correct? If I then use this function on a snapshot of the wave and plug my "wavenumbers" into the formula for the group speed (according to the discretization scheme), I get the following:

I like this because the group speeds on the staggered grid are non-negative. Does this look right? Also, since I'm using the formula for the group speed, does that not make this exercise a circular reasoning? Because even when the numerical solution is exactly the same on both grids, my computation of the group speeds would still get me a different result on both grids.


问:
How to constrain the every optimized vector component to be nonnegative?
I am building a gradient descent model based on portfolio optimization. Currently, I have finished the model and am able to run it smoothly without any problem. However, there's one issue that I couldn't fix: the optimal vector components have negative values. As my portfolio does not support short-selling stocks, this optimal solution may not be a convincing solution to my portfolio.

Model Description

My cost function for the portfolio problem is defined as

$$F(x)=\frac{\beta_1}{2}(x^T\boldsymbol{\Sigma}x)-\mu^Tx+\frac{\gamma}{2}(e^Tx-1)^2+\lambda(e^Tx-1)+\frac{\beta_2}{2}\Vert{x}\Vert^2_2+\rho\sum_{i=1}^n\max(0,-x_i)^2$$

where $x$ and $\boldsymbol{\Sigma}$ are my minimizer vector and positive definite covariance matrix. $\frac{\gamma}{2}(e^Tx-1)^2+\lambda(e^Tx-1)$ is the augmented lagrangian method for the constraint $e^Tx=1$ and $\rho\sum_{i=1}^n\max(0,-x_i)^2$ is the penalty imposing $x_i>0]$ for all components of $x$.

As can be seen here, my objective is to find the optimum solution for $x$ using the Steepest Descent method. Hence, I start to build my Python code as below:

```<code>import numpy as np
from numpy.linalg import norm

def penalty_f(v):
    return  np.sum([np.power(max(0., -xi),2) for xi in v])

def penalty_df(v):
    return -2.*np.array([max(0, -xi) for xi in v])

def aug_lag(lam, v):
    new_lam = lam + gam*(np.sum(v)-1.)
    return new_lam

def lipschitz(v):
    lips = beta_1*np.sqrt(np.trace(cov@cov)) + beta_2*np.sqrt(len(v)) + gam*len(v)
    return min(1,1/lips)

def project_f(v, lam_k):
    func = 0.5*beta_1*v.T@cov@v - b*mean.T@v + gam/2*(np.sum(v)-1.)**2 + lam_k*(np.sum(v)-1.) + beta_2/2*v@v + rho*penalty_f(v)
    return func

def project_df(v, lam_k):
    par_func = beta_1*cov@v - b*mean + gam*np.ones_like(v)*(np.sum(v)-1.) + lam_k*np.ones_like(v) + beta_2*v + rho*penalty_df(v)
    return par_func

def gradient_descent(f, df, ini_v, tolerance, lam, MAX_ITER=10000, output_fname='output.txt'):
    vec = ini_v
    negative = 0
    aug_lam = lam
    with open(output_fname, 'w') as out:
        for i in range(MAX_ITER):
            alpha = lipschitz(vec)
            f_value = f(vec, aug_lam)
            gradient = df(vec, aug_lam)
            direction = np.negative(gradient)
            vec += alpha * direction
            aug_lam = aug_lag(aug_lam, vec)
            print("No. of zeros: ", len(vec)-np.count_nonzero(vec))
            print(norm(gradient,2), "," , f_value)
            msg = f'{f_value},{norm(gradient,2)}, {str(list(vec))}\n'
            out.write(msg)

            # stopping criteria
            if norm(gradient,2) < tolerance:
                print ("The optimum vector for", {df}, " is at ", vec,"at iteration ", i+1)
                for i in vec:
                    if i < 0:
                        negative += 1
                print("No of negative: ", negative)
                print("Vector sum: ",np.sum(vec))
                print("Gradient: ", norm(gradient,2))
                break

            if i == MAX_ITER:
                print ('Higher no. of iterations is needed for', {df})
                print ("Vector: ", vec)
                print("Vector sum: ",np.sum(vec))
                print("Gradient: ", norm(gradient,2))

    return vec
</code>```

Here is the result I obtained using 33 stocks:

```<code>The optimum vector for {<function project_df at 0x000001FA50FF2D30>}  is at  [-0.03442667 -0.01839447  0.03305939 -0.04319195  0.00353058  0.13003235
  0.09405175  0.00022958  0.08010633 -0.02831581 -0.04066535 -0.02875569
  0.04953858 -0.04019486  0.22866438 -0.03474848 -0.00029308  0.09792178
  0.11669689 -0.08569448  0.00405406  0.05544077  0.12273198  0.09615974
 -0.04374827  0.16994986  0.00670517  0.04772238  0.01729612  0.05218545
  0.03484764 -0.03438397 -0.00811171] at iteration  1230
No of negative:  13
Vector sum:  1.0000000001114249
Gradient:  9.961948573544883e-07
</code>```

When I assess 30 or more stocks, negative components start to appear. Is there anything I can improve or take note of when performing portfolio optimization?

P.S. The data to the mean and covariance is in the link https://www.dropbox.com/scl/fi/rycj948t4bnq5u60m13ow/Covariance.xlsx?dl=0&rlkey=d8u18ntuxk7wjcl1eup8gmxoa

```<code>import pandas as pd

df1 = pd.read_excel('Covariance.xlsx', sheet_name=0, header=None)
df2 = pd.read_excel('Covariance.xlsx', sheet_name=1, header=None)
np_cov_1 = df1.values
mean_1 = df2.values.reshape(len(df2))
ini_vec_1 = np.array([1. / (len(df2)) for i in range(len(df1))])
</code>```


问:
Sufficient condition for real roots of a polynomial of order $n>5$ with arbitrary real coefficients
I ask for help in solving the problem. I am developing an optimization program that selects the coefficients of a polynomial of order $n> 5$ so that all its zeros are just real numbers. And I cannot find sufficient conditions suitable for this, which would include its coefficients.

In the article A Sufficient Condition for All the Roots of a Polynomial To Be Real  conditions are given that seem to indicate whether the zeros of the polynomial are valid. But there are also counterexamples when the condition is satisfied, but the polynomial has zeros with a nonzero imaginary part. A modification of the method is intended only for polynomials with positive coefficients. I am working with polynomials that have arbitrary real coefficients.

Therefore, I need sufficient conditions for the validity of zeros for such a case. Who can suggest anything?


问:
Cellular automata - can scaling factors that influence probability of spread be applied to rate of spread?
This is a follow-up question to my earlier question Stochastic cellular automata - algorithm limited by 1 cell per timestep.  I am considering blending two approaches to a cellular automata model of how mold spreads in a petri dish.

A probability based approach

A quick review: I am trying to model the spread of mold in a petri dish. The petri dish can be thought of as a grid of 1mm x 1mm squares, each called a Cell. Each square has a value for the kind of food that is present there, and the probability $\rho$ of mold spreading into that Cell is based on some numerical value that can be derived from the food that is present.

We seed the dish by saying that a certain Cell has gotten some mold on it. We then iterate over a series of timesteps. Let's say each timestep is 1 minute. Within each timestep, we iterate over every Cell that has mold on it.  And for every one of those Cells, we iterate over the neighbors of the Cell. Based on the fuel for each neighbor, we decide whether or not it is likely in that step for the mold to spread to that neighbor Cell.

For simplicity, we can imagine this problem in 1 dimension. We seed the top left Cell to be considered moldy, and it can only spread in the x direction. A visual would look like this:

The effective rate of spread (RoS) emerges from $\rho$.  So if $\rho$ is 0.5, after 60 timesteps, we should have 30 moldy cells, creating a RoS of about 30mm/hr.

A Rate of spread based approach

For reasons explained in my prior question, I am switching to a RoS based approach.  Rather than relying on $\rho$, I have data available based on the type of food in each cell which gives me the RoS directly.  With that resource, I can switch to a priority-queue based model, which has been working well.

The complication - more influences

There are more factors than just food source which affect the rate of spread.  Temperature is a good example.  Some research has been done on the effect of temperature on $\rho$.  Let's say, for example, there is a temperature coefficient that goes like:

$$ \phi_T(T) = a * \frac{1}{b*e^{-(T - 75)^{3}} + 0.5}$$

Where a and b are constants.   (Google graph for the curious).  And its effect on $\rho$ would look like:

$$ \rho_{with T}(\rho, \phi_T) = 1 - (1 - \rho) ^ {\phi_T} $$

(For the curious, here's a graph of $\rho_{with_T}$ as a function of $\rho$ for $\phi_T$ = 0.5 )

The relationship between $\rho$ and RoS is linear.  The relationship between $\rho$ and $\rho_{with_T}$ is very nonlinear.

The question

Is it possible to say that the effect of T on RoS is the same as it is on $\rho$?  For example, based on the above equation, can we say

$$ RoS_{with T}(RoS, \phi_T) = 1 - (1 - RoS) ^ {\phi_T} $$

Or is this an extrapolation that we cannot make?  My instincts tell me we can make this assumption, but I can't think of why or why not.

Bonus question

If the answer to the above question is yes, does it still hold up when we expand the problem back to 2 dimensions?


问:
Fitting data with a Voigt function
I have some data, (xrd data), that I would like peak fit with a pseudo-Voigt function, a combination of a Gaussian and a Lorentzian function. These are the functions

$G(x) = I \exp\left( -\frac{4\ln(2)(x-\Theta)^2}{w_G}\right)$

$L(x) = I \frac{w_{L}^2}{w_{L}^2+\left(x-\Theta\right)^2}$

Where $I$ is the intensity, $\Theta$ is the peak position, $w_G$ and $w_L$ are the full width at half maximum for the Gaussian and the Lorentzian respectively, together they form the Voigt function

$V(x) = \eta L(x) + (1-\eta)G(x)$

where $\eta$ is a mixing parameter

The unknowns are $I$, $\Theta$, $w_G$, $w_L$, and $\eta$. The data is given by $x$.

I want to fit the data with the method of least squares, MSE, but since the problem is nonlinear, I find it to be quite a pain.
for would i have read I need to solve the Hessian for only the five unknown parameters, but I don't quite understand the math of what this solves.

So the question is, what is the math to solve this problem and if able to, some pseudo-code?

Title: Fitting data with a Voigt function


问:
How to determine the minimum grid length | Numerical Plasma physics
I am trying to understand the attached python 3X code.

The following dispersion relation is given

$$\epsilon (k, \omega) = 1 - \frac 1 2 \left[ \frac{\omega_p^2}{(\omega-kv_0)^2} + \frac{\omega_p^2}{(\omega+kv_0)^2}\right]$$

What I am trying to understand is how to find the minimum grid length $L_{min}$ (as a function of $\frac{v_0}{w}$). $L_{min}$ indicates the needed minimum grid length to support such unstable modes.

To do so the code below should be checked in order to understand the plasma behavior either side of $L_{min}$ threshold, ie $L_{min} < L$ and $L < L_{min}$. To do so we should: 1) Adjust the number of simulation particles to grid points to improve the statistics, 2) Fix the number of particles per cell (i.e. npart/ngrid) and keep it $>> 1$, in order to reduce numerical noise and the runtime needed to observe the instability can be estimated from the maximum growth rate.

The theory is OK but I am struggling to actually find the minimum grid length $L_{min}$ once given the code. Might you please shed some light and guide me through? Thank you.

```<code>#! /usr/bin/python
#
#  Python script for computing and plotting single charged particle
#  trajectories in prescribed electric and magnetic fields.
#  Roughly equivalent to boris.m matlab program

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.widgets import Slider, Button, RadioButtons
from mpl_toolkits.mplot3d import Axes3D
import os
import os.path
import sys
from sys import exit
from time import sleep

# ===================================
#
# Function to integrate particle trajectory
# in given E, B fields
#
# ===================================

def integrate(E0, B0, vz0):
   global dt, v0, x0, xp, yp, zp, qom, larmor, nsteps
   wc=qom*B0 # cyclotron frequency
   larmor=vperp/wc
   print ("Cyclotron frequency =",wc)
   print ("Perpendicular velocity v_p=",vperp)
   print ("Larmor radius=",larmor)

   norm = 1.  # choose whether to normalise plot axes dimensions to Larmor radius
   trun=5*2*np.pi/wc  # total runtime
   dt=.1/wc  # timestep - adjust to current B-field

   nsteps=int(trun/dt)  # timesteps
   E=np.array([0.,E0,0.])  # initial E-field
   B=np.array([0.,0.,B0])  # initial B-field
   u=np.array([0.,0.,0.])  # intermediate velocity
   h=np.array([0.,0.,0.])  # normalized B-field
   xp[0]=x0[0]
   yp[0]=x0[1]
   zp[0]=x0[2]
   v0[2]=vz0 # z-component

   v=v0+.5*dt*qom*(E+np.cross(v0,B)) # shift initial velocity back 1/2 step
   x=x0

   for itime in range(1,nsteps):
     x=x+dt*v
     xp[itime]=x[0] /norm
     yp[itime]=x[1] /norm
     zp[itime]=x[2] /norm
     tp[itime]=itime*dt
#
# Boris mover: solves dv/dt = q/m*(E + vxB) to 2nd order accuracy in dt
#
     qomdt2 = dt*qom/2
     h = qomdt2*B
     s=2*h/(1+np.dot(h,h))
     u = v + qomdt2*E
     up=u+np.cross(u+np.cross(u,h),s)
     v=up+qomdt2*E

#     vxp[itime] = v[0]

# ===================================

# Make 2D plots of particle orbit
#
# ===================================

def plot_track2D():
  global xp,yp,nsteps,ax1

  fig = plt.figure(figsize=(8,8)) # initialize plot
  xmin=np.min(xp)
  xmax=np.max(xp)
  ymin=np.min(yp)
  ymax=np.max(yp)
  fig.add_subplot(221) # 1st subplot in 2x2 arrangement
  plt.cla()
  plt.grid(True, which='both')
  plt.xlim( (xmin, xmax) )
  plt.ylim( (ymin, ymax) )
  plt.xlabel('<span class="math-container">$x$</span>')
  plt.ylabel('<span class="math-container">$y$</span>')
  plt.plot(xp[0:nsteps],yp[0:nsteps],c='b')

  fig.add_subplot(222) # 2nd subplot

#  fig.add_subplot(223) # 2nd subplot
#  fig.add_subplot(224) # 2nd subplot

  plt.draw()
  plt.savefig('./particle_orbit.png') # Save plot to file

# ===================================
#
# Make 3D plot of particle orbit
#
# ===================================

def plot_track3D():
  global xp,yp,zp,nsteps,ax1
  xmin=np.min(xp)
  xmax=np.max(xp)
  ymin=np.min(yp)
  ymax=np.max(yp)
  zmin=np.min(zp)
  zmax=np.max(zp)
  ax1.cla()

  plt.ion()
  plt.grid(True, which='both')
  ax1.set_xlim( (xmin, xmax) )
  ax1.set_ylim( (ymin, ymax) )
  ax1.set_zlim( (zmin, zmax) )
  ax1.set_xlabel('<span class="math-container">$x $</span> [m]')
  ax1.set_ylabel('<span class="math-container">$y $</span> [m]')
  ax1.set_zlabel('<span class="math-container">$z $</span> [m]')
#ax1.set_aspect(1.)
  ax1.scatter(xp,yp,zp,c=tp,marker='o') # tracks coloured by elapsed time since start
  plt.draw()

# =============================================
#
#  Main program
#
# =============================================

print ("Charged particle orbit solver")
plotboxsize   = 8.
animated = True

x0=np.array([0.,0.,0.])     # initial coords
vz0=0.
v0=np.array([-1e2,0.,vz0]) # initial velocity
vperp = np.sqrt(v0[0]**2+v0[2]**2)
E0=0.
B0=.1

e=1.602176e-19 # electron charge
m=9.109e-31 # electron mass
qom=e/m  # charge/mass ratio

wc=qom*B0 # cyclotron frequency
larmor=vperp/wc
print (wc,vperp,larmor)

trun=5*2*np.pi/wc  # total runtime
dt=.1/wc  # timestep - adjust to current B-field

nsteps=int(trun/dt)  # timesteps
B1=np.array([0.,0.,0.1])  # gradient B perturbation

#wc=qom*np.linalg.norm(B) # cyclotron frequency

#nsteps=2
tp = np.zeros(nsteps)  # variables to store particle tracks
xp = np.zeros(nsteps)
yp = np.zeros(nsteps)
zp = np.zeros(nsteps)
vxp = np.zeros(nsteps)
vyp = np.zeros(nsteps)
vzp = np.zeros(nsteps)

# Compute orbit
integrate(E0, B0, vz0)

# 2D orbit plotter
plot_track2D()

exit(0) # Quit script before 3D plot - comment out to continue!

# Start 3D interactive mode with sliders for B, E and v0

plt.ion() # Turn on interactive plot display
fig = plt.figure(figsize=(8,8))
# Get instance of Axis3D
ax1 = fig.add_subplot(111, projection='3d')

# Get current rotation angle
print (ax1.azim)

# Set initial view to x-y plane
ax1.view_init(elev=90,azim=0)
ax1.set_xlabel('<span class="math-container">$x $</span>[microns]')
ax1.set_ylabel('<span class="math-container">$y $</span>[microns]')
ax1.set_zlabel('<span class="math-container">$z $</span>[microns]')
plot_track3D()

#filename = 'a0_45/parts_p0000.%0*d'%(6, ts)
#plot_from_file(filename):
axcolor = 'lightgoldenrodyellow'
axe0 = fig.add_axes([0.1, 0.95, 0.3, 0.03])#, facecolor=axcolor) # box position, color & size
axb0  = fig.add_axes([0.5, 0.95, 0.3, 0.03])#, facecolor=axcolor)
axv0  = fig.add_axes([0.1, 0.9, 0.3, 0.03])#, facecolor=axcolor)

sefield = Slider(axe0, 'Ey [V/m]', -5.0,5.0, valinit=E0)
sbfield = Slider(axb0, 'Bz [T]', -1.0, 1.0, valinit=B0)
svz = Slider(axv0, 'vz [m/s]', 0.0, 1.0, valinit=0.)

def update(val):
    E0 = sefield.val
    B0 = sbfield.val
    vz0 = svz.val

    integrate(E0,B0,vz0)
    plot_track3D()
    plt.draw()

sefield.on_changed(update)
sbfield.on_changed(update)
svz.on_changed(update)

resetax = fig.add_axes([0.8, 0.025, 0.1, 0.04])
button = Button(resetax, 'Reset', color=axcolor, hovercolor='0.975')
def reset(event):
    global ax1
    sefield.reset()
    sbfield.reset()
    svz.reset()
    ax1.cla()
    ax1.set_xlabel('<span class="math-container">$x $</span>[microns]')
    ax1.set_ylabel('<span class="math-container">$y $</span>[microns]')
    ax1.set_xlim( (0., 10.) )
#    ax1.set_ylim( (-sigma, sigma) )
    ax1.grid(True, which='both')
    plt.draw()
button.on_clicked(reset)

#plt.show()
plt.show(block=False)
</code>```


问:
Solution fails to converge with different collocation point selection
I'm trying to learn about collocation. As an example, I am using Numerical Recipes, 20.7.12, but changing the basis. To wit, I'm trying to solve

$$y'' + y' - 2y + 2 = 0, -1 \le x \le 1, y(-1) = y(1) = 0$$

which has exact solution $y(x) = 1 - \frac{e^{x} \sinh(2) + e^{-2x}\sinh(1)}{\sinh(3)}$.
Now, Numerical Recipes is using the Chebyshev basis, but just to generate some understanding of the method I tried to use the following basis:

$$y(x) = \sum_{n=1}^{N} a_n \sin\left(\frac{n\pi}{2}(x+1)\right)$$

The collocation method becomes

$$ \sum_{n=1}^{N} \left[\left(2 + \frac{n^2\pi^2}{4} \right) \sin\left(\frac{n\pi}{2}(x_m+1) \right) - \frac{n\pi}{2} \cos\left(\frac{n\pi}{2}(x_m+1) \right)\right]a_n = 2$$

where we choose $N$ points for $x_m$.

My understanding was that I had lots of freedom to choose $x_m$. So for example, I could choose $x_m = -1.0 + 2(m-1)/(N-1)$, or $x_m \sim \mathrm{unif}(-1,1)$, or $x_m = \cos(\frac{2m-1)\pi}{2N})$ (the Chebyshev nodes). The only thing to worry about was convergence rate or conditioning of the resulting matrix.

Is this belief correct? In the following code, uniformly spaced nodes converge (albeit slowly) but choosing the Chebyshev nodes for the collocation points appear to diverge!

```<code>using LinearAlgebra
using Plots
Plots.GRBackend()

function exact(x::Float64)
    return 1 - (exp(x)*sinh(2.0) + exp(-2*x)*sinh(1.0))/sinh(3.0)
end

function approximate(x::Float64, v::Vector{Float64})
    y = 0.0
    for n in 1:length(v)
        k = n*π/2
        y += v[n]*sin(k*(x+1))
    end
    return y
end

function ode(N::Int64)
    if N < 2
      error("N >= 2 is required.")
    end
    M = Matrix{Float64}(undef, N, N)
    for m in 1:N
      #=Uncomment and the solution will diverge!=#
      #=arg = cos((2*m-1)*π/(2.0*N))=#
      arg = -1.0 + 2.0*(m-1)/(N-1)
      for n in 1:N
        M[m,n] = (2+n*n*π*π/4)*sin(n*π*(arg+1)/2) - n*π*cos(n*π*(arg+1)/2)/2
      end
    end

    c = cond(M,1)
    println("Condition number of matrix is $c")
    B = lu(M)
    b  = Vector{Float64}(undef, N)
    for n in 1:N
      b[n] = 2.0
    end
    v = B\b
    p = plot([exact, x -> approximate(x, v)], xlims=(-1.0, 1.0))
    savefig(p, "solution.png")
end
</code>```

Now, if you run this code you'll see that the condition number is worse for the Chebyshev nodes at a given $N$ than for the uniformly spaced nodes. However, I'm not sure that's a sufficient explanation, because even before the condition number renders to solve worthless, then collocation fails to converge.


问:
Is a 2021 MacBook Air sufficient for modeling and simulations research in fluid dynamics?
I code mostly in Matlab, simulating various simple fluid dynamics models.

One particular simulation has the potential to become sort of complex, though I believe the set of coupled equations will still be pretty simple / greatly simplified.

Is a 2021 MacBook Air sufficient for simulations?

My 10-year-old MacBook Pro runs Matlab code perfectly fine, but it's so old now that soon I won't be able to use new versions of Matlab, etc.

My concern about the MacBook Air is that it doesn't have a fan, so I wonder if its processor chip will be damaged during high-demand simulations.

My gut feeling is to stick with what has worked well, and pay up for a new MacBook Pro.


问:
Why are fast Givens rotations mentioned so little in the recent literature?
Fast Givens rotations seem like a nice upgrade from standard Givens rotations. They take fewer multiplications to apply and avoid the calculation of a square root.

They are, however, given very little focus in the recent literature. They appear as an exercise in Golub and Van Loan. They get a short section in Gentle but, as far as I can tell, they are unmentioned entirely in Davis. Since Davis is very focused on efficient calculation this is surprising.

Is there some issue with fast Givens rotations? Do they not give performance advantages on modern computer architectures?

References:

Davis, T., 2006. Direct methods for sparse linear systems. Philadelphia: Society for Industrial and Applied Mathematics.

Gentle, J., 2007. Matrix Algebra. New York: Springer.

Golub, G. and Van Loan, C., 2013. Matrix computations. 4th ed. Baltimore: Johns Hopkins Univ Press.


问:
Poisson equation fast solver on large grids
Below is kind of a broad question.

Was given Poisson equation with BC, on meshgrid of NxM.

Rewrote as system of N*M linear equations,

Converted to sparse format (CSR/CRS/Yale), solved using Conj Grad method

As for 100x100 grid computation time is ~0.3s and
if I increase the grid to 500x500, code runs for 11 seconds already.

Was initially looking for something fast that can be used as a subroutine for Navier-Stokes solver, i.e. will run the Poisson solver a big number of times.

Question - what is the optimal computation time and how to chose the grid size properly. Or, perhaps, shall I change the method completely, or maybe something wrong with it?

Below posted my conj function for CSR input format:

```<code>double *cj(double *val, int *col, int *row, int valSize, int colSize, int rowSize, double *b, int ansSize){
    double *x;
    x = (double *) malloc(ansSize * sizeof(double));
    double *r;
    r = (double *) malloc(ansSize * sizeof(double));
    double *p;
    p = (double *) malloc(ansSize * sizeof(double));
    double *Ap;
    Ap = (double *) malloc(ansSize * sizeof(double));
    double rsold, alpha, rsnew, sum;
    int colStart = 0;
    int i, j, l, k, rowStart, rowEnd, colEnd;
    for (i = 0; i < ansSize; i++){// fill x with 0
        x[i] = 0;
    }
    for (i = 0; i < ansSize; i++){
        rowStart = row[i];
        rowEnd = row[i+1];
        sum = 0;
        for (j = rowStart; j < rowEnd; j++){
            k = col[j];
            sum = sum + val[j]*x[k];
        }
        r[i] = b[i] - sum;
    }
    for (i = 0; i < ansSize; i++){
        p[i] = r[i];
    }
    rsold = 0;
    for (i = 0; i < ansSize; i++){
        rsold = rsold + r[i] * r[i];
    }
    for (i = 0; i < ansSize; i++){
        for (l = 0; l < ansSize; l++){
            rowStart = row[l];
            rowEnd = row[l+1];
            sum = 0;
            for (j = rowStart; j < rowEnd; j++){
                k = col[j];
                sum = sum + val[j]*p[k];
            }
            Ap[l] = sum;
        }
        sum = 0;
        for (j = 0; j < ansSize; j++){
            sum = sum + p[j] * Ap[j];
        }
        alpha = rsold / sum;
        for (j = 0; j < ansSize; j++){
            x[j] = x[j] + alpha * p[j];
        }
        for (j = 0; j < ansSize; j++){
            r[j] = r[j] - alpha * Ap[j];
        }
        rsnew = 0;
        for (j = 0; j < ansSize; j++){
            rsnew = rsnew + r[j] * r[j];
        }
        if (sqrt(rsnew) < epsilon){
            return x;
        }
        for (j = 0; j < ansSize; j++){
            p[j] = r[j] + (rsnew/rsold)  * p[j];
        }
        rsold = rsnew;
    }
    free(r);
    free(p);
    free(Ap);
    return NULL;
}
</code>```


问:
robust numerical calculation with large (almost) offsetting terms
I am attempting to evaluate an expression that defines a probability of a particular event. With 3 different events, I can (under some statistical assumptions) characterize the probability of event 1 (for example) occurring as

$$\frac{e^{v_1+2 v_2} (v_2-v_1)}{\left(e^{v_2}-e^{v_1}\right)^2 \left(e^{v_2}-e^{v_3}\right)}+\frac{e^{v_1+2 v_3} (v_3-v_1)}{\left(e^{v_3}-e^{v_1}\right)^2 \left(e^{v_3}-e^{v_2}\right)}+\frac{e^{2 v_1}}{\left(e^{v_1}-e^{v_2}\right) \left(e^{v_1}-e^{v_3}\right)},$$

where $(v_1, v_2, v_3)$ are arbitrary real numbers that are associated with the different events. (I have general expressions for an arbitrary number of events, but I believe this case is sufficient to understand the problem.)

The practical numerical issue I have been struggling with is that the expression is not numerically robust when two or more elements in $(v_1, v_2, v_3)$ are numerically similar in value, as the denominator in the expression involves pairwise differences. In the case that they are identical I can calculate limits. However, when they similar the associated numerical error could result in an evaluation that is outside of the unit interval as at least two or more of the fraction terms in the expression will become very large in absolute terms (and theoretically almost offsetting).

Any help in how I should approach this would be much appreciated. Thank you.