问:
Quantifying Superiority of Golden Section vs Binary Search vs Brent on Continuous Function
This concerns algorithms for searching for the minimum of a continuous unidimensional unimodal function.  The simplest method is three points, and bi-sectioning.  The golden-section search seems more natural.  Although Brent seems a little better yet.

I am trying to find a reference as to how much more efficient these searches are relative to one another.

I need something like an "expected number of function samples to reach tolerance xtol."  One can assume a uniform distribution of the minimum, and some distribution of the brackets, perhaps.

does someone have a good reference?


问:
Best Possible Convex bounds for optimization problems
Suppose we have a primal problem

$$
p^{*}=\min_x f(x), \\\text{s.t.}~~ h_i(x) \leq 0,
$$

where $f(.)$ and $h_i(.)$ are possibly non-convex.

Then its Lagrangian is

$$\mathcal{L}(x,z_i)= f(x) + \sum_i z_i h_i(x)$$

and the dual problem is

$$d^{*}=\max_{z_i\geq 0} \min_x \mathcal{L}(x,z_i)$$

Now it is well-known that $p^{*}\geq d^{*}$ and that the dual problem is convex. Is this the best possible convex lower bound in general? Are there instances of well-known problems where we know that there exists a better lower bound which also comes from a convex optimization problem.

ContentLicense: CC BY-SA 4.0


问:
How to prevent BFGS from getting stuck on astronomically large gradient?
I have implemented BFGS myself from scratch in order to solve minimization problems.  Part of BFGS, as I understand it, is that the approximation to the Hessian is supposed to be positive definite, which means it is supposed to be invertible.   I have been following the algorithm outlined here: https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm#Algorithm .  This is how we update our approximation to the Hessian:

$$ B_{k+1} = B_k + \frac{ y_k y_k ^T }{y_k^T s_k} – \frac{ B_k s_k s_k^T B_k}{s_k^T B_k s_k} $$

where
$ y_k = \nabla f(x_k) - \nabla f(x_{k-1}) $ and

$s_k = -1 * B_k^{-1} \nabla f (x_k)$

I made a few slight modifications to the standard BFGS algorithm in my implementation:

1)  I check to see if $y_k^T s_k$ is negative.  If it is, then this has the potential to make the next update $B_k$ negative, too.  So, I switch to an alternative update when that term is negative:

$$ B_{k+1} = B_k + \frac{ y_k y_k ^T }{y_k^T y_k} – \frac{ B_k s_k s_k^T B_k}{s_k^T B_k s_k} $$
2)  I have a check for when the direction it is about to step $d_k$ might be too large or too small compared to the magnitude of the gradient:  if either $||d_k|| \ge 100*||\nabla f|| $ or $||\nabla f^T * d_k|| \le 10^{-4} * ||\nabla f||*||d_k||$.  If either one of those are true, then I just use the gradient step, $d_k = \nabla f $.

Now I am running into an issue where the gradient somehow becomes very large:  $\nabla f$ becomes very large.  So, when I go to update $s_k$, it becomes very large:  $s_k = -1 * B_k^{-1} \nabla f$.  Then, this causes the next iterate of $B_k$ to be very large.  Eventually, the matrix becomes so large that np.linalg.inv(B_k) throws an error and will not be inverted.

How do I protect against this happening?  Do I need another check, and why isn’t my existing preventative measure working for this?


问:
Computing Multivariate Cumulative Normal Distribution
I am trying to compute the CDF of a Multivariate Cumulative Normal Distribution in the 1000th dimension (I have a 1000 vectors and their covariance matrix). I haven't been able to find a fast way to compute this, numpy and scipy all take take too long.

Are there any libraries (platform and language doesn't matter) that can compute this in a reasonable amount of time? If not, are there any numerical methods that can achieve this?


问:
Finding the smallest root of a function on $[0, \infty)$
I would like to find the smallest real root of a 1-D real-valued function $f(x)$ on the domain $x\in [0,\infty)$. In this problem, I can make the following guarantees on $f$:

$f$ does have a root at some finite, positive $x$.

$f$ is differentiable at all finite, positive $x$ and $f'(x)$ is either known or can be approximated accurately.

$f$ may have multiple roots (and usually does).

At its smallest root, $f$ changes sign from positive to negative.

Further, I can estimate

an approximate lower bound on the spacing between roots, call it $\Delta$.

an upper bound on the location of the smallest root, call it $x_{\max}$. If the root is beyond this point, I can compute its approximate location analytically.

My current algorithm for this problem is to bracket the root via a naive linear search. Starting from $x=0$, I march along in increments of $\Delta$ until $f(x)$ and $f(x+\Delta)$ differ in sign. With the root bracketed, I defer to a standard root-finding algorithm.

With the knowledge I have of $f$ and the estimates I am willing to make about its root structure, can I bracket the smallest root more quickly than a simple linear search?


问:
Transforming an arbitrary linear program into one with an interior point
Primal-dual interior point methods for linear programming require interior primal and dual starting points. I am looking for a good reference containing a description for modifying a given linear program (without solving an auxiliary LP and without doubling the number of variables or constraints) to one with a primal and dual interior point. The linear program is given in the form  $\min c^Tx$ s.t. $Ax\ge b$ (or $Ax=b$, $x\ge 0$) and may possibly be primal or dual infeasible.</p>


问:
PETSc SNES for user defined state
How to use PETSc SNES (scalable nonlinear equation solver), when the solution is not a vector but a user defined state?

I am implementing a non-linear mechanics problem (geometrically exact shell 5-paramter model), the update

$$X^{(k+1)} = X^{(k)} + \Delta_X$$

is not additive and requires the use of exponential maps.

More specifically, the state I mentioned above stores $ X $ :
$$ X = \left \{ \Phi , \Lambda  \right \} $$
$ \Phi $ (vector of size 3) is mid-surface coordinates and have direct updates i.e  $$ \Phi _{n+1} = \Phi_{n} + \Delta \Phi $$
While $ \Lambda $, the rotation matrix (correpondiong to director at mid surface), is updated through exponential map (the update increment is a vector of size 2, which is obtained after exploiting the one to one connection between tangent space of SO3 and S2 and fixing a reference vector).
The residual at  each point of mid surface(node) is a vector of size 5 and consists of $ \Lambda $ , $ \Phi $ and their derivatives.

I am using Newton-Raphson and many times the solution couldn't converge using a complete step size. Consequently, I thought of resorting to PETSc SNES (line search).

At any state, I know the residual vector and stiffness matrix.
I could not figure how to make sense of solution vector that I get from SNES solver (or to send data to SNES), given that I have my own update procedure.

Title: PETSc SNES for user defined state


问:
Can I MPI_Test for an Ibarrier
In this code, all processes post a barrier, sleep a while for good measure, then first Test and then Wait for the barrier. The Test says no and the wait succeeds. A test should be like a non-blocking wait, so it should succeed. What am I missing?

```<code>  printf("[%d] posting barrier, then sleeping\n",procid);
  MPI_Request final_barrier;
  MPI_Ibarrier(comm,&final_barrier);
  sleep(2);

  int all_done_flag=0;
  MPI_Test(&final_barrier,&all_done_flag,MPI_STATUS_IGNORE);
  printf("[%d] all done: %d\n",procid,all_done_flag);

  MPI_Wait(&final_barrier,MPI_STATUS_IGNORE);
  printf("[%d] concluded\n",procid);
</code>```

Title: Can I MPI_Test for an Ibarrier


问:
Simplest meaningful PDE/FEM calculation for mechanical stress due to heat
W have a complicated structure on which we do some FEM calculations regarding electrical potentials and heat distribution. The equations have the form

$\nabla\kappa\nabla u = f + g\rvert_{N}$

where $u$ is the temperature or electrical potential, $\kappa$ is a material parameter and $f$ and $g$ are heat or current densities that enter the investigated object over the whole volume given as a density-per-volume or over the surface $N$, given as density-per-area, respectively.

Now, we'd like to calculate the mechanical stress induced by heating of the material. In literature, this is usually described with vector valued functions and using nontrivial tensors. We wonder whether it is possible to formulate a meaningful calculation of the mechanical stress in the same form as the equation above, in order to re-use our code that we wrote to solve it. We are aware that this is not possible for all materials and larger stresses/displacements.

Is it possible to formulate a meaningful calculation of the mechanical stress in the same form as the equation above? For which type of materials and other restrictions would this form be valid?


问:
Discrete sine and cosine transform for mixed derivatives
Using sine and cosine transforms to solve Poisson's equation with Dirichlet boundary conditions seem quite standard nowadays (see, e.g., here or Table 2 in this paper). In the case of Poisson's equation, the multiplication with $-k^2$ preserves the symmetry of the Fourier coefficients such that the same transform can be used for going back to real space.

But what happens if the differential operator in Fourier space changes the symmetry? E.g., in the case of the transport equation, multiplication with $ik$ would transform even symmetry into odd, such that another transform would be needed for the backtransformation to real space. This should happen for all first or mixed derivatives. Are there any papers/codes around for these cases?


问:
Accelerate computation speed by using a different syntax
I was reading a book about Matlab ("Accelerating Matlab with GPU computing - A Primer with Examples" by Jung Suh and Youngming Kim, 2013. Chapter 1.7 Examples).

I read an example where it said that:

by using elemental calculation, matrix calculation, and function call
  instead of for loops, the speed can be dramatically accelerated.

In the case I read, it reduced calculation time from 50s to 0.5s.

My questions:

which syntax is more "resources saving" and "faster": function call, for, while, if, switch, and (i.e. search elements in a matrix or set) etc.?

How fast is the element-wise calculation, addition, and matrix calculation etc.?

For the convince of understanding, a short example was given below:

```<code>%%% First case
tic;
B = 0;
A = 1:1e8;
for i=1:1000
    B = B + A(i);
end
toc;

%%% Second case
tic;
B = 0;
A = 1:1e8;
B = sum(A);
toc;

%%% Third case
C = ones(1e8,1);
tic;
B = 0;
A = 1:1e8;
B = A*C;
toc;
</code>```

with the following timing results

Elapsed time is 0.463159 seconds.

Elapsed time is 0.501524 seconds.

Elapsed time is 0.532470 seconds.

In the above example, all three methods had different but comparable speed.


问:
Equal area algorithm to find shock location
I am looking to solve 1D burgers equation with various random initial conditions. What is the best algorithm to find the exact solution?

One method that is covered in literature is the equal area principle (R. J. LeVeque, Numerical Methods for Conservation Laws). Is this the standard procedure employed in numerics papers on limiters (for ex: this paper) for the exact solutions?

Can anyone provide a pseudo-code or algorithm to implement this method to find the shock location and hence the exact solution?


问:
Avoiding divergent solutions with `odeint`? shooting method
I am trying to solve an equation in Python. Basically what I want to do is to solve the equation:

$$
\frac{1}{x^2}\frac{d}{dx}\left(Gam \frac{dL}{dx}\right)+L\left(\frac{a^2x^2}{Gam}-m^2\right)=0
$$

This is the Klein-Gordon equation for a massive scalar field in a Schwarzschild spacetime. We know $m$ and $Gam=x^2-2x$.

The initial/boundary condition that I know are $L\rvert_{2+\epsilon}=1$ and $L\rvert_{\infty}=0$.

Notice that the asymptotic behavior of the equation is

$$
L(x\to\infty)\to \frac{e^{\pm(m^2-a^2)x}}{x}
$$

Then, if $a^2>m^2$ we will have oscillatory solutions, while if $a^2 < m^2$ we will have a divergent and a decay solution. What I am interested in is  the decay solution; however when I am trying to solve the above equation transforming it as a system of first-order differential equations and using the shooting method in order to find the $a$ that can give me the behavior that I am interested about, I am always having a divergent solution (for $0<a^2<m^2$). I suppose that it is happening because odeint is always finding the divergent asymptotic solution.

Is there a way to avoid or tell to odeint that I am interested in the decay solution? If not, do you know a way that I could solve this problem? Maybe using another method for solving my system of differential equations? If yes, which method?

Basically what I am doing is to add a new system of equations for $a$

$$\frac{d^2a}{dx^2}=0,\\
\frac{da}{dx}(2+\epsilon)=0,\\
a(2+\epsilon)=a_0)$$

in order to have $a$ as a constant. Then I am considering different values for $a_0$ and asking if my boundary conditions are fulfilled.

EDIT 1.

Hello, let me try to explain a little more my issue

I am incorporating the value at infinity considering the assimptotic behavior, it means that I will have a relation between the field and its derivative. I will post the code for you if it is helpful:

```<code>from IPython import get_ipython
get_ipython().magic('reset -sf')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from math import *
from scipy.integrate import ode
</code>```

These are initial conditions for Schwarzschild. The field is invariant under reescaling, then I can use $L(2+\epsilon)=1$

```<code>def init_sch(u_sch):
    om = u_sch[0]
    return np.array([1,0,om,0]) #conditions near the horizon, [L_c,dL/dx,a,da/dx]
</code>```

These are our system of equations

```<code>def F_sch(IC,r,rho_c,m,lam,l,j=0,mu=0):
    L = IC[0]
    ph = IC[1]
    om = IC[2]
    b = IC[3]

    Gam_sch=r**2.-2.*r

    dR_dr = ph
    dph_dr = (1./Gam_sch)*(2.*(1.-r)*ph+L*(l*(l+1.))-om**2.*r**4.*L/Gam_sch+(m**2.+lam*L**2.)*r**2.*L)
    dom_dr = b
    db_dr = 0.
    return [dR_dr,dph_dr,dom_dr,db_dr]
</code>```

Then I try for different values of "om" and ask if my boundary conditions are fulfilled. p_sch are the parameters of my model. In general what I want to do is a little more complicated and in general I will need more parameters that in the just massive case. Howeve I need to start with the easiest which is what I am asking here

