问:
Efficiently computing the properties of a Chebyshev series
Suppose we have some function $f(x)$ defined as a Chebyshev series up to order $M$:
$$ f(x) = \sum_{n=0}^{M} c_n T_n(x). $$

For a given coefficients vector $\mathbf{c}$, and $x \in [-1,1]$ I'm looking to determine the following:

Is $f(x) > 0$ everywhere?

The maximum value of the gradient $f'(x)$

The number of stationary points i.e. where $f'(x) = 0$

These properties are of course not difficult to determine, but computational efficiency is an important factor for me in this case.

I'd love any suggestions of efficient ways to answer these qustions starting only from the coefficients, and no prior calculation of $f(x)$ values.

Thanks for any help!


问:
Numerical error in implementation of iterative algorithm
I am trying to implement (in Python for now) low thrust orbit propagation for spacecraft using universal variables. For a given central body with the gravitational parameter $\mu$ and an orbit with the semi-major axis $a$ and the initial position $\vec{r}_0$ and velocity $\vec{v}_0$ at $t = t_0$ the position for a given time is given by:
$$
\vec{r} = \vec{r}_0 \ f(s) + \vec{v}_0 \ g(s)
$$
and the velocity by
$$
\vec{v} = \vec{r}_0\ \dot{f}(s) + \vec{v}_0\ \dot{g}(s)
$$

Where
$$
f(s) = 1-\left(\frac{\mu}{|\vec{r}_0|}\right) s^2 c_2(\alpha s^2)
$$
$$
g(s) = t-t_0-\mu s^3c_3(\alpha s^2)
$$

With $\alpha=\frac{\mu}{a}$, $c_n$ the n-th Stumpff function and $s$ the solution of
$$
t-t_0=|\vec{r}_0| s \ c_1(\alpha s^2) + |\vec{r}_0| s^2 \ c_2(\alpha s^2) + \mu s^3 \ c_3(\alpha s^3)
$$

I am calculating $s$ by using Newton's Method and everything works fine for scenarios without thrust. The orbit is elliptical and closed, parabolic for higher initial velocities, if I choose a $\vec{r}_0$, $\vec{v}_0$ and $t_0$ and propagate $t$ forward.

For cases with thrust where I would change the velocity in every iteration

$$
\vec{v}_{k} = \vec{v}_{k-1} + \vec{a}_{k} \cdot \Delta t
$$

it is necessary to update $\vec{r}_0$ and $\vec{v}_0$ every time, propagate $\Delta t$ forward and repeat instead of choosing initial conditions and just propagating $t$ forward. In preparation for that I chose not to add any velocity so both approaches should yield the same results.

Here is my Problem:
Even after the first iteration the differences are large and the error increases fast

<img src="https://i.stack.imgur.com/opE2X.png" alt="blue: fixed initial conditions, propagating $t$; red: propagating initial conditions, fixed $\Delta t$" />

blue: fixed initial conditions, propagating $t$

red: propagating initial conditions, fixed $\Delta t$

I thought I could solve this by using a more advanced integration algorithm like Runge-Kutta and tried to calculate every iteration in as few steps as possible but I was not able to transform the equations so that I could use Runge-Kutta (since it is not an ODE) and reducing the steps did not help at all.

Can someone help to fix this or give any hints to why this error is so large?

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import newton
import matplotlib.animation as animation
import cmath as mth

plt.close("all")

#defining some constants
gamma           = 4*np.pi**2*6.6108**3

#stumpff functions
f_c0            = lambda x: mth.cos(x)
f_c1            = lambda x: mth.sin(x)/x
f_c2            = lambda x: (1.-mth.cos(x))/(x*x)
f_c3            = lambda x: (x-mth.sin(x))/(x*x*x)

#f
f_uvf           = lambda mu,r0,s,a,sq:     1.-mu/r0*s**2*f_c2(sq)
#g
f_uvg           = lambda t,tau,mu,s,a,sq:  t-tau-mu*s**3*f_c3(sq)
#df/dt
f_uvdf          = lambda mu,r,r0,s,a,sq:   -(mu/(r*r0))*s*f_c1(sq)
#dg/dt
f_uvdg          = lambda mu,r,s,a,sq:      1.-mu/r*s**2*f_c2(sq)

def getPos(t,tau,r0,v0,mu,a):

    #absolute values of initial position and velocity
    r0a     = np.sqrt(r0[0]*r0[0]+r0[1]*r0[1])
    v0a     = np.sqrt(v0[0]*v0[0]+v0[1]*v0[1])

    #semi-major axis from gravitational paramteter and initial conditions
    a       = (2./r0a-v0a**2./mu)**-1.

    #defining function to numericly find s
    def f_root(s):
        sqrtx = mth.sqrt(mu/a*s**2)
        return r0a*s*f_c1(sqrtx)+r0a*s**2*f_c2(sqrtx)+mu*s**3*f_c3(sqrtx)-(t-tau)

    #using newton algorithm to evaluate the function
    s       = newton(f_root,2,maxiter=100,tol=1e-10)

    sq = mth.sqrt(mu/a*s**2)

    #calculating the functions
    f       = f_uvf(mu,r0a,s,a,sq)
    g       = f_uvg(t,tau,mu,s,a,sq)

    #calculating new position
    rx          = r0[0]*f + v0[0]*g
    ry          = r0[1]*f + v0[1]*g

    #absolute value of new position
    ra          = np.sqrt(rx**2+ry**2)

    df          = f_uvdf(mu,ra,r0a,s,a,sq)
    dg          = f_uvdg(mu,ra,s,a,sq)

    #calculating new velocity
    vx          = r0[0]*df + v0[0]*dg
    vy          = r0[1]*df + v0[1]*dg

    return (rx,ry),(vx,vy)

#setting initial values
initial_r   = (0.,5.)
initial_v   = (30.,0.)

#setting simulation parameters
dt          = .001
frames      = 1000

#setting up plot
scale = 5.5

fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False,xlim=(-scale, scale), ylim=(-scale, scale))
ax.grid()

#creating list with points in time to propagate t
t   = np.linspace(0,frames*dt,frames+1)

#creating list for position and velocity for fixed initial conditions while propagating t (method 1)
data_r = [initial_r]
data_v = [initial_v]

#same for propagating initial conditions, fixed delta t (method 2)
data2_r = [initial_r]
data2_v = [initial_v]

for i in range(len(t)-1):
    #evaluating with method 1
    data = getPos(t[i],0,data_r[0],data_v[0],gamma,3.5)
    #evaluating with method 2
    data2= getPos(dt,0,data2_r[i],data2_v[i],gamma,3.5)

    #adding data to lists
    data_r.append(data[0])
    data_v.append(data[1])

    data2_r.append(data2[0])
    data2_v.append(data2[1])

#formatting data
data_r = list(zip(*data_r))
data_v = list(zip(*data_v))
data2_r = list(zip(*data2_r))
data2_v = list(zip(*data2_v))

#plot data
ax.plot(data_r[0],data_r[1],'b')
ax.plot(data2_r[0],data2_r[1],'r')

#creating line for animation
line, = ax.plot([0], [0], 'o-', lw=2)

#animation loop
def animate(i):

    #accelerate animation
    di = i*10

    #animating mehtod 1
    #line.set_data([0,data_r[0][di]],[0,data_r[1][di]])

    #animating mehtod 2
    line.set_data([0,data2_r[0][di]],[0,data2_r[1][di]])

    return line,

#creating animation
ani = animation.FuncAnimation(fig, animate, interval=500, frames=frames, repeat=False)

plt.show()

RAW Paste Data
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import newton
import matplotlib.animation as animation
import cmath as mth

plt.close("all")

#defining some constants
gamma           = 4*np.pi**2*6.6108**3

#stumpff functions
f_c0            = lambda x: mth.cos(x)
f_c1            = lambda x: mth.sin(x)/x
f_c2            = lambda x: (1.-mth.cos(x))/(x*x)
f_c3            = lambda x: (x-mth.sin(x))/(x*x*x)

#f
f_uvf           = lambda mu,r0,s,a,sq:     1.-mu/r0*s**2*f_c2(sq)
#g
f_uvg           = lambda t,tau,mu,s,a,sq:  t-tau-mu*s**3*f_c3(sq)
#df/dt
f_uvdf          = lambda mu,r,r0,s,a,sq:   -(mu/(r*r0))*s*f_c1(sq)
#dg/dt
f_uvdg          = lambda mu,r,s,a,sq:      1.-mu/r*s**2*f_c2(sq)

def getPos(t,tau,r0,v0,mu,a):

    #absolute values of initial position and velocity
    r0a     = np.sqrt(r0[0]*r0[0]+r0[1]*r0[1])
    v0a     = np.sqrt(v0[0]*v0[0]+v0[1]*v0[1])

    #semi-major axis from gravitational paramteter and initial conditions
    a       = (2./r0a-v0a**2./mu)**-1.

    #defining function to numericly find s
    def f_root(s):
        sqrtx = mth.sqrt(mu/a*s**2)
        return r0a*s*f_c1(sqrtx)+r0a*s**2*f_c2(sqrtx)+mu*s**3*f_c3(sqrtx)-(t-tau)

    #using newton algorithm to evaluate the function
    s       = newton(f_root,2,maxiter=100,tol=1e-10)

    sq = mth.sqrt(mu/a*s**2)

    #calculating the functions
    f       = f_uvf(mu,r0a,s,a,sq)
    g       = f_uvg(t,tau,mu,s,a,sq)

    #calculating new position
    rx          = r0[0]*f + v0[0]*g
    ry          = r0[1]*f + v0[1]*g

    #absolute value of new position
    ra          = np.sqrt(rx**2+ry**2)

    df          = f_uvdf(mu,ra,r0a,s,a,sq)
    dg          = f_uvdg(mu,ra,s,a,sq)

    #calculating new velocity
    vx          = r0[0]*df + v0[0]*dg
    vy          = r0[1]*df + v0[1]*dg

    return (rx,ry),(vx,vy)

#setting initial values
initial_r   = (0.,5.)
initial_v   = (30.,0.)

#setting simulation parameters
dt          = .001
frames      = 1000

#setting up plot
scale = 5.5

fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False,xlim=(-scale, scale), ylim=(-scale, scale))
ax.grid()

#creating list with points in time to propagate t
t   = np.linspace(0,frames*dt,frames+1)

#creating list for position and velocity for fixed initial conditions while propagating t (method 1)
data_r = [initial_r]
data_v = [initial_v]

#same for propagating initial conditions, fixed delta t (method 2)
data2_r = [initial_r]
data2_v = [initial_v]

for i in range(len(t)-1):
    #evaluating with method 1
    data = getPos(t[i],0,data_r[0],data_v[0],gamma,3.5)
    #evaluating with method 2
    data2= getPos(dt,0,data2_r[i],data2_v[i],gamma,3.5)

    #adding data to lists
    data_r.append(data[0])
    data_v.append(data[1])

    data2_r.append(data2[0])
    data2_v.append(data2[1])

#formatting data
data_r = list(zip(*data_r))
data_v = list(zip(*data_v))
data2_r = list(zip(*data2_r))
data2_v = list(zip(*data2_v))

#plot data
ax.plot(data_r[0],data_r[1],'b')
ax.plot(data2_r[0],data2_r[1],'r')

#creating line for animation
line, = ax.plot([0], [0], 'o-', lw=2)

#animation loop
def animate(i):

    #accelerate animation
    di = i*10

    #animating mehtod 1
    #line.set_data([0,data_r[0][di]],[0,data_r[1][di]])

    #animating mehtod 2
    line.set_data([0,data2_r[0][di]],[0,data2_r[1][di]])

    return line,

#creating animation
ani = animation.FuncAnimation(fig, animate, interval=500, frames=frames, repeat=False)

