问:
How can this attempt to implement a Milstein numerical approximation for an Ito process of multiple components be fixed?
I have been reading Kloeden and Platen's Numerical Solution of Stochastic Differential Equations, and have more or less been trying to systematically complete the various exercises therein as I go along. I have been effectively stuck on a single problem for a long while, and despite returning to attack it periodically, I have not been successful. At this point, naturally, I am pretty desperate for any insight that might lead me to finally moving past it. To that end, I shall try my best to explain my problem succinctly.

I am given the following Ito process:
$(1)$ $dX_t = -\frac{1}{2}X_tdt + X_tdW_t^1 + X_tdW_t^2$,

which I have determined has the exact solution given by $(2)$  $X_t = X_0e^{((-\frac{1}{2} - 1^2)t + 1\cdot(W_t^1 + W_t^2))}$,

where here I write the drift coefficient $ a = -\frac{1}{2}$ distinct from the diffusion coefficient $b = 1.0$ for clarity later. Also, $X_0 = 1$. I will add that this text makes the (I think) somewhat regrettable convention choice of writing $a(x,t)$ shorthand as simply $a$. So a regular sight is to see something like "$a(x,t) = ax$" where the latter $a$ is of course just the constant $-\frac{1}{2}$ I mentioned. Similarly, "$b(x,t) = bx$" in this problem. (Also, this all corresponds to the SDE introduced in question 10.2.1 of Kloeden.)

Now, my problem is a later one, 10.3.3, in which I am to implement the following generalized Milstein numerical scheme for a $k$-dimensional Ito process with $m$ components. In the case here, k = 1 and m = 2:

$(3)$ $Y_{n+1}^k = Y_n^k + \underline{a}^k\Delta + \sum_{j=1}^m b^{k,j}\Delta W^j + \sum_{{j_1,j_2}=1}^m\underline{L}^{j_1}b^{k,j_2}J_{(j_1,j_2)}$,

with $\Delta$ = the step size = $h$ in my code, $\Delta W^j$ = the $n$-th Gaussian increment of the Wiener process $W^j$, and so in the implementation, $\Delta W^j$ = $\sqrt{h}\cdot W_n$, with $W_n$ the $n$-th standard normal random number, which in my case is generated using Matlab's $randn$ function. $\underline{a} = a - \frac{1}{2}bb'$, $\underline{L}^{j_1} = \sum_{{i}=1}^kb^{i,j_1}\frac{\partial}{\partial{x^i}}$, the usual Stratonovich operator for obtaining the coefficient functions.

Finally, the Stratonovich integrals themselves, $J_{(j_1,j_2)}$, are approximated in the following way:

$J^p_{(j_1,j_2)} = \Delta\Big(\frac{1}{2}\xi_{j_1}\xi_{j_2} + \sqrt{\rho_p}(\mu_{j_1,p}\xi_{j_2} - \mu_{j_2,p}\xi_{j_1})\Big) + \frac{\Delta}{2\pi}\sum_{{r}=1}^p\Big[\frac{1}{r}\big(\zeta_{j_1,r}\big(\sqrt{2}\xi_{j_2} + \eta_{j_2,r}\big) - \zeta_{j_2,r}\big(\sqrt{2}\xi_{j_1} + \eta_{j_1,r}\big)\big)\Big]$,

where $\rho_p = \frac{1}{12} - \frac{1}{2\pi^2}\sum_{{r}=1}^p\frac{1}{r^2}$.

$p \in \mathbb{N} $, with larger $p$ producing a better approximation, $\xi_j = \frac{1}{\sqrt{\Delta}}\Delta W^j = W_n^j$. Again, I think the notation here is somewhat ambiguous, but basically it seems to me it would be clearer to write  $\xi_n^j = W_n^j$ to make it absolutely clear that it is the $n$-th value of $\xi_j$. In any case, $\xi_j,\mu_{j,p},\zeta_{j,r},$ and $\eta_{j,r}$ are independent $N(0;1)$ Gaussian random variables.

Now, I honestly don't know at this point if I have been too explicit to the point of pedantic in defining things, or not enough-- I really don't think I am personally at the level of competence in this subject yet that I have any sense at all of what "should be taken for granted". But if I've left anything out or haven't been clear enough, I am quite happy to come back and remedy that as needed.

So after all that, my problem is that my implementation doesn't seem to work and after weeks of writing and re-writing it, I can't seem to make it work. And sadly, I have literally no other colleagues available I can talk to about this. So I am desperately hoping someone here has the knowledge and willingness to look over what I have done, and hopefully spot and point out whatever mistake I have made. As it stands, this ought to be a strong order 1.0 method, but my implementation isn't displaying any error convergence behavior with decreasing step size at all. And honestly, it's driving me mad. I do not know what I am missing. Anyway, without further ado, here is what I have (for use in Matlab):

```<code>clear;

t0 = 0;
T = 1;
X0 = 1.0;
a = -0.5;
b = 1.0;

p = 2;

trajectories = 1000;
error = [];

for N = [16,32,64,128];
trajsum = 0;
for i1 = 1:trajectories;

    h = (T-t0)/N;

    W1 = randn(1,N+1);
    W1(1) = 0;
    W2 = randn(1,N+1);
    W2(1) = 0;

    delW1 = zeros(1,N);
    delW1(1,1:N) = W1(1,2:N+1)*sqrt(h);
    delW2 = zeros(1,N);
    delW2(1,1:N) = W2(1,2:N+1)*sqrt(h);

    sumdelW1 = zeros(1,N);
    sumdelW2 = zeros(1,N);

    for i = 1:N;
        tempsum1 = 0;
        tempsum2 = 0;
        for j = 1:i;
            tempsum1 = tempsum1 + delW1(1,j);
            tempsum2 = tempsum2 + delW2(1,j);
        end
        sumdelW1(1,i) = tempsum1;
        sumdelW2(1,i) = tempsum2;
    end

    milE1 = delW1/sqrt(h);
    milE2 = delW2/sqrt(h);

    milstein = zeros(1,N+1);
    milstein(1,1) = X0;

    exact2 = zeros(1,N+1);
    exact2(1,1) = X0;

    for i = 1:N;
        milstein(1,i+1) = milstein(1,i) + (a - b^2/2)*milstein(1,i)*h + b*milstein(1,i)*(delW1(1,i) + delW2(1,i)) + (1/2)*b^2*milstein(1,i)*((delW1(1,i))^2 + (delW2(1,i))^2) + b^2*milstein(1,i)*(J12(h,p,milE1(1,i),milE2(1,i)) + J21(h,p,milE1(1,i),milE2(1,i)));
        exact2(1,i+1) = X0*exp((a - (b^2))*h*i + b*(sumdelW1(1,i) + sumdelW2(1,i)));
    end
    trajsum = trajsum + abs(milstein(1,N+1) - exact2(1,N+1));
end
trajsum = (1/trajectories)*trajsum;
error = [error trajsum];
end

error

figure(1);
subplot(1,1,1);
plot([0:h:1],exact2(1,1:N+1),'k-');
hold on;
plot([0:h:1],milstein(1,1:N+1),'r--');
</code>```

Moreover, I use the following two functions to simulate the stochastic integrals.

J12 given by:

```<code>function a = J12(h,p,E1,E2);
N1 = randn(1,p);
N2 = randn(1,p);
C1 = randn(1,p);
C2 = randn(1,p);
U1 = randn;
U2 = randn;

pp = 0;

for i = 1:p;
    pp = pp + 1/i^2;
end

pp = (1/12) - (1/(2*(pi^2)))*pp; % This completes getting rho-P.

a = h*((0.5)*E1*E2 + sqrt(pp)*(U1*E2 - U2*E1));

y = 0;
for k = 1:p;
    y = y + (1/k)*(C1(1,k)*(sqrt(2)*E2 + N2(1,k)) - C2(1,k)*(sqrt(2)*E1 + N1(1,k)));
end

a = a + (h/(2*pi))*y;

end
</code>```

And J21 given by:

```<code>function a = J21(h,p,E1,E2);
N1 = randn(1,p);
N2 = randn(1,p);
C1 = randn(1,p);
C2 = randn(1,p);
U1 = randn;
U2 = randn;

pp = 0;

for i = 1:p;
    pp = pp + 1/i^2;
end

pp = (1/12) - (1/(2*(pi^2)))*pp; % This completes getting rho-P.

a = h*((0.5)*E2*E1 + sqrt(pp)*(U2*E1 - U1*E2));

y = 0;
for k = 1:p;
    y = y + (1/k)*(C2(1,k)*(sqrt(2)*E1 + N1(1,k)) - C1(1,k)*(sqrt(2)*E2 + N2(1,k)));
end

a = a + (h/(2*pi))*y;

end
</code>```

Anyway my apologies for the essay. Of course I am happy to clean it up if it doesn't conform to standards. As to the problems, regardless of the choice of p, or N (the number of steps) that I divide the interval into, checking the "error" variable makes it clear that the error is not showing the order of convergence this method is supposed to have, nor indeed any convergence at all. I am sincerely grateful to anyone willing to give this a look. Thank you.


问:
About the the stability of using an explicit scheme on the heat equation
Before I get to the heat equation I'd like to talk about the advection equation.
Descritize that with FD in time and BD in space:
\begin{equation}
\dfrac{u^{n+1}_i - u^{n}_i}{\Delta t} + v \dfrac{u^{n}_{i} - u^{n}_{i-1}}{\Delta x}= 0
\end{equation}
there are higher order terms which were cancelled out so let's bring those back in ie. trucnation errors:
\begin{equation}
\dfrac{\partial u}{\partial t} + \dfrac{\Delta t}{2}\dfrac{\partial^2 u}{\partial t^2}  + v \dfrac{\partial u}{\partial x} -v\dfrac{\Delta x}{2}\dfrac{\partial^2 u}{\partial x^2} = 0
\end{equation}
Move things around a lil:
\begin{equation}
\dfrac{\partial u}{\partial t}   + v \dfrac{\partial u}{\partial x}  = -v\dfrac{\Delta x}{2}\dfrac{\partial^2 u}{\partial x^2} + \dfrac{\Delta t}{2}\dfrac{\partial^2 u}{\partial t^2}
\end{equation}
with some vodo magic
\begin{equation}
\dfrac{\partial u}{\partial t}   + v \dfrac{\partial u}{\partial x}  = -v\dfrac{\Delta x}{2}\dfrac{\partial^2 u}{\partial x^2} + v^2\dfrac{\Delta t}{2}\dfrac{\partial^2 u}{\partial^2 x}
\end{equation}
and some compression
\begin{equation}
\dfrac{\partial u}{\partial t}   + v \dfrac{\partial u}{\partial x}  = \dfrac{v}{2} \left(- \Delta x + v \Delta t    \right)\dfrac{\partial^2 u}{\partial^2 x} + \text{higher order terms}
\end{equation}
When the convection equation is discretized we don't solve the original equation but the equation written above. If the RHS coefficient  $-\Delta x + v \Delta t$

Positive $-\Delta x + v \Delta t > 0 \implies \dfrac{v \Delta x}{\Delta t}<1$ we get the CFL condition and a positive diffusion coefficient which is why explicit schemes are explicit. A low CFL would imply more diffusion.

negative $-\Delta x + v \Delta t > 0$ implies a negative coefficient for second order derivative which is why the solution blows up.

Now I am trying to apply this rational to the heat equation. Descritize that with FD in time and CD in space:
\begin{equation}
\dfrac{u^{n+1}_i - u^{n}_i}{\Delta t} = D  \dfrac{u^{n}_{i+1} -2  u^{n}_{i} + u^{n}_{i-1}}{\Delta x^2}
\end{equation}
the truncated equation:
\begin{equation}
\dfrac{\partial u}{\partial t} + \dfrac{\Delta t}{2}\dfrac{\partial^2 u}{\partial t^2}  = D \dfrac{\partial^2 u}{\partial x^2} +  D \dfrac{\Delta x^2}{12}\dfrac{\partial^4 u}{\partial x^2}
\end{equation}
and finally
\begin{equation}
\dfrac{\partial u}{\partial t}  = D \dfrac{\partial^2 u}{\partial x^2} -D^2 \dfrac{\Delta t}{2}\dfrac{\partial^4 u}{\partial x^4} +  D \dfrac{\Delta x^2}{12}\dfrac{\partial^4 u}{\partial x^4}
 \end{equation}
compress a lil:
\begin{equation}
\dfrac{\partial u}{\partial t}  = D \dfrac{\partial^2 u}{\partial x^2}
+\left (-D^2 \dfrac{\Delta t}{2} +  D \dfrac{\Delta x^2}{12} \right)
   \dfrac{\partial^4 u}{\partial x^4}
 \end{equation}
I wouldn't want a negative coefficient next to the fourth order derivative because that would mean that I would be solving an equation that blows up (I guess) so I want to guarantee that:

$-D^2 \dfrac{\Delta t}{2} +  D \dfrac{\Delta x^2}{12} > 0 \implies \dfrac{D \Delta t }{\Delta x^2}< \dfrac{1}{6} $

which is not consistent with the Von Neumann stability condition for this scheme which is:

$\dfrac{D \Delta t }{\Delta x^2}< \dfrac{1}{2} $

any hints where I went wrong?

ContentLicense: CC BY-SA 4.0