```<code>p_sch = (1,1,0,0) #[rho_c,m,lam,l], lam and l are for a more complicated case
ep = 0.2
ep_r = 0.01
r_end = 500
n_r = 500000
n_omega = 1000
omega = np.linspace(p_sch[1]-ep,p_sch[1],n_omega)
r = np.linspace(2+ep_r,r_end,n_r)
tol = 0.01
a = 0

for j in range(len(omega)):
    print('trying with $omega =$',omega[j])
    omeg = [omega[j]]
    ini = init_sch(omeg)
    Y = odeint(F_sch,ini,r,p_sch,mxstep=50000000)
    print Y[-1,0]
    #This is my boundary condition that I am implementing. Basically this should be my conditions at "infinity
    if abs(Y[-1,0]*((p_sch[1]**2.-Y[-1,2]**2.)**(1/2.)+1./(r[-1]))+Y[-1,1]) < tol:
        print(j,'times iterations in omega')
        print("R'(inf)) = ", Y[-1,0])
        print("\omega",omega[j])
        omega_1 = [omega[j]]
        a = 10
        break
    if a > 1:
        break
</code>```

Basically what I want to do here is to solve the system of equations giving different initial conditions and find a value for "a=" (or "om" in the code) that should be near to my boundary conditions. I need this because after this I can give such initial guest to a secant method and try to fiend a best value for "a". However, always that I am running this code I am having divergent solutions that it is, of course, a behavior that I am not interested. I am trying the same but considering the scipy.integrate.solve_vbp, but when I run the following code:

```<code>from IPython import get_ipython
get_ipython().magic('reset -sf')
import numpy as np
import matplotlib.pyplot as plt
from math import *
from scipy.integrate import solve_bvp

def bc(ya,yb,p_sch):
    m = p_sch[1]
    om = p_sch[4]
    tol_s = p_sch[5]
    r_end = p_sch[6]

    return np.array([ya[0]-1,yb[0]-tol_s,ya[1],yb[1]+((m**2-yb[2]**2)**(1/2)+1/r_end)*yb[0],ya[2]-om,yb[2]-om,ya[3],yb[3]])

def fun(r,y,p_sch):
    rho_c = p_sch[0]
    m = p_sch[1]
    lam = p_sch[2]
    l = p_sch[3]

    L = y[0]
    ph = y[1]
    om = y[2]
    b = y[3]

    Gam_sch=r**2.-2.*r

    dR_dr = ph
    dph_dr = (1./Gam_sch)*(2.*(1.-r)*ph+L*(l*(l+1.))-om**2.*r**4.*L/Gam_sch+(m**2.+lam*L**2.)*r**2.*L)
    dom_dr = b
    db_dr = 0.*y[3]
    return np.vstack((dR_dr,dph_dr,dom_dr,db_dr))

eps_r=0.01
r_end = 500
n_r = 50000
r = np.linspace(2+eps_r,r_end,n_r)
y = np.zeros((4,r.size))
y[0]=1
tol_s = 0.0001
p_sch= (1,1,0,0,0.8,tol_s,r_end)

sol = solve_bvp(fun,bc, r, y, p_sch)
</code>```

I am obtaining this error: ValueError: bc return is expected to have shape (11,), but actually has (8,).
ValueError: bc return is expected to have shape (11,), but actually has (8,).


问:
Geometric interpretation of lemma
I am currently studying eigenvalue problems. I already worked through the minimax-principles and seen why $\lambda_{h, m} \geq \lambda_m$, when comparing the eigenvalues of a discretization and the continuous ones.

The next lemma tries to establish conditions when we can get convergence between these two, finally arriving at the theorem that eigenvalues converge twice as fast as eigenvectors do. The lemma states the following with $V$ being the space we are looking for eigenvectors in, $V_h$ being a discretation of that and $V_m$ being the space which is spanned by the first $m$ eigenvectors:

Define
$$C_{h,m} := \inf_{v \in V_m} \frac{\| P_h v \|}{\|v\|}, \quad 1 \leq m \leq N.$$
If there exists a $m$ such that $C_{h,m} > 0$, we get
$$\lambda_m \leq \lambda_{h,m} \leq C_{h,m}^{-2} \, \lambda_m.$$

There is a remark stating that the $1-C^{-1}_{h,m}$ is a measure for how close $V_m$ and $P_h V_m$ are and that it is the norm of the inverse of the operator $P_h$.  Can you help me how to arrive at this conclusion?

My thoughts so far were that the $C_{h,m}$ are measuring the minimal $\cos(\theta)$ between $V_m$ and $P_h V_m$. As we want to get $C_{h,m} \to 1$ later on for $h \to 0$, it means that $P_h V_m$ would become parallel to $V_m$. What is the connection to the inverse of the operator $P_h$?

Title: Geometric interpretation of lemma


问:
How to numerically transform a 2D Fourier spectrum with arbitrary frequency shift to center frequency?
Suppose $F(u,v)$ is the center frequency Fourier representation of some $f(x,y)$ in 2D.

$$
f(x,y)=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}F(u,v)e^{2\pi i (xu+yv)}dudv
$$

In principle for the same $f(x,y)$, the representation could also be $F(u-u_0, v-v_0)$, where $u_0$ and $v_0$ are arbitrarily shifted frequencies.

Now assume the fact is we already have a representation $F(u-u_0, v-v_0)$ and coefficients(arranged as a matrix) and also we know the shift $u_0$ and $v_0$, how can we derive the coefficients of $F(u,v)$ corresponding to the same $f(x,y)$?


问:
How to model a non-linear least-squares problem for triangles
I have a non-linear least-squares problem to solve and with my current modeling, the solver is either very slow or does not converge to a correct solution.

For the problem, I need to minimize energy on a set of triangles, which is defined by

$$\sum_f E_d(f) + \sum_e E_s(e)$$
with faces given by their vertices $f=((u_1,v_1), (u_2,v_2), (u_3,v_3))$
and shared edges given by the vertices of both adjacent triangles $e=((u_1, v_1)_{T_{1}}, (u_2, v_2)_{T_{1}}, (u_2, v_2)_{T_{2}}, (u_2, v_2)_{T_{2}})$, where 1,2,3 are local indices for the triangle/edge and the vertices are from a common set $X = \{(u_1, v_1), \dots, (u_{|V|}, v_{|V|})\}$.

The edge indices are defined such that for the inner edges of a manifold mesh holds $(u_1, v_1)_{T_{1}}=(u_1, v_1)_{T_{2}},\ (u_2, v_2)_{T_{1}}=(u_2, v_2)_{T_{2}}$.

I am using the ceres solver using the line search algorithm and modeling the problem as two residual blocks of size $|F|$ and $|E|$, defined by

$$
\vec{E}_f(\vec{X}) = (E_f(f_1), \dots, E_f(f_{|F|})) \\
\vec{E}_s(\vec{X}) = (E_s(e_1), \dots, E_s(e_{|E|}))
$$

This works for many inputs but is converging slowly and sometimes to local minima.

To model it using more blocks, I tried to create one block per triangle and one block per edge, leading to one block per face and one block per edge having only one residual

$$
E_{f_i}(\vec{X}) = E_f(f_i) \\
E_{e_j}(\vec{X}) = E_s(e_j) \\
$$

where $E_{f_i}$ and $E_{e_j}$ are only depending on the 6 respectively 8 values which are relevant for their energy, which should result in a sparse problem which is hopefully fast to solve.

As far as I understood the problem modeling in the ceres solver, modeling the problem this way should be possible and probably the faster way to solve it, but I get totally different results and the result of the modeling using separate blocks seems completely wrong.

What can be my problem here and what is the best way to model such a problem for fast solving?

(The problem is based on the methods in the autocuts parametrization paper. The solving process there is iterated a few times sharpening a smoothing term in $E_s$, but the question is only about a single iteration.)


问:
A good 2D finite difference for the continuity equation
How could I go about solving the continuity equation below in 2D?

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho u\right)=0$$

I saw that a similar question was posted here: A good finite difference for the continuity equation, but that only discussed the 1-D case.

In it, one of the answers defines incoming, and outgoing fluxes, as
$$\Phi_{i+1/2} = \dfrac{u_{i+1/2}+|u_{i+1/2}|}{2}\rho_{i}
              + \dfrac{u_{i+1/2}-|u_{i+1/2}|}{2}\rho_{i+1}\\
\Phi_{i-1/2} = \dfrac{u_{i-1/2}+|u_{i-1/2}|}{2}\rho_{i-1}
              + \dfrac{u_{i-1/2}-|u_{i-1/2}|}{2}\rho_{i}$$

With these, one could use the finite difference scheme

$$\dfrac{\partial(\rho)}{\partial t}=-\dfrac{\partial(\rho u)}{\partial x} = -\dfrac{\partial \Phi}{\partial x} \approx -\dfrac{\Delta \Phi}{\Delta x} \approx -\dfrac{\Phi_{i+1/2}-\Phi_{i-1/2}}{\Delta x}$$

The problem is that when I use a naive generalization of this, namely defining

$$\Phi_{i+1/2,j;x} = \dfrac{u_{i+1/2,j;x}+|u_{i+1/2,j;x}|}{2}\rho_{i,j}
              + \dfrac{u_{i+1/2,j;x}-|u_{i+1/2,j;x}|}{2}\rho_{i+1,j}\\
\Phi_{i,j+1/2;y} = \dfrac{u_{i,j+1/2;y}+|u_{i,j+1/2;y}|}{2}\rho_{i,j}
              + \dfrac{u_{i,j+1/2;y}-|u_{i,j+1/2;y}|}{2}\rho_{i,j+1}$$

and taking

$$\dfrac{\partial (\rho u_x)}{\partial x}  \approx \frac{\Phi_{i+1/2,j;x}-\Phi_{i-1/2,j;x}}{\Delta x}$$
and
$$\dfrac{\partial (\rho u_y)}{\partial y}  \approx \frac{\Phi_{i,j+1/2;y}-\Phi_{i,j-1/2;y}}{\Delta y}$$
my density diverges. What scheme could I use to solve this equation properly in 2D?


问:
What is the mathematical and physical principle behind of RBE2 element？
I am writing a 3d linear finite element code to solve the standard linear elasticity equation on a tetrahedron mesh of a gearbox. Notice that, the two rectangular plates above the gearbox are fixed, which means their displacement are zeros.

The standard linear elasticity model is
$$
-\mu\Delta \bf u - (\lambda + \mu)\nabla \text{div} \bf u = \bf f
$$

where $\bf u$ is the displacement, $\mu$ and $\lambda$ are the Lamé parameters. Then I get the weak formulation as follows:

$$
\int_\Omega \sigma(u):\varepsilon \rm(v) dx = \int_\Omega fv\rm d x
$$

where $\sigma$ is the stress and $\varepsilon$ is the strain. I use linear finite element to discrete above weak formulation and get

$$
KU=F
$$

where $K$ is the stiffness matrix, $U$ is the displacement vector, and $F$ is the loads.

At last, I want to use RBE2 element to add some constraints to the above finite element model. I add an extra node $A$ to the tetrahedron mesh, and connect it to a boundary node sets $\{B_i\}_{i=1}^N$ of the tetrahedron mesh. Then I apply a load on node $A$, include a force $F$ and a moment $M$.

I have tried my best to search a document which can explain the details of the mathematical and physical of RBE2 element clearly but got nothing. So I need your suggestion and help, thanks very much.

One can take a look at this document, which explain RBE3 element very clearly but the RBE2 element is still not clear for my case.

Edit:
For RBE2 element, the $A$ is an independent node, which displacements (transition and rotation) are unknown. Furthermore, the displacements of $\{B_i\}_{i=1}^N$ are also unknown. I just know the concentrated load $F$ and moment $M$ on $A$, so my question is how to apply the $F$ and $M$ into $KU=F$?


问:
Which solvers for BVP in python are the best? Is there something better that scipy.integrate.solve_bvp?
I am trying to solve a boundary value problem with Python. I have been using scipy.integrate.solve_bvp but the result that it is giving me is completely wrong. Basically my code is as follows:

```<code>from IPython import get_ipython
get_ipython().magic('reset -sf')
import numpy as np
import matplotlib.pyplot as plt
from math import *
from scipy.integrate import solve_bvp

r_in = 0.0
r_end = 500.
l = 0.
rho_c = 1.
m = 1.

def system_DE(r,y,p_self,l=l,rho_c=rho_c):
     R = y[0]
     ph = y[1]
     a = y[2]
     al = y[3]

     om = p_self[0]
     u = p_self[1]

    dR_dr = ph
    da_dr = a*(((2.*l+1.)*r/2.)*(-om**2.*a**2.*R**2./al**2.+ph**2.+l*(l+1.)*a**2.*R**2./r**2.+m**2.*a**2.*R**2.)-(a**2.-1.)/(2.*r))
    dal_dr = al*(da_dr/a-l*(l+1.)*(2.*l+1.)*a**2.*R**2./r-(2.*l+1.)*m**2.*a**2.*r*R**2.+(a**2.-1.)/r)
    dph_dr = -2.*ph/r+(l*(l+1.)*(2.*l+1.)*a**2.*R**2./r+(2.*l+1.)*m**2*a**2*r*R**2-(a**2-1)/r)*ph-om**2.*a**2.*R/al**2.+l*(l+1.)*a**2.*R/r**2.+m**2.*a**2.*R

     return [dR_dr,dph_dr,da_dr,dal_dr]

def bc(ya,yb,p_sch,l=l,rho_c=rho_c,r_in=r_in,r_end = r_end,m=m):

    om = p_self[0]
    u = p_self[1]

    if l==0.:
        return np.array([ya[0]-rho_c,ya[1]-rho_c*r_in,yb[0]*((m**2.-om**2.)**(1./2.)+1./(r_end))+yb[1],ya[2]-1.,ya[3]-u,yb[3]*yb[2]-1.])
    else:
        return np.array([rho_c*r**l,rho_c*l*r**(l-1.),1.,u])

n_r = 1000
r = np.linspace(r_in,r_end,n_r)
y = np.zeros((4,r.size))
y[0,0] = rho_c
y[0,1] = rho_c

y[2,:] = 1.

p_self=[0.8,0.85] #[om,u]

sol = solve_bvp(system_DE,bc, r, y, p=p_self, tol=0.01,max_nodes=100000)

x_plot = np.linspace(r_in, 200  , 5000)
L = sol.sol(x_plot)[0]
Ph = sol.sol(x_plot)[1]
a = sol.sol(x_plot)[2]
alpha = sol.sol(x_plot)[3]

plt.subplot(211)
plt.plot(x_plot, L)
plt.xlabel("x")
plt.ylabel("y")
plt.subplot(212)
plt.plot(x_plot, a)
plt.plot(x_plot, alpha)
plt.xlabel("x")
plt.ylabel("y")
plt.show()
</code>```