plt.show()
```

ContentLicense: CC BY-SA 4.0


问:
Will penalty-augmented stiffness matrix cause numerical issues in eigenvalue analysis?
In the finite element method, we often construct the constraints of the system by adding penalty-function terms ( which often are many many magnitudes, up to $10^6$ order bigger than the largest stiffness matrix term) to the stiffness matrix by using the so-called Penalty Method (section 9.2)

On the other hand,

the eigenvalue analysis, as implemented by ARPACK, makes use of Arnoldi iteration scheme, which involves working on the stiffness matrix many times over in the iteration.

My hunch is that, the out-of-magnitude penalty term in stiffness matrix will make the whole Arnoldi iteration unstable, when it is iterated many times over, resulting in inaccuracy when we compute the eigenvalue analysis. I also don't know how to estimate the inaccuracy level .

Due to lack of competency in rigorous math, I am unable to prove my hunch. Am I right? Can anyone furnish a proof to prove or disprove my hunch?

PS:
My hunch is somewhat validated by a note in the textbook reference above ( section 9.7),

But still, I'll be interested to see a rigorous mathematical argument.

ContentLicense: CC BY-SA 3.0


问:
bound error for iterative method for solving linear system
$A$ is square and positive definite, and let $r_k = Ax_k - b$. Also let $M = \frac{1}{2}(A+A^T)$. I want to show that
$$\frac{||r_{k+1}||_2}{||r_k||_2} \le \left(1-\frac{\lambda_\min(M)^2}{\lambda_\max(A^TA)} \right)^{1/2}$$

We define $x_{k+1} = x_k + \alpha_k r_k$, where $\alpha_k$ is chosen to minimize $||r_{k+1}||_2$.

I think I found $\alpha_k$, to be specific I think that
$$\alpha_k = -\frac{r_k^TA^TAx_k}{r_k^TA^T Ar_k}$$

An earlier part of the problem has us prove that
$$\frac{v^TAv}{v^Tv}\ge \lambda_\min(M) > 0$$

I am having trouble using those two facts to get the desired result. I tried expanding $||r_{k+1}||_2^2 = r_{k+1}^Tr_{k+1}$ by using its definition and the definition of $\alpha$, but I didn't get really anywhere. Any ideas? Thanks!


问:
NURBS-based isogeometric analysis using Nutils
How can one construct non-uniform rational B-splines (NURBS) in Nutils, and use them for NURBS-based isogeometric analysis?</p>


问:
Spherical Advection Discretization (boundary nodes)
Consider the spherical advection problem: describing the conservation of a property $u$ in a closed spherical domain.

$$
\frac{\partial u}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2u\ \boldsymbol{v} \right)=0
$$

In the domain $0<r<1,\ t\geq0$, with the following initial/boundary conditions:

$$
\left.\frac{\partial u}{\partial r}\right|_{r=a}=0,
\\
u(r,t=0)=u_i
$$
Where $\boldsymbol{v}=\boldsymbol{v}(r,t)$ and $u_i$ are known. ($\boldsymbol{v}\leq0$). Other posts (here, and here) suggested discretizing the conservative form.
Forward explicit (downwind) discretization leads to:
$$
r\to i\Delta r,\ t\to n\Delta t,\ u(r,t)\to u_i^n
\\
\frac{\partial u}{\partial t}
\to
\frac{u^{n+1}_i-u^{n}_i}{\Delta t}
\\
\frac{\partial(r^2u\ \boldsymbol{v})}{\partial r}
\to
\frac{
((i+1)\Delta r)^2 u_{i+1}^n \boldsymbol{v}_{i+1}^n-
(i\Delta r)^2 u_{i}^n \boldsymbol{v}_{i}^n
}{\Delta r}
$$
And the equation becomes:
$$
\frac{u^{n+1}_i-u^{n}_i}{\Delta t}=
-\frac{1}{(i\Delta r)^2}
\frac{
((i+1)\Delta r)^2 u_{i+1}^n \boldsymbol{v}_{i+1}^n-
(i\Delta r)^2 u_{i}^n \boldsymbol{v}_{i}^n
}{\Delta r}
$$

My problems are:

I cannot use this equation at $r=1$ (last node) because i have no boundary condition to fix the ghost node at $r=1+\Delta r$.

I'm unsure how to write the equation at $r=0$. I can either use L'Hôpital to write
$$
\frac{\partial u}{\partial t}+
\lim_{r\to0}\left(\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2u\ \boldsymbol{v} \right)\right)
=
\frac{\partial u}{\partial t}+
\lim_{r\to0}\left(\frac{1}{2r}\frac{\partial^2 }{\partial r^2}\left(r^2u\ \boldsymbol{v} \right)\right)
=\\
\frac{\partial u}{\partial t}+
\frac{1}{2}\frac{\partial^3 }{\partial r^3}\left(r^2u\ \boldsymbol{v} \right)
=0
$$
And discretize the above or I can use the forward derivative and write
$$
\frac{u_1^{n}-u_0^{n}}{\Delta r}=0\therefore u_0^{n}=u_1^{n}
$$

What is the best way to discretize the problem and how do I address the first and last nodes?

Edit:

I found another post (here) with almost the same problem. One of the answers suggested an interpolation for the stencil at the last node, which i tried with unsuccessful results. The solution returns values greater than $1$ at $r=1$ and does not converge at all as the nodes increase (unstable).
The image below shows the solutions for 21, 51 and 101 spatial nodes. The CFL number is below $1$. Solution should stabilize at $0.91$, and be less than $1$ for every instant.


问:
Metropolis algorithm and thermal sine-Gordon model
I try to simulate thermal version of 1D $(x, t)$ sine-Gordon field model. I am interested in finding thermal static solution that minimizes functional of energy $E$:

$$E = \int dx \left( \frac{1}{2} \phi' ^2 + 1 - \cos \phi \right)  \ ,$$

where $\phi' = \partial_x \phi $

What is very confusing that acceptance ratio of Metropolis algorithm is too high - more than $0.95$, so almost every new proposed configuration is accepted. On each Metropolis step I change field value at one spatial point and calculate difference in energy. To propose new configurations uniform sampling is used with step parameter $\delta = 0.5$, i.e.

$$\phi_{new} = \phi_{old} + r \ ,$$

where $r$  is random number from $\phi_{old} - \delta$ to $\phi_{old} + \delta$.

It seems to me that if acceptance ratio is too high then algorithm does not work correctly. However, such algorithm is perfectly applied to the ground state harmonic oscillator (via path integral Monte Carlo).

Here is my code:

```<code>double E_part(double dx, double phi, double phi_plus, double phi_minus) {
return dx * ( 0.5 * pow ( 0.5 * ( phi - phi_minus ), 2.0 ) + 0.5 * pow ( 0.5 * ( phi_plus - phi ), 2.0 ) + 1 - cos ( 0.5 * ( phi_minus + phi ) ) - cos ( 0.5 * (phi + phi_plus ) ) );
}

double Metropolis(int N, int Steps, double dx, double Temperature, double* dphi, double* phi, double delta, double* EnergyArray) {

double Beta = 1.0 / Temperature;

double Epart;

double Enewpart;

double phinew;

double phi_plus, phi_minus;

int AcceptanceCounter = 0;

 //Initial energy

double r; //random number from 0 to 1;

for (int k=0; k<Steps; k++) {

    for (int j=0; j<N; j++) {

        int pos = RandomUniformInt(1 , N-2); //random position to change;

        phinew = phi[pos] - delta + RandomUniform() * 2 * delta;

        phi_plus = phi[pos + 1];
        phi_minus = phi[pos - 1];

        Epart = E_part(dx, phi[pos], phi_plus, phi_minus);
        Enewpart = E_part(dx, phinew, phi_plus, phi_minus);

        r = RandomUniform();

        if ( Enewpart * Beta - Epart * Beta < 0.0 || exp( Epart * Beta - Enewpart * Beta ) >= r ) {
            phi[pos] = phinew;
            dphi[pos - 1] = 0.5 * ( phi[pos] - phi[pos - 1] );
            dphi[pos] = 0.5 * ( phi[pos + 1] - phi[pos] );
            AcceptanceCounter++;

        }
    }

   EnergyArray[k] = Energy(N, dx, dphi, phi);

}
return ( AcceptanceCounter / static_cast<double>(Steps * N) );  //AcceptanceRatio is returned
}
</code>```

I don't know what might be wrong.


问:
Adjoint of the MATLAB $\tt dwt3$ (3D wavelet transform) operator
How do I compute the adjoint of MATLAB's dwt3 operator?

In other words, how do I compute the adjoint of the linear operator that takes a 3D complex array x as input and returns the quantity dwt3(x,'db4') as output.

This adjoint computation is an important step in many iterative reconstruction algorithms.

One might guess that idwt3 is the adjoint of dwt3, but the following code seems to show that this is not the case (unless my code contains an error):

```<code>sz = [384,300,144];

% Randomly generate a 3D complex array
x1 = randn(sz) + 1i*(randn(sz));

% Compute wavelet transform of x1
Wx1 = dwt3(x1,'db4');

% Randomly generate x2 which has the same size and shape as Wx1 (which is a Matlab structure)
x2 = Wx1;
for i1 = 1:2
    for i2 = 1:2
        for i3 = 1:2
            arrSz = size(x2.dec{i1,i2,i3});
            x2.dec{i1,i2,i3} = randn(arrSz) + 1i*randn(arrSz);
        end
    end
end

% Compute inverse wavelet transform of w2
WTransx2 = idwt3(x2);

% Check to see if the adjoint property <x1,WTransx2> = <Wx1,x2> is satisfied
prod1 = x1(:)'*conj(WTransx2(:));
prod2 = 0;
for i1 = 1:2
    for i2 = 1:2
        for i3 = 1:2

            Wx1_arr = Wx1.dec{i1,i2,i3};
            x2_arr = x2.dec{i1,i2,i3};
            prod2 = prod2 + Wx1_arr(:)'*conj(x2_arr(:));

        end
    end
end

max(abs(prod1 - prod2))
</code>```


问:
$(max(0, f(x)))^2$ or $(max(0, exp(f(x))))^2$ for soft constraints with Gauss-Newton
I need some kind of inequality constraint in my optimization problems (rude version of SVM for example or skeleton based mesh fitting). However hard constraints is not suitable for me because sometimes small constraint violations are allowed.

I think about $(max(0, f(x)))^2$ as soft constraint function where $f(x)$ is constraint. I.e. i want to keep $f(x) \le 0$. $f(x)$ is linear function.

At first glance this function is ok because it provide some kind of barrier for constraint and it can be expressed as sum of squares as Gauss-Newton solvers need.
But I'm not sure that Gauss-Newton can hold such kind of function. Because expression inside squares (max(0, f(x))) is not smooth.

I think about $(max(0, exp(f(x))))^2$ as "barrier". This function is smooth (because $max(0, exp(f(x)))$ is smooth).

So my question is:

Is ok to use $(max(0, exp(f(x))))^2$ or $(max(0, f(x)))^2$ within Gauss-Newton for smooth constaint? If not which kind of function i should use?

PS Sorry for my poor English and math skills.


问:
Cavity Flow CFD Boundary conditions and strange waves
So I have a PDE that I use to describe how material flows through a volume(2D or 3D).
$$\frac{\partial C}{\partial t} + \vec{u} \cdot \nabla C = (D' + D )\nabla^2C$$
Now using finite differences I get to the following:
$$C^{i+1}_{j,k} =C^i_{j,k}+\Delta t \biggr[ (D+D')\biggr[\frac{C^i_{j+1,k}-2C^i_{j,k}+C^i_{j-1,k}}{(\Delta x)^2}+\frac{C^i_{j,k+1}-2C^i_{j,k}+C^i_{j,k-1}}{(\Delta y)^2}\biggr]-u_x \frac{C^i_{j,k} -C^i_{j-1,k}}{\Delta x} - u_y \frac{C^i_{j,k} -C^i_{j,k-1}}{\Delta y}\biggr]$$
where D is diffusion coefficient and superscript is temporal, subscript is spatial. (Just imagine a grid). We are given $u_x$ and $u_y$ through some already validated code. Now, we use a boundary condition where $C=0$ at the edges. I know this is a silly boundary condition but it's irrelevant because the problem I'm about to show you persists even if we use small times.(where the concentration doesn't reach the edge)

[Look in the comments...It's telling me I can't post 2 links... Where do I have the second?]

here is the code. It's relatively simple.

```<code>import numpy as np
import matplotlib.pyplot as plt
import time
%matplotlib inline
from scipy import ndimage
import scipy.io as sio
a =sio.loadmat('C:\\Users\\Tsila elbag\\Documents\\IPythonNotes\\New folder\\Cavity10k_1.mat')
b = a['VelField']
N =50
M =300
ux = b[0][0][0]
uy = b[0][0][1]
#ux = np.zeros((50,50))
#uy = np.ones((50,50))
dx = 2/N
dy =2/N
D= .001
Dprime =0
dt = .0001
C = np.zeros((N,N))
#C[np.floor((.25+1)/dx):np.floor((.75+1)/dx),np.floor((.25+1)/dy)\
#  :np.floor((.75+1)/dy)]=1
C[np.floor((.25+1)/dy)\
  :np.floor((.75+1)/dy),.1/dx:np.floor((.25)/dx)]=1
f, axarr = plt.subplots(2, sharex=True)
axarr[0].contourf(C)
Chist= C
Ctemp = C
for T in range(M):
    for j in range(N):
        for k in range(N):
            if k==49  or j ==49 or k ==0 or j ==0:
                Ctemp[j,k]=0
            else:

                Ctemp[j,k]= C[j,k]+dt*((D+Dprime)*(((C[j+1,k]-2*C[j,k]+\
                C[j-1,k])/(dx**2))+((C[j,k+1]-2*C[j,k]+C[j,k-1])/(dy**2)))-\
                ux[j,k]*(C[j,k]-C[j-1,k])/(dx)-uy[j,k]*(C[j,k]-C[j,k-1])/dy)
    #Chist.append(list(C))
    C = Ctemp