问:
Global reconstruction defined elementwise in a-posteriori error estimator
This question is a follow-up of this previous one. In "Error Control for Discontinuous Galerkin Methods for First Order Hyperbolic Problems" by Georgoulis et al., an error estimator is provided in equation (4.7).

It contains the term $$||f -b \cdot \nabla \hat u - cu_h - \Pi(f-b \cdot \nabla \hat u-c u_h) ||_T$$

where :

$\Pi$ is the local orthogonal $L^2$ projection onto the classical polynomial DG finite element space $S_h= \{v \in L^2: v_{|T} \in \mathbb{P}^{p}(T) \}$

$u_h$ is the finite element solution, which I have on each cell $T$

$\hat{u}$, called reconstruction, is defined elementwise (see definition 4.3 or the beginning of Section 4) by the relation $$\int_Tb \cdot \nabla \hat{u} v_h = \int_T b \cdot \nabla u_h v_h - \int_{\partial_{-T}} (b \cdot n)[u_h]v_h^{+}$$
for all $v_h \in S_h$ (classical DG fe space)

It's written that the reconstruction space is $\{w_h \in C({\bar{\Omega}}): w_{h|T} \in \mathbb{P}^{p+2}(T)\}$

My big issue is with this reconstruction part!

As far as I understand, I should construct $\hat{u}$ cell-by-cell by using that definition. Of course, as I know $u_h$ and $v_h$, it turns out that a linear system for each cell has to be solved.

However, it seems to me that $\hat{u}$ is in a continuous finite element space, instead of a discontinuous one like $v_h$. So, by expanding the l.h.s I would have

$$(b \cdot \nabla \varphi_j,\phi_i)_T$$

where $\varphi_j$ is a basis function in the continuous finite element space, while $\phi_i$ is a basis function of the discontinuous one $S_h$.

I am a bit confused on how to do this. If I create a different DoFHandler object on the same triangulation, distribute the dofs corresponding to the continuous space, I would not even end up with a square matrix since I don't have the same number of degrees of freedom in these two spaces. So I think I am clearly missing the point. Any help is extremely appreciated.


问:
Comparing minimas of two different functions
The goal is to find vectors $x_u$ and $y_i$, both of the same length $f=64$, and to do this the following loss function is minimized:

$$\sum_{u, i} (1 + \alpha \cdot r_{ui})(p_{ui} - x_{u}^{T}y_i)^2$$

The vectors obtained are user-factors and item-factors, used for the purpose of rating calculation in a Recommender System. My goal is to find the optimum value of $\alpha$, and the initial thought was to try different values for $\alpha$, for which I would calculate the min of the above function; in the end, I would pick $\alpha$ for which the function has the lowest value. Is this a good idea? How could one find an optimal $\alpha$ value?

As for the $r_{ui}$ and $p_{ui}$ values, these would be the same for each user-item pair in each run of the minimization. Given that the function is used for a Recommender System, one way to identify good $\alpha$ would be to withhold a test set, and then check how the recommendations (top $x_u^{T}y_i$ values) are recovered.


问:
Weird "oscillatory" modes appearing in FEM simulations
I am using COMSOL to solve a mathematical model involving thermoelectric hydrodynamic (TEMHD) flow. I am running a very large parameter sweep and using the solutions obtained to make some plots. However, I get a rather strange phenomenon appearing from time to time, in which the solution appears to exhibit strange oscillatory behaviour, and I have tried several things such as mesh refinement, reducing the relative tolerance, and changing the type of solver.

Below is one of the plots of the velocity inside an open channel as part of a parameter sweep.

However, the solution with neighbouring parameter values much more closely resembles this:

There is no mathematical reason as far as I am aware as to why this should happen, and there seems to be no regularity as to when these quirky plots arise. Even then, this only seems to happen a few percent of the time but is still rather frustrating.

Does the phenomenon I have illustrated have a name? How can it be resolved?


问:
Attempt on 2d Advection with FDM - With Code
I tried to implement the 2d advection problem with a velocity field, that is not constant in space. My problem is, that the "mass" of my shifted density gets "eroded" or just disappears, i.e. the norm $\|u_t\|$ goes to zero as time progresses. At least I want to minimize its decay somehow.

The initial condition has non zero pixels in each corner where it is basically a plateau disc function.

The velocity field points always to the center of the domain, so that the points have to meet in the center.

Here is my code.

```<code>    # -*- coding: utf-8 -*-

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation

n = 32
X = np.linspace(-2,2,n)
Y = np.linspace(-2,2,n)

XX,YY = np.meshgrid(X,Y)

#------------- velocity field --------
# a disc in each corner
u0 = np.zeros((n,n))
u0 += (XX-1.25)**2 + (YY-1.25)**2 < 0.25
u0 += (XX+1.25)**2 + (YY+1.25)**2 < 0.25
u0 += (XX-1.25)**2 + (YY+1.25)**2 < 0.25
u0 += (XX+1.25)**2 + (YY-1.25)**2 < 0.25

# flow field points always to origin
flow = np.zeros((2,n,n))
flow[0,:,:] = YY
flow[1,:,:] = XX
flow *= -1 #

# normalizing flow field - is it necessary?
norm = np.linalg.norm(flow, axis = 0, ord=2)
flow[0,:,:] =  np.divide(flow[0,:,:], norm, out=np.zeros_like(norm), where=norm>10**-7)
flow[1,:,:] =  np.divide(flow[1,:,:], norm, out=np.zeros_like(norm), where=norm>10**-7)

dx = (XX.max() - XX.min())/n
dt = dx/5

u = np.copy(u0)
U = [u0.copy()]

frames = 100

for i in range(frames):

    #upwind scheme
    ux_minus = np.vstack((np.zeros((1,n)), u[1:-1,:] - u[0:-2,:],np.zeros((1,n))))
    uy_minus = np.hstack((np.zeros((n,1)), u[:,1:-1] - u[:,0:-2],np.zeros((n,1))))

    # downwind scheme
    ux_plus = np.vstack((np.zeros((1,n)), u[2::,:] - u[1:-1,:], np.zeros((1,n))))
    uy_plus = np.hstack((np.zeros((n,1)), u[:,2::] - u[:,1:-1], np.zeros((n,1))))

    pulse_x = (flow[0,:,:] < 0) * flow[0,:,:]*ux_plus + (flow[0,:,:]>0)*flow[0,:,:]*ux_minus
    pulse_y = (flow[1,:,:] < 0) * flow[1,:,:]*uy_plus + (flow[1,:,:]>0)*flow[1,:,:]*uy_minus

    u -= (dt/dx)*(pulse_x + pulse_y)

    U.append(u.copy())

def update(i):
    matrice.set_array(U[i])

fig, ax = plt.subplots()
plt.show()
matrice = ax.matshow(U[0])

plt.colorbar(matrice)

ani = animation.FuncAnimation(fig, update, frames=frames, interval=1)

</code>```


