问:
Calculating a submanifold of minimal values
I have a (numerical, noisy) function $f$ that I'd like to optimize (in 3 < dimension < 10).

The function $f$ has a fairly prominent valley/ravine structure due to physical symmetries present in the input data. To be more precise: there is a submanifold $V$ of low dimension (let's say 1 or 2) along which function values are close to the global minimum $x \in V$. I would like to approximate the submanifold $V$.

The problem of such valleys/ravines is well-known in the literature, any first course on numerical optimization will mention the Rosenbrock function. However, I keep coming up short when looking for methods of computing/parametrizing the valley itself. I guess I could look at the gradient flow and treat the valley as the slow manifold of that, but due to the fact that my valley potentially has lots of local minima due to noise, I don't expect this to work particularly well.

What is the most common way of dealing with this kind of problem?


问:
'Items in a bucket' optimizaiton problem?
Sorry for question's silly tittle. I don't quite know how to name my problem, that's why I'm here.

Imagine you wish do add the collumns of you data ('1','2','3'...'30k' etc.) in such a way as to obtain the numbers in the Goal column.

You don't have to combine all the columns, combine only as many as necessary for you to obtain the column 'Goal'. In fact, my objective is to know which columns lead me to the column Goal. I may have to combine 2, 3, 4, 30k columns, doesn't matter as long as I find the 'shortest' solution possible, namely using as few collumns as possible.

My analogy to this problem would be something like puttings items in a blender until I reach the desired outcome. The order in which I add the collumns doesn't matter.

I mean to minimize the mean abstolute error (MAE) between the combination of columns and the Goal.

I have been reading about integer programming, but I don't quite know where exactly I should go to. Could anybody point me to which class of algorithms I can use to approach this type of problem?

ContentLicense: CC BY-SA 4.0


问:
Can I build an electrostatic simulation of conductive objects using just the method of images?
I wish to build a simulation of self-assembling wires. Think lots of ball bearings in a 20KV electric field and they form emergent structures.

If I use the method of images to figure out what the induced charges on various parts of these ball bearings due to all the other charges, is this enough to specify a general simulation that should computationally reproduce the experiment?

Are there any python, R, C#, or other libraries that already do this? Note, that I want to do it in 3D.


问:
Can we sparse solve a few eigenvalues specified by index range?
I need to solve a few eigenvalues of a large sparse matrix specified by their index range. These indices are according to the whole eigenspectrum sorted in algebraic (not absolute value) ascending order.

For example in python, I notice this is made available by the subset_by_index option in scipy.linalg.eigh for non-sparse matrices, but not in the sparse counterpart eigsh. In MKL, we have non-sparse routines like heevx with the similar option but I don't find clear info in the sparse docs.

Is this possible at all in the sparse context?


问:
How to implement large rotations in total lagrangian formulation (nonlinear FEM)?
I have developed an Octave script to solve the nonlinear Euler-Bernoulli beam equations with linearized von Karman-strains, i.e. higher-order terms are dropped. The simulation results agree with analytical results in small and moderately large displacements in static analysis. However, when doing a dynamic simulation with large displacements I get incorrect results.

I used this site as a reference for assembling the stiffness matrix and the tangent stiffness matrix. I should note that I think the site has a typo for the stiffness coefficients $K_{21}$. I think the 0.5 should be dropped from the $K_{21}$ expression and also for the tangent stiffness matrix there should be $T_{21}=K_{21}$ instead of $T_{21}=2K_{21}$ (that is, if you don't symmetrize the matrices).

Now as the static linear analysis seems to work fine, I tried a dynamic nonlinear analysis. For this I used a simple pendulum with one beam element (length L=1m) starting from rest and a horizontal position. For the first few time steps everything seems fine until the beam element starts elongating (see figures below). At t=0.5 the length of the beam is almost 1.5m, which should not be possible. The internal force vector is initially pointing vertically upward but starts oscillating along with the inertial force vector and is never really aligned with the length of the beam as you would excpect for a tension force. Also if I do a static analysis and rotate the beam 90 degrees with Dirichlet BCs, I get large moments acting on the beam even though we are simulating rigid body rotation.

Pendulum at t=0.32s, Pendulum at t=0.5s

So my question is, is it possible to simulate large rotations in the total Lagrangian description (assuming small strains) and using the von Karman strains? Or is there some other reason why the simulation is not working? Or is my time integration scheme not suitable for this problem?

For the dynamic analysis I used the Newmark method for time integration (trapezoid rule), with boundary conditions set to u(0)=0 and w(0)=0. The stiffness matrices are integrated with gauss quadrature with reduced integration for the nonlinear terms, i.e. $K_{12}$, $K_{21}$, $T_{12}$, $T_{21}$ and $K_{22}^{NL}$ and full integration for the linear terms.


问:
Eigen-decomposition one eigenpair by one eigenpair?
Is it possible to conduct an Eigen-decomposition of a matrix one eigenpair by one eigenpair?

And related to this question, what is the time complexity of truncated eigendecomposition?

I am trying (hard) to find a way to reduce the time complexity of my algorithm which involves the calculation of eigendecomposition of a positive definite matrix. The full eigendecomposition is $O(n^3)$ in general. However, I may not need full eigendecomposition. I only need to stop when the $\lambda_i$, i.e. the $i$-th largest eigenvalue (and its eigenvector), is less than $\epsilon>0$.

Is there a way to do that? We might assume we know the number of eigenvalues greater or equal to $\epsilon$, if necessary, e.g. something like truncated eigendecomposition.


问:
Cyipopt fails to converge for NLP problem which fmincon() can solve
I'm currently trying to implement a python script for solving a constrained nonlinear optimization problem with ~800 variables and 2 constraints, one linear and one nonlinear. There already exists a Matlab implementation of this script and I'm trying to use cyipopt to implement an equivalent Python implementation that produces similar results.

The Matlab code has analytic solutions for both the objective and constraint and jacobian, which I've rewritten entirely in Python, and uses L-BFGS for hessian approximation. Calling the Python jacobian functions from fmincon() results in the same solution so I think these functions are correct, and I've also used the derivative checkers in both fmincon() and cyipopt to make sure. However, my Python implementation doesn't converge and fails to match the performance of the Matlab code. In particular, the objective doesn't decrease smoothly. Here's some example output from Ipopt:

```<code>This is Ipopt version 3.14.4, running with linear solver MUMPS 5.2.1.

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

Total number of variables............................:      801
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      267
                     variables with only upper bounds:        0
Total number of equality constraints.................:        2
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.6958872e-01 5.68e-14 1.04e-02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.1548097e-01 3.25e+00 7.09e+00   0.6 3.87e+00    -  9.90e-01 4.34e-01f  1
   2  1.7668628e-01 3.05e+00 8.15e+00   0.9 2.38e+00    -  1.00e+00 9.72e-01f  1
   3  1.7381539e-01 1.78e-01 3.39e+00   0.2 1.09e+00    -  7.37e-01 1.00e+00h  1
   4  1.7270930e-01 4.47e-03 2.42e-01  -0.5 2.10e-01    -  1.00e+00 1.00e+00h  1
   5  1.7178940e-01 4.69e-03 4.70e-02  -1.2 7.62e-02    -  1.00e+00 1.00e+00f  1
   6  1.7055753e-01 5.01e-03 8.89e-03  -1.9 1.21e-01    -  1.00e+00 1.00e+00h  1
   7  1.6692096e-01 1.04e-02 1.15e-03  -2.9 7.68e-02    -  1.00e+00 1.00e+00h  1
   8  1.6323948e-01 3.18e-01 8.33e-03  -3.3 5.42e-01    -  9.42e-01 1.00e+00h  1
   9  1.6005166e-01 1.43e-01 6.64e-03  -9.3 2.50e-01    -  5.61e-01 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.5786474e-01 1.22e-01 4.18e-03  -3.7 4.04e-01    -  7.71e-01 1.00e+00h  1
  11  1.5672747e-01 6.04e-02 5.11e-03  -4.2 1.61e-01    -  9.61e-01 1.00e+00h  1
  12  1.5472463e-01 2.40e-02 2.34e-03  -4.8 1.23e-01    -  9.96e-01 1.00e+00h  1
  13  1.5533721e-01 9.66e-04 6.81e-03  -5.0 1.84e-01    -  9.99e-01 1.00e+00H  1
  14  1.5377085e-01 3.02e-02 3.08e-03  -5.3 1.65e-01    -  1.00e+00 1.00e+00h  1
  15  1.5361137e-01 2.09e-02 4.51e-03  -5.7 1.11e-01    -  1.00e+00 1.00e+00h  1
  16  1.5296208e-01 6.95e-03 1.95e-03  -6.0 8.23e-02    -  1.00e+00 1.00e+00h  1
  17  1.5273048e-01 1.90e-03 1.01e-03  -7.0 3.83e-02    -  1.00e+00 1.00e+00h  1
  18  1.5311695e-01 4.91e-05 5.20e-03  -8.0 2.13e-01    -  1.00e+00 1.00e+00H  1
  19  1.5262106e-01 7.92e-03 6.41e-03  -5.9 4.57e-01    -  1.00e+00 2.30e-01h  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  1.5261389e-01 1.19e-02 1.27e-02  -3.9 4.90e-01    -  1.00e+00 1.02e-01h  2
  21  1.5310130e-01 4.88e-02 4.09e-03  -3.6 2.90e-01    -  4.10e-01 1.00e+00h  1
  22  1.5330300e-01 2.75e-02 3.72e-03  -3.9 1.29e-01    -  1.00e+00 1.00e+00h  1
  23  1.5201654e-01 1.23e-02 1.40e-03  -3.9 9.69e-02    -  1.00e+00 1.00e+00h  1
  24  1.5227903e-01 3.92e-04 6.45e-03  -4.8 2.11e-01    -  9.86e-01 1.00e+00H  1
  25  1.5495353e-01 1.86e-04 1.26e-02  -4.5 3.70e-01    -  1.00e+00 1.00e+00H  1
  26  1.5557422e-01 3.93e-05 4.07e-03  -4.6 3.53e-01    -  1.18e-01 1.00e+00H  1
  27  1.5281467e-01 7.97e-02 4.35e-03  -4.3 3.09e-01    -  1.00e+00 1.00e+00h  1
  28  1.6236141e-01 3.83e-02 1.24e-02  -4.4 4.15e-01    -  7.25e-01 8.80e-01H  1
  29  1.5299068e-01 8.52e-04 6.17e-03  -4.4 2.48e-01    -  7.17e-02 1.00e+00F  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  1.5139425e-01 2.53e-02 3.94e-03  -4.4 9.96e-01    -  2.38e-01 1.55e-01h  3
  31  1.5138493e-01 5.50e-04 3.76e-03  -4.8 2.40e-01    -  9.94e-01 1.00e+00H  1
  32  1.5200809e-01 2.14e-04 5.16e-03  -4.9 1.32e-01    -  9.96e-01 1.00e+00H  1
  33  1.5074554e-01 1.86e-02 1.91e-03  -5.0 1.77e-01    -  9.98e-01 1.00e+00h  1
  34  1.5130813e-01 1.49e-04 3.53e-03  -5.6 1.56e-01    -  1.00e+00 1.00e+00H  1
  35  1.5036412e-01 1.28e-02 1.63e-03  -5.7 1.18e-01    -  1.00e+00 1.00e+00h  1
  36  1.5026330e-01 6.65e-04 8.92e-04  -7.2 3.49e-02    -  1.00e+00 1.00e+00h  1
  37  1.5076495e-01 1.34e-04 5.36e-03  -7.7 1.20e-01    -  1.00e+00 1.00e+00H  1
  38  1.5066012e-01 1.69e-02 6.61e-03  -7.9 1.11e-01    -  1.00e+00 1.00e+00h  1
  39  1.5005302e-01 1.25e-02 1.81e-03  -8.0 1.05e-01    -  5.44e-01 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  1.5001109e-01 3.87e-03 2.08e-03  -8.6 4.98e-02    -  1.00e+00 1.00e+00h  1
  41  1.4987818e-01 2.98e-03 2.39e-03  -9.7 5.14e-02    -  1.00e+00 1.00e+00h  1
  42  1.5014976e-01 2.70e-05 3.82e-03  -9.8 1.02e-01    -  1.00e+00 1.00e+00H  1
  43  1.4997585e-01 1.01e-02 5.16e-03 -10.4 6.66e-02    -  1.00e+00 1.00e+00h  1
  44  1.4970903e-01 9.51e-03 2.36e-03 -10.5 8.90e-02    -  8.37e-01 1.00e+00h  1
  45  1.4940598e-01 2.37e-03 9.31e-04 -11.0 4.20e-02    -  1.00e+00 1.00e+00h  1
  46  1.4968497e-01 2.68e-05 4.63e-03 -11.0 1.21e-01    -  1.00e+00 1.00e+00H  1
  47  1.4995823e-01 1.29e-06 7.01e-03 -11.0 3.32e-01    -  1.00e+00 1.00e+00H  1
  48  1.4974993e-01 1.46e-02 6.86e-03  -9.0 2.03e+00    -  1.00e+00 1.18e-01f  2
  49  1.4948599e-01 1.72e-02 5.95e-03  -9.1 7.57e-01    -  2.53e-01 8.53e-02h  3
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  1.5186145e-01 6.77e-04 8.93e-03  -9.1 2.87e-01    -  7.46e-01 1.00e+00H  1
  51  1.5120058e-01 1.35e-03 4.47e-03  -9.1 3.78e-01    -  1.19e-01 1.00e+00H  1
  52  1.5550082e-01 1.54e-04 8.87e-03  -9.1 3.78e-01    -  5.42e-01 1.00e+00H  1
  53  1.5020757e-01 1.12e-04 4.29e-03  -9.1 4.29e-01    -  1.00e+00 1.00e+00H  1
  54  1.4918169e-01 1.28e-02 1.97e-03  -9.1 1.68e+00    -  2.17e-01 6.23e-02h  3
  55  1.5133875e-01 4.49e-05 8.71e-03  -9.5 2.28e-01    -  9.87e-01 1.00e+00H  1
  56  1.4882545e-01 2.34e-02 2.51e-03  -9.6 1.72e-01    -  4.77e-01 1.00e+00h  1
  57  1.4946850e-01 1.48e-02 6.52e-03 -10.4 1.07e-01    -  1.00e+00 1.00e+00h  1
  58  1.4859822e-01 5.15e-03 2.33e-03 -10.5 6.48e-02    -  9.27e-01 1.00e+00h  1
  59  1.4902635e-01 7.60e-05 4.53e-03 -11.0 1.06e-01    -  1.00e+00 1.00e+00H  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  60  1.4843429e-01 4.65e-03 1.51e-03 -11.0 7.41e-02    -  1.00e+00 1.00e+00h  1
  61  1.4839694e-01 1.37e-03 2.49e-03 -11.0 3.13e-02    -  1.00e+00 1.00e+00h  1
  62  1.4830507e-01 3.51e-04 8.59e-04 -11.0 2.56e-02    -  1.00e+00 1.00e+00h  1
  63  1.4849382e-01 1.36e-05 5.56e-03 -11.0 1.22e-01    -  1.00e+00 1.00e+00H  1
  64  1.4823266e-01 3.86e-03 1.66e-03 -11.0 9.19e-02    -  1.00e+00 1.00e+00h  1
  65  1.4821480e-01 5.60e-04 1.20e-03 -11.0 1.53e-02    -  1.00e+00 1.00e+00h  1
  66  1.4818681e-01 1.50e-04 2.13e-04 -11.0 1.25e-02    -  1.00e+00 1.00e+00h  1
  67  1.4816045e-01 1.46e-04 5.02e-04 -11.0 1.61e-02    -  1.00e+00 1.00e+00h  1
  68  1.4850582e-01 5.85e-05 5.29e-03 -11.0 1.57e-01    -  1.00e+00 1.00e+00H  1
  69  1.5001793e-01 1.58e-06 1.08e-02 -11.0 1.41e-01    -  1.00e+00 1.00e+00H  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  70  1.4861627e-01 2.68e-07 5.43e-03 -11.0 2.13e-01    -  2.10e-01 1.00e+00F  1
  71  1.4842174e-01 1.21e-02 4.02e-03 -10.6 3.61e-01    -  1.00e+00 2.50e-01f  3
  72  1.4816338e-01 4.71e-03 2.55e-03 -10.7 6.93e-02    -  1.00e+00 1.00e+00h  1
  73  1.4841277e-01 1.20e-04 5.26e-03 -10.7 1.25e-01    -  1.00e+00 1.00e+00H  1
  74  1.4923205e-01 4.65e-05 7.02e-03 -11.0 9.58e-02    -  1.00e+00 1.00e+00H  1
  75  1.4793006e-01 1.50e-02 1.22e-03 -11.0 1.02e-01    -  9.33e-01 1.00e+00h  1
  76  1.4790107e-01 1.98e-03 2.30e-03 -11.0 4.36e-02    -  1.00e+00 1.00e+00h  1
  77  1.4795403e-01 6.21e-06 2.74e-03 -11.0 3.90e-02    -  1.00e+00 1.00e+00H  1
  78  1.4783354e-01 3.00e-03 1.85e-03 -11.0 5.85e-02    -  1.00e+00 1.00e+00h  1
  79  1.4775646e-01 1.00e-03 1.66e-03 -11.0 2.36e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  80  1.4772077e-01 2.46e-04 8.16e-04 -11.0 1.37e-02    -  1.00e+00 1.00e+00h  1
  81  1.4772337e-01 1.36e-06 1.67e-03 -11.0 2.97e-02    -  1.00e+00 1.00e+00H  1
  82  1.4800983e-01 3.43e-07 4.69e-03 -11.0 1.73e-01    -  1.00e+00 1.00e+00H  1
  83  1.4785970e-01 3.85e-03 5.48e-03 -10.2 1.53e+00    -  1.00e+00 7.61e-02f  3
  84  1.4778198e-01 6.41e-03 6.58e-03 -10.3 3.42e+00    -  6.76e-02 2.92e-02h  3
  85  1.4773098e-01 1.02e-02 5.56e-03 -10.3 2.22e+00    -  6.49e-02 4.72e-02h  3
  86  1.5775339e-01 1.78e-04 7.17e-03 -10.3 1.27e+00    -  1.35e-01 1.00e+00H  1
  87  1.5716495e-01 1.63e-03 5.41e-03 -10.3 1.70e+00    -  1.70e-01 3.17e-02h  4
  88  1.5247575e-01 1.33e-04 4.44e-03 -10.3 3.35e-01    -  5.46e-01 1.00e+00H  1
  89  1.5188072e-01 1.15e-03 3.96e-03 -10.3 1.80e+00    -  7.48e-01 3.87e-02h  5
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  90  1.5109270e-01 2.99e-03 3.54e-03 -10.3 1.65e+00    -  4.22e-01 5.10e-02h  4
  91  1.4947946e-01 5.89e-02 4.20e-03 -10.3 1.50e+01    -  4.46e-02 2.31e-02h  2
  92  1.4898332e-01 5.71e-02 3.70e-03 -10.3 6.05e-01    -  3.30e-01 1.30e-01h  3
  93  1.4893398e-01 4.01e-02 6.60e-03 -10.3 1.73e-01    -  1.00e+00 1.00e+00h  1
  94  1.4813809e-01 2.83e-02 2.88e-03 -10.3 2.64e-01    -  1.00e+00 1.00e+00h  1
  95  1.4953181e-01 5.14e-04 7.13e-03 -10.3 1.61e-01    -  1.00e+00 1.00e+00H  1
  96  1.4757049e-01 2.19e-02 1.86e-03 -10.3 1.12e-01    -  1.00e+00 1.00e+00h  1
  97  1.4838902e-01 1.46e-02 5.65e-03 -10.8 1.05e-01    -  1.00e+00 1.00e+00h  1
  98  1.4742559e-01 8.28e-03 5.50e-04 -10.9 7.16e-02    -  1.00e+00 1.00e+00h  1
  99  1.4739592e-01 1.63e-04 7.52e-04 -11.0 1.75e-02    -  1.00e+00 1.00e+00h  1
</code>```

Compared to fmincon():

```<code>                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    1.695878e-01    5.684e-14    1.102e-02
    1       2    1.691645e-01    1.296e-04    4.729e-03    2.229e-02
    2       3    1.665545e-01    7.175e-03    9.554e-03    1.857e-01
    3       4    1.635913e-01    4.118e-02    6.511e-03    5.518e-01
    4       5    1.626062e-01    3.235e-02    1.080e-02    6.491e-01
    5       6    1.617645e-01    2.861e-02    3.397e-03    7.418e-01
    6       7    1.607960e-01    1.297e-02    2.636e-03    6.037e-01
    7       8    1.595868e-01    4.010e-02    4.019e-03    1.089e+00
    8       9    1.589008e-01    2.304e-02    4.341e-03    7.246e-01
    9      10    1.581489e-01    3.965e-03    1.667e-03    2.636e-01
   10      11    1.575713e-01    1.046e-02    1.000e-03    6.139e-01
   11      12    1.573578e-01    3.968e-03    1.179e-03    4.697e-01
   12      13    1.570832e-01    4.880e-03    1.247e-03    6.278e-01
   13      14    1.567794e-01    5.258e-03    1.000e-03    6.424e-01
   14      15    1.567468e-01    7.916e-03    2.665e-03    6.705e-01
   15      16    1.565624e-01    1.995e-03    1.452e-03    1.778e-01
   16      17    1.565058e-01    9.510e-04    1.309e-03    1.582e-01
   17      18    1.564559e-01    8.874e-04    1.000e-03    2.408e-01
   18      19    1.556438e-01    2.632e-03    7.432e-04    2.010e-01
   19      21    1.542575e-01    1.763e-05    1.908e-03    6.447e-01
   20      23    1.538128e-01    8.862e-06    2.772e-03    3.140e-01
   21      24    1.535477e-01    1.020e-03    1.317e-03    1.241e-01
   22      25    1.531595e-01    2.294e-03    7.816e-04    2.747e-01
   23      26    1.529536e-01    5.902e-04    8.447e-04    1.609e-01
   24      28    1.524555e-01    4.137e-07    1.424e-03    4.367e-01
   25      30    1.521461e-01    2.148e-06    1.428e-03    2.705e-01
   26      31    1.520693e-01    3.572e-04    9.764e-04    7.670e-02
   27      33    1.518785e-01    1.542e-06    1.821e-03    3.113e-01
   28      34    1.516523e-01    7.795e-04    1.052e-03    1.813e-01
   29      35    1.515454e-01    3.468e-04    6.490e-04    1.431e-01
   30      36    1.514033e-01    6.813e-04    7.628e-04    2.445e-01

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   31      37    1.511966e-01    1.321e-03    1.065e-03    3.590e-01
   32      38    1.509549e-01    2.310e-03    9.453e-04    5.227e-01
   33      39    1.507060e-01    2.390e-03    1.893e-03    5.238e-01
   34      40    1.505274e-01    1.191e-03    4.680e-04    2.007e-01
   35      41    1.504018e-01    3.995e-04    7.486e-04    1.168e-01
   36      42    1.501441e-01    8.662e-04    5.371e-04    3.376e-01
   37      44    1.498647e-01    3.364e-06    1.694e-03    5.798e-01
   38      45    1.496842e-01    9.911e-04    7.416e-04    3.829e-01
   39      46    1.496415e-01    3.685e-04    5.740e-04    1.727e-01
   40      47    1.495934e-01    9.745e-04    1.368e-03    3.101e-01
   41      48    1.494941e-01    6.501e-04    6.743e-04    3.020e-01
   42      49    1.493509e-01    9.861e-04    3.953e-04    3.974e-01
   43      50    1.492265e-01    7.372e-04    6.141e-04    3.550e-01
   44      51    1.490722e-01    1.345e-03    7.354e-04    4.686e-01
   45      52    1.489949e-01    9.143e-04    5.345e-04    2.969e-01
   46      53    1.489130e-01    8.080e-04    7.258e-04    2.955e-01
   47      54    1.488745e-01    1.176e-03    1.445e-03    2.992e-01
   48      55    1.488323e-01    4.312e-04    5.706e-04    1.628e-01
   49      56    1.487758e-01    3.440e-04    8.374e-04    1.993e-01
   50      57    1.486999e-01    1.017e-03    6.874e-04    3.923e-01
   51      59    1.486946e-01    8.701e-06    3.005e-03    4.956e-01
   52      60    1.485959e-01    3.722e-04    9.674e-04    7.810e-02
   53      61    1.485395e-01    3.153e-04    5.199e-04    1.871e-01
   54      62    1.484944e-01    3.883e-04    9.311e-04    2.667e-01
   55      63    1.484078e-01    6.886e-04    1.208e-03    4.148e-01
   56      64    1.482546e-01    2.396e-03    1.258e-03    9.236e-01
   57      65    1.481228e-01    2.799e-03    1.541e-03    9.661e-01
   58      66    1.479890e-01    7.611e-04    5.626e-04    4.157e-01
   59      67    1.478692e-01    1.020e-03    9.555e-04    5.353e-01
   60      68    1.477201e-01    1.205e-03    2.336e-03    6.695e-01

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   61      69    1.474560e-01    3.556e-03    1.588e-03    1.439e+00
   62      70    1.471856e-01    2.773e-03    1.294e-03    1.377e+00
   63      71    1.469532e-01    1.536e-03    1.162e-03    8.180e-01
   64      72    1.467492e-01    1.390e-03    8.390e-04    7.719e-01
   65      73    1.464826e-01    2.369e-03    2.728e-03    1.328e+00
   66      74    1.461599e-01    2.225e-03    1.697e-03    1.276e+00
   67      75    1.458142e-01    2.839e-03    7.377e-04    1.729e+00
   68      76    1.455273e-01    1.941e-03    1.540e-03    1.555e+00
   69      77    1.451421e-01    2.727e-03    2.188e-03    2.072e+00
   70      78    1.446933e-01    6.136e-03    3.128e-03    2.527e+00
   71      79    1.442774e-01    4.200e-03    1.005e-03    1.977e+00
   72      80    1.440018e-01    1.392e-03    1.141e-03    4.495e-01
   73      81    1.436448e-01    3.002e-03    2.341e-03    1.335e+00
   74      82    1.432150e-01    4.996e-03    3.013e-03    1.999e+00
   75      83    1.427922e-01    3.273e-03    7.698e-04    1.622e+00
   76      84    1.424515e-01    1.755e-03    1.547e-03    1.999e+00
   77      85    1.420710e-01    1.985e-03    1.608e-03    2.162e+00
   78      86    1.414360e-01    5.327e-03    1.684e-03    3.757e+00
   79      88    1.410962e-01    8.909e-05    5.970e-03    7.093e+00
   80      89    1.405116e-01    3.316e-03    1.177e-03    8.055e-01
   81      90    1.402913e-01    1.138e-03    1.549e-03    3.040e-01
   82      91    1.398874e-01    3.111e-03    2.327e-03    1.897e+00
   83      93    1.391672e-01    3.531e-06    3.701e-03    5.891e+00
   84      95    1.383758e-01    5.642e-05    3.542e-03    6.218e+00
   85      96    1.378290e-01    4.166e-03    1.036e-03    2.995e+00
   86      97    1.375162e-01    1.780e-03    1.041e-03    1.597e+00
   87      98    1.370376e-01    5.208e-03    1.925e-03    3.465e+00
   88      99    1.366982e-01    7.099e-03    4.238e-03    3.143e+00
   89     100    1.362886e-01    2.651e-03    2.194e-03    2.074e+00
   90     101    1.360006e-01    2.295e-03    1.189e-03    1.023e+00

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   91     102    1.358196e-01    9.869e-04    1.742e-03    1.240e+00
   92     103    1.355212e-01    2.761e-03    1.672e-03    1.891e+00
   93     105    1.353705e-01    1.162e-04    4.479e-03    4.375e+00
   94     106    1.348992e-01    2.477e-03    1.632e-03    3.747e-01
   95     107    1.347110e-01    1.267e-03    1.224e-03    3.268e-01
   96     108    1.345712e-01    8.739e-04    1.756e-03    8.004e-01
   97     109    1.343092e-01    3.165e-03    2.251e-03    1.691e+00
   98     111    1.341921e-01    4.984e-03    2.419e-03    1.037e+00
   99     112    1.339379e-01    1.967e-03    9.087e-04    2.052e+00
</code>```

I don't really have any experience with either the theory or practice of optimization, so I'm wondering if anyone can help spot the issue here? Is there any way I can tweak the options I'm passing to IPOPT to solve the problem or is IPOPT simply less suited to this type of optimization problem?

For reference, here's the python code:

class initial_problem():
    def __init__(self, cell_pairs, dC):
        self.cell_pairs = cell_pairs
        self.dC = dC

    def objective(self, x):
        x = x.reshape(x0.shape, order='F')
        q = x[:,0:2]
        p = x[:,2]

        b = np.matmul(self.dC, np.multiply(q.T,p).T)
        delta_p = np.matmul(self.dC, p)

        t1 = np.divide((b - (np.multiply(r1.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r1.T,delta_p).T), axis=1)).T
        t2 = np.divide((b - (np.multiply(r2.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r2.T,delta_p).T), axis=1)).T

        E = 0.5 * np.mean(np.power(np.sum(t1 * tau_1, axis=1), 2) +
                          np.power(np.sum(t2 * tau_2, axis=1), 2))
        return E

    def gradient(self, x):
        # Jacobian of objective
        x = x.reshape(x0.shape, order='F')
        q = x[:,0:2]
        p = x[:,2]

        b = np.matmul(self.dC, np.multiply(q.T,p).T)
        delta_p = np.matmul(self.dC, p)

        t1 = np.divide((b - (np.multiply(r1.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r1.T,delta_p).T), axis=1)).T
        t2 = np.divide((b - (np.multiply(r2.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r2.T,delta_p).T), axis=1)).T

        ip1 = np.sum(t1 * tau_1, axis=1)
        ip2 = np.sum(t2 * tau_2, axis=1)

        drX1 = sx_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],0], q[self.cell_pairs[:,1],1], r1[:,0], r1[:,1])
        drX2 = sx_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],0], q[self.cell_pairs[:,1],1], r2[:,0], r2[:,1])
        drY1 = sy_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],0], q[self.cell_pairs[:,1],1], r1[:,0], r1[:,1])
        drY2 = sy_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],0], q[self.cell_pairs[:,1],1], r2[:,0], r2[:,1])

        dE = np.multiply(ip1 * tau_1[:,0], drX1.T).T + np.multiply(ip1 * tau_1[:,1], drY1.T).T + \
             np.multiply(ip2 * tau_2[:,0], drX2.T).T + np.multiply(ip2 * tau_2[:,1], drY2.T).T
        dE = dE / self.dC.shape[0]

        rows = np.concatenate([self.cell_pairs[:,0],self.cell_pairs[:,0]+self.dC.shape[1],self.cell_pairs[:,0]+2*self.dC.shape[1],
                               self.cell_pairs[:,1],self.cell_pairs[:,1]+self.dC.shape[1],self.cell_pairs[:,1]+2*self.dC.shape[1]])

        dE = np.bincount(rows, weights=np.ravel(dE,order='F'))
        return dE

    # Define constraints
    def constraints(self, x):
        # linear constraint
        Aeq = np.concatenate((np.zeros(x0.shape[0]*2), np.ones(x0.shape[0])))
        lincon = np.dot(Aeq, x)

        # nonlinear constraint
        x = x.reshape(x0.shape, order='F')
        q = x[:,0:2]
        p = x[:,2]

        b = np.matmul(self.dC, np.multiply(q.T,p).T)
        delta_p = np.matmul(self.dC, p)

        l1 = np.linalg.norm(b - (np.multiply(r1.T,delta_p).T), axis=1)
        l2 = np.linalg.norm(b - (np.multiply(r2.T,delta_p).T), axis=1)

        nonlincon = 0.5*(np.mean(l1) + np.mean(l2)) - scale

        cons = np.append(nonlincon, lincon)
        return cons

    def jacobian(self, x):
        # Jacobian of nonlinear constraints
        x = x.reshape(x0.shape, order='F')
        q = x[:,0:2]
        p = x[:,2]

        b = np.matmul(self.dC, np.multiply(q.T,p).T)
        delta_p = np.matmul(self.dC, p)

        t1 = np.divide((b - (np.multiply(r1.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r1.T,delta_p).T), axis=1)).T
        t2 = np.divide((b - (np.multiply(r2.T,delta_p).T)).T,np.linalg.norm(b - (np.multiply(r2.T,delta_p).T), axis=1)).T

        drX1 = rx_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,1],0], r1[:,0])
        drX2 = rx_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],0], q[self.cell_pairs[:,1],0], r2[:,0])
        drY1 = ry_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],1], r1[:,1])
        drY2 = ry_grad(p[self.cell_pairs[:,0]], p[self.cell_pairs[:,1]], q[self.cell_pairs[:,0],1], q[self.cell_pairs[:,1],1], r2[:,1])

        dE = 0.5 * (np.multiply(t1[:,0],drX1.T).T + np.multiply(t1[:,1],drY1.T).T) + \
             0.5 * (np.multiply(t2[:,0],drX2.T).T + np.multiply(t2[:,1],drY2.T).T)

        rows = np.concatenate([self.cell_pairs[:,0],self.cell_pairs[:,0]+self.dC.shape[1],self.cell_pairs[:,0]+2*self.dC.shape[1],
                               self.cell_pairs[:,1],self.cell_pairs[:,1]+self.dC.shape[1],self.cell_pairs[:,1]+2*self.dC.shape[1]])

        dE = np.bincount(rows, weights=np.ravel(dE/self.dC.shape[0],order='F'))

        # Add jacobian of linear constraints
        lin_dE = np.concatenate((np.zeros(2*x0.shape[0]), np.ones(x0.shape[0])))

        jac = np.concatenate([dE, lin_dE])
        return jac