axarr[0].set_title('Original') axarr[2].contourf(C)
axarr[2].set_title('{} seconds later'.format(dt*M))
</code>```

So the question is why do I have these waves and if anyone knows effective boundary conditions I would be grateful. Thanks in advance.

OwnerDisplayName: Aviv Moshe Albeg


问:
Nelder-Mead algorithm with constraints
Im currently implementing a Nelder-Mead algorithm in MATLAB that can deal with constraints. In detail, the constraints should be added to the objective function with a Lagrange Multiplier($\lambda$) and the new function is then minimized by Nelder-Mead. Looks something like this.

$$L(x_1, \ldots,x_n,\lambda_1, \ldots,\lambda_n) = f(x_1,\ldots,x_n)+ \lambda\left(c-g(x_1,\ldots,x_n)\right)$$

this works for constraints that are equations, $g(x_1,\dots,x_n)= c$

Questions:

How does it work with constraints that are upper and lower bounds? cant put them into the objective function like that

Are the multipliers to be treated as additional variables by the nelder mead or are they constants you need to set?

Thank you guys.

Title: Nelder-Mead algorithm with constraints


问:
Numerical integration using interval arithmetic, nowadays
Is there now a package for rigorous numerical integration that uses interval arithmetic and has access to a well-developed library of special functions?

By "well-developed", I mean something that, at the very least, includes the error function $\text{erf}$, the incomplete gamma function $\Gamma(s,x)$ and the like. In other words, I would like to integrate expressions involving these functions. VNODE-LP doesn't seem to satisfy this requirement, as it relies on older interval-arithmetic libraries lacking such basic special functions. (Besides, perhaps it is just me, but VNODE-LP actually seems harder to install than three years ago; I am not sure it is being adequately maintained.)

PS. An ideal library would have both arbitrary-precision and double-precision modes, the latter one fast. That isn't really an issue for what I am doing now, but it sounds like a basic desideratum, just as, say, an intuitive interface written in C++ (as opposed to plain C) would be nice.


问:
Eigenvalues of a Laplacian operator on an irregular mesh
I have the following setting:

An irregularly-shaped domain, expressed as a mesh of points

A Laplacian operator, together with boundary conditions

I am looking for the eigenvalues of that operator, i.e. functions (expressed on the mesh, though clearly approximating an underlying continuous function) satisfying the boundary conditions and
$$
\triangledown^2 f = \lambda f
$$
for some function $f$ and constant $\lambda$. I don't know $\lambda$ a priori, in fact part of the problem is to find the possible values!

What I got so far for the special case, where the boundary conditions are a constant (e.g. 0s):

The Laplacian operator, when acting on a discrete mesh of points, becomes just a linear operation - e.g. a matrix transformation

So I can just look for the eigenvalues of that

That works pretty well, but doesn't generalise so well, I think. E.g. now I would like to solve the problem with boundary conditions such that on the boundary
$$
\frac{ \partial f }{ \partial n } = 0
$$
and now I  don't think I can just express that as a matrix eigenvector problem - or maybe I can?

For context, I am solving the problem of drum vibrations on an irregular-shaped domain. Given a domain, the solutions to the problem above give the shape of vibrations of the drum, along with the associated frequencies (which can be deduced from the eigenvalues). Here the boundary conditions is $f=0$ on boundary, and it solves nicely - I can produce a pretty picture like this, for example (here my domain is a semi-circle):

However, now I would like to extend this to water waves in a cyllindroid container. The boundary condition becomes that of partial derivatives vanishing at the boundary described above - and I feel out of ideas. All help will be highly appreciated.


问:
Looking for C++ function for performing optimization of parameters for multivariable function
I am adapting a Java program from C++ and need a C++ function to perform the same task as the Java BOBYQAOptimizer() function. Can anyone recommend a C/C++ library with equivalent or similar functions and good documentation on how to use them? I've found one that might be suitable ( http://dlib.net/optimization.html ), but it's not entirely clear if its optimizers can be used for this task.

The program is intended to find the set of parameters that give the maximum statistical likelihood of a phylogenetic tree of HIV virus data across thousands of patients. The parameters are infection rates between demographic groups. (Different versions use either 4 or 8 groups, giving 16 or 64 parameters, though in practice some are held constant at 0.)


问:
Quasi-Newton Optimization with parallel function evaluation
I have a function of many variables (~200-2000) which I am optimizing with some success using L-BFGS. While the function is expensive to evaluate, the gradient can be computed with not much additional effort.

I would like to know if there are any techniques which could allow this optimization to be sped up by parallel evaluation of the function. On this problem, L-BFGS averages ~1.1-1.3 evaluations per iteration, so I think to a good approximation, I can ignore the line search component of the algorithm. From this perspective then, on every iteration, L-BFGS chooses a tentative input evaluates the function and its gradient, updates its estimation of the inverse Hessian and repeats.

I imagine it is possible to enhance this approach by choosing multiple tentative inputs and using the multiple resulting gradients to form a better estimation of the Hessian. I imagine there is diminishing returns of some sort, the point selected by L-BFGS probably has the highest expected information, but surely there is some additional information about the Hessian that one can infer by using more function evaluations per iteration.

Has this type of approach been investigated? I know there are many resources for parallel gradient-free optimization, as well as parallel gradient estimation, but I have not come across anything like this yet.


问:
Modified convection diffusion reaction equations
I am trying to solve convection diffusion type equations by a modification:
$$-\epsilon \Delta u+ b \nabla u+cu=f_1,$$
with the following boundary conditions:
$$u=g~ on ~\Gamma_D,~~~~~~~~~~~~~~ \epsilon \frac{\partial u}{\partial n}=g_N ~~on ~~\Gamma_N.$$

The idea is  to solve another equation:
$$-\nu\Delta u+ \nabla u+\frac{c}{b}u=f_2,$$
$$u=g~ on ~\Gamma_D,~~~~~~~~~~~~~~ \epsilon \frac{\partial u}{\partial n}=g_N ~~on ~~\Gamma_N.$$

where $\nu=\frac{\epsilon}{b}$
and then use $\nu+\delta$ instead of $\nu.$
so we should solve
$$$$
$$-(\nu+ \delta) \Delta u+ \nabla u+\frac{c}{b}u=f_3,$$
$$u=g~ on ~\Gamma_D,~~~~~~~~~~~~~~ \epsilon \frac{\partial u}{\partial n}=g_N ~~on ~~\Gamma_N.$$
Here
$(\nu+\delta)=\frac{h}{2}*coth(\frac{h}{2*\nu})$, where $h$ is the mesh size.

Now the question is that what the relations are between $f_1,f_2,f_3.$ I used $f_2=f_3$ and $f_2=\frac{f_1}{b}$. Is it true?

Also if $b$ be a vector I used $||b||_\infty$, but the results are acceptable just for some examples. Am I wrong in this way?


问:
Dielectric in an external electric field
I want to numerically solve for electric field of an object in an external electric field. The object has linear displacement field such that its permitivity ($\epsilon$) is constant. The object's shape is not important but I would prefer an ellipse. I do not know how should I impose curl and divergence equations into simulation as well as boundary conditions.

What I need is electric field or charge density. Does anyone know how can I calculate charge density of a dielectric ellipse in an external electric field?


问:
trigonometric interpolation with non-equidistant sampling
I want to interpolate a periodic function on a non-equidistant grid and have implemented it the using the Lagrangian formula described in wikipedia. For an odd number of data points, this takes the form

$$p(x) = \sum_{k=0}^{2K} y_k\,t_k(x)$$

with

$$t_k(x) = \prod_{m=0,m\ne k}^{2K} \frac{\sin\frac12(x-x_m)}{\sin\frac12(x_k-x_m)}$$

which bears strong resemblance to the usual Lagrangian interpolation formula, but with sines and a the factors 1/2 appearing. This formula makes perfect sense and in fact, I used the equivalent of the barycentric form described in the wikipedia page of Lagrange polynomials (I cannot post two links):

$$ p(x) = \frac{\sum_{k=0}^{2K} \frac{w_k}{\sin{\frac{1}{2}\left(x-x_k\right)}}y_k}{\sum_{k=0}^{2K} \frac{w_k}{\sin{\frac{1}{2}\left(x-x_k\right)}}} $$

with

$$w_k = \frac{1}{l'\left(x_k\right)}$$

and

$$l\left(x\right) = \prod_{k=0}^{2K} \sin{\frac{1}{2}\left(x-x_k\right)}$$ .

This works magnificently, as can be seen in the following figure, where I use 15 (red) points of

$$f\left(x\right) = \sin{\left(x\right)} + \frac{1}{4}\cos{\left(2x\right)} + 2 \sin{\left(6x\right)} + 3 \cos{\left(7 x\right)}$$ ,

defined on a non-equidistant grid (in this case with more points on the edges of the region $2\ldots \pi$). As Nyquist states, this is indeed enough information to interpolate $f\left(x\right)$, which I illustrate on 100 (blue) points, also non-equidistant (but in this case with more points in the middle of the region).

However, for an even number of points, there is an additional degree of freedom. In the Wikipedia entry for trigonometric interpolation, it is stated that this is usually employed to make sure the highest modal content only has a cosine, and no sine.

This, they state, is implemented by using additional factors

$$\alpha_k=\sum_{m=0,m\ne k}^{2K-1} x_m$$

that modify the formulas to

$$t_k(x) = \frac{\sin\frac12(x-\alpha_k)}{\sin\frac12(x_k-\alpha_k)}\prod_{m=0,m\ne k}^{2K-1} \frac{\sin\frac12(x-x_m)}{\sin\frac12(x_k-x_m)}$$

This, however, I have tried, and I have found it does not work. Not only does the interpolated function not go through the correct points, but also it blows up to high absolute values.

Therefore, I have investigated it in more detail, and have found in fact that this definition of $\alpha_k$ contains singularities.

Take for example the case with 4 points

$$x_k = \frac{k \pi}{4} \ , k = 0\ldots 3$$

applied to a simple sine. Note that according to Nyquist, 3 points would be enough, so there is indeed a redundancy of 1 point.

Let's again look at the interval $x= 0\ldots 2 \pi$. As the only non-vanishing components of $y_k$ are 1 and 3 (with $y_1 = 1$ and $y_3=-1$), the necessary factors $t_k$ are given by

$$t_1 = \sqrt{2} \sin{\frac{x}{2}} \cos{\frac{x}{2}} \left( \cos{\frac{x}{2}} + \sin{\frac{x}{2}} \right ) \frac{\sin{\frac{1}{2}\left(x-\alpha_1\right)}}{\sin{\frac{1}{2}\left(\frac{\pi}{2} - \alpha_1\right)}} \\
t_3 = \sqrt{2} \sin{\frac{x}{2}} \cos{\frac{x}{2}} \left( \cos{\frac{x}{2}} - \sin{\frac{x}{2}} \right ) \frac{\sin{\frac{1}{2}\left(x-\alpha_3\right)}}{\sin{\frac{1}{2}\left(\frac{3\pi}{2} - \alpha_3\right)}}$$

where it can be seen that there are extra degrees of freedom $\alpha_k$.

However, now using their definition:

$$\alpha_1 = \frac{5}{2} \pi \\
\alpha_3 = \frac{3}{2} \pi $$ ,

this clearly leads to singularities in the denominators of $t_1$ and $t_3$. Does anybody know why this happens? There must either be something wrong in the definition of $\alpha_k$ or in my understanding of it.

I have tried to find other resources but the only thing I could find limited the discussion to an odd number of points (which indeed works fine.)

I have also tried to follow the arguments in the wikipedia page and I am not able to capture why they reach that expression for $\alpha$, repeated here:

$$\alpha_k=\sum_{m=0,m\ne k}^{2K-1} x_m$$

By the way, as a side note, it can be seen that throwing away the fourth point $\frac{3}{2}\pi$, the remaining three points indeed lead to the correct interpolation:

$$p\left(x\right) =  \sum_{k=0}^{2} y_k\,t_k(x) = t_1 = \frac{\sin{\left(\frac{x}{2}\right)}}{\sin{\left(\frac{\pi}{4}\right)}}\frac{\sin{\left(\frac{x}{2}-\frac{\pi}{2}\right)}}{\sin{\left(-\frac{\pi}{4}\right)}} = 2 \sin{\left(\frac{x}{2}\right)}\cos{\left(\frac{x}{2}\right)} = \sin{\left(x\right)}$$

EDIT

In fact, the interpolating formula for $p\left(x\right)$, using 4 points, can be written readily as

$$ p\left(x\right) = 2 s c \left[ \left(s + c \right) \frac{s c_1 - c s_1}{c_1-s_1} - \left(c-s\right) \frac{s c_3 - c s_3}{s_3 + c_3} \right]$$

with $c = \cos{\left(\frac{x}{2}\right)}$, $s = \sin{\left(\frac{x}{2}\right)}$, $c_i = \cos{\left(\frac{\alpha_i}{2}\right)}$ and $s_i = \sin{\left(\frac{\alpha_i}{2}\right)}$ as a short-hand notation.

With basic trigonometric identities $\sin{\left(x\right)} = 2 \sin{\left(\frac{x}{2}\right)} \cos{\left(\frac{x}{2}\right)}$, $\cos{\left(x\right)} = \cos^2{\left(\frac{x}{2}\right)} - \sin^2{\left(\frac{x}{2}\right)}$, this can be rewritten as

$$p\left(x\right) = \sin{\left(x\right)} \left[  1 + \frac{1}{2} \cos{\left(x\right)} \left( \frac{s_1+c_1}{s_1-c_1} + \frac{s_3-c_3}{s_3+c_3}  \right)\right] $$

So clearly, to avoid having a term $\propto \sin{\left(2 x \right)}$, the condition

$$\frac{s_1+c_1}{s_1-c_1} + \frac{s_3-c_3}{s_3+c_3} = 0$$

needs to hold. But, again, with $\alpha_1 = \frac{5}{2} \pi$ and $\alpha_3 = \frac{3}{2} \pi$, we get

$$ s_1 = -\frac{\sqrt{2}}{2} \ , \ c_1 = -\frac{\sqrt{2}}{2}  \\
s_3 = \frac{\sqrt{2}}{2} \ , \ c_3 = -\frac{\sqrt{2}}{2}$$

i.e. singularities in both denominators...

EDIT 2
Above condition can be rewritten to

$$\frac{\cos{\frac{1}{2}\left( \alpha_1 - \alpha_3 \right)}}{\left(\sin{\left(\frac{\alpha_1}{2}\right)}-\cos{\left(\frac{\alpha_1}{2}\right)}\right) \left(\sin{\left(\frac{\alpha_3}{2}\right)}+\cos{\left(\frac{\alpha_3}{2}\right)}\right)} = 0$$

which is satisfied for $\alpha_1-\alpha_3 = (1+2 r) \pi$.

For our simple test case, this is indeed correct, but at the same time they give infinities in the denominator, which becomes

$$\left(\sin{\left(\frac{\alpha_1}{2}\right)}-\cos{\left(\frac{\alpha_1}{2}\right)}\right)^2 = 0^2$$

But this is probably not always the case. Indeed, I have run interpolations where there were no problem. Would it just be due to an unfortunate choice of sampling points?

MY OWN, SLIGHTLY UNSATISFYING ANSWER

After some more thought I think that

One of the formulas on Wikipedia was slightly wrong. I have corrected it there, but the result is still $\alpha_k = \sum_{j=0, j\ne k}^{2K-1} x_j$.

This formula means nothing more than set the phase of the contribution of this point to the highest frequency to zero, i.e. it should be a cosine. Logically, another option is to set it to $\pi$ in order to get sines at the highest frequency. Or anything else. There is no clear preference

When forcing the phase to be some value, it can often happen that the denominators in the formulas contain zeros, leading to infinities. This is the case when $\sin{\frac{1}{2}\left(x_k - \alpha_k\right)} = 0$. It is very important, therefore, to have the $\alpha_k$'s satisfy the condition $x_k = \alpha_k \ne m \pi$ with $m$ an integer. This can impede the preference for the phase of the highest frequency.

Sometimes it is clear that the phase of the highest frequency does not matter as the amplitude should be zero, as in the case for my simple sine. There is, however, imho no criterion that can be used there, as it requires information about the function values $f_k$ and not only the values of the points $x_k$.

In particular, in the case for a simple sine, a phase of 0 is not possible, as this leads to infinities. An other phase has to be used, but it does not matter which: Due to the fact that $\alpha_1 -\alpha_3 = \pi$ and $y_1 = - y_3$, the global contribution to the highest frequency is zero, as required.

It is better just to take an odd number of points.


问:
How to reuse permutation-orderings within scipy's SuperLU-wrapper?
i'm solving sparse linear equations within scipy 0.18 which internally resorts to SuperLU (after umfpack got removed due to license-issues).

Current, i'm doing a complete re-factorization in each iteration of my algorithm, but now i want to re-use the permutation-ordering (at least the column-permutation; but maybe even row-permutation) calculated in the first iteration.
How can i achieve this in a nice way using scipy?

I got no previous usage-history of SuperLU (and also not much practical experience with sparse-algebra packages), but the user-guide makes me believe, that this should be done through the options-parameter.

Fact
  Specifies whether or not the factored form of the matrix A is supplied on entry, and if not,
  how the matrix A will be factorized base on the previous history, such as factor from scratch,
  reuse Pc and/or Pr, or reuse the data structures of L and U. fact can be one of:

– DOFACT: the matrix A will be factorized from scratch.

– SamePattern: the matrix A will be factorized assuming that a factorization of a matrix
  with the same sparsity pattern was performed prior to this one. Therefore, this
  factorization will reuse column permutation vector perm c.

– SampPattern SameRowPerm: the matrix A will be factorized assuming that a factorization
  of a matrix with the same sparsity pattern and similar numerical values was performed
  prior to this one. Therefore, this factorization will reuse both row and column
  permutation vectors perm r and perm c, both row and column scaling factors Dr and
  Dc, and the distributed data structure set up from the previous symbolic factorization

As an alternative, there is also the following within options:

ColPerm
  Specifies how to permute the columns of the matrix for sparsity preservation.

– NATURAL: natural ordering.

– MMD ATA: minimum degree ordering on the structure of AT A.

– MMD AT PLUS A: minimum degree ordering on the structure of AT + A.

– COLAMD: approximate minimum degree column ordering

– MY PERMC: use the ordering given in perm c input by the user.

Questions

Which part of SuperLU should i use to achieve my goal? One of the above? Something else?

How can it be done within scipy?

The ColPerm-parameter is missing MY_PERMC in scipy's docs (maybe it's not usable)


问:
Using Line Search Method for Constrained Optimization
Suppose we have a $f(x)$ to be minimized (we only know that $f(x)$ is three-differentiable), and a feasible, convex set of $S$ such that all $x$ belong to $S$. Using line search method, how we can "mathematically" prove that a descending direction pointing from previous $x$ towards a new $x$ within the $S$ exists?

And, how we can find the best $a\geq 0$ such that $x + ad$ belongs to $S$? I mean if we have a direction, how can we optimize the $d$ to have the least iterations?

1.I'm neither a beginner nor an expert, somewhere in between.
2.I need all the iterates to be feasible.

Many thanks in advance.


问:
What's the optimal method to solve for the top eigenvectors of a very large, real, symmetric matrix of limited rank?
Consider a real symmetric matrix of dimension N~10^5 and rank m~2000.  What is the most efficient algorithm for determining the top m eigenvectors?  If the answer isn't obvious, are there existing routines for this in languages like MATLAB or Python?</p>


问:
Useful Quantity for Heat Equation?
I'm interested in testing some algorithms on the heat equation, and I'd like to assess their accuracy. When evolving a Hamiltonian system, one has the energy to check the validity/correctness of the algorithm, which turns out to be (arguably) the most important quantity to check. Is there a similar quantity for the heat equation - i.e. a quantity which would very well characterize the error of the algorithm?

I could compute the time-evolution of the "energy" of the heat equation easily by the following

$$u_t(x,t)=\hat{H}u(x,t),~~~\hat{H}=\frac{\partial^2}{\partial x^2}$$

The "energy" would then be the expectation value of the $\hat{H}$ operator, which is simple to evaluate in Fourier space. However, I don't know if this "quantity" would even be useful.

Title: Useful Quantity for Heat Equation?


问:
Algorithm for Adaptive Mesh Refinement
I am trying to implement Adaptive Mesh Refinement. I am not a Mathematics/Computational Science person so I will try to write the algorithm in a simpler way. I will be grateful if experts can comment on the correctness of steps:

Aim: To solve Laplace Equation with Dirichlet boundary conditions i.e.
$\nabla^{2} u = 0 $ and $u = 1$ on $\partial \Omega$. Assume 3 levels of static refinement (say for 2-D mesh on $[0,1]$ I refine cells with $x^{2}+y^{2}>0.7$ and $x^{2}+y^{2}<0.9$, where $x, y$ are coordinates).

Initialize $u = 0$ on all levels.

Solve on coarsest level (level 1) using 5-pt stencil Jacobi iterations. The initial value of $u$ at level 2 is used as Coarse-Fine interface value.

Exchange ghost cells for parallel implementation at coarsest level.

Interpolate coarse grid cells at level 1 to fill ghost cells at level 2.

Solve on level 2. The initial value of $u$ at level 3 is used as Coarse-Fine interface value (and also the interpolated values in ghost cells).

Exchange ghost cells for parallel implementation at level 2.

Repeat same at level 3.

Now restrict from level 3 to 2 and level 2 to 1 (so that coarser grids see some part of the solution at finer grids).

Repeat these steps.

Are these steps correct ?

(In 3-D my mesh looks like the following)

Title: Algorithm for Adaptive Mesh Refinement


问:
Multigrid for Robin boundary conditions
I have a question regarding the treatment of Robin boundary conditions in a multigrid solver. I am solving the Poisson equation in $\Omega=(0,1)^2$ with Robin boundary conditions on the boundary,
$$- \Delta u = f \text{ in } \Omega$$
and
$$\frac{\partial u}{\partial n} - au= g \text{ on } \partial \Omega.$$
I use a finite difference discretization with central differences (a standard five point stencil) for the interior points, and eliminate the boundary conditions so that, for example, at a point $x_C$ on the left boundary with east, south and north neighbors $x_E, x_S, x_N$ the resulting equation is

$$\frac{1}{h^2}\left( 4u(x_C)-u(x_S)-u(x_N)-2u(x_E)\right) - 2\frac{a}{h}u(x_C) = \frac{2}{h}g(x_C)+f(x_C).$$

In chapter 5 of the book 'Multigrid' of Trottenberg et. al., it is explained that if this approach is followed for a Neumann problem ($a=0$) and the full-weight restriction operator $R$ is modified on the boundary points, e.g., on a vertical boundary the modified stencil is
$$ \frac{1}{16}\left[\begin{array}{cc} 2 &2 \\ 4&4 \\2&2 \end{array}\right]_{h}^{2h} $$
then a relaxation method can be applied directly on the equations for the boundary points with the eliminated boundary condition (no need of transferring the equation and boundary condition separately). I am using the interpolation operator $I=4R^T$ and Galerkin coarse grid systems.

My question is: Is this multigrid approach correct (modified FW operator on the boundaries + relaxation on the equation with eliminated boundary conditions and Galerkin coarse systems) also for more general Robin boundary conditions, provided a good relaxation method is used?

Thanks in advance.

Title: Multigrid for Robin boundary conditions


问:
Under what circumstances does Elemental's distributed SVD not work?
I am playing around with Elemental's distributed singular value decomposition and am running into two particular issues.

Building the test at tests/lapack_like/SVD.cpp, and running with

```<code>mpirun -np 4 ./SVD --height 3 --width 3
</code>```

I get

```<code>SVD residual was unacceptably large
</code>```

for the distributed Complex<float> test (but sequential is OK).

Additionally, when I run with --height 50 --width 52, I run into some problem

```<code>H and A must have the same height
</code>```

The Elemental documentation for SVD is outdated, and I don't know what requirements there are on the dimensions of a matrix in order to take its SVD.

I have Elemental built with pulled OpenBLAS, Scalapack, Metis and Parmetis. I could not manage to get libFLAME linked, even using the reference BLAS/LAPACK implementations.


问:
Finite Difference advection-reaction-diffusion in spherical coordinates, problem with Diffusion
I have a problem with the following equation:

$$
V_r\frac{\partial C_A}{\partial r} + \frac{V_\theta}{r}\frac{\partial C_A}{\partial \theta} = \frac{2}{Pe_A} \left[ \frac{\partial ^2 C_A}{\partial r ^2} + \frac{2}{r} \frac{\partial C_A}{\partial r} - k C_A\right]
$$

I'm trying to solve this with finite difference technique. My discretization :

$$
V_r \frac{A_{i+1,j}-A_{i-1,j}}{2\Delta r} + \frac{V_\theta}{r} \frac{A_{i,j+1}-A_{i,j-1}}{2\Delta \theta} - \frac{2}{Pe} \left [ \frac{A_{i+1,j}-2A_{i,j}+A_{i-1,j}}{\Delta r ^2}+\frac{2}{r} \frac{A_{i+1,j}-A_{i-1,j}}{2\Delta r} - kA_{i,j}\right] = 0
$$

I first tried it with this, but because I used central-difference in doesn't work for $Pe(cell)>2$, I had strange oscillations at the boundary. So I thought I use a upwind scheme now.

$V_r$ and $V_\theta$ are given through polynoms as a function ($r, \theta$).

So for $\theta<90 °$ I Use second order upwind for the convection term $V_r > 0$.

And for $\theta >90 °$ in the other direction.
$V_\theta$ is always $<0$, so I only use this form for this.

My problem is, for high Peclet number and low reaction constants, the concentration is rising up where it shouldn't. This happens in a region, where the convection term gets stronger.

What am I supposed to do ? Is there a difference about using upwind in polar coordinates ?

Edit :

I want go more in detail. This is my coordinate system :

The problem appears at the zero degress line. I also have seen that the problem is strongly dependent on $V_\theta$. If I multiply $V_\theta$ with 0.9 for example there is no problem. The concentration is not rising. I've looked at the speed vectors in more detail. The problem is vanishing when the speed vectors near the zero degree line are horizontal or show away from it. I have Neumann boundary condition there ( also at 180°). ($\frac{\partial C_A}{\partial \theta}=0 $). Now I also have seen that for the regions where this "pseudo-source" exist the boundary condition isn't met perfectly. There is a gradient. So I am thinking maybe this acts like a "pseudo-source ". I also have drichilet boundary condition for r=1 and r-->infinite. ( My calculation begins at r=1).

My grid number is the following ( for example I have three radial and five angular steps ). My notation is c(i,j), i is the radial direction and j the angular direction. So say for i=1, it is c(1,j) with j=0...4 from phi=0...180°. Let's have a look for zero degress with constants l_i :
$c10*l1+...c12*l2+c14*l4=0$.

I need no ghoistpoints here. I just sad c12=c14 and rearranged the equation to $c10*l1...+c12(l2+l4)$.
For 180 ° I need ghoistpoints because of second order upwind scheme I have c(i,j-2) which is outside my domain, but at 180 ° there is no problem.

Should I implement Neumann on another way ? Or do you have any other suggestions ?


问:
Chinese character component-based layout engine algorithms
I am looking for a PDF describing Chinese character component-based layout engine algorithms. These should allow for recursive encodings of radical and non-radical components, and the encodings should account for the structural layout of characters based on these layouts.

Thank you.


问:
How to apply Dirichlet boundary condition for lowest order Nedelec element over tetrahedral domain?
The Dirichlet boundary condition is $n \times \vec{E} = n \times \vec{g}$. I don't know how to calculate the boundary value. Any kind of help would be appreciated.</p>


问:
Simulating a Hill-Type muscle model
I'm trying to do biologically based locomotion animation.

I have a Hill-Type Muscle Model (see reference 1) that I've implemented in code. I can find the initial lengths of the tendons, then do a step and find the new velocities and forces.

My question is, how do I use this information to simulate the muscle?

Specifically, we will call our muscle a Muscle Tendon Complex (MTC) that has four parts:

The contractile element (CE) which models the active force production

The parallel elastic element (PEE) which is in parallel with the CE

The serial elastic element (SEE) which is in series to the CE, and

The serial damping element (SDE) which is in parallel to the SEE

Each of these have a length $l$, where

$$l_{SEE}=l_{SDE}, l_{PEE}=l_{CE} \text{ and } l_{MTC}=l_{SEE}+l_{CE}$$

and each has a velocity $i$. We also know a relationship between the force $F$ that each tendon produces:

$$F_{CE} + F_{PEE} = F_{SEE} + F_{SDE}$$

And the output force $F_{MTC}$ of the muscle is

$$F_{MTC}=F_{SEE}+F_{SDE}$$

Then the muscle has some activation value $0.001 \leq q \leq 1$.

What I know so far is that you do the following:

Start the muscle not moving ($i_{MTC}=0$ and $i_{CE}=0$), then choose a reasonable starting length for your tendons. They use $l_MTC=0.264$ so I'll use that too. Then run a solver to find the $l_{CE}$ such that the force generated by the muscle while not moving is 0 to find the initial $l_{CE}$.

Given the current velocities and lengths, calculate the force generated by each element and a new velocity $i_{CE}$. if $i_{CE}>0$, recalculate it using other equations. This part I can do.

Calculate $F_{MTC}=F_{SEE}+F_{SDE}$ which is the output of our muscle tendon complex model.

"Integrate $i_{CE}$ to get $l_{CE}$ along with the ..."

"Multi-body simulation that calculates the generated motion from $F_{MTC}$, which results in new $l_{MTC}$, $i_{MTC}$, and $q$. Now repeat from step 2".

4 and 5 are direct quotes from the paper (including the "..."), but I'm not sure what they mean or what I'm supposed to do. Is there something trivial here I am missing?

References

Haeufle, D. F. B., et al. "Hill-type muscle model with serial damping and eccentric force–velocity relation." Journal of biomechanics 47.6 (2014): 1531-1536.

Title: Simulating a Hill-Type muscle model


问:
Generalization error and Sample Complexity estimation for Least Squares
I am wondering how to draw a sample complexity plot similar to the following figure which shows the estimated number of samples to incur no more than 10 percent generalisation error on average for the least squares algorithm for classification.

The problem assumes that we have $m$ patterns $\mathbf{x}_1,...,\mathbf{x}_m$ sampled uniformly at random from $\{−1,1\}^n$, where each label is defined as $y_i:=x_{i,1}$ (s.t. the label of pattern $\mathbf{x}$ is just its first coordinate). And the goal is to estimate the sample complexity as a function of the dimension ($n$) of the dataset.

The m (examples) versus n (dimension of dataset) plot also includes error bars indicating the standard deviation for the estimates of m. I'm not interested in that, just the function estimation.


问:
Strange convergence behavior of WENO5 for Hamilton--Jacobi equations
I have the following question. I have a code function that computes right-biased and left-biased approximations of the derivative of a function using WENO5 for Hamilton--Jacobi equations as described in

```<code>@Article{Jiang2000,
  Title                    = {Weighted {ENO} Schemes for {Hamilton--Jacobi} Equations},
  Author                   = {Jiang, Guang-Shan and Peng, Danping},
  Journal                  = SIAM Journal of Scientific Computing,
  Year                     = {2000},
  Month                    = jan,
  Number                   = {6},
  Pages                    = {2126--2143},
  Volume                   = {21},
  DOI                      = {10.1137/S106482759732455X},
}
</code>```

I refer to this WENO as WENO5JP below.

I unit-test my implementation of WENO5JP by checking the order of convergence while approximating derivatives of simple smooth functions.
I obtain right- and left-biased approximations of the derivative, average them, and compare this result to the analytic derivative. Then I measure error in $L_{\infty}$-norm, and then I compute the observed order of accuracy.

However, I get strange behavior when for one function I get between 5 and 6.5 order of accuracy, while for another function I get between 4 and 6 order of accuracy.

Computation of the observed order of accuracy. I use the following standard formula to estimate the order of accuracy:

$$
p_k = \frac{\ln \left( \|E_{k-1}\| / \|E_k\| \right)}
  {\ln \left(N_k / N_{k-1} \right)},
$$
where $p_k$ is the observed order for the grid size $N_k$, and here I use $N_k / N_{k-1} = 2$.

Details. Below are the details of convergence behavior that I get.

I apply WENO5JP to a function $y = \sin x$ on domain $[0; 2\pi]$ and obtain the order of accuracy as shown in the table below

```<code>|------------|-------|
| Resolution | Order |
|         11 |   nan |
|         41 |  6.68 |
|         81 |  6.24 |
|        161 |  6.12 |
|        321 |  6.15 |
|        641 |  5.12 |
|       1281 |  0.60 |
|       2561 | -1.36 |
|------------|-------|
</code>```

Linear least squares applied to the function $E = p \ln N + \ln C$ gives estimate of the order $p$ as $p = 6.06$. And this is the plot in a log-log scale "Error vs grid size":

I apply WENO5JP to the function $ \exp x \, \sin(5x)$ from Trefethen's book on spectral methods (p. 56) and obtain the following behavior

```<code>|------------|-------|
| Resolution | Order |
|         11 |   nan |
|         21 |  4.16 |
|         41 |  5.64 |
|         81 |  4.90 |
|        161 |  4.47 |
|        321 |  4.11 |
|        641 |  4.21 |
|       1281 |  4.45 |
|       2561 |  5.92 |
|       5121 | -0.84 |
|------------|-------|
</code>```

Linear least squares applied to the function $E = p \ln N + \ln C$ gives estimate of the order $p$ as $p =4.58$. And this is the plot in a log-log scale "Error vs grid size":

As you can see, WENO5JP exhibits some strange behavior, where for one example it gives more than the 5th order of accuracy, while in another example the order of accuracy drops to 4. Last rows of each table show that the rounding error dominates over the truncation error of the method at these grid step sizes.

My question is: is it normal for WENO5JP to have that volatile order of accuracy?

As far as I can say from the literature on WENO for conservation laws, people usually don't get strictly 5th order even for simple linear problems.

However, it is still not clear to me whether I can say that my implementation is correct or not due to this strange convergence behavior.

UPDATE: you can find the code here: https://github.com/dmitry-kabanov/weno5jp-issue


问:
What does it take to prove that a multigrid algorithm scales linearly with system size?
It is my understanding that the multigrid solution techniques are generally the preferable method to solve large Poisson problems. Now assume I have written a multigrid solver that is tailored to my special needs that actually works quite well, i.e. the residual norm decreases by a factor ~10 with every successive V-cycle.

How would I actually show/prove/test for textbook efficiency? Can optimal scalability only be shown for a fixed continuous problem that is discretized onto increasingly finer grids? If so, are there any special properties such a problem should have?


问:
Preconditioners and discrete Lagrange multipliers
There is a huge literature on efficient preconditioners for saddle-point problems. In computational physics, the case where the Lagrange multipliers enjoy a weak formulation (say, the Stokes equations) is almost always considered. In such cases, very fast solvers can be designed.

Unfortunately, in structural mechanics, the Lagrange multipliers are often introduced on the discrete level. Consider for instance multi-points constraints between degrees of freedom (this feature is very common in industrial software). I know this is not the good way of doing things, but this is how it is often done.

I know about the ancient work of St-Georges, Notay, Warzée, or the recent work of Rees and Greif, Gould, Hribar and Nocedal.
But I wonder if someone has experienced other strategies (such as the ones available in the FieldSplit preconditioners of PETSc for instance) for preconditioning saddle point systems with those horrible discrete Lagrange multipliers? I am indded seeking a way of using efficient preconditioners (say multigrid) on these saddle-point systems.

Finally, I wonder if the best way to use efficient preconditioners with Lagrange multipliers wouldn't be the domain decomposition methods, such as FETI or BDD, that can naturally handle multipliers.


问:
What is a good algorithm to solve a discrete continuity equation in Cylindrical coordinates?
The equation is:

$\partial f/\partial t + \nabla \cdot (v f) = 0$    $, \;\; f \in [0,1] $

and $v$ is a velocity known at every grid cell. A more precise constraint is that $f$ is either 0 or 1 but I can also live with it being in the domain [0,1]. The geometry is 2D cylindrical R-Z coordinates.

For the Initial condition, $f(R, Z)$ is known for all zones in the $R-Z$ domain and is not a simple function. Note that $v$ is evolving in time and space i.e. $v$ is really $v(R, Z, t)$. I do not want the solution to be diffusive.


问:
Finding eigenvectors and eigenvalues of large matrices in Python's numpy
I am running a PCA analysis on a data set using Python's (v2.7.10) NumPy.  I validated that my program works by running the PCA analysis on a smaller dataset and then confirming that I get similar results using sklearn.decomposition.

When I run my full dataset (with 60000 variables and 93 samples), I experience some strangeness.  After I construct my covariance matrix (which is 60000 x 60000), I compute the eigenvalues and eigenvectors using numpy.linalg.eig().  When I inspect the eigenvalues and eigenvectors, all the entries are exactly 0.

This leads me to believe that there is something strange with the behavior of numpy.linalg.eig() due to the large matrix size.

QUESTION : What is going on?  What assumptions are breaking down? How do I fix it?


问:
How to obtain projections from sinogram in ART reconstruction technique?
I'm kind new in the Computed Tomography field and I'm trying to understand and implement ART technique. Said it, I started to read the book The Mathematics of Medical Imaging - A Beginners Guide by Timothy G. Feeman. In the chapter 9, author describes the ART method. Well, If I understood fine (correct me If I'm wrong, please), in ART we are just trying to solve this system: Ax=p. One approach to solve it is to use Kaczmarz method (I read it and it is ok to understand it).

Now, let's understand each variable in the equation above: x is the image matrix that I want to found. A is the projection matrix and p is the sinogram. Is this reasoning right? If yes, how I obtain A from p? Can anyone show me an example with some numbers? (Ex: let A = [...], p = [....], etc.) A very simple example for the sake of simplicity would be fine.

Thanks in advance, guys!


问:
Are linear programming algorithms faster than quadratic programming algorithms?
I have an objective function that I can write either in quadratic programming (QP) such as
$$\sum_{i=1}^N \sum_{j=1}^N C_{ij}^2$$
or as an LP problem
$$\sum_{i=1}^N \sum_{j=1}^N |C_{ij}|$$
which can be done by splitting $C_{ij}$ into two parts $C_{ij} = C_{ij}^+ + C_{ij}^-$ and imposing $C_{ij}^{+}>0$ and $C_{ij}^{-}>0$. This doubles the number of optimisation variables.

Both approaches will have additional linear constraints which will be the same.

My question is the following: Is it possible to argue that the LP algorithm will always be faster than the QP algorithm given that we use the state of the art algorithms for each ?

Update: I realise that the solutions will not be the same. But assuming that I want to compare both approaches, I simply want to know which one will in general be faster. In practice I find LP to be faster. But I would like to know if this is a hard rule or does it simply come down to which algorithm I am using and the nature of the problem and constraints.


问:
Open Source Quadratic Programming with Piecewise Linear Objective
I am looking for an open source solver to solve the a quadratic programming problem with an additional piecewise linear objective, as show below. The problem is fairly small ($\mathbf{x}$ is a vector of dimension 120).

\begin{equation*}
 \max_{x}: \mathbf{x}^{T} \mathbf{q_0}
 - \mathbf{x^{T}} P_0 \mathbf{x} - f(\mathbf{x})
\end{equation*}

Note that above $P_0$ is positive semidefinite and $f$ is piecewise linear for each variable in $\mathbf{x}$ but not convex. For each $x_i$, $i \in \left\{1, ..., 120\right\}$, $f(x)$ would be a vector of coordinate points, e.g.

```<code>[(x_0, y_0), (x_1, y_1), ..., (x_N, y_N)]
</code>```

For any given $x_i$, $f(x_i)$ looks something like the following.

Since $f$ does not depend on any cross terms in $\mathbf{x}$, the above optimization could also be written as

\begin{equation*}
 \max_{x}: \mathbf{x}^{T} \mathbf{q_0}
 - \mathbf{x^{T}} P_0 \mathbf{x} - (f(x_1) + ... + f(x_{120}))
\end{equation*}

The are also linear constraints of the form $A\mathbf{x} < \mathbf{b}$ are added.

I have formulated the problem in Gurobi which is fairly straightforward but was hoping to compare this to an open source solver. Looking around it seems like GLPK would work for this, but I have no experience using this library so was hoping for a bit of info related to this or alternative solutions.

A sample problem for gurobi using the matlab API is shown below.

```<code>N = 2;