问:
Variational loss of hp-Variational Physics Informed Neural Networks for 2D-Poisson Equation in Tensorflow
I am trying to reproduce the results from the hp-VPINN paper (https://arxiv.org/pdf/2003.05385.pdf) on tensorflow (v1) for Poisson's equation, particularly the two-dimensional Poisson equation.

In one dimension the variational loss is given by
$$
    L^p = \sum_{e=1}^{N_{el}}\frac{1}{K^{(e)}}\sum_{k=1}^{K^{(e)}}\left|\mathcal{R}_k^{(e)}\right|^2
$$

where $N_{el}$ is the number of elements and $K^{(e)}$ is the number of test functions used (in this case corresponding to different degrees of Legendre polynomials)

where $\mathcal{R}_k^{(e)}$ is
$$
    \mathcal{R}_k^{(e)} =  \int_{x_{e-1}}^{x_e} \frac{du_{NN}(x)}{dx} \frac{dv_k^{(e)}(x)}{dx}dx  - \left. \frac{du_{NN}(x)}{dx}v_k^{(e)}(x)\right|_{x_{e-1}}^{x_{e}} - F_k^{(e)}
$$

with $u_{NN}$ being the network's output and and $v_k^{(e)}$ the test function. The boundary term vanishes since the test functions are such that they are zero on the element boundaries
and $F_k^{(e)} = \int_{x_{e-1}}^{x_e} f(x)v_k^{(e)}dx$.

To calculate it Gauss-Lobatto quadrature is used giving the following code

```<code>import tensorflow as tf
import numpy as np
u_nn_element = tf.reshape(
    tf.stack(
        [
            tf.reduce_sum(
                w_quadrature
                * du_nn_quad_element
                * dtest_quad_element_dx[i]
            )
            for i in range(n_test_functions)
        ]
    ),
    (-1, 1),
)
f_element = jacobian * np.asarray(
    [sum(w_quadrature * f(x_quadrature) * t(x_quadrature) for t in test_functions]
)
residual_element = u_nn_element - f_element
loss_element = tf.reduce_mean(tf.square(residual_nn_element))
total_varloss += loss_element
</code>```

When combined with the boundary loss (mean squared error for the boundary points at $[-1, 1]$)

In two dimensions the variational loss is given by

$$
    L= \sum_{e_x=1}^{N_{el_x}}\sum_{e_y=1}^{N_{el_y}}\frac{1}{K_1K_2}\sum_{k=1}^{K_1K_2}\left|\mathcal{R}_k^{(e_xe_y)}\right|^2
$$

with a residual given by

$$
    {}^{(2)} \mathcal{R}_{k_1k_2}^{(e_xe_y)} = - \int_{x_{e-1}}^{x_e}\int_{y_{e-1}}^{y_e} \left(\frac{\partial u_{NN}(x)}{\partial x} \frac{d\phi_{k_1}^{(e_x)}(x)}{dx}\phi_{k_2}^{(e_y)}(y) + \frac{\partial u_{NN}}{\partial y}\phi_{k_1}^{(e_x)}(x)\frac{d\phi_{k_2}^{(e_y)}(y)}{dy} \right) dxdy \\
- F_k^{(e)}
$$

with the test functions $\phi_{k_1} \phi_{k2}$

My code attempt has been

```<code>l = [
    tf.reduce_sum(
        -wy_quadrature
        * wx_quadrature
        * (
            du_nn_dx * dtest_x_quad_element[i] * y_test_quad_element[j]
            + du_nn_dy
            * dtest_y_quad_element[j]
            * x_test_quad_element[i]
        )
    )
    for i in range(n_test_functions_x)
    for j in range(n_test_functions_y)
]
u_nn_element = tf.reshape(
    tf.stack(l),
    (-1, 1),
)

# element boundaries
jacobian_x = (x1 - x0) / 2
jacobian_y = (y1 - y0) / 2

f = [
    sum(xw_quad * yw_quad * f_quad_element * t1 * t2)
    for t1 in self.test_functions
    for t2 in self.test_functions
]
self.f = jacobian_x * jacobian_y * np.asarray(f)
self.f = self.f.flatten()
self.f = self.f[:, None]

residual_nn_element = u_nn_element - f_element
loss_element = tf.reduce_mean(tf.square(residual_nn_element))
self.varloss_total += loss_element
</code>```

though my variational error doesn't converge and the net only learns the boundary values. How is the double integral supposed to look in code?


问:
Calculating the number of Flops of SPH density calculation
I would like to calculate the number of floating point operations (Flops) my code is performing in my machine. To do so, I would like to be sure I am counting the operations in the inner-most loop correctly.

My understanding is that the relevant number, which are the operations inside the vectorized loop, total $25 N^2$, $13 = 3\times 1 + 4 \times 2$ outside the b-spline kernel and $12 = 6+6$ operations inside the kernel.

```<code>int compute_density_3d_ref(int N,double h,
                           double* restrict x, double* restrict y,
                           double* restrict z, double* restrict nu,
                           double* restrict rho){
  const double inv_h = 1./h;
  const double kernel_constant = w_bspline_3d_constant(h);
  #pragma omp parallel for
  for(size_t ii=0;ii<N;ii+=1){
    double xii = x[ii];
    double yii = y[ii];
    double zii = z[ii];
    double rhoii = 0.0;
    #pragma omp simd reduction(+:rhoii) aligned(x,y,z,nu)
    for(size_t jj=0;jj<N;jj+=1){
      double q = 0.;

      double xij = xii-x[jj]; // One operation
      double yij = yii-y[jj]; // One operation
      double zij = zii-z[jj]; // One operation

      q += xij*xij; // Two operations
      q += yij*yij; // Two operations
      q += zij*zij; // Two operations

      q = sqrt(q)*inv_h; // Two operations

      // Two operation plus the number of operations in w_bspline
      rhoii += nu[jj]*w_bspline_3d_simd(q);//*w_bspline_3d(q); // box->w(sqrt(dist),h);
    }
    rho[ii] = kernel_constant*rhoii;
  }

  return 0;
}
</code>```

The functions called are:

```<code>double w_bspline_3d_constant(double h){
  return 3./(2.*M_PI*h*h*h);
}

#pragma omp declare simd
double w_bspline_3d(double q){
  double wq = 0.0;
  double wq1 = (0.6666666666666666 - q*q + 0.5*q*q*q); // Six operations
  double wq2 = 0.16666666666666666*(2.-q)*(2.-q)*(2.-q); // Six operations

  if(q<2.)
    wq = wq2;

  if(q<1.)
    wq = wq1;

  return wq;
}
</code>```

I am counting ${\rm sqrt(q)} = \sqrt{q}$ as a single operation under the suposition that the CPU can dispatch one (vector) sqrt instruction per cycle, suposing AVX256. I am also counting additions and multiplications both as floating point operations on equal footing.

I would like to know if I going about this calculation right or I am over-estimating (or even under-estimating) the number of operations involved in ths computation.

Thanks in advance.


问:
Weird behavior in for solving TISE in harmonic oscillator potential using the shooting method
I was solving the time independent Schrödinger equation using the shooting method for harmonic oscillator potential. This is the code that I wrote for that with the results (code is written in julia):

m, omega, hbar, del_x = 1.0, 1.0, 1.0, 1e-3

function get_eigenfunction(xgrid, V_x, E_n, del_x, m=1; hbar=1, s=1e5)
    psi_vals = [0, s]  # Initial Conditions

    for i = 2:(length(xgrid) - 1)
        psi_ii = 2psi_vals[i] * ((m * del_x^2) * (V_x[i] - E_n) / hbar^2 + 1) - psi_vals[i-1]
        append!(psi_vals, psi_ii)
    end

    normalize!(psi_vals)
    return psi_vals
end

E_100 = (100 + 0.5)*hbar*omega
x1 = collect(-15:del_x:15)
V_x1 = [1/2 * m * omega * x^2 for x in x1]
phi2_x1 = get_eigenfunction(x1, V_x1, E_100, del_x).^2
scaling = maximum(V_x1) / maximum(phi2_x1)
plot(x1, V_x1)
plot!(x1, scaling*phi2_x1 .+ E_100)
plot!(x1, [E_100 for x in x1])

```
<p><img src="https://i.stack.imgur.com/qU3kr.png" alt="Plot Result 1" /></p>
<p>Now, when I increase the range of x, to sample more value of potential then the result is quite weird as I get an insane overshoot at the right boundary point:</p>
<pre class="lang-matlab prettyprint-override"><code>
E_100 = (100 + 0.5)*hbar*omega
x1 = collect(-20:del_x:20)
V_x1 = [1/2 * m * omega * x^2 for x in x1]
phi2_x1 = get_eigenfunction(x1, V_x1, E_100, del_x).^2
scaling = maximum(V_x1) / maximum(phi2_x1)
plot(x1, V_x1)
plot!(x1, scaling*phi2_x1 .+ E_100)
plot!(x1, [E_100 for x in x1])

</code>```

Is there any way to prevent this problem?


问:
2d advection-diffusion: cell Péclet number and numerical stability
I am studying the numerical resolution of 2d advection-diffusion problems with finite element methods.

$$\frac{\partial u}{\partial t} + \beta\cdot\nabla u = \nabla\cdot(\nabla u) \, .$$

It is said in this article on the COMSOL blog that it has been mathematically proven that the condition on the cell Péclet number

$$\mathrm{Pe}:= \frac{\Vert \beta\Vert h}{2c} > 1\, ,$$

leads to numerical instabilities.

Could you please point me to any reference rigorously justifying such a statement?


问:
What is the limit involving `Sum`, `Subsets`, and `RankedMax` as `t` approaches infinity?

Motivation

Suppose we have a countably infinite $A$ with order and group structures and suppose $F_1,F_2,\cdot\cdot\cdot$ are an infinite sequence of finite sets (denoted $\left\{F_n\right\}_{n=1}^{\infty}$) such that $\bigcup\limits_{n=1}^{\infty}F_n=A$, and $F_1\subset F_2\subset \cdot\cdot\cdot$

We state that $F_t$ is equidistributed in the interval $[\inf(A),\sup(A)]=[a,b]$ when for any sub-interval $[c,d]$ of $[a,b]$ we have there exists $\frac{d-c}{b-a}$ such for every arbitrarily small positive $\epsilon$ there exists a sufficiently large integer $N$ where for all $t\ge N$

$$\left|\frac{\left|F_t\cap[c,d]\right|}{|F_t|}-\frac{d-c}{b-a}\right|\le\epsilon$$

The discrepancy measures how much a structure deviates from being equidistributed in $[a,b]$. While there already have been measures for the discrepancy such as

$$\sup_{a\le c \le d \le b}\left|\frac{\left|F_t\cap[c,d]\right|}{|F_t|}-\frac{d-c}{b-a}\right|$$

The discrepancy only measures the maximum deviance of $F_t$ from equidistribution in $[a,b]$ and always converges to zero as $t$ grows larger.

Hence, we need an alternative definition to the discrepancy that converges to different values, depending on the arrangement of $F_t$

Core Question

I wish to compute my measure of deviance from equidistribution (for particular $F_t$). In my new definition, a set is equidistributed when the differences between consecutive pairs in a set are all the same.

For instance $\left\{\frac{1}{3},\frac{2}{3},\frac{3}{3}\right\}$ is equidistriubted since $\frac{2}{3}-\frac{1}{3}=\frac{3}{3}-\frac{2}{3}$
and $\left\{\frac{1}{4},\frac{2}{4},\frac{3}{4},\frac{4}{4}\right\}$ is also equidistriubted since $\frac{2}{4}-\frac{1}{4}=\frac{3}{4}-\frac{2}{4}=\frac{4}{4}-\frac{3}{4}$

We then measure the deviance from equidistribution using the following notations:

Definition:

If:

$\Delta{F_t}$ is a multiset that arranges the elements in $F_t$ from least to greatest and takes the difference between consecutive pairs. (For example if $F_3=\left\{\frac{1}{3},\frac{2}{3},\frac{3}{3}\right\}$ then $\Delta F_3=\left\{\frac{2}{3}-\frac{1}{3},\frac{3}{3}-\frac{2}{3}\right\}=\left\{1/3,1/3\right\}$)

$||F_t||=\sum\limits_{x\in \Delta F_t}x$ normalizes the elements in $\Delta F_t$ so the sum of all the elements is $1$ (similar to a probability distribution.)

$\mathcal{P}_{n}(\Delta F_t/||F_t||)$ takes all subsets of $\Delta F_t/||F_t||$ (with a cardinality of $n$) such that all elements (including repeated elements) are treated distinctively. (For example if $\Delta F_3/||F_3||=\left\{\frac{1}{3},\frac{1}{3},\frac{1}{3}\right\}$ then $\mathcal{P}_{3}\left(\Delta F_3/||F_3||\right)=\left\{\left\{\frac{1}{3},\frac{1}{3},\frac{1}{3}\right\}\right\}$ and $\mathcal{P}_{2}\left(\Delta F_3/||F_3||\right)=\left\{\left\{\frac{1}{3},\frac{1}{3}\right\},\left\{\frac{1}{3},\frac{1}{3}\right\},\left\{\frac{1}{3},\frac{1}{3}\right\}\right\}$)

$S_{n,t}=\bigcup\limits_{X\in P_{n}(\Delta F_t/||F_t||)\setminus\emptyset}\left\{\prod\limits_{x\in X}x\right\}$

${\max}^{r}(S_{n,t})$ represents the $r$-th largest element in $S_{n,t}$. For example if $\Delta F_3/||F_3||=\left\{1/6,1/2,1/2\right\}$ then ${\max}^{1}\left(\Delta F_3/||F_3||\right)=\frac{1}{2}$, ${\max}^{2}\left(\Delta F_3/||F_3||\right)=\frac{1}{2}$, and ${\max}^{3}\left(\Delta F_3/||F_3||\right)=\frac{1}{6}$

And $$N_{n,t}=\frac{2|\Delta F_t/||F_t|||^{n-1}}{\binom{\left|\Delta F_t/||F_t||\right|}{n}+1}$$

Starting with $n=1$, we arrange the values in $S_{1,t}$ taking its greatest value, its greatest value plus the second greatest, its greatest plus second plus third greatest and continuing till we add the greatest up to the least greatest value. We then take the mean of these values giving $d_{1,t}$ and repeating this process for $n=2$ up to $n=|F_t|$. This causes the mean to converge to different values depending on the arrangement of $F_t$.

$$d_{n,t}=\frac{1}{|S_{n,t}|}\sum\limits_{p=1}^{|S_{n,t}|}\sum\limits_{r=1}^{p}{\max}^{r}\left(S_{n,t}\right)$$

If $F_t$ was completely equidistributed in $[a,b]$ we would get different values for $d_{n,t}$ for different $n$. Since want every $d_{n,t}$ to be equal, we multiply them by $$N_{n,t}=\frac{2|\; \Delta F_t/||F_t||\;|^{n-1}}{\binom{\left|\;\Delta F_t/||F_t||\;\right|}{n}+1}$$

And add them up:

$$\sum_{n=1}^{|F_t|}N_{n,t}d_{n,t}$$

Hence for every $F_t$ we want a $d(F_t,A)$ such for every arbitrarily small positive $\epsilon$ there exists a sufficiently large integer $N$ where for all $t\ge N$

$$\left|d(F_t,A)-\sum_{n=1}^{|F_t|}{N_{n,t}}d_{n,t}\right|\le\epsilon$$

or if we expand:

$$\left|d(F_t,A)-\sum_{n=1}^{|F_t|}\frac{N_{n,t}}{|S_{n,t}|}\sum\limits_{p=1}^{|S_{n,t}|}\sum\limits_{r=1}^{p}{\max}^{r}\left(S_{n,t}\right)\right|\le\epsilon$$

For a given $A$, the $F_t$ that gives $|d(F_t,A)-1|$ the smallest value is the most equidistrubuted $F_t$

Examples to compute

If $A=\left\{\frac{1}{n}:n\in\mathbb{N}\right\}$ and $F_t=\left\{\frac{1}{n}:n\in\mathbb{N},n\le t\right\}$, I did the following to compute $d$:

```<code>F[t_] := F[t] = 1/Range[1, t]
F1[t_] := F1[t] = Total[Differences[F[t]]]
F2[t_] := F2[t] = Sort[Differences[F[t]]/F1[t]]
d1[t_, n_] := d1[t, n] = Subsets[F2[t], {n}];
S[q_, t_] :=
 S[q, t] = (2 (Length[F2[t]])^(q - 1))/(Binomial[Length[F2[t]], q] + 1)
d[t_] := d[t] =
  Sum[S[q, t] Mean[
     Table[Sum[
       RankedMax[
        Table[Times @@ d1[t, q][[s]], {s, 1, Length[d1[t, q]]}],
        n], {n, 1, p}], {p, 1, Length[d1[t, q]]}]], {q, 1,
    Length[F[t]]}]
</code>```

But once t goes beyond $10$ it takes more time to compute, (if it converges it would do so extremely slow).

I tried getting someone to compute this on Mathematica but no one was able to answer. The moderator sent the question to computational stack exchange but someone stated that answering through programming is impossible.

Is wish for a way we can mathematically simplify this for any $F_t$ but as a user pointed out, generalizing this for any $F_t$ is impossible. However, I wish to compute $d$, where $A=\left\{\frac{1}{n}:n\in\mathbb{N}\right\}$, for the following $F_t$:*

$F_t=\left\{\frac{1}{n}:n\in\mathbb{N},n\le t\right\}$

$F_t=\left\{{1}/{\left[2^t/m\right]}:m\in\mathbb{N}, m\le 2^t \right\}$ (where $[\cdot]$ is the nearest integer function)

Here, I hope 1. and 2. will give different results for $d$.

I also wish to calculate $d$, where $A=\mathbb{Q}\cap[0,1]$, for the following:

$F_t=\left\{\frac{m}{n!}: m,n\in\mathbb{N}, m\le n!\le t\right\}$

$F_t=\left\{\frac{m}{n}:m,n\in\mathbb{N},m\le n\le t\right\}$

I also hope, in this case, that 1. and 2. give different results with $d$.

Edit: Perhaps combinatorics identities such as in this answer helps.

For any set $S$ (with distinct elements):

$$\sum_{X \in \mathcal{P}(S)} (\prod_{x\in X} x) = \prod_{x\in S} (1+x)$$

However, this time we are applying to $\mathcal{P}_{n}(S)$ instead of $\mathcal{P}(S)$

OwnerDisplayName: Arbuja


问:
Scheduling Algorithm Goal is to Fix Max Number of Appointments
First post.  Hi everyone.

I’m trying to develop an algorithm to later code that schedules appointments in a way such that the number of hours occupied by the appointments in a given day is maximized.

Appointments can be an hour or half hour long and all appointments should have a 15 minute buffer between them.

I’m getting stumped by the complexity of trying to come up with rules for this.

It’s currently being done by humans.

Someone texts, says they want a half hour or full hour, specifies a time that works best for them and they are fit into the schedule like a puzzle by a human being.

Preference is given to fitting the most appointments per day into a given 8 hour day.  Hour appointments are preferred but half hour is just fine.

The texters have times that work for them and a daily schedule is built up by hand in a very qualitative way.

I’m trying to automate this and it’s hard. Really hard.

There are cases where the wrong time appointment could be booked and it would eat up 2 hours of time.  For example:

Say you have 1pm to 3pm open without appointments.  Then say you have 4pm to 6pm open. A pair of 2 hour time slots. Ideally, you would want to fit a full hour and a half hour into each of those time slots.  However, if someone books their own half hour appointment at 1:45, , through the new self booking system, it makes it such that you can’t book any other appointments except the poorly placed 1:45 half hour in the entire 2 hour open period!  (Remember there is a 15 min buffer between each appointment). Now imagine someone does the same thing again booking 4:45 also.  You now have only 2 half hours taking up 4 hours of the day.  Bad!

The human mind just fits these all in as best it can and easily optimizes for the maximum amount of booked appointment time in a day. I’m really struggling to develop algorithms to do the same.  Any ideas?

Thank you!


问:
Stable iterative solver for complex symmetric linear systems
I am interested in the iterative solution (preferably Krylov-type solvers) of a problem $\boldsymbol{A}x=b$, with $x,b\in\mathbb{C}^{n\times1}$   and $\boldsymbol{A}\in\mathbb{C}^{n\times n}$. $\boldsymbol{A}$ is symmetric, invertible, and its real and imaginary parts are indefinite. This comes from solving a harmonic, damped Finite Element problem.

Now, because the idea is to write a few iterations explicitly by hand, the challenge is that I would like to use a:

simple algorithm (i.e. without the need for orthogonalization during each step, as in GMRES) which converges well with a diagonal preconditioner.

in which the estimates of $x_n$ are guaranteed to "monotonically decrease" the residual $r_n = b - \boldsymbol{A}x_n$. By this I mean that, for example, $||r_n||_2$ or $r_n^H\boldsymbol{A}r_n$ are guaranteed to decrease monotonically per iteration. Sort of like the Conjugate Gradient method for hermitian matrices.

Are there any such methods?

I tried applying the Conjugate A-Orthogonal Conjugate Residual (COCR) method, by Sogabe and and Zhang, 2007[1], with simple diagonal preconditioning. I have a working implementation for coarse meshes, but for less coarse meshes, the algorithm fails to converge.

I then read about the Quasi-Minimal Residual COCR (QMRCOCR) method, by Gu et. al., 2014[2], which seems to be much more stable, still relatively simple to implement and the error seems to decrease monotonically. But this last remark is my interpretation of the results, not any rigorous convergence analysis. And a Symmetric Successive Over-Relaxation (SSOR) preconditioner is always used in the article, so it cannot be known whether the convergence will be as well-behaved if a simple diagonal preconditioner were used instead.

A possible solution that comes to mind is solving $$\boldsymbol{A}^H\boldsymbol{A}x=\boldsymbol{A}^Hb$$ where the superscript H denotes the hermitian adjoint. Since $\boldsymbol{A}$ is invertible, then $\boldsymbol{A}^H\boldsymbol{A}$ is necessarily hermitian (it is also positive-definite) Therefore, the simple Conjugate Gradient method could be used.

Now, everywhere I read about this mentions stability issues. I cannot find the reason why. The only concrete thing I have found is that the condition number of $\boldsymbol{A}^H\boldsymbol{A}$ will be equal to the condition number of $\boldsymbol{A}$ squared. But in my FE code, I always apply left- and right-diagonal preconditioning of A, thereby obtaining relatively good condition numbers.

So why exactly is solving $\boldsymbol{A}^H\boldsymbol{A}x=\boldsymbol{A}^Hb$ with an iterative solver bad?

After some experiments, the answer to the above question is: it just is.

Here are some images of the the curves $||r_n||/||b||$ for solving my problem using the Conjugate Gradient method (CGM) of $\boldsymbol{A}^H\boldsymbol{A}$, of the COCR method of the complex linear system and of the QMRCOCR method. I only used simple Jacobi preconditioning. It is a simple example for a small $\boldsymbol{A}$ of shape $6930 \times 6930$.

Figure 1: Clearly, CGM of $\boldsymbol{A}^H\boldsymbol{A}$ is bad in this case.

Figure 2: COCR and QMRCOCR both converge in 1132 and 1136 iterations, respectively. Clearly QMRCOCR is more stable.

Figure 3: The issue for me is that the residuals of neither COCR or QMRCOCR are guaranteed to decrease in the first iterations. For CGM it does decrease.

So is there a way to guarantee at least 2 first steps of a decreasing residual for COCR or QMRCOCR? I thought the first step of these types of methods was always similar to a gradient descent. Maybe for complex matrices, things work differently?

References:

Sogabe, T.; Zhang, S. A COCR method for solving complex symmetric linear systems. J. Comput. Appl. Math. 2007, 199 (2), 297–303. DOI: 10.1016/j.cam.2005.07.032.

Gu, X.; Huang, T.; Li, L.; Li, H.; Sogabe, T.; Clemens, M. Quasi-Minimal Residual Variants of the COCG and COCR Methods for Complex Symmetric Linear Systems in Electromagnetic Simulations. IEEE Trans. Microwave Theory Techn. 2014, 62 (12), 2859–2867. DOI: 10.1109/TMTT.2014.2365472.


问:
Help with restart functionnality in sef-made GMRES solver in python
I am new to this forum and to computational science in general.

I started to learn numerical liner algebra on my own and would like to code a GMRES solver in python (no preconditioner for the time being) with the ability to restart. My end goal is to code a Flexible GMRES [2].
However, I am facing some difficulties in understanding the restart part of GMRES(m).

From [1], the algorithm for the GMRES(m) is:

Start: Choose $x_0$ and compute $r_0 = f - Ax_0$ and $v_1 = \frac{r_0}{\left \| r_0 \right \|} $

Iterate: for $j= 1,2,3,...m$ do: [note: $m$ should be the number of restart iterations]

$h_{i,j} = (Av_j, v_i), i=1,2,...,j$,

$\hat{v}_{j+1}=Av_j-\sum_{i=1}^j h_{i,j}v_i$,

$h_{j+1,j} =\left \| \hat{v}_{j+1} \right \| $,

$v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1,j}}$