The idea is that with the parameters that I am giving to the code I should obtain something similar to the figure [1] (up) on this article: https://arxiv.org/pdf/0908.2435.pdf, however if I run the code I am obtaining something completely wrong. Is there another solver that I could use? Is there something wrong with my code?

Thanks for your time. Regards,

Luis P.


问:
How do I configure PETSc to run long double precision or some other precision that is greater than default?
I'm installing PETSc for complex numbers with the C99 standard. I'd like to have it installed using a higher precision than the default (double, I presume) since that will most likely make the LU preconditioner run better. I've tried several different commands that I could google, but it doesn't work. I've searched the manual but haven't found anything on the subject.

I tried invoking ./configure with the

precision=longdobule

and

precision=__float128

with

download-f-blas-lapack

and

download-f2cblaslapack

in all four different combinations. The one time I used download-f2cblaslapack and precision=__float128 it did configure successfully but at the next stage when I compiled the libraries I got the following error for every library being compiled:

error: unknown type name ‘PetscScalar’

What do I have to write to have it configured at a higher precision?


问:
Review of modern homotopy methods and practical techniques
I'm hoping someone can recommend recent literature concerning homotopy methods for solving systems of nonlinear equations.  Already by the time of Layne Watson's 1986 paper there were a lot of methods, and I'm sure there have been developments since then as well as practical shifts in the field due to hardware advancements.

I'm looking to apply globally convergent methods to problems in nonlinear elasticity with ~$10^6$ degrees of freedom.  The Jacobians are sparse but local methods like Newton-Raphson or quasi-Newton with line search or trust region tend to get trapped at points where the Jacobian is singular.  I've had some success with a crude implementation of a homotopy method but I'm looking to get better performance and better familiarity with modern methods.

Edit: I'm mainly interested in the benefits and drawbacks of the three methods for curve tracking that Watson describes: the ODE algorithm, the normal flow algorithm, and the augmented Jacobian algorithm.  Watson's papers are very nice but they are a bit outdated for assessing my problem.  For example, it seems that the limited memory version of BFGS had not been developed yet, which made it impossible to fairly test a sparse version of the augmented Jacobian algorithm.  For my problem, I am only interested in the solution, i.e. $x$ s.t. $F(x,\lambda)=0$ for $\lambda=1$, but I am very concerned with stability and efficiency.

The linear algebra steps seem quite straightforward, and I can test different preconditioners, etc.


问:
What is the finite-difference representation of the Laplacian operator with periodic boundary conditions?
I am using a central-difference scheme to solve the eigenvalue problem

$$\frac{d^2}{dx^2}u = \lambda u$$

on a unit interval with periodic boundary conditions. My understanding is the eigenvectors $u_i$ have the analytic form

$$u_i(x) = e^{in\pi x}$$

I know that, for example, with dirichlet boundary conditions $u_0 = u_{N+1} = 0$ and $N = 6$, the central-difference representation of the Laplacian is

$$
\frac{1}{h^2}\begin{bmatrix}
    2       & -1 & 0 & 0 & 0 & 0 \\
    -1       & 2 & -1 & 0 & 0 & 0 \\
    0 & -1 & 2 & -1 &  0 & 0\\
    0 &0 & -1 & 2 & -1 &  0\\
    0 & 0 & 0 & -1 & 2 & -1 \\
    0 & 0 & 0 & 0 & -1 & 2 \\
\end{bmatrix}
$$
But I'm not sure about with periodic boundary conditions. My best guess so far is

$$
\frac{1}{h^2}\begin{bmatrix}
    2       & -1 & 0 & 0 & 0 & -1 \\
    -1       & 2 & -1 & 0 & 0 & 0 \\
    0 & -1 & 2 & -1 &  0 & 0\\
    0 &0 & -1 & 2 & -1 &  0\\
    0 & 0 & 0 & -1 & 2 & -1 \\
    -1 & 0 & 0 & 0 & -1 & 2 \\
\end{bmatrix}
$$

But this doesn't give me the expected eigenvectors.

My question: Is my form above correct?


问:
Automatic differentiation of barycentric rational functions
By a barycentric rational interpolant we understand a function of the form

\begin{align*}
r(t) := \frac{\sum_{i=0}^{n-1} \frac{w_i y_i}{t-t_i} }{ \sum_{i=0}^{n-1} \frac{w_i}{t-t_i}}
\end{align*}
In Schneider and Werner, (proposition 12) a stable algorithm for computation of the derivatives of $r^{(k)}$ is presented. (This algorithm is reviewed here  in more modern notation, in equation 5.1a.) However, this algorithm requires two passes through the data vectors $\{w_{i}\}_{i=0}^{n-1}$, $\{y_{i}\}_{i=0}^{n-1}$ and $\{t_{i}\}_{i=0}^{n-1}$.

Since I only require the first derivative $r'$, I am curious as to whether I can get the first derivative via automatic differentiation, in one pass, perhaps increasing the speed of evaluation. So I have two questions: What is the automatic differentiation formula for $r$, and by using it, will I lose the stability guaranteed by Werner and Schneider's algorithm?


问:
Preconditioned residual converges, but true residual doesn't
I'm using Albany w/ Trilinos to solve an elasticity problem with thermal expansion mismatch. I'm using block GMRES with MueLu preconditioning. It works for problem size of several million dofs, but when the problem size gets bigger, the true residual won't converge even though the preconditioned residual converges. I've tried different parameters for MueLu (sym Gauss-Seidel, Jacobi, ...), but it won't work. Any suggestions on the convergence problem?

For the problem size that converges with Block GMRES with MueLu, I tried Block GMRES with ILU preconditioner, and the true residual norm stays 1. Does that mean I have a singular preconditioner with ILU? Why did that happen?

Thanks.


问:
Convergence of a very large non-linear least squares optimization
(note: I also posted this question on stackoverflow  before finding this community here, which seems a better place for it)

I'm trying to solve the following problem: I have a lot (~80000) surface patches of an organ that's growing. I measure each of its areas over time (18 time-points) and want to fit a growth curve to it (bi-logistic model, eg. just the sum of two logistic functions bcs. there are two 'growth spurts' happening in the observed period).

I have box constraints to ensure that the exponential terms don't explode and a linear constraint that one growth spurt has to happen after the other. Also, in order to enforce some sort of spatial continuity in the fitted parameters, I add a penalty term to the objective function (least squares) given as the squared sum of differences between (some of the) parameters of neighbouring patches, so the individual model fits are not independent anymore.

I have gradients and the hessian of the whole thing, and I solve it at increasingly finer scales, starting with the overall surface area, subdividing, mapping the parameters of the global model fit to each patch, running the solver again, subdividing ... until I get my desired resolution.

So, optimization using IPOPT works, but it is terribly slow, and since I'm a bit of a noob in all things concerning optimization theory I was wondering if I'm doing something terribly stupid in the settings. I'm using ma86 as linear solver compiled against openBLAS running on a machine with 256G RAM and 56 cores as well as metis reordering. The other settings I use are:

```<code>% nlp scaling
solv_options.ipopt.nlp_scaling_method = 'gradient-based';
solv_options.ipopt.nlp_scaling_max_gradient = 1;
solv_options.ipopt.nlp_scaling_min_value = 1e-16;

solv_options.ipopt.bound_mult_init_method='constant';

% Barrier Parameter
solv_options.ipopt.mu_strategy = 'adaptive';
solv_options.ipopt.mu_oracle = 'quality-function';
solv_options.ipopt.fixed_mu_oracle = 'average_compl';
solv_options.ipopt.adaptive_mu_globalization = 'kkt-error';
solv_options.ipopt.corrector_type = 'affine';

% linear solver
solv_options.ipopt.max_soc=0;
solv_options.ipopt.accept_every_trial_step='yes';
solv_options.ipopt.linear_system_scaling = 'none';
solv_options.ipopt.neg_curv_test_tol = 0;
solv_options.ipopt.neg_curv_test_reg = 'yes';
solv_options.ipopt.max_refinement_steps=0;
solv_options.ipopt.min_refinement_steps=0;

% ma86 settings
solv_options.ipopt.linear_solver='ma86';
solv_options.ipopt.ma86_order='auto';
solv_options.ipopt.ma86_scaling='mc64';
solv_options.ipopt.ma86_small=1e-10;
solv_options.ipopt.ma86_static=1;
solv_options.ipopt.recalc_y='yes';
</code>```

and this gives me something like this:

```This is Ipopt version 3.12, running with linear solver ma86.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:     2560
Number of nonzeros in Lagrangian Hessian.............:   112280

Total number of variables............................:     8960
                     variables with only lower bounds:     2560
                variables with lower and upper bounds:     3840
                     variables with only upper bounds:     2560
Total number of equality constraints.................:        0
Total number of inequality constraints...............:     1280
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:     1280

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  9.8260736e-01 0.00e+00 1.19e-02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.6288666e-01 0.00e+00 1.24e-02  -4.5 1.14e-02   0.0 9.91e-01 1.00e+00f  1
   2  9.1582880e-01 0.00e+00 1.16e-02  -4.4 2.72e-02  -0.5 1.00e+00 1.00e+00f  1
   3  8.2635857e-01 0.00e+00 1.01e-02  -4.6 1.39e-01  -1.0 9.99e-01 1.00e+00f  1
   4  7.8943781e-01 0.00e+00 9.40e-03  -4.9 2.70e-02  -0.5 1.00e+00 1.00e+00f  1
   5  7.2123624e-01 0.00e+00 8.12e-03  -5.3 8.70e-02  -1.0 1.00e+00 1.00e+00f  1
   6  6.9535003e-01 0.00e+00 7.56e-03  -6.1 2.20e-02  -0.6 1.00e+00 9.06e-01f  1
   7  6.6635914e-01 0.00e+00 7.00e-03  -6.6 6.21e-02  -1.1 1.00e+00 5.40e-01f  1
   8  6.5683787e-01 0.00e+00 6.78e-03  -7.6 2.14e-02  -0.6 1.00e+00 3.81e-01f  1
   9  6.4238130e-01 0.00e+00 6.53e-03  -7.7 1.53e-01  -1.1 1.00e+00 2.90e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  6.3440625e-01 0.00e+00 6.33e-03  -8.7 2.20e-02  -0.7 1.00e+00 3.30e-01f  1
  11  6.2431009e-01 0.00e+00 6.16e-03  -8.3 3.08e-01  -1.2 1.00e+00 2.04e-01f  1
  12  6.1872460e-01 0.00e+00 6.05e-03  -8.8 2.29e-02  -0.7 1.00e+00 2.32e-01f  1
  13  6.0753815e-01 0.00e+00 5.86e-03  -8.9 2.73e-01  -1.2 1.00e+00 2.31e-01f  1
  14  6.0575477e-01 0.00e+00 5.82e-03  -9.9 2.60e-02  -0.8 1.00e+00 7.44e-02f  1
  15  6.0089103e-01 0.00e+00 5.71e-03 -11.0 1.11e-02  -0.4 1.00e+00 4.48e-01f  1
  16  5.9426852e-01 0.00e+00 5.58e-03 -11.0 3.14e-02  -0.8 1.00e+00 2.68e-01f  1
  17  5.9175418e-01 0.00e+00 5.52e-03 -11.0 1.17e-02  -0.4 1.00e+00 2.26e-01f  1
  18  5.8756155e-01 0.00e+00 5.44e-03 -11.0 3.97e-02  -0.9 1.00e+00 1.65e-01f  1
  19  5.8651702e-01 0.00e+00 5.42e-03  -9.5 1.26e-02  -0.5 1.00e+00 8.89e-02f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  5.8278151e-01 0.00e+00 5.34e-03 -10.3 5.18e-02  -0.9 1.00e+00 1.40e-01f  1
  21  5.8100614e-01 0.00e+00 5.30e-03 -11.0 1.37e-02  -0.5 1.00e+00 1.42e-01f  1
  22  5.7827149e-01 0.00e+00 5.25e-03 -11.0 7.12e-02  -1.0 1.00e+00 9.77e-02f  1
  23  5.7564914e-01 0.00e+00 5.19e-03  -9.4 1.52e-02  -0.6 1.00e+00 1.97e-01f  1
  24  5.7340095e-01 0.00e+00 5.15e-03 -10.0 1.06e-01  -1.0 1.00e+00 7.66e-02f  1
  25  5.7134734e-01 0.00e+00 5.10e-03 -11.0 1.78e-02  -0.6 1.00e+00 1.45e-01f  1
  26  5.6897098e-01 0.00e+00 5.06e-03 -11.0 1.82e-01  -1.1 1.00e+00 7.65e-02f  1
  27  5.6722484e-01 0.00e+00 5.02e-03 -11.0 2.11e-02  -0.7 1.00e+00 1.16e-01f  1
  28  5.6402838e-01 0.00e+00 5.25e-03 -11.0 3.29e-01  -1.1 1.00e+00 9.57e-02f  1
  29  5.6299118e-01 0.00e+00 5.26e-03 -10.5 2.38e-02  -0.7 1.00e+00 6.53e-02f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  5.5753178e-01 0.00e+00 5.63e-03 -11.0 1.38e-01  -1.2 1.00e+00 1.62e-01f  1
  31  5.5730217e-01 0.00e+00 5.63e-03 -11.0 2.37e-02  -0.8 1.00e+00 1.38e-02f  1
  32  5.5363398e-01 0.00e+00 5.71e-03 -11.0 8.55e-02  -1.2 1.00e+00 1.06e-01f  1
  33  5.5223091e-01 0.00e+00 5.70e-03 -11.0 2.43e-02  -0.8 1.00e+00 7.99e-02f  1
  34  5.5163441e-01 0.00e+00 5.70e-03  -9.9 1.30e-01  -1.3 1.00e+00 1.65e-02f  1
  35  5.4699694e-01 0.00e+00 5.67e-03 -10.9 2.56e-02  -0.9 1.00e+00 2.48e-01f  1
  36  5.4559780e-01 0.00e+00 5.67e-03 -11.0 2.30e-01  -1.3 1.00e+00 3.71e-02f  1
  37  5.4248376e-01 0.00e+00 5.63e-03 -11.0 2.58e-02  -0.9 1.00e+00 1.60e-01f  1
  38  5.3955718e-01 0.00e+00 5.63e-03 -10.7 6.78e-01  -1.4 1.00e+00 6.97e-02f  1
  39  5.3606672e-01 0.00e+00 5.57e-03 -10.9 2.79e-02  -1.0 1.00e+00 1.73e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  5.3481638e-01 0.00e+00 5.55e-03 -11.0 1.32e-02  -0.5 1.00e+00 1.33e-01f  1
  41  5.3248177e-01 0.00e+00 5.50e-03 -11.0 3.23e-02  -1.0 1.00e+00 1.10e-01f  1
  42  5.3143109e-01 0.00e+00 5.48e-03 -10.4 1.37e-02  -0.6 1.00e+00 1.05e-01f  1
  43  5.2975978e-01 0.00e+00 5.45e-03 -11.0 3.92e-02  -1.1 1.00e+00 7.50e-02f  1
  44  5.2855042e-01 0.00e+00 5.42e-03 -11.0 1.50e-02  -0.6 1.00e+00 1.13e-01f  1
  45  5.2670676e-01 0.00e+00 5.38e-03 -11.0 4.80e-02  -1.1 1.00e+00 7.85e-02f  1
  46  5.2428812e-01 0.00e+00 5.31e-03 -11.0 1.64e-02  -0.7 1.00e+00 2.12e-01f  1
  47  5.2221881e-01 0.00e+00 5.27e-03 -11.0 5.84e-02  -1.2 1.00e+00 8.45e-02f  1
  48  5.2084516e-01 0.00e+00 5.23e-03 -10.8 1.79e-02  -0.8 1.00e+00 1.14e-01f  1
  49  5.1928549e-01 0.00e+00 5.19e-03 -11.0 7.06e-02  -1.2 1.00e+00 6.06e-02f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  5.1787195e-01 0.00e+00 5.15e-03 -11.0 1.95e-02  -0.8 1.00e+00 1.10e-01f  1
  51  5.1639782e-01 0.00e+00 5.12e-03 -10.9 8.49e-02  -1.3 1.00e+00 5.45e-02f  1
  52  5.1606363e-01 0.00e+00 5.11e-03 -11.0 2.13e-02  -0.9 1.00e+00 2.45e-02f  1
  53  5.1178863e-01 0.00e+00 5.01e-03 -11.0 1.04e-01  -1.3 1.00e+00 1.50e-01f  1
  54  5.1173718e-01 0.00e+00 5.01e-03 -11.0 2.33e-02  -0.9 1.00e+00 3.61e-03f  1
  55  5.0944130e-01 0.00e+00 4.96e-03 -11.0 1.40e-01  -1.4 1.00e+00 7.70e-02f  1
  56  5.0727659e-01 0.00e+00 4.89e-03 -11.0 2.55e-02  -1.0 1.00e+00 1.43e-01f  1
  57  5.0638483e-01 0.00e+00 4.85e-03 -11.0 1.11e-02  -0.5 1.00e+00 1.29e-01f  1
  58  5.0426908e-01 0.00e+00 4.78e-03 -11.0 2.79e-02  -1.0 1.00e+00 1.32e-01f  1
  59  5.0333858e-01 0.00e+00 4.74e-03 -11.0 1.22e-02  -0.6 1.00e+00 1.26e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  60  5.0161958e-01 0.00e+00 4.69e-03 -11.0 3.05e-02  -1.1 1.00e+00 1.02e-01f  1
  61  5.0037437e-01 0.00e+00 4.64e-03 -11.0 1.34e-02  -0.6 1.00e+00 1.58e-01f  1
  62  4.9971029e-01 0.00e+00 4.62e-03 -11.0 3.34e-02  -1.1 1.00e+00 3.72e-02f  1
  63  4.9914235e-01 0.00e+00 4.59e-03 -11.0 1.48e-02  -0.7 1.00e+00 6.69e-02f  1
  64  4.9664256e-01 0.00e+00 4.51e-03 -11.0 3.70e-02  -1.2 1.00e+00 1.32e-01f  1
  65  4.9621567e-01 0.00e+00 4.50e-03 -11.0 1.62e-02  -0.7 1.00e+00 4.72e-02f  1
  66  4.9453149e-01 0.00e+00 4.44e-03 -11.0 4.19e-02  -1.2 1.00e+00 8.39e-02f  1
  67  4.9304626e-01 0.00e+00 4.38e-03 -11.0 1.79e-02  -0.8 1.00e+00 1.53e-01f  1
  68  4.9163087e-01 0.00e+00 4.34e-03 -11.0 6.89e-02  -1.3 1.00e+00 6.69e-02f  1
  69  4.9031174e-01 0.00e+00 4.28e-03 -11.0 1.97e-02  -0.8 1.00e+00 1.28e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  70  4.8795240e-01 0.00e+00 4.22e-03  -8.8 2.01e-01  -1.3 9.38e-01 1.04e-01f  1
  71  4.8721723e-01 0.00e+00 4.18e-03  -9.8 2.16e-02  -0.9 1.00e+00 6.73e-02f  1
  72  4.8560412e-01 0.00e+00 4.10e-03 -11.0 9.08e-03  -0.5 1.00e+00 3.33e-01f  1
  73  4.8395604e-01 0.00e+00 4.03e-03 -11.0 2.38e-02  -0.9 1.00e+00 1.42e-01f  1
  74  4.8308417e-01 0.00e+00 3.99e-03 -11.0 9.99e-03  -0.5 1.00e+00 1.69e-01f  1
  75  4.8236398e-01 0.00e+00 3.96e-03 -10.8 2.63e-02  -1.0 1.00e+00 5.83e-02f  1
  76  4.7919513e-01 0.00e+00 3.79e-03 -10.9 1.11e-02  -0.6 1.00e+00 5.72e-01f  1
  77  4.7844811e-01 0.00e+00 3.76e-03 -11.0 3.12e-02  -1.0 1.00e+00 5.79e-02f  1
  78  4.7801338e-01 0.00e+00 3.74e-03 -11.0 1.22e-02  -0.6 1.00e+00 7.40e-02f  1
  79  4.7631089e-01 0.00e+00 3.67e-03 -11.0 4.30e-02  -1.1 1.00e+00 1.23e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  80  4.7516941e-01 0.00e+00 3.61e-03 -11.0 1.35e-02  -0.7 1.00e+00 1.81e-01f  1
  81  4.7405503e-01 0.00e+00 3.57e-03 -11.0 7.24e-02  -1.1 1.00e+00 7.54e-02f  1
  82  4.7394863e-01 0.00e+00 3.57e-03 -11.0 1.49e-02  -0.7 1.00e+00 1.57e-02f  1
  83  4.7209146e-01 0.00e+00 3.51e-03  -8.9 1.76e-01  -1.2 9.65e-01 1.15e-01f  1
  84  4.7199185e-01 0.00e+00 3.50e-03 -10.1 1.66e-02  -0.8 1.00e+00 1.36e-02f  1
  85  4.7051045e-01 0.00e+00 3.41e-03 -11.0 6.70e-03  -0.3 1.00e+00 4.74e-01f  1
  86  4.7028360e-01 0.00e+00 3.40e-03 -11.0 1.92e-02  -0.8 1.00e+00 2.87e-02f  1
  87  4.6892068e-01 0.00e+00 3.31e-03 -11.0 7.44e-03  -0.4 1.00e+00 4.02e-01f  1
  88  4.6754658e-01 0.00e+00 3.25e-03 -11.0 2.37e-02  -0.9 1.00e+00 1.62e-01f  1
  89  4.6723856e-01 0.00e+00 3.23e-03 -11.0 8.24e-03  -0.4 1.00e+00 8.47e-02f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  90  4.6625331e-01 0.00e+00 3.18e-03 -11.0 3.07e-02  -0.9 1.00e+00 1.08e-01f  1
  91  4.6621348e-01 0.00e+00 3.18e-03 -11.0 9.16e-03  -0.5 1.00e+00 1.01e-02f  1
  92  4.6511327e-01 0.00e+00 3.17e-03 -11.0 4.12e-02  -1.0 1.00e+00 1.11e-01f  1
  93  4.6363418e-01 0.00e+00 3.14e-03 -11.0 1.02e-02  -0.5 1.00e+00 3.45e-01f  1
  94  4.6314545e-01 0.00e+00 3.14e-03 -11.0 6.17e-02  -1.0 1.00e+00 4.60e-02f  1
  95  4.6235630e-01 0.00e+00 3.12e-03 -11.0 1.13e-02  -0.6 1.00e+00 1.70e-01f  1
  96  4.6143809e-01 0.00e+00 3.11e-03  -9.7 1.03e-01  -1.1 1.00e+00 7.93e-02f  1
  97  4.6069624e-01 0.00e+00 3.10e-03 -10.6 1.38e-02  -0.6 1.00e+00 1.47e-01f  1
  98  4.5824023e-01 0.00e+00 3.68e-03  -8.6 3.14e-01  -1.1 8.52e-01 1.81e-01f  1
  99  4.5822224e-01 0.00e+00 3.68e-03  -9.9 2.27e-02  -0.7 1.00e+00 3.30e-03f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  4.5747364e-01 0.00e+00 3.71e-03 -11.0 7.25e-03  -0.3 1.00e+00 3.33e-01f  1
 101  4.5651156e-01 0.00e+00 3.77e-03 -11.0 2.65e-02  -0.8 1.00e+00 1.62e-01f  1
 102  4.5624508e-01 0.00e+00 3.79e-03 -11.0 8.36e-03  -0.3 1.00e+00 1.09e-01f  1
 103  4.5460694e-01 0.00e+00 3.91e-03 -10.5 3.10e-02  -0.8 1.00e+00 2.54e-01f  1
 104  4.5415701e-01 0.00e+00 3.92e-03 -11.0 9.64e-03  -0.4 1.00e+00 1.69e-01f  1
 105  4.5379967e-01 0.00e+00 3.95e-03 -11.0 3.54e-02  -0.9 1.00e+00 5.13e-02f  1
 106  4.5282822e-01 0.00e+00 3.99e-03 -11.0 1.10e-02  -0.4 1.00e+00 3.33e-01f  1
 107  4.5236345e-01 0.00e+00 4.03e-03 -11.0 4.03e-02  -0.9 1.00e+00 6.16e-02f  1
 108  4.5158337e-01 0.00e+00 4.06e-03 -11.0 1.24e-02  -0.5 1.00e+00 2.45e-01f  1
 109  4.4983066e-01 0.00e+00 4.19e-03  -8.9 4.54e-02  -1.0 1.00e+00 2.15e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 110  4.4976166e-01 0.00e+00 4.20e-03 -10.3 1.39e-02  -0.5 1.00e+00 2.01e-02f  1
 111  4.4733643e-01 0.00e+00 4.34e-03  -8.3 4.72e-02  -1.0 9.30e-01 2.78e-01f  1
 112  4.4674741e-01 0.00e+00 4.36e-03  -9.7 1.51e-02  -0.6 1.00e+00 1.60e-01f  1
 113  4.4600496e-01 0.00e+00 4.39e-03  -8.6 4.25e-02  -1.1 1.00e+00 8.06e-02f  1
 114  4.4596146e-01 0.00e+00 4.39e-03 -10.0 1.64e-02  -0.6 1.00e+00 1.09e-02f  1
 115  4.4157951e-01 0.00e+00 4.49e-03  -8.4 4.50e-02  -1.1 1.00e+00 4.49e-01f  1
 116  4.4137682e-01 0.00e+00 4.48e-03  -9.6 1.60e-02  -0.7 1.00e+00 4.90e-02f  1
 117  4.3873388e-01 0.00e+00 4.47e-03  -7.5 3.70e-02  -1.2 9.19e-01 2.63e-01f  1
 118  4.3815153e-01 0.00e+00 4.45e-03  -8.7 1.61e-02  -0.7 1.00e+00 1.33e-01f  1
 119  4.3569743e-01 0.00e+00 4.41e-03  -7.5 3.98e-02  -1.2 1.00e+00 2.34e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 120  4.3444019e-01 0.00e+00 4.36e-03  -8.6 1.67e-02  -0.8 1.00e+00 2.71e-01f  1
 121  4.3145745e-01 0.00e+00 4.27e-03  -7.3 4.24e-02  -1.3 1.00e+00 2.74e-01f  1
 122  4.2975412e-01 0.00e+00 4.19e-03  -8.0 1.81e-02  -0.8 1.00e+00 3.55e-01f  1
 123  4.2433966e-01 0.00e+00 3.98e-03  -6.9 4.56e-02  -1.3 1.00e+00 4.92e-01f  1
 124  4.2414838e-01 0.00e+00 3.96e-03  -7.5 1.92e-02  -0.9 1.00e+00 4.01e-02f  1
 125  4.2194620e-01 0.00e+00 3.87e-03  -6.9 4.80e-02  -1.4 1.00e+00 1.96e-01f  1
 126  4.1778190e-01 0.00e+00 3.63e-03  -7.7 2.08e-02  -0.9 1.00e+00 8.34e-01f  1
 127  4.1737428e-01 0.00e+00 2.94e-03  -7.1 4.72e-02  -1.4 1.00e+00 3.66e-02f  1
 128  4.1315312e-01 0.00e+00 2.54e-03  -7.8 2.22e-02  -1.0 1.00e+00 8.46e-01f  1
 129  4.1173197e-01 0.00e+00 2.51e-03  -7.1 6.31e-02  -1.5 1.00e+00 1.27e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 130  4.0915489e-01 0.00e+00 2.46e-03  -7.7 2.32e-02  -1.0 1.00e+00 5.13e-01f  1
 131  4.0515277e-01 0.00e+00 2.38e-03  -6.9 9.97e-02  -1.5 1.00e+00 3.50e-01f  1
 132  4.0343660e-01 0.00e+00 2.34e-03  -7.3 2.38e-02  -1.1 1.00e+00 3.43e-01f  1
 133  3.9516682e-01 0.00e+00 2.17e-03  -6.8 1.68e-01  -1.6 1.00e+00 7.24e-01f  1
 134  3.9098427e-01 0.00e+00 2.06e-03  -7.2 2.34e-02  -1.1 1.00e+00 9.07e-01f  1
 135  3.9058959e-01 0.00e+00 2.05e-03  -7.0 5.50e-02  -1.6 1.00e+00 3.97e-02f  1
 136  3.8628795e-01 0.00e+00 1.93e-03  -7.7 2.39e-02  -1.2 1.00e+00 9.42e-01f  1
 137  3.8288703e-01 0.00e+00 1.84e-03  -7.2 6.76e-02  -1.7 1.00e+00 3.49e-01f  1
 138  3.8113960e-01 0.00e+00 1.79e-03  -7.8 2.37e-02  -1.2 1.00e+00 4.02e-01f  1
 139  3.7895865e-01 0.00e+00 1.74e-03  -7.3 8.98e-02  -1.7 1.00e+00 2.23e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 140  3.7630588e-01 0.00e+00 1.65e-03  -7.9 2.56e-02  -1.3 1.00e+00 6.04e-01f  1
 141  3.7040736e-01 0.00e+00 2.44e-03  -7.4 4.40e-01  -1.8 1.00e+00 5.52e-01f  1
 142  3.6744301e-01 0.00e+00 2.20e-03  -8.1 5.47e-02  -1.3 1.00e+00 6.87e-01f  1
 143  3.6668338e-01 0.00e+00 2.11e-03  -9.7 1.27e-02  -0.9 1.00e+00 4.33e-01f  1
 144  3.6563026e-01 0.00e+00 2.04e-03 -10.0 1.45e-01  -1.4 1.00e+00 2.33e-01f  1
 145  3.6473497e-01 0.00e+00 1.94e-03 -11.0 1.73e-02  -1.0 1.00e+00 4.85e-01f  1
 146  3.6423323e-01 0.00e+00 2.14e-03 -11.0 2.77e-01  -1.5 1.00e+00 1.00e-01f  1
 147  3.6397151e-01 0.00e+00 2.15e-03 -10.1 2.02e-02  -1.0 1.00e+00 1.32e-01f  1
 148  3.6196968e-01 0.00e+00 2.50e-03  -8.5 1.56e-01  -1.5 1.00e+00 3.92e-01f  1
 149  3.6176499e-01 0.00e+00 2.49e-03  -9.6 1.65e-02  -1.1 1.00e+00 1.00e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 150  3.6069639e-01 0.00e+00 2.45e-03  -8.4 3.49e-02  -1.6 1.00e+00 2.18e-01f  1
 151  3.6041434e-01 0.00e+00 2.43e-03  -9.4 1.63e-02  -1.1 1.00e+00 1.30e-01f  1
 152  3.5921303e-01 0.00e+00 2.37e-03  -8.5 3.78e-02  -1.6 1.00e+00 2.34e-01f  1
 153  3.5869485e-01 0.00e+00 2.33e-03  -9.4 1.64e-02  -1.2 1.00e+00 2.27e-01f  1
 154  3.5699552e-01 0.00e+00 2.25e-03  -7.9 4.06e-02  -1.7 1.00e+00 3.20e-01f  1
 155  3.5674149e-01 0.00e+00 2.23e-03  -8.7 1.77e-02  -1.2 1.00e+00 1.07e-01f  1
 156  3.5471902e-01 0.00e+00 2.11e-03  -7.5 4.33e-02  -1.7 1.00e+00 3.71e-01f  1
 157  3.5335745e-01 0.00e+00 1.99e-03  -8.3 1.91e-02  -1.3 1.00e+00 5.64e-01f  1
 158  3.5245420e-01 0.00e+00 1.94e-03  -7.5 4.53e-02  -1.8 1.00e+00 1.67e-01f  1
 159  3.5021073e-01 0.00e+00 1.73e-03  -8.3 2.05e-02  -1.3 1.00e+00 9.16e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 160  3.4610138e-01 0.00e+00 1.04e-03  -7.3 4.72e-02  -1.8 1.00e+00 7.94e-01f  1
 161  3.4469412e-01 0.00e+00 9.99e-04  -8.0 2.07e-02  -1.4 1.00e+00 6.40e-01f  1
 162  3.4399767e-01 0.00e+00 1.02e-03  -7.7 4.72e-02  -1.9 1.00e+00 1.42e-01f  1
 163  3.4260897e-01 0.00e+00 9.60e-04  -8.1 2.22e-02  -1.4 1.00e+00 6.21e-01f  1
 164  3.4129508e-01 0.00e+00 9.37e-04  -7.4 4.92e-02  -1.9 1.00e+00 2.68e-01f  1
 165  3.4016827e-01 0.00e+00 9.15e-04  -8.0 2.34e-02  -1.5 1.00e+00 5.07e-01f  1
 166  3.3825082e-01 0.00e+00 8.81e-04  -7.2 5.10e-02  -2.0 1.00e+00 3.98e-01f  1
 167  3.3610567e-01 0.00e+00 8.39e-04  -7.7 2.44e-02  -1.5 1.00e+00 1.00e+00f  1
 168  3.3460677e-01 0.00e+00 8.12e-04  -7.1 5.92e-02  -2.0 1.00e+00 3.31e-01f  1
 169  3.3266040e-01 0.00e+00 7.74e-04  -7.7 2.51e-02  -1.6 1.00e+00 9.50e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 170  3.3037125e-01 0.00e+00 7.34e-04  -7.1 8.30e-02  -2.1 1.00e+00 5.31e-01f  1
 171  3.2873890e-01 0.00e+00 7.02e-04  -7.7 2.53e-02  -1.6 1.00e+00 8.56e-01f  1
 172  3.2622239e-01 0.00e+00 6.57e-04  -7.2 5.31e-02  -2.1 1.00e+00 6.32e-01f  1
 173  3.2470308e-01 0.00e+00 6.28e-04  -7.8 2.53e-02  -1.7 1.00e+00 8.69e-01f  1
 174  3.2372187e-01 0.00e+00 6.27e-04  -7.5 5.26e-02  -2.2 1.00e+00 2.63e-01f  1
 175  3.2283458e-01 0.00e+00 5.92e-04  -8.2 2.60e-02  -1.7 1.00e+00 5.20e-01f  1
 176  3.2183252e-01 0.00e+00 6.17e-04  -7.5 5.47e-02  -2.2 1.00e+00 2.71e-01f  1
 177  3.2046827e-01 0.00e+00 5.47e-04  -8.1 2.70e-02  -1.8 1.00e+00 8.08e-01f  1
 178  3.1975003e-01 0.00e+00 6.03e-04  -7.1 5.85e-02  -2.3 1.00e+00 2.04e-01f  1
 179  3.1900035e-01 0.00e+00 5.19e-04  -7.6 2.77e-02  -1.8 1.00e+00 4.53e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 180  3.1615561e-01 0.00e+00 4.90e-04  -6.8 6.32e-02  -2.3 1.00e+00 8.53e-01f  1
 181  3.1524761e-01 0.00e+00 4.49e-04  -7.3 2.68e-02  -1.9 1.00e+00 5.85e-01f  1
 182  3.1326122e-01 0.00e+00 4.87e-04  -7.0 6.24e-02  -2.4 9.98e-01 6.26e-01f  1
 183  3.1296679e-01 0.00e+00 4.09e-04  -7.7 2.63e-02  -1.9 1.00e+00 2.01e-01f  1
 184  3.1055716e-01 0.00e+00 4.34e-04  -6.8 6.57e-02  -2.4 1.00e+00 8.22e-01f  1
 185  3.1020674e-01 0.00e+00 3.60e-04  -7.5 2.70e-02  -2.0 1.00e+00 2.54e-01f  1
 186  3.0929569e-01 0.00e+00 4.68e-04  -7.0 6.86e-02  -2.5 1.00e+00 3.18e-01f  1
 187  3.0819030e-01 0.00e+00 3.23e-04  -7.8 2.94e-02  -2.0 1.00e+00 8.02e-01f  1
 188  3.0763044e-01 0.00e+00 4.47e-04  -7.2 8.21e-02  -2.5 1.00e+00 2.02e-01f  1
 189  3.0725472e-01 0.00e+00 3.07e-04  -8.0 3.11e-02  -2.1 1.00e+00 2.85e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 190  3.0578606e-01 0.00e+00 4.05e-04  -7.5 2.33e-01  -2.6 1.00e+00 5.11e-01f  1
 191  3.0495222e-01 0.00e+00 3.43e-04  -8.1 2.91e-02  -2.2 1.00e+00 6.59e-01f  1
 192  3.0413633e-01 0.00e+00 1.17e-03  -7.6 4.93e-01  -2.6 1.00e+00 2.84e-01f  1
 193  3.0331078e-01 0.00e+00 9.60e-04  -7.7 3.14e-02  -2.2 1.00e+00 6.81e-01f  1
 194  3.0239991e-01 0.00e+00 8.65e-04  -6.7 7.67e-02  -2.7 9.99e-01 4.00e-01f  1
 195  3.0137954e-01 0.00e+00 6.77e-04  -7.5 3.36e-02  -2.3 1.00e+00 8.52e-01f  1
 196  3.0056119e-01 0.00e+00 6.24e-04  -7.5 8.05e-02  -2.7 1.00e+00 3.28e-01f  1
 197  3.0014178e-01 0.00e+00 5.65e-04  -8.3 3.58e-02  -2.3 1.00e+00 3.62e-01f  1
 198  2.9923633e-01 0.00e+00 5.17e-04  -8.2 8.53e-02  -2.8 1.00e+00 3.59e-01f  1
 199  2.9888005e-01 0.00e+00 4.75e-04  -8.9 3.83e-02  -2.4 1.00e+00 3.15e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 200  2.9848396e-01 0.00e+00 4.57e-04  -8.1 8.98e-02  -2.8 1.00e+00 1.58e-01f  1
 201  2.9821853e-01 0.00e+00 4.30e-04  -8.8 4.15e-02  -2.4 1.00e+00 2.29e-01f  1
 202  2.9738245e-01 0.00e+00 3.96e-04  -7.2 9.62e-02  -2.9 9.99e-01 3.39e-01f  1
 203  2.9665614e-01 0.00e+00 3.36e-04  -7.8 4.41e-02  -2.5 1.00e+00 6.30e-01f  1
 204  2.9578199e-01 0.00e+00 3.07e-04  -6.8 1.02e-01  -2.9 1.00e+00 3.87e-01f  1
 205  2.9465429e-01 0.00e+00 2.39e-04  -7.4 4.54e-02  -2.5 1.00e+00 1.00e+00f  1
 206  2.9301415e-01 0.00e+00 1.99e-04  -6.7 1.07e-01  -3.0 1.00e+00 8.01e-01f  1
 207  2.9292261e-01 0.00e+00 1.94e-04  -7.3 4.47e-02  -2.6 1.00e+00 8.17e-02f  1
 208  2.9132376e-01 0.00e+00 2.05e-04  -6.8 1.11e-01  -3.0 1.00e+00 7.08e-01f  1
 209  2.9068030e-01 0.00e+00 1.46e-04  -7.4 4.72e-02  -2.6 1.00e+00 5.83e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 210  2.8989341e-01 0.00e+00 2.15e-04  -6.8 1.19e-01  -3.1 1.00e+00 3.71e-01f  1
 211  2.8905846e-01 0.00e+00 1.25e-04  -7.4 5.00e-02  -2.7 1.00e+00 7.99e-01f  1
 212  2.8817010e-01 0.00e+00 2.35e-04  -6.9 1.32e-01  -3.1 1.00e+00 4.36e-01f  1
 213  2.8741077e-01 0.00e+00 1.25e-04  -7.5 5.21e-02  -2.7 1.00e+00 7.41e-01f  1
 214  2.8705020e-01 0.00e+00 2.96e-04  -7.2 1.57e-01  -3.2 1.00e+00 1.65e-01f  1
 215  2.8639044e-01 0.00e+00 1.29e-04  -8.0 5.54e-02  -2.8 1.00e+00 6.44e-01f  1
 216  2.8567025e-01 0.00e+00 2.13e-04  -7.4 4.30e-01  -3.2 1.00e+00 3.16e-01f  1
 217  2.8519078e-01 0.00e+00 1.22e-04  -8.1 6.05e-02  -2.8 1.00e+00 4.79e-01f  1
 218  2.8508336e-01 0.00e+00 1.93e-04  -7.5 3.80e-01  -3.3 1.00e+00 4.58e-02f  1
 219  2.8459190e-01 0.00e+00 1.16e-04  -8.3 6.70e-02  -2.9 1.00e+00 4.75e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 220  2.8382163e-01 0.00e+00 1.63e-04  -7.2 3.30e-01  -3.3 1.00e+00 3.29e-01f  1
 221  2.8341130e-01 0.00e+00 1.28e-04  -7.8 7.27e-02  -2.9 1.00e+00 3.93e-01f  1
 222  2.8284739e-01 0.00e+00 1.56e-04  -7.0 2.21e-01  -3.4 1.00e+00 2.45e-01f  1
 223  2.8183563e-01 0.00e+00 1.19e-04  -7.7 7.96e-02  -3.0 1.00e+00 9.45e-01f  1
 224  2.8117611e-01 0.00e+00 1.40e-04  -6.9 4.77e-01  -3.4 1.00e+00 2.96e-01f  1
 225  2.8033847e-01 0.00e+00 1.04e-04  -7.6 8.59e-02  -3.0 1.00e+00 7.96e-01f  1
 226  2.8005033e-01 0.00e+00 9.73e-05  -9.1 3.26e-02  -2.6 1.00e+00 6.19e-01f  1

Number of Iterations....: 226

                                   (scaled)                 (unscaled)
Objective...............:   2.8005032523583379e-01    2.8005032523583379e-01
Dual infeasibility......:   9.7272674831558608e-05    9.7272674831558608e-05
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   3.2315149852880829e-07    3.2315149852880829e-07
Overall NLP error.......:   9.7272674831558608e-05    9.7272674831558608e-05

Number of objective function evaluations             = 227
Number of objective gradient evaluations             = 227
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 228
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations             = 226
Total CPU secs in IPOPT (w/o function evaluations)   =   3403.326
Total CPU secs in NLP function evaluations           =    214.977

EXIT: Optimal Solution Found.
```