q_0 = [-3.31e3, -5.07e3];
P_0 = [-0.90e-04  -0.63e-04;
       -0.63e-04  -0.90e-04];

x = [-2.0000e9, -1.0000e9, -0.8000e9, -0.6000e9, -0.4000e9, -0.2500e9,...
    -0.1000e9, -0.0500e9, -0.0250e9, -0.0100e9, -0.0010e9, 0, 0.0010e9,...
    0.0100e9, 0.0250e9, 0.0500e9, 0.1000e9, 0.2500e9, 0.4000e9,...
    0.6000e9, 0.8000e9, 1.0000e9, 2.0000e9];

y = [2.938964e6, 0.753364e6, 0.488104e6, 0.280144e6, 0.129472e6,...
    0.053766e6, 0.011007e6, 0.003751e6, 0.001435e6, 0.000465e6,...
    0.000037e6, 0, 0.000037e6, 0.000465e6, 0.001435e6, 0.003751e6,...
    0.011007e6, 0.053766e6, 0.129472e6, 0.280144e6, 0.488104e6,...
    0.753364e6, 2.938964e6];

params.outputflag = 0;

%formulate model
model.obj = zeros(N, 1);
model.Q = sparse(P_0);
model.A = sparse(zeros(1, N));
model.sense = '=';
model.rhs = 0;
model.ub = repmat(1e9, N, 1);
model.lb = repmat(-1e9, N, 1);
model.modelsense = 'max';