initial_minimization = cyipopt.Problem(
    n = x0.size,
    m = 2,
    problem_obj = initial_problem(self.cell_pairs,self.dC),
    lb = np.concatenate((-np.inf*np.ones(x0.shape[0]), -np.inf*np.ones(x0.shape[0]), 0.001*np.ones(x0.shape[0]))),
    ub = np.concatenate((np.inf*np.ones(x0.shape[0]), np.inf*np.ones(x0.shape[0]), 1000*np.ones(x0.shape[0]))),
    cl = np.array([0, np.mean(p0)*p0.shape[0]]),
    cu = np.array([0, np.mean(p0)*p0.shape[0]])
)

initial_minimization.add_option('max_iter',2000)
initial_minimization.add_option('tol',1e-5)
initial_minimization.add_option('nlp_scaling_method','none')
initial_minimization.add_option('print_timing_statistics','yes')
initial_minimization.add_option('hessian_approximation','limited-memory')

res = initial_minimization.solve(x0.ravel(order='F'))
```
<p>Compared to the original matlab code:</p>
<pre class="lang-matlab prettyprint-override"><code>energyFunc = @(x) reducedEnergy( x, sparse(d0), t1, t2, r1, r2, bCells );
nonlinFunc = @(x) mycon( x, sparse(d0), Lscale, r1, r2, bCells );

lb = -inf*ones(size(x0));
lb(:,3) = .001;
ub = inf*ones(size(x0));
ub(:,3) = 1000;

Aeq = [zeros(1,size(d0,2)),zeros(1,size(d0,2)),ones(1,size(d0,2))];
beq = size(d0,2)*mean(p);

optimset = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point','MaxFunEvals',2e6,...
                        'MaxIter',2e3,'TolFun',1e-5,'GradObj','on','GradConstr','on',...
                        'HessianApproximation','lbfgs');

[x,ERes] = fmincon(energyFunc,x0,[],[],Aeq,beq,lb,ub,nonlinFunc,optimset);

function [ clin, ceq, dlin, deq ] = mycon( x, d0, Lscale, r1, r2, bCells )

    NB = size(d0,1);
    NC = size(d0,2);
    q = x(:,1:2);
    p = x(:,3);
    clin = [];

    B = d0*bsxfun(@times,q,p);
    deltaP = d0*p;

    B1 = B - bsxfun(@times,deltaP,r1);
    B2 = B - bsxfun(@times,deltaP,r2);

    L1 = sqrt(sum(B1.^2,2));
    B1 = bsxfun(@rdivide,B1,L1);

    L2 = sqrt(sum(B2.^2,2));
    B2 = bsxfun(@rdivide,B2,L2);

    ceq = .5*(mean(L1) + mean(L2)) - Lscale;
    if (nargout > 2)
        dlin = [];
        deq = .5*(bsxfun(@times,B1(:,1),drX1) + bsxfun(@times,B1(:,2),drY1)) + ...
              .5*(bsxfun(@times,B2(:,1),drX2) + bsxfun(@times,B2(:,2),drY2));
        rows = [bCells(:,1);bCells(:,1)+NC;bCells(:,1)+2*NC;bCells(:,2);bCells(:,2)+NC;bCells(:,2)+2*NC];
        deq = accumarray(rows,deq(:)/NB,[3*NC,1]);
    end
end

function [ E, dE ] = reducedEnergy( x, d0, t1, t2, r1, r2, bCells )

    q = x(:,1:2);
    p = x(:,3);

    B = d0*bsxfun(@times,q,p);
    deltaP = d0*p;

    B1 = B - bsxfun(@times,deltaP,r1);
    B2 = B - bsxfun(@times,deltaP,r2);

    L1 = sqrt(sum(B1.^2,2));
    B1 = bsxfun(@rdivide,B1,L1);
    L2 = sqrt(sum(B2.^2,2));
    B2 = bsxfun(@rdivide,B2,L2);

    E = .5*mean( dot(B1,t1,2).^2 + dot(B2,t2,2).^2 );

    if (nargout == 2)
        NB = size(d0,1);
        NC = size(d0,2);
        [ drX1, drY1, drX2, drY2 ] = fitDual.AFN.reducedLocalGrads( q, p, bCells, r1 ,r2 );

        IP1 = dot(B1,t1,2);
        IP2 = dot(B2,t2,2);

        dE = bsxfun(@times,IP1.*t1(:,1),drX1) + bsxfun(@times,IP1.*t1(:,2),drY1) + ...
             bsxfun(@times,IP2.*t2(:,1),drX2) + bsxfun(@times,IP2.*t2(:,2),drY2);
        dE = dE / NB;
        rows = [bCells(:,1);bCells(:,1)+NC;bCells(:,1)+2*NC;bCells(:,2);bCells(:,2)+NC;bCells(:,2)+2*NC];
        vals = dE(:);
        dE = accumarray(rows,vals',[3*NC,1]);

    end
end

<span class="math-container">```
```


问:
Change in Variables applied to biharmonic equation

Background

I want to solve the following biharmonic equation:

$$\frac{ \partial^4 s }{ {\partial \xi}^4 }+\frac{ \partial^4 s }{ {\partial \xi}^2{\partial t}^2 }+\frac{ \partial^4 s }{ {\partial t}^4 }=0\label{1}\tag{1}$$

on the unit square with the four Edges E1, E2, E3, and E4 (Fig. 1 left) and its prescribed boundary conditions (Fig. 1 right). The values of the boundary conditions
$\frac{ \partial \xi }{ {\partial s} }\bigg\rvert_{E1}$ and $\frac{ \partial \xi }{ {\partial s} }\bigg\rvert_{E2}$ can be set by me arbitrarily.

Additionally, $s$=$s$($\xi$,t) is known and is plotted in (Fig. 2).

Although the data in (fig.2) lies in the space s($\xi$,t), like the solution of Eq \eqref{1}, it's not the desired solution of Eq \eqref{1}. It is the starting point for further consideration, which becomes clear immediately.

Problem

Now I am going to quote the actual task, which I am struggling with:

... is solved by applying a change in variables. In fact we solve $\Delta\Delta$ $f$= 0 where $f$: $s$ $\in$ [0,1] $\mapsto$ [0,1] is a monotonic functions whichs maps a unit inverval onto a unit interval. The boundary conditions for $s$ are transferred to corresponding boundary conditions for $f$. After solving $\Delta\Delta$ $f$= 0, we find $f$ values at interior grid points and the corresponding $s$ values are found using $f^{-1}$. In practice, we define $f$: $s$ $\in$ [0,1] $\mapsto$ [0,1] such that $f$(0.5$\times$($s_{E3}(\xi)$+$s_{E4}(\xi)$)) $\equiv$ $\xi$

(Source: S. P. Spekreijse "Elliptic generation Systems", page 19; 1997) https://reports.nlr.nl/bitstream/handle/10921/1402/TP-1996-735.pdf?sequence=1)

Attempted solution

In the following, I want to show you how I understand the change in variables, what my approach is and which detail I don't grasp. If anything along the way is wrong, let me know, please! I am an engineer and my mathematical understanding is limited. Thus, some steps might be right, but it might also be the case, that everything is just wrong.

I am supposed to transform (Eq. 1) into:

$$\frac{ \partial^4 f }{ {\partial s}^4 }+\frac{ \partial^4 f }{ {\partial s}^2{\partial t}^2 }+\frac{ \partial^4 f }{ {\partial t}^4 }=0\label{2}\tag{2}$$

or

$$\frac{ \partial^4 \xi }{ {\partial s}^4 }+\frac{ \partial^4 \xi }{ {\partial s}^2{\partial t}^2 }+\frac{ \partial^4 \xi }{ {\partial t}^4 }=0\label{3}\tag{3}$$

(Eq. 2) and (Eq. 3) are the same because in the quoted task the relation $f$(0.5$\times$($s_{E3}(\xi)$+$s_{E4}(\xi)$)) $\equiv$ $\xi$ is defined, which implies $f$ $\equiv$ $\xi$. This change in variable is illustrated in (fig 3.).

In order to solve (Eq. 3) the boundary conditions in (fig. 1 right or rewritten in fig. 4 left ) are to be transformed from the space $s=s(\xi,t)$ to the space $\xi=\xi(s,t)$. This transformation is shown in (fig. 4).

Once this transformation has been completed, we have a boundary value problem with prescribed boundary conditions. By applying "fast Fourier transformation" or the "general coupled approach for biharmonic functions" this boundary value problem can be solved.

I have calculated the transformed values for $\xi$ from the the given relation $f$(0.5$\times$($s_{E3}(\xi)$+$s_{E4}(\xi)$)) $\equiv$ $\xi$.

Looking at at the boundary conditions in (fig. 4), I don't know how I can deduce the normal derivatives (Neumann Boundary Conditions) in (fig. 4 right), i.e.:

$$\frac{ \partial \xi }{ {\partial s} }\bigg\rvert_{E1}, \frac{ \partial \xi }{ {\partial s} }\bigg\rvert_{E2} \quad\mathrm{and}\quad \frac{ \partial \xi }{ {\partial t} }\bigg\rvert_{E3}, \frac{ \partial \xi }{ {\partial t} }\bigg\rvert_{E4}$$

from

$$\frac{ \partial s }{ {\partial \xi} }\bigg\rvert_{E1}, \frac{ \partial s }{ {\partial \xi} }\bigg\rvert_{E2} \quad\mathrm{and}\quad \frac{ \partial s }{ {\partial t} }\bigg\rvert_{E3}=0, \frac{ \partial s }{ {\partial t} }\bigg\rvert_{E4}=0$$

in (fig. 4 left).

I believe I can solve (Eq. 3) once I know the answer to this question.


问:
TVD slope / flux limiters formulation
Even if the formulation is the same the TVD slope limiter can be applied:

to state reconstruction at the interface, in 1D FV formulation, we reconstruct the $Q^*_{j+1/2}$ and the $Q^*_{j-1/2}$ in the formulation:
$$
Q_J^{n+1}=Q_J^n+\frac{\Delta t}{\Delta x}\left[F\left( Q^*_{j+1/2}\right)- F\left( Q^*_{j-1/2}\right)\right]
$$
where: $$Q^*_{j-1/2} = Q_{j-1}+\phi\left(Q_j-Q_{j-1}\right)  $$
$$Q^*_{j+1/2} = Q_{j+1}+\phi\left(Q_{j+1}-Q_{j}\right)  $$
(formulation took from this book)

For mixing up different numerical fluxes of different order. And hence approximating directly the fluxes as: $$ F \left( u_{i + \frac{1}{2}} \right) = f^{low}_{i + \frac{1}{2}}  - \phi\left( r_i \right) \left( f^{low}_{i + \frac{1}{2}}  - f^{high}_{i + \frac{1}{2}}  \right)$$
$$F \left( u_{i - \frac{1}{2}} \right) = f^{low}_{i - \frac{1}{2}}  - \phi\left( r_{i-1} \right) \left( f^{low}_{i - \frac{1}{2}}  - f^{high}_{i - \frac{1}{2}}  \right)  $$
(formulation took from this book or this book)

I see a strong difference that is in the nature of the numerical fluxes raising from the two cases. In the first case the numerical fluxes are built by a left and a right state raising from a linear approximation of the two state. In the second case I can have a purely Godunov method, and hence a purely piecewise approximation of left and right state, but the TVD flux helps in the treatment of the discontinuity and avoid the oscillations in proximity of that. Hence 2 questions:

Is it the only the TVD property in common with the two methods just cited?

Is it in case possible to combine the TVD reconstruction (first case) with TVD flux (second case) even with different slope/flux limiters?

Title: TVD slope / flux limiters formulation


问:
Solving huge dense square symmetric linear system
I have a linear system of the type

$A x = y$

where A is a dense, square, symmetric, positive definite matrix, $x$ a vector of unknown parameters, and $y$ is a vector of observed quantity.
I know that using a Cholesky decomposition based algorithm should be the fastest way to find the values of $x$ that fulfil the equation.
The only problem is I would like to do so for a matrix $A$ with dimension 500000x500000, and $x$ and $y$ 500000x1.
Is there any scientific computing library to do so? Which kind of hardware should I use?


问:
Slender body theory for plates? (hydrodynamics, Stokes flow regime)
A basic question on low Reynolds number hydrodynamics. When a slender beam is moved through a fluid, the dynamics may be calculated using slender body theory (e.g. Slender-body theory for slow viscous flow J. B. Keller and S. I. Rubinow, J. Fluid Mech. 1976; https://en.wikipedia.org/wiki/Slender-body_theory). In effect, it allows one to calculate drag coefficient on a material element of the beam and thus avoid having to solve for the dynamics of the fluid around the beam. Essentially, drag per unit length ends up being proportional to viscosity of the ambient fluid, and, obviously, the dimensions work out ( viscosity is in Pa*s, multiplied by velocity leaves N/m which is the required drag per unit length of the beam). My question is if one can do something similar for a 2D plate? It appears to not be possible, since, if nothing else, dimensions don't work out correctly. Is there a simple reason for why this approach fails for the case of a plate?</p>


问:
How to discretize PDEs using Chebyshev spectral method into a system of differential algebraic equations (DAEs)?
Let's take the heat equation. We have a time derivative and spatial derivative. How to discretize the spatial derivative using Chebyshev spectral method and convert it into DAEs? Like in the form of $ F(\dot{x}(t),x(t),t) $

$$ \frac{\partial u}{\partial t} = \frac{\partial^2u}{\partial x^2} $$

Assuming boundary conditions $u(0,t) = u(L,t) = 0$ and initial conditions $u(x,0) = \sin(x)$.

I tried to apply Chebyshev spectral method assuming $ \frac{\partial^2u}{\partial x^2} = 0 $ and the code to find u is given below.

```<code>N = 16; [D,x] = cheb(N); D2 = D^2;
D2N = D2(2:N,2:N);
f = sin(x(2:N));
u = D2N\f;
u = [0; u; 0];
</code>```

Does this mean that the values we get in u are in the form of x(t)? or my whole approach is wrong? How do I use the solver ode15s from this?


问:
Absolute value constraint in quadratic programming optimization
$$
argmin(x,y)=x^2+y^2+2y
$$

$$
s.t.\ \  y=|x-10|
$$

How can I convert the absolute value constraint to the constraint matrix (GX<=h, AX=b) in cvxopt?


问:
Finite difference solver for the 2D Poisson's equation with an integral boundary condition
I wanted to attempt an implementation of a finite-difference-based solver for the 2D elctrostatic Poisson equation when metallic objects are present. Also, I hope to take as input, the location of charges and their magnitude.

For the sake of simplicity, I consider a rectangular domain, inside which I attempt to solve the equation.

Something like so

(The red dots denote charges)

I have the following three points that I wanted to get some advice on:

1) Condition at points on metallic objects:

For the objects (which could be shells or filled), I intend to use the fact that metals are equipotential. To bring this into a matrix equation, what I thought was this: given an object, pick a point in/on it as a reference, and set the potentials of all other points to be equal to the potential of this point.

I am a little lost as to what condition I should use for the point chosen as reference since the Laplace equation does not hold for it. Currently, I am assuming that the Laplace equation holds at these reference points.
Is there an alternative way to convey that the object is equipotential?

2) Boundary condition:

I have studied some bit of surface integral equations; if I assume that the potential inside the domain is $\phi_{1}$, and outside, it is $\phi_{2}$, then the following boundary condition holds:

$$\oint_{S2}[\phi_{2}\nabla g(r,r')-g(r,r')\nabla\phi_{2}]\cdot\hat{n_{2}}dS=
-\oint_{S1}[\phi_{2}\nabla g(r,r')-g(r,r')\nabla\phi_{2}]\cdot\hat{n_{1}}dS$$ if $r'$ is outside the region enclosed by the two surfaces and the r is inside it or on the boundary.
The integration is done w.r.t to the variable $r$.

Here, $g(r,r')$ is Green's function for the 2D Poisson equation in free space, that is, it satisfies the equation

$$\nabla^{2}g(r,r') = \delta(\lvert r-r'\rvert)$$

For this problem, $g(r,r') = \frac{1}{2\pi}log(\lvert r-r'\rvert)+C$

(It is assumed that $\phi_{1}$ is zero outside the domain, and that $\phi_{2}$ is zero inside the domain. On the common boundary ($S1$), they have the same value)

$\hat{n_{1}}$ is the normal to the surface $S_{1}$ that points into the region outside the volume enclosed by $S_{1}$ and $S_{2}$.

However, I would appreciate any suggestions for alternative boundary conditions, since this integral makes my otherwise sparse matrix dense.

What I basically want to enforce through this integral is that there is no object or source or anything outside the boundary; i.e. the domain is an artificial limiting of the space under consideration.

Does anyone know of any boundary condition for electric potential that performs the same task, but is something like the Robin BC in electromagnetism?

3) Matrix solver:

Since this problem will generate a sufficiently large (but somewhat sparse) matrix, should I use a direct solver, or should I opt for an iterative solver? From some preliminary implementations I have done for this problem, the matrix is usually $10000 \times 10000$ or somewhere of that order; the integral condition makes the matrix dense.

I would greatly appreciate help on this, since I have spent a month on this problem. Also, I am averse to using FEM on this problem, since that is too cumbersome to implement.

In case anyone wants my code for the problem (I have written it in Fortran 2003), I can provide a Google Drive link for that.)


问:
Chemical advection of a fluid in a porous media
I am trying to solve an equation for chemical advection of a fluid in a porous media with this equation:

$\begin{align}
\frac{\partial (\phi(x,t) C(x,t))}{\partial t} = - \nabla . (\vec{q(\phi, t, x)} C(x,t))
\end{align}
$

with $\phi$ beeing the porosity, C the concentration of an element, and $\vec{q}$ the fluid flux.

I have another system of equations to compute $\phi$ and $\vec{q}$ at the next time step, so I only have 1 unknown for this equation: C.

What I first did was to simplify the system that way:

$\begin{align}
\frac{\partial C}{\partial t} = -   \frac{\vec{q}}{\phi} \nabla  ( C)
\end{align}
$

Which is a quasi-linear advection equation and easy to solve. The thing is that I would like to solve this using the conservative form without simplifying my system because mass conservation is important.

I would like to use a WENO scheme for the spatial discretization and a SSP Runge-kutta method for the time because I have sharp gradients in composition. I've already did it with the classical Burgers' equation, but I don't know what to do with $\phi$ in the time derivative.

Here are my 2 ideas:

either developing the time derivative that way:

$\begin{align}
\frac{\partial (\phi(x,t) C(x,t))}{\partial t} = \phi(x,t) \frac{\partial  C(x,t)}{\partial t} + C(x,t) \frac{\partial \phi(x,t)}{\partial t}
\end{align}
$

where I can approximate the second term because I know the value of $\phi$ at the next timestep. My gut feeling tells me this is a bad idea.

My second idea is to make a change of variable:

$\begin{align}
\phi C = \tilde{C}
\end{align}
$

knowing that the flux is equal to this, with $\vec{v}$ the velocity of the fluid:
$
\vec{q} = \vec{v} \times \phi
$

My equation then would becomes:

$\begin{align}
\frac{\partial (\tilde{C}(x,t))}{\partial t} = - \nabla . (\vec{v} \tilde{C})
\end{align}
$

That I should know how to solve. Then I could retrieve $C$ by dividing $\tilde{C}$ by $\phi$ at the next timestep.

Would that works? Or do you have other suggestions?


问:
"Unspoken Tradeoff" in Convergence Rates for "Quasi-Newton vs. Gradient Descent"
Is there an "Unspoken Tradeoff" in Convergence Rates for "Quasi-Newton Methods vs. Gradient Descent"?

As a quick summary:

Gradient Descent based algorithms try to find the minimum of the function they are optimizing by repeatedly using information obtained from the First Derivatives

Quasi Newton methods (e.g. BFGS) tries to find the minimum of the function by repeatedly using/approximating information obtained from the Second Derivatives

Thus, it is almost a logical conclusion to believe that for the exact same (convex) function, Gradient Descent would on average require fewer calculations per iteration when compared to Quasi Newton Methods.

But in the spirit of "tradeoffs" - are there any theoretical results that show Quasi-Newton Methods gain stronger convergence rates compared to Gradient Descent based algorithms ?

For instance, since information about the second derivative gives the optimization algorithm additional knowledge about the function's local curvature - in theory, we would expect this additional knowledge would "better lead" the optimization algorithm to the minimum when compared to information only about the first derivative (especially in higher dimensional functions). Or is this additional second derivative information gained by Quasi-Newton Methods generally negligible compared to the first derivative information, and as a result, there is no general "strength in convergence" obtained by Quasi-Newton Methods when compared to Gradient Descent?

Thus: Are there any theoretical results that compare the "strength of convergence" for Quasi-Newton Methods vs. Gradient Descent Based Algorithms?


问:
Solving linear system and obtaining operator norm
I need to solve a linear system of the form $(\mathrm{Id} + \mathbf{J})\mathbf{x} = \mathbf{b}$ for $\mathbf{x}$ and I also need to compute the operator norm of $\mathbf{J}$ (i.e. the largest singular value). The matrix $\mathbf{J}$ has no special structure other than being real-valued; the vector $\mathbf{b}$ is also real-valued.

I'm currently doing this in two steps, treating the two problems as essentially independent. I'm wondering what is the best numerical method for this problem.


问:
Generate numerically stable equation automatically
I remember there is a software (and its website interface) that can generate numerically stable equation from an expression automatically.

To see the problem it tries to solve, let’s take a look at a simple example of $f(x) = \log(1 + \exp(x))$.
If $x$ is a relatively large number (e.g. 800), $\exp(x)$ could be very large and reach inf in the floating point representation.
So it is better to split how to compute the function: (1) if $x < a$ (say $a = 50$), then it computes $f(x) = \log(1 + \exp(x))$ directly, (2) if $x \geq a$, then it computes $f(x) = x$.
The software I mention in the previous paragraph can automatically generate such equation splitting.

My question is: how does it (the software) work? And if you can put the reference to the website, it would be even better.
Thanks!


问:
Convergence of Conjugate Gradient Method
Can anyone help me prove that, if the residual $r_0 = b-Au_0$ can be written as a linear combination of $k$ eigenvectors of $A$, then the CG method will converge to a solution $x^{*}= A^{-1}b$ in at most k steps. With $N=I_n$ and $A$ is symmetric and positive definite.


问:
Solving 1D Advection Equation Using Midpoint-Rule and Finite-Differences
I want to solve the advection equation ($v_0 \in \mathbb R$)
$$
\frac{\partial f}{\partial t} + v_0 \frac{\partial f}{\partial x} = 0
$$

using 2nd order Runge Kutta like the midpoint rule for the time integration and an upwind (backwards finite differences) for the spatial derivation.

Here is a general Runge-Kutta implementation in python

```<code>def rk2a( f, x0, t ):
    n = len( t )
    x = np.array( [ x0 ] * n )
    for i in range( n - 1 ):
        h = t[i+1] - t[i]
        k1 = h * f( x[i], t[i] ) / 2.0
        x[i+1] = x[i] + h * f( x[i] + k1, t[i] + h / 2.0 )

    return x
