问:
Jacobian Matrix of 2D element mapped to 3D
Note: I previously posted this question to MathStackExchange, but got no attention there. So I'm rewritting and trying over here.

Problem summary

Given a common¹ set of shape functions defined at master element space $\{\psi_j(\xi,\eta)\}$ and associated to each corner $j$ of that element;

We can compute the affine transformation that maps parametric space coordinates $(\xi,\eta)$ into 'real' space coordinates $(x,y,z)$ as an interpolation of corner coordinates $\mathbf{p_j}$:

$$\bf x(\xi,\eta) = \sum_j\psi_j p_j$$

In similar fashion, we can compute the gradient² of that map as a rectangular matrix and using the partial derivatives of shape functions:

$${\rm grad}_\xi \mathbf{x}=
\sum_j
\left[
\begin{aligned}
\frac{\partial \psi_j}{\partial \xi}p_{jx}
&&
\frac{\partial \psi_j}{\partial \eta}p_{jx}
\\
\frac{\partial \psi_j}{\partial \xi}p_{jy}
&&
\frac{\partial \psi_j}{\partial \eta}p_{jy}
\\
\frac{\partial \psi_j}{\partial \xi}p_{jz}
&&
\frac{\partial \psi_j}{\partial \eta}p_{jz}
\end{aligned}\right]$$

(I define the Jacobian Matrix down below)

Prove the Jacobian matrix of that transformation can be computed through the $QR$ decomposition of ${\rm grad}_\xi\bf x$, where $[R]_{2\times 2}$ is the Jacobian matrix, and $[Q]_{3\times 2}$ is a unitary matrix whose columns are orthonormal vectors that exist in the plane of the mapped element.

Context:

So last semester I took a course in intro to Finite Element approximations. In the course, the professor instructed us to compute the Jacobian Matrix of the geometrical map from the parametric space (of the master element) to 'real' space through the $QR$ decomposition of the 'gradient' of that map, in which

$R$ is the Jacobian matrix we seek and

$Q$ is a unitary matrix whose columns are orthonormal vectors in the plane of the mapped element.

He justified it by explaining that, in mapping 2D-to-2D, or 3D-to-3D elements, the ${\rm grad}_\xi\bf x$ exactly matches the Jacobian matrix, but that is not the case when mapping 2D master elements to 3D. For those, the linear geometrical map is an affine transformation from their 2-dimensional parametric space (of a triangle or a quadrilateral) into the actual coordinates of the element, which exist in $\mathbb{R}^3$.

Although I trust that he is right and I kinda have an intuition for it, he didn't prove that $R$ is indeed the so-called Jacobian matrix.

When I asked him for a proof, he seemed annoyed and brushed me off, which gave me the impression that this is supposed to be trivial and I was asking a dumb question. Yet, none of my classmates managed to prove it either.
Now the course is over I decided to go back to the problem that's been bugging me for a few months.

Traditional Jacobian Matrix definition

(for mapping from 2D isoparametric space to 2D 'real' space as an example)

The Jacobian Matrix $\mathbf{J}$ is the linear transformation that maps line segments $d\xi$ and $d\eta$ in the parametric space $(\hat\Omega)$ to the real space $\mathbb R^2$

$$
\begin{bmatrix}
dx \\ dy
\end{bmatrix}
=
\begin{bmatrix}
\dfrac{\partial x}{\partial \xi}
&&
\dfrac{\partial x}{\partial \eta}
\\
\dfrac{\partial y}{\partial \xi}
&&
\dfrac{\partial y}{\partial \eta}
\end{bmatrix}
\begin{bmatrix}
d\xi \\ d\eta
\end{bmatrix}
=
\mathbf{J}
\begin{bmatrix}
d\xi \\ d\eta
\end{bmatrix}
$$
$$
\mathbf{J} =
\begin{bmatrix}
\dfrac{\partial x}{\partial \xi}
&&
\dfrac{\partial x}{\partial \eta}
\\
\dfrac{\partial y}{\partial \xi}
&&
\dfrac{\partial y}{\partial \eta}
\end{bmatrix}
$$

Notes:

By 'common' I mean the classic FEM linear shape functions. Valued 1 at their corresponding node and 0 on all other nodes.

Pardon me if terms 'gradient' and 'Jacobian Matrix' I used aren't exactly their usual formal definitions, but I'm matching the terms used by the professor. I've put some effort into trying to keep it unambiguous;

I'm a last year engineering student so, if you could keep your proof within the realms of my elementary maths, I'd appreciate it. But I'll take any help I can get;

This is not homework, so feel free to answer a complete proof. I'm just honestly trying to learn.

There's a question in math.stackexchange that seems to hint me in the right direction, but I haven't been able to complete this proof using it either.

There's a closely related question here in which the poster is asking for the solution my professor gave us. They discuss it a bit there (applied to triangles), but it doesn't solve my problem. I'll gladly reply to them once I can actually prove it.

I'll take any help anyone can give me. If you can't prove it but help me prove it, by sharing some knowledge or pointing to a reference, I'll be very much thankful.


问:
Path constraints for state variables - fmincon, ODE45

My problem lies in constraining state variables (look for 23,24,25). I am currently using ODE45 to solve equations and fmincon to find best control variable.How would you go about solving this?  Here is my code so far:

```<code>function ster_6

optionsopt = optimset('Display','iter-detailed','Algorithm','sqp', 'TolFun', 10^(-8), 'TolX', 10^(-8), 'MaxFunEvals', 10^6, 'MaxIter', 50);

u_lb = [ 0.0];

u_ub = [ 50.0]

u0 = (u_lb+u_ub)/2;

[rozw, dokladnosc] = fmincon(@model_procesu,u0,[],[],[],[],u_lb,u_ub,[],optionsopt)

end

function J = model_procesu(u)

tspan_1 = [0.0  140.0];
x0_1 = [ 1.5    0.0    0.0    7.0 ];

options_ode = odeset('RelTol',1e-6,'AbsTol',1e-6);

[tsol_1,xsol_1] = ode45(@(t,x) penicillin_production_problem( t,x,u(1) ), tspan_1, x0_1, options_ode);

J = -(xsol_1(end,2) * xsol_1(end,4));

plot(tsol_1,xsol_1)
title('Trajektorie zmiennych stanu')

ylabel('Stężenie')

xlabel('Czas')
legend('x1_1(t)','x2_1(t)','x3_1(t)','x4_1(t)')

end

function dx = penicillin_production_problem(t,x,u)

dx = zeros(4,1);

h1 = 0.11 * (x(3)/(0.006*x(1)+x(3)));

h2 = 0.0055 * (x(3)/(0.0001 + x(3) * (1 + (10 * x(3)))));

dx(1) = h1 * x(1) - u * ( x(1) / (500*x(4)) );

dx(2) = h2 * x(1) - 0.01*x(2) - u * ( x(2) / (500*x(4)));

dx(3) = (-((h1 * x(1)) / 0.47)) - (h2 * (x(1) / 1.2)) - (x(1) * ((0.029*x(3)) / (0.0001 + x(3)))) + ((u/x(4)) * (1-(x(3)/500)));

dx(4) = u/500;

end

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

ContentLicense: CC BY-SA 4.0


问:
Using the lower triangular portion of a matrix to return a Symmetric Positive Definite Matrix?
I've recently posted this question.
To summarise, I'm dealing with supposedly Symmetric Positive Definite(SPD) matrices, but due to machine-precision they end up not being SPD.
In a comment, a user suggested to use the lower triangular portion of the matrix.

I would like to be sure in what way we can use that? What algorithm is he possibly alluding to? I thought of the Cholesky Decomposition...

With that method, or another which uses the lower triangular part of the matrix, will I be able to take any 'almost' SPD matrix, and return an SPD matrix?


问:
Number of words that a processor can handle for PMMM
The PMMM which stands for parallel matrix-matrix multiplication essentially accelerates the algorithm of the matrix-matrix multiplication of two matrices $A$ and $B$ both of size $n$ so that $C:= AB$. Let us say we have $4$ processors and I wish to know how many words are stored per processor.

To my knowledge, since we have $4$ processors then we can convert $A, B$, and $C$ to block matrices. Each block matrix has $4$ sub-blocks since we have $4$ processors so processor $i$ should store the matrix multiplications that will occur between sub-blocks from matrix $A$ and $B$ where ($i=1,2,3,4$). Now each sub-block must have $n^{2}/4$ elements. Taking processor $1$ for example, we have :
$$
C_{1,1}=A_{1,1}B_{1,1}+A_{1,2}B_{2,1}
$$
We may see that we should have a total of $n^{2}/4+4n^{2}/4=5n^{2}/4$ words per processor. Is my assertion correct? It might appear that the value obtained is illogical but let's say $n$ is divisible by $4$.


问:
Why does including the pressure in this FVM for Stokes 2nd Problem lead to wrong solutions?
I'm trying to learn how to use finite volume methods and I want to solve a more general case of Stokes' second problem i.e. an infinite half-plane oscillating harmonically with no-slip boundary conditions.

I've been following the code from this lecture and changed the boundary conditions to:

```<code>u[1:-1, 1]  = u[1:-1, 2]                #left
u[1:-1, -1] = u[1:-1, -2]               #right
u[-1,1:]    = 2*np.cos(n*dt) - u[-1,1:] #top
u[0, 1:]    = u[1,1:]                   #bottom