Form the approximate solution:

$x_m = x_0 + V_m y_m$, where $y_m$ minimizes $\left \| \beta e_1 - \bar{H}_m y \right \|$, $y \in R^m$

Restart:

Compute $r_m = f-A x_m$;

if satisfied then stop, else compute $x_0:=x_m$ and $v_1:=\frac{r_m}{\left \| r_m \right \|}$ and go to point 2.

Note: Point 2 is known as Arnoldi iterations.

I have the following code for a GMRES solver without restart.

```<code>def GMRES (A, b, max_iterations, x0=None, tol = 1e-8):
    """
    Parameters
    ----------
    A : np.array
        (m x m) numpy array. Coefficient matrix
    b : np.array
        (m x 1) numpy array. Source term vector
    max_iterations : int
        maxium number of iterations to perform
    x0 : np.array
        (m x 1) numpy array. Initial guesss for the iterative procedure
    tol : double
        tolerance. The default is 1e-8.

    Returns
    -------
    x : np.array
        (m x 1) numpy array.
    residualVector : list
        values of residual throughout the iterative procedure

    """
    # Algorithm status variable
    status = -1

    # Size of solution
    m = len(b)

    if np.all(x0==None):
        x0 = np.zeros((m, 1))

    # Compute the residual
    r0 = b - (A @ x0)

    # Compute beta
    beta = norm(r0)

    # List to store the residual values
    residualVector = [beta]

    for i in np.arange(1, (max_iterations + 1)):

        # Compute Arnoldi iterations
        Q, h = arnoldi_iter(A, r0, i)

        # Solve least square problem
        y = np.linalg.lstsq(h, beta*e1(i), rcond=None)[0]

        # Computes the error of the iteration
        rn = error (r0, h, y, i)

        # Appends the value of error to the list
        residualVector.append(rn)

        # Check convergence
        if (rn < tol):
            # Solver terminated successfully
            status = 0
            break

    if (i == max_iterations):
        # Solver did not converge
        status = max_iterations

    # Form solution
    x = x0 + np.matmul(Q[:,:-1], y )

    return x, residualVector, status
</code>```

The code for the Arnoldi Iterations is:

```<code>def arnoldi_iter(A, b, N):
    """
    Computes a basis of the (n+1)-Krylov Subspace of A

    Krylov subspace A = <b, bA, A²b, A³b, ... A^nb>

    Parameters:
        A: m x m Matrix (numpy array)
        b: m x 1 vector (numpy array)
        N: Number of iterations

    returns:
        Q: The columns of an orthogonal basis of A's Krylov subspace. m x (n+1) Matrix
        h: A on basis Q, It is an upper-Hessenber Matrix of size (n+1) x n
    """

    # Set m size of A
    m = len(b)

    # Initializes the Matrices to the required dimensions
    Q = np.zeros((m, N+1))

    h = np.zeros((N+1, N))

    # Computes the first q-vector
    q = b/norm(b)

    # Sets the first column of Q to the transpose of the q-vector
    Q[:, [0]] = q

    # Loop
    for n in np.arange(N):

        # Calculate v-vector (New candidate vector)
        v = A @ q

        for j in np.arange(n + 1):
            qj= Q[:, j].reshape(m, 1) # make vector (m, 1) 2-D array instead of (m,). Causes problems with subtraction when using (m,1) - (m,) vectors

            h[j, n] = ( qj.conj() ).T @ v

            v = v - ( h[j, n] * qj)

        h[n + 1, n] = norm(v)

        tol = 1e-15 # avoid division by 0

        if (h[n+1, n] > tol):
           q = v/h[n + 1, n]

           Q[:, [n + 1]] = q

        else:
            return Q,h
    return Q,h
</code>```

The function for computing the error is:

```<code>def error (b, h, y, i):
    """
    Parameters
    ----------
    b : np.array
        residual vector.
    h : np.array
        Hessenberg matrix at the i'th iteration
    y : np.array
        Solution value at i'th iteration.
    i : int
        iteration number.

    Returns
    -------
    err :
        DESCRIPTION.

    """
    # Compute norm of the b vector
    beta = norm(b)

    # Compute the error

    err = norm(np.matmul(h, y) - ( beta*e1(i) ))

    return err
</code>```

And finally, the function to return the first column vector of the identity matrix of size n+1 is:

```<code>def e1(n):
    """
    Compute the first column vector of the identity matrix with size (n+1) x 1

    Parameters
    ----------
    n : iteration counter

    Returns
    -------
    tmp : np.array
        first column vector of the identity matrix with size (n+1) x 1

    """
    tmp = np.zeros((n + 1 , 1))
    tmp[0] = 1
    return tmp
</code>```

I have been trying out the above algorithm as:

```<code>N = 30

b = np.random.normal(1, scale=5, size=(N, 1))

# Create non-singular matrix. Singular matrices are not inverteble ...
matrixQualityCheck = False

while(not(matrixQualityCheck)):
    A = np.random.normal(1,scale=5, size=(N, N))
    determinant = np.linalg.det(A)

    if(determinant != 0):
        matrixQualityCheck=True

# Solve system exactly
x = np.linalg.solve(A,b)

# Solve with GMRES
x2, RESVEC, status = GMRES(A, b, N, tol=1e-12)

# Compute the norm of the difference between exact and GMRES solution
normDiff = np.abs(norm(x - x2))

if (  normDiff < 1e-10  ):
    print("Same result")
else:
    print("solvers are different")
    print(normDiff)
</code>```

The above code appears to work correctly.

Now, for the restart part of the code I have:

```<code>def restart_GMRES(A, b, max_iterations, restart_iterations, x0=None, tol = 1e-8):

    for i in np.arange(restart_iterations):
        if (i == 0):
            # First iteration computes normal GMRES with x0 by the user
            x, RESVEC, status = GMRES(A, b, max_iterations, x0=x0, tol=tol)
        else:
            x, RESVEC, status = GMRES(A, b, max_iterations, x0=x, tol=tol)

        if (status == 0):
            # solver finished with success
            return x, RESVEC, [status, i]
    return x, RESVEC, [status, i]
</code>```

Since I have a status variable inside the code that is telling me that the solution convergence to the correct value I am using it to exist the restart algorithm. However, the results do not appear to be Ok and I do not know why...

I would really appreciate some help!

Sorry for the long post.

References

[1] - Saad, Y., & Schultz, M. H. (1986). GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7(3), 856–869. https://doi.org/10.1137/0907058

[2] - SAAD, Y. (1993). A FLEXIBLE INNER-OUTER PRECONDITIONED GMRES ALGORITHM. SIAM Journal on Scientific Computing, 14(2), 461–469.


问:
Multigrid method: linear solver and modified residual
I am trying to better understand the FAS multigrid algorithm for Euler equation in FV discretization. The usage of the modified residual (the residual with forcing) inside the different cases:

Unsteady computations : referring to the Jameson's paper I understand that it is used in the dual time step computation that allows for a faster convergence. This fact seems confirmed by this paper, where a similar approach is used.

Steady computations: if I have well understood (seen that there is not the time integration) is used in the linear solver process, hence for example with a Krylov based method.

I wanted to ask if these two points are correct, and, in the second case, if we have to rely compulsorily on a Krylov based method (e.g GMRES). From my understanding another solver (e.g Jacobi) would not take into account the modified residual with forcing and so would not take into account the history of the different multigrid levels residual, making the multigrid algorithm ineffective.


问:
Varying Young modulus in FEM simulation
I'm working on a project for which I have inherited some FEM code. This implemented FEM calculates, given some force field, displacements on a discretised square using square elements and assumes a homogeneous isotropic and linear material.

I want to do away with the assumption homogeneity. I would like to assume a spacial dependent Young's modulus $E(x) = r(x) E_0$, where $E_0$ is some fixed quantity and $r(x)$ is a scaling quantity which can be different for every element.