% piecewise tcosts
for i = 1:2
    model.pwlobj(i).var = i;
    model.pwlobj(i).x   = x;
    model.pwlobj(i).y   = x * q_0(i) - y;
end

result = gurobi(model, params);
</code>```


问:
"Solution path" for quadratic program as regularizer changes
I am solving a quadratic program with regularization parameter $\alpha\geq0$ to get the solution to a problem of the form
$$
p(\alpha):=
\arg\min_{p\in\mathbb{R}^n}\ [\alpha(v^\top p)+f(p)],
$$
where $f(\cdot)$ is a strictly convex, continuous, and piecewise quadratic function of $p\in\mathbb{R}^n$.

I'm hoping to say something about the behavior of $p$ and/or the optimal objective value as $\alpha$ changes, in the style of solution path papers like this one.  Sadly, the paper I linked does not seem to apply since the function next to $\alpha$ is not bounded below.

Does anyone have pointers for theorems that characterize the "solution path" ($p$ as a function of the regularizer $\alpha$) of problems in this form?


问:
The Schrodinger equation for time-dependent Hamiltonian after one timestep, taking exponential or use ode solver?
The Schrodinger equation for time-dependent Hamiltonian is

$$i\hbar\frac{d}{dt}\psi(t) = H(t)\psi(t) \, .$$

Assuming I knew $\psi(t)$, I want to know $\psi(t+\Delta t)$.

Should I take exponential or use ode solver?

I know that the mathematical solution after one timestep should be (or not?)

$$\psi(t+\Delta t) = e^{-\frac{i}{\hbar}H(t)\Delta t}\psi(t) \, .$$

It gives a different result from using ode solver unfortunately. Which is the most correct one?


问:
Calculating fixed point inside limit cycle
So I'm working with a rather complicated dynamical system. Instead of writing it all out. It's probably easier if you just clone my git repository.

```<code>git clone https://gitlab.com/mdornfe1/vortex_acoustic
</code>```

The dynamics of the system is encoded in the function flow inside vortex_acoustic.py. If I integrate the system from rest I find oscillations for certain values of the parameters. The system is 13 dimensional, but here is a plot of the first variable vs. time

So it seems like the system is approaching a limit cycle attractor. From my understanding of dynamical systems there should be an unstable fixed point inside that limit cycle. I want to try to calculate it so I can analyze its stability. My idea is that the fixed point should be close to the mean of the oscillating quantities, so I use that as a starting value for a numerical root solver. Here's the code for finding the fixed point, but you can also just run stack_exchange.py.

```<code>import pickle, numpy as np, matplotlib.pyplot as plt
import vortex_acoustic as va
from scipy import optimize
from constants import *