v[1:,0] = v[1:,1]      #left
v[1:,-1] = v[1:,-2]    #right
v[-1,1:-1] = 0.0       #top
v[1,1:-1] = 0.0        #bottom
</code>```

(My thinking was that this would imply no gradients at the left, right, and bottom boundary so the flow would just continue out of/into the volume and to oscillate the top by setting the value through a ghost cell.)

This leads to vortex formations however

But we know by the analytical solution
$$ u(y,t) = e^{-y\sqrt{\frac{\omega}{2\nu}} } \cos\left(\omega t - y\sqrt{\frac{\omega}{2\nu} } \right)$$
that there will be no velocity component in the y-direction. And it's not like this vortex is just a small deviation.

However if I simply exclude the projection step/set the pressure to 0 everywhere I recover the analytical solution. Why does this happen? From what I understand the pressure is essentially used to enforce that the flow field is divergenceless, however $u(y)$ obviously also has 0 divergence so why does the projection method fail here?


问:
How to distinguish primary hosts (stars) and orbiting satellites (planets) and tertiary bodies (moons) by their mass and trajectory?
I posted this question in the astronomy stackexchange. There are no responses, and it was suggested that I pose the question here. The "too long, didn't read" was taken from a comment, and my question is copy-pasted below.

TLDR:

You're asking about how some kind of classifier might work that
accepts a group of trajectories or state vectors for  bodies over an
extended period of time, and identifies orbital hierarchies; i.e. "Who
orbits who orbits who?"

Suppose one has run a gravitational simulation of N bodies (has the mass, vector positions, vector velocities, etc for each body), but knows nothing a priori about the system (ie, one does not know anything regarding composition, or which bodies are stars or planets). I would like to use the simulation results to categorize which bodies are primary hosts and which ones are satellites (and if possible, which "tertiary" bodies are orbiting the satellites). Is there a way to create a hierarchy that identifies whether a body is a primary host or satellite or tertiary? I would like to find a method of doing this such that the method can be applied to cases for N > 10^3 bodies.

For context, consider a 9-body system like a simplified example of our solar system (8 planets + 1 star). One can identify the Sun as the only primary host and the planets as the satellites because the ratio of the mass of the Sun to the mass of Jupiter (most massive planet of our solar system) is 1048. One could also use the approach of noticing that the center of mass would be inside of the Sun (or very close to in a point-mass simulation). But, the center of mass approach only works for nearly circular stable orbits around a single host; for unstable orbits, a body can deviate from its tangential path about the center of mass.

If we consider a binary star system for which both stars are of comparable mass, then I think that the ratio of the mass of any single host star to the mass of the most massive planet orbiting this binary star system can decrease; in other words, for each primary host body added to a system, the combined gravitational force at certain points in the orbit required to keep a satellite body in orbit should decrease if the system is stable, thereby decreasing the ratio of the mass of any host body to the mass of the most massive satellite body. (I suppose the same logic could be extended to host planets and satellite moons, though I suspect the ratio of masses would decrease; would this logic also apply to unstable orbits?)

Is this line of thinking correct? More importantly, what is the right way to think about this?


问:
Pade-like approximation, but force poles to be negative
Are there techniques to form a Pade approximation (or Pade-like approximation), except force the poles of the rational function to be negative?

I am trying to use Pade approximations to extrapolate a curve based on known derivatives of the curve at a point. The Pade approximation usually does much better than a Taylor series. However, occasionally the Pade approximation creates a spurious pole in the direction I'm extrapolating in, and in this case the approximation is much worse because of the spurious pole.

Has this been studied? I looked around and couldn't find much practical information on this. What do people normally do to avoid spurious poles in Pade approximations?

I would be happy with losing some accuracy in the approximation, if this ensures that the poles are negative (i.e., behind the "starting point" for the curve).


问:
Equivalence between zero sum games and linear program
It is well known that you can use the algorithm for finding the equilibrium of a Zero-sum game to solve a linear program. In particular, you can take a LP and reduce it to a zero-sum game, and use the famous algorithm for zero sum games proposed by Grigoriadis  and Khachiyan.

Is this equivalence used in practice, or the simplex algorithm (or any interior point methods) is always preferable? Why?

Notably, the approach of zero sum games offers an approximate solution in polynomial time, and converges in $O(1/\epsilon^2)$ iterations, contrary to the simplex algorithm which find the optimal solution in (worst-case) exponential time.


问:
Interpolation is correct for grid, but wrong for random points
I am using scipy interpolation to construct a map of two coordinates to two other coordinates --- effectively one set of coordinates are 2D cartesian coordinates (called X and Y), and the other are similar to polar coordinates (radial variable r and angle coordinate t) but with no analytical expressions in term of the cartesian coordinates. The variable r labels different contours, and t is the angle around the contour.

First:

```<code>import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d, RectBivariateSpline
from skimage import measure
</code>```

I start with a data set containing r as a function of X, Y at specific points. X_data has nX points, Y_data has nY points and r_data is a nX by nY array.

```<code>X_data = np.linspace(1.0, 3.0, 1000)
Y_data = np.linspace(-3.0, 3.0, 1000)
XX, YY = np.meshgrid(X_data, Y_data)
X_centre = 2.0
Y_centre = 0.0
r_data = ((XX-X_centre)**2 + (YY-Y_centre)**2)**0.5 + 0.25*np.sin(YY)
</code>```

Note that this is just a mockup of the original data, which has no analytical form. I build a 2D array of X and Y on a regular grid of (r, t):

```<code>nR, nT = 250, 250
t_array = np.linspace(0.0, 2.0*np.pi, nT)
r_array = np.linspace(0.01, 1.0, nR)
X_grid = np.zeros([ nR, nT])
Y_grid = np.zeros([ nR, nT])
</code>```

Now I want to find contours of r in r_data, using the following function, which uses measure.find_contours from skimage to find contours, then processes that output to get the X, Y coordinates along a complete contour of r=constant. Then I calculate the angle along the contour (Contour_t).

```<code>def GetContour(r, r_data, X_data, Y_data, X_centre, Y_centre):
    Contours = measure.find_contours(r_data, r) # gives X,Y indices that form contour

    Contour_X, Contour_Y = np.array([]), np.array([])
    for contour in Contours: #combine lists of partial contours for full contour
        Contour_X = np.concatenate((Contour_X, X_data[0] + contour[:,0]/len(X_data)*(X_data[-1] - X_data[0]) ))
        Contour_Y = np.concatenate((Contour_Y, Y_data[0] + contour[:,1]/len(Y_data)*(Y_data[-1] - Y_data[0]) ))

    # re-arrange contour so angle is taken from Y_data = 0 on right hand side of X_centre (r=0)
    PossibleInds = np.argsort(np.abs(Contour_Y)) #indices for the best candidates for the start
    j = 0
    while True:
        if Contour_X[PossibleInds[j]] >= X_centre:
            if Contour_Y[PossibleInds[j]] >= Y_centre:
                break
        j += 1
    Ind = PossibleInds[j]

    if Ind!=0: # if it's not already in the right order (by chance)
        Contour_X = np.concatenate(( Contour_X[Ind-1:], Contour_X[:Ind]))
        Contour_Y = np.concatenate(( Contour_Y[Ind-1:], Contour_Y[:Ind]))

    Contour_l = np.zeros(len(Contour_X)) # length along the contour
    Contour_l[1:] = np.sqrt(np.diff(Contour_X)**2 + np.diff(Contour_Y)**2 )
    Contour_l = np.cumsum(Contour_l)
    Contour_t = Contour_l/Contour_l[-1] * 2.0*np.pi  # angle along contour

    return Contour_X, Contour_Y, Contour_t
</code>```

Now find each contour for every element in r_array, then plot them over the data for a sanity check:

```<code>for i, r in enumerate(r_array):
   X, Y, Contour_t = GetContour(r, r_data, X_data, Y_data, X_centre, Y_centre)
   fX = interp1d( Contour_t, X) # interpolate X as a function of angle along contour
   fY = interp1d( Contour_t, Y)
   X_grid[i,:] = fX(t_array)    # now save X at fixed values of t given by t_array
   Y_grid[i,:] = fY(t_array)
   plt.plot(X_grid[i,:], Y_grid[i,:])

plt.contourf(X_data, Y_data, r_data.T)
plt.show()
</code>```

Great. I have found the contours successfully. Now I use RectBivariateSpline to construct a map of X(r,t) and Y(r,t):

```<code>X_map = RectBivariateSpline( r_array, t_array, X_grid )
Y_map = RectBivariateSpline( r_array, t_array, Y_grid )
</code>```

I can check this works by plotting the maps using something like:

```<code>r_test = np.linspace(0.0, 1.0, 50)
t_test = np.linspace(0.0, 2.0*np.pi, 50)
rr, tt = np.meshgrid(r_test, t_test, indexing='ij')
plt.contourf(rr, tt, X_map(r_test, t_test))
plt.show()
plt.contourf(rr, tt, Y_map(r_test, t_test))
plt.show()
</code>```

and it gives exactly what I would expect. However, I actually want to randomly generate points in r, t, then use the maps to get the X, Y coordinates for the randomly generated points.

```<code>N = 100000
r_pts = np.random.uniform(0.01, 1.0, size = N)
t_pts = np.random.uniform(0.0, 2.0*np.pi, size = N)
X_pts = X_map(r_pts, t_pts, grid=False)
Y_pts = Y_map(r_pts, t_pts, grid=False)

plt.hist2d(X_pts, Y_pts, 100)
plt.show()
</code>```

Now all the points are clustered at the centre, not up to r=1.0, and they definitely are not uniform in angle!

QUESTION(S): I suspect I am messing up the indexing, but I have played around and can't see an obvious mistake. Can you? Why does it look fine when I plot the maps, but using the maps on randomly generated points gives the wrong result?

I have tried various interpolation methods (interp2d, RegularGridInterpolator, map_coordinates) and these have the same issue...

Thanks!


问:
Method of Lines for Wave Equation with Periodic Bounday Condition
I am trying to solve, using the method of lines in Matlab, the wave equation IVP:

\begin{cases}
      \frac{\partial^2 v}{\partial t^2} = \frac{\partial^2 v}{\partial x^2} \\
      v(x, 0) = e^{-x^2} & -5\lt x\lt 5 \\
      \frac{\partial v}{\partial t}(x, 0) = 0 \\
      v(-5, t) = v(5, t)
\end{cases}

I rewrote the wave equation as two coupled equations:

\begin{equation}
\frac{\partial u}{\partial t} = w\\
\frac{\partial w}{\partial t} = \frac{\partial^2u}{\partial x^2}
\end{equation}

The second equation is just the heat equation, which I think can be represented by

\begin{equation}
\frac{d}{dt}\begin{pmatrix}
w_1\\w_2\\\dots\\w_n
\end{pmatrix} = \frac{1}{\Delta x^2}\begin{pmatrix} -2 & 1 & 0 & \dots & 1\\ 1 & -2 & 1 & \dots & 0 \\ 0 & 1 & -2 & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots \\1 & 0 & 0 & \dots & -2\end{pmatrix}\begin{pmatrix}
w_1\\w_2\\\dots\\w_n
\end{pmatrix}
\end{equation}
where the matrix A is n x n and tridiagonal and the 1's in the top right and bottom left corners provide the periodic boundary condition; or
\begin{equation}
\frac{d}{dt}\vec{w} = A\vec{w}
\end{equation}
and
\begin{equation}
\frac{d}{dt}\vec{u} = \vec{w}
\end{equation}

I am a little unsure how to plug this into, for example, ode45 in Matlab.

I have done the following:

```<code>function solns = coupledODEs(t,U)
global A n

len = length(U);
u = U(1:len/2);
w = U(len/2+1:len);

solns(1:len/2) = A*w; %W'
solns(len/2+1:len) = u;  %U'
solns = transpose(solns);
</code>```
```<code>global A n

xi = -5;
xf = 5;
n = 1000;   % no of grid points Xn
h = (xf-xi)/(n+1.); % grid spacing in x direction

T = 50.; % length of time for solution

r = 1/(h^2);

A = 2*-r*diag(ones(1, n)) + r*(diag(ones(1, n-1), 1) + diag(ones(1, n-1), -1));
A(n, 1) = r;
A(1, n) = r;

% set up mesh in x-direction:
x = xi:h:xf;
x = x';
x1 = x(1:n);
U0 = exp(-x1.*x1);
W0 = 0*x1;

ICs = [U0; W0];
[t, results] = ode45('coupledODEs',[0,T], ICs);

resultsize = size(results);

U = results(:, 1:resultsize(2)/2); %capture the first half of the results, which returns U?