</code>```

And that makes sense and if I'd have to solve something like $y'=f=x$ I'd know how to pass $f$ but I'm confused because now we have $f=-v_0\frac{\partial f}{\partial x}$.

Since I want to use an upwind scheme for the spatial derivative, we have

$$
\frac{\partial f}{\partial x}|_{i-1/2} = \frac{q^{n}_i - q^{n}_{i-1}}{x_i - x_{i-1}}
$$

whereas $n$ denotes the $n$-th timestep and $i$ denotes the $i$-th grid point.

Now to implement that, I'd have to write down f( x[i], t[i] ) and f( x[i] + k1, t[i] + h / 2.0 ) resp. provide a general f.

Now while I think I can write down f( x[i], t[i] ), I have no idea how I'd write down f( x[i] + k1, t[i] + h / 2.0 ).

So the question is: How do I use a upwind scheme in the midpoint rule?

Edit: I assume equispaced grid.


问:
Branch and Bound vs Evolutionary Algorithms
On an informal level, I have heard the following statements being made:

The "Branch and Bound" Algorithm has the ability to provide exact solutions to optimization algorithms, but can take a long time to run.

Evolutionary Algorithms can only provide approximate answers, but take less time to run.

My Question: Are there any references which quantify this above statement using mathematical formulas? Would any of these references be able to provide detailed comparisons regarding the theoretical behavior of both types of algorithms (e.g. strength of convergence, varying levels of convexity, etc.)?

Thanks!


问:
3D web visualization with arbitrary cutting plane
How can I visualize 3D data volumes with arbitrary cutting plane on web pages? I could not find a library that allows so.</p>


问:
Transparent Boundary Conditions for Finite Difference ADI PR 2D TDSE solution
I want to put (non-dirichlet) boundary conditions inside the code I wrote to solve the 2dim TDSE using the alternating direction implicit Peaceman - Rachford method.
$$
(1 + iB\Delta t/2 ) \psi^{n+1/2} = (1 - iA\Delta t/2) \psi^n
$$
$$
(1 + iA\Delta t/2 ) \psi^{n+1} = (1 - iB\Delta t/2) \psi^{n+1/2}
$$, where $A$ contains operators only on the coordinate $x$ and $B$ contains operators only on the coordinate $y$ and $\psi = \psi(x,y,t)$.

This PR method, as applied to my particular case, is detailed in https://www.sciencedirect.com/science/article/pii/S0021999196955886 and in the Appendix of http://www.iap.tu-darmstadt.de/tqp/papers/BauMul99.pdf.

I have come across http://rrp.infim.ro/2003_55_4/carjan.pdf, from now on called Carjan 2003 which details a particular approach to implement Transparent Boundary Conditions for the 1 dim TDSE.
Because the ADI basically solves a one-dim problem at each half timestep (from $t_n$ to $t_{n+1/2}$ and from $t_{n+1/2}$ to $t_n$) via Crank Nicolson in 1D, I believe I can implement BCs from 1dim TDSE studies. Am I correct?

If I am correct and I can implement them to my 2D problem, I am having trouble understanding how to treat the left boundary of a 1dimensional line (x-coordinate). Carjan 2003 only explains how to treat the right boundary (page 560, or 6/25 from the pdf). I have also went to the source of these Transparent Boundary Conditions described in Carjan 2003: the original (quoted as [11, 12] in the above link) is a treatment of the propagation of a laser beam in space [nothing related to quantum mechanics, but there is a nice linkage between the two, made in the publication from link above]. Still, nothing related to the left boundary, only the $x = x_{max}$ boundary (RHS one).

In the paper, the right boundary is treated by firstly assuming that the waves near it are of the form $\psi_0 \exp(i k_x x)$, where $i^2=-1$, and $k_x$ and $\psi_0$ complex constants to be determined. Then the continuity equation for the wavefunction $\psi$ is applied (1-dim problem) and integrated across the simulation domain:
$$
\rho = \int_{a}^{b} \vert \psi \vert^2 dx
$$
$$
\dfrac{\partial \rho}{\partial t} = \dfrac{i \hbar}{2 m} \left[ \psi^* \dfrac{\partial \psi}{\partial x} - \psi \dfrac{\partial \psi^*}{\partial x} \right] \vert_{a}^{b} \equiv -F_b + F_a
$$

They say that $F_b$ signifies the "energy" flux leaving through the right boundary and $F_a$ is the "energy" flux entering through the left boundary. It is then wanted to obtain the $k_x$ such that $-F_b$ is positive. They say that this implies that the contribution to the overall change in energy from this RHS boundary will be always negative (radiative energy can only flow out of the problem region). They obtain a $k_x$ s.t. respect the above, in particular the condition is: real part of $k_x$ shall be positive to ensure that radiation only flows outward at the RHS boundary.

I have looked here https://www.youtube.com/watch?v=Ls5HS2MLXpg and I understand that the continuity equation arises due to a fluid flowing from the LHS into a volume (my line domain) and going out on the right (for no creation/destruction of the fluid).
I can also understand that the derivation can be applied the other way around: the fluid coming into the volume from the right and leaving through the left boundary.

My 2D TDSE problem, my aim, for which I want to implement BC's, has a wavefunction $\psi(x,y,t)$ which spreads from its initial condition (localized at a central point) across the whole simulation domain as time passes. It delocalizes over time. I need to remove (absorb, not reflect, as much as possible) the wavefunction which going out of the boundaries of the simulation box. So the $x_{min}$ and $x_{max}$ boundaries shall both have such transparent boundary conditions. They both shall remove any incoming wavefunction.

My question is: in the derivation for the left boundary, how shall I proceed to do it? I am particularly concerned about the starting assumption: in the paper, for the RHS boundary, they said: waves near the boundary are of the form $\psi_0 \exp(i k_x x)$. Now, when treating the LSH boundary, shouldn't the waves be of the form: $\psi_0 \exp(- i k_x x)$. [Tangent question: where has the time dependence disappeared in their ansatz?] Then the continuity equation, after it's been integrated as (see the boundaries, changed):
$$
\rho = \int_{b}^{a} \vert \psi \vert^2 dx
$$
shouldn't it read:
$$
\dfrac{\partial \rho}{\partial t} = \dfrac{i \hbar}{2 m} \left[ \psi^* \dfrac{\partial \psi}{\partial x} - \psi \dfrac{\partial \psi^*}{\partial x} \right] \vert_{b}^{a} \equiv -F_a + F_b
$$.
But if I assumed $exp(-i k_x x)$ for the waves near the LHS boundary, the derivation, from now on, is not similar to the derivation for the RHS boundary presented in the papers. Again, I need to make $F_a$ > 0, but this time i have that -1 in the exponential. What I want to emphasize is that I believe the derivation shall be more or less identical w.r.t. which boundary we consider, but I don't get this.
In the end, eq. (3.7) from Carjan 2003 shows the result of putting a transparent boundary condition at the RHS boundary: it changes the last equation from the Crank-Nicolson scheme (i.e. the one for the node prior to the very last one, i.e. the node before the boundary one). When I work with the LHS boundary, the first equation from the Crank-Nicolson scheme will have to change. Thus my concerns because a mistake in the derivation will ruin everything.

Thank you!


问:
Geodesic approximation algorithms for minimal geodesic curvature

Introduction

I am building an application in which, given a surface and a pair of points on it, an analytic expression of some geodesic arc between them is needed, preferably the one with the shortest length. Furthermore, this parameterized curve, needs to have minimal geodesic curvature.

Is there a suggested tool that specializes in minimizing the geodesic curvature of the generated approximation?

I am looking for geodesic curvature less than $10^{-8}$, ideally on any given NURBS representation of a smooth surface.

What I have tried so far:

The geodesic element in Grasshopper. This produces a geodesic curve with around roughly $10^{-1}$ geodesic curvature. (reference surface is a cone)

The shortpath algorithm implemented in Grasshopper. This produces a geodesic curve with geodesic curvature roughly of the order $10^{-3}$ to $10^{-5}$. (reference surface is a cone)

A remark:
When I assign a value to the geodesic curvature of a curve (say $10^{-5}$), I am referring to the absolute maximum of a dense sampling of the geodesic curvature along said curve.


问:
How is ADMM Separable?
I'm learning about ADMM by reading Boyd's paper Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.

The paper says that ADMM is an improvement over the Method of Multipliers because the latter requires the minimization of the augmented Lagrangian:

$L_\rho(x,y) = f(x) + y^T (Ax − b)+(\rho/2)\lVert Ax − b\rVert^2_2$

which we cannot minimize in parallel because the added quadratic penalty term makes the function non-separable.

However, the paper says that ADMM requires minimization of the augmented Lagrangian functions:

$x^{k+1} = \underset{x}{\operatorname{argmin}} L_{\rho}(x,z^k,y^k)$

$z^{k+1} = \underset{z}{\operatorname{argmin}} L_{\rho}(x^{k+1},z^k,y^k)$

where the augmented Lagrangian is defined as

$L_{\rho}(x,z,y) = f(x) + g(z) + y^T(Ax + Bz − c)+(\rho/2)\lVert Ax + Bz − c\rVert^2_2$

The paper says that ADMM is an improvement because it has decomposability unlike the Method of Multipliers. However, shouldn't the augmented Lagrangian for ADMM also be non-separable because it also contains a quadratic penalty term $(\rho/2)\lVert Ax + Bz − c\rVert^2_2$ at the end? What allows it to be more decomposable?

Title: How is ADMM Separable?


问:
Well-conditioned pseudospectral for computing eigenvalues to (partial) differential equations
I am working on writing a Chebyshev pseudospectral method (see for example "Chebsyhev and Fourier spectral methods" by John Boyd) to solve for the eigenvalues of differential equations of the form

$ \left(a_2(x) \frac{d^2}{dx^2} + a_1(x)\frac{d}{dx} + a_0(x)\right)u - \lambda u = 0$

There are "regular" boundary conditions (i.e. $u$ should be regular on the boundary; see the edits below). Using Chebyshev pseudospectral methods I can turn the above differential equation into a matrix eigenvalue equation

$ Mu=\lambda u$

My problem is that the Chebyshev pseudospectral differentiation matrices are very poorly conditioned, e.g. the condition number of the pseudospectral matrix to $d^2/dx^2$ scales as $N^4$, where $N$ is the number of collocation points. This makes computing the eigenvalues to the operator $M$ numerically unstable when I use many collocation points, and I need to use many points to properly resolve rapidly oscillating eigenfunctions $u$.

I could apply a preconditioner matrix $P$ to $M$ to get a generalized eigenvalue equation

$PMu - \lambda Pu=0$

but I'm not sure that would actually make the eigenvalue solver more stable. Does anyone know applying a preconditioning matrix could make computing the eigenvalues more numerically stable? And if so, could they post several references that shows this is true/provides some examples?

EDITS:

I'd like to consider equations where $a_2(x)$ could vanish at the boundaries (see third bullet point), so there are two choices of boundary conditions: "regular" and "irregular" boundary condtions. What is nice about using pseudospectral methods is that they naturally pick the "regular" solution. All of which is to say that I want to use pseudospectral methods to solve this ODE.

The matrix M is not generally self-adjoint or symmetric for Chebyshev pseudospectral methods, so I'm not worried about losing that property in the system

The functions $a_i(x)$ that I'm looking at are fairly complicated, but roughly they take the form

$a_2(x) = x^2(1-x)$

$a_1(x) = c_0 + c_1x + c_2x^2 + c_3x^3$

$a_0(x) = d_0 + d_1x + d_2x^2$

The coefficients $c_i,d_i$ are generally complex (not real).

I don't need all of the eigenfunctions, but I would like to have more than the first few lowest one. To make things worse, I want to have high accuracy (say a relative error of $10^{-10}$ or so), so that is driving the need to have many Chebyshev coefficients even for the lower eigenmodes.


问:
Problem in parameterizing a CVXPY program
I am trying to parameterize a CVXPY program as I need to repeatedly solve the problem, but I noticed that when my parameters are complex numbers, CVXPY models the problem in each iteration. For example, suppose I want to solve the following problem:

$$
\underset{\mathbf x}{\text{minimize}} | \mathbf A \mathbf x - \mathbf b| +  |\boldsymbol{\gamma} \mathbf x|,
$$
where $\mathbf A \in \mathbb C^{3\times 4}$, $\mathbf x \in \mathbb C^{4 \times 1}$, $\mathbf b \in \mathbb C^{3 \times 1}$, and $\boldsymbol{\gamma} \in \mathbb C^{1 \times 4}$.

I don't need to solve this problem repeatedly, but I am forcefully doing so to check whether I can consider $\boldsymbol \gamma$ as a parameter or not?

Below is the code:

```<code>import cvxpy as cp
import numpy as np