.... the timing is ok at this resolution, but at higher ones, the linear solves become quite large and extremely slow. So I guess my question is: can I expect something like this (at a higher resolution, with ~500 000 parameters) to take a couple of days to solve or are there any settings I could tweak to increase convergence / reduce the number if iterations needed. The initial guesses are quite good actually btw, and a local optimum is actually all I'm looking for.


问:
Unstable convergence of a Poisson equation
What could be the reason that the solution of a Poisson equation is smooth when obtained by an iterative solver, only if the maximum residual is set to a high value (e.g. 0.1)? When the maximum residual is set as a typical lower value (e.g. 1e-6), then the solution oscillates - often (but not exclusively) near Neumann boundary conditions. The solution with the high residual is underpredicted, and with the low residual it oscillates around expected values.

The matrix is non-symmetric sparse, weakly diagonally-dominant (a diagonal entry equals the negative sum of other row entries). The system includes both Dirichlet and Neumann boundary conditions. Similar outcomes happen using various iterative solvers: BiCGStab, CGS, QMR and GMRES. The preconditioner is always a diagonal preconditioner, since I have to use matrix-free solving approach on the GPU. Using no preconditioner yields similar (though a bit smoother) results.


问:
Finite difference Neumann boundary conditions: uneven weighting of edge nodes?
Originally asked this on math.stackexchange, but I figure it's also appropriate here. I'm reading through some finite difference code for a diffusion equation and came across something odd for the boundary conditions, and I was wondering about the validity of the method. I'm familiar with the finite difference method, but I've never seen this certain method before.