</code>```

Visually the wave appears to be passing through the sides like a periodic boundary allows, but I'm not sure how to confirm that it is actually correct. I am not terribly familiar with Matlabs ODE functions. Am I doing this correctly? Thank you in advance for any assistance.


问:
Uses for Automatic Differentiation in Investment Banking
I've finally wrapped my head around the advantages and disadvantages of AD/AAD compared to FD and SD, but could someone please explain to me where you'd use AD/AAD in an investment bank/banking generally?

If this question has already been answered, feel free to point me in the right direction.


问:
A JAVA solver for ODEs with boundary conditions (BVP)
I need to solve a system of linear first order ODEs with boundary conditions in JAVA. I was wondering if any of you know of a JAVA package with the capability of solving a boundary value problems? (i.e. a solver using finite differences or finite elements)

The type of equation I want to solve is of the following type:
$$
n''(x) = -D(x)n(x) + n'(x)/x \\
n'(0) = 0 \\
n(1) = 0
$$
Through substitution this equation can be transformed into two first order differential equations.

I know there is the ode solver of Apache commons but it can only solve initial value problems. I have also read about the numerical methods but would preferably use existing code.


问:
Unstable Algorithms which become stable when hardware provides Kulisch exact dot product instruction
In John Gustaffson's book The End of Error, he discusses Ulrich Kulisch's exact dot product, which (in double precision) requires a 2100 bit fixed point register which rounds only once after the computation of the dot product is complete.

I have found only one statement about how this instruction can be used in numerical analysis, namely Kulish claims "it is the EDP which makes residual correction effective . . " However, I could not find a reference demonstrating this claim.

What numerical methods would benefit from an exact dot product? And more interestingly, are there algorithms which become stable in the presence of a hardware dot product which are unstable otherwise?


问:
Expression for numerical amplification factor for Euler time integration and CD2 scheme for one degree hyperbolic function

This is the expression to be derived but I am not getting the exact expression as given in the image. Is the given expression given for the implicit Euler method?

Here $N_c = c (\delta t)/h$ and $\tan(\beta_2) = N_c \sin (kh)$.


问:
For a hyperbolic PDE, is there any proof that the BDF2 method is stable for integrating them?
I would like to ask a question on the stability of BDF2 applied to hyperbolic PDEs.

Say I have a hyperbolic equation as $\frac{\partial c}{\partial t} + {\bf U} \cdot {\bf \nabla}c=0$. This system is coupled to the Navier-Stokes equation (the governing equations for $\bf U$) : $\frac{\partial {\bf U}}{\partial t} + {\bf U} \cdot {\bf \nabla}{\bf U}= -\nabla p + \frac{1}{Re}\nabla^2{\bf U} + c^2 {\bf j}$, where $\bf j$ is in the vertical direction. I have applied BDF2 to these equations and the solutions seem stable.

The thing is that some people add an artificial diffusion term in the $c$ equation in order to secure the numerical stability (with some explicit time-integration method). My experience is that the BDF2 method is stable. Is there any proof that BDF2 is always stable in this case?

Thank you very much. Any help is appreciated.


问:
Constrained optimization for non-linear equations in octaveGNU
I have installed Optim1.6.1 package. I would like to solve a system of equations in non linear finite element analysis using constraints as u=1 at certain nodes. u=0 at certain nodes. Typically I find kmat and res as below

[kmat,res]=systemmatrices(u,load,boundarycondition,material,constraints)
The optimisation syntax in octave is as follows
[x, obj, info, iter, nf, lambda]=sqp(x0, @g, @h, upperbound, lowerbound, iter, tol)

Here, 'g' seems to be the objective function which is given as a function of variables and equations. Similarly 'h' is considered the system of equations in symbolised form such as x1, x2 etc. I see several examples in this template. I'm not clear on how to input the stiffness matrix, residual matrix, since the matrices are computed in a separate script as a function of several variables.

In general what is the objective function in a non linear FEA problem with mechanical loading? It the energy residual? Can I find a simple example of a simple mesh?

I feel quite stuck with this problem for quite a long time. If these are too basic please suggest me some reading material.


问:
Transient advection equation with stabilized FEM
I am interested in solving the transient advection equation

$\left\{\begin{array}{ll}\partial_{t} u+\beta \cdot \nabla u=f & \text { in } \Omega, t>0 \\ u=0 & \text { on } \partial \Omega^{-}, t>0 \\ u(\cdot, 0)=u_{0} & \text { in } \Omega,\end{array}\right.$

with a stabilized FEM formulation (SUPG, GLS or similar), but I am having issues with the preconditioner using a backwards Euler method. I am wondering if one should actually use an explicit method given that it is an advection problem, if there are special time stepping algorithms for stabilized FEM, or if one should be careful when discretizing the stabilized formulation in time (references will be appreciated)


问:
How to interpolate stress at unknown points from the stress values available based on geometrical position for constant load?
I am working on a combined contact, bending, and torsion problem. I have data on geometrical points and their instantaneous stress components. However, based on the available data, I have to interpolate the stress values at points where the stress values are not available, i.e., I would like to interpolate the stress values at points where the stresses are not calculated by the computational model. Any inputs on the interpolation functions are highly appreciated.

I did try interpolation function and scattered data function available in MATLAB. I am getting an error as the points are not monotonically increasing and the points are not under a convex hull. I did use mldivide in MATLAB to obtain the fit based on the available data. However the accuracy of the interpolated stress is a question. I am still contemplating whether the use of Airy's stress equation and Prandtl stress equation will be appropriate to this condition. Any inputs on interpolation functions or directions on how to arrive at interpolation function will be highly appreciated.

The mesh is made up of mid noded quadrilateral element. It is 8 noded higher order element. As I mentioned earlier the computational model calculates only instantaneous stress values based on which point is in contact. For other nodes which are not in contact/or will be in contact in future there is no stress value. The idea is to obtain realistic stress histories of all the contact points.


问:
Solving PDEs: What is the best way to deal with non-banded/dense jacobians?
I have a system of PDEs describing atmospheric chemistry and transport. I use finite-differences to make my system of PDEs into a system of ~10,000 ODEs. I then integrate the ODEs forward in time with the CVODE BDF method (from Sundials). For my current discretization of my PDEs, the jacobian is banded (band width = ~200), so I have been using CVODE's direct banded solver.

However, I'd like to couple additional PDEs to my system, which will cause my jacobian to no longer be banded. Of the options below, what is likely the best way to dealing with this new system of equations, which no longer has a banded jacobian.

Use GMRES to do the linear solves, with the CVODE BDF method

Try to work out the sparsity pattern of my new jacobian, and use a direct linear solver for that sparsity pattern (again sticking with the CVODE BDF method)

Use a split implicit/explicit ODE solver, like ARKODE, instead of CVODE BDF. The new equations I'm adding are not stiff ODEs. It seems like I should be able to solve my stiff ODEs with the implicit method and a direct banded matrix solver, and my new non-stiff ODEs with an explicit method (no jacobian required).

Other?

Thanks for any feedback.


问:
Hanging nodes in deal.ii tutorials: how is the continuity constraint imposed?
While looking at step6 of deal.ii tutorials, I decided to try to understand how the constraints coming from hanging nodes are imposed. So I started by watching video lecture 16 by prof. Bangerth

As far as I understood: some of the basis functions $\varphi_i$ are now discontinuous, and that's of course a problem as we want to compute $L^2$ products of gradients on local cells. Also, our space is no more a subspace of $H^1$. In practice, what they do is to ensure that every linear combination $\sum_i V_i \varphi(x)$ of such basis functions is indeed continuous, by playing with the cofficients and imposing continuity at the hanging nodes.

What I can't understand is the sentence at 10:10, when basically he says:

I need that the value of the function $v_h$ at vertex 1 is one half the value at vertex 0 plus one half the value at vertex 2

I know this follows by imposing continuity, but I can't see the computation behind this. I can't see how the continuity is imposed, mathematically.


问:
Do the class of Construction Heuristic that does or does not change previous decisions have a name?
I'm writing a paper where I am discussing different types of construction heuristics.

One type does not change previous decisions when adding new elements to the solution. I'd call them "immutable", as the (partial) solution is immutable during construction. This behavior is known in many (possibly most) construction heuristics.

The other type can change/adapt previous decisions (at least to some extend) when adding new elements. I'd call them "mutable", as the partial solution can be changed during construction.

My question is: Does anybody know a name for those types or for this characteristics of construction heuristics (including a citation) that is already "known" in the scientific community? Or am I the first to discuss this feature (probably not)?


问:
Strange Picard iteration
I am interested in solving the equation
$$
\begin{aligned} \nabla \cdot\left(\nabla \phi-\frac{\nabla \phi}{|\nabla \phi|}\right) &=0 & & \text { in } \Omega \\ \phi &=0 & & \text { on } \Gamma\left(\phi^{0}\right) \text { and } \\ \nabla \phi \cdot \hat{n} &=\operatorname{sign}\left(\phi^{0}\right) & & \text { on } \partial \Omega
\end{aligned}
$$
from the paper A high-order elliptic PDE based level set reinitialisation method using a discontinuous Galerkin discretisation. They apply a Picard iteration method (Section 2.2) like
$$
\nabla \cdot \nabla \phi^{m}=\nabla \cdot \frac{\nabla \phi^{m-1}}{\left|\nabla \phi^{m-1}\right|} \quad in~ \Omega
$$
where they solve for $\phi^m$ with $\phi^{m-1}$ from the previous iteration. However, the Picard iteration that I am used to seeing would have been something like
$$
\nabla \cdot\left( \nabla \phi^{m} - \frac{\nabla \phi^{m}}{\left|\nabla \phi^{m-1}\right|}\right)=0\quad in~ \Omega
$$
where one simply assigns the nonlinear terms with values from the previous iterations. This approach does not work however, whereas the one they describe does. I am wondering why. It might have to with the fact that the diffusion term at the current iteration
$$
\left(1 - \frac{1}{\left|\nabla \phi^{m-1}\right|}\right)
$$
is zero for $\left|\nabla \phi^{m-1}\right|=1$? And that might require the different strategy they employ? If so, does anybody know of a theory that back that approach?

Title: Strange Picard iteration


问:
Reproducing a paper's result for Topological Insulators
For the past weeks I have been trying to reproduce Agarwala's results but I've been unsuccessful. From this paper I am trying to reproduce the first and last columns of Fig.2, by implementing eq.2; The image below is the crystalline case from the paper:

The code is made of two parts mostly: the building of the system and the method that checks which sites are connected to which. The latter method I have checked, tested and improved it and it appears to be 100% functional. The building of the system, although appears to be correct, is likely the part where the issue arises.

I have tried several things and changed the code a whole lot but the results are always incompatible with the reference. My output is:

The results are clearly different but I cannot find any issues in the code and bear in mind that since it is a crystalline structure the result are supposed to be, at least, very similar to one another.
Here's my github with the code I implemented and two images: one for the crystal case and another for the random case.

Thanks in advance

PS: I was going to post this on physics exchange but apparently since it is about writing a code I am supposed to ask here.


问:
Viable algorithms for efficiently solving a block matrix system with non-uniform sparsity structure
I am trying to use Newton's method to get a stationary solution for a system of equations of the following form:

$$
\begin{Bmatrix}
\frac{\partial x}{\partial t} \\
0
\end{Bmatrix} = \begin{Bmatrix}
f(x, y) \\
g(x, y)
\end{Bmatrix}
$$

For vectors $x$ and $y$.

Due to the nature of the equations I'm trying to solve, the Jacobian can be expressed in the form:

$$
J = \begin{bmatrix}
\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\
\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}
\end{bmatrix} = \begin{bmatrix}
A & B \\
C & D
\end{bmatrix}
$$

Such that $A$ and $C$ are dense matrices, while $B$ and $D$ are sparse. It should also be noted that the elements in the diagonals of $A$ and $C$ are much larger than those in $B$ and $D$ (by ~100 times).

Newton's method works perfectly when I solve the linear problem using Gaussian elimination, and the eigenvalues of the Jacobian are such that no severe numerical problems are faced with that algorithm.

I need to use an iterative method for larger problems, however, and the linear problem isn't appropriately solved with GMRES, IDR(s) or biconjugate gradient stabilized method when the matrix is taken as a whole, so I suppose another technique or precaution could be necessary.

I have tried preconditioning with an incomplete LU decomposition and with a diagonal preconditioner, and failed with both.
When individually solving for blocks $A$ and $D$, however, one notices that the linear system $A\Delta x = - f(x, y)$ exhibits the same behaviour as the complete Jacobian, while $D \Delta y = - g(x, y)$ is easily solved by any of the Krylov subspace algorithms I cited above.

Is there any algorithm that could help solve the linear system?


问:
Can we solve can overdetermined coupled equations in matlab?
```<code>for k=1:len
         thetaHattt(k)=atan(((Ahat+Rhat)*(1/k11(k)*Bhat-Chat)+Ahat*(Chat-1/k22(k)*Bhat))/((1/k22(k)*(Ahat+Rhat)*(1/k11(k)*Bhat-Chat)+1/k11(k)*Ahat*(Chat-1/k22(k)*Bhat))));
        if k<len&&k>1
            chatt(k)=Ahat/(k11(k-1)-k11(k))*(sin(thetaHattt(k-1))-sin(thetaHattt(k))-k11(k-1)*cos(thetaHattt(k-1))+k11(k)*cos(thetaHattt(k)));
            Ahatt(k)=(k11(k-1)-k11(k))*chatt(k)/((sin(thetaHattt(k-1))-sin(thetaHattt(k))-k11(k-1)*cos(thetaHattt(k-1))+k11(k)*cos(thetaHattt(k))));
            Bhatt(k)=k11(k)*(chatt(k)+Ahatt(k)*cos(thetaHattt(k)))-Ahatt(k)*sin(thetaHattt(k));
            Rhatt(k)=(k22(k)*chatt(k)+k22(k)*Ahatt(k)*cos(thetaHattt(k))-Bhatt(k)-Ahatt(k)*sin(thetaHattt(k)))/(sin(thetaHattt(k))-k22(k)*cos(thetaHattt(k)));
        end
    end