m = 4
n = 3
A = np.random.randn(n,m)+1j*np.random.randn(n,m)
b = np.random.randn(n)+1j*np.random.randn(n)

x = cp.Variable((m,),complex=True)
gamma = cp.Parameter((1,m),complex=True)
obj = cp.norm(A@x-b)+cp.abs(gamma@x)
prob = cp.Problem(cp.Minimize(obj))

for iIter in range(5):
    gamma.value = np.random.randn(1,m)+1j*np.random.randn(1,m)
    prob.solve(verbose=True)

</code>```

What I see from verbose is that CVXPY first models my problem in each iteration and then solves it. On the other hand, if I assume $\boldsymbol \gamma \in \mathbb R^{1\times 4}$, it bypasses the modeling from the next iteration. Why is this so? Are we not allowed to create complex-valued parameters?


问:
Solving the non-linear Hamiltonian using Scipy's root finding method
I am a complete novice to computational physics and am finding difficulty in implementing a code to iteratively solve for a $2\times2$ nonlinear Hamiltonian using Scipy's root solver. I can't seem to find open-source codes that I can take inspiration from and learn by myself. Here is the problem:

Consider a single unit cell consisting of two wavefunctions (corresponding to a dimer). I need to solve for the nonlinear Hamiltonian associated with this by giving initial guesses for my wavefunction ($x$ values).

If anyone here can help me explain, and/or provide me with a link to an example code that solved something similar/using similar methods, I will appreciate it tremendously.


问:
Generate polynomial basis through a sequence of SVD
I need help to understand how to use the result given by an algorithm for constructing an orthonormal polynomial basis over $L^{2}(X)$, where $X\subset\mathbb{R}^2$, with respect to the inner product $(f,g) = \iint_X f(x,y)g(x,y)\,dxdy$. We denote $\mathcal{P}_{n} = \{\text{Polynomials }P(x,y): \text{ total degree of } P \leq n\}\subset L^{2}(X)$. The algorithm is as follows.

(1) Consider the vector corresponding to the constant polynomial $1$. Normalize it. The resulting vector provides an orthonormal basis of $\mathcal{P}_{0}$.

(2) Suppose we have constructed an orthonormal basis of $\mathcal{P}_{k}$, consisting of $K$ vectors. Multiply each of the vectors in this basis by $x$ and $y$, respectively, and put the resulting vectors, together with the vectors in the basis of $\mathcal{P}_{k}$, in a matrix having $3K$ columns. Calculate the SVD of the matrix. The SVD yields an orthonormal basis of the column space of the matrix, which is also an orthonormal basis of $\mathcal{P}_{k+1}$.

(3) Apply the procedure described in step (2) for each $k$ in $\{0, 1, \ldots , n − 1\}$ in order to obtain an orthonormal basis of $\mathcal{P}_n$.

My thoughts: One need to implement an SVD that uses projection with respect to the given inner product. Also, I assume that a quadrature rule is given to compute the inner product.

Questions:
(1) What does it mean to multiply with x and y in step (2)? The SVD is not performed symbolically. Should I sample $X$ and multiply by those values?

(2) Executing the algorithm above, how do I use the result to evaluate the resulting orthonormal polynomial basis at a point in $X$? How do I take the derivative of this basis?


问:
solve Ax=b for outrigger A matrix python
I implement Crank-Nicolson 2D finite-difference method.

I get a matrix A which is banded with 1 band above and below the main diagonal, but also contains 2 additional bands , further apart from the main diagonal, so it is NOT penta-diagonal.

A picture showing the structure is below. My matrix is the RHS one. The LHS is easy, it's the penta-diagonal one.

I couldn't find up until now a way to solve Ax = b with A being the RHS matrix from the photo in python.
I could barely find a name for it, in these lecture notes https://ocw.mit.edu/ans7870/2/2.086/F12/MIT2_086F12_notes_unit5.pdf it is called an 'outrigger' matrix (page 403).

At the moment I am using spsolve from from scipy.sparse.linalg, into which I feed two arguments, namely sparse.csc_matrix(A) and sparse.csc_array(b), where A and b have been defined initially as A = sparse.dok_matrix((size, size), dtype=np.complex64) and b = sparse.dok_array((size, 1), dtype=np.complex64), then populated with values by iterating element by element through them.
It is extremely slow and I was wondering maybe someone more experienced knows a way to exploit the structure appearing in A.

Thank you!


问:
Help with CVXPY and Disciplined Convex Programming
I'm trying to recreate Figure 1 in this paper. This requires maximizing equation (19), which I have convinced myself is concave, but I am having trouble implementing it in CVXPY. Here is the code I think is relevant to this question:

```<code>import cvxpy as cp
import numpy as np