T = 5000
dt = 0.05
tn_transients = int(0.9 * T / dt)

def calc_fixed_point(seed, params):
    return optimize.fsolve(va.flow_star, seed, params)

def calc_closest_fixed_point(y, params):
    seed = np.mean(y[tn_transients:, :], 0)
    seed[Nq+1:2*Nq+1] = 0
    y_star = calc_fixed_point(seed, params)

    return y_star, seed

y = np.load('simulation.npy')
with open('params.p', 'r') as f:
    params = pickle.load(f)

plt.plot(y[:,0])

y_star, seed = calc_closest_fixed_point(y, params)
</code>```

The function calc_closest_fixed_point takes the output of the simulation and the params passed to the simulation. It calculates the mean of the fluctuating quantities after the transients has died out. It sets the derivatives to 0 and uses that quantity seed as a starting point for fsolve. The function calc_fixed_point then finds the fixed point of the system y_star. I thought this approach would show that seed and y_star are very close to each other. But for this example I'm seeing they're quite different. Furthermore the actual value isn't inside the limit cycle.

```<code>print(seed)
>>>array([ 0.18307736, -0.11207637,  0.00382286,  0.        ,  0.        ,
       -0.02016129,  0.08881281, -0.07852504,  0.07439292,  0.0109214 ,
       -0.00894938, -0.02896082,  0.00833759])
print(y_star)
>>>array([  1.06466138e+00,  -6.65605305e-01,   2.15920825e-02,
        -5.18569988e-25,   1.44087666e-24,   7.55288134e-03,
         5.17875853e-01,  -4.41152924e-01,   4.72080951e-01,
         4.86813549e-02,  -6.04678149e-02,  -1.42533599e-01,
         5.90028044e-02])
</code>```

Also seed is pretty obviously not a fixed point. If you calculate the function flow at seed the derivatives are pretty far from 0.

```<code>va.flow(seed,0,*params)
>>>array([-0.00543102,  0.        ,  0.        , -0.00113827, -0.00042674,
       -0.00973662, -0.00775962,  0.01449168,  0.01369892,  0.00331889,
        0.00502761, -0.00026469, -0.00076228])
</code>```

Anyone know what's going? I'm correct in assuming there is a fixed point inside the limit cycle right?


问:
Decreasing - increasing - stabilising $l_{2}$ norm
Let $\bar{x}$ denote the analytical solution of a PDE. Let $x^{(k)}$ be the solution at the $k^{th}$ iteration. The initial guess for the solution is $x^{(0)} = 0$. Let $r_{0} = ||\bar{x}-x^{(0)}||_{2}$ i.e. the $l_{2}$ norm of the error. Let $r_{k} = ||\bar{x}-x^{(k)}||_{2}$ be the $l_{2}$ norm at the $k^{th}$ iteration.

Question: Is there a possibility that $\frac{r_{k}}{r_{0}}$ first decreases then increases and then stabilises to a particular value ?

Additional info: When I examine this stable value of $\frac{r_{k}}{r_{0}}$ for increasing mesh sizes, it does give me $O(h^{2})$ convergence for a uniform cell-centered grid with Dirichlet boundaries solved using Finite Difference Method (using Jacobi iterations).


问:
Solving HJ using Hopf-Lax formula (Burgers case)
Hi I want to approximate $u(x,t)$ solution of:

$\begin{cases}u_t+H(u_x)=0 \\u(x,0)=u_0(x)= e^{−(x−3)^2}
\end{cases}$,  where

$H(u_x)=\frac{u_x^2}{2}$ and my spatial domain is $\Omega=[0;10]$.

Now i want to apply Hopf-Lax formula (Evans' chapter 3 theorem 4), who says

$u(x_j,t_n+\Delta t)=\min_{a\in \mathbb{R}} \left\{ \Delta t H^*(a)+u(x_j-a\Delta t,t_n)\right\} $

in my case $H^*(p)=H(p)=\frac{p^2}{2}$ and $u(x_j-a\Delta t,t_n)$ needs to be interpolated but... which different values I can assign to $a$?

I tested it with $a=\left\{ \min_j u(x_j,t_n),\max_j u(x_j,t_n)\right\}$ but it seems only partially working.

Any help?


问:
C1 continuous spline on regular 2D-grid with quadratic 1D cuts
I want some scalar spline function defined on regular 2D grid $F(x,y)$ with continuous first derivative which is easy to intersect with arbitrary ray/line ${\vec l}(t) = (c_x t,c_y t,c_z t)$.

Finding the intersection means finding roots of equation $f_{c_x,c_y}(t)-c_zt=0$ where  $f_{c_x,c_y}(t)=F(c_x t, c_y t )$ are values of the spline along the line ${\vec l}(t)$.

Since it is easy to find roots of quadratic polynominal I want $f_{c_x,c_y}(t)$ to be piecewise quadratic polynominal for any ${c_x,c_y}$.

What is the problem / what is not a solution:

consider bi-quadratic B-spline on rectangular grid created by tensor-product of 1D quadratic B-splines. The result is bi-quadratic which means it is $f_{c_x,c_y}(t)$ is 4-th order polynominal.

function composed of Quadratic Bezier-triangles. While it is certainly just quadratic function along any direction, it is hard ensure continuous derivatives at the boundary between triangles.


问:
Numerical solution of Dirac equation (eigenvalue problem)
Suppose we have equation of the form:

$$H \Psi = E \Psi $$

where $H$ is Dirac Hamiltonian (also my question can be answered by people who are not familiar with Dirac Hamiltonian but familiar with numerical solution of different eigenvalue problems):

$$ H = \vec{\alpha} \vec{\nabla} + V(r) = H_0 + V(r)$$

and $\vec{\alpha}$ and $\beta$ are $4 \times 4$ Dirac matrices:

$$\vec{\alpha} = \left( \begin{matrix} 0  & \vec{\sigma} \\ -\vec{\sigma} & 0 \end{matrix} \right) \ , \ \ \beta = \left( \begin{matrix} 1 & 0 \\ 0 & -1 \end{matrix} \right)$$

$\vec{\sigma}$ are $2 \times 2$ Pauli matrices.

Here $H_0$ is free Dirac Hamiltonian (without potential, only kinetic term is preserved).

I think that the best way to start solving such an equation numerically is to introduce plane-wave basis (which diagonalizes free Hamiltonian $H_0$):

$$\Psi = \sum_n c_n \psi_n $$

Then we have a problem like:

$$H \sum_n c_n \psi_n = E \sum_n c_n \psi_n $$

My question is: what are the most commonly used ways to solve this equation numerically for arbitrary potential? Just to compute matrix elements $H_{nm}$ in this plane-wave basis, to find eigenvalue $E$ and coefficients $C_n$ from this computation? Or maybe there are some better methods?


问:
Irregular data access on CUDA
I am trying to parallelize a program with CUDA that has irregular memory access. To explain what I am talking about, I have attached two figures. In the first figure, I have plotted the indices of the output and input arrays that needs to be accessed during each loop iteration. Obviously I have plotted only a few iterations (the total number is actually in millions). For the loop iterations in between the black horizontal lines, I assign to the same element in the output array. It can be seen that the number of times this happens for different elements in the output array is varying.

In the next picture which is a zoom in of the first one, it can be seen that I need to access four elements of the input array for each iteration. These can be quite close together or quite far apart. I have only plotted the case when the number of elements in the input array are in the hundreds. This can be in tens of thousands.

As a first naiive approach for parallelising this, I can read all the elements in the input array into the shared memory. Then I can distribute the iterations over the threads. This way I only need to to reads/writes from the global memory and if I am careful in how I assign the work, these can be coalesced too.

Now, the problem comes when the input array don't fit into the shared memory. As you can see, for some of the iterations, I need to deal with elements quite far apart.

So, my question is this, for general, large sized (can't fit all into the shared memory) problems, what sort of algorithms or techniques can I apply to make this efficient? I have read around a bit regarding irregular programs but didn't find much that is directly applicable to my problem.

P.S. I already have this working very well with MPI. Doing this on CUDA is a bit of a curiosity thing.

Title: Irregular data access on CUDA


问:
Vectorize function integration
I need to implement the following in python:

For a given discrete time series $Z_t$ ($t={0...T}$), find the smallest $t$ such that:

$$c\sum_{s=0}^t e^{[k(Z_t-Z_s)+m(t-s)]}  \geq \frac{p^*}{1-p^*} $$

where $c,k,m$ are constants and

$$P^*=\int_0^{0.5} \frac{(1-2y)e^{-d/y}}{(1-y)^{1+d}y^{1-d}}dy=\int_{0.5}^{p^*} \frac{(2y-1)e^{-d/y}}{(1-y)^{1+d}y^{1-d}}dy $$

(d is another constant)

I need to implement this for each row of an $N\times w$ array $Z$.

I implemented the sum in equation 1 above utilizing the fact

$$ e^{[k(Zt-Zs)+m(t-s)]} = \frac{e^{kZt+mt}}{e^{kZs+ms}} $$

hence

$$\sum_{s=0}^t e^{[k(Z_t-Z_s)+m(t-s)]} =e^{kZ_t+m_t}*\sum_{s=0}^t\frac{1}{e^{kZ_s+m_s}} $$ which is reflected in use of cumsum in code below

```<code>def f(z):
    return ((1-2*z)*np.exp(-d/z))/(((1-z)**(1+d))*(z**(1-d)))