</code>```

In the above code,there are 5 unknowns Ahat,Bhat,Chat, Rhat and thetaHattt. All the five equations are given. K11 and K22 are vectors of length 10 and are known. The question is with the known K11 and K22 can we somhow find the values of the 5 unknowns given the constraint equations as above. chatt is same as Chat, Ahatt is same as Ahat. I am giving different variables just to show the equations. Any comment please.

Edit: One can further put a constraint that rate of change of thetaHatt with respect to the variable (k) is constant.


问:
Tolerance is not satisfied in eigsh function
I am diagonalizing a sparse matrix of size $153600\times 153600$ via eigsh function. I did not set the tolerance to specific value, so it would be the default value about $2\times 10^{-16}$. However, if I run the same code 3 times, the difference between eigenvalues are much bigger (about $6 \times 10^{-8}$) than this!
The eigenvalues for 3 trials are shown in the below figure. The real eigenvalue is $10^{-3}$ times the displayed one.

Question 1. Why the tolerance is not satisfied?

Question 2. How can I improve tolerance? Does it work if I just set the tol value smaller than machine epsilon? (I am worried that numpy module does not support high precision computation than machine epsilon.)

Question 3. The documentation for eigsh states that the iterative algorithm stops when eigenvalues reach within the tolerance limit. However, I think eigenvectors have a larger discrepancy for each trials than eigenvalues. Can I control the errors for eigenvectors?


问:
Does the time-dependent 1D advection-diffusion with point sources have an analytical solution?
I am looking for the analytical solution of 1-dimensional advection-diffusion equation with several point sources, Q, along the axial length of a cylinder through which the fluid flow occurs. Neumann boundary condition is specified at the outlet and Dirichlet boundary condition at the inlet.

$$\frac{\partial C}{\partial t} =  D\frac{\partial ^2 C}{\partial x^2}-v\frac{\partial C}{\partial x} + Q $$
with initial condition
$$c(x,0) = C_0$$
Neumann boundary condition at the right boundary
$$\frac{\partial C}{\partial x}=0\text{ at }t>0.$$

Dirichlet boundary condition at the left boundary
$$c(0,t) = C_L$$

The method to obtain a numerical solution for this problem is answered here and the Mathematica code is also posted.

I would like to compare the numerical solution with the analytical solution.
Are there some references in which I can find an analytical solution for this problem?
Could someone please suggest a reference or share provide an answer?


问:
Fast integration scheme for path integral of Gaussian over a cubic curve
I need to numerically compute an integral of the following form:
$$\int_0^1 \frac{1}{2\pi\sigma^2}\exp\left(-\frac{\|(q_0t^3 + q_1t^2 + q_2t + q_3) - a\|^2}{2\sigma^2}\right)\|3q_0t^2 + 2q_1t + q_2\|\,dt.$$
Here, $q_0,q_1,q_2,q_3,a$ are constant vectors in $\mathbb{R}^2$, and $\sigma$ is real. This comes from the line integral of a Gaussian over a cubic curve.

Obviously I can integrate this in Python using something like scipy.integrate.quad, but this is extremely slow, and this is an integral I need to compute many times over. Does anyone know of a fast way I can approximate an integral of this form?

One idea: if I plot the integrand for some particular values of $q_i$'s, I get something like this:

Of course, this looks like a Gaussian. But I'm unable to see the reformulation of the integral that would put it in the form of a Gaussian integral. Of course, if this really were a Gaussian integral in disguise, then computing the definite integral is easy with some approximations of the error function.

EDIT: This certainly cannot be reduced to a Gaussian integral. A particular choice of $q_i$'s and $a$ gives the following asymmetric plot of the integrand:


问:
Ritz projection error estimate for bilinear Lagrange element?
For second order elliptic problem
\begin{align}
-\Delta u=f,\quad in ~~\Omega\\
u=0,\quad on~~ \partial\Omega,
\end{align}
we have for the Ritz projection for $P_1$ conforming element
\begin{align}
\|u-R_hu\|_{W^{0,p}}+h\|u-R_hu\|_{W^{1,p}}\leq Ch^2\| u \|_{W^{2,p}}，\quad u\in W_0^{1,p}(\Omega)\cap W^{2,p}(\Omega).
\end{align}
Does this error estimate also hold for bilinear $Q_1$ element?


问:
Is today's enterprise computing sytem architecture with CPU, GPU and DPU is still Von Neumann architecture?
The current enterprise server rack has, in most cases, more than 1 CPU socket per board and a bunch of graphics card together with a data processing unit on a single board in a node. Several such nodes in a rack are closely stitched together for data pathways and tightly integrated into the cluster for a tailored application as deemed fit by the system architect. My question is Are we still in the Von-Neumann architecture era? The reason I raise this question is 2 fold

We no longer have a centralized processing unit, for processing a myriad of data, but rather separate discrete processing units other than CPU for handling specialized function including GPU, DPU and in some cases dedicated accelerators. Architecturally also, in some cases, GPU can have separate access to DRAM and SSD. So more than one processing unit (assuming entire socket as a processor which is not the case though). How come this classifies under Von-Neumann architecture?

If the so called neuromorphic computing ever comes into dominance would the current heterogenous compute elements would still be there or We will be completely shifting towards the neural engine cores? In other words would we still be using the traditional multiply and accumulate based registers or we will start moving to the compute in memory thing?

Are we really facing a architectural change to incorporate NVME drives, most importantly intel pmem,   into traditional architecture? Clearly pmem shows architectural choke points for utilizing its full potential

regards,
eigenvalue


问:
Discretizing Multi-species Ion Exchange Equations by Finite Volume Method

<p>
I'm solving a system of multispecies ion exchange equations (diffusion+drift fluxes) in 1-d spherical domain using finite volume method to obtain the ion concentrations at the next time step. After polishing the generic equation with so many assumptions, I ended up with this equation for species A
\[\frac{{\partial {C_A}}}{{\partial t}} = D\frac{{{\partial ^2}{C_A}}}{{\partial {r^2}}} - {\mu _{_A}}{C_A}\frac{{{D_A}\frac{{{\partial ^2}{C_A}}}{{\partial {r^2}}} + {D_B}\frac{{{\partial ^2}{C_B}}}{{\partial {r^2}}} + 2{D_C}\frac{{{\partial ^2}{C_C}}}{{\partial {r^2}}}}}{{{\mu _A}{C_A} + {\mu _B}{C_B} + 4{\mu _D}{C_C}}} + \frac{2}{r}{D_A}\frac{{\partial {C_A}}}{{\partial r}} - \frac{2}{r}{\mu _A}{C_A}\frac{{{D_A}\frac{{\partial {C_A}}}{{\partial r}} + {D_B}\frac{{\partial {C_B}}}{{\partial r}} + 2{D_C}\frac{{\partial {C_C}}}{{\partial r}}}}{{{\mu _A}{C_A} + {\mu _B}{C_B} + 4{\mu _D}{C_C}}}\]
</p>
 D is the diffusivity and mu is the ion mobility and they are all constants. There are two more similar equations for B and C. My question is: how can I discretize the second term in the RHS? should I assume that CA, CB and CD are known at initial time and discretize second order derivatives? knowing that the transient term in LHS can be solved using forward Euler's methods for time first order derivative (n is time index): </p> \[\frac{{\partial {C_A}}}{{\partial t}} = \frac{{{C_A}^{n + 1} - {C_A}^n}}{{\Delta t}}\]


问:
constrained zero-sum two person game
Finding the saddle point of a constrained zero-sum two-person game is equivalent to a resolution of primal-dual programs (with bi-linear objective function).

I am looking for a free solver to compute the saddle point in the constrained zero-sum two-person game (in primal-dual programs with linear constraints and bilinear objective function).

Primal-dual programs
$$
\begin{align}
\max_x \qquad &x^tAy \\
\text{s.t.} \qquad & Bx \leq b \\
& x_i\geq 0 \\
& x_1+x_2 +\dots+x_m = 1, \\ \\
\min_y \qquad &x^tAy \\
\text{s.t.} \qquad & Dy \leq d \\
& y_i\geq 0 \\
& y_1+y_2 +\dots+y_n = 1.
\end{align}
$$

Title: constrained zero-sum two person game


问:
Do the projection (with Jacobian) and marginalisation (inversion of matrix and remove a row/column and reinversion) commute?
I try to check the equality or the inequality between 2 Fisher matrices.

The goal is too see if the projection (with Jacobian) and marginalisation (inversion of matrix and remove a row/column and reinversion) commute.

Each of these 2 matrices is computed slightly differently. These 2 matrices are Fisher matrices.

Actually, this is the computation of changing parameters between initial parameters for each row/colum and final parameters for the final matrix. That's why in both computations, I am using the Jacobian J formulating the derivatives between initial and final parameters :

The formula is : $F_{\text{final}} = J^T F_{\text{initial}} J$

The first matrix has size 5x5 and the second one has 4x4 size. They are identical except the 4th row/column.

First : I inverse the 5x5 matrix (which gives a 5x5 covariance matrix). Then, I "marginalise", that is to say, I remove the 4th row/column of this covariance matrix. Then, I inverse again to get a 4x4 matrix.

Finally, I perform a projection with a Jacobian 4x4 with formula : $F_{\text{final}} = J^T F_{\text{initial}} J$ : so I get at the end a 4x4 matrix

For the second matrix to build : I am doing directly projection on 5x5 second matrix (which I recall is identical to the 5x5 except for the 4th row/column).

I perform this projection with the Jacobian 5x5. Then I get the second projected matrix 5x5. Finally, I remove the 4th row/column on this 5x5 matrix to get a 4x4 matrix new projected matrix.

I wonder under which conditions I could have equality between the 2 matrices 4x4. I don't know if my method is correct.

To show you a practical example, I put below a small Matlab script that tries to follow all the reasoning explained above :

```<code>clear;
clc;