def tanh_atom(x):
    return (cp.exp(2 * x) - 1) / (cp.exp(2 * x) + 1)

def functional(X, T, dt, G=None):
    if G is None:
        G = np.identity(X.shape[0]-1)
    Q = np.linalg.inv(np.matmul(G, G.T))
    dXk_dTk = cp.diff(X, axis=0) / dt
    Fk = tanh_atom(X)[:-1]
    diff = dXk_dTk - Fk
    S = dt * cp.quad_form(diff, Q)
    return -.5 * S

def norm_pdf_atom(x, mu, sd):
    coef = 1 / (sd * (2*np.pi)**.5)
    z = (x - mu) / sd
    return coef * cp.exp(-.5 * cp.power(z, 2))

def merit(obs, X, T, dt):
    term1 = functional(X, T, dt)
    term2 = norm_pdf_atom(X[0][0], 0, .16**.5)
    term3 = norm_pdf_atom(obs, X[-1][0], .16**.5)
    return term1+term2+term3

t0 = 0
t_end = 5
N = 5
obs = 1.5
T = np.reshape(np.linspace(t0, t_end, N), (N, 1))
dt = T[:,0][1] - T[:,0][0]
dt = cp.Parameter(value=dt, nonneg=True)
T = cp.Parameter(np.shape(T), value=T, nonneg=True)
X = cp.Variable(np.shape(T))
constraints = [X >= 0, X <= 10]
prob = cp.Problem(cp.Maximize(merit(obs,X,T,dt)), constraints)
prob.solve(verbose=True)
</code>```

Two issues

CVXPY tells me that term1 is not DCP siting this expression:

```<code>(exp(Promote(2.0, (5, 1)) @ var2) + Promote(-1.0, (5, 1))) /
(exp(Promote(2.0, (5, 1)) @ var2) + Promote(1.0, (5, 1)))
</code>```

For both term2 and term3 CVXPY tells me that the expression below is not DCP. It actually suggests using DQCP with qcp=True in the solve() call but this also gives an error. Both these terms are just maximizing a Gaussian distribution so I really feel like this shouldn't be an issue. Is there a better way to do this other than creating the pdf function?

```<code>exp(-0.5 @ power((var2[0, 0][0] + -0.0) / 0.4, 2.0))
</code>```

In both scenarios I can't figure out why they are not DCP, as they both are concave, unless I messed up my reasoning for that. I've tried using as many CVXPY atoms as possible but still don't know what's going on. Any help would be much appreciated.

Note: I individually tested each term in merit just by commenting out the other terms and running the script to see if there were errors. I have also recreated this situation using numpy by changing every place in the script with cp to np. Then I would randomly generate column vectors X and feed them to merit and plot the best X. Although this is a very crude way to find the optimal X, it gave me results close enough to the desired figure to make me think the math behind the code is correct.

Edit:
I found some errors in my original question so I have made changes as necessary.


问:
Expected residual at double root
In a generic rootfinding problem $f(x) = 0$, we assume that the probability that the root $x$ is a floating point representable is zero. Hence, the best floating point approximation $\hat{x}$ to $x$ gives $f(\hat{x}) = f(x(1+\delta)) = x\delta f'(x) + \mathcal{O}(\delta^2) \approx x\delta f'(x)$ where $|\delta| \le \mu_{M}$, where $\mu_{M}$ is the unit roundoff.

Hence, any rootfinding iteration can be terminated when $|f(\hat{x})| \le \mu_{M}|x f'(x)|$. However, this termination criteria fails in the presence of a double root.

What is a sensible termination criteria for a rootfinding problem in the presence of a double root? (The assumption being that $f \in C^{2}$ and $f''$ is available.)

@LutzLehman makes a valid point that this may be impossible in general. Nonetheless I have succeeded in every case I have tested for the case of polynomials by use of the expected residual

$$
\mu_{M}\left[\sum_{k=0}^{n} |c_k||x|^k + |xp'(x)|  \right]
$$

Maybe using a running error estimate in the evaluation of $f$ might succeed in the general case?

Title: Expected residual at double root


问:
Is it possible to impose periodic boundary conditions on Chebyshev differentiation matrix?
This is a repost from another community, it was recommended that I post it here instead.

I have been using the Chebyshev collocation method for solving PDEs. In particular to solve a 3D generalized Poisson's equation. One of the problems I am facing is that one of my dimensions is periodic and would require periodic BCs. I have been following with the Trefethen's book "Spectral Methods in Matlab", but so far do not see anything like that. I am wondering if it's not possible even though we can apply FFTs to these derivative matrices to speed up things.

A little background, I am solving for 3D fluid equations to study ionospheric turbulence, basically I am writing a potential solver to solve for a current closure (charge conservation) and this is my equation:
$$\textbf{a}\nabla_{\perp} \cdot \left(n \nabla_{\perp} \phi \right) + \textbf{b} \frac{\partial}{\partial z} \left( n \frac{\partial \phi}{\partial z}\right) = \textbf{S} $$

Originally, I was able to solve this above equation in 2D using spectral methods with Fourier basis function with periodic BCs, but with the addition of the third non-periodic dimension I thought of rewriting everything using Chebyshev. Any suggestions or examples using periodic BCs with chebyshev matrices would be appreciated. Thanks.


问:
Simultaneously partition all mesh elements (nodes, edges and cells) among processors
The numerical method that I'm implementing stores various fields in various elements of an unstructured triangular mesh. Data may be stored in cells, nodes and even edges of the mesh. For a parallel implementation I need to distribute all the mesh elements among processors.

There are well known approaches to the problem of distribution cells or nodes, but I failed to find how to distribute all the elements at once.

Of course one can distribute all mesh elements independently by partitioning the corresponding adjacency graphs. But this approach gives no guarantee that the cells of i-th processor and the nodes of the same processor are close (that is cells of i-th processor are mostly composed of the nodes from the same processor).

So I came up with an approach to partition only mesh nodes and then try to partition edges and cells optimally, i.e. minimize cross-processor edges and cells while maintaining a fair distribution. Anyway, a proper implementation is quite tough and has some corner cases.

Do you know any existing solutions to the problem above?


问:
Passing boundary conditions to solver
Quite broad question,
Currently building own Poisson solver subroutine for CFD solver.

Works smooth, the goal is to generalise the input and make it flexible.

Description of also:

Memory allocation.

Filling source right hand side $-\Delta u(x,y) = f(x,y)$ function of Poisson equation.

Filling sparse matrix with inner nodes with coefficients of discretisation (1/dx2, 2/dx2, 1dy2, 2/dy2)

Entering into matrix BC's (here I do hardcoding - L R T B separately and corners).

Converting to Sparse CSR format.

Calling Conjugate Gradient solver function.

Writing output to a file.

Ending, freeing memory.

If I figure our how to generalise input will likely get rid of the step 5 and input directly into the dense matrix of CSR format. Because initially I used two matrices

```<code>matrixOfValues
matrixOfLocations
</code>```

Both are matrices of size $(N_x N_y)\cdot 5$ of int/double (location/value),
5 is for the node with 4 neighbours in x and y directions.

The purpose of the intermediate storage was to not mess up, however most of the time is consumed by ConjGrad solver anyway.

Any recommendations on passing the BC's to the solver more generally without hardcoding?
On top of that, there are three different types of BC, namely: Dirichlet, Neumann and Robin.

Feels like throwing math function as an argument is not a good idea, I could pass an array (say Western wall nodes as a start, then move counterclockwise, and passing all of the boundary node coefficients, i.e. a triangle with neighbours and, perhaps, separate for corners), but why need to invent the bicycle if this should've been done earlier in the most effective way.

Any relevant literature you can recommend?

Thank you!

Title: Passing boundary conditions to solver


问:
Semi-lagrangian method for compressible fluid (non divergence free)
I am looking for a semi-Lagrangian method for advection with a non divergence free velocity field.

The equation is

\begin{align}
\frac{\partial C(x,t)}{\partial t} &= - \nabla \cdot (\vec{v}(x,t) C(x,t)) \\
C(x, 0) &= C_0(x)
\end{align}

I've solved that using the formulation in this paper [1], which describes a quasi-monotone and conservative scheme, where they use a property of the Jacobian determinant to come up with

\begin{align}
\frac{\partial J(x,t)}{\partial t} &= J \nabla \cdot \vec{v}(x,t) \\
J(x, 0) &= 1
\end{align}
where $J$ is the determinant of the Jacobian.

Then they replace the first 2 equations with

\begin{align}
\frac{D J(x,t)}{D t} &= 0 \\
C(x, 0) &= C_0(x)
\end{align}

I have managed to implement it and it works fine. My worries is that the paper is a bit old (2002) and I would like to know if there is another way to achieve that with a semi-Lagrangian scheme? I think the main drawback of this approach is that conservation is achieved globally and not locally, and that the use of the Jacobian prevents the use of dimensional splitting.

I've not found easier alternatives in the literature I have read. I have read a lot about the combination with WENO/ENO scheme that seems to be promising, but I didn't find a case for non divergence free velocity field.

Is there some more up to date techniques than the one I used?


Title: Generalized Eigenvalue Problem using MATLAB

Body:
I'm trying to solve a generalized eigenvalue problem. I have two matrices $H$ and $S$ such that:

$$
HX=λSX
$$

I need to find the eigenvalues $\lambda$. The matrices $H$ and $S$ are real, asymmetric, contain positive and negative decimal numbers, and (can be) singular. I am using MATLAB and have tried two different approaches:

I used eig(H,S) command which uses QZ Algorithm to solve for eigenvalues.
I formed an equation in $\lambda$ by using $|H-\lambda S|=0$ and then solved for $\lambda$.
Here are the matrices: (both $H$ and $S$ are $24\times 24$)

```<code>H=[0,0,0,0,192,1917.04064,10332.51505,40092.51227,125681.1486,338350.2206,811892.8294,1779728.921,3625982.355,6953387.916,12670976.81,22100000,37132930.27,60353006.25,95276316.19,146559937.4,220274060.5,324208411.9,468219308,664618721.8;
0,0,0,0,192,1893.475124,10051.90014,38308.22391,117609.6433,309187.9535,722364.2569,1537115.973,3030677.025,5606681.841,9824567.083,16426180.74,26355891.41,40770017.02,61031174.08,88683247.04,125403139.6,172926323.5,232944447.3,306974887;
0,0,0,0,192,1846.924351,9507.779872,34923.83828,102685.9158,256815.0258,566753.8301,1130503.089,2072244.588,3531888.862,5644710.273,8510616.36,12154554.65,16481872.51,21234784.98,25958063.51,29983201.26,32440241.94,32304892.86,28485362.66;
0,0,0,0,192,1778.534558,8732.953666,30281.20319,83090.50097,191433.6927,383369.9941,681577.7806,1089000.57,1571702.814,2043718.6,2360415.659,2327150.205,1728131.214,376286.8361,-1820993.208,-4791939.112,-8243585.409,-11638086.16,-14220529.21;
0,0,0,0,192,1689.989727,7773.203992,24831.36406,61515.88289,124689.9821,212051.9399,303583.89,356337.2627,308129.1939,94242.56924,-323226.8121,-920239.2345,-1583058.32,-2104105.201,-2214813.175,-1659107.74,-295993.1841,1796496.926,4254923.427;
0,0,0,0,192,1583.470126,6683.555129,19072.50505,40639.69537,66711.57632,81796.22718,60590.29663,-21861.2573,-171429.7577,-355062.5907,-492927.1661,-475398.2758,-208930.4549,320950.1747,997476.1882,1566265.895,1699768.472,1129529.762,-194208.6509;
0,0,0,0,192,1461.598622,5523.845413,13484.72812,22648.14037,23924.65349,4088.400074,-44523.67789,-110266.4937,-154796.0153,-124116.6074,18352.09343,250026.0416,467393.1673,513093.6957,258167.4211,-292568.7362,-936296.041,-1313188.428,-1070939.742;
0,0,0,0,192,1327.376095,4354.022806,8472.235861,8915.430969,-1775.330646,-26379.77787,-53287.95266,-55123.02039,-5131.095738,92963.75123,184430.9087,180646.3185,21137.78348,-252747.0452,-475735.34,-436097.8969,-33470.52706,583046.8515,1024123.527;
0,0,0,0,192,1184.107549,3229.595901,4321.651839,-104.2768198,-12455.61648,-25344.84529,-20852.91355,14642.78081,67162.133,88643.95935,29331.92183,-103721.9759,-215885.8219,-175668.2952,63339.97258,366040.1121,463007.7228,161145.8379,-426631.9431;
0,0,0,0,192,1035.320732,2197.653181,1181.837405,-4776.384006,-12614.95411,-10733.55999,10443.54925,39434.88701,40860.0242,-11919.57845,-90391.52336,-109544.5177,-4764.383231,169447.4417,238928.2291,63214.86651,-273067.0285,-452649.3406,-197417.4769;
0,0,0,0,192,884.6792681,1293.804398,-933.634351,-6033.445815,-7367.269353,3768.214436,22179.4179,22745.61092,-13409.71654,-59759.11695,-51183.66478,39553.74618,132447.0485,93359.81585,-98263.14966,-255824.9417,-144286.083,212696.4921,445204.5408;
0,0,0,0,192,735.8924507,540.3053094,-2124.343367,-5057.868211,-969.9812354,11405.43117,15940.52452,-5622.481002,-37623.47011,-29881.32575,35855.80045,87872.36995,28692.89466,-112944.5884,-154561.6584,24454.28757,254280.9736,202369.5839,-178659.2813;
0,0,0,0,192,592.6239049,-54.49033533,-2567.50857,-2988.296704,3941.876267,11526.38431,2233.953421,-22136.82236,-24139.16168,20459.96219,60268.45865,14532.0698,-88098.13253,-93277.95948,64910.42847,192595.8473,49803.13745,-244897.5404,-254777.9707;
0,0,0,0,192,458.4013779,-495.3380524,-2477.176593,-714.5839156,6392.684212,6873.986266,-9352.41825,-21134.06087,3844.754647,42020.47316,20505.10932,-58980.138,-70793.54333,51327.19794,141815.5672,6229.868714,-208390.0472,-130756.2533,224730.6889;
0,0,0,0,192,336.5298736,-798.086388,-2066.223865,1211.650343,6609.005295,832.2069828,-14542.95731,-9394.727768,23888.17206,29260.48864,-28610.53289,-63533.47519,18266.13986,110340.0964,20413.44929,-159630.1792,-99650.79156,191428.6434,224840.2275;
0,0,0,0,192,230.010273,-986.6727322,-1517.645248,2561.671612,5397.538212,-4248.815769,-13851.25499,3944.496017,28647.88724,2246.348802,-50351.46819,-19936.09886,77156.40402,55808.14847,-103790.0999,-116346.81,120768.0917,206007.3722,-114289.9057;
0,0,0,0,192,141.4654422,-1089.461027,-968.7547191,3346.025684,3627.717551,-7474.097775,-9965.085871,13492.9252,22500.21035,-20544.91549,-44251.15549,26525.3978,78409.54708,-27774.42335,-127872.1365,18879.12551,194650.9903,7368.517829,-279197.6607;
0,0,0,0,192,73.07564879,-1135.312342,-508.5655728,3706.949976,1946.428626,-9020.436512,-5497.089494,18315.25514,12841.39554,-32819.98602,-26305.72752,53592.96081,48907.83539,-81344.97478,-84369.59991,116249.2281,137093.641,-157744.2065,-212101.4947;
0,0,0,0,192,26.52487642,-1149.801347,-185.5323873,3822.418606,714.8963246,-9523.125445,-2036.182875,19914.17359,4805.440196,-36987.23947,-9962.723932,63039.66842,18780.08487,-100646.6018,-32907.74918,152629.087,54418.22672,-222018.2854,-85848.08848;
0,0,0,0,192,2.959359616,-1151.972632,-20.71532046,3839.781056,79.90093794,-9599.042127,-227.8620289,20156.93483,538.572741,-37623.95402,-1118.549352,64493.6094,2112.7619,-103642.0699,-3710.543262,158327.7533,6151.495274,-232190.8611,-9731.39151;
1,0,-1,0,1,0,-1,0,1,0,-1,0,1,0,-1,0,1,0,-1,0,1,0,-1,0;
0,0,4,0,-16,0,36,0,-64,0,100,0,-144,0,196,0,-256,0,324,0,-400,0,484,0;
0,0,0,24,192,840,2688,7056,16128,33264,63360,113256,192192,312312,489216,742560,1096704,1581408,2232576,3093048,4213440,5653032,7480704,9775920;
0,1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400,441,484,529]