I use the $E_0$ as the Youngs modulus in creating the element stiffness matrices using the material matrix $D$
$$
D= \frac{E_0}{1-v^2} \left(\begin{matrix} 1 & v & 0 \\ v & 1 & 0 \\ 0 & 0 & \frac {1}{2}(1-v) \end{matrix}\right).
$$
In the assembly of the stiffness matrix I do the following. Suppose we are adding the element matrix $K_e$ corresponding to element $e$, I then add $r(e) K_e$ to the global stiffness matrix. Here $r(e)$ is just $r(x)$ evaluated on some point of the element $e$.

Is this a valid approach for modelling non-homogeneous materials? If not, can someone point me to a source that does explain how to model such materials in detail?

If relevant: I have a mathematics background, but am new in using FEM in structural problems.

Title: Varying Young modulus in FEM simulation


问:
Looking for a library for solving convex-convex quadratic fractional programming problems
Is there a library of some programming language (MatLab, C, Python, R, etc.) that includes a function for efficient calculation of convex-convex quadratic fractional programming problems of the form

$$
\begin{align*}
\sup& &\frac{x^TQx}{x^TPx}\\
\text{s.t.}& &A^Tx \leq b\\
&&x\geq 0
\end{align*}
$$
where $Q\in\mathbb{R}^{n\times n}$ is a given symmetric positive semidefinite matrix, $P\in\mathbb{R}^{n\times n}$ is a given symmetric positive definite matrix, $A \in \mathbb{R}^{m\times n}$ is a given matrix, $b \in\mathbb{R}^{m\times 1}$ is a given column vector, and $m,n \in \{1,2,...\}$ are given.

The function should run fast when $m,n \leq 500$.

EDIT: I'll be satisfied with an implementation of the algorithm described in the following paper:

R. Yamamoto and H. Konno: An efficient algorithm for solving convex-convex quadratic fractional programs. J. Optim. Theory Appl. 133 (2007), 241–255. DOI:10.1007/s10957- 007-9188-y


问:
Floquet theory for periodic delay differential equations: current numerical routines
I would like to determine the stability of a system of periodic delay differential equations (a seasonal host-parasite model). I've tried to implement the method described in Lemma 2.5 in this paper:

Liang, X., Zhang, L. & Zhao, XQ. Basic Reproduction Ratios for Periodic Abstract Functional Differential Equations (with Application to a Spatial Model for Lyme Disease). J Dyn Diff Equat 31, 1247–1278 (2019).:

Lemma 2.5 Assume that (E, E+) is an ordered Banach space with E+ being normal and Int(E+) $\neq$ ∅,
which is equipped with the norm ||·||E . Let L be a
positive bounded linear operator. Choose v0 ∈
Int(E+) and define an =
||Lvn-1||E , vn = $\frac{Lv_{n-1}}{
> a_{n}}$, ∀n ≥ 1. If $\lim_{(n→+∞)}a_n$ exists, then r(L) =
$\lim_{n→+∞}a_n$.