% Big_31 Fisher :
FISH_Big_1_SYM = sym('sp_', [4,4], 'real');
% Force symmetry for Big_31
FISH_Big_1_SYM = tril(FISH_Big_1_SYM.') + triu(FISH_Big_1_SYM,1);

% Big_32 Fisher :
FISH_Big_2_SYM = sym('sp_', [5,5], 'real');
% Force symmetry for Big_32
FISH_Big_2_SYM = tril(FISH_Big_2_SYM.') + triu(FISH_Big_2_SYM,1);

% Jacobian 1
J_1_SYM = sym('j_', [4,4], 'real');
% Jacobian 2
J_2_SYM = sym('j_', [5,5], 'real');

% Remove 3th row/column
J_2_SYM(4,:) = [];
J_2_SYM(:,4) = [];

% Projection
FISH_proj_1 = J_1_SYM'*FISH_Big_1_SYM*J_1_SYM;
size(FISH_proj_1)

% Invert Fisher_2
COV_Big_2_SYM = inv(FISH_Big_2_SYM);

% Remove 3th row/column
COV_Big_2_SYM(5,:) = [];
COV_Big_2_SYM(:,5) = [];

% Re-inverse
FISH_Big_2_SYM_new = inv(COV_Big_2_SYM);

% Projection 2x2
FISH_proj_2 = J_2_SYM'*FISH_Big_2_SYM_new*J_2_SYM;
size(FISH_proj_2)

% Test equality between 2 matrices
isequal(FISH_proj_1,FISH_proj_2)
</code>```

The problem with this script is even I have small matrices (5x5 or 6x6), the code is very
slow with Matlab Symbolic at execution and I can't conclude up to now.

I hope to have been clear, don't hesitate me to ask for further information if my explanations are not enough precise.


问:
How to numerically solve PDE that governs the free vortex wake model?
Crossposted at Math SE

I am reading a paper on the free vortex wake model for a helicopter rotor blade, which is described by the following PDE:

$$\frac{\partial \vec{r}}{\partial \psi} (\psi, \zeta) + \frac{\partial \vec{r}}{\partial \zeta} (\psi, \zeta) = \frac{1}{\Omega} \vec{V} \left( \vec{r}(\psi, \zeta), t\right)$$

where $\zeta$ is the wake age angle, and $\psi$ is the azimuthal blade location (both measured in radians).

Also, it is known that $d \psi=d t \times \Omega$, where $\Omega $ is a constant that refers to the rotational speed of a rotor blade (rad/s).$V$ refers to the velocity, and $ r(\psi ,\zeta )$ is the location of a Lagrangian marker.

Since it is a technical paper, it only gives a brief summary of the method used, and I really have no idea how I should write a code to solve for $ r(\psi ,\zeta )$ with $V$. The paper uses a predictor-corrector 2nd backward scheme as shown here ($\frac{\partial \vec{r}(\psi, \zeta)}{\partial \psi}$ is computed using a 2-step backward Euler's method and $\frac{\partial \vec{r}(\psi, \zeta)}{\partial \zeta}$ is computed using central difference) :

\begin{array}{r}
\vec{r}^{m}\left(\psi+\Delta \psi_{i}, \zeta+\Delta \zeta\right)= \\
\left\{\frac{\vec{V}}{\Omega}-\left(-\frac{1}{2 \Delta \zeta}+\frac{23}{\phi}\right) \vec{r}^{m}\left(\psi+\Delta \psi_{i}, \zeta\right)-\left(\frac{1}{2 \Delta \zeta}-\frac{21}{\phi}\right) \vec{r}^{m}(\psi, \zeta+\Delta \zeta)\right. \\
+\left(\frac{1}{2 \Delta \zeta}+\frac{21}{\phi}\right) \vec{r}^{m}(\psi, \zeta)+\frac{3}{\phi} \vec{r}^{m}\left(\psi-\Delta \psi_{i-1}, \zeta+\Delta \zeta\right)+\frac{3}{\phi} \vec{r}^{m}\left(\psi-\Delta \psi_{i-1}, \zeta\right) \\
\left.-\frac{1}{\phi} \vec{r}^{m}\left(\psi-\Delta \psi_{i-1}-\Delta \psi_{i-2}, \zeta+\Delta \zeta\right)-\frac{1}{\phi} \vec{r}^{m}\left(\psi-\Delta \psi_{i-1}-\Delta \psi_{i-2}, \zeta\right)\right\}\left\{\frac{1}{2 \Delta \zeta}+\frac{23}{\phi}\right\}^{-1}
\end{array}

and
\begin{aligned}
\vec{V}=& 4 V_{\infty}+V_{i n d}\left(\vec{r}^{m-1}(\psi, \zeta)\right)+V_{i n d}\left(\vec{r}^{m-1}(\psi+\Delta \psi, \zeta)\right) \\
&+V_{i n d}\left(\vec{r}^{m-1}(\psi, \zeta+\Delta \zeta)\right)+V_{i n d}\left(\vec{r}^{m-1}(\psi+\Delta \psi, \zeta+\Delta \zeta)\right) \\
\phi=& 46 \Delta \psi_{i}+4 \Delta \psi_{i-1}-2 \Delta \psi_{i-2}
\end{aligned}

V_ind has to be computed by integrating the Bio-Savart equation:
$d \vec{V}_{k}=F_{\nu} \frac{\Gamma}{4 \pi} \frac{d \vec{l} \times \vec{r}}{|\vec{r}|^{3}}$

$\Gamma$ is a constant value referring to circulation and $F_{\nu}$ is a correction factor preventing singularity.

$F_{\nu}=\frac{|\vec{r}|^{2}}{\left(|\vec{r}|^{2 n}+r_{c}^{2 n}\right)^{1 / n}}$, n is usually taken to be 2 and $r_{c}$ is a generic core radius (meters).

Could someone give me some guidance on how I can solve for $r$ ins the original PDE by writing code in Python or in MATLAB? Also, since $\psi$ is a function of time, is this an initial boundary value problem or just a regular boundary value problem?


问:
Confusion about preconditioner for incompressible Navier-Stokes equation with implicit-explicit method
Consider the time-dependent Navier-Stokes equation

$$u_t + (u \cdot \nabla) u - \Delta u  + \nabla p = f$$
$$\operatorname{div}(u)=0$$

Looking at deal.ii tutorials, I've notice that there are extensive discussions about the choice of a preconditioner in the time independent case, while for the time dependent one I've only found this one the code gallery which fits quite well my testcase

Applying an IMEX scheme, I got

$$
\frac{M}{\Delta_t} U_{n+1} + A U_{n+1} + B^t P_{n+1} = G(t_n,U^n,f)
$$
$$
B U_{n+1} = 0
$$

Therefore at each time-step I need to solve the saddle point problem

\begin{bmatrix}
\ \frac{M}{\Delta_t} + A & B^t\\B & 0
\end{bmatrix}

Of course, I need to apply a preconditioner and, following the discussion, I'd chose

$$P = \begin{bmatrix} \frac{M}{\Delta_t} + A & B^t \\ 0 & -S^{-1} \end{bmatrix}$$ where $S$ is the Schur complement of $\frac{M}{\Delta_t}+ A$

Now, let's go to the hot-spot: I need to understand what can be a good way to invert the top left block and the bottom one $S$.

Since the top left block $\frac{M}{\Delta_t} +  A$ is symmetric and positive definite, then in order to "invert" I would just use Conjugate Gradient.

For $S$ I would use the pressure mass matrix, which is usually called $M_p$ in deal.ii website. Essentially, it is composed by $L^2$ products of basis functions for the pressure. This is symmetric and positive definite too, so I can use CG again.


问:
How to treat nonlinear radiation term in heat equation using Finite-element method?
I am trying to solve the time-dependant non-linear heat equation with radiation. This equation is coupled to the radiative transfer equation but for the purpose of my question, this does not matter much.
I am using a Crank-Nicolson scheme for the time stepping scheme as well as a Taylor Development of $T^4$.
The Taylor Formula for timestep $n+1$ reads:

$$
T_{n+1}^4 \approx T_{n}^4 + 4T_{n}^3(T_{n+1}-T_{n})
$$

and using the Crank-Nicolson scheme I end up with a term like this in my equation
$$
\frac{T_{n+1}^4+T_{n}^4}{2}
$$
which simplifies with the Taylor's Formula in :

$$
\frac{T_{n+1}^4+T_{n}^4}{2} =-T_{n}^4+2T_{n}^3T_{n+1}
$$

Now I am trying to solve my equation using a finite element method, so over a cell $K$ of my mesh. I can write my temperature as a linear combinaison of trial function $w \in H^2(K)$, such that $T_{n} = \sum_{n_{dof}} T_{n,i}w_i$, with $n_{dof}$ the number of degrees of freedom on my cell $K$. By producting the weak form of my equation, I have to calculate terms like :

for $1 \leq i \leq n_{dof}$ :
$$
\int_K T_{n}^3T_{n+1}w_i dK
$$

Now my problem is how to calculate properly this term. If I replace $T_{n}$ by its linear combinaison of trial function I provided, I end up with

$$
\sum_{n} T_{n+1,n} \int_{K} \left(\sum_{k,l,m} T_{n,k}T_{n,l}T_{n,m}w_kw_lw_mw_nw_i\right) dK
$$

and so if I do that for $1 \leq i \leq n_{dof}$, I end up with a linear algebra expression $MT_{n+1}$, with

$$
M = \left( \int_{K} \left(\sum_{k,l,m} T_{n,k}T_{n,l}T_{n,m}w_kw_lw_mw_nw_i\right) dK\right)_{1 \leq n,i \leq n_{dof}}
$$

Which is quite ugly to work with, but is it the only way to approach this term $T_n^3T_{n+1}$ ?

The only other method I could think of is using an orthonormal basis for $w_i$'s such that $T_{n}^3 = \sum_i^{n_{dof}} T_{n,i}^3w_i$. If I do so the matrix M is like this :

$$
M = \left( \int_{K} \left(\sum_{k} T_{n,k}^3w_kw_nw_i\right) dK\right)_{1 \leq n,i \leq n_{dof}}
$$

Which looks more manageable, but I am not quite sure this is optimal. So my question is : in the finite element method framework, what would be the best way to treat a term like $T_{n}^3T_{n+1}$ ?


问:
PLASMA usage (Linear algebra routines that supports multithreading)
I have been looking for linear algebra libraries that support multithreading. I have found PLASMA which looks promising. It is from the same group that developed LAPACK.

http://icl.cs.utk.edu/projectsfiles/plasma/html/InstallationGuide.html

https://icl.utk.edu/projectsfiles/plasma/pdf/users_guide.pdf

However, it doesn't seem to be very well maintained/documented which is surprising. For instance, I was unable to find the user-forum to ask questions.

I have been able to install the library correctly. I created a mybuild directory and used cmake to install.

debasish@debasish-HP-Z840:~/Desktop/softwares/plasma-20.9.20/mybuild$ ls

CMakeCache.txt  CMakeFiles  cmake_install.cmake  libplasma_core_blas.so  libplasma.so  Makefile  plasmatest

Next, I am able to run tests as following:

debasish@debasish-HP-Z840:~/Desktop/softwares/plasma-20.9.20/mybuild$ ./plasmatest dgesv

Status      Error       Time    Gflop/s       n    nrhs    nb    ib  padA  padB  mtpf

pass   1.28e-19     0.1492    17.8622    1000    1000   256    64     0     0     1

However, I am unsure how to use these libraries in my own code.

Sorry, in the absence of any documentation I am unable to be more specific but I thought I will give it a try here. In general, I am not good at using libraries without comprehensive documentation.


问:
TVD Lax-Wendroff with non-constant velocity
I am dealing with a linear advection equation with a non-constant velocity, where I would like to apply a TVD Lax-Wendroff scheme in 1D.

The equation is the following:

\begin{equation}
    \frac{\partial \: C(x, t)}{\partial t} =  -  \frac{\partial\:(\vec{v(x,t)} \:C(x, t))}{\partial x}
\end{equation}

Which describes the advection of a concentration by a fluid, with $C$, the concentration, and $\vec{v}$, the fluid velocity.

I understand how TVD Lax-Wendroff works in 1D for a linear advection equation, where I can write the flux like a combinaison of the upwind and Lax-Wendroff fluxes for a constant velocity (see (1) or (2) p177, for example).

However, I am currently stuck on how to apply this method with a non-constant velocity.

Using a flux conservative form, in 1D, I can write, with $v_x C = F$:

\begin{equation}
\frac{C^{n+1}_i - C^{n}_i}{\Delta t} = - \frac{1}{\Delta x}(F^{n}_{i+1/2} - F^{n}_{i-1/2})
\end{equation}

for a constant and positive v, the flux can be expressed as:

\begin{equation}
F^n_{i+1/2} = v u_i + \frac{\phi(r_{i+1/2})}{2} (v - \frac{v^2\Delta t}{\Delta x}) (u_{i+1} - u_{i})
\end{equation}

with $\phi$ beeing the flux limiter function.

The first term is the upwind flux, and the second term is the anti-diffusive Lax-Wendroff flux.

My idea is to use staggered velocities.

I can easily replace $v u_i$ by $v_{i+1/2} u_i$ but I don't know how to discretize the velocity in the anti-diffusive term, as I don't quite understand how it was derived.

I would write something like that for the right flux:

\begin{equation}
F^n_{i+1/2} = v_{i+1/2} u_i + \frac{\phi(r_{i+1/2})}{2} (v_{i+1/2} - \frac{v_{i+1/2}^2\Delta t}{\Delta x}) (u_{i+1} - u_{i})
\end{equation}

Does it make sense for you?

It doesn't seem to work, as it creates lots of diffusion when I am trying to implement it.


问:
Is there a graphical interpretation or explanation of automatic differentiation compared to numerical differentiation
I have been looking at automatic differentiation for solving differential equations lately. I understand the basic ideas of using Dual numbers and such for finding derivatives, etc. However, I feel like I am not understanding the geometric interpretation of automatic differentiation and dual numbers.

Let's consider a simple differential question or system of differential equations.

$$
\frac{dy}{dt} = f(y, t)
$$

So let me set a little context first for how I understand the application of automatic differentiation to numerically solving differential equations. If I make any mistakes, please correct me. So if you have a relatively low dimensional system with reasonable time-scale separation, then you could get away with some second order explicit or RK type solvers. That should provide good accuracy and fast compute time. In such cases, you don't need use automatic differentiation, because you can compute the derivative f(y,t) directly--if the solver is explicit.

It seems that AD become more relevant for larger systems of equations where there is less time-scale separate and more stiffness. In those cases, implicit solvers are the best approach to take. But when using an implicit solver you have to solve a root finding problem using something like Newton's method or Anderson Acceleration, etc.

So in Newton's method we have to compute the jacobian of the original differential equation $f(y,t)$. AD become relevant or useful for computing the jacobian of the original RHS of the differential equation when finding the roots for the implicit equation.

It seems like the older approach was to use numerical differentiation or symbolic methods to represent the jacobian in the root finding implementation. And in the case of numerical differentiation, this is susceptible to floating point errors and is generally ill conditioned. And in the symbolic case, this is generally slow to solve.

So that is my understanding of why AD is used.

BUT, in many of the explanation of AD, the authors will explain the reason for using AD, but then will dive deep into algebraic constructions of Dual numbers and Dual number arithmetic. So I don't have a sense of the geometric intuition behind dual numbers. I get the algebraic proofs. I understand that Dual number arithmetic is an application of the chain rule. But I was hoping that someone could explain some sort of graphical intuition behind what dual numbers are and why they can solve for derivatives and jacobians. Is there some graphical explanation for what a number represented by a value, and derivative are? How do these compose graphically to produce the overall derivative of the function?

Any help is appreciated.


问:
Can a symmetric matrix be found using matrix-vector operation while maintaining matrix-vector operatin?
Consider $A\in\mathbb{R}^{n\times n}$ its not a special matrix and in the worst case all of its entries are non-zero. I am looking for a way to compute $AA^{T}$ using matrix-vector operation. The total number of flops are $(2n-1)\frac{n(n+1)}{2}$ because all I have to do for a symmetric matrix like $C=AA^{T}$ is to compute the diagonal entires which are $C(i,i)=A(i,:)^TA(i,:)$. What I want now is to compute the lower triangular part and then when I am done I just say that the upper triangular part is the same as the lower triangular part. The problem is can I do it in matrix-vector multiplication or will it force me to perform unnecessary multiplication (like multiplying elements in the upper part)? It is clear that a vector-vector computation would work but I am interested to know if a matrix-vector computation would work or not.


问:
Procedure to convert continuous equations of motion to discrete version
Let's take a mobile robot and suppose we know its continuous equations of motion, for example this car-like simple model.

Now if I am simulating this robot in a continuous 2D plane, then coding the simulation is simple. But suppose i am now in a 2D grid world. Movement is only from one cell to another cell. Now how do i simulate this motion given that i only know its continuous version.

I remember seeing purely discrete equations like those for a cellular automata are derived once and for all, and then these purely discrete equations are naturally used in the grid world just like continuous equations are used naturally for the continuous world.

Is there a good example that teaches how we can go from ODEs to CA-like equations?


问:
How to solve nonlinear second order ODE in Matlab?
I am working on simulating a car suspension system using Matlab.

Specifically, I have to derive equation of motion using the Lagrange method and then use ode 45 to solve it. However, while using odetoVectorField I keep getting the error: System contains a nonlinear equation in diff(theta_rod_1(t), t). The system must be quasi-linear: highest derivatives must enter the differential equations linearly.

Is there anyway to get around this?

Here is part of the code:

```<code>%Defining relations between variables

th_M = asin( ( L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2)) / (LMx - 2*(L1+L2) ) );

L_spring_1 = sqrt(L_rod^2 + L2^2 - 2*L_rod*L2*cos(pi/2 + th_M - th_rod_1) ) ; %%Law of cosine

th_spring_1 = pi/2 - th_M - acos((L2^2 - L_rod^2 + L_spring_1^2)/(2*L2*L_spring_1));

L_spring_2 = sqrt(L_rod^2 + L2^2 - 2*L_rod*L2*cos(pi/2 - th_M - th_rod_2) );

th_spring_2 = pi/2 + th_M - acos((L2^2 - L_rod^2 + L_spring_2^2)/(2*L2*L_spring_2));

dx = LMx/2 * th_M;

dth_M = (L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2)/((1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx));

ddx = (LMx*(L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2))/(2*(1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx));

%% Kinematics

% Define the displacements of the double spring in Cartesian coordinates

xm1 = (x0+L1+L2)*cos(th_M) - L_rod*sin(th_rod_1); %x-pos of the left wheel, the cos(th_M) projects the length onto the x-direction

ym1 = 0;

xm2 = (x0+LMx-L1-L2)*cos(th_M) + L_rod*sin(th_rod_2);
ym2 = 0;

xM = x0;

yM = L_rod*cos(th_rod_1) + LMy/2*cos(th_M) - (LMx/2 -L2-L1)*sin(th_M); %yM is defined as the height of mass M center of gravity

%Find the velocities by differentiating the displacements with respect to
%time (using diff function)

vxm1 = ((L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))*(L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2)*(L1 + L2 + x0))/((1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx)^2) - L_rod*cos(th_rod_1)*Dth_rod_1
vym1 = 0;