lhs=integrate.quad(f,0,0.5)

def rhs(p):
    return integrate.quad(-f,0.5, p)
p_star= fsolve(rhs-lhs,0.75) # will depend on time series only indrectly when d will be optimized

for i in np.arange(N):
    z=Z[i,:]
    main = np.exp (kz+m*np.arange(w))
    cumsum_t=np.cumsum(1/main)
    final_sum= main*cumsum_t
    t_solution= # index i where final_sum[i]> p_star/(1-p_star)) # not implemented yet
</code>```

Is there any way I can vectorize this for $Z$? In this case $N$ is ~400,000 so vectorization would really help. f(z) will be fine with vector inputs but I don't think rhs(p) will be fine as it uses integration.

Title: Vectorize function integration


问:
Does applying the Newton-Raphson iteration for matrix reciprocal refine a matrix inverse from LU/GE?
This is a follow-up to this answer.

Suppose you have a possibly very ill-conditioned matrix $A$, and you compute its inverse with LU/GE to get $X_{\text{lu}}\approx A^{-1}$.

The Newton-Raphson iteration
$$ f: X \mapsto 2X-XAX $$
generally converges to the matrix inverse of $A$.

Does computing applying the Newton-Raphson iteration to $X_{\text{lu}}$, that is, computing $f(X_{\text{lu}})$, $f(f(X_{\text{lu}}))$ make the resulting matrix inverse more accurate, or less accurate?

I was unable to find a reference directly, so I'm asking if someone knows for sure.

The absolute error will be proportional to $\|(f')^{-1}\|$ (e.g., Accuracy and Stability of Numerical Algorithms, Chapter 25), but $f'$ vanishes at $X=A^{-1}$. So presumably this means that applying $f$ could only help when the errors came from the growth factors in the LU factorization step? But doesn't that mean that applying $f$ hurts at other times since it relies on computing a root of a function with a vanishing Jacobian, which would make it less accurate than LU? For example, in the one-dimensional case with a root of multiplicity 2, Newton's method would have error $O(\sqrt{\epsilon})$, but LU/GE always has error $O(\rho \epsilon)$.

I tried it on all $16\times 16$ matrices in MatrixDepot.jl, and it seems like it almost always hurts accuracy, never helps by more than $\frac13$, and often fails. (It's not very clear what matrices to test.)

```<code>module MatrixInverseNewtonRaphson

using MatrixDepot

F(A, X) = 2X - X*A*X
R(A, X) = norm(A*X - eye(A), 1)

function go(name, m)
    A = full(matrixdepot(name, m))
    B = inv(A)
    B1 = F(A, B)
    X = 0.5A' / norm(A*A', 2)
    for i in 1:2m; X = F(A, X); end
    (R(A, B1) / max(eps(eltype(A)), R(A, B)), R(A, B), R(A, X), name)
end

function main(m=16)
    e = []
    for name in matrixdepot("all")
        try
            push!(e, go(name, m))
        catch err
            println(err)
        end
    end
    sort(e)
end