S=[1,0.998458667,0.993839419,0.986156496,0.975433581,0.96170373,0.945009268,0.925401657,0.902941342,0.87769756,0.84974813,0.819179209,0.786085032,0.750567618,0.712736454,0.672708162,0.630606134,0.586560159,0.540706014,0.493185053,0.444143767,0.393733335,0.342109153,0.289430364;
1,0.98618496,0.945121551,0.877944359,0.786509494,0.673343309,0.541572595,0.394838187,0.237194368,0.07299685,-0.093217576,-0.256856394,-0.413398249,-0.558517877,-0.688205612,-0.798878172,-0.887477663,-0.951556077,-0.98934292,-0.999794139,-0.982620967,-0.938297899,-0.868049586,-0.773816993;
1,0.961939766,0.850656228,0.67462034,0.447232036,0.18580022,-0.089774795,-0.35851611,-0.599967012,-0.795748145,-0.930956556,-0.99530012,-0.983880972,-0.897568346,-0.742932397,-0.531744087,-0.280079168,-0.007094493,0.266430219,0.519674138,0.733360219,0.891222577,0.981244656,0.996573933;
1,0.926320082,0.716137789,0.400425549,0.025706666,-0.352800347,-0.679318759,-0.90573287,-0.998678335,-0.944458724,-0.751063831,-0.446992295,-0.077052048,0.304242576,0.640704064,0.882751506,0.994716832,0.960100849,0.784004562,0.492377492,0.128193756,-0.254880591,-0.600395776,-0.857436738;
1,0.880202983,0.549514582,0.087165765,-0.396067449,-0.784405265,-0.984804259,-0.949250027,-0.686261152,-0.258848199,0.230583238,0.664768308,0.939678856,0.989447956,0.802151229,0.422663852,-0.058091262,-0.524928056,-0.86599522,-0.999575095,-0.89366274,-0.573634124,-0.116166194,0.369134463;
1,0.824724024,0.360339432,-0.230362851,-0.740310987,-0.990741662,-0.893865914,-0.483643725,0.096120716,0.642189852,0.963138082,0.946456378,0.597992543,0.039901255,-0.532177495,-0.917700386,-0.981521616,-0.701268527,-0.175184388,0.412310981,0.85526993,0.998412337,0.79155935,0.307223688;
1,0.761249282,0.15900094,-0.51917058,-0.949437402,-0.926346503,-0.460923818,0.224590651,0.802862762,0.997766752,0.716235686,0.092701052,-0.575098468,-0.968287643,-0.899118079,-0.400618342,0.289177228,0.840890257,0.991076981,0.668023024,0.025987114,-0.62845768,-0.98281303,-0.867873748;
1,0.691341716,-0.044093263,-0.75230874,-0.996111568,-0.624998222,0.131936882,0.807425162,0.984476513,0.553794202,-0.218754445,-0.856262349,-0.965185319,-0.4782834,0.303870785,0.8984405,0.93838801,0.399053054,-0.386623964,-0.933631603,-0.904292986,-0.316719326,0.46637042,0.96156198;
1,0.616722682,-0.239306267,-0.911893888,-0.885465021,-0.180278837,0.663100925,0.998177599,0.568096607,-0.297461473,-0.934999082,-0.855808809,-0.120594327,0.707062296,0.992717038,0.517399932,-0.354532491,-0.954696389,-0.823033344,-0.060470273,0.748446566,0.98363822,0.464817436,-0.410311308;
1,0.539229548,-0.418462989,-0.990524765,-0.649777453,0.289766361,0.96227862,0.74801177,-0.155578523,-0.915796843,-0.832070912,0.0184424,0.851960286,0.90036192,0.119043216,-0.771978681,-0.951590646,-0.254272907,0.677367717,0.984786282,0.384684007,-0.569920316,-0.999319756,-0.507805164;
1,0.460770452,-0.575381181,-0.991007746,-0.337872993,0.679643962,0.964192705,0.208899055,-0.771683681,-0.920037132,-0.07616817,0.849845048,0.859335144,-0.057932562,-0.91272237,-0.783178435,0.190991406,0.959184828,0.692936648,-0.320615363,-0.98839682,-0.590232736,0.444473211,0.99983298;
1,0.383277318,-0.706196995,-0.924615899,-0.002571609,0.92264462,0.70982912,-0.378521817,-0.999986774,-0.38802268,0.702546189,0.926562719,0.007714758,-0.920648935,-0.713442468,0.373756304,0.999947095,0.392757778,-0.698876799,-0.928485029,-0.012857704,0.918628896,0.717036943,-0.368980903;
1,0.308658284,-0.809460128,-0.808351431,0.310451397,0.999998222,0.306864073,-0.810565945,-0.807239861,0.312243405,0.999992888,0.305068772,-0.811668881,-0.80612542,0.314034304,0.999983998,0.303272386,-0.81276893,-0.805008112,0.315824085,0.999971552,0.301474921,-0.813866089,-0.803887941;
1,0.238750718,-0.88599619,-0.66181517,0.569978496,0.93398072,-0.124001362,-0.993191548,-0.350249028,0.825947135,0.74463997,-0.47038048,-0.969247324,0.007563492,0.972858903,0.456978031,-0.754651237,-0.81732508,0.364377338,0.991315782,0.10897737,-0.939278931,-0.557484408,0.673079326;
1,0.175275976,-0.938556665,-0.504288846,0.761777225,0.771331339,-0.491385519,-0.943587492,0.160609082,0.999889319,0.18990407,-0.933318077,-0.517080544,0.752054483,0.78071471,-0.478373418,-0.948409446,0.145906635,0.999557301,0.204490127,-0.927872888,-0.529757779,0.742165265,0.789925261;
1,0.119797017,-0.971297349,-0.352514068,0.886837082,0.564994942,-0.751467664,-0.745042111,0.572960019,0.882319914,-0.361561431,-0.968947876,0.1294073,0.999953093,0.110175495,-0.973555702,-0.343433634,0.891271052,0.556976861,-0.757822719,-0.738546663,0.580871344,0.877719972,-0.370574875;
1,0.073679918,-0.989142539,-0.2194398,0.956805927,0.360434564,-0.903692348,-0.49360252,0.830955162,0.616051936,-0.74017385,-0.725123833,0.633319721,0.818449723,-0.512713105,-0.894003042,0.380972963,0.950143155,-0.240960024,-0.985650985,0.095714657,0.999755481,0.051609146,-0.992150366;
1,0.038060234,-0.997102837,-0.113960168,0.988428136,0.18919978,-0.97402616,-0.263343106,0.95398036,0.335960537,-0.928406887,-0.406631304,0.897453922,0.474945916,-0.861300817,-0.540508536,0.820157054,0.602939275,-0.774261035,-0.661876387,0.723878695,0.716978371,-0.669301966,-0.76792595;
1,0.01381504,-0.999618289,-0.041434573,0.998473449,0.069022474,-0.996566352,-0.096557681,0.993898456,0.124019175,-0.990471796,-0.151385989,0.986288989,0.178637233,-0.981353228,-0.2057521,0.975668281,0.232709893,-0.969238489,-0.259490029,0.962068758,0.286072066,-0.954164565,-0.312435708;
1,0.001541333,-0.999995249,-0.004623985,0.999980994,0.007706592,-0.999957238,-0.010789127,0.999923978,0.013871559,-0.999881217,-0.016953859,0.999828954,0.020035998,-0.999767189,-0.023117946,0.999695925,0.026199675,-0.99961516,-0.029281155,0.999524896,0.032362357,-0.999425133,-0.035443251;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
</code>```

Results from eig(H,S):

```<code>
          -10.4013941701241 +                     0i
             495.8709816739 +                     0i
          -2075.89359089271 +                     0i
           3810.96559828315 +                     0i
          -12007.5792049748 -      7859.80373595807i
          -12007.5792049748 +      7859.80373595807i
           13744.8803145627 -      6829.57196964173i
           13744.8803145627 +      6829.57196964173i
           40714.6215488508 +                     0i
           161777.282995074 +                     0i
          -222054.135160397 +                     0i
           456554.085296176 +                     0i
           1046529.30310225 +                     0i
           2209556.44199227 +                     0i
          -3698516.26538781 +                     0i
           5938712.03152599 +                     0i
           29139261.3696759 +                     0i
           35187640.9660181 +                     0i
          -106153683.549455 +                     0i
            625134760.07866 +                     0i
                        Inf +                     0i
                        Inf +                     0i
                        Inf +                     0i
                        Inf +                     0i
</code>```

Results from putting determinant equal to 0:

```<code>6.0880681896251522577558901777469
493.13352335963732050050918987025
3805.0426185157246029406709075097
14617.451723278368450735731815741
39943.815388339366263101056654306
89135.406154470667286843490819588
173881.32991449447500370097212608
308208.37416682144561534408214656
508470.65137071189157395633863007
793718.91013854156465411829146852
1184484.6868564064238468530882583
1708985.9187321214627561917325303
2400030.5803698046017748760892131
3509269.6989698602635588940827738
6117137.3935227069102232614222001
11673551.342361595565898213914818
30121894.392585463384868409705274
104001718.25765847906945083161619
953411363.34772056816450115047618
2832353346.2556813059459515513392
</code>```

The problem is that both of the approaches are giving different eigenvalues. The second method(determinant one) gives correct results as they match with the values present in the research paper that I'm trying to implement. So, I have the following questions:

Why is eig(H,S) not giving correct eigenvalues?

Evaluating the determinant involves symbolic computation and the size of my matrices could be as large as $100\times100$ and this takes very long for MATLAB to solve. Is there any efficient way to evaluate determinants of large(non-sparse) symbolic matrices?


问:
The condtion for Augmented Langrangian Multiplier
I am currently learning the usage of Augmented Lagrangian Multiplier to achieve my equality constraint. I have learnt from the https://en.wikipedia.org/wiki/Augmented_Lagrangian_method that I have two multipliers, $\mu$ and $\lambda$.

According to wikipedia, it shows me how I can update my $\lambda$ with any value of $\mu$. However, this gets my curiosity. Since it did not mention the condition for setting initial $\lambda$ and $\mu$, can both of them be negative value initially, and also during the iterations?


问:
Preserving conservation properties across time-integration schemes
I am interested in deriving an accurate, implicit discretization of a nonlinear advection-diffusion equation
$$
\partial_t f = \frac{1}{v^2} \partial_v J(f,v) \equiv F(f,v) \tag{$\star$}
$$
with flux
$$
J(f,v) = C(f,v) f + D(f,v) \partial_v f
$$
where $f(v,t)$ is the unknown function of a radial coordinate $v\in [0,\infty)$, and $C$ and $D$ are linear integral operators on $f$. In the physics context where this comes up, the operator $F$ satisfies two conservation laws,
$$
\int_0^\infty v^2 F \,dv = 0  \qquad \int_0^\infty v^4 F\,dv=0
$$
and there is also a positivity requirement, $f(v,t) > 0$. Historically, there has been a good deal of effort spent on designing finite-difference schemes such that a first-order update of ($\star$) satisfies these constraints to a high degree of accuracy.

The issue is that in practice, $F$ is not the only operator driving $\partial_t f$, and for some of these other operators - call them $A(f,v,t)$ - I would like to use at least second-order time integration schemes. The naive approach would be to use operator splitting, keeping the conservative first-order discretization for $F$ while using a higher-order scheme for $A$. But the overall accuracy of the scheme can never be better than first order in time, limited by the backward Euler update for $F$. So I am looking for something fancier.

I have heard of "strong stability preserving" methods. My rough understanding is that given a semi-discrete scheme that obeys certain constraints under a forward Euler update, then SSP methods can preserve these constraints for other time-stepping schemes. Are conservation properties like the above included in the class of constraints that an SSP scheme can preserve? My thought is to find a second-order (or higher) SSP scheme that is suitable for $A$ without ruining the conservation properties of $F$. Does this make sense, or is my understanding wrong?

Edit:

In response to comments, I'll write out the discretization of the operator $F$ (note I have made a minor change to the definition of the flux above). The nonlinearity is handled by simple fixed-point iteration, and I use a superscript "$*$" to mark quantities that are treated as "constant" and then re-evaluated at each stage of the iteration.

The equation is implicitly differenced in time and the divergence of the flux $J$ is computed from a central difference of values at cell boundaries (fractional indices)
$$
\frac{f^{n+1}_k - f^n_k}{\Delta t} = \frac{J^{n+1}_{k+\frac12} - J^{n+1}_{k-\frac12}}{v_k^2 \Delta v}
$$
The flux at cell boundaries is constructed as
$$
J^{n+1}_{k+\frac12} = C^*_{k+\frac12} f^{n+1}_{k+\frac12} + D^*_{k+\frac12} \frac{f^{n+1}_{k+1} - f^{n+1}_k}{\Delta v}
$$
The value of $f$ at the cell boundary is interpolated using the values from adjacent cell centers:
$$
f^{n+1}_{k+\frac12} = (1 - \delta^*_{k+\frac12}) f^{n+1}_{k+1} + \delta^*_{k+\frac12} f^{n+1}_k
$$
where the interpolation weight is
$$
\delta^*_{k+\frac12} = \frac{1}{w^*_{k+\frac12}} - \frac{1}{\exp(w^*_{k+\frac12}) - 1} ,\qquad w^*_{k+\frac12} = \Delta v \frac{C^*_{k+\frac12}}{D^*_{k+\frac12}}
$$
The "coefficients" $C^*_{k+\frac12}$ and $D^*_{k+\frac12}$ are evaluated at each iteration using midpoint quadrature. The $C$ quadrature needs only $f^*_k$ (cell centers), but the $D$ quadrature uses $f^*_{k+\frac12}$ (cell edges). Since the interpolation weights $\delta^*_{k+\frac12}$ are themselves functions of $D^*_{k+\frac12}$, an inner iteration is performed to ensure the $\delta$'s and $D$'s are consistent.

The comment by @DavidKetcheson mentions that if the above scheme were linear, then it could in principle be used in an RK scheme and still have the desired conservation properties. However, after writing things out, I think the dependence of the interpolation weights $\delta_{k+\frac12}$ on the unknowns $f_k$ is nonlinear, so there is no guarantee that this scheme can be trivially made to work with a higher-order time-integrator without sacrificing conservation properties. Is that the end of the story?


问:
How can I couple MCNP with FLUKA? To attain reactor dynamics and simulate an ADS experiment? (Under Subcritical fast system)
I'm new in Particle Physics based Computing. I want to simulate an ADS experiment (Accelerator Driven System) in subcritical fast system. I was advised to simulate the experiment with FLUKA with MCNP being coupled with it.

On the otherhand, I saw CERN now uses Flair more than FLUKA. is FLUKA and Flair same and if someone can mention how I can couple Flair with MCNP that will be helpful!

I don't know how to actually couple both of them with each other. If anyone can help me out or suggest something that would be wonderful?

(Edited) Note: I'm trying to analyze Nuclear waste Neutralization using proton beam through ADS concept. And I don't know to couple MCNP with FLUKA!

FLUKA: http://www.fluka.org/


问:
Stiff ODE solver in the web browser
I'd like to make a web application that lets people play with solving ODE systems, changing parameters with sliders etc. but instead of doing the computations on the server side, solving the equations on their own computer via the web-browser.

There's some really nice examples here: http://blog.littleredcomputer.net/jekyll/update/2016/04/03/diffeq-javascript.html using this javascript library
https://www.npmjs.com/package/odex

But the ODE system I want to play with is pretty stiff, so would almost certainly need an implicit, and probably adaptive timestep, solver like CVODE. Does anyone know whether any libraries exist for that?

Title: Stiff ODE solver in the web browser