vxm2 = L_rod*cos(th_rod_2)*Dth_rod_2 - ((L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))*(L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2)*(L1 + L2 - LMx - x0))/((1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx)^2);

vym2 = 0

vxM=0;

vyM = ((L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2)*(L1 + L2 - LMx/2))/(2*L1 + 2*L2 - LMx) - L_rod*sin(th_rod_1)*Dth_rod_1 + (LMy*(L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))*(L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2))/(2*(1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx)^2);

%% Lagrange

% Kinetic energy

T = 0.5*M*ddx^2 + 0.5*J0*dth_M^2 + 0.5*m1*vxm1^2 + 0.5*m2*vxm2^2 ;

%

T = simplify(T);

% Potential energy

V = m1*g*ym1 + m2*g*ym2 + M*g*yM + 0.5*k1*(L_spring_1-L_spring_rest)^2 + 0.5*k2*(L_spring_2-L_spring_rest)^2 ;

V = simplify (V);

%Dissipation energy

R = 0.5*b*(vxm1*sin(th_spring_1))^2 + 0.5*b*(vxm2*sin(th_spring_2))^2;

Lagrange = T-V;

%Generalized coordinates:

q   = [th_rod_1, th_rod_2];

dq  = [Dth_rod_1, Dth_rod_2];

ddq = [DDth_rod_2, DDth_rod_2];