But I'm not sure it's working properly (I've probably made a mistake somewhere in the implementation). The problem is it doesn't always seem to correspond to equilibrium points I find when simply simulating the model for a very long time. I looked into different ways of approximating Floquet multipliers/monodromy matrices etc. for delay equations, just so I could double check my work, and there seem to be many different approaches (e.g., semi-discretization, finite element analysis, Chebyshev polynomials, etc.)

My question, for anyone much more familiar with this literature than I am, is: Is there a currently preferred method? Are there any that are easy and quick to implement (I'm not a very sophisticated programmer or mathematician)? I would be especially interested in an example with published code. Some advice would be greatly appreciated.


问:
Can you provide an example of an algorithm that has strong scalability but not weak scalability? Can you also provide an example the other way around?

Strong Scaling and Amdahl's Law

Strong scaling is a measure of how, for a fixed overall problem size, the time to solution decreases as more processors are added to a system. An application that exhibits linear strong scaling has a speedup equal to the number of processors used.

Strong scaling is usually equated with Amdahl's Law, which specifies the maximum speedup that can be expected by parallelizing portions of a serial program. Essentially, it states that the maximum speedup S of a program is:

$$S=\frac{1}{(1−P)+P/N}$$

Here P is the fraction of the total serial execution time taken by the portion of code that can be parallelized and N is the number of processors over which the parallel portion of the code runs.

In reality, most applications do not exhibit perfectly linear strong scaling, even if they do exhibit some degree of strong scaling. For most purposes, the key point is that the larger the parallelizable portion P is, the greater the potential speedup. Conversely, if P is a small number (meaning that the application is not substantially parallelizable), increasing the number of processors N does little to improve performance. Therefore, to get the largest speedup for a fixed problem size, it is worthwhile to spend effort on increasing P, maximizing the amount of code that can be parallelized.

Weak Scaling and Gustafson's Law

Weak scaling is a measure of how the time to solution changes as more processors are added to a system with a fixed problem size per processor; i.e., where the overall problem size increases as the number of processors is increased.

Weak scaling is often equated with Gustafson's Law, which states that in practice, the problem size scales with the number of processors. Because of this, the maximum speedup S of a program is:

$$S=(1-P)+NP$$

Here P is the fraction of the total serial execution time taken by the portion of code that can be parallelized and N is the number of processors over which the parallel portion of the code runs.

the scaled speedup is calculated based on the amount of work done for a scaled problem size (in contrast to Amdahl’s law which focuses on fixed problem size). Or say, it is not the problem size that remains constant as we scale up the system but rather the execution time.

Can you provide an example of an algorithm that has strong scalability but not weak scalability? Can you also provide an example the other way around?


问:
How to compute the Eigenvalue and Eigenstates of Quantum well with Effective mass using finite difference method in Python?
I want to compute the eigenvalues and eigenstates of a quantum well with different effective masses of electron in the barrier and in the quantum well. As can be seen [1]: https://github.com/mholtrop/QMPython/blob/master/Finite%20Well%20Bound%20States.ipynb [1], the mass of the electron is constant, therefore the Hamiltonian can be easily diagonalized and eigenvalues/eigenstates can be evaluated.

However, in my situation, the mass of the electron takes different values in the barrier and the well. Could any one suggest how I could obtain the matrix which satisfies the effective mass condition, and hence can be used to compute eigenvalues and states using similar method.

$$a_i\psi_{i-1}+b_i\psi_{i}+c_i\psi_{i+1}=E\psi_{i}$$

$$a_{i+1}=c_i=\frac{-\hbar^2}{2m^*(\delta z)^2} \text{ and } b_i=\frac{\hbar^2}{m^*(\delta z)^2}+V_i$$

In the above equations, it is apparent that the Hamiltonian can be diagonalized as mentioned in the link above.

$$a_{i+1}=c_i=\frac{-\hbar^2}{2m_{i+\frac{1}{2}}^*(\delta z)^2} \text{ and } b_i=\frac{\hbar^2}{2(\delta z)^2}\bigg(\frac{1}{m^*_{i+\frac{1}{2}}}+\frac{1}{m^*_{i-\frac{1}{2}}}\bigg)+V_i$$

My question is: how can I diagonalize my Hamiltonian with $m^*$(z) being different in the barrier and the well using the finite difference.

The complete Schrodinger equation  can be written as,

$$\frac{-\hbar}{2}\frac{\partial}{\partial z}\frac{1}{m^*(z)}\frac{\partial}{\partial z}\psi(z)+V(z)\psi(z)=E\psi(z)$$

$\qquad\qquad\qquad\qquad\qquad\quad$

Here, mass depends on z growth axis of quantum well (in picture its x instead of z) $m^*(z)$ in case of finite quantum well the value mass in barrier remains same in both regions 1 and 3 and is different in the well region 2.

I would like to know how can I use the same approach as mentioned [1] to solve  the
Schrodinger equation with a spatially dependent mass.

OwnerDisplayName: celerion


问:
Fast evaluation of trigonometric polynomials
Suppose you have a trigonometric polynomial of the form
\begin{equation*}
x(t) = \sum_{k = 0}^N a_k \cos(2 \pi k f_0 t).
\end{equation*}
Using Clenshaw algorithm, one can evaluate this polynomial in $O(N)$, just as an ordinary polynomial via Horner's rule. Now:

Is it possible to evaluate $x(t)$ in less than $O(N)$ in general? Also negative or partial results are appreciated.

What if $x(t)$ has a particular form? For example, if $a_k = \frac{(-1)^k}{k}$, $x(t)$ is the truncated Fourier series of a sawtooth wave.

In case of negative answer, is there at least some fast algorithm that approximates $x(t)$?

Thank you in advance.


问:
Is this a legit way to sample a random matrix spectrum?
In order to undergird a theoretical model concerning many body physics, I want to have exponentially large eigenvalue spectra from the random matrix GOE ensemble. its properties are mainly

(i) a semicircle distribution,

(ii) a level repulsion between the eigenvalues.

These are met in some physical systems, so it is natural to test a model using those random matrices.

However, my PC can diagonalize only 40000x40000 matrices before reaching RAM limitations so I tried to come up with a clever way of creating larger spectra by means of the corresponding PDF and CDF. I want to  make use of the inverse transform sampling method which transforms uniformly distributed numers into eg semicircular ones, which then matched the (i) condition, but obvioulsy had rather a Poissoinian level attraction, because the uniform input was uncorrelated in the first place. So my idea was to sample random numbers which inherit the desired level spacing from the (ii) property by means of the corresponding CDF and feed these numbers into the inverse sampling described above.
Here you can see the side by side comparison of the two method in the case of a 10000x10000 matrix, left being a diagonalized matrix spectrum, and right following from the algorithm described.
is this a proper method to generate arbitrarily large "spectra"?
I found that for larger N the right semicircular histogram becomes much smoother, maybe it converges only for asymptomatical N?

can I by chance incorporate another feature of random matrix theory, the Tracy-Widom distribution?

can I perform a suitable test to judge the quality of this method?

or am I doing a horribly naive mistake in my wish of finding a neat and clever RAM-workaround?


问:
Good non oscilliatory derivatives for an exsisting grid
I'm calculating the entropy production of a shockwave by utilizing the equations:
\begin{equation}
\sigma = J'_q\frac{\partial}{\partial x}\left(\frac{1}{T}\right) +\frac{1}{T}\frac{4\eta}{3}\left(\frac{\partial v}{\partial x}\right)^2
\end{equation}
and
\begin{equation}
\sigma = \frac{\partial \rho_s}{\partial t} +\frac{\partial J_s}{\partial x}
\end{equation}
I have access to all the variables on a spatial $[x_i,...x_N]$ and temporal grid $[t_i,...t_K]$. Currently I've used central discretization on both the spatial and temporal grid. The method works for the first equation, but for the second there seems to be oscillations before and behind the shock wave. Are there any better discretization methods I can use? Should I use a different method in the temporal dimension?


问:
How to identify the most different short sequence from a corpus that does not occur?
I would like to generate the most unique string (20 long, primer) of a defined length that does  not  occur in a much larger string (6.4 billion long, human genome). Of course, there could be many unique strings, but I want to avoid string with any mismatches/deletions/insertions that could occur in the larger stringer. My problem is easier than regular text because there only 4 possible elements (A, T, C, G) in the string, but I also need to worry about the reverse complement of the unique string that might occur in the larger string.

I have tried identifying these unique string by brute force search for mismatch, but I was wondering if there is better method. For bioinformatics experts, I did a BLAST search against the genome from a randomly generated string.


问:
What is the best method to do a MC Integration of a multidimensional integral where the integration limits depend upon other variables?
What is the best method to do a Monte Carlo Integration of a multidimensional integral where the integration limits depend upon other variables?

I am interested in getting a numerical value of a 5 dimensional integral, the limits look like this:

$$\int_{0}^{a^2} dS_{12} \int\limits_{S_{12}}^{a^2} dS_{123} \int\limits_0^{(a^2-S_{123})(S_{123}-S_{12})/S_{123}} dS_{34} \int\limits_0^1 d\tilde{S}_{13}\int_\limits{-\pi/2}^{\pi/2} d\tilde{S}_{134}$$

I guess the first two integrals could be done by choosing any random number between their respective limits, but what about the three last ones? I could find no source that could tell me about that, the only reference I could find that barely talks about this kind of stuff is this (page 41 in the pdf file): http://www.hep.fsu.edu/~harry/papers/MonteCarloTheoryPractice.pdf

My first attempt at solving this integral was with what the author of the previous link calls "the obvious way", which they say is wrong, and it seems like this is the case (I am calculating the decay width of a four-body decay and comparing it with other simulations and the value given by different collaborations).

Sorry if it seems to be a rather elementary question, this is the first time that I am doing a Monte Carlo integration problem.


问:
Numerical Partial Differentiation Check
In my computer vision course, we are working on extracting a 3D surface from a chain of 2D images taken under several conditions. This procedure is known as Photometric stereo. Prior to extracting the surface $f(x,y)$ numerically, we have to perform what's called a test of integrability :
$$
\frac{\partial^{2}f}{\partial x\partial y}=\frac{\partial^{2}f}{\partial y\partial x}
\tag{1}$$
This sanity test can be tested numerically by discretizing the partial derivative so that it accepts small values (and thus we have to accept that to recover $f$, we recover it starting from small values). Numerically this translates to :
$$
\frac{\partial\displaystyle\left(\frac{\alpha(x,y)}{\gamma(x,y)}\right)}{\partial y}-\frac{\partial\displaystyle\left(\frac{\beta(x,y)}{\gamma(x,y)}\right)}{\partial x}\approx0
\tag{2}$$
where :
$$
\frac{\partial f}{\partial x}=\frac{\alpha(x,y)}{\gamma(x,y)}\qquad\text{and}\qquad\frac{\partial f}{\partial y}=\frac{\beta(x,y)}{\gamma(x,y)}
\tag{3}$$
essentially where $(\alpha(x,y),\beta(x,y),\gamma(x,y))$ are the coordinates of the measured value of the unit-normal vector at some point on the surface $f$.

In short, how can we numerically check if $(2)$ is valid? what algorithm or procedure should we do to discretize the partial differentiation?

I am not familiar with advanced concepts such as discretizing partial differential equations. However, I am familiar with numerical integration techniques.

Title: Numerical Partial Differentiation Check


问:
Is there a way to generate a sample $(X_i, Y_i, Z_i)$ from custom distribution?
I'm newbie here.

I'm wondering if it's possible to generate $(X_i, Y_i, Z_i)$ from my own distribution function? I know that there is a way to make own class for 1D variable. But what about 3D case?


问:
Digital beamforming: how payload manipulation can change beam direction without manipulating the carrier?
I'm interested in how digital beamforming works and I can not find an answer for a lot of time.
I googled, asked teammates, and couldn't get it.

Let me describe my question.

From my understanding, the beam means the area in the space where waves from a number of antenna elements are in superposition which means, in turn, sum of amplitudes. Those waves must be in-phase for that and the "phase" here belongs to the carrier HF wave. Payload wave (envelope) is NOT changed here and it is obviously not needed for such beamforming.

If it's so I can imagine how analogue beamforming acts by an appropriate shifting of phase of HF carrier.

In digital beamforming, HF carrier phase is not manipulated and all antenna elements should radiate the beam in the straight direction in this case. The theory says, that by applying some vector to the payload (BaseBand signal) the beam can be steered.

This is what I can not understand. How manipulating of payload can change the direction of the beam without manipulating the HF carrier (so, no phase shifting of the carrier)?

As I said earlier the beam from my understanding is the area where carriers from different antenna elements are in-phase.


问:
The implicit form of a NURBS curve
I am trying to evaluate and analyse a NURBS curve to generate a mechanism. I understand that the general form of a NURBS curve is commonly written as a parametric equation in the form of $f_{par}(t)$.

Now, the criteria to generate the above mechanism from a curve are well-established. However, they are based on the implicit form of the curve $f_{imp}(x, y, z)$ as a (quite high order) polynomial. These criteria describe, for example, the relationships among the coefficients of $f_{imp}(x, y, z)$, in addition to some properties of the curve (degree, etc).

My thought is, if I could represent the input NURBS curve in its implicit form, I could conveniently apply the above criteria to the NURBS curve (e.g., in an application, return error if a designer specifies an invalid set of parameters), and hopefully redefine these criteria within the context of NURBS specification. Am I correct to say this?

I also understand there exist implicitization techniques, to represent parametric curves in their implicit equations. However, so as to not reinvent the wheel: knowing that NURBS is pretty much an industry standard and widely researched, I was just wondering if this has been done before in a paper or article?

What I have so found so far:

It is mentioned in Les Piegl's The NURBS Book that "Theoretically, a precise conversion of a NURBS to piecewise implicit form is possible, using techniques known as implicitization", but no details is given (implicitization not being their interest)

Busé's 2014 paper "Implicit matrix representations of rational Bézier curves and surfaces" [link] (but not NURBS)

Thank you!

PS: previously asked in MathOverflow, but I was suggested to ask here.

Title: The implicit form of a NURBS curve


问:
State change of input-output system
Edited

Given a computer model $F:\mathbb{R}^3 \to \mathbb{R}$, with inputs $x, w$ and $z$, and output $y=F(x,w,z)$, where for any input, we are able to evaluate the output, my goal is to tune the inputs to achieve a certain state change (change in output).

More precisely:

Goal: Given an input $x_0, w_0, z_0$ and an output $y_0$, I want to increase my output to a state $y_1$.

Question: How to change my inputs to achieve the new state $y_1$.

My approach:
\begin{align}
dy &= \frac{\partial F}{\partial x}(x_0, w_0, z_0)dx + \frac{\partial F}{\partial w}(x_0, w_0, z_0) dw + \frac{\partial F}{\partial z}(x_0, w_0, z_0)dz \\
y_1-y_0 &= \frac{\partial F}{\partial x}(x_0, w_0, z_0) (x_1 - x_0) + \frac{\partial F}{\partial w}(x_0, w_0, z_0) (w_1 - w_0)+ \frac{\partial F}{\partial z}(x_0, w_0, z_0) (z_1 - z_0)
\end{align}
Since $x_0, w_0, z_0, y_0, y_1$ are given, and assuming we are able to compute from the model $F$ the numerical approximation of the derivative of $F$ with respect to each input at the point $(x_0, w_0, z_0)$, the new values $x_1, w_1, z_1$ then should lie on the plane given by the second equation above.

I am also aware that I made a linear approximation that is only exact when $F$ is linear with respect to the inputs.

I am going in the right direction in formulation the problem ? If not, any suggestions?

Title: State change of input-output system


问:
Help with debugging block GMRES
I have written block version of GMRES by referring [1] and MATLAB implementation of gmres. I need to write it for complex matrices. My block implementation when run on single RHS is giving correct answer but when run with block RHS is converging but LHS is not matching. Can anyone please help me to debug or point in direction how can this be debugged? Are there any standard set of matrix and vectors where intermediate vectors like Krylov subspace basis or H matrix are know (or can be easily computed)?

[1] Iterative Methods for Sparse Linear Systems, Saad Yousef

Pasting here my blockGMRES and wrapper testcode I am using.

```<code>
function [x, error, iter, flag] = myblockgmres1( A, x, b, max_it, tol )

% input   A        complex nonsymmetric positive definite matrix
%         x        complex initial guess vector block
%         b        complex right hand side vector block
%         max_it   INTEGER maximum number of iterations
%         tol      complex error tolerance
%
% output  x        complex solution vector block
%         error    real error norm
%         iter     INTEGER number of iterations performed
%         flag     INTEGER: 0 = solution found to tolerance
%                           1 = no convergence given max_it

   iter = 0;
   flag = 0;
   numRHS = min(size(b));

   bnrm2 = norm( b );
   if  ( bnrm2 == 0.0 ), bnrm2 = 1.0; end

   r =  b-A*x ;
   error = norm( r ) / bnrm2;
   if ( error < tol ) return, end

   [n,n] = size(A);
   m = max_it;
   p = numRHS;
   V(1:n, 1:numRHS) = zeros(n,numRHS);
   H = [];
   QH = [];
   cs = [];
   sn = [];
   e1 = zeros(n, numRHS);
   for i=1:numRHS
    e1(i,i) = 1 + 0i;
   end

   r = b-A*x;
   %calculate QR of initial residue r0
   [r0,R] = qr(r, 0);
   V(:,1:numRHS) = r0(:,1:numRHS);
   g = e1*R(1:numRHS, 1:numRHS);

     for j = 1:max_it,
         Avj = A*V(:,(j-1)*p+1:j*p);
         for i = 1:j
             H((i-1)*p+1:i*p, (j-1)*p+1:j*p) = V(:,(i-1)*p+1:i*p)'*Avj;
             Avj = Avj - V(:,(i-1)*p+1:i*p)*H((i-1)*p+1:i*p, (j-1)*p+1:j*p);
         end
         [Q2, R2] = qr(Avj, 0);
         H(j*p+1:(j+1)*p, (j-1)*p+1:j*p) = R2(1:numRHS, 1:numRHS);
         V(:, j*p+1:(j+1)*p) = Q2(:, 1:numRHS);
         %QH is H which gets reduced through Given's
         QH = [QH; zeros(p, (j-1)*p)];
         QH = [QH zeros((j+1)*p, p)];
         QH(:, (j-1)*p+1:j*p) = H(:, (j-1)*p+1:j*p);

         %apply Givens rotation in all previous iterations
         %to new block of columns
         for k = (j-1)*p+1:j*p
             %columns-wise
             for l = 1:k-1
                %reverse order of rows in that block
                for i = p:-1:1
                        % apply Givens rotation
                        col = k;
                        row1 = l + i - 1;
                        row2 = l + i;
                        rr = QH(row1, col);
                        hh = QH(row2, col);
                        QH(row1, col)  = cs(l, i)*rr - sn(l, i)*hh;
                        QH(row2, col) = conj(sn(l, i))*rr + cs(l, i)*hh;
                end
             end
             %new column reduction
             for i = p:-1:1
                col = k;
                l = k;
                row1 = l + i - 1;
                row2 = l + i;
                rr = QH(row1, col);
                hh = QH(row2, col);
                rnorm = norm(rr);
                hnorm = norm(hh);
                rhnorm = sqrt(rnorm*rnorm + hnorm*hnorm);
                sn(l, i) =  -1*(hh*rr) / (rnorm*rhnorm);
                cs(l, i) =  real(rnorm/rhnorm);
                temp = conj(sn(l, i))*g(row1, :) + cs(l, i)*g(row2, :);
                g(row1, :) = g(row1, :)*cs(l, i) + sn(l, i)*g(row2, :);
                g(row2, :) = temp;
                QH(row1, col) = cs(l, i)*rr - sn(l, i)*hh;
                QH(row2, col) = 0.0;
             end
         end
         'iter=', j
         noexit = 0;
         %F-norm
         for i=1:numRHS
             error(i)  = norm(g(j*p+1:(j+1)*p, i));
         end
         norm(error)
         if(norm(error) > tol)
           noexit = 1;
         end
         %solve triangular systtem
         if ( noexit == 0 || j == max_it),
            y = zeros(size(b));
            y = QH(1:j*p,1:j*p) \ g(1:j*p,:);
            %y = y*norm(b);
            y1(1:j*p, :) = y(1:j*p, :);
            x = x + V(:,1:j*p)*y1;
            break;
         end
     end

   if ( error > tol ) flag = 1; end;
</code>```

Testcode to run blockgmres, here you can replace A, RHS and LHS with your test matrices and vectors. If required I can attach test matrix and vector I am using.

```<code>clear all;
A = dlmread('SystemMatrix_1.00e+006_real.txt') + 1i*dlmread('SystemMatrix_1.00e+006_imag.txt');
RHS1 = dlmread('System_1.00e+006_RHS_block1_real.txt') + 1i*dlmread('System_1.00e+006_RHS_block1_imag.txt') ;
RHS0 = dlmread('System_1.00e+006_RHS_block0_real.txt') + 1i*dlmread('System_1.00e+006_RHS_block0_imag.txt') ;

LHS0 = zeros(size(RHS0));
LHS1 = zeros(size(RHS0));
LHS0_ = dlmread('System_1.00e+006_LHS_block0_real.txt') + 1i*dlmread('System_1.00e+006_LHS_block0_imag.txt') ;
LHS1_ = dlmread('System_1.00e+006_LHS_block1_real.txt') + 1i*dlmread('System_1.00e+006_LHS_block1_imag.txt') ;

err = 0;
iter = 0;

%[LHS0,flag,relres,iter,resvec] = gmres(A, RHS0, length(RHS0), 1e-4);
%[LHS1,flag,relres,iter,resvec] = gmres(A, RHS1, length(RHS0), 1e-4);

%[LHS0, err, iter] = mygmres(A, LHS0, RHS0, length(RHS0), 1e-4);
%[LHS1, err, iter] = mygmres(A, LHS1, RHS1, length(RHS0), 1e-4);

%call block GMRES
LHS = [LHS0 LHS1];
RHS = [RHS0 RHS1];
%[LHS0, err, iter] = myblockgmres1(A, LHS(:,1), RHS(:,1), length(RHS0), 1e-4);
%[LHS1, err, iter] = myblockgmres1(A, LHS(:,2), RHS(:,2), length(RHS0), 1e-4);
[LHS, err, iter] = myblockgmres1(A, LHS, RHS, length(RHS0), 1e-4);
LHS0 = LHS(:,1);
LHS1 = LHS(:,2);

close all
plot(abs(LHS0_))
hold on
plot(abs(LHS0))
legend('expected LHS0','LHS0')
fprintf('relative error = %f\n', norm(LHS0-LHS0_)/norm(LHS0_))
fprintf('relative error = %f\n', norm(LHS1-LHS1_)/norm(LHS1_))
<span class="math-container">```
```

Title: Help with debugging block GMRES


问:
Regularisation of ill-conditioned matrix-vector problem
I have a linear* problem which arises from an integro-differential system, and writes:
$$
(\mathbf{I}+\lambda \mathbf{A})x = b
$$
where $\mathbf{A}$ is a real full matrix, size $n\times n$, but is not symmetric and has some zero eigenvalues. I know what these eigenvalues correspond to.

When $\lambda$ increases (but is still much smaller than 1) the condition number of $\mathbf{I}+\lambda \mathbf{A}$ (let's call it $\mathbf{K}$) deteriorates extremely rapidly (to $10^{10}$ and more). Even with a direct solver (luckily $n=100$ is enough for me) the results become meaningless.

I have no certainty, but I believe this is linked with the zero eigenvalues of $\mathbf{A}$.
Is there a way I could exploit my knowledge of the corresponding eigenvectors to regularize the numerical resolution?

(* note that this is actually the fixed-point linear problem that I am solving, the original problem has $\lambda=\lambda(x)$)

EDIT:
$\mathbf{A}$ itself has indeed a very large condition number (numerically evaluated by numpy with L2-norm), $10^{19}$, with $\sigma_{\min} = 1.8e-07$  and $\sigma_{\max} = 3.2e11$ ($\sigma_{\min}$ should analytically be 0 of course, since as said above, $\mathbf{A}$ has a 0 eigenvalue). Its largest eigenvalue has modulus about 600, here is the complete spectrum for $n=102$


问:
FEM : energy minimization VS PDE solving

Engineering FEM

When I studied engineering, I learned the traditional approach for finite elements for elasticity. The point was to solve the PDE $-div(\sigma)=f$ as:

Multiply your PDE with a test function, integrate, find the weak form: $\int \sigma(u): \epsilon(v) - \int v f = 0$

Interpolate $u=\sum_i u_i N_i$, set $v=\sum_j N_j$

Assemble the linear system and solve it to find $u_i$

Computer graphics FEM

Then I enrolled in a PhD in computer graphics, and the way they do FEM for simulation is totally different. The point is to minimize the integrated energy with Newton Raphson

Given the integrated energy $\int \psi$, compute the gradient $G = \nabla_u  \int \psi$ and the hessian $H = \nabla_{u} \nabla_{u} \int \psi$

Discretize the domain, for each element compute $H_{elt}$ and $G_{elt}$

Assemble the full linear system, and apply Newton Raphson to find displacements s.t. $\int \psi$ is minimal.

The computer graphics way has the advantage to be easy to implement. Whatever energy you chose, the point is to find a minimum.

The engineering way is harder to implement, but has a strong math background (Lax theorems to ensure convergence).

So here are my questions

Are both ways equivalent ?

The computer graphics way is an optimization problem, the engineering way is a numerical analysis problem. How do we apply the convergence theorems (stability, consistency, etc) in the energy minimization case ?


问:
3D Cooley-Tukey FFT
To compute the $N$-point DFT
$$ X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn} $$
where $N = N_1 N_2$, we can write the indices as $n = N_2 n_1 + n_2$ and $k = k_1 + N_1 k_2$, (effectively packing the data into a 2D array) and then

$$
X[k_1 + N_1 k_2] = \sum_{n_1=0}^{N_1 - 1} \sum_{n_2 = 0}^{N_2 - 1} x[N_2 n_1 + n_2] W_N^{kn}
$$
and
$$
W_N^{kn} = W_N^{(N_2 n_1 + n_2)(k_1 + N_1 k_2)} = W_N^{N_2 n_1 k_1} W_N^{N_1N_2n_1k_2}W_N^{n_2 k_1} W_N^{N_1 n_2 k_2} = W_{N_1}^{n_1 k_1} W_N^{n_2 k_1} W_{N_2}^{n_2 k_2}
$$

so
$$
X[k_1 + N_1 k_2] = \sum_{n_1=0}^{N_1 - 1} \left( \sum_{n_2 = 0}^{N_2 - 1} x[N_2 n_1 + n_2] W_{N_1}^{k_1 n_1} W_N^{n_2 k_1} \right) W_{N_2}^{k_2 n_2}
$$

This means we can compute the DFT of the columns of the 2D matrix, then multiply by twiddle factors $W_N^{n_2 k_1}$, then compute the DFT of the rows.

Does this extend to 3 dimensions? i.e. if $N =N_1 N_2 N_3$ can we pack the data into a 3D array, compute the DFT across each axis and multiply by twiddle factors? If so how would the indexing work and what would the twiddle factors look like.

Thanks for any help

Title: 3D Cooley-Tukey FFT


问:
How do you correctly implement Scipy's FFT procedures to produce a low-pass filter - image processing
I'm following this low-pass filter example in the text "Image Operators: Image Processing in Python 1st Edition" by Jason M. Kinser, but can't seem to duplicate their results.

The text's code:

```<code>import imageio
import scipy.fftpack as ft
import mgcreate as mgc

amg = imageio.imread('data/reschartsmall.png', as_gray=True)
bmg = 1 - amg/amg.max()
cmg = ft.fftshift(ft.fft2(bmg))
V,H = cmg.shape
circ = mgc.Circle((V,H), (V/2,H/2), 10)
dmg = ft.ifft2(ft.fftshift(circ*cmg))
imageio.imsave('figure.png', -dmg.real)
</code>```

where the Circle function is given by

```<code># mgcreate.py
def Circle(size, loc, rad):
    b1,b2 = np.indices(size)
    b1,b2 = b1-loc[0], b2-loc[1]
    mask = b1*b1 + b2*b2
    mask = (mask <= rad*rad).astype(int)
    return mask
</code>```

The following is the original image and the final image after the low-pass filter code above is applied:

I load in the image from
https://raw.githubusercontent.com/joefoxva1/CDS468/master/Lectures10_11_12/USAF1951w.jpg
but otherwise, I'm using the same code. This is the image I'm producing:

My image appears to be inverted (apparently fixed by saving the image as dmg.real instead of -dmg.real) and nonuniform (the two FFT procedures together are creating a band in the center of the image). I don't understand why I'm observing different behavior of the same code. My guesses were that maybe the Scipy fftpack functions have changed since the publishing of the text or that the image I'm importing isn't the same somehow (doesn't seem likely, the code produces similar behavior with other images).


问:
Is there any reliable free/open source tool for structured mesh smoothing?
I have been using Pointwise for grid generation and found the quality of smoothed grids to be stunning.

I am not aware of any free/open source alteranative that offers the same capabilities for structured grids (I am not interested in unstructured meshes).

I am aware of Meshlab but never used it (it doesn't even have undo functionality!).

Here is what I mean by smoothing (Using the elliptic solver to enhance the mesh quality): See the video starting from 8:19 : https://youtu.be/ED3J8SO7kOg?t=498

Could you please suggest any candidate free/open source software/libraries for that purpose?


问:
Is it possible to use a fixed point iteration for solving this nonlinear system?
Consider the following differential equation
\begin{align}
\frac{\partial f(u)}{\partial x} &= g(x), \ \ x\in [x_{L},x_{R}] \label{Eq2.2} \\
u(x_{L}) &= g_{1}
\end{align}
where $f(u)$ is a nonlinear function. Let $\{x_{i}\}_{i=1}^{N+1}$ be a partition of $\Omega$ with $N$ elements such that
\begin{equation}
x_{L} = x_{1}<x_{2}<\cdots<x_{N}<x_{N+1} = x_{R}
\end{equation}
where $D_{k} = [x_{k},x_{k+1}]$. The weak formulation of the problem is
\begin{equation}
\int_{D_{k}} \Bigg[ -f(u)\frac{d \phi_{i}^{k}}{d x}\Bigg] = \int_{D_{k}} g(x)\phi_{i}^{k}-\Big[(f(u))^{*}\phi_{i}^{k}  \Big]_{x_{k}}^{x_{k+1}}  \label{WeakFormulationNL}
\end{equation}
where the numerical flux $(f(u))^{*}$ is given by
\begin{equation}
(f(u))^{*} = \{f(u)\}+\frac{C}{2}[\hspace{-0.6mm}[ u ]\hspace{-0.6mm}]
\end{equation}
where $\{\cdot\}$ is the average, $[\hspace{-0.6mm}[ \cdot ]\hspace{-0.66mm}]$ is the jump and $C$ is given by $C = \max{ |f(u)| }$. After the space discretization we get a system of the form
\begin{equation}
-r(u_{h}) = b-w(u_{h})
\end{equation}
which can be rewritten as follows
\begin{equation}
r(u_{h}) + b-w(u_{h}) = 0
\end{equation}
The entries of the operators $r(u_{h})$ and $b_{k}$ in the element $D_{k}$ are given by
\begin{align}
(r_{k})_{i} &= \int_{D_{k}} f(u)\frac{d \phi_{i}^{k}}{d x} \\
(b_{k})_{i} =& \int_{D_{k}} g(x,t)\phi_{i}^{k}
\end{align}
where $w(u_{h})$ in the element $D_{k}$ is the term $\Big[(f(u))^{*}\phi_{i}^{k}  \Big]_{x_{k}}^{x_{k+1}}$.
Let $F(u_{h}):= r(u_{h}) + b-w(u_{h})$, then we have a nonlinear system of the form $F(u_{h}) = 0$. I want to rewrite this in the form $G(u_{h}) = u_{h}$ in order to solve the system using a fixed point iteration.

Here's my attempt: let $\lambda \neq 0$ then
\begin{align*}
F(u_{h}) =& 0 \implies \\
r(u_{h}) + b-w(u_{h}) =& 0  \implies \\
r(u_{h}) + b-w(u_{h})+\lambda u_{h}  =& \lambda u_{h} \implies \\
\frac{1}{\lambda}\Big[r(u_{h}) + b-w(u_{h})+\lambda u_{h} \Big]=&  u_{h}
\end{align*}
Then we define
\begin{equation}
G(u_{h}) = \frac{1}{\lambda}\Big[r(u_{h}) + b-w(u_{h})+\lambda u_{h} \Big]
\end{equation}
My first question is:
Do you think is possible to use the fixed point iteration to solve the nonlinear system?

if it is possible:
Do you think this approach is correct? or Do you know a better way to solve the system? (Without Newton's method)


问:
Fourier integral over elements
Suppose I have a triangular element with vertices ${\vec{r_1},...,\vec{r_3}}$ and a function $f(\vec{r})$. I want to calculate the fourier integral over this triangle such that:

$$F(k_x,k_y)=\int \int f(x)e^{-j2\pi(k_xx+k_yy)}dxdy \, .$$

Normally, for any integral inside such a triangle element, I would do the master element transformation and calculate the fourier transform as:

$$F(k_x,k_y)=\int_{0}^{1} \int_{0}^{1-\eta} f(x(\xi,\eta),y(\xi,\eta))e^{-j2\pi(k_xx(\xi,\eta)+k_yy(\xi,\eta))}|J|d\xi d\eta \, ,$$

where $|J|$ is the Jacobian of the transform from $(x,y)$ to $(\xi,\eta)$. In this case, I do not know if this approach is correct. Is there anyone that has used such an approach to calculate the Fourier transform of a "meshed" object before? I would also be very happy if someone could share an article on the subject.

Title: Fourier integral over elements


问:
Any literature that extensively discusses the ability to strongly inter-/extrapolate computationally on little empirical data?
Any literature that extensively discusses the ability to strongly inter-/extrapolate computationally on little empirical data?

This topic has fascinated me, but I find that it seems a bit novel.

Particulary, instead of gathering large amounts of empirical observations, is there something that could guide interpolating a lot of empirical result from very limited but "skillfully" selected or manipulated data?


问:
Pass forward intermediate results during iterative optimization
To investigate a counter-current flow heat exchanger while considering temperature dependent physical properties (such as specific heat $c_\textrm{p,i}$, heat conductivity $\lambda_\textrm{i}$, viscosity $\eta_\textrm{i}$ ...) I am trying to solve a system of equations by finding its roots using scipy.optimize.fsolve().

The temperature dependance can be summarized as

$$
\vartheta_\textrm{i,mean} = \frac{\vartheta_\textrm{i,in} + \vartheta_\textrm{i,out}}{2} \\
\eta_\textrm{i} = f(i, \vartheta_\textrm{i,mean}) \\
c_\textrm{p,i} = f(i, \vartheta_\textrm{i,mean}) \\
\lambda_\textrm{i} = f(i, \vartheta_\textrm{i,mean}) \\
\alpha_\textrm{i} = \frac{\textrm{Nu}_\textrm{i} \cdot \lambda_\textrm{i}}{d_\textrm{h}} \\
$$

and is implemented with regression polynomials.

Based on these physical properties, dimensionless numbers are calculated (Reynold's number $Re$, Prandtl's number $Pr$ and Nusselt's number $Nu$:

$$
\textrm{Re}_\textrm{i} = \frac{\dot{m}_\textrm{i} \cdot d_\textrm{h}}{A_\textrm{c} \cdot \eta_\textrm{i}} \\
\textrm{Pr}_\textrm{i} = \frac{\eta_\textrm{i} \cdot c_\textrm{p,i}}{\lambda_\textrm{i}} \\
\textrm{Nu}_\textrm{i} = \textrm{Pr}_\textrm{i}^{1/3} \cdot \textrm{Re}_\textrm{i}^{0.62} \\
$$

As heat is transferred from a hot to a cold fluid stream, I need to evaluate the beforementioned equations for each fluids which are denoted by the general index $i$:

$$
i = \{1, 2\}
$$

The above equations are then combined with each other to calculate the thermal transmittance $k$:

$$
k = \frac{1}{ \frac{1}{\alpha_\textrm{1} } + \frac{s_\textrm{w}}{\lambda_\textrm{w}} + \frac{1}{\alpha_\textrm{2}}}
$$

The transferred heat $\dot{Q}$

$$
\dot{Q} = \dot{m}_\textrm{1} \cdot c_\textrm{p,1} \cdot (\vartheta_\textrm{1,out} - \vartheta_\textrm{1,in}),
$$

can then be used to define the system of equations to be finally solved:

$$
0 = \vartheta_\textrm{2,mean} - \vartheta_\textrm{1,mean} - \Delta\vartheta_\textrm{mean} \\
0 = \dot{m}_\textrm{2} \cdot c_\textrm{p,2} \cdot (\vartheta_\textrm{2,in} - \vartheta_\textrm{2,out}) - \dot{Q} \\
0 = k \cdot A_\textrm{hex} \cdot \Delta\vartheta_\textrm{mean} - \dot{Q} \\
$$

My Python implementation looks like this (with pseudo-code elements):

```<code>import numpy as np
from scipy.optimize import fsolve

# -- stream 1
m_flow_1 = 30
theta_1_in = 40.9
theta_1_out = 70    # initially guessed value for solver

# -- stream 2
m_flow_2 = 100
theta_2_in = 105.1
theta_2_out = 90    # initially guessed value for solver

# -- heat exchanger configuration
area_heatex = 2312.5e-6
lambda_wall = 15
wall_thickness = 0.75e-3

def alpha_1(*args):
    # calculate alpha for fluid stream 1 here
    pass

def alpha_2(*args):
    # calculate alpha for fluid stream 2 here
    pass

def heat_cap_1(theta):
    # calculate heat capacity for fluid stream 1 here
    pass

def heat_cap_2(theta):
    # calculate heat capacity for fluid stream 2 here
    pass

def nusselt_1(*args):
    # calculate Nusselt's number for fluid stream 1 here
    pass

def nusselt_2(*args):
    # calculate Nusselt's number for fluid stream 2 here
    pass

def thermal_transmittance(alpha_1, alpha_2, thickness, thermal_conductivity):
    a = 1 / alpha_1
    b = thickness / thermal_conductivity
    c = 1 / alpha_2
    return 1 / (a + b + c)

def theta_mean(theta_in, theta_out):
    return (theta_in + theta_out) / 2

def delta_theta_mean(theta_2_in, theta_2_out, theta_1_in, theta_1_out):
    return theta_mean(theta_2_in, theta_2_out) - theta_mean(theta_1_in, theta_1_out)

def func(z):
    theta_1_out = z[0]
    theta_2_out = z[1]

    theta_1 = theta_mean(theta_1_in, theta_1_out)
    nu_1 = nusselt_1(m_flow_1, theta_1)
    alph_1 = alpha_1(nu_1, theta_1)

    theta_2 = theta_mean(theta_2_in, theta_2_out)
    nu_2 = nusselt_2(m_flow_2, theta_2)
    alph_2 = alpha_2(nu_2, theta_2)

    k = thermal_transmittance(alph_1, alph_2, wall_thickness, lambda_wall)
    delta_t = delta_theta_mean(theta_2_in, theta_2_out, theta_1_in, theta_2_out)

    cp_1 = heat_cap_1(theta_1)
    cp_2 = heat_cap_2(theta_2)

    q_dot = m_flow_1 * cp_1 * (theta_1_out - theta_1_in)

    f = np.empty(2)
    f[0] = m_flow_2 * cp_2 * (theta_2_in - theta_2_out) - q_dot
    f[1] = k * area_heatex * delta_t - q_dot

    return f

guessed = np.array([theta_1_out, theta_2_out])
solution = fsolve(func, guessed, full_output=True)
print(solution)
</code>```

The function func containing the system of equations is passed to fsolve in order to find a valid solution/root.

As you can see from the above equations and the code implementation, new physical properties need to be calculated based on the temperature values calculated by the solver. However, this results in wrong physical properties as the values from each iteration step do not seem to represent the iteratively calculated temperatures but rather their relative deviation (?).

Is there a way to access the iteration results to calculate the required physical properties and dimensionless numbers inside the iteration function?


问:
How do people deal with resized grid steps while numerically integrating using discrete Fourier Transform?
I am trying to simulate light propagation on python using FFT following the Fresnel diffraction equation given on Wikipedia:

The problem with this is that the output matrix from the DFT would be given in coordinates $(p,q)$, which would inversely scale with my variable propagation distance $z$. To use the new $(x,y)$ coordinates to account for the values outside of the FFT I would need to factor out at least the $z$ value from my coordinate grid to get a 1-to-1 match, but doing so would increase my grid steps with increasing $z$.

What would be the normal solution to something like this? My current idea would be to convert the $h(x,y)$ section to $h(p,q)$ and then apply the evolution again(and use $z/2$ for each evolution), which should take me back to the same grid size. However, that would require applying FFT twice, which I would prefer to prevent if possible.

Any suggestion would be appreciated.

<s>EDIT: Having thought about it a bit more, the two step fft might not work (I need to think about this more). The reason being that $p$ and $q$ being dependent on $z$ seems problematic on integration. </s>

EDIT2: Never mind, it would be fine in this case simply because each integration is fixed with fixed distance $z$. Therefore I can treat it like a constant.


问:
2D axisymmetric Magnetostatics Finite Element Method solver
I am new to the field of magnetostatics. I wish to write a 2D finite element code for obtaining the magnetic field inside a solenoid coil. I have started with a 2D code and have followed the method as described in the book by J.Bastos and N Sadowski. I managed to get the solutions and verified it from FEMM (a open source software package for magnetics and electrostatics). I have also built the code for 2d axisymmetric problem using the method described in the book by Nathan Ida. However, when I am trying to verify my code using FEMM (as there are no validation problems in the book), I am not able to match the results.

The axisymmetric problem I am trying is a square domain with side =60mm.The axis of symmetry is the bottom side of the square. A small square at the centre with 20mm side carries a current density of $10^6$ $A/m^2$. All the edges of the outer square boundaries are at A=0 where A is the magnetic vector potential. The material of the inner square with current density is one with $\mu_r$= 1 and the outer square has a $\mu_r$= 10.I am solving for the vector potential inside the bigger square.The equations I am solving are

$\vec{\nabla} \times \vec{H}=\vec{J}$

$\vec{\nabla}\cdot\vec{B}=0$

$\vec{B}=\mu \vec{H}$

I am solving for the magnetic vector potential $A$ where

$\vec{\nabla} \times \vec{A}=\vec{B}$

My questions are

I have been stuck on this for a while. Is there any paper or book which has a verification problem for the same which clearly brings out the method?

Is there a open source code for axisymmetric case that I can refer to? Or any experts that can help me in the algorithm for axisymmetric FEM for magnetic vector potential?

P.S.: I am not sure if I have explained the problem clearly. Let me know if any clarifications are required


问:
Solute transport around a solid obstacle
I am a newbie in CFD and single/multiphase flow and transport in general. As part of my quest to learn, I am trying to model solute transport around a solid object in the center of a 2D domain. The governing equation is the advection-diffusion equation, given as:

where C is the concentration, v is the flow velocity, D is the diffusion coefficient. At a constant flow velocity, the above reduces to:

This is my code so far:

```<code>clc; clear; close all;

%%%%%%%%%% Specify inputs
CircleDiam = 40;        %in pixels
Co    = 0;          %initial conc. in the domain [mol/L]
Cin   = 0.01;       %conc. of injected fluid [mol/L]

Lx    = 0.1/100;        %Length of domain [m]
nx    = 200;        %spatial gridpoints in x
dx    = Lx/(nx-1);      %Length step size [m]
Ly    = 0.05/100;       %Width of domain [m]
ny    = 100;        %spatial gridpoints in y
dy    = Ly/(ny-1);      %Length step size [m]

T     = 5/(3600*24);    %Simulation time [days]
nt    = 8000;           %shifts
dt    = T/nt;       %Time step [days]

% Flow
u = 103.68;         %Velocity in x direction [m/day]
v = 0;              %Velocity in y direction [m/day]
De = 8.64e-04;      %Dispersion coeff. [m2/day]
betaX  = u*dt/dx;
betaY  = v*dt/dy;
gammaX = De*dt/(dx^2);
gammaY = De*dt/(dy^2);

%%%%%%%%%% Create image with solid object
radius  = CircleDiam/2;

% obtain full output grids from grid vectors
[Colgrids, Rowgrids] = meshgrid(1:nx, 1:ny);

% create a logical mask for the circle by specifying the center and diameter of the circle.
centerX = 0.5 + (nx/2);
centerY = 0.5 + (ny/2);

% obtain image from: ( (y-y0)^2 + (x-x0)^2 )  <= r^2, where (y0,x0) is the centre point of circle
SolidImg = (Rowgrids - centerY).^2 + (Colgrids - centerX).^2 <= radius.^2;

% change from logical to numeric labels. Also, transpose matrix to conform with the conc. matrix
P   = double(SolidImg');
%figure, imshow(~SolidImg, [], 'InitialMagnification','fit'); box on;

% Gridblocks
x = 0:dx:Lx;
y = 0:dy:Ly;
t = 0:dt:T;
[X,Y] = meshgrid(x, y);

%specify initial conditions
C        = zeros(nx, ny, nt+1);
C(:,:,1) =  Co;              %Initial condition

%iterate finite difference equations
for k = 1:nt
    for j = 2:ny-1
        for i = 2:nx-1

            if P(i,j)==1    %Solid pixels
                %C(i,j,k+1) = C(i,j,k);
                C(i,j,k+1) = 0;
            else
                C(i,j,k+1) = C(i-1,j,k)*(betaX/2+gammaX) + C(i+1,j,k)*(gammaX-betaX/2)...
                    + C(i,j-1,k)*(betaY/2+gammaY) + C(i,j+1,k)*(gammaY-betaY/2)...
                    + C(i,j,k)*(1-2*gammaX-2*gammaY);
            end

        end
    end

    % Insert boundary conditions
    C(:,1,k+1) = C(:,2,k+1);        % bottom bc
    C(:,end,k+1) = C(:,end-1,k+1);  % top bc
    C(1,:,k+1) = Cin;           % left bc
    C(end,:,k+1) = C(end-1,:,k+1);  % right bc

%     C_all = squeeze(C(:,:,k));
%     pcolor(C_all'); shading interp; colorbar; colormap(jet); hc = colorbar;
%     ylabel(hc,'Concentration [mol/L]');
%     caxis([0 max(max(max(C)))])
%     title( sprintf('Time = %f seconds', k) )
%     pause(0.001);
%
end

</code>```

My result (see the first image) does not look like what I would expect. Ideally, I expect to see something like the above (second) image, and so I am wondering where I might be getting it wrong. I have assumed that both velocity and diffusion coefficient are zero at the solid boundary. Is there an additional boundary condition to be implemented at the solid boundary?

I would very much appreciate if someone can help me out here or point me in the right direction.

PS:

The axes are the same in both image. The difference is that mine is in pixels while the other is in cm.

Thank you in anticipation.

EDIT:

Previously, I had velocity in x = velocity in y. This is wrong because the flow is supposed to be coming from one edge only. I have now corrected this to be velocity in x = +103.68 m/day, and velocity in y = 0.

Previously, the boundary condition at the bottom edge was C(:,2,k+1) = C(:,1,k+1);. It should have been C(:,1,k+1) = C(:,2,k+1);.

The above changes means that I now have a much closer solution (see first image) to the problem than my previous version.


问:
Computationnal Mechanics : Three Bodies contact problem with ALM
I am facing an issue with convergence of a contact problem consisting of three successive bars, as presented below.

A force is applied on the right hand of the first bar so there is contact between both the first and the second bar and the second and the third bar. Another bar is linked to the second bar so there is no ill-conditionned stiffness matrix.

I implemented a Python code to solve this problem using the Augmented Lagrange Multiplier Method with Uzawa's algoritm to update the Lagrangian Multipliers, as presented in Yastrebov, "Numerical Methods in Contact Mechanics".

As it stands now, I am able to solve this problem using my code, generating the expected results.

Nonetheless, I am not satisfied with the convergence the ALM provides for this problem, and I am looking for a way to improve it. Also, the robustness of ALM for this problem is poor, in such way that a sligth modification in penalty prevents the convergence.

I am looking for a way to improve both robustness and convergence of this problem.

One way I found that worked for other problems but not for this one, is to update the contact forces before updating the displacements, and looping again on  to converge Lagrangian Multipliers. This way, the system anticipates the contact forces. For a reason I don't understand, this method works for other problems, reducing the total amount of iterations, but it prevents this particular problem from converging.

Edit :

The algorithm is constituted of two loops :
First loop regards the update of the displacements
Second loop is inner to the first one and regards the update of the Lagrange Mulipliers.
I'll use the superscipt $i$ to refer to the i-th iteration of the first loop and $j$ to refer to the j-th iteration of the second loop.

Given an initial displacement $u = 0$, the penalty $\epsilon$, the Lagrange Multiplier $\lambda=0$ and the gap function of a contact $G = x_2 - x_1$ and given the contact force : $F_c = max(\lambda - \epsilon \times G, 0)$

First Loop :

Computing the Residuals $R^i = R(u^i) = \sum\left(F_i(u^i) + F_e(u^i) + F_c(u^i) \right) $ : sum of internal, external and contact forces.

Computing a descent direction $v^i$

Second loop

Compute the step $\rho^j$ minimizing the energy $W(u+v^i*\rho^j)$

Update the Lagrange Multiplier : $\lambda^{j+1} = \lambda^j - \epsilon \times G$

If ($G - \rho^j\times v^i > 0$ and $\epsilon \times (G - \rho^j\times v^i)< Tol^j$) or $j = j^{max}$ go to 4 . Else $j = j+1 $and go to 3.

$u^{i+1} = u^i + \rho^j\times v^i$ and $x^{i+1} = x^i+u^{i+1}$

If $R(u^{i+1}) < Tol^i$ or $i = i^{max}$ end. Else $i = i+1$

I do not mind doing more iterations of the inner loop. Yet I would like to minimize the number of iterations from the first loop.