The equation I am solving is

$$ \frac{\partial f}{\partial t}=C(x)\frac{\partial^2 f}{\partial x^2}$$

with boundary condition

$$ \frac{\partial f}{\partial x}=0\ \text{on}\ x=0,1 $$

across time, $t$, and position, $x$, with position-dependent diffusion constant, $C(x)$. It's being solved in an explicit way using the forward difference for time and the second-order central difference for space

$$ \frac{f^\text{new}_i-f^\text{old}_i}{\Delta t}=C_i\frac{f^\text{old}_{i-1}-2f^\text{old}_{i}+f^\text{old}_{i+1}}{(\Delta x)^2}$$

along discretised space $\{x_1,\dots,x_n\}$. Rearranging the above equation gives

\begin{align}
f^\text{new}_i=&f^\text{old}_i + \frac{C_i\Delta t}{(\Delta x)^2}\Big(f^\text{old}_{i-1}-2f^\text{old}_{i}+f^\text{old}_{i+1}\Big)\\
=&f^\text{old}_{i-1}(S_i)+f^\text{old}_{i}(1-2S_i)+f^\text{old}_{i+1}(S_i)
\end{align}

where $S_i=\frac{C_i\Delta t}{(\Delta x)^2}$.

This is simply represented in the code. For the boundary at $x_1$ (and at $x_n$, but I'll only talk about the first as they are similar), there are some strange things going on. Usually I would do a Taylor series expansion of function $f$ at points $x_2$ and $x_3$

$$ f_2=f_1+\frac{\partial f}{\partial x}\bigg\rvert_i \Delta x + \frac{1}{2} \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i (\Delta x)^2 + \mathcal{O}((\Delta x)^3)$$
$$ f_3=f_1+\frac{\partial f}{\partial x}\bigg\rvert_i (2\Delta x) + \frac{1}{2} \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i (2\Delta x)^2 + \mathcal{O}((\Delta x)^3)$$

which, due to boundary condition $\frac{\partial f}{\partial x}=0$ can be rearranged to

$$ \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{2}{(\Delta x)^2}\Big(f_2 - f_1\Big)\tag{$A$}$$
$$ \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{2(\Delta x)^2}\Big(f_3 - f_1\Big)\tag{$B$}.$$

Usually, I'd average $A$ and $B$ to determine how the boundary conditions are applied in the eventual matrix solution

$$ \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{(\Delta x)^2}\Big(-\frac{5}{4}f_1 + f_2 + \frac{1}{4}f_3 \Big).$$

However, the code that I'm reading does something different and takes the boundary as a weighted sum of $A/5+4B/5$, as shown below

$$\frac{A+4B}{5} \Rightarrow \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{(\Delta x)^2}\Big(-\frac{4}{5}f_1 + \frac{2}{5}f_2 + \frac{2}{5}f_3 \Big).$$

I don't think this as being inherently wrong, as the LHS is still $\frac{\partial^2 f}{\partial x^2}\big\rvert_i$, but I'm not sure. Is this a valid step to make? Is it ok to "weight" results from different nodes, and if so, in what situations would weighting or not weighting be useful?


问:
Inconsistency in optimize.minimize
I am trying to fit a time-dependent curve at each time step. I do so in minimizing along $x_c$ the quadratic error between the curve and a reference solution $ 1/(1 + \exp\left(\sqrt{S}(x-x_c)\right) $:

$$ \frac{1}{N+2} \sum_{i=0}^{N+1} \left(y_i - \frac{1}{1 + \exp\left(\sqrt{S}(x_i-x_c)\right)}\right)^2 $$

Because the calculations are easy, I calculated the gradient so as to use a conjugate gradient method implemented in optimize.minimize in scipy.

So: fairly easy stuff.

Below is my MWE and here is my data file:

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

def energie(centre, x, y, S):
    erreur = 0.
    for i in range(0, y.size):
        erreur += (y[i] - 1./(1.+np.exp(np.sqrt(S)*(x[i]-centre))))**2.
    return erreur/y.size

def gradenergie(centre, x, y, S):
    grad = 0.
    for i in range(0, y.size):
        grad += - (y[i] - 1./(1.+np.exp(np.sqrt(S)*(x[i]-centre)))) * np.sqrt(S) * \
        np.exp(np.sqrt(S)*(x[i]-centre)) * (1.+np.exp(np.sqrt(S)*(x[i]-centre)))**-2.
    return 2.*grad/y.size

T = np.arange(0., 36., 0.5/np.sqrt(2))
x = np.arange(301) * 0.5
q = np.load('data/data.npy')
centres = np.zeros(q.shape[0])

for n in range(q.shape[0]):
    # initial guess
    centreq = centres[n-1]
    # call to optimize.minimize
    resultatq = optimize.minimize(energie, centreq, method = 'CG', jac = gradenergie, args = (x, q[n , : ], 0.1))
    # write to centres
    centres[n] = resultatq.x

dcentresqdt = np.gradient(centres, T[1]-T[0])
plt.plot(dcentresqdt[2 : ])
plt.show()
</code>```

which produces the following graph:

with obvious irregularities. Any idea why?

--

Edit:

in the code centre or centreq stands for $x_c$, each q[n, :] stands for the $(y_i)$, and the $(x_i)$ are the same at each time step, they are in the array x.

Indeed I plot the derivative of centres with respect to time, but if I plot centres itself I have the same kind of problem.

--

Edit2:

I realized my problem is so standard that there is actually a scipy function that does the job in no time (101 seconds for Nelder Mead vs 0.39 for the version below):

```<code>from scipy.optimize import curve_fit

def kink(x, centre):
    S = 0.1
    return 1./(1.+np.exp(np.sqrt(S)*(x-centre)))
</code>```

and later on in the code:

```<code>for n in range(q.shape[0]):
    centreq = centres[n-1]
    popt, pcov = curve_fit(kink, x, q[n , : ], p0 = centreq)
    centres[n] = popt
</code>```

Title: Inconsistency in optimize.minimize


问:
Weighted Set Cover in practice, beyond the greedy algorithm
According to the wikipedia page for Set Cover, the greedy algorithm for weighted set cover achieves the polynomial-time approximation bound. There are other techniques for solving Set Cover, such as randomized rounding of a linear programming relaxation, which I believe achieve the same asymptotic bound.

Is there information comparing these asymptotically-equivalently-performing algorithms in practice?


问:
conjugate gradient for Newton's method with non positive definite Hessian matrix
I want to minimize a non-linear function $f(x)$ using Newton's method. At each optimization step, I compute a descent direction $d$ to update $x$ using a second-order approximation of $f(x)$:

$$ \nabla^2f(x) \ d = -\nabla f(x) $$

So I need to solve a linear system, but the Hessian matrix $\nabla^2f(x)$ might not be positive definite which could give me an invalid descent direction. I want to run conjugate gradient to solve the linear system. In the conjugate gradient method, one of the steps involves this operation $v^T \nabla^2\ f(x)\ v$ (computing $\alpha$ according to Shewchuk's notes, page 50 https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf , I'm using $v$ rather than $d$ because in my convention $d$ is the descent direction to update $x$). I could recycle this operation to know if the Hessian is not positive definite (if such operation is negative).

I implemented my algorithm like that, so as soon as I detect that operation is negative, I stop the CG solver and return the solution iterated up to that point. I've actually seen it works pretty well in practice, but I have no rigorous justification for doing it. I remember hearing / reading people saying that each CG iteration takes care of different eigenvalues in order of "importance", so by stopping early you end up exploring the positive definite part of the Hessian. Can anybody point me to a more rigorous justification of this? Is this a valid way to enforce positive definiteness using CG for non-linear optimization?


问:
Understanding design patterns
I read some books and scoured the web about design pattern,but almost all ways of expressing design patterns are the same. They define what is it,draw UML and give an example,but what I want is to suggest me any site or book that writes some bad code and refactor using one of the design patterns?</p>

Title: Understanding design patterns


问:
What Derivative-free optimization method should I use when my initial guess is very good?
I am trying to minimize a function where my initial guess is quite close to the minimum.
I'm trying to minimize

$$f(q) = \text{angle}(qw_1q*, v_1) + \text{angle}(qw_2q*, v_2) + \text{angle}(qw_3q*, v_3)$$

with $q$ being a quaternion and it's conjugate $q*$ and $w$ and $v$ being vectors. What algorithms are there that can find my minimum without requiring expensive derivatives? Speed is my biggest concern.


问:
Basic approach for numerical solution of PDE
I'm looking for some guidance on how to write a program to numerically solve a PDE. As an example for comparison in 1D:

$$\frac{d^2u}{dx^2} = f\;\;\;\;u(0) = 0\;\;\;\;u(1) = 0$$

We could try 6 equally spaced nodes between 0 and 1, interpolate within each element with a quadratic, and ensure continuity of $u$ and $du/dx$ at the boundaries, and to solve within each element we can integrate the differential equation over each element. I believe this constant-weight integration makes this FVM. Is that correct?

In any case, there were 3 DOFs for each of the 5 elements, 1 integral equation per element, 2 compatibility conditions for each of the 4 internal boundaries, and 2 boundary conditions given as part of the problem. So that's 15 equations in 15 variables which can then be solved. When trying to make the equivalent 2D method, I cannot work out how to generate the right number of compatibility equations. For example:

$$\nabla^2{u} = f\;\;\;\;u(0,y) = u(1,y) = u(x,0) = u(x,1) = 0$$

Imagine we mesh with a grid of 5 by 5 squares with each square divided into two triangles. Each triangular element can be represented with a quadratic in $x$ and $y$ (6 DOFs). The DE is integrated over for each element. Continuity of $u$ and $\nabla{u}$ could be ensured at the boundaries, though already this might be over/under-constrained.

There are 50 elements, 65 internal boundaries and 20 external boundaries. There are 300 DOFs. Compatibility of $u$ and $\nabla{u}$ on a single internal boundary yields 5 equations, but either the sets of 5 equations per internal boundary are not independent of each other or the system is over-constrained because already there are too many equations.

As far as I understand, representing the DOFs with the value of $u$ at the corners and edge midpoints of the triangles is a way of intrinsically satisfying the continuity of $u$ condition and reducing the number of DOFs, but it does not help intrinsically satisfy the other compatibility conditions. But I mention this in case proper representation of the DOFs is key to helping understand the problem.

If I am thinking about this in roughly the right way, how do I determine the compatibility conditions and generate the right number of equations to solve for the given number of DOFs? If I'm not thinking about the problem in the right way but you can help me funnel this into a more specific, generally-relevant question(s), please do.


问:
Simultaneous update to barycenters
Suppose a tiling is given in 2D (an embedding of a planar triangulated graph), with all faces convex.

Now suppose one moves each point, one by one, to the barycenter of its neighbors. I think that with such an update one can guarantee that the planarity is preserved at each step (after each point is moved). But could planarity be guaranteed in case whole configuration moves simultaneously?

I understood that once such simultaneous position update converges (in case the positions of the points of the outer convex face are fixed), it can be guaranteed that the new tiling is also planar (Tutte's embedding). My question is: given no fixed positions, what happens if one performs a single simultaneous coordinate update from the initial configuration that is already planar (all faces convex)?

In case no guarantees can be made, are there special cases with which the simultaneous update would preserve planarity, i.e., special kinds of initial configurations, etc?

In mathematical terms: The initial configuration is $Y_0$, and can be written as $Y_0 = T_0Y_0$, where $T_0$ is a transition matrix (each row/column having weight coefficients expressing each coordinate as a convex combination of its neighbors). The update is achieved as $Y_1 = T_bY_0$, where $T_b$ is a transition matrix that moves a configuration to the barycenter (a simultaneous update). Obviously, $Y_1 = T'Y_0$, where $T'$ is also a transition matrix $T' = T_bT_0$. But is that sufficient to guarantee that $Y_1 = T''Y_1$, for some transition matrix $T''$, i.e., each coordinate is a convex combination of its neighbors?

Title: Simultaneous update to barycenters


问:
Should I expect computational gains using a second-order splitting method here?
I am trying to solve a three-dimensional baroclinic transport problem. The hydrodynamic (three-dimensional shallow water) equations are:
\begin{align}
\nabla \cdot \vec{v} = 0, \tag{1} \\
\frac{\partial u}{\partial t} + \nabla \cdot (u \vec{v} ) + \frac{1}{\rho_0} \frac{\partial P}{\partial x} - f_{\text{coriolis}}v - \frac{1}{\rho_0}(\nabla \cdot T_{\text{mx}}) - S_x = 0, \tag{2} \\
\frac{\partial v}{\partial t} + \nabla \cdot (v \vec{v} ) + \frac{1}{\rho_0} \frac{\partial P}{\partial y} + f_{\text{coriolis}}u - \frac{1}{\rho_0}(\nabla \cdot T_{\text{my}}) - S_y = 0, \tag{3} \\
\frac{\partial H}{\partial t} + \frac{\partial (\bar{u}H)}{\partial x} + \frac{\partial (\bar{v}H)}{\partial y} = 0 \tag{4}.
\end{align}
The constituent (salinity or temperature) transport equation is:
\begin{align*}
\frac{\partial c}{\partial t} = \nabla \cdot (\vec{v} c) + \nabla \cdot (\bar{\bar{D}} \nabla c) = 0 \tag{5}.
\end{align*}
The problem is to solve for water velocity $\vec{v} = [u,v,w]^T$, river depth $H$, and constituent concentrate $c$ over a fluid domain $\Omega \subset \mathbb{R}^3$ and timescale $[0, t_{\text{end}}]$ of interest. The first equation is usually called the incompressible continuity equation. The second and third equations are called the horizontal momentum equations. The fourth equation is called the depth-averaged continuity equation. The fifth equation is the transport equation.

The other important relationship to note is that the fluid pressure $P$ (which appears in the horizontal momentum equations) is related to the concentration $c$. Specifically, $P$ depends on the fluid density $\rho$, and $\rho$ is related to $c$ through an algebraic equation of state. So there is two-way coupling between the hydrodynamic equations and the transport equation.

The solution procedure used to run the simulation one time-step is as follows. Note that spatial discretization is done using a first-order SUPG FEM, and temporal discretization is done using a second-order (implicit) BDF method.

Using $u,v,w,H,P$ from the old time-step, we solve discrete weak version of the horizontal momentum equations (2)-(3) and  the depth average continuity (4) monolithically to get $u, v, H$ at the new time-step. Since the horizontal momentum equations (2)-(3) are nonlinear, Newton's method with line-search is used.

Using the new $u, v, H$ from step 1, a discrete weak version of incompressible continuity equation (1) is solved to get $w$ at the new time-step. This is a linear problem, so this is computationally easy compared to step 1.

Using the new $u, v, w, H$ from steps 1 and 2, a discrete weak version of the transport equation (5) is solved to get $c$ at the new time-step. This is also linear, so this is computationally easy compared to step 1. Also, this new $c$ value is used in the calculation of a new $P$ value for use in step 1.

After these three steps, I have all the variables we are looking for at the next discrete time value for all spatial coordinates in our grid. So I repeat these three steps to march in time.

This solution procedure might be seen as a first-order operator splitting schemes, since the hydrodynamics processes (equations (1)-(4) and steps 1-2) are solved separately from the transport process (equation (5) and step 3). Should one expect to gain any sort of advantage (e.g., in computational efficiency) by implementing  second-order splitting schemes (e.g., from this webpage: http://www.asc.tuwien.ac.at/~winfried/splitting/)?

I have ran some computational experiments testing a few second-order schemes where the hydrodynamic equations (1)-(4) are split from equation (5), i.e., the coefficients $a_i$ correspond fractional step coefficients for hydrodynamic processes while the coefficients $b_i$ correspond to fractional step coefficients for transport processes; however, I do not seem to see any computational gains at the moment.

The experiments are run on a lock-exchange problem that is similar to the one found in [Shin et al. 2014]. Initially, the left-half of the flume (with length 2 meters, width 0.2 meters, and depth 0.2 meters) is filled with saltwater, the right-half is filled with freshwater, so there is a sharp interface in the middle at $t = 0$. The barrier is instantly removed, so the dense saltwater flows beneath the less-dense freshwater. The fixed uniform mesh that I use to run my experiments is relatively coarse (grid spacing is 0.2 meters), but I don't think this matter since I am only interested in temporal convergence. I roughly measure computationally efficiency by plotting the discrete $L_2$ error versus the wall-clock time of the simulation. Below is an example of one such plot (split1 through split5 correspond to different implementations of second-order splitting schemes, while split0 corresponds to first-order splitting). Note that the "reference solution" used in this discrete $L_2$ error measurement is actually an approximate solution to the semi-discrete system, which is computed using a very small time-step size of 0.01 seconds. Also, I am specifically looking at the error in salt concentration $c$.

Looking at points which are vertically aligned, you can see that the black points (first-order splitting) achieves an $L_2$ error which is roughly the same as the $L_2$ errors computed using second-order splitting (in fact, they are a little better).

References:

Shin, J. O., S. B. Dalzieland, and P. F. Linden. 2004. Gravity currents produced by lock exchange. Journal of Fluid Mechanics 521:1–34.


问:
Cressman interpolation and objective analysis
I have read this question and answer – Interpolation of scattered data to a regular grid in python and I am doing something similar as I have temperature values of the atmosphere at different heights defined on a  coarse grid  and I want to interpolate this to a regular grid in Python and I also need to weight the matches using the Cressman weighting scheme as shown in this old paper – Cressman weighting scheme.

In this regard, I have a couple of questions. Looking at equation (1)

$$C_h =  -W  E_h$$

I understand what $W$ is and how to go about calculating it but it is not clear to me what the author means when he writes:

$E_h$ is the error of the interpolated value of the first guess field at the location of the observation

An error is usually a difference between two values. If that is the case which two values are being referred to here? As a first guess, I have my temperature values at irregular grid intervals and I can easily find matches within a radius of influence for a particular grid point. But which two values are contributing to the error that is being talked about here?

I have also read papers (cf. Vertical and horizontal Cressman radii) in which it is mentioned that the radius of influence that Cressman talks about in the above paper sometimes can be treated as a "deep layer" in the atmospheric sciences. In other words, the interpolated value at a grid point can include matches from not only the horizontal radius of influence on a given horizontal surface but also matches from above and below the surface in question. Is that accurate?

OwnerDisplayName: gansub


问:
How to analyze the dispersion and dissipation of a certain FEM?
In Finite Difference method or Finite volume textbook. Dispersion/Dissipation can be analyzed by set $u = u_0\exp{\omega t +\mathbf{kx}}$。

However, I cannot find something about this kind of analysis in FEM without convert it to a system of difference equations. And FEM is a little more complex because of there are two function spaces and some time they are not continuous.

Is there example to do it in the framework of FEM? I saw someone used the Eigen value distribution of assembled matrix.


问:
Solving a PDE implicitly by iteration in python
Connected to this question here on Computational Science, I've posted a follow-up question on how to solve a PDE using an implicit scheme like Crank-Nicholson in general in this question on SO.

But I guess it is better to ask here it is correct to solve a PDE with an implicit step in this way in general and ask the more specific coding issues on SO. So I'll try to split up the question I asked on SO.

So I've got a non-linear advection-diffusion(-reaction) PDE which I want to solve for one single timestep using Crank-Nicolson. In this example my 1D-PDE is the following:

$$\frac{T^{n+1} - T^n}{\Delta t} \approx \frac{1}{2} \left(\left(\alpha(T^n) \frac{T^n_{i-1} - 2T^n_{i} + T^n_{i+1}}{\Delta x^2} + v(T^n) \frac{T^n_{i\pm1} - T^n_{i}}{\Delta x}\right) + \left( \alpha(T^{n+1}) \frac{T^{n+1}_{i-1} - 2T^{n+1}_{i} + T^{n+1}_{i+1}}{\Delta x^2} + v(T^{n+1}) \frac{T^{n+1}_{i\pm1} - T^{n+1}_{i}}{\Delta x}\right)\right)$$

Where $\alpha$ is the temperature dependent diffusion coefficient and $v$ is the temperature dependent advection speed. The advection is discretized using an upwind scheme, since when there is advection $|Pe|\geq2$ can be assumed.

My current program code returns the derivatives of the PDE for given input temperatures, summing up the specific advection and diffusion derivatives. Thus my solver can't "see" the diffusion coefficient and the advection speed. Using a matrix notation to solve the PDE is thus not applicable.

Furthermore the derivatives and value (temperature) arrays have the same shape as the calculation grid, meaning that for the simplest case of the 1D PDE with, for example 10 nodes, the derivatives and value arrays are vectors with 10 cells (of course neglecting the additional ghost cells I use for fast stencil solving and boundary conditions).

The derivative terms in my PDE are thus in my program code replaced with the following vector derivative where each element in $derivative$ corresponds to the derivative of its respective temperature element in the value array:
$$derivative(T^n) = \alpha(T^n) \frac{T^n_{i-1} - 2T^n_{i} + T^n_{i+1}}{\Delta x^2} + v(T^n) \frac{T^n_{i\pm1} - T^n_{i}}{\Delta x}
$$

So in pseudo-code equation my PDE looks like this:
$$\frac{T^{n+1} - T^n}{\Delta t} \approx \frac{1}{2} \left(derivative(T^{n}) + derivative(T^{n+1})\right)$$

To solve this, I reorganize this equation so that
$$T^{n+1} - T^n - \frac{\Delta t}{2} \left(derivative\left(T^{n}\right) + derivative(T^{n+1})\right) = 0$$

which in python looks like

```<code>def crank_nicolson(y, yprev, h):
    return (y - yprev - h / 2 * (diff(y) + diff(yprev)))
</code>```

with the stepsize $h = \Delta t$, $y=T^{n+1}$, $yprev=T^{n}$ and $diff = derivative(T)$.

This can easily be fed to any root finding algorithm like the scipy.optimize.root module or using any other self-written iterative solver for example like:

```<code>from scipy import optimize

def diff(T):  # simply example differential function
    return np.sqrt(T)

def iterate_cranknicolson(y0, h, rtol, diff):  # simple iterative solver
    y = y0
    diff_old = diff(y)  # save forward diff
    f = np.ones_like(y)
    i = 0
    while np.any(f > rtol):
        y_old = y
        y = y0 + h / 2 * (diff_old + diff(y))  # at the first step: euler forward predictor
        f = (y / y_old - 1)  # simple stop criterion
        i += 1
    return y

# some test runs
T0 = np.arange(10, 20).astype(float)  # set starting temperature

cn_hybr = optimize.root(crank_nicolson, T0, args=(T0, 1.), method='hybr', tol=1e-6)
cn_it = iterate_cranknicolson(T0, 1., 1e-6, diff)
</code>```

TL,DR. The question(s) summmarized:

Is this approach to iteratively solve an implicit discretization of a PDE correct? With "this approach" I mean especially reorganizing the PDE to $T^{n+1} - T^n - \frac{\Delta t}{2} \left(derivative\left(T^{n}\right) + derivative(T^{n+1})\right) = 0$ and finding the root of this equation.

Is there any better way to do it without being able to use matrix notation?

Is the derivative of $T^n$ the first column of the Jacobian matrix, so that I could use a modified Newtonian iteration, as proposed here and does this work without having a matrix notation at all?

(not really important, thus secondary: Is it ok to use upwind scheme for advection in Crank-Nicolson or should I stick to central differences?)

Thanks alot in advance! Please tell me if I should refine parts of my question and/or add more detail.


问:
Derive total energy balance equation from Chapman-Enskog analysis of lattice Boltzmann equation
I'm interested to derive the total energy balance from Chapman-Enskog analysis of lattice Boltzmann equation (LBE). I know, I should go to the second moment of LBE (zeroth moment gives mass conservation and first moment gives momentum conservation). In fact, I derived the main equation for the second moment at convective ($t_{1}$) and diffusive ($t_{2}$) time scales as:

From Chapman-Enskog analysis up to the second order:

Convective time scale: $\frac{\partial f_{i}^{eq}}{\partial t_{1}} + \vec{e_{i}} \cdot \nabla_{1} f_{i}^{eq} = -\frac{f_{i}^{(1)}}{\tau}$

Diffusive time scale: $\frac{\partial f_{i}^{eq}}{\partial t_{2}} + (1-\frac{1}{2 \tau}) \Big [ \frac{\partial f_{i}^{(1)}}{\partial t_{1}} + \vec{e_{i}} \cdot \nabla_{1} f_{i}^{(1)} \Big ] = -\frac{f_{i}^{(2)}}{\tau}$

Where $f_{i}^{eq}$ is the equilibrium distribution function (expanded Maxwell-Boltzmann up to third order):

$f_{i}^{eq} = \rho \omega_{i} (1 + \frac{\vec{e_{i}} \cdot \mathbf{u}}{c_{s}^{2}} + \frac{(\vec{e_{i}} \otimes \vec{e_{i}} - \mathbf{I})\cdot \mathbf{u} \otimes \mathbf{u}}{2 c_{s}^{4}} + \frac{(\vec{e_{i}} \otimes \vec{e_{i}} \otimes \vec{e_{i}} - 3 \vec{e_{i}} \otimes \mathbf{I}) \cdot \mathbf{u} \otimes \mathbf{u} \otimes \mathbf{u}}{6 c_{s}^{6}})$

$f_{i}^{(1)}$ and $f_{i}^{(2)}$ are the first and second order terms of distribution function Taylor's expansion.

We know the total energy and its flux could be derived as:

$\rho e = \rho (\frac{3}{2} U^{2} + \frac{1}{2} |\mathbf{u}|^{2}) = \sum_{i} \frac{1}{2} \vec{e_{i}} \cdot \vec{e_{i}} f_{i}^{eq}$

$\rho e \mathbf{u} = \rho (\frac{3}{2} U^{2} + \frac{1}{2} |\mathbf{u}|^{2}) \mathbf{u} = \sum_{i} \frac{1}{2} \vec{e_{i}} \cdot \vec{e_{i}} \vec{e_{i}} f_{i}^{eq}$

Finally, we have:

Convective time scale :$\frac{\partial}{\partial t_{1}} (\rho e) + \nabla_{1} \cdot (\rho e \mathbf{u}) = -\frac{1}{\tau} \sum_{i} \frac{1}{2} f_{i}^{(1)} \vec{e_{i}} \cdot \vec{e_{i}}$

Diffusive time scale: $\frac{\partial}{\partial t_{2}} (\rho e) + (1-\frac{1}{2 \tau}) \frac{\partial}{\partial t_{1}} (\sum_{i} \frac{1}{2} \vec{e_{i}} \cdot \vec{e_{i}} f_{i}^{(1)}) + (1-\frac{1}{2 \tau}) \nabla_{1} \cdot (\sum_{i} \frac{1}{2} \vec{e_{i}} \cdot \vec{e_{i}} \vec{e_{i}} f_{i}^{(1)}) = -\frac{1}{\tau} \sum_{i} \frac{1}{2} f_{i}^{(2)} \vec{e_{i}} \cdot \vec{e_{i}}$

The terms $\sum_{i} \frac{1}{2} f_{i}^{(1)} \vec{e_{i}} \cdot \vec{e_{i}}$ and $\sum_{i} \frac{1}{2} f_{i}^{(2)} \vec{e_{i}} \cdot \vec{e_{i}}$ are equal to zero because energy is a conserved quantity.

As a result:

Convective time scale :$\frac{\partial}{\partial t_{1}} (\rho e) + \nabla_{1} \cdot (\rho e \mathbf{u}) = 0$

Diffusive time scale: $\frac{\partial}{\partial t_{2}} (\rho e) + (1-\frac{1}{2 \tau}) \nabla_{1} \cdot (\sum_{i} \frac{1}{2} \vec{e_{i}} \cdot \vec{e_{i}} \vec{e_{i}} f_{i}^{(1)}) = 0$

Right now, my question is how can I simplify the term $(1-\frac{1}{2 \tau}) \sum_{i} \frac{1}{2} \vec{e_{i}} \cdot \vec{e_{i}} \vec{e_{i}} f_{i}^{(1)}$ to prove that its divergence is equal to $\nabla \cdot \mathbf{q} - \nabla \cdot (\sigma \cdot \mathbf{u})$, where $\mathbf{q}$ is the heat flux and $\sigma$ is the Cauchy stress tensor:

$\nabla_{1} \cdot ((1-\frac{1}{2 \tau}) \sum_{i} \frac{1}{2} \vec{e_{i}} \cdot \vec{e_{i}} \vec{e_{i}} f_{i}^{(1)}) = \nabla \cdot \mathbf{q} - \nabla \cdot (\sigma \cdot \mathbf{u})$

Finally, we could find the total energy balance equation as:

$\frac{\partial}{\partial t} (\rho e) + \nabla \cdot (\rho e \mathbf{u}) = -\nabla \cdot \mathbf{q} + \nabla \cdot (\sigma \cdot \mathbf{u})$


问:
Identifying the time window at which a meteorological factor influences tree growth the most
This is a problem I am encountering in a research project. I have thought of a naive approach to solve it which relies on many assumptions, is not guaranteed to work, and would be tedious and computationally expensive to implement. Hence, I am wondering if there would be cleverer or more general approaches. Below I will try to explain my problem as clearly as possible.

Tree growth at day d is thought of as a function of meteorology that has happened in the past. For instance, tree growth on a given day d may depend, on the cumulative temperature between day d-1 and day d-20 but probably not on the temperature a year before. Likewise for all the other meteorological variables.

I have daily meteorological data and daily tree growth data measured in number of cells created.

For any given meteorological variable, I would like to know for which period in the past the variable influences tree growth the most (in terms of daily creation of cells).

Naive approach:

Let's choose any meteorological variable, like temperature.

I have daily mean temperature data $T(d)$ and daily cell creation data $C(d)$. These are discrete functions but I can extend them to all time by interpolation.

Solving the problem is equivalent to finding values $0 < \lambda \leq \mu$ which maximise the correlation

$$K(\mu,\lambda) = cor(\: \frac{1}{2 \lambda}\int_{(d - \mu) - \lambda}^{(d - \mu) +\lambda} T(t) \: dt,\quad C(d)\:)$$

$2\lambda$ would then be length of the time interval which most influences growth and $\mu$ the number of days in the past at the center of the interval.

We can refine this further. There is probably a sweet spot in the past for which temperature particularly influences present growth. Conversely, temperatures further away in time from that sweet spot influence present growth less. To take this into account we can weigh the integration with a normal distribution.

Solving the problem then becomes maximising

$$\tilde{K}(\mu, \sigma, \lambda) = cor( \: \frac{\int_{(d-\mu) - \lambda}^{(d-\mu) +\lambda} T(t)\:F_{d-\mu, \sigma}(t) \: dt}{\int_{(d-\mu) - \lambda}^{(d-\mu) +\lambda} F_{d-\mu, \sigma}(t) \: dt},\:  C(d)\:)$$

for $\sigma \geq 0, \: \mu \geq \lambda \geq 0 $,

where $F_{d-\mu, \sigma}$ is the density function of a normal distribution of mean $d-\mu$ and standard deviation $\sigma$.

Lastly, I am coding in R and the correlation methods available are "Pearson" (the usual one), "Spearman" and "Kendall". Pearson mostly measures linear relationships. However, I would not expect growth to be linearly dependent with the temperature. There could be a polynomial relation with optimum temperatures causing maximum growth and temperatures far from that optimum causing less or no growth. For that reason I was going to use Spearman correlation which better captures non-linear relationships.

With the above approach I need to assume that the meteorological influences (e.g. temperature and rainfall) are not dependent, but they probably are. I am treating each meteorological factor separately for now.

Going back to questions:

Do you have any thoughts on how I could tackle this problem?

Have you spotted significant flaws in this method that could be addressed?

Do you know of better methods to solve this problem?

Do you know of similar problems that have already been solved, which I could look into for a solution?

Thanks a lot for your help!


问:
How can Navier--Stokes equations have asymmetric solutions such as Karman vortex streets
The Navier--Stokes equations are axially symmetric, so with symmetric boundary conditions, how can features such as Karman vortex streets develop?

I understand that in reality symmetry does never completely hold due to minimal fluctuations of geometries, densities, etc., but in the linked article, vortex streets develop even in numerical simulations


问:
What does the Jackson Kernel measure?
A certain filter I'm writing uses two different kernels.  The Fejer kernel (which is common) and the Jackson kernel:

$$ \Delta_T(x) = T \,\left(  \frac{\sin \pi T x}{\pi T x}\right)^2  \quad\text{and}\quad
J_T(x) = \frac{3T}{4} \left(  \frac{\sin \pi T x/2}{\pi T x/2}\right)^4
 $$

So I had to remind myself how these two kernels behaved and what their properties were:

For a fairly broad class of functions, convolving with the Fejer kernel is the same as integrating against a triangle:

$$ \int_{-T}^T \left( 1 - \frac{|\xi|}{T} \right) \hat{f}(\xi) \,e^{2\pi i \, x \xi} \,d\xi = (f * \Delta_T)(x) $$

Notably, the discrete Fejér kernel (over $S^1$) is a limit of the continuous Fejer kernel (over $\mathbb{R}$):

$$ F_N(x) = \frac{D_0(x) + \dots + D_{N-1}(x)}{N}  = \frac{1}{N} \left(\frac{\sin \left(Nx/2\right)}{\sin (x/2)}\right)^2 = \Delta_{1/T}(x)$$

Here it emulates Cesaro summabilition on the Fourier series:

$$ \big(f*F_N\big)(x) = \frac{ S_0(f) (x) + \dots + S_{N-1}(f)(x)}{N} = \frac{1}{N} \sum_{n=0}^{N} \sum_{|x| \leq n} \widehat{f}(x)$$

The similarities continue that the $S^1$ Fejér kernel is the periodization of the $\mathbb{R}$ Fejér kernel:

$$ \sum_{n \in \mathbb{Z}} \Delta_{\frac{1}{N}}(x + n) = F_N(x)  $$

These are taken from Fourier Analysis by Elias Stein.  My question is what similar properties should occur for the Jackson kernel?

Googling the Jackson kernel doesn't return much. E.g. this note on trigonometric polynomials and modulus of continuity.

Title: What does the Jackson Kernel measure?


问:
Maintain sorted ring buffer
I would like to insert elements into a ring (circular) buffer one at a time and maintain a permutation array which keeps track of the sorted elements in ascending order. To do this, I have adapted the "Insertion Sort" algorithm to use circular array indexing:

```<code>/* sort ringbuf entries with insertion sort algorithm */
void
ringbuf_isort(const ringbuf *rbuf, int *idx)
{
  const int n = ringbuf_n(rbuf);
  int i, j;

  for (i = 1; i < n; i++)
    {
      int id = idx[i];
      double v = rbuf->array[(rbuf->head + id) % rbuf->size];

      for (j = i; j > 0; j--)
        {
          if (rbuf->array[(rbuf->head + idx[j-1]) % rbuf->size] < v)
            break;

          idx[j] = idx[j-1];
        }

      idx[j] = id;
    }
}
</code>```

This algorithm works correctly - the idx array contains the permutation of the buffer elements which tracks their sorted status. The problem is that the algorithm is extremely slow, due to the modular arithmetic needed for the circular indexing, i.e. (head + idx) % size.

In fact, I have found it is faster to simply copy the ring buffer into a linear array and use a quicksort algorithm each time an element is added, rather than use the code above. For a linear array, I have confirmed the insertion sort algorithm outperforms quicksort when only a single element is changed, so I am hoping I can find a better solution for the ring buffer.

Does anyone know how to optimize the code above to speed up the modular arithmetic (or avoid it altogether?)

Title: Maintain sorted ring buffer


问:
An invertible matrix that minimizes the norm of the product with a given matrix
Given a fat matrix $B \in \mathbb{C}^{n \times m}$ (where $m > n$) with full row rank, I would like to find (numerically) a full-rank matrix $A$ that minimizes the Frobenius norm of the product $A B$. Formally,

$$\underset{A \in \mathbb{C}^{n \times n}}{\text{minimize}} \quad \frac{1}{2} \|AB\|_F^2 \quad \text{subject to} \quad \det (A) \neq 0$$

The value of $m$ is typically an order of magnitude larger than $n$. The sizes ($n,m$) I am interested in may be on the order of hundreds.

I have found the following discussion, which I guess could be generalized to the above case. I wonder if there is a simpler solution in this particular scenario.


问:
Showing date in Paraview's Annotate Time Filter
I have a time dependent data set where each frame corresponds to a certain date. I need to show the time annotation, much like what Annotate Time Filter does, but show the date instead of time as a float value. How can I do that?

In python I can convert the time index to a date with datetime for instance:

```<code>init_date = datetime.datetime(2018, 1, 1)
time_step = 3600.
itime = 20
date = init_date + datetime.timedelta(seconds=itime*time_step)
print(date.strftime('%Y-%m-%d %H:%M'))
</code>```


问:
Stochastic conjugate directions to improve convergence in narrow valleys
My question concerns a specific statement in this paper:

N. N. Schraudolph and T. Graepel, "Conjugate Directions for Stochastic Gradient Descent," in Int. Conf. Artificial Neural Networks, Berlin, 2002 pp. 1351–1356. The preprint is available here.

Schraudolph and Graepel address the problem that in stochastic settings of ill-conditioned problems the sampled estimate of the gradient is very unlikely to point in the right direction due to the ill-conditioning of the (estimated) Hessian.

In Chapter 3, the authors claim that directions associated with large eigenvalues of the Hessian (directions of large curvature) can be identified by multiplying the stochastic estimates of Hessian $H$ and gradient $g$.

I don't understand this argument. What point am I missing?


问:
Solving an ODE from many starting points
I have a situation where I'm interested in solving the same ODE many times, from different initial locations. More precisely, I'm interested in solving an ODE of the form

\begin{align}
\frac{dx}{dt} &= f(x) \\
x(0) &= a \in \mathbb{R}^n
\end{align}

where I know that this ODE has a $T$-periodic orbit $\mathcal{A}$, i.e for $a \in \mathcal{A}$, we have that

\begin{align}
x(0) &= a \implies x(T) = a
\end{align}

My problem is such that at each iteration $k$, I am given a point $a_k \in \mathcal{A}$, and need to output where the ODE takes $a_k$ after time $t$, i.e.

\begin{align}
\frac{dx}{dt} &= f(x) \\
x(0) &= a_k \\
\implies b_k &\triangleq x(t)
\end{align}

and I want to report $b_k$.

My question is: is there a better way of doing this than using an ODE solver at every step $k$? For example, if possible, I would prefer to incur a one-off cost at the start of my algorithm which allows me to map directly from $a_k$ to $b_k$, rather than having to repeatedly solve the same problem.

Additional details:

The system I'm interested in is a two-dimensional Hamiltonian system, given in coordinates as

\begin{align}
\frac{dq}{dt} &= |p|^{\alpha - 1} \cdot \text{sign} (p) \\
\frac{dp}{dt} &= -|q|^{\beta - 1} \cdot \text{sign} (q)
\end{align}

where $\alpha, \beta > 1$ such that $1/\alpha + 1/\beta = 1$.

The standard integrator for such systems is the leapfrog integrator, which I would use if I was solving for a single trajectory. As written, the system is non-smooth - I would rather not worry about that right now, and if it makes the question easier to answer, please feel free to assume that everything involved is smooth.

When $\alpha = \beta = 2$ things are much easier, and I can handle this situation fine - I'm asking about the other cases.