end
</code>```


问:
Lambdify error when sympy array contains number
This example is a dummy example, but it describes exactly whats goes on in my problem.

As a result from some sympy calculation and manipulation i got two arrays

```<code>array_1 = ([f_1(x,y)      0         f_3(x,y) ])
array_2 = ([g_1(x,y)   g_2(x,y)     g_3(x,y) ])
</code>```

where f_n(x,y) and g_n(x,y) are complex functions. When I print them to se what kind of creature they are, they both are on the form

```<code>Matrix([[#, #, #]])
</code>```

as expected from a sympy matrix.

Furthermore I use lambdify to convert the arrays by

```<code>lambda_array_1 = lambdify([x, y], array_1, 'numpy' )
lambda_array_2 = lambdify([x, y], array_2, 'numpy' )
</code>```

Then I let

```<code>x = 1
y = 1
x_array = np.array([x], dtype = complex)
y_array = np.array([y], dtype = complex)

pythonarray_1 = lambda_array_1(x_array , y_array )
pythonarray_1 = np.reshape(pythonarray_1, 3)

pythonarray_2 = lambda_array_2(x_array , y_array )
pythonarray_2 = np.reshape(pythonarray_2, 3)
</code>```

When I then print the two pythonarrays I get

```<code>pythonarray_1 = array([array([# + #j]), array([# + #j]), array([# + #j])], dtype=object)

pythonarray_2 = array([# + #j, # + #j, # + #j ])
</code>```

where I need the arrays on the form as in pythonarray_2.

What is wrong with pythonarray_1, how can I fix it, and why does it happen?


问:
Porting from MPI to GPU
Assuming that the underlying algorithm can be ported to multi GPU, what aspects should one consider while porting from MPI (on multiple nodes) to multi GPU (again on multiple nodes)?

Making use of Hyper Q ? Of MPS ? Peer to peer communication ?

For example, if we consider an MPI code that solves for the Navier Stokes equations using ADI, on multiple nodes with domain decomposition, and solves for the pressure and velocity by solving the schur complement matrices and interface unknowns induced by the domain decomposition. We could solve for internal unknowns by a simple tridiagonal solver (say, Thomas algorithm) which is highly parallel and use the schur complement method to solve for the interface unknowns and repeat this throughout the time stepping.

The schur complement would involve solving directional systems (MPI ranks in depending on the cartesian decomposition in one direction) and hence would involve communication between the ranks in one direction.

The issue of communication arises when porting to the GPU, as we would need to communicate between the GPU's on different nodes which are logically in the same direction (considering 1:1 ratio of MPI ranks: GPU).

Any suggestions or ideas are welcome.

Thank You.

Title: Porting from MPI to GPU


问:
Precision of ratio constraint in a linear program
I am trying to solve a LP in which one of my constraints is of the form

$$\frac{A(x)}{B(x)} = 1$$

I transform this constraint into a linear one

$$A(x) - B(x) = 0$$

However when CPLEX solves the problem, this constraint is of the order of ~1.e-5, but the ratio is of the order of ~1.4, which is a large deviation for my needs. Unfortunately the kappa value (condition number of the optimal basis) of the problem is kappa ~10^10, and this doesn't make things easier. I suspect this is the source of the problem.

My question: is it possible (given the value of kappa) to increase the precision of this constraint such as the ratio to be accurate within 2 decimal digits (at least) ~ 1.00something ?


问:
vectorizing optimization or root finding
I need to find the roots of a function. I am currently using scipy.optimize.fsolve

```<code>    def f(z):
        return ((1-2*z)*np.exp(-d/z))/(((1-z)**(2+d))*(z**(2-d)))
        # d is a constant ndarray of shape(300000,)
        # hence f(z) is an ndarray of shape (300000,)

    def trap(func,a,b,num): #trapezoidal approach to integration
        xlinear=np.linspace(a,b,num)
        dx=np.diff(xlinear)
        fx=func(xlinear[:,None])
        f12=fx[:-1,:]+fx[1:,:]
        return np.dot(dx,f12)/2.0
        #also returns an array of shape (300000,)

    def integral(p):
        return trap(f,0,p,1000)

    root= fsolve(f,0.75) # 0.75 because the answer should lie in (0.5,1)
</code>```

The issue is that I was expecting fsolve to return an array in line with returns of functions f and trap. However, somehow I am only getting a array of shape (1,)


问:
Sensitivity Analysis -- Total Variation for a function with categorical arguments?
I have an application in sensitivity analysis of complex system models with moderately nonlinear interactions between arguments

Arguments potentially dozens or hundreds in number

Arguments mostly real-valued but some categorical with a handful of levels (unordered settings) per variable

The outputs of interest are real-valued and highly interdependent, and respond well to the usual dimensionality reduction techniques

Assume for present purposes that the model output is a single real number.

I have the compute resource to evaluate the model several hundred or a few thousand times, but not millions of times, and I will typically do this on random or Sobol'-gridded points in the design space, i.e. I will have the output value at a few thousand points randomly scattered in a few dozen dimensions.

I want to do the usual thing of assigning importance rankings to inputs and combinations of inputs using techniques from Saltelli et al -- see this paper to get a flavor. I want to estimate the total variation of the expectation of the model output conditional on one or a few inputs at a time.

I plan to do this by fitting a low-order smoothing spline through the space of one input, or a few possibly interacting inputs, plus the model output, and calculating the total variation of the smoother by numerical integration. The interacting sets will typically be three or four variables, so sparsity in high dimensions will never be an issue. I also know how to cross-validate and regularise the splines for optimal fit, and correct for randomness in the total-variation measure.

What I don't know is how to proceed when one or more of the inputs of interest is categorical. I can see what to do for a single binary-valued input, but anything more complex defeats me. I've wiki'ed and googled this up and down without success.

So to summarise the question:

Is there an extension of the concept of total variation to scalar real-valued functions with (possibly multiple) (more than binary) unordered categorical arguments?

and

Is there an explanation digestible by a mere engineer?

Thanks in advance, I hope, but I'm actively looking in multiple places and I'll update this post if I find anything useable elsewhere.


问:
Algorithm for face based data-structure - CFD
Good morning

I'm trying to develop an unstructured CFD code to solve Euler equations in a finite-volume (cell-centered) context (learning purposes). I was able to build from a cgns file some basic information for the code. What I have now are:

Node index and its 'x' 'y' coordinates.

Element index and nodes (indexes) that build these elements.

Boundary conditions.

Since the code is 2D every element is actually a cell and every face is actually a line forming this cell.

The basic data-structures I build right now are:
+ Node array where I put the node index and choose between the x and y coordinates.
+ Element/node connectivity (anti-clockwise manner) where (for Quad):

```<code>connec(1,nelem) = node 1
connec(2,nelem) = node 2
connec(3,nelem) = node 3
connec(4,nelem) = node 4
</code>```

These data-structures were very simple to build but the next one is very obscure for me. I need a face based data-structure where each face "knows" the nodes that forms it the left element index and the right element index.

The end result would be:

```<code>faces(1,face_index) = node 1 of the face.
faces(2,face_index) = node 2 of the face.
faces(3,face_index) = left element.
faces(4,face_index) = right element.
</code>```

I searched a lot but could not find a reference that explains me the algorithms to build this face data-structure. Lohner's book about Finite-elements have very nice explanations about node-based datas-tructures but not about face-based ones.

Anyone has a good reference that explains these algorithms or even the algorithm itself ?

Thank you !
cheers


问:
How to compute $\Delta f/ f$ in computational neuroscience?
I want to compute $\Delta f/f$ (as explained in this link, reference 1), but I am not able to do integration as I have distinct values and no function to generate the continuous time based values.

I have a $1\times N$ time series data where $N$ is the total number of neuron response for one neuron, one for each millisecond. How do I calculate $\Delta f/f$ for this single neuron? If I do $(f_{i+1} - f_{i})/f_{i}$, it gives me very small values.

References

Jia, Hongbo, et al. "In vivo two-photon imaging of sensory-evoked dendritic calcium signals in cortical neurons." Nature protocols 6.1 (2011): 28-35. doi:10.1038/nprot.2010.169


问:
Geometric Multigrid for Conform and Non–″ Elements: Restriction Operators
First of all, let me set up some notations.

Suppose we have a hierarchy of meshes $\mathcal{T}_0 \subset \mathcal{T}_1 \subset \dots \subset \mathcal{T}_n$. For simplicity I restrict myself here to $\mathcal{T}_H := \mathcal{T}_0 \subset \mathcal{T}_1 =: \mathcal{T}_h$:

$\mathbb{U}_H$ and $\mathbb{U}_h$ are spaces of FE–interpolants defined on the coarse and fine mesh, respectively. $\mathbb{R}_H$ and $\mathbb{R}_h$ are corresponding spaces of vectors of coefficients.

$\mathcal{P} \; : \; \mathbb{R} \longleftrightarrow \mathbb{U}$ is a natural isomorphism. For one, for linear Lagrange shapes we have $\mathcal{P}_H \; : \; \mathbb{R}_H \ni \langle \, \xi_0, \, \xi_1, \, \dots, \, \xi_{16} \, \rangle \, \mapsto \, \sum_{i = 0}^{16} \xi_i \, S^H_i(x, y) \in \mathbb{U}_H$, where $S^H_i$ is a usual “hat–function.” $S^H_{15}$ looks like this:

i. Prolongation, $\mathbb{U}_H \subset \mathbb{U}_h$

We want to define a prolongation operator $P \, : \, \mathbb{R}_H \longrightarrow \mathbb{R}_h$ in order to travel between meshes in MG–routines. It is easy in this case:

$$
P \, : \, \mathbb{R}_H \stackrel{\mathcal{P}_H}{\longrightarrow} \mathbb{U}_H \subset \mathbb{U}_h \stackrel{\mathcal{P}_h^{-1}}{\longrightarrow} \mathbb{R}_h,
$$

so for linear Lagrange elements $P$ is a sparse matrix:

You can see that coarse DOFs, $i = 0, 1, \dots, 16$, are multiplied by unity (i.e. they remain the same), and fine DOFs are just means of DOFs associated with the ribs they lie on (since the shapes are linear functions).

Let’s fix some $\xi \in \mathbb{R}_H$. I obtain the following results:

$\mathcal{P}_H \, \xi$ (green) and $\mathcal{P}_h \, P \, \xi$ (blue)

So you can see that everything works fine.

ii. Prolongation, $\mathbb{U}_H \not\subset \mathbb{U}_h$

Let’s stick to Crouzeix–Raviart elements. In this case DOFs are associated with ribs, so we should enumerate them:

As one may notice, we do not have $\mathbb{U}_H \not\subset \mathbb{U}_h$ since $u_H \in \mathbb{U}_H$ may be discontinuous at midpoints of fine ribs.
So we have to use some tricks in order to define $P$.

Namely, we take a space $\mathbb{V} := \{ \, \text{functions that are linear on each fine triangle of } \mathcal{T}_h \, \}$, $\mathbb{U}_H, \mathbb{U}_h \subset \mathbb{V}$; then we define a mapping $\phi \, : \, \mathbb{V} \longrightarrow \mathbb{U}_h$ as follows:

$$
(\phi \, v)(x) = \frac{1}{2} \left( v|_{t_1}(x) + v|_{t_2}(x) \right)
$$

for each pair of fine triangles $t_1$ and $t_2$ that share a rib ($x$ is a midpoint of the rib). Now we are able to define our operator:

$$
P \, : \, \mathbb{R}_H \stackrel{\mathcal{P}_H}{\longrightarrow} \mathbb{U}_H \subset \mathbb{V} \stackrel{\phi}{\longrightarrow} \mathbb{U}_h \stackrel{\mathcal{P}_h^{-1}}{\longrightarrow} \mathbb{R}_h.
$$

Pattern of sparse matrix $P$,
 $\mathcal{P}_H \, \xi$ (red), and $\mathcal{P}_h \, P \, \xi$ (yellow)

So you can see that everything works fine here, too.

iii. Restriction Operators: Here Comes a Problem

In the corresponding literature it is common to see that the rule of thumb of defining a restriction operator is that it should be “somewhat transpose of $P$.” And by “somewhat” they usually mean “multiply by a scalar.”

For one, in [1, p. 436] the prolongation / restriction operators are defined via interpolation / weighting (without information extraction about FE–spaces at hand). In this sense I prefer approaches of [2] and [3] more.

The question is, given $P$ (from (i) and (ii)), how can one canonically define the corresponding restriction $R$?

For the case (i), if I fix some vector $\xi \in \mathbb{R}_h$ and blindly multiply it by $P^T$, I will get the following (I used 3 mesh levels here for the sake of clarity):

$\mathcal{P}_2 \, \xi$ (green), $\mathcal{P}_1 \, P_{1 \rightarrow 2}^T \, \xi$ (blue), and $\mathcal{P}_0 \, P_{0 \rightarrow 1}^T \, P_{1 \rightarrow 2}^T \, \xi$ (red)

So you can see that such “restrictions” are not good at all. And I doubt that multiplication by a scalar will do here. For the case (ii) I have similar problems.

I bet the solution is “hidden” in definitions of the following inner products:

[3, p.160]

[2, p.110]

Please help me to clarify this so I will be able to implement restrictions correctly.

References:

Yousef Saad

Iterative Methods for Sparse Linear Systems (2nd Edition)

2003

Maxim A. Olshanskii, Eugene E. Tyrtyshnikov

Iterative Methods for Linear Systems

2014

Susanne C. Brenner, L. Ridgway Scott

The Mathematical Theory of Finite Element Methods (3rd Edition)

2008


问:
A test suite of large systems of nonlinear equations
I am looking for a kind of modern test set of large nonlinear problems. The only option I managed to find so far is rather dated: http://folk.uib.no/ssu029/Pdf_file/Testproblems/testprobRheinboldt03.pdf </p>


问:
Equation of Motion Moving Cylinder, Displacement due to spring force
I have to simulate a cylinder moving in transverse direction due to an external force (no damping). The equation of motion is: $M\frac{d^2Y}{dt^2}+kY = F_Y$

My confusion is regarding $kY$, how to get displacement when the cylinder is not at $(0,0)$ but at an arbitrary position say, $x_c=0.5, y_c=0.5$. The cylinder should oscillate freely about the ($0.5,0.5$) point.

I would like to solve this equation using Euler method, what will be the $Y$ in this scenario. My version is:

```<code>Y = Y0 + dt * V0;
V = V0 + dt * (-K * Y0 + FY) / M;
</code>```

V0 and Y0 are the values at old/previous time step.


问:
Scaling step size in adaptive runge-kutta method
I'm developing my own generic Runge-Kutta solver, and I'm currently implementing the adaptive step-size routine. I say generic because I want to be able to test different RK implementations by only passing the solver a Butcher tableau, such as the following.

\begin{array}{c|cccc}
c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\
c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\
\vdots &\vdots &\vdots &\ddots &\vdots
\\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\
\hline
&b_{1}&b_{2}&\dots &b_{s}\\
&b_{1}^*&b_{2}^*&\dots &b_{s}^*\\
\end{array}

Currently, I'm only focused on adaptive-explicit methods. As far as I understand it, a non-adaptive method would only have only the first row of $b$, and an explicit method would be zero on the diagonal and upper triangle.

Essentially what I'm doing within a single iteration of my solver is passing a timestep span [$t_0, t_1$], and initial guess of the step ($h$), and an initial value of $y_0$, which could be a vector, and the function $\frac{dy}{dt} = f(t, y)$.

I calculate $y(t_0 + h)$ and $y^*(t_0 + h)$ using coefficients of the Butcher tableau, calculate the element-wise error between them, and then calculate a scaling factor of the step using the maximum absolute error. The scaling factor is used to determine if $h$ is too small, a good size, or too big, and is adjusted accordingly. Currently, I'm calculating the scaling factor like this (described in this youtube video):

$$ s = \frac{\epsilon \cdot h}{2 \cdot (t_1 - t_0) \cdot max(|y-y^*|)}$$

where $t_0$ and $t_1$ are the beginning and end of the time step of interest, and epsilon is some error parameter set previously (e.g. 1e-8).

The resulting value of $h$ would be scaled using $s$:
$$h = sh$$

If s is larger than 2, then I assume the error between $y(t_0 + h)$ and $y^*(t_0 + h)$ is small and I can continue with the next step. If s is less than 1, then I adjust h without moving forward and try the same timestep again. If $s$ is between 1 and 2, then I move on as well, but I don't adjust h because it's good enough.

What as not readily obvious to me at the start, but now becomes so as the $tspan$ increases ($t_1-t_0$ goes up), is that the acceptable h is kept really far down. This also seems to be greatly influenced by the RK method used.

For example, the following system takes about 2 seconds on my laptop with 100 pts between [0, $4\pi$]:

\begin{array} {rl}
y_1'(t) &= y_2 \\
y_2'(t) &= -y_1 \\
y_1(0) &= 0 \\
y_2(0) &= 0.1 \\
\end{array}

Decreasing the number of sampling points to 10 within the same range ([0, $4\pi$]) now takes 16 seconds because I think the value of s is too small, due to the $t_1-t_0$ term.

It may be that I'm only interested in a very few number of points (or even the final value), but this is causing the solution time to increase, even for this relatively small problem.

Can someone point me to a reference that goes into more detail about this scaling ($s$) parameter? Or is there another way to solve this approach this problem?


问:
Eigenvalues with absolute values close to 1
I have a generalized eigen-problem: $A\psi = \lambda B \psi$ with $A$ and $B$ are large (>1000) complex non-Hermitian matrices. I know that eigenvalues with largest and smallest absolute values can be found using iterative approaches like ARPACK. But I am only interested in those complex eigenvalues (not all of them) with absolute values closest to 1. It there any fast and parallel approach to do this?</p>


问:
Implementing initial conditions into the solution domain of a 1-D advection-diffusion equation
I have the following PDE.

$\frac{\partial c}{\partial t} + v\frac{\partial c}{\partial x} - D\frac{\partial^2 c}{\partial x^2} = 0
\\$.

I have discretized it such that i now have
$\frac{dC}{dt} =  aC_{i+1} +bC_{i} +cC_{i-1} \\$    (1)

I have been given the following initial conditions:

$\\C=C_0\quad\forall  \ 0<t<\tau\\$

$\\C=0\quad\forall  \ t>\tau\\$

i have been asked how i would go about: ' describe how to implement the initial conditions in terms of the solution domain and in terms of any extra vectors'.

This question is from a degree course and is about injecting a trace of dye.

I have discretised the convective component using simple upwind and the diffusive component using central differences. I'm struggling to understand what the question is asking me.

I understand that evaluation at i=1 poses a problem in that it requires C at $\ C_0\\$. This is outside of the solution domain and would require a ghost node, but this is outside the solution domain.

What are peoples thoughts?

For the 1st initial condition i was thinking of ignoring the diffusive component and for the 2nd initial condition ignoring both the diffusive and convective terms since there is no dye left. Then i would have two equations that describe the initial conditions.


问:
1D FEM with nonhomogeneous Dirichlet BC
Consider the problem $-u''=f$, $u(a)=A$, $u(b)=B$. Everywhere, FEM seems to be introduced by using the case $(a,b,A,B) = (0,1,0,0)$, in which case the term $u'(a)v(a)-u'(b)v(b)$ resulting from integration by parts vanishes. If we divide $[a,b]$ into $a=x_0 < ... < x_n = b$, we'll end up solving an $(n-1)\times (n-1)$ system, fine and dandy.

How do you solve the above problem without convenient boundary conditions? I have tried searching and the few results I find are for higher dimensions, or I simply don't understand the answers. I have looked up several books on the matter and I can't find an answer.

Naively I would form a stiffness matrix $K=(K_{ij})_{i,j=0}^n$ and the load vector $b$ per the "usual" way, then add... something, somehow. In a case with BCs $au'(0) = \kappa_0(u(0)-g_0)$, $-au'(L) = \kappa_L(u(L)-g_L)$ the solution was to add $\kappa_0$ to $K_{00}$, $\kappa_L$ to $K_{nn}$, $\kappa_0 g_0$ to $b_0$ and $\kappa_L g_L$ to $b_n$. I don't know if that's at all applicable (I don't have any information about first derivatives) but those are the only special cases I see covered. On the other hand, $(n+1)\times (n+1)$ doesn't sound right since we know the values at the boundary and so we shouldn't want degrees of freedom there...

I really need help with this, I'd appreciate any advice, thanks.

Title: 1D FEM with nonhomogeneous Dirichlet BC