DL_Dq = jacobian(Lagrange,q')';  % Operator ' is transpose

DL_Ddq = jacobian(Lagrange,dq);

DDL_DtDdq = jacobian(DL_Ddq',[q, dq]) * [dq, ddq]';

DR_Ddq = jacobian(R,dq')';    %Damping term of Lagrange_Euler equation

EoM = DL_Dq - DDL_DtDdq    ;  %Euler-Lagrange equation with damping

%- DR_Ddq

EoM = EoM == 0;

%simplify(EoM(1));

%%
syms theta_rod_1(t) theta_rod_2(t)

EoM=subs(EoM,[th_rod_1, Dth_rod_1, DDth_rod_1, th_rod_2, Dth_rod_2, DDth_rod_2], ...
[theta_rod_1, diff(theta_rod_1), diff(theta_rod_1,2), theta_rod_2, diff(theta_rod_2), diff(theta_rod_2,2)]);

%% Compute with value

% m1 m2 M g k1 k2 b J0  x0 L1 L2 LMy LMx L_rod L_spring_rest

x0=0;
L1 = 1;
L2 = 3;
LMx = 10;
LMy = 5;
L_spring_rest = 2;
L_rod = 3;
m1 = 1;
m2 = 1;
M = 5;
g = 9.8;
k1 = 500;
k2 = 700;
b = 100;
J0 = 50;

eqn1 = subs(EoM(1))

eqn2 = subs(EoM(2));

%The two equations are nonlinear second-order differential equations.

%To solve these equations, convert them to first-order differential equations

[V,S] = odeToVectorField(eqn1,eqn2);
</code>```

OwnerDisplayName: Luận Trần Công


问:
Maximizing $l_1$-normalized entropy using CVXPY
Suppose that $x = (x_1, ..., x_n)$ is a vector of variables and I would like to maximize the Shannon entropy of $\frac{|x|}{||x||_1}$ (i.e. the vector of absolute values of $x_i$, normalized to have $l_1$ norm = 1) subject to some linear constraints. However, this function is not directly expressible as a composition of atoms recognized by CVXPY as convex/concave (although CVXPY recognizes entropy itself). Are there any tricks to write this program in a way suitable for CVXPY?


问:
Efficient solver of a Integer programming
I am solving an Integer programming using MATLAB, yet the efficiency is low.
Here is the problem:

Suppose $v$ is a $N \times 1$ vector. For $v_i \in v$, $v_i \in \{0,1\}$.
$D$ is a 0-1 matrix, which means for each $d_{ij} \in D$, $d_{ij} \in \{0,1\}$.

The translation of this problem is to make element $1$ in $v$ to be  as few as possible when the constraints are met.
\begin{equation}
\begin{split}
&\min \left\|\boldsymbol{v}\right\|_2 \\
&s.t.\  \begin{array}{c}
D\boldsymbol{v} \ge \boldsymbol{1}\\
\end{array}
\end{split}
\label{eq:core_op}
\end{equation}

Is there a good solver for this kind of problem in the field of optimization?


问:
Understanding inf-sup conditions for classical saddle point problems
I'm studying the inf-sup conditions for saddle point problems. I'm referring to the usual one $$\begin{cases}Au + B^t p = f \\Bu=g \end{cases}$$ In the book I'm using (Ern - Guermond: Theory and practice of FEM) they say (pg. 100)

The problem is well posed iff

$$ \exists \alpha >0: \inf_{u \in \ker(B)} \sup_{v \in \ker(B)} \frac{a(u,v)}{||u||_X ||v||_X} \geq \alpha$$

$$\forall v \in \ker(B), \forall u \in \ker(B), a(u,v)=0 \implies v=0$$
[... and then we have also the inf-sup on $B$, which is not a problem for me right now]

I understand how they derive them, but then, looking at my lecture notes, my professor replaced the second condition with another inf-sup

$$ \exists \alpha >0: \inf_{v \in \ker(B)} \sup_{u \in \ker(B)} \frac{a(u,v)}{||u||_X ||v||_X} \geq \alpha$$

I've searched in lots of books and lecuture notes, but I still cannot understand why this one above is equivalent to the second condition. In particular, it's clear to me that the last "inf-sup" must be strictly greater than $0$, but I don't see why the constant $\alpha$ must be the same.

What am I missing?


问:
Invert a huge sparse operator;
please help me with this question, I want to invert a huge sparse (non-circulant) this below in a $Ax=y$ equation:

$$(\lambda I+ \beta D+ \sigma C)x=y$$
where
I is an Identity Matrix,D is a Diagonal Matrix,C is a circulant Matrix and $\lambda$, $\beta$, $\sigma $ are positive constants.

And all have the same size of 10000x10000 (in fact it could be much larger because I calculate for processing an image the size of 1000 x 1000, so my matrix will be 1000000x1000000).

I did with scipy solver and the inversion is applicable only for image size of 100x100, which means 100000x100000 matrix.

Title: Invert a huge sparse operator;


问:
About the condition $\ker(B_h) \subset \ker(B)$ in mixed finite elements formulation
I'm studying mixed finite elements. The problem is a classical saddle-point one: we seek for $(u,p)$ in $V \times Q$:
$$A u + B^t p = f$$
$$Bu = g$$
where $A: V \rightarrow V', B:V \rightarrow Q'$  and $V,Q$ are Hilbert spaces. For this problem, we have the inf-sup conditions on $A$ and $B$ that essentially tell us that the problem is well posed.

As far as I understand, the issue one the discrete level (in contrast to what happens for Poisson equation, and other coercive problems) is that if we consider $V_h \subset V$ and $Q_h \subset Q$ as finite dimensional spaces in which we seek the solution, then we can no more apply the continuous properties to the discrete problem because for instance $$\ker(B_h) \not \subset \ker(B)$$ in general.

So far so good. Then, while doing the analysis, my professor wrote the following diagriam

and he said that we would like to join with dotted arrow and we would like that $B_h$ maps entire $V_h$ into the entire $Q_h'$ and we would like this diagram to commute. He then said that it's trivial if $\ker(B_h) \subset \ker(B)$.

I really can't understand why the inclusion of the kernels implies the commutativity of the diagram. Any help is highly appreciated, and even an example where $B$ is the usual divergence operator is perfectly fine to me.

EDIT: Maybe I got the way to go: the diagram commuting means that $$ \Pi_{Q_h'} \circ B =  B_h  \circ \Pi_{V_h}$$

Now, if the $\ker(B_h) \not \subset \ker(B)$, this means there exists a function $v_h \in V_h \subset V$ such that $$B_h v_h = 0$$ and $$B v_h \ne 0$$

If it commutes, it should hold:$$ B_h \Pi_{V_h} v_h=\underbrace{B_h v_h}_{0}  = \Pi_{Q_h'} \underbrace{B v_h}_{\ne 0} \ne 0$$

which is clearly false. In the first equality I used the fact that $v_h$ is of course not affected by the projection. What do you think ? I think that something more can surely be said. Please feel free to add more things about this


问:
How to obtain smallest eigenvalues with Arnoldi iteration
I understand that the Arnoldi iteration produces a basis which tends to include in its span the eigenvectors corresponding to eigenvalues of large magnitude (hence the analogy between the last vector of the Krylov subspace, $A^{n - 1} b$, and the result of $n - 1$ iterations of power iteration method made in this Wikipedia article).

However, most solver implementations (for example this) allow for the user to choose between obtaining the eigenvalues of largest or smallest magnitude, or even the ones of smallest real part.

How would that be implemented?


问:
Finite Difference Approximation for the Laplacian in 2D that produces a nonsymmetric matrix
Consider the following PDE

\begin{align}
-\Delta u  &= f \ \ \text{en} \ \ (0,1)\times (0,1) \label{P1}  \\
u &= 0  \ \ \text{en} \ \partial ((0,1)\times (0,1)) \label{P2}
\end{align}

if we use 5-point finite difference approximations we get a system
\begin{equation}
Ax = b  \label{Sist_Planteamiento}
\end{equation}

where the matrix $A$ is symmetric.

Questions:

Do you know a finite difference approximation that produces a system

\begin{equation}
Ax = b
\end{equation}
where $A$ is nonsymmetric?

Do you know books or papers where I can study systems like the system in question 1?

Thanks in advance


问:
Solving multiple linear regression in parallel
I am working on a problem where I need to solve approximately 500 Million Linear Regressions (OLS).
What would be the most efficient way to do this (e.g. using GPU or a some framework that can do this efficiently)?

Note that in my case each individual OLS problem is "small" in size e.g. each input dataset has 100 observations with (at most) 5 features.

Here's what I have tried:

I use a 25 core machine and parallelize my computation across all cores using Python3 multiprocessing module. To solve OLS, I use numpy.linalg.solve as I found that to be more efficient than sklearn.linear_model.LinearRegression.


问:
How to determine the orientation of convex/concave hexahedra?
I am writing a code that checks the orientation of a list of vertices (along with face connectivity) describing both convex and concave hexahedra. The face connectivity table stores the list of vertex indices that constitute each face, as described on p.27 of this document. I flag a hexahedron as being positively oriented if all faces are oriented (i.e. clockwise or counter-clockwise) such that the face normal points outward. The direction of the face normal can be determined by applying the right-hand rule in the plane of the face: four fingers curl in the direction of the face vertex orientation, and the thumb points in the direction of the face normal.

I also store the vertex list in the same order as in the linked document. However, the issue I have is that there are certain geometric operations (e.g. reflection about a plane) in my code that flip the orientation of a given hexahedron. I am currently able to check the orientation robustly for convex hexes by computing the centroid and ensuring that the face normals are outward-pointing with respect to it. However, this is not a robust procedure in the case of concave hexes because - (i) the centroid may lie outside the hex, or (ii) the interior side of a face may not be 'visible' to the centroid.

Can anyone suggest a robust procedure for checking the orientation of arbitrary hexahedra? Any advice would be greatly appreciated.


问:
Getting euclidean distance between vector A and C without anyway of retrieving them when their distances with a common vector B is known
Motivation:

My plan is to get the overall euclidean distance matrix for all the vectors in N number of dataset. Each dataset is basically an array of n-dimensional points. For e.g: A dataset can be like [[1,2,3,...n],[1,2,3,..n],...[1,2,3,n]]. However, the datasets are not gonna be shared to a single entity for which I cannot compile them and thus will not know all the points to calculate pairwise euclidean distance for some similarity calculation and clustering analysis and modeling. However, I will be only notified what are the common vectors in the datasets and distance of the vectors in any dataset with respect to those common vectors without knowing or passing the point's coordinates. If the dataset could be shared in the same model, I wouldn't have to face this distributed calculation problem.

Details:

A and C are two n-dimensional vectors from two different dataset. They have a common vector B. I want to calculate euclidean distance between A and C without exposing A and C, rather using B. Lets say ED(A,B) is calculated in model_1 and ED(B,C) is calculated in model_2, I can easily use ED(A,B) i.e AB and ED(B,C) i.e BC to compute ED(A,C). However, if someone know B they can easily find out A and C. Even if B is a randomly created common point, if it available, A and C can also be found. Is it possible to use B such that A and C cannot be retrieve in anyway. I have looked into differential privacy, PATE, SMPC techniques but they have their limitations in preserving privacy.

Work done so far :

https://www.online-python.com/Pq7OhM8FJe

Question

Given the euclidean distance between vector A and vector B, and that between vector B and vector C, how to calculate the euclidean-distance between point A and C such that there is no way to retrieve A and C even if B is known?

Any computation technique I can follow?

Thanks in advance.


问:
Efficient multidimensional numerical integration in R and C++
I'm trying to perform a 4-dimensional numerical integration in R using a function I wrote in C++ code which is then sourced in R using the Rcpp package.

Below there is my code:

```<code>// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#define EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
#include <Eigen/Eigen>
#include <Rcpp.h>
#include <math.h>
#include <iostream>
#include <cstdlib>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/gamma.hpp>

using namespace std;
using namespace Rcpp;
using namespace Numer;

// [[Rcpp::depends(BH)]]

// [[Rcpp::export]]
double gaussian_free_diffusion(double x,double x0, double sigma, double t) {
  double pi = 2 * acos(0.0);
  double a1 = (1/sqrt(2.0 * pi * sigma * t));
  double b1 = exp(-pow((x - x0), 2.0)/(2.0 * sigma * t));
  double res = a1 * b1;
  return res;
}

// [[Rcpp::export]]
double integral_CppFunctionCUBATURE(NumericVector Zt_pos, const double& xA0, const double& xB0, const double& yA0,
   const double& yB0, const double& t1, const double& sigma){

  double xAt_pos = Zt_pos[0];
  double xBt_pos = Zt_pos[1];
  double yAt_pos = Zt_pos[2];
  double yBt_pos = Zt_pos[3];

  double temp_pbxA = gaussian_free_diffusion(xAt_pos, xA0, sigma,t1);
  double temp_pbxB = gaussian_free_diffusion(xBt_pos, xB0, sigma, t1);
  double temp_pbyA = gaussian_free_diffusion(yAt_pos, yA0, sigma,t1);
  double temp_pbyB = gaussian_free_diffusion(yBt_pos, yB0, sigma, t1);

  return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);

};
</code>```

integral_CppFunctionCUBATURE is the function I use with cubintegrate. In R, I would then do:

sourceCpp('./myCppcode.cpp')

```<code>xA0<-3
xB0<-2
yA0<-5
yB0<-10
t<-500
sigma<-1

cubature::cubintegrate(integral_CppFunctionCUBATURE,lower=rep(-1000,4), upper=rep(1000,4),
                       xA0=xA0, xB0=xB0, yA0=yA0,yB0=yB0, t1=t, sigma=sigma,method='cuhre')$integral
</code>```

Result: 0.9999978

The multidimensional integration takes about 4 sec, but I would like to speed it up as much as possible, since when I run the multidimensional integration on several xA0 values this can take quite some time. Do you have any suggestion on how I could improve the speed of the code?

And also, would you suggest an alternative way to perform fast 4-dimensional integration in C++/Rcpp?


问:
Plotting the motion of a positive charge in a cylindrically symmetric magnetic field
I want to plot the motion of a positive charge in a cylindrically symmetric magnetic field.

I am assuming a cylinder around the z-axis, with the magnetic field going in clockwise direction. The B-field has magnitude of 6 T and the distance R from the z-axis is 3 m. The charged particle is launched in positive direction along the z-axis and has the energy 2 MeV.

I am uncertain of how to simulate this B-field correctly. I was thinking to create the B-field in cylindrical coordinates, cylinder from 0 to 2pi:

theta = numpy.linspace(0, 2*numpy.pi, 360)
x = R*numpy.cos(theta)
y = R*numpy.sin(theta)

Bx = B0*(numpy.cos(numpy.arctan2(y,x)
By = B0*(-numpy.sin(numpy.arctan2(y,x)))
Bz = 0
```

And then create a vector B=[Bx, By, Bz] from which I would calculate the acceleration using Lorentz force for a timespan t.

But I think I am going in circles with this. Is there another way to create a cylindrically symmetric magnetic field?


问:
Binarization for optimization problems
I have a nonlinear mixed-integer optimization problem, and because of very high complexity when solving it using methods like Branch and Bound, I resorted to solve it using alternating method and continuous relaxation of the binary variables. So the issue is that the resulting optimum values are fractional. How to properly get binary values from the optimum fractional variables, such that the constraints in the problem remain satisfied?

For example, if $x$ is originally a binary variable $\{0,1\}$ and is linearly relaxed to be $0 \leq x \leq 1$, and I solved the optimization problem based on this relaxed constraint. If the problem is finished and the optimum solution is in the continuous range $[0,1]$ not binary, and I want to properly binarize it, keeping the problem still feasible (i.e. the binarization should not affect the problem feasibility).
The Constraints are that in the attached image

Title: Binarization for optimization problems


问:
Numerical Soultion to Background equations of cosmology
I am trying to solve the background equations of cosmology numerically using Runge-Kutta Dormand Prince method with simplified assumption $8\pi G=1$ and $c=1$. The equations are
$$\ddot a = - \frac{1}{6}(\rho + 3p) a $$

and
$$\dot {\rho} = -3 \frac {\dot a}{a}(\rho +p) $$

for matter$(p=0)$ and radiation $(p=\frac{1}{3}\rho)$.

I should get $a\propto t^{2/3}$ for matter and $a\propto t^{1/2}$ for radiation. But I am getting same variation for scale factor $a$ in both the cases. I am a beginner. I am unable to figure out where I have gone wrong . Please help me.

The equations I pass to my code are as follows

For matter

```<code>       void fcn(double t, double *y, double *f){

                double p,hub;
                p=0;
                hub=y[1]/y[0];
                f[0] = y[1];
                f[1] = -(y[0]*y[2]+3*p*y[0])/6.0;
                f[2] = -3*hub*y[2]-3*hub*p;

             }
</code>```

For radiation

```<code>           void fcn(double t, double *y, double *f){
           /* a=y[0],adot=y[1],rho=y[2]  */
                     double p,hub;
                     hub=y[1]/y[0];
                     p= (2.0/3.0)*y[2];
                     f[0] = y[1];
                     f[1] = -(y[0]*y[2]+3*p*y[0])/6.0;
                     f[2] = -3*hub*y[2]-3*hub*p;
                 }
</code>```

The plot I have obtained for $a$ vs $t$ is
<img src="https://i.stack.imgur.com/PeRnx.png" alt="$a$ vs $t$ plot" />

I am just trying to see what is the type of functional dependence we get , not putting the initial conditions as per cosmology. Can this be because of wrong initial conditions?

One more thing is that variation of $\rho$ with time should be same in both cases ($\rho \propto t^-2$) but I am getting different result for this also.

My plot for $\rho$ vs $t$ is as below

<img src="https://i.stack.imgur.com/T7emw.png" alt="$\rho$ vs $t$ plot" />

My code was just working excellent to solve other coupled differential equations.

```<code>         #include<stdio.h>
         #include<math.h>
         #include<stdbool.h>
         #include  <stdlib.h>

         #define N 3

         double ddopri5(void fcn(double, double *, double *),
                        double *y);
         double alpha;
         void fcn(double t, double *y, double *f);
         double eps;

         int main(void){

               double y[N];
               //eps = 1.e-9;
               printf("Enter epsilon:\n");
               scanf("%lg", &eps);
               y[0]=1.0;
               y[1]=12.00;
               y[2]=30.00000567845;
               ddopri5(fcn, y);

          }

         void fcn(double t, double *y, double *f){
                double p,hub;
                 p=0;
                 hub=y[1]/y[0];
                 f[0] = y[1];
                 f[1] = -(y[0]*y[2]+3*p*y[0])/6.0;
                 f[2] = -3*hub*y[2]-3*hub*p;

          }

           double ddopri5(void fcn(double, double *, double *),
                         double *y){
                 double t, h, a, b, tw, chi;
                 double w[N], k1[N], k2[N], k3[N], k4[N], k5[N],
                        k6[N],k7[N], err[N], dy[N];
                 int i;
                 double errabs;
                 int iteration;
                 FILE *fpw;
                 fpw=fopen("case1.dat", "w");

                 iteration = 0;
                 //eps = 1.e-9;
                  h = 1.0e-6;
                  a = 0.0;
                  b = 30;
                  t = a;
                 while(t < b -eps){
                      fcn(t, y, k1);
                      tw = t+ (1.0/5.0)*h;
                      for(i = 0; i < N; i++){
                         w[i] = y[i] + h*(1.0/5.0)*k1[i];
                       }
                     fcn(tw, w, k2);
                      tw = t+ (3.0/10.0)*h;
                      for(i = 0; i < N; i++){
                       w[i] = y[i] + h*((3.0/40.0)*k1[i]+
                              (9.0/40.0)*k2[i]);
                            }

                         fcn(tw, w, k3);
                         tw = t+ (4.0/5.0)*h;
                         for(i = 0; i < N; i++){
                           w[i] = y[i] + h*((44.0/45.0)*k1[i] -n
                                 (56.0/15.0)*k2[i] +
                                 (32.0/9.0)*k3[i]);
                              }
                    fcn(tw, w, k4);
                    tw = t+ (8.0/9.0)*h;
                     for(i = 0; i < N; i++){
                       w[i] = y[i] + h*((19372.0/6561.0)*k1[i] -
                             (25360.0/2187.0)*k2[i] +
                                (64448.0/6561.0)*k3[i] -
                                   (212.0/729.0)*k4[i]);
                         }
                   fcn(tw, w, k5);
                   tw = t + h;
                   for(i = 0; i < N; i++){
                         w[i] = y[i] + h*((9017.0/3168.0)*k1[i] -
                                (355.0/33.0)*k2[i] +
                                 (46732.0/5247.0)*k3[i] +
                                  (49.0/176.0)*k4[i] -
                                   (5103.0/18656.0)*k5[i]) ;
                            }
                  fcn(tw, w, k6);
                  tw = t + h;
                  for(i = 0; i < N; i++){
                     w[i] = y[i] + h*((35.0/384.0)*k1[i] +
                     (500.0/1113.0)*k3[i] + (125.0/192.0)*k4[i] -
                      (2187.0/6784.0)*k5[i] + (11.0/84.0)*k6[i]);
                    }
                  fcn(tw, w, k7);
                  errabs = 0;

                 for(i = 0; i < N; i++){
                   dy[i] = h*((35.0/384.0)*k1[i] +
                         (500.0/1113.0)*k3[i] +
                         (125.0/192.0)*k4[i] -
                          (2187.0/6784.0)*k5[i] +
                           (11.0/84.0)*k6[i]);

                    err[i] =  h*((71.0/57600.0)*k1[i]  -
                            (71.0/16695.0)*k3[i]
                            + (71.0/1920.0)*k4[i] -
                            (17253.0/339200.0)*k5[i] +
                        (22.0/525.0)*k6[i] - (1.0/40.0)*k7[i]);
                   errabs+=err[i]*err[i];
                }

               errabs = sqrt(errabs);
               if( errabs < eps){
                    t+= h;
                    for(i = 0; i < N; i++){
                         y[i]+=dy[i];
                    }
                 }
              chi=errabs/eps;
              chi = pow(chi, (1.0/6.0));
              if(chi > 10)    chi = 10;
              if(chi < 0.1)   chi = 0.1;
              h*=  0.95/chi;
              if( t + h > b ) h = b - t;
              iteration++;

             fprintf(fpw ,"%.25lf %.25lf %.25lf %.25lf \n",
                      t,y[0],y[1],y[2]);
             if(iteration > 30000) break;
            }

             fclose(fpw);
             return 0;
           }
</code>```

I also tried with Scipy . I am sharing my Scipy code below

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

      def odes(x,t):
      #assign each ode to a vector
      A=x[0]
      B=x[1]
      C=x[2]

      #define each ODE
      dAdt=B
      dBdt=-(1/6)*A*C
      dCdt=-(3*C*B)/A

      return [dAdt,dBdt,dCdt]

      #define initialcondition
      x0=[1,5,13]

      #declare a time vector
      t=np.linspace(0,30,1000)
      x=odeint(odes,x0,t)

      A=x[:,0]
      B=x[:,1]
      C=x[:,2]

      #plot the results
      plt.semilogy(t,A)
      plt.semilogy(t,B)
      plt.semilogy(t,C)
      plt.show()
</code>```


问:
Typo in a-priori error estimate in a Discontinuous Galerkin paper
I'm looking at this famous paper which is available in the link below:

Franco Brezzi, LD Marini, Endre Süli, Discontinuous Galerkin methods for first-order hyperbolic problems,  Mathematical Models and Methods in Applied Sciences Vol. 14, No. 12 (2004) pp. 1893-1903, doi:10.1142/S0218202504003866, ResearchGate.

I am wonrdering if I've found a typo or (more probably) I am missing something.

I will use the enumeration of the paper.
Let me quickly write the setting: they introduce a norm in eq. $(31)$ $$|||\delta|||^2 = ||\delta||_{0,\Omega}^2 + \sum_{e \in E_j} ||c_e^{\frac{1}{2}} [\delta]||_{0,e}^2$$ What is relevant is not what $c_e$, but just the definition of the norm.
In order to derive an a-priori bound for $\delta$ in this norm, they start from $$C_s |||\delta||| \leq a_h(n ,\delta)+ b_h^s(n,\delta)$$

and they bound both $a(n,\delta)$ and $b_h^s(n,\delta)$ in the following way:
$$a_h(n,\delta) \leq C h^{k+1}||\delta||_{0,\Omega}||u||_{k+1,\Omega} \quad (45)$$

$$b_h^s(n,\delta) \leq C h^{k+\frac{1}{2}} ||u||_{k+1,\Omega} \Bigl( \sum_{e \in E_h}
 ||c_e^{\frac{1}{2}}[\delta]||_{0,e}^2\Bigr)^{\frac{1}{2}} \quad (49)$$

So they bound the two terms in the sum $(43)$ and obtain $$C_s |||\delta||| \leq C h^{k+\frac{1}{2}} ||u||_{k+1,\Omega} |||\delta||| \quad (50)$$ My problem is exactly in this last bound: I don't see how it's possible to have only $h^{k+\frac{1}{2}}$ in the latter bound, as summing $(45)$ and $(49)$ I have a term with $h^k$ and a term with $h^{k+\frac{1}{2}}$. Also, the constant $C$ should not depend on $h$, since they wrote:

$C$ will denote a generic positive constant which depends only on the degree $k$ of the polynomials, on the minimum angle of the mesh, and on the maximum value of the stabilizing functions $c_e$.

Any help is highly appreciated.

EDIT:

Maybe I spot the tricky part (to me). In computing $(45)$ they're using the inverse inequality $$|\delta|_{1,T} \leq h_T^{-1} ||\delta||_{0,\Omega}$$ which I know it holds under the assumption $h = \max_T h_T \leq 1$

But if $h \leq 1$, then for sure I can bound $h^{k + 1}$ with $h^{k+\frac{1}{2}}$ and I would be able to prove this. Does it sound reasonable to you?


问:
Extracting a mid-plane for thick shell analysis
I have a complex part that contains features of the form shown in the figure below.

Because of the cost of 3D finite element simulation of the part, I want to try an analysis with 2D thick shells.

However, I am running into a problem when trying to extract the mid-surface of the part.  Automated tools provide me with a mid-surface that is shown by the red dashed line in the figure.  That leads to interpenetration of the top and bottom surfaces in bending.

I would like a mid-surface that is shown by the green dashed line.

Can anyone point me to automated mid-surface extraction algorithms that can do what I seek?


问:
Identifying an unknown P.D.E. from solution data
I have a black-box simulation that produces the time evolution of a probability density function p(x, t) in 1 dimension from arbitrary initial conditions p(x, 0). The underlying simulation occurs on a discrete lattice of sites x=0, 1, 2, ... but when p(x, t) is a sufficiently smooth function of x, we expect that the behavior is approximately captured by an unknown P.D.E.
(Additionally, even with non-smooth initial conditions solutions become smoother and smoother over time.)

The system is local, translationally-invariant and time-independent. However, the P.D.E. is almost certainly non-linear.
What is the best strategy to identify the unknown P.D.E.? The simulation can be run on many different initial conditions at low cost.

A few sub-questions:

is there a way to determine even just the order of a PDE (i.e. the maximum number of derivatives needed to describe the solutions, in space and in time) from solution data?

are there ways to produce useful initial conditions to run through the simulation?


问:
What is the most common loss function used with collocation methods for differential equations
I was looking at the Cheney and Kincaid book (6th edition) on numerical methods, with respect to collocation method for differential equations. Now for linear systems of ODES, collocation is just a linear solve. But for nonlinear systems of equations there is a residual function that we minimize.

The residual function they show in the Cheney and Kincaid book is p. 619 and shows.

$$
\mathcal{L} = \sum_{j=0}^{J}c_jv_j - b
$$

This seems a bit odd, since I was expect something like the usual squared loss instead of just the sum of the differences. With this formulation, the positive errors can cancel out the negative errors.

Could someone just confirm that this is indeed the correct loss or residual function to use, or let me know if there a different residual function that is more appropriate. I checked a few other books. Most of them talk about using collocation on a linear system, which like I said is just a linear solve.


问:
elementwise function of low rank approximation
I recently ran into an interesting result.

I have a matrix $D$ containing pairwise distances between all points in a dataset. This matrix is converted to a similarity matrix $S$ via an RBF kernel.

Because calculating the full distance matrix is expensive, I use randomized linear algebra techniques to approximate $S$. As $S$ is positive definite, there are some specialised techniques to find good approximations cheaply. For my application, the approximation performs very well.

In my research, I also tried approximating $D$ and converting the approximation to an approximate $S$ afterwards. I expected that this wouldn't perform well, because I thought a low rank approximation of $D$ would not approximate small elements well and these small elements become large elements in $S$. However, the experiments indicate that it performs very well.

Can anyone help me explain these results? Am I wrong in thinking that the low rank approximation will not approximate small elements well?


问:
Sample Average Approximation vs. Numerical Integration
To calculate the expected value of objective functions, we have two choices:

Sample Average Approximation (SAA):
$$
\frac{1}{N}\sum_{i=1}^N f(x,\xi^i).
$$

Numerical Integration (e.g., Monte Carlo method):
$$
\frac{1}{N} \sum_{\xi^i\in\Xi^N}^N f(x,\xi^i)p(\xi^i),
$$
where $p(\cdot)$ is a probability density function and $\Xi^N\subset\Xi$ is the approximation of the support $\Xi$ of the probability function.

The value of SAA converges to the true expectation ($\mathbb{E}[f(x,\xi)]$) as $N\to\infty$, as does the numerical integration, due to the law of large numbers. However, from the numerical viewpoint, what is the difference between these two methods?

More specifically, in stochastic optimization or expected residual minimization of the least squared mean for a system of equations, how do these methods differ in convergence rate?