问:
Computing solutions with singularities using MATLAB ODE45
I am new to solving numerically ODES and thus it is difficult for me to judge the reliability/trustworthiness of the results that I have produced for the following problem:

I am dealing with a 2nd order ODE (I can add a non-linear term to it, too) on the interval $[0,1]$ with regular singularity at $0$ and $1$.
I converted the ODE into a 1st order one and interpreted it as a IVP by setting an initial condition at $x=1$ motivated by symmetry arguments and let ODE45 integrate from $x=1$ to $x=0$.
In this way I can produce solutions which is nice.

But here comes the problem: I have to use these solutions later on, e.g., I can also compute the $L^2$-norm of the solutions on the interval $[0,1]$ with respect to some measure and what I get as a result seems to depend sensitively on

(1.) the cut-off at $x=0$ (I let ODE45 run until "eps", the machine precision but I could also take smaller and bigger values)

(2.) the smallness of AbsTol and RelTol (both around eps)

(3.) the tspan size, i.e., the number of slices into which I cut the interval, like $10^6$.

Depending on how I vary the values in (1.)-(3.) I get results which differ from one another and not to a small extent.
I believe that the code is correctly implemented, rather I guess the problem comes from the regular singularity which leads me to the question how trustworthy ODE45 is for such a problem.
Does anybody have experiences with such a problem and can recommend what I could do to get more reliable results?
Thanks!

Update:

The equation reads $2x~(1-x)~f''(x)+(3-4x)~f'(x)+a~f(x)+b~f^n(x)=0;~~ a,b\in\mathbb{R}, n\in\mathbb{N}.$

"not to a small extent" means a change in the result for the normalization from e.g $10^{-2}$ to $10^{20}$ if I solve until $x=eps*100$ vs. $x=eps$.


问:
Finite difference aproximation - Darcy law
I am solving following problem:

Filtration of water can be described in bi-dimensional case by
$$- \partial_x(K(x,y) \partial_x u ) - \partial_y (K(x,y) \partial_y u ) = 0, $$
where $u$ - water level
$K$ - hydraulic conductivity
On uniform mesh of $N$, $K$ is defined as
\begin{align}
K(x,y) & =  \begin{cases}
    (x -0.5)^2 + (y - 0.5)^2 - 0.04 &  \geq 0\\
     &  <0
  \end{cases}\\
\end{align}

\begin{align}
- \partial_x K \partial_x u  \ \ \text{is defined as} \ \  \frac{K(i+ 1/2) \partial_x u_{i+1/2} - K(i - 1/2) \partial_x u_{i - 1/2} }{h},
\end{align}

where $(i + 1/2)$ refers to mid-points of $(i + 1)th$ interval.

I have extended discretization into 2D case and obtained following system (not sure, if its correct)
\begin{align}
\begin{bmatrix}
   c_{i}       & c_{i+1} & 0 & d_{i+1} & 0 & 0 & 0 & 0 \\
    c_{i-1}       & c_{i} & c_{i+1} & 0 & d_{i+1} & 0 & 0& 0 \\
0       & c_{i-1} & c_{i}  & c_{i+1} & 0 & d_{i+1}&  0 & 0\\
    d_{i-1} & \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \ddots \\
    0  & 0  & 0 & \dots &\dots & &  c_{i-1} & c_{i}
\end{bmatrix}
*
\begin{bmatrix}
    u_{1,1}  \\
    u_{2,1}  \\
    u_{N-1, 1}  \\
    \vdots  \\
    u_{N-1,N-1} \\
\end{bmatrix} = 0
\end{align}
where coefficients are calculated as
\begin{align}
c_i  &= + K(i -1/2,j) + K(i +1/2,j) + K(i ,j+1/2) + K(i,j-1/2)\\
c_{i + 1} &= - K(i+1/2,j) \\
c_{i- 1} &= - K(i - 1/2,j) \\
d_{i - 1} &= -  K(i,j - 1/2) \\
d_{i + 1} &= - K(i,j +1/2) \\
\end{align}

System has to satisfy following BC:
\begin{align}
\partial_x u &= 0 \ \  &\text{on} \{0,1 \} \times (0,1)\\
u &= 10 \ \  &\text{on} \{0 \} \times (0,1) \\
u &= 0 \ \  &\text{on} \{1 \} \times (0,1)
\end{align}

Can anybody give me some advices how to handle them ?

ContentLicense: CC BY-SA 3.0


问:
Spectral Collocation (or Weighted Residual) Methods to solve Stiff ODEs?
I have a system of ODEs which is (at least moderately) stiff.

Consider the class of spectral collocation methods https://en.wikipedia.org/wiki/Spectral_method or the related class of weighted residual methods.

Basically, take a system of ODEs like
$$
y'(t) = f(t, y(t))
$$
where $f(\cdot,\cdot)$ is nonlinear in my case and leads to a stiff ODE, and I want to solve it on $t \in [0,M]$.  Then approximate the solution $$y(t) \approx \sum_{i=0}^N a_i T_i(t)$$
where $T_i(\cdot)$ are a basis (e.g. Chebyshev polynomials), $N$ is the maximum order, and the basis is over $[0,M]$ if appropriate.  Note that given $\{a_i\}$, $y'(t)$ can be calculated by differentiating the polynomial basis.

Now, evaluate the ODE given a $\{a_i\}$ at the roots of the polynomial basis (which can be shown to be optimal) and find the residual.  Choose $\{a_i\}$ to minimize this residual, or just solve a system of equations with the residual = 0.

Question: how well does spectral collocation with a Chebyshev basis handle systems that seem stiff when solved with finite differences?  The worry is that the same issues that cause finite differences to have issues (e.g. convergence to a very flat function with no slope) could cause major issues with a finite dimensional approximation of the function.

If Chebyshev has problems, then is there a better basis which might have less issues.

A little more context: I have implemented spectral collocation with Chebyshev polynomials to solve a stationary Kolmogorov Forward equation in CDFs.  These, of course, need to be strictly positive and weakly monotonically increasing and asymptotically go towards $1$ if it is a proper probability distribution.  I am seeing weird behavior where the pdf drops below $0$ (i.e., it becomes decreasing) or strange jumps at the corners.  It is very hard for me to tell if these are artifacts in the solution scheme or if I have just chosen parameters for the model where no solution exists.  Increasing the number of Chebyshev basis functions doesn't eliminate them, but that could be because they can't, or only slowly, converge.

Even more context:  For what it is worth, in the simplest form of my problem there are two sets of ODEs (which are ultimately coupled):

Parameters: $\lambda_{\ell}, \lambda_h, r,$ and $\chi$.  In some ways $F_{\ell}'(0), F_h'(0),$ and $g$ are also parameters for the purposes here.  There are two discrete states which leads to the system: $i \in \{\ell,h\}$

The variable $z\in[0,\infty)$ in reality, but these converge pretty quickly so choosing some $z \in [0,\bar{z}]$ for collocation methods is reasonable.

System (1) The following comes from a stationary Kolmogorov Forward Equation (given a $\gamma(z)$ function such that $\gamma(0) = 0$ and $\lim\limits_{z\to\infty} \gamma(z) = g$.  The latter condition leads to a vanishing derivative term in the ODE (i.e., a singular mass matrix if put in that canonical form.)
$$
\begin{align}
0 &= g F_{\ell}'(z) + \lambda_h F_h(z) - \lambda_{\ell}F_{\ell}(z) + g (F_{\ell}'(0) + F_h'(0))(F_{\ell}(z) + F_h(z)) - g F_{\ell}'(0)\\
0 &= (g - \gamma(z))F_h'(z) + \lambda_{\ell}F_{\ell}(z) - \lambda_h F_h(z) - g F_h'(0)
\end{align}
$$

Forget that this is a linear ODE, as the actual one is nonlinear and much trickier.  However, it can still be written in a $M(z, F(z))\cdot F'(z) = \Phi(z, F(z))$ term for some operators and a mass matrix.

The trouble seems to come out of the term on the derivative going to $0$ (i.e., a singular mass matrix).

It might seem odd that there is a $F_{\ell}'(0)$ and a $F_h'(0)$ in the ODE as parameters but this is not a mistake.  You can easily show that any solution to the ODE will have the $F_i'(0)$ matching these constants by construction.

The initial condition is $F_i(0) = 0$.  Furthermore, from any $F_i'(0)$, one can analytically find $\lim\limits_{z\to\infty}F_h(z)$.  With $F_{\ell}(\infty) + F_h(\infty) = 1$, this means we can write it as a BVP if we wish.

System (2) The following comes from a Hamilton-Jacobi-Bellman equation

\begin{align}
0 &= 1  - (r + \lambda_{\ell})w_{\ell}(z) - g w'_{\ell}(z) + \lambda_h w_h(z)\\
0 &= 1 - (r + \lambda_h)w_h(z) - (g - \frac{\chi}{2} w_h(z))w_h'(z) + \lambda_h w_{\ell}(z) + \frac{\chi}{4} w_h(z)^2
\end{align}
subject to, $w_{\ell}(0) = w_h(0) = 0$ and defined on $z \in [0,\infty)$.  You can analytically find $w_{\ell}(\infty)$ and $w_h(\infty)$ which can help with various methods to set it as a BVP.

One can show that $\lim\limits_{z\to\infty}(g - \frac{\chi}{2} w_h(z)) = 0$ which leads to a vanishing derivative term in the ODE.

Coupled: In reality, the $\gamma(z)$ from the first solution is
$$
\gamma(z) \equiv \frac{\chi}{2}w_h(z)
$$

Leaving Out: A whole bunch of equilibrium conditions in terms of integrals of $w_i(z)$ and $F_i(z)$ which pin down $g, F'_{\ell}(0)$ and $F'_h(0)$.  Evaluating these concurrently with the solution to the ODEs is the reason that spectral methods are preferable since I can add in more equations and solve everything at the same time.

... is what I have written useful?  Maybe not, but you can see the asymptotic singularity in the mass matrix.  Is this stiff?  I think so.  When you solve the second system with finite differences, for example, it has all the hallmarks of a stiff system.


问:
Solving the transient advection-diffusion equation analytically
Can the following 1D transient advection diffusion and corresponding initial and boundary conditions:

$$
\frac{\partial c}{\partial t} + v\frac{\partial c}{\partial x} - D\frac{\partial^2 c}{\partial x^2} = f \quad\forall x\in(0,1)
\\
c(x=0,t) = 0 \\
c(x=1,t) = 0 \\
c(x,t=0) = 0
$$
be solved such that we first move the advection and diffusion operators to the RHS, integrate with respect to time, and obtain the following second order ODE:
$$
c + vt\frac{\partial c}{\partial x} - Dt\frac{\partial^2 c}{\partial x^2} = tf
$$
while assuming that $v,\;D,\;and\;f$ are constants (i.e., they do not vary spatially or temporally)? Because when I tried this and compared the solution to my SUPG and Operator Splitting formulation (assuming the same $h$-size with backward euler), it seems these numerical discretization approach steady state much quicker than the analytical solution.

For instance, if $h$-size = 1/100, $\Delta t$ = 0.01, $v$ = 1,  $D$ = 1/150, and $f$ = 1, SUPG/OS comes close to the steady-state solution at $t$ = 1.0, but for the analytical solution to be close, the value of $t$ needs to be somewhere around 10 to have roughly the same closeness to the steady-state solution.

I want to know if this is an acceptable approach if I want to perform a numerical convergence study comparing SUPG/OS with the analytical solution. I would rather compare the solutions at the exact time levels as opposed to seeing how close they get to the steady-state solution after $T$ amount of time.


问:
How can I use ODE events in MATLAB?
I need to have a better understanding about how to define ODE events.
What I know is that if I have my ODE defined as

```<code>[t,q,te,qe,ie] = ode45(odefun,tspan,q0(:,1),options);
options = odeset('Events',@events);
function [value,isterminal,direction] = events(t,q)
    value = ...;
    isterminal = ...;
    direction = ...;
end
</code>```

I need to define the value as a function of q and tin such a way that it goes to zero. So, for example, if I want to find when v = 0 I just need to write

```<code>value = q(2);
isterminal = 1;
direction = 0;
</code>```

But is writing

```<code>value = [q(2),q(2)];
isterminal = [1,1];
direction = [1,-1];
</code>```

equivalent to the previous formulation?

And what if I want to stop the integration when q(2) < toll? I mean, is it possible to stop the integration if q(2) is very close to xero but not exactly zero? (I am thinking about something like the Karnopp model used to model the sticking condition in this picture http://vibrationacoustics.asmedigitalcollection.asme.org/data/journals/jvacek/28878/s_011506vaj1.jpeg)

Also, in the vector te and qe I have the time instants and the solution values when the events are encountered. If I have 2 different events to be satisfied (lets call them type 1 and type 2) but I want to plot only the events of type 1, how can I tell Matlab?

Finally, what checks can I perform to be sure the ODE correctly identifies all the events?

Title: How can I use ODE events in MATLAB?


问:
Permissive Math Library for Parameter Statistics in C++
I am looking for recommendations for a C++ math library, with a permissive licence, well suited to calculating a wide variety of statistics on segmentations of timebased parameter data.

I would be particularly interested in anything that can quantify properties of the curve shapes as well as the raw data set.


问:
Need suggestions on how to implement this time stepping for wave equation
I have the following system of equations obtained by implementing Sympletic Euler time scheme to wave equation. I want to model this in Fenics. Here 'u' is the displacement and 'p' is corresponding velocity.

$$ \dot{u} = p $$
$$ u(x, 0) = u_i$$
$$ p(x, 0) = 0$$
$$ u(x, t) = 0 \quad \text{on} \quad \partial \Omega$$

I read in a documentation that the velocity 'p' should be obtained first and thereafter 'u' should be computed. Also, that it must be ensured that we never have to explicitly compute and store the inverse of the mass matrix M.

Here, $M = uvdx$ and $K = \nabla u \nabla v dx$ .

Is this the correct way to implement the above system of equations? Can anyone kindly guide on how this can be corrected / refined.

```<code>u = interpolate(ui,V)
V = FunctionSpace(mesh,"CG",2)

u_p = function(V)
u_new = function(V)
p = interpolate(Constant(0.0),V)

while t =< T
    bc.apply(M,u)        # apply boundary conditions
    solve(M, u_p.vector(), u.vector)
    p_new = p.vector() - dt*K*u_p.vector()
    t += dt              # solve for next timestep
    u_new.vector()[:] = u.vector() + dt*p_new.vector()
    u.assign(u_new)      # update old vector
    p.vector()[:] = p_new
</code>```


问:
Algorithm to determine if two polygons intersect
I'm working on an algorithm which should check if two polygons, described by their vertex coordinates, are: one inside the other, are intersecting or are separated

image below describe this three cases:

i'm thinking about how to do it but i'm not having any idea. Any suggestion?

solved:
i solved using Shamos-Hoey Algorithm and Winding Number method


问:
Allocation under constraint (optimization ?)
I have a problem of allocating x between p entities, each entity has a constant $\rho_i$ such that the allocated part should be proportionnal.

Without loss of generality we can divide the ressources by x and normalise the $\rho_i$ by $\Sigma_i\rho_i$ to see the unconstrained problem has a simple solution:

$$x_i = n \frac{\rho_i}{\Sigma_i\rho_i} $$

My problem is to solve a constrained problem, such that:

$$ a_i \leq x_i \leq b_i $$

I have an iterative algorithm, starting at $x_i = n \frac{\rho_i}{\Sigma_i\rho_i} $ at each iteration we reallocate the ressources that excess the constraint. Theoretically this works because $\Sigma b_i > 1$ and $\Sigma a_i < 1$.

I have the feeling that it could be reformulated as an optimization problem but can't figure out what the objective function would be.


问:
Software Implementing Jiles-Atherton
I have been doing much research on the internet to see if there was a FEM simulator that natively implemented the Jiles-Artherton model whether it be commercial or open source

So far, my conclusion is that there is none. Before I make a final verdict though, I would like to check with the community that this statement is correct. Is there any software package that natively implements the Jiles-Artherton model?

I am looking into this because I am working on finding a focal point to my PhD thesis research and I was thinking about creating a simulator for non-linear magnetic simulations. It would include a number of different features such as SPICE simulation that couples with the FEA simulation and parallel computing (both SPICE and FEA simulation). This would also include implementing a Jiles-Artherton model natively into the software.

I am though concerned that my project is not too much on the researchy side of things and worried this may not be funded since alot of the commercial software packages implement what I would like to do. So, I would also like to ask the community for some feedback, if possible, on the current research topic. Everything is still in an incubation since I start in the coming spring semester. I am asking this because both my professor and I are relatively new to the field.

Title: Software Implementing Jiles-Atherton


问:
Can I make this numerical integration continuously differentiable?
Suppose I have the discrete values $f(x_i)$ for every discrete value $x_i$ greater than some $\varepsilon$, and I want to numerically calculate the following integral:
\begin{equation}
 n = \int_\varepsilon^\infty f(x)Z(x)dx
\end{equation}
with the known function $Z(x)$.

We can approximate $f(x)=f(x_i)$ for $x\in[x_{i-},x_{i+}]$ with $x_{i-}=x_i-\dfrac{\Delta x}{2}$ and $x_{i+}=x_i+\dfrac{\Delta x}{2}$,
and calculate the above integral like:
\begin{equation}
 n \approx f(x_k)\int_{\varepsilon}^{x_{k+}} Z(x)dx + \sum_{j=k+1}^{\infty} f(x_j)\int_{x_{j-}}^{x_{j+}} Z(x)dx
\end{equation}
with $\varepsilon\in[x_{k-},x_{k+}]$. First consider the case of $Z(x)=\sqrt{x-\varepsilon}$, for which the above expression gives a continuously differentiable
approximation with respect to $\varepsilon$.

Next, consider another function $Z(x)=c$. We can check and see that the above approximation is no longer differentiable when $\varepsilon$ moves continuously from one
$[x_{i-},x_{i+}]$ to the next $[x_{(i+1)-},x_{(i+1)+}]$. In order to have a continously differentiable approximation to $n$, we can use a piecewise linear interpolation for
$f(x)$. That is,
\begin{align}
 & f(x) \approx \dfrac{f(x_{j+1})-f(x_j)}{\Delta x}(x-x_j) + f(x_j)  \qquad\quad   x\in [x_{j-},x_{j+}]~,   \\
 & n \approx \int_{\varepsilon}^{x_{k+}} f(x)Z(x)dx + \sum_{j=k+1}^{\infty} \int_{x_{j-}}^{x_{j+}} f(x)Z(x)dx
\end{align}
results in a continuous derivative $\dfrac{dn}{d\varepsilon}$.

Is there some piecewise interpolation for $f(x)$ in order to make the above mentioned numerical approximation continuously differentiable
for the function $Z(x) = \dfrac{1}{\sqrt{x-\varepsilon}}$?


问:
Time discretization of wave equation
I am trying to model the seismic wave equation and have therefore been reading about discretization schemes and their stability. I recently came across an insightful paper on 'Galerkin FEM methods for wave equation' http://www.math.tamu.edu/~bangerth/publications/2010-waves.pdf. In this document and many other documents, for the newmark beta scheme, the final discretised form of wave equation is in terms of acceleration, velocity and displacement. That is how I have seen it in other documents as well.

However, I have also seen another form of Newmark scheme of wave equation, i.e., the following: (Pg 130)

http://www.lmn.pub.ro/~daniel/ElectromagneticModelingDoctoral/Books/Computational%20EM/Bondeson-%20Computational%20%20Electromagnetics.pdf

Now, I am really confused as to how a wave equation can have 2 different forms of Newmark discretisation. I' ll be very grateful if you could shed some light on this or point me in the right direction.

The end purpose is to model it in fenics and DEAL 2

Title: Time discretization of wave equation


问:
understanding interior penalty jump of basis function
If cardinal basis functions are used (i.e. $\psi_{ij}=1$ iff $i=j$, and 0 otherwise) in interior penalty methods for elliptic equations such as SIPG & NIPG, shoudn't the jump in basis functions when computing surface integrals be [$\psi$] = 0 ? If that is the case the second term (consistency term) is always zero, right? How about the penalty term that has $[\psi]$ ? I am sure I am misunderstanding something here. For the poisson equation
$-\nabla. (a\nabla u)=0\  \text{in} \  \Omega$ with $u=0$

Using the symbol v for the basis functions $\psi$ ,

$$a_{\epsilon}(u,v)=l(v)$$ where  $$a_{\epsilon}(u,v)=\sum_{K\in T_h}\int_K a\nabla u\cdot \nabla v-\sum_{e\in \Gamma_h\cup \Gamma_D}\int_e\{a\nabla u\cdot n_e\}[v]+ \sum_{e\in \Gamma_h\cup \Gamma_D}\int_e\epsilon\{a\nabla v\cdot n_e\}[u] + \sum_{e\in \Gamma_h\cup \Gamma_D}\frac{\sigma_0}{|e|}\int_e[u][v]$$


问:
FFT-based Image Rotation Algorithms More Accurate Than Chirp-Z?
We're currently using a Chirp-Z based implementation:

R. W. Cox and R. Tong, "Two- and three-dimensional image rotation using the FFT," IEEE Trans. Image Processing, vol. 8, no. 9, pp. 1297–1299, Sep. 1999

for rotation of astronomical images. I was wondering if there were any algorithms/implementations out there that could also be used to rotate images using FFTs, and were more accurate than the one we're using now. Perhaps any progress made in the past few years in the area?

ContentLicense: CC BY-SA 4.0


问:
Inverse quadratic method only makes one iteration - python
I'm tryin' to implement this method with python3 but i don't solve the problem. This is my code:

```<code>#!/usr/bin/python
# -*- coding: utf-8 -*-

import sympy
import numpy as np

x = sympy.symbols('x')

expr = x ** 2 - 2
fx = sympy.lambdify(x, expr, 'numpy')

def quad_inv(f, x0, x1, x2, tol, maxiter):
    fx0 = f(x0)
    fx1 = f(x1)
    fx2 = f(x2)
    x3 = None
    iters = 0
    error = tol + 1
    if np.sign(tol) < 0:
        raise ValueError('La tolerancia debe ser estrictamente positiva.')

    while (error >= tol) & (iters <= maxiter):
        a = fx1 * fx2 / ((fx0 - fx1) * (fx0 - fx2)) * x0
        b = fx0 * fx2 / ((fx1 - fx0) * (fx1 - fx2)) * x1
        c = fx0 * fx1 / ((fx2 - fx0) * (fx2 - fx1)) * x2
        x3 = a + b + c
        x0 = x1
        x1 = x2
        x2 = x3
        error = abs((x3 - x0) * (x3 - x1) * (x3 - x2))
        print(error)
        iters += 1
    return[x3, iters]

 sol = quad_inv(fx, 1., 1.2, 1.6, 0.0001, 300)
    print('La raíz aproximada con el Método Inversa Cuadrática es:\n{0}'.format(sol[0]))
    print('Se han realizado un total de {0} iteraciones.'.format(sol[1]))
</code>```

I only get 1 iteration. I don't find anything about the error cote with this method.

```<code>La raíz aproximada con el Método Inversa Cuadrática es:
1.4195804195804194
Se han realizado un total de 1 iteraciones.
</code>```


问:
Reference Request: Variational Problem
I want to solve approximately the following variational problem:

Given a function $c:[-1,1]^2\rightarrow [0,1]$, constants $p_1...p_n\in \mathbb{R}^+$, $\alpha_1...\alpha_n\in \mathbb{R}$, and $\beta_1...\beta_n\in \mathbb{R}$ let

$$V =\sup_{f_i,g_i}\left\{\sum_{i=1}^n \;\;p_i\cdot \int_0^1\int_0^1  f_i(x,y)\cdot g_i(x,y)\; dx\; dy \;\;: \sum_{i}\alpha_i\cdot f_i = c = \sum_i \beta_i\cdot g_i \right\}$$

Where the supremum is taken over all square integrable functions
$f_i,g_i:[-1,1]^2\rightarrow [-1,1]$.

Observe that since all functions have domain $[-1,1]^2$ and range
$[0,1]$, the value of each integral $\int_{0}^1\int_{0}^1 f(x,y)\cdot g(x,y) \;dx\; dy$ lies between $0$ and $1$. Therefore $V$ lies between $0$ and $\sum_{i} p_i$.

Assumptions:

All functions $f_i$ and $g_i$ and $c$ can be assumed to be square integrable.

If it simplifies the problem, $c:[-1,1]\rightarrow [0,1]$ can be assumed
to be $c(x,y)=1$ if $x=y$ and $c(x,y)=0$ otherwise.

Questions:

What numerial methods are available to approximate the value of $V$ up to an $\varepsilon$ additive factor? In other words, I want to find a $V'$ such that $|V-V'|\leq \varepsilon$, where $\varepsilon$ is a given precision parameter.

What is the ratio of convergence of such methods?

Title: Reference Request: Variational Problem


问:
Coriolis test problem
I have a eulerian hydrodynamic code, which I have recently implemented the Coriolis force in the integrator. Is there a standard test problem for the Coriolis force?

Since this isn't a Lagrangian code there is no tracking of fluid particles. Right now I am thinking about periodic boundary conditions and setting the entire box to constant density with an initial velocity and watching the entire box evolve. Thankfully since the force is only dependent on velocity and the box's orbital frequency, the medium should evolve homogeneously.

It would be nice to track an individual particle instead of just the velocity evolution of the entire box. If I created an over density, on a static background density, with a non-zero velocity it would be almost analogous. There are some obvious problems with that which make me unsure if it would be a clean test.

Title: Coriolis test problem


问:
problem with understanding the fluid boundary conditions of a 1D probelm
I am having problems understanding the boundary conditions of the problem described in this paper on researchgate

Essentially the problem consists of a one dimensional fluid chamber in contact with a piston on the right hand side. The problem makes use of the ALE mesh moving technique. The ALE is moved in such a way that the elements remain equidistant, so essentially the mesh velocity varies linearly at all times.

My understanding is that on the left hand side, the velocity of the fluid and the mesh are set to zero whilst the pressure and density are linearly interpolated. On the right hand side, the fluid and mesh velocities are set equal to the piston velocity whilst pressure and density are interpolated.

I would appreciate some guidance on those.


问:
Family of linear one-step method convergence question
This question is not homework it is recommended exercises to prepare for the final exam.

Consider the family of linear one-step methods defined by $$y_n = y_{n-1} + h(\theta f_n + (1 - \theta)f_{n-1})$$ where $0\leq \theta \leq 1$.

2) For what values of $\theta$ in the interval $0\leq \theta \leq 1$ is the associated method convergent?

Attempted proof:

Let's first prove convergence of forward Euler i.e. when $\theta = 0$, the forward Euler method says that $$u_{n+1} = u_n + hf(t_n,u_n)$$ and Taylor expansion with remainders of each component gives $$ y_j(t_{n+1}) = y_j(t_n) + hf_j(t_n,y(t_n)) + \frac{1}{2}h^2y_j''(t_n + \theta_j h)$$ $$= y_j(t_n) + hf_j(t_n,y(t_n)) + h\tau_{nj}$$ for some unknown collection of numbers $\theta_j$ between $0$ and $1$. This defines the local truncation error $\tau_n$, the amount by which the exact solution $y$ fails to satisfy the method after one step. Subtraction gives a difference equation for the error $$e_{n+1} = u_{n+1} - y(t_{n+1}) = e_n + h[f(t_n,u_n) - f(t_n,y(t_n))] + h\tau_n$$ Assume $f$ is Lipschitz constant $L$ and the local truncation error satisfies a bound $$||\tau_n|| \leq \tau$$ for all $n$ (which depends on the exact solution) then
$$||e_{n+1}|| \leq ||e_n|| + hL||e_n|| + h\tau = (1 + hL)||e_n||+ h \tau$$

$$\leq (1 + hL)^{2}||e_{n-1}|| + (1 + (1 + hl))h\tau $$
$$\leq \ldots$$
$$\leq (1 + hL)^{n+1}||e_0|| + \frac{(1 + hL)^{n+1} - 1}{(1 + hL) - 1}h\tau$$
since $(1 + hL) \leq e^{hL} = 1 + hL + \frac{1}{2}(hL)^2 + \ldots$, this gives $$||e_n|| \leq e^{nhL}||e_0|| + \frac{1^{nhL} - 1 }{hL}\tau \leq e^{LT}||e_0|| + \frac{e^{LT} - 1}{L}\tau$$ for $0\leq t_n = nh\leq T$. This shows stability. Now the local truncation error satisfies $$h||\tau_n|| = ||\frac{1}{2}h^2y'(t_n + \theta_n h)|| \leq \frac{1}{2}Mh^2$$ assuming $y$ is $C^2$ and $||y''|| \leq u$. Consistency gives a local bound on $\tau$ and stability allows us to conclude convergence $$|u_n - y(t_n)|| \leq e^{LT}||u_0 - y_0|| + \frac{e^{LT} - 1}{LT}\frac{T}{2}Mh \leq O(|u_0 - y_0|) + O(h)$$ Convergence implies that $u_n\rightarrow y_n$ on the whole interval $[0,T]$ as $h\rightarrow 0$ and $u_0\rightarrow y_0$. A similar proof can be done for any one-step method. A little extra work is required for implicit methods such as implicit Euler $$u_{n+1} = u_n + hf(t_{n+1},u_{n+1})$$ because the right hand side depends on the new unknown $u_{n+1}$. Assume you can solve the implicit equation for $u_{n+1}$ at each step so we can estimate the errors as for explicit Euler $$||e_{n+1}|| \leq e_n + hL||e_{n+1}|| + h||\tau_n||$$ collect each terms on the left, and for $hL \leq \frac{1}{2}$ we can then divide by $1 - hL$ to get $$||e_{n+1}|| \leq \frac{1}{1 - hL}(||e_n|| + h||\tau_n||)$$ using the additional inequality $$\frac{1}{1 - hL}\leq 1 + 2hL \ \ \text{for} \ \ hL \leq 1/2$$ gives stability as for Euler's method, and the consistency computation follows.


问:
Why not solve molecular dynamics trajectories exactly?
The basis of Molecular Dynamics simulation is the integration of newton's equations of motion:

$$\frac{F}{m} = \frac{d^2r}{dt^2}$$

I understand that there are various methods based on Taylor's expansion, to integrate forward in time and calculate the new position. Doing this always leaves an error term in the expansion, however.

To solve for position exactly, I can integrate both sides:

$$r = \frac{Ft^2}{m} + c_1 t + c_2$$

If $r(t=0) = 0$ then $c_2 = 0$.

If it's possible to determine $c_1$ too (?) we have an exact solution, and so I can calculate the next position at $t+dt$ with perfect accuracy, unlike with Taylor expansions.

Can anyone show me where my understanding is wrong? Thank you.


问:
How to construct a subring of a polynomial ring in Magma
This is a question about the computer algebra system Magma. I have been looking for a place to ask this type of question on the SE network and scicomp.SE was suggested to me; hopefully it finds a home here.

Suppose one has constructed a polynomial algebra $A$ over a ring $R$ in Magma. How does one construct the sub-$R$-algebra of $A$ generated by a given list of elements of $A$?

This seems to me to be a very basic operation so I can't believe there isn't a way to do it, but I haven't so far found it in the handbook. (I see functionality to construct subalgebras of matrix algebras and of endomorphism rings of abelian varieties, but not polynomial rings.)

EDIT 11/3/15:

1) As the functionality of Magma may vary with the ring $R$, take $R=\mathbb{Z}$.

2) hardmath suggests in comments to construct the desired subring as the quotient of another polynomial ring. I take him to mean as the image of a map from another polynomial ring. I tried to implement this in a test case as follows:

```<code>A := PolynomialRing(IntegerRing(), 3);
s1 := ElementarySymmetricPolynomial(A,1);
s2 := ElementarySymmetricPolynomial(A,2);
s3 := ElementarySymmetricPolynomial(A,3);
B := PolynomialRing(IntegerRing(),3);
f := hom< B ->A | s1,s2,s3>;
S := Image(f);
</code>```

I get a Runtime error: No constructor provided for this type of object message when Magma tries to implement Image(f).


问:
Numerically calculating the divergence of a set of oriented points
Say I have a set of oriented points at locations $\vec{v_i}$ with each some direction $\vec{n_i}$, in practice they represent the normals of some surface that has been non uniformly sampled. How would I calculate the divergence of this vector field? Standard finite-difference methods would obviously not work because there is no grid.

Can I  perhaps assume a grid and use nearest neighbours to infer the vector field at each point (I think it's called "splatting")? The discretized vector field would then have at each point a weighted average of nearby points. This would then allow me to use standard FD-methods to calculate (or rather estimate) the divergence at each gridpoint. Is this a common method, or are there better methods?

Edit:
Asked it on MO, got an answer.


问:
Do collocated grid arrangements definitely result in the checkerboard effect?
I understand the checkerboard effect due to the use of collocated grid arrangements in FVM. However, I wanted to know whether this problem is definitely bound to effect the results? For instance, I have modelled some problems on 1D structured grids with a collocated arrangement and have produced smooth results. Perhaps I am mistaken in how exactly the results should be effected? </p>


问:
Resolving a stiff hyperbolic problem with Neumann boundary conditions
I am trying to numerically resolve the equation for an Euler-Bernoulli beam that is inextensible, unshearable, and subject to planar deformations:

$$\rho I(s) \frac{\partial^2 \theta}{\partial t^2}(s,t)
= E \frac{\partial}{\partial s} \left[I(s) \left( \frac{\partial \theta}{\partial s}(s,t) - \kappa_p(s,\theta,t) \right) \right]\,,\,\,\,\,\, t>0\,,\,\,\,\,\, 0<s<L\,,$$

$$\frac{\partial \theta}{\partial s}(0,t) = \frac{\partial \theta}{\partial s}(L,t) = 0\,,$$

$$\theta(s,0) = \frac{\partial \theta}{\partial t}(s,0) = 0\,,$$

where $\kappa_p(s,\theta,t)$ is a preferred curvature. The ends of the beam are force-free and moment-free. Here, $E \gg 1$ so that $\frac{\partial \theta}{\partial s}(s,t) \approx \kappa_p(s,\theta,t)$. I cannot directory approximate $\theta(s,t)$ with $\kappa_p(s,\theta,t)$ since $\kappa_p(s,\theta,t)$ depends on $\theta(s,t)$. I have been reading into what may be a good numerical method to use to resolve this problem, however, I am not finding very good solutions since this equation has a stiff part in both the linear term and the nonlinear forcing term. An implicit method would free me from the minisulce time-step coming from $E$, however, the phase will not be resolved properly. I can use an implicit method for the linear terms, however, using something like multi-level approximation in time for the source term also introduces errors. Also, an implicit method may not resolve the phase well if I take big time-steps. What would be the best way to resolve this problem? If the area moment of inertia, $I(s) = 0$, at the end points, would it be a good idea to simply replace the geometry of the beam where the endpoints have a finite, but small $I(s)$ so that the PDE isn't singular or is there a nice way that I can numerically treat this singularity?

Aside: To show that the system is hyperbolic, I can rewrite the PDE as a system of equations and then compute the eigenvalues of the linear part.

$$ \rho I(s) \frac{\partial \theta}{\partial t}
= \frac{\partial \omega}{\partial s} $$

$$ \frac{\partial \omega}{\partial t} = E \left[I(s) \left( \frac{\partial \theta}{\partial s}(s,t) - \kappa_p(s,\theta,t) \right) \right] $$

Now, if $I(s)$ is strictly positive, then the eigenvalues are $\lambda_{\pm} = \pm \sqrt{E/\rho}$ which is a real quantity in this problem. Thus, the problem is hyperbolic. For non-negative $I(s)$, one recovers the same eigenvalues when taking the limit. Note that the stiffness for the problem is:

$$ S = \frac{|\lambda_+|}{|\lambda_-|} = 1 $$

however, depending on the numerical method, the variation in $I(s)$, and if I am not careful about the treatment of the singular points (from $I(0) = I(L) = 0$), then I will get numerical stiffness.


问:
Modeling First Order Parabolic PDE (Battery Storage Model)
I'm trying to solve the following first order parabolic partial differential equation,

\begin{equation*}
X \frac{\partial V}{ \partial Q} = -\frac{1}{2} \sigma^2 \frac{\partial^2 V}{\partial X^2} + rV + p \min (X, 0)
\end{equation*}

subject to conditions
\begin{equation*}
\begin{split}
V &= 0 \hspace{27pt} \text{as} \hspace{20pt} X \rightarrow \infty \\
V &= pQ \hspace{20pt} \text{as} \hspace{20pt} X \rightarrow -\infty \\
\frac{1}{2} \sigma^2 \frac{\partial^2 V}{\partial X^2} - rV &= 0 \hspace{20pt} \text{if} \hspace{20pt} X < 0 \hspace{10pt} \text{and} \hspace{10pt} Q = 0 \\
\frac{1}{2} \sigma^2 \frac{\partial^2 V}{\partial X^2} - rV &= 0 \hspace{20pt} \text{if} \hspace{20pt} X > 0 \hspace{10pt} \text{and} \hspace{10pt} Q = Q_{max}
\end{split}
\end{equation*}

What I understand is that Matlab pdepe solver solves
\begin{equation*}
c \Big(x, t, u, \frac{\partial u}{\partial x} \Big) \frac{\partial u}{\partial t} = x^{-m} \frac{\partial }{\partial x} \Big( x^{m} f \Big( x, t, u,  \frac{\partial u}{\partial x} \Big) \Big) + s \Big( x, t, u, \frac{\partial u}{\partial x} \Big)
\end{equation*}

So if redefine my model,
\begin{equation*}
X \frac{\partial V}{ \partial Q} = -\frac{1}{2} \sigma^2 \frac{\partial^2 V}{\partial X^2} + rV + p \min (X, 0)
\end{equation*}
I can let
\begin{equation*}
\begin{split}
c &= X \\
f &= -\frac{1}{2} \sigma^2 \frac{\partial V}{\partial X} \\
s &= rV + p \min (X, 0)  \\
m &= 0
\end{split}
\end{equation*}

My question is: how do I get the boundary condition supposedly in the form
\begin{equation*}
p(x,t,u)  + q(x,t) f \Big( x,t,u, \frac{\partial u}{\partial x} \Big) = 0
\end{equation*}

and how do I capture conditions as stated above? The domain should be divided into two regions, one (X < 0) for solving upwards, and the other (X > 0) in which solve downwards. Along the boundary X = 0 the solutions must be matched.

I'm a bit lost at this point. Any insight would really help. Thanks!

Edit:
If I let the following be p,f, and q for the boundary condition
\begin{equation*}
\begin{split}
p &= -rV \\
f &= -\frac{1}{2} \sigma^2 \frac{\partial V}{\partial X} \\
q &= 1
\end{split}
\end{equation*}
would this be correct in capturing the above mentioned conditions?

This is the code in Matlab:

```<code>function pde
global sigma r p
sigma = 0.5;
r = 0.01;
p = 1;

x_inf = 25;
N = 10;
M = 20;
delta_x = x_inf/N;
q_max = 1;
delta_q = q_max/M;

m = 0;
x = linspace(0, delta_x, 2*N);
t = linspace(0, delta_q, M);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1);

% A surface plot is often a good way to study a solution.
surf(x,t,u)
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,q,v,DvDx)
global sigma r p
c = x;
f = -(1/2)*sigma*DvDx;
s = r*v + p*min(x,0);
% --------------------------------------------------------------
function u0 = pdex1ic(~)
u0 = 0;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,vl,xr,vr,q)
global r
pl = -r*vl;
ql = 1;
pr = -r*vr;
qr = 1;
</code>```

Any insight at all on the steps to solving this (whether in Matlab with pdepe or not) will be very much appreciated!


问:
CFL Condition and Convection Diffusion Equation in 2D
I am solving the convection-diffusion equation in 2D using Finite Differences with the $\theta$ scheme. The velocity of the fluid and the diffusion coefficient is low in my case (in the range of $10^{-15}$). The boundary conditions are: open boundary on two sides (Robin boundary), and Dirichlet boundary $u=0$ on the other two.

I have some questions

Is the CFL-Number of any importance when solving the Convection Diffusion Equation in 2D using the $\theta$ scheme and Finite Differences?

How does the diffusion coefficient factor into the CFL-Condition?

<s> I know the implicit case is supposed to be stable for all time steps and step sizes, but I get ugly oscillations. Is it only important for the explicit case?</s>
(edit: I have no more oscillations)


问:
Center of mass in systems with periodic boundary conditions
I have a question about the calculation of center of mass (COM) in systems with
periodic boundary conditions. There is a
method introduced here:

https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions

But by this method the COM can't be calculated properly, see the image bellow, because the system is extended over the whole box and not only the lower right corner but also the upper right one is considered as periodic image.
Is there a straightforward way for calculating COM in systems like this?


问:
Boundary Conditions for the given PDE
I'm working on the Black-Scholes equation, but I'm pretty new to financial modeling. Right now, I am trying to understand the Black-Scholes PDE. I understand that the Black-Scholes equation is given by
\begin{equation*}
\frac{\partial C}{\partial t} + \frac{1}{2}\sigma^2 S^2\frac{\partial^2 C}{\partial S^2} + rS \frac{\partial C}{\partial S} - rC = 0
\end{equation*}
with initial condition
\begin{equation*}
C(S,T) = \max (S-K, 0)
\end{equation*}
and boundary conditions
\begin{equation*}
C(0,t) = 0 \hspace{35pt} C(S,t) \rightarrow S \text{ as } S \rightarrow \infty
\end{equation*}
and $C(S,t)$ is defined over $0 < S < \infty$, $0 \leq t \leq T$.

This can be further transformed and simplified into a heat diffusion equation as described here.

If we make the following change of variable
\begin{equation*}
u = e^{-r\tau}C \hspace{20pt} \text{ or } \hspace{20pt} C = ue^{r\tau}
\end{equation*}
and
\begin{equation*}
S = e^x \hspace{25pt} \text{ and }  \hspace{25pt} t = T- \tau \hspace{20pt}
\end{equation*}
we get the transformed heat equation

\begin{equation*}
\frac{\partial u}{\partial \tau} = \frac{\partial^2 u}{\partial x^2} + (k-1)\frac{\partial u}{\partial x} - ku
\end{equation*}

Where $k = \frac{2r}{\sigma^2}$. The following matlab code implements this. My question is, what exactly is the form of the boundary conditions for the the transformed equation? I can't seem to understand the parameters (related to the boundary conditions) given in the Matlab code. Any related literature would be highly appreciated.

Title: Boundary Conditions for the given PDE


问:
Clustering state-space models according to similarity
I have a plant that can be modeled with nonlinear equations as
$$f(x,u,p)=0$$
where $x$ is the state vector, $u$ the control vector, and $p$ the model parameters.

In order to control this system, I compute the discrete state-space model
$$x_{k+1}=A_p x_k+B_p u_k$$
where $A_p$ and $B_p$ depend on the parameters $p$.

Now, I consider $N$ different $p$ vectors referring to possible malfunctions of the plant. Thus, I get $N$ possible state-space models. I am searching for a methodology to cluster these state-space models according to their "similarity" and end up with a subset of $M(<N)$ cluster state-space models.

A trivial example is when the plant has two identical components sharing some task, when one malfunctions ($p_1$), I get the state-space model $(A_{p_1},B_{p_1})$. When the second malfunctions ($p_2$), I get $(A_{p_2},B_{p_2})$. However, the two state-space models are the same and I need only to consider one of them.

Any references to papers or techniques are appreciated.


问:
Is it worth switching to timesteppers provided by PETSc if I can't write down a Jacobian for my problem? Case study with "the amoeba" toy problem
I am considering using petsc4py instead of scipy.integrate.odeint (which is a wrapper for Fortran solvers) for a problem involving the solution of a system of ODEs. The problem has the potential to be stiff. Writing down its Jacobian is very hard.

So far, I have been able to produce reasonable speed gains by writing the RHS functions in "something like C" (using either numba or Cython). I'd like to get even more performance out, hence my consideration of PETSc.

Introducing a toy version of the problem: The Amoeba

Consider the following figure, which defines the basics of an amoeba:

So, an amoeba consists of $N$ nodes and edges, constrained to a 2D space. The position of the nodes describing the amoeba evolve with time.

Let the position of the $n$th ($n \in \{0, 1, ..., N-1\}$) node be described a position vector, $\textbf{x}_n$:

$$\textbf{x}_n = \left[\begin{array}{cc} x_{0,n} \\ x_{1,n} \end{array}\right]$$

For the purposes of simplification, assume that the amoeba lives in a world with first order equations of motion, so given an applied force $\textbf{F}$ on a node, the nodal velocity $\dot{\textbf{x}}$ is produced:

$$ \textbf{F} = c \dot{\textbf{x}}$$

where $c$ is a constant I pick.

We are very interested in tracking the position of each node,so immediately we produce $2n$ ODEs (one ODE for each spatial component) of interest (assume that edges are elastic springs):

$$ \frac{\textrm{d}\textbf{p}_n}{dt} = \text{elastic edge forces} + \textbf{F}_{n, active}$$

$\textbf{F}_{n, active}$ is the force generated by the amoeba at node $n$. One can imagine a wide range of rules that could be devised to determine how $\textbf{F}_{n, active}$ is determined, and I suggest that we consider the following rules:

1) There are red balls and blue balls inside the amoeba.
^
2) Both red and blue balls have "switched-on" (symbolized $r^*$/$b^*$) and "switched-off" (symbolized $r$/$b$) state.

3) Both red and blue balls can be found either in the body of the amoeba where they are always switched off, or at the nodes of the amoeba, where they can be switched on or switched off.

4) The total number of red balls regardless of location of status is a constant $R$ and similarly for blue balls, a constant $B$

5) When red or blue balls are in the body, they are free to jump onto any node they want, but when they are on nodes, some of them might be induced to move to neighbouring nodes due to some sort of "diffusion"

6) Switched on red balls at a node cause a force that attempts to reduce the area of the amoeba

7) Switched on blue balls at a node cause a force that attempts to increase the area of the amoeba

8) Switched on red balls tend to switch off blue balls that are switched on, and vice versa, switched on blue balls tend to switch off red balls that are switched on

9) ...

Okay, you get my point. I can come up with a model that is fairly convoluted in terms of how the various variables of interest interact.

Due to the large number of equations involved, it is already tedious to think about writing downa Jacobian. However, since we have access to a computer, it may be that we can write functions governing a particular interaction that do not have neat analytical forms (let alone whether or not their derivatives have neat analytical forms), so we might have a mess of piecewise functions needed to approximate them if we were to go about still trying to produce a Jacobian...

Things get even worse if my problem involves not just one amoeba, but multiple amoebae (which is when things ARE interesting), with new interaction rules enforcing that amoebae cannot violate each others' areas, and so on.

All the toy examples I see of PETSc time stepping problems have Jacobians defined, so I wonder if I would even get a speed gain going from switching to it, if perhaps one of the reasons why I have a high computational cost is due to not being able to provide a Jacobian function?


问:
Duality/Lagrangian condition and Variational Inequality of a cost functional
Given the functional $J_A(y,u) = \frac{1}{2}||y-y_d||^2_{L^2(w,\omega)} + \frac{\lambda}{2}||u||^2_{L^2(w^{-1}, \omega)}$ where w is a function belongs to Muckenhoupt class.

Given the optimization problem: min $J_A$ subject to

-$div(A\triangledown{y}) = u$ in $\omega$. Assume matrix $A$ satisfies these inequalities for almost every $x\in \omega$,

$w(x)|\zeta|^2 \leq C(\zeta^{T}\cdot A)$,
$\ \zeta^{T}\cdot A \leq D(w(x)|\zeta|^2)$ for some constants $C$ and $D$

$y = 0$ on $\partial \omega$,

and the control constraints $u\in U_A$ where $U_A$ is a nonempty, closed and convex set of $L^{2}(w^{-1}, \omega)$.

I'm interested in finding the dual problem + Lagrangian condition for this optimization problem. I also tried quite hard to prove the following inequality but still without any success: $(DJ_A(\bar{u}), (u-\bar{u}))\geq 0$ where $\bar{u} = $ optimal control solution to the optimization problem above.

Hope someone can help me with either questions above. It's quite interesting, but also challenging, so I would really appreciate any help you can give.


问:
How to get mathematical model from a data set with MATLAB
I have some values imported from excel about the annual sell of a product, with 3 var:

month,
price,
sold items

I have seen how to plot a chart using curve fitting tool but I need to define a math expression that can describe the relation between these variables.

How can I do this in matlab? Is there any tool that can achieve this automatically. If you look at the chart should exist a way to use an approximation that can mime the behavior(maybe with spline functions or n-point-hermit), is enough one approximated expression of the model that for the 3 var replicates a similar graph if plotted.


问:
Order of accuracy of linearised vs non-linear system
Does the order of accuracy of a combination of schemes applied to solve a system of non-linear equations, match those of the same schemes applied to the linearised version of the system?

In other words, does linearisation alter the order of accuracy of a scheme?


问:
Not sure if sparse quadrature routine is working correctly?
I have written a sparse quadrature routine to integrate multidimensional Gaussian integrals. I'm trying to show convergence plots in my thesis and demonstrate that the code works correctly. I am in the process of benchmarking this using multidimensional integrals with closed forms. With normal numerical schemes one would assume the relative error to decreases continuously a with higher number of quadrature points. However the smolyak algorithm selects only those tensor products obeying the following criteria:

$$\ell-M+1\leq||i||_1\leq\ell$$

Where $M$ is the number of dimensions in the integral, $\ell$ is the order  and $\|i\|_1$ the summation of the number of quadrature points used for each of the $M$ one-dimensional sets. This results in a sparse tensor product. As we can see when large $\ell$ is used the lower order tensor products are dismissed, particularly when $\ell>>M$. Hence would one still expect the error to decrease continuously with a higher order $\ell$? Within my code I am seeing that for higher $\ell$, the error plateaus, as $M$ is increased this plateaus moves such that the accuracy is decreased. For example when $M$ is 1, the error plateaus at 1E-16 i.e. double precision limit as expected. But when $M=5$ the plateau is at around 1E-12, this seem quite large and I've tested this for several different multidimensional integrals. Is this correct? For $\ell=35$ we have

$$31\leq||i||_1\leq35$$

So we are losing a lot of tensor products. Am I being to pedantic? I will only every require the first four significant figures due to input data.


问:
DIIS method to accelerate SCF convergence for stretched geometries
I am implementing from scratch an Hartree-Fock calculation in the STO-3G basis set to perform Born-Oppenheimer molecular dynamics. I have a Restricted Hartree-Fock procedure that can reproduce very well the total energies of Ref. [1] and of Gaussian09 with pre-optimized molecular geometries.

Now, during a Born-Oppenhaimer MD simulation, I will have geometries that can be quite far from equilibrium. I saw that my SCF cycle does not converge for stretched geometries: the SCF cycle for the $CO$ molecule with a bond length of 2.132 bohr converges, while for a bond length of 3.132 bohr I reach the maximal number of allowed iterations (500).

In Table 10.4 of Ref. [2] it is clear that an implementation of the DIIS (Direct Inversion of the iterative Subspace) method speedup the SCF convergence. In addition to this speedup, it seems that the convergence for a stretched molecule can be easily reached with this method. For this reason I implemented the DIIS algorithm following Ref. [3].

I believe I implemented the DIIS method quite correctly, since I was able to obtain a speedup of SCF convergence for pre-optimized molecular geometries. In the following table you can find the number of simple SCF and SCF-DIIS steps I have, for different molecules:

```<code>Molecule    SCF     DIIS
H2O         16      9
CO          53      19
HeH+        8       7
CH4         11      8
FH          10      7
O2          41      18
N2          110     22
</code>```

In every molecule I tested I get a speedup, which is quite impressive where a lot of normal SCF cycles were needed. This is a good result, but unfortunately the convergence of stretched molecules is not improved. By doubling the distances (in bohr) of the pre-optimized structure, the DIIS algorithm (as well as the normal SCF) fails to converge.

There is some step to add to the original DIIS algorithm of Ref. [3] that I am not considering? How I can make the SCF cycle convergent for stretched molecules, provided that these calculations converge in Gaussian09?

EDIT

Gaussian09 brakes down as my program does with the option SCF=NoDIIS. With the option SCF=DIIS the SCF converges even for very distorted molecules.

[1] A. Szabo and N. Ostlund, Modern Quantum Chemistry, Dover, 1996.

[2] T. Helgaker, P. Jørgensen and J. Olsen, Molecular Electronic-Structure Theory, Wiley, 2000.

[3] P. Pulay, Improved SCF Convergence Acceleration, Journal of Computational Chemistry, 1982.


问:
Interior nodes of a closed graph?
Does anybody know if any graph partitioner library such as Metis, Scotch, or Zoltan can (besides splitting a domain), differentiate between internal (i) and boundary (b) nodes?

Title: Interior nodes of a closed graph?


问:
Memory requirement to find eigenvalues and -vectors of large sparse matrix
How can I estimate how much memory will be needed to find eigenvalues and eigenvectors of a given large sparse matrix?

I have a real symmetric matrix with roughly $5 \times 10^4$ rows and columns, and an average of $10$ nonzero elements per row. I would like to find the smallest eigenvalue and the corresponding eigenvector, using the built-in Eigensystem function in Mathematica (which treats the matrix as sparse and uses an ARPACK Arnoldi algorithm). Is there a simple way of estimating how much memory this will take?


问:
Nitsche' method coercivity again
I know that for a dependent-mesh norm Nitsche' method for boundary conditions is coercive. But how to prove this for the usual norms for $H^1$?</p>

Title: Nitsche' method coercivity again


问:
End conditions on cubic splines interpolation
I'm currently working on cubic splines interpolation. My question has to do with the end conditions. Instead of using natural or clamped cubic splines, I am asked to use the following condition :

$S_{1} = S_{2} $ and $S_{n} = S_{n-1} $

but I am not pretty sure how this is applied.

Could anyone help me clear this out?

EDIT

I am editting due to the first comment. On which points do I appply this condition? $S_{1}(x_{2}) = S_{2}(x_{2}) $ is a condition that I've already used and $S_{1}(x_{1}) = S_{2}(x_{2}) $ leads to a contradiction...So is it something like $S_{1}(x_{1}) = S_{2}(x_{1}) $ ?


问:
Solving a nonlinear problem with CDF
I'm trying to solve this problem:

$\begin{cases}
\partial_t E=-k\left([f(\rho)-i.\left[\delta+\frac{1}{2}a\left[\dfrac{\nabla^2_{\bot}}{4}+1-\rho^2\right]\right]]E - 2CP\right)\\
\partial_t P=-\gamma_{\bot}[(1+i\delta)P + ED]\\
\partial_t D=-\gamma_{\parallel}(D-\chi(\rho)-\frac{1}{2}(E^*P+EP^*))
\end{cases}$

Here $f(\rho)$ and $\chi(\rho)$ are radial functions. and the Laplacian is only en x and y (transversal Laplacian).

$\rho$ is the radial coordinate in a cylindrical coordinate system.

with

$E=E_x,E_y$ ;and   $E_x , E_y \in \mathbb{C}$

$P=P_x,P_y$ ; and  $P_x , P_y \in \mathbb{C}$

so $E_x,\ E_y,\ P_x,\ P_y$ are complex functions. which amounts to a total of 9 equations.

I don't know much about CDF, I tried doing some explicit forward method, since I saw a Navier-Stokes problem being solved that way, but its unstable:

from the equation for the evolution of $P$:
$$P^{j+1}_{n,m}=-\gamma_{\bot}[(1+i\delta)P^j_{n,m}+E^j_{n,m}D^j_{n,m}]\Delta t + P^j_{n,m}$$

from the equation for the evolution of $D$:
$$D^{j+1}_{n,m}=-\gamma_{\parallel}\Delta t (D^j_{n,m}-\chi(\rho)_{n,m}-\frac{1}{2}((E^j_{n,m})^*P^j_{n,m}+E^j_{n,m}(P^j_{n,m})^*))+D^j_{n,m}$$

from the equation fo the evolution of $E$:

$\partial_t E=i\dfrac{ka}{8}\nabla^2_{\bot}E-k([f(\rho)E-i[\delta+\frac{1}{2}a[1-\rho^2]]E-P)\\$

then: $$\Re(E)^{j+1}_{n,m}= -\frac{ka}{8}\Delta t
\left(\frac{\Im(E)_{n+1,m}^{j} - 2\Im(E)_{n,m}^{j} + \Im(E)_{n-1,m}^{j}}{\Delta x^2}+\frac{\Im(E)_{n,m+1}^{j} - 2\Im(E)_{n,l}^{m} + \Im(E)_{n,m-1}^{j}}{\Delta y^2}\right)
-\left(f(\rho)_{n,m}\Re(E)^j_{n,m} + (\delta+\frac{a}{2}[1-\rho^2])\Im(E)^j_{n,m} - \Re(P)^j_{n,m}\right)k\Delta t + \Re(E)^{j}_{n,m}\\$$

witch holds for $\Re(E_x)$ and $\Re(E_y)$

$$\Im(E)^{j+1}_{n,m}=
\frac{ka}{8}\Delta t
\left(\frac{\Re(E)_{n+1,m}^{j} - 2\Re(E)_{n,m}^{j}+\Re(E)_{n-1,m}^{j}}{\Delta x^2} + \frac{\Re(E)_{n,m+1}^{j} - 2\Re(E)_{n,m}^{j} + \Re(E)_{n,m-1}^{j}}{\Delta y^2}\right)
- \left(f(\rho)_{n,m}\Im(E)^j_{n,m} - (\delta+\frac{a}{2}[1-\rho^2])\Re(E)^j_{n,m} - \Im(P)^j_{n,m}\right)k\Delta t + \Im(E)^{j}_{n,m}\\$$
witch holds for $\Im(E_x)$ and $\Im(E_y)$

I have been reading up on CDF for several weeks, but still feel kind of lost.

Any suggestions about winch method should I use, and why, or maybe some book.

Title: Solving a nonlinear problem with CDF


问:
Detecting and joining series of line segments that run along each other
Given: Several circular series of map GPS coordinates for several bus routes. The GPS coordinates are not all equal when they run along the same road. The number of GPS coordinates for a single bus route runs from 140-600.

Problem: When the points are downloaded drawn, the raw routes don't follow each other perfectly and especially when zoomed-out you can't see some of the routes. I want to somehow show the multiple routes that run along a road, mostly likely dashed with multiple colors. But in order to do that, I need to detect the segments of the route that run really close to each other.

What would be perfect is finding the segments of the routes and mutating their GPS coords to be shared among them. Once that's done I can easily render it.

So my question is: are there any good algorithms for this purpose, or some combination/tweaks of other algorithms I can use? It all has to be done autonomously in an app, not hand-done/hand-checked.

If you need any other information please ask.


问:
Is this the correct way to calculate the Madelung constant?
I'm trying to compute the Madelung constant of the ZnS lattice.
The method is as follows:
The lattice is a face-centered-cubic with basis $S^-: (0,0,0)$ and $Zn^+: (1/4,1/4,1/4)$.
The Madelung constant is
$$\alpha=\sum_j \frac{q_j}{r_j}$$
where $q_j$ is the charge (in this case divided by the number of neighboring blocks) of $j$
and $r_j$ the distance to a fixed ion.
Evjen's method tells us that we need to form blocks of neutral charge to obtain a good
convergence.

I've chosen to start with two cubes and the origin is placed in the center of the face that
is common to both cubes (so I have a block of size 1x1x2).
Then, two get more precision I add two cubes in each direction (so in the second iteration
I have a block of size 3x3x4).

The problem is that this method converges to $\alpha\approx 5.35$, when it should be $\alpha\approx 1.64$.
Is this the correct block expansion?  I've tried with other block and I get the same value.

I haven't included the code, because I think that the problem is in the block chosen.


问:
convergence of a method
I  want to show convergence of a finite element  method for a higher order equation.I have a coupled equations that solved together and gives two variables as answer $[u, v]$.
$$w+\Delta u=0$$
$$\Delta w=f$$
for finding the answer
$$[u,w]=K^{-1}F,$$
Unfortunately, I am really confused about that and have a lot of questions such as:
What space should I used?  Should I show convergence for both of variables together or separately?
If I have a bilinear form what do I need to show (boundedness, coercively and so existence, uniqueness )? </p>

<p>Does anyone know a good reference that I can find my answers? </p>

Title: convergence of a method


问:
Numerically solving geodesic differential equations with a priori knowledge of the Riemann curvature tensor
The geodesic differential equations are given as

\begin{align}
\frac{d^2 x^j}{ds^2} + \Gamma^{\phantom{h}j}_{h\phantom{j}k}\frac{dx^h}{ds}\frac{dx^k}{ds} = 0,
\end{align}

where the $\Gamma^{\phantom{h}j}_{h\phantom{j}k}$ is the Christoffel symbol of second kind.

According to T. A. Moore - A general relativity workbook,

the Riemann curvature tensor quantifies the relative variation of
  initially parallel geodesics. Since such geodesics only variate relative to each other in a curved space, the Riemann tensor will be zero everywhere in flat space.

If one can conclusively distinguish flat space vs. curved space, can such information have any usefulness when attempting to solve numerically the geodesic differential equations


问:
Meshing two parts together
I have two parts, which need to be meshed. One part is located inside another. To simplify the problem, imagine two cubes: one inside another.

Geometries are stored in STL files. From CAD point of view, one file contains an internal cube, another contains cube with a hole for an internal cube. All dimensions fits (there are no intersections). My task is to create a conformal tetrahedral mesh suitable for mechanical FEM analysis - without duplicated nodes on shared faces. Each element should contain information to which part it belongs. File format does not matter but Abaqus is preferred.

I decided to use gmsh, but I cannot mesh both parts together.
I tried the following commands:

File -> Open Internal cube

Add volume to internal cube

File -> Merge External cube

Add volume to external cube

This produces mesh but elements in internal cube are duplicated. It looks like gmsh does not treat hole (in external cube) as a void. I also tried to use Coherence in gmsh but it failed with error:

Error   : Self intersecting surface mesh, computing intersections (this could take a while)

Error   : Surface mesh is wrong, cannot do the 3D mesh

I assume, that my geometries are fine:

Internal cube: http://pastie.org/10673562

External cube: http://pastie.org/10673564

How to create such a mesh which satisfy my requirements?

Title: Meshing two parts together


问:
Poisson equation in frequency domain
I need some help in numerically solving the nonlinear Poisson's equation for electrons in frequency domain. The steady-state equation is:
\begin{equation}
\nabla.(\epsilon\nabla\varphi) = q\left(n_i\exp\left(\dfrac{\varphi-\Phi_n}{V_\text{T}}\right) -N_D\right)
\,.
\end{equation}

In small-signal analysis, we have $\varphi = \varphi_{\text{DC}} + \delta\varphi$, with $\delta\varphi <<\varphi_{\text{DC}}$.
This yields:
\begin{equation}
\nabla.\Big[\epsilon\nabla(\varphi_{\text{DC}} + \delta\varphi)\Big] =
q\left(n_i\exp\left(\dfrac{\varphi_{\text{DC}} + \delta\varphi-\Phi_n}{V_\text{T}}\right) -N_D\right)
\end{equation}
which I can further simplify to:
\begin{equation}
\nabla.(\epsilon\nabla\varphi_{\text{DC}}) + \nabla.(\epsilon\nabla\delta\varphi) =
q\left(n_i\exp\left(\dfrac{\varphi_{\text{DC}} + \delta\varphi-\Phi_n}{V_\text{T}}\right) -N_D\right)
\,\,.
\end{equation}

I do not know how to proceed from here. That is, how should I treat the exponential term on the RHS and what are the basic assumptions in doing so?

Title: Poisson equation in frequency domain


问:
Which topics are hot in computational science?
I am accepted to computational science master program. The master program just started to accept students in our university and today each group of 3 students was assigned to an advisor teacher (we are only 13 students). Our advisor asked us to research for valuable topics in Computational Science, so I thought why don't I ask to experienced scientists in this forum.

My group is combined of 1 electronics engineer (me), 1 mathematician and 1 computer engineer. So, our strengths are math and software programming.


问:
Divide and Conquer division algorithm explained (as used in GMP bignum)
I am trying to understand the divide and conquer division algorithm that is used in the GMP bignum arithmetic library.

The code is very optimised and that makes it somewhat hard to understand. the doc does not really offer an explanation of how it is done, (and the references do not point to an actual algorithm, not one that I could see anyway).

I understand the concept of divide and conquer for sorting, but I am not sure I understand how it can be used in a division.

How can $x \over y$ be calculated using divide and conquer?

Can someone give a small example of how it is done?


问:
Extracting sinograms from tomography projections
my name is Lorenzo and I am a postdoc researcher in Italy. My job is related to tomography imaging using synchrotron radiation. This is a new field for me and I start to face with Matlab to write some code.

I am completely new with tomography imaging and Matlab shows new challenges for me. My actual problem is to create sinograms from a stacking of parallel image projections. For whose are not in the field, a sinogram is a map detecting the projection on the detector of the position of a feature in the sample, as function of angle between the X-Ray beam and the sample.

What I get from an experiment is a series of 2D radiography at different angles, you can see it as a rectangular "volume" where dimensions are the number of rows and columns in the single projection, respectively, and volume is given from the number of angles. A sinogram is just a transversal cut of this volume. This means reading the volume not from lateral side but from the top, so I will create a new array of images, whose dimensions are the number of columns and projection, and array lenght the number of rows in the projections. To give numbers in this experiment I have 4000 projections with size 2048x1370 pixels, so for my computer skills this is a huge computational problem.

I need your help to perform some operations faster. My code in a first part allocates an array to contain all images, this is a 34 Gb array but I have 130 Gb RAM, so no problem. The code to perform this operation uses a loop of imread:

```<code>for i=2:num_proj
    filename=strcat(path_im,list_proj(i).name);
    image=imread(filename);
    imArray(:,:,i)=image;
end
</code>```

This is not the fastest way off course, now it takes 332 seconds to create this array. I have found several solutions to improve this and I will do it.

Second step is to divide for a flat field (image taken without the sample). My code takes the flat and using imdivide divides each image in the array for the flat image:

```<code>for i=1:size(imArray, 3);
    imArray(:,:,i)=imdivide(imArray(:,:,i), flat);
end
</code>```

this step seems to be fast however it is called 4000 times. Do you have suggestion? Is there a better way to perform it?

Now my greatest problem is how to take in the fastest way the horizontal cut in the projection volume? My basic idea to read the volume from the "top". is reported below, it does its job, however it takes long time:

```<code>for i=1:size(sinogram, 3)
        for k=1:size(sinogram, 1)
            sinogram(k,:,i)=imArray(i,:,k);

        end
end
</code>```

Can you help me to speed up also this operation? I hope my question is clear, otherwise please ask and I will try to explain it better.


问:
scipy.linalg.sparse.eigsh does not work for generalised eigenvalues
I asked this question over at StackOverflow and someone told me that I'd get a better answer here. So here's my problem:

I'm working on a machine learning project which involves doing a Principal Component Analysis on some labeled data and using those labels to extract more valuable information from the data.

To do that, I'm calculating a scatter matrix for each class, and for each pair of classes I need to solve a generalised eigenvalue problem for their scatter matrices, as follows:

$$S_i  v = w  (S_j + \beta I) v,$$

where $\beta$ is a multiplier and $I$ is the identity matrix. Now, this is the code in python:

```<code>        jeigenvalues = eigsh(scatter_j, k=10, return_eigenvectors=False, maxiter=100)
        print('eigenvalues made')
        beta = betaMult*mean(jeigenvalues)
        print(beta)
        print(scatter_j+beta*eye(shape(x_data)[1]))
        w, v = eigsh(scatter_i,M=scatter_j+beta*eye(shape(x_data)[1]),k=int(numberOfEVs/45), maxiter=100)
        print(i,j,'done')
</code>```

numberOfEVs is 90 in my current code (so that it's divisible by 45).

But the problem is, at the line where I use the eigsh for the aforementioned formula, it never gives me an answer. It keeps eating more and more memory without even completing a single iteration (I set its maxiter input to 1, and it still didn't give an answer). When I don't give the eigsh function the M argument (which is the matrix on the right side of the generalised EV problem and it is assumed to be "I" when not specified), it works correctly. But when M is provided, it becomes unresponsive.

Any ideas?

EDIT: The scatter matrices have rather small entries, mostly around 10^-5. I've also tried multiplying the left hand side by the inverse of the RHS matrix, and again it's having the same issue (goes on for a long time without an answer). Is the smallness of these entries the issue? How can I solve it, then?


问:
Implementation of a direct solver in Fortran 90?
My question may be elementary, but it is quite essential as I am getting confused.
Here I am supposed to solve the following equation:

$Ax=B$

From my understanding I have options of using either direct solvers or iterative solvers. I have tried using iterative solver GMRES(It is quite simple in terms of implementation on Fortran 90) but now I would like to use a direct solver for certain reasons and from resources and reviews I have shortlisted three of them:

UMFPACK

SUPERLU

MUMPS

I think the first two are developed using C and could be called from Fortran 90; The third one is in Fortran 90.

Now, my question is how do I implement this or call function from my Fortran 90 code, ie. is there a reliable procedure to install or call the libraries?

There looks to be tons of literature beyond my scope of understanding.

To simplify, I have a sparse symmetric matrix $A$ , stored in skyline format and rhs vector $B$, how do I calculate $x$ using one of the direct sparse solvers in Fortran 90?


问:
3D Diffusion Equation in Fourier space
I'm solving the 3D Diffusion equation
$$u_t=k(u_{xx}+u_{yy}+u_{zz})$$

in MATLAB using Fourier techniques.  I assume a 3D Fourier expansion $(e^{-ipx},e^{-imy},e^{-imz})$of the solution.

Physical space: $u(x,y,z,t)$.
Fourier Space:  $c(m,n,p,t)$.

Substitution and differentiation result in:
 $$c(m,n,p)^{N+1}-c(m,n,p)^{N} = -\frac{k\Delta T}{2} (p^2+m^2+n^2)(c(m,n,p)^{N+1}+c(m,n,p)^{N})$$
after applying a Crank-Nicolson scheme.

I'm using fftn() and ifftn() to forward my coefficients in time and bring them back to physical space.  However I achieve universal decay from the initial condition to zero with no heat flux in any direction, for all time.  Typical time step: 0.00001.  Typical k=0.005.

Is the problem with my application of fftn() or the stability of my finite-difference?

Edit:  I've taken an initial condition of $50 \sin(2x)$. I merely took the boundaries of that initial condition and imposed them as the boundary conditions for all time. No good I'm guessing?

Edit 2: I forgot the wavenumber squared on the rhs.  Thanks James!

Title: 3D Diffusion Equation in Fourier space


问:
Calculate fluid speed in a conduit with obstacles
I need some advices to start to solve some problems about fluid motion in a conduit.

The problem domain is shown below: a conduit with some obstacle of different geometric properties.

The fluid is a "real" and viscous, the viscosity is very high in relation to speed (low speed, very high viscosity). The fluid is in pressure and comes from the blue arrow, I need to know the "speed or amount" of fluid in each output (red arrows). I don't need a real-time simulation.

How can I approach this kind of problem?


问:
Difference of convex functions optimization problem in R
I am seeking of any already written R package which could help in an optimization technique which is called  Difference of convex functions. This technique is sketched here and could be very useful for reinforcement learning processes.

Does anyone know something, that could be helpful?
I have found package called cccp, but it is all about conic optimization and I don't understand now how could it be used for such a problem.

I would like to solve the following (a bit simplified) minimization problem, which decomposes in difference of two convex functions:

$$\min_\mathbf{w} S(\mathbf{w}) =\min_\mathbf{w} \sum_{i=1}^{n}  \min ( \dfrac{|X_i\mathbf{w}|}{\phi}, 1) =  \min_\mathbf{w}~ (S_1(\mathbf{w}) - S_2(\mathbf{w}))$$

$$S_1(\mathbf{w}) =  \sum_{i=1}^{n}  \dfrac{|X_i\mathbf{w}|}{\phi} $$

$$S_2(\mathbf{w}) =  \sum_{i=1}^{n}    \left( \dfrac{|X_i\mathbf{w}|}{\phi} - 1   \right)^+$$

where $(x)^+ = \max(x, 0)$

If you have any suggestions or useful related links with any examples, they would be very appreciated.

PS. I have found that difference of convex functions optimization problem is closely related with "concave-convex optimization procedure" (CCCP). If anyone has thoughts of how the aforementioned problem could be solved by CCCP, please share!


问:
How to obtain the reduced model from a subspace projection method?
I have a system of ordinary differential equations (ODEs). It is a large system that has dozens of equations and hundreds of parameters. I wish to reduce its size so it becomes computationally more efficient to simulate, and I'm willing to tolerate some error in the reduced model. So the end result should be a system of ODEs that has a smaller number of equations or parameters.

I can supposedly use methods of subspace projection to obtain a reduced model. For example Principal component analysis / Proper orthogonal decomposition (POD) with the method of snapshots. With this method a basis is obtained.

This leads to the following POD algorithm:

Input: the data in the matrix $W$ consisting of the snapshots obtianed by simulating the full systems.

Calculate the singular value decomposition so that:
$$W = UVD' .$$

Output: The first $l$ columns of $U$ form the reduced subspace and the basis is then $\Phi = [u_1, ... , u_l]$.

If the coefficients of my linear ODE system are in $A$, I can use $\Phi' A \Phi$ to project the system to a reduced subspace and simulate it there.

Having obtained the simulation results $Y$ in the reduced space, I can do $\Phi Y$ to project the results back to the original space.

Now if $A$ models a biological system with meaningful equations, the same cannot be said for $\Phi' A \Phi$. But is there any way $\Phi$ can be used to remove equations or variables from $A$ to obtain a simpler 'real world' system?

Edit: I'm not very capable of reading math, so if possible, please explain like I am five.


问:
3D plot from Fourier coefficients in Matlab
I would like to construct a surface on a 2D mesh, given a set of Fourier coefficients for the original function $V(x,y)$.

Given some $\mathbb{K}\subset \mathbb{Z}^2$ I have a partial double sum of a Fourier series for $V(x,y)$, which could be some potential, for example

$V(x,y)=\frac{1}{\sqrt{x^2+y^2}}$

then

$f(x,y)=\sum_{(j,k)\in \mathbb{K}} v_{j,k} \exp(2 \pi i \cdot (jx+ky))$.

I am going to use Matlab for this to begin with so my question:

How can I generate a finite set of Fourier coefficients given $V$ in Matlab?

Once I have these coefficients what would the best method to plot them in 3D/contour plot?


问:
Error norm in maple
I have the function $cos(x+t)$; the following PDEs have $cos(x+t)$ as the analytical solution of them:

$$u_t = u_{xx} -\sin(x+t) + \cos(x+t) \\
u_x(0,t)=-\sin(t) \\
u_x(1,t)=-u(1,t)^4-\sin(1+t)+\cos^4(1+t)\\
u(x,0) = \cos(x)$$

$$ v_t = v_{xx} -\sin(x+t)+\cos(x+t) \\
v_x(0,t)=-\sin(t)\\
v_x(1,t)=-\sin(1+t) \\
v(x,0) = \cos(x)$$

Now I computed the error norm $\| E \| = \sqrt{\sum_{i=0}^{1/h} h\cdot | u_{numerical}(i*h)-\cos(i*h+t_{final})|^2 }$ for a set of values of spacestep $h$, the set is: $\{ 1/16 , 1/32 , 1/64 , 1/128 , 1/256 \}$.

Now the plots of $\log_{10} \| E \| $ vs $\log_{10} h$ have been plotted, and I also calcualted the slope of this linear squares fit, for the second PDE above the slope is approximately $2$ for $t_{final} =5$, but for the first PDE above the slope is approximately $0.15$, how did this difference came about?

I am adding the code below in maple 2015.

```<code>hvals:=[1/16, 1/32, 1/64, 1/128, 1/256]:
tval:=5:
dopd:= (x, y, hv, tv)->  rhs
                           ( pdsolve
                             ( x,
                               y,
                               numeric,
                               time = t,
                               range = 0 .. 1,
                               spacestep = hv,
                               timestep = hv,
                               errorest=true
                             ):-value
                                ( t=tv,
                                  output=listprocedure
                                )[3]
                           ):
  pdes:=[ diff(u(x, t), t) = diff(u(x, t), x, x)-sin(x+t)+cos(x+t),
          diff(v(x, t), t) = diff(v(x, t), x, x)-sin(x+t)+cos(x+t)
        ]
  bcs:= [ [ D[1](u)(0,t)=-sin(t),
            D[1](u)(1,t)=-sin(1+t),
            u(x,0)=cos(x)
          ],
          [ D[1](v)(0,t)=-sin(t),
            D[1](v)(1,t)=-v(1,t)^4-sin(1+t)+(cos(1+t))^4,
            v(x,0)=cos(x)
          ]
        ]:
  for p from 1 by 1 to numelems(hvals) do
      h:= hvals[p]:
      for k from 1 by 1 to 2 do
          sol[k]:= dopd( pdes[k], bcs[k], h, tval);
    od:
      E[p]:= sqrt
             ( h*add
                 ( (sol[2](j*h)-cos(j*h+tval))^2,
                   j=0..1/h
                 )
             ):
  od:
  printf( "\n\n      Step Value     Norm Value\n\n");
  seq( printf( " %12a,   %12a\n",
               hvals[j], E[j]
              ),
       j = 1..numelems(hvals)
     );
  pA:=Array(1..2):
  pA[1]:= plot( [ seq( [ hvals[j], E[j] ],
                       j=1..numelems(hvals)
                     )
                ],
                labels=["space\\timestep", "Norm"],
                labelfont=[times, bold, 16],
                labeldirections=[horizontal, vertical],
                title="Norm vs step (lin-lin)",
                titlefont=[times, bold, 20]
              ):
  pA[2]:= plot( [ seq( [ log10(hvals[j]), log10(E[j]) ],
                       j=1..numelems(hvals)
                     )
                ],
                labels=["space\\timestep", "Norm"],
                labelfont=[times, bold, 16],
                labeldirections=[horizontal, vertical],
                title="Norm vs step (log10-log10)",
                titlefont=[times, bold, 20]
             );
  plots:-display(pA);
with(CurveFitting):
LeastSquares([[log10(hvals[1]),log10(E[1])],[log10(hvals[2]),log10(E[2])],[log10(hvals[3]),log10(E[3])],[log10(hvals[4]),log10(E[4])],[log10(hvals[5]),log10(E[5])]],v);
</code>```

By the comments that I got I am adding a txt file with the computed errors, a plot of the computed solution of the PDE by maple and an explanation of the default meothd being used by maple:

"The first mode of operation uses the default method, which is a centered implicit scheme, and is capable of finding solutions for higher order PDE or PDE systems. The PDE systems must be sufficiently close to a standard form for the method to find the numerical solution.
The default method uses a second order (in space and time) centered, implicit finite difference scheme to obtain the computed solution. The number of points in the stencil of the finite difference scheme is one greater than the order of each equation.
  Note: A first order accurate accurate boundary condition (the default) can be used for second order accurate method, and still provide a solution of second order accuracy. Use of boundary conditions with an order of accuracy that is too low degrades the accuracy of the solution, while use of boundary conditions with an order of accuracy that is too high may negatively impact the stability of the problem. For first and second order accurate methods, it is generally recommended to use one or two more boundary points than the differential order of the PDE in the space variable. For third and fourth order accurate methods it is generally recommended to use two or three more boundary points than the differential order of the PDE in the space variable.
  PDE that are greater than first order in time are solved by converting to an equivalent first order system. Note: The conversion to first order is the expected one, namely that for any dependent variables that are higher than first order, variables representing all lower order time derivatives, and corresponding equations, are appended to the system. This choice of additional variables may not always be desirable, particularly when error control is in use, so often better performance can be obtained by choosing the first order system yourself. For more information on this topic, see pdsolve/numeric/firstorder.
Requested values that do not lie on a mesh point are computed using a Lagrange interpolation that uses 4 space points and 3 time points. The time points are chosen (whenever possible) so that only one time value is larger than the requested time.
"

In the next file there's the txt file with the errors.
txt file with errors computed

I hope you can help me Kirill or anyone else.

Thanks in advance.

P.S
I am adding the plot of the absolute value between the computed solution and the actual solution.

Title: Error norm in maple


问:
Computing eigenvectors from the QR algorithm
I've seen a few other posts  on this topic but none have full answers.

I'm trying to implement some eigen-decomposition algorithms. I've managed to get the Explicit QR algorithm and the Implicit (Francis) algorithm to compute the eigenvalues. My question however relates to computing the eigenvectors.

My core question is:

Do I need to use the Inverse Iteration Algorithm to get the eigenvectors after completing the QR algorithm?

I understand additionally that I'll need to use the transformation matrices from my QR algorithm (balancing, the hessenberg reduction and the QR decomposition). But I'm more interested right now in whether there is some way to get the eigenvectors from the QR algorithm itself (as wikipedia implies there is).


问:
How to find active region in treepartition Matlab object
I am trying to identify a non-linear plant using the identification toolbox of Matlab.

In particular I would like to identify a non-linear arx model, where the nonlinearities are expressed by means of tree partitions.

According to Mathworks (http://it.mathworks.com/help/ident/ref/treepartition.html):

TREEPARTITION is a nonlinear function y = F(x), where y is scalar and
x a 1-by-m vector.                      F is a piecewise linear
(affine) function of x:

F(x) = x*L+ [1,x]*C_a + d when x belongs to D_a, where

L is a 1-by-m vector, C_k is a 1-by-(m+1) vector,

and D_k is partition of the x-space.                                                                   The active partition D_a is determined as an intersection of
half-spaces by a binary tree as follows:         first a tree with N
nodes and J levels is initialized. A node at level J is a  terminating
leaf               and a node at a level j < J has two descendants at
level j+1. All levels are complete, so N = 2^(J+1)-1.        The
partition at node r is based on [1,x]*B_r > 0 or <= 0 (move to left or
right descendant),                 where B_r is chosen to improve the
stability of least-square computation on the partitions at the
descendant  nodes. Then at each node r the coefficients C_r of best
linear approximation of unknown regression function   on D_r are
computed using penalized least-squares algorithm.

[...]

When the value of the mapping F, defined by the treepartition object, is computed at x, an adaptive algorithm selects the active node k of the tree on the branch of partitions which contain x.

In this description it is not clear what is the underlying algorithm used to identify the current active node. Indeed, at a first glance, I tought that all the intermediate nodes were decisional nodes and the terminating leafs were containing the unknown coefficients that have to be used for a particular PWA representation.

That approach is not correct, indeed as bolded in the citation, every single node has its own set or parameters that best fit the unknown regression function in D_r.

Is there a clearer reference and/or explanation about this?

Thanks


问:
PDE discretization (via finite difference sheme) question
So after posting this question and reading all your comments I would like to make this new question (update).

If you consider the three equations presented here:
$$\frac{\partial \rho}{\partial t}
+\frac{\partial (\rho u)}{\partial z}=0$$
$$\frac{\partial C_i}{\partial t}
+\frac{\partial (u C_i)}{\partial z}=Reactionrate$$
$$\rho=MW_{mean}\frac{P}{RgT}, MW_{mean}=\sum^{N_{comp}}_{i=1}(Y_i*mw_i), Y_i=\frac{C_i}{\sum^{N_{comp}}_{i=1}{C_i}}$$
continuity equation, mass balance for each species & an equation to calculate density, what would the discretized form (in $z$ direction, with backwards finite differences) be? I mean, if we consider (for solution purposes) that the variables are $C_i$, $u$ & $\rho$ (Case_1 a, b & c) then it is totally different than considering that the variables are $C_i$, $(uC_i)$, $\rho$ & $(\rho u)$ (Case_2)(as mentioned by Wolfgang Bangerth at the previous question, the term $\frac{\partial (\rho u)}{\partial z}$ implies conservation of quantity $(\rho u)$ in direction $z$, so should we treat this quantity like one?). Next I present the discretized form of both this cases (please correct me if I am wrong, or more precisely, tell me which is the correct one):

Case_1a

$$\frac{\partial \rho_j}{\partial t}
+\frac{(\rho_j u_j - \rho_{j-1} u_{j-1})}{\Delta z}=0$$
$$\frac{\partial C_{i,j}}{\partial t}
+\frac{(u_j C_{i,j} -u_{j-1} C_{i,j-1})}{\Delta z}=Reactionrate$$
$$\rho_j=MW_{mean,j}\frac{P}{RgT}, MW_{mean,j}=\sum^{N_{comp}}_{i=1}(Y_{i,j}*mw_i), Y_{i,j}=\frac{C_{i,j}}{\sum^{N_{comp}}_{i=1}{C_{i,j}}}$$

Case_1b

$$\frac{\partial \rho_j}{\partial t}
+\rho_j\frac{(u_j - u_{j-1})}{\Delta z}
+u_j\frac{(\rho_j - \rho_{j-1})}{\Delta z}=0$$
$$\frac{\partial C_{i,j}}{\partial t}
+u_j\frac{(C_{i,j} - C_{i,j-1})}{\Delta z}
+C_{i,j}\frac{(u_j - u_{j-1})}{\Delta z}=Reactionrate$$
$$\rho_j=MW_{mean,j}\frac{P}{RgT}, MW_{mean,j}=\sum^{N_{comp}}_{i=1}(Y_{i,j}*mw_i), Y_{i,j}=\frac{C_{i,j}}{\sum^{N_{comp}}_{i=1}{C_{i,j}}}$$

Case_2

$$\frac{\partial \rho_j}{\partial t}
+\frac{(\rho u)_j - (\rho u)_{j-1}}{\Delta z}=0$$

$$\frac{\partial C_{i,j}}{\partial t}
+\frac{(u C_i)_j - (u C_i)_{j-1}}{\Delta z}=Reactionrate$$

$$\rho_j=MW_{mean,j}\frac{P}{RgT}, MW_{mean,j}=\sum^{N_{comp}}_{i=1}(Y_{i,j}*mw_i), Y_{i,j}=\frac{C_{i,j}}{\sum^{N_{comp}}_{i=1}{C_{i,j}}}$$

First of all, should there be any difference detween the results in Case_1a & Case_1b given the same initial $(t=0)$ and boundary $(z=0)$ conditions? What is the correct formulation (mathematically speaking)? Do you think that this set of equations has one unique solution? Because I am almost sure that it does not. What assumption could justify a change like this:

Case_1c
$$\frac{\partial \rho_j}{\partial t}
+\rho_j\frac{(u_j - u_{j-1})}{\Delta z}=0$$
$$\frac{\partial C_{i,j}}{\partial t}
+u_j\frac{(C_{i,j} - C_{i,j-1})}{\Delta z}=Reactionrate$$
$$\rho_j=MW_{mean,j}\frac{P}{RgT}, MW_{mean,j}=\sum^{N_{comp}}_{i=1}(Y_{i,j}*mw_i), Y_{i,j}=\frac{C_{i,j}}{\sum^{N_{comp}}_{i=1}{C_{i,j}}}$$
which I believe that does have one uniquie solution?

And referring to the Case_2, this way we have four(*) unknown variables ($C_i$, $(uC_i)$, $\rho$ & $(\rho u)$) which means that we need one more equation? Would an equation like:
$$\frac{(uC_i)_j}{C_{i,j}}=\frac{(\rho u)_j}{\rho_j}$$
work? and with given initial and boundary conditions provide one and only one solution to this problem?

*actually the variables are more than four. If $N_c=Number of Components$ then the variables are $N_c*C_i+N_c*(uC_i)+\rho+(\rho u)$ times the number of discrete points in z direction but the use of the last equation should give the same number of equations and variables along with the boundary conditions.


问:
ADMM for Linear Program over graph
I want to use ADMM to solve a LP defined over a graph.
According to

Distributed optimization and statistical learning via the alternating direction method of multipliers
  S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, 2011

the first primal update ($x^{k+1}$) is given by the solution of

$$ \begin{bmatrix} \rho I & A' \\ A & 0 \end{bmatrix}  \begin{bmatrix} x^{k+1} \\ \nu \end{bmatrix} + \begin{bmatrix} q - \rho(z^k - u^k) \\ -b \end{bmatrix} = 0$$

Thus, to solve it for $x^{k+1}$ one has to invert the matrix and therefore needs the full vectors $z^k$ and $u^k$.

Because my A is (almost) the incidence matrix, I would like to compute $x^{k+1}_i$ only as a function of $z^k_j$ and $u^k_j$, for $j$ such that $A_{ij} \neq 0$ (in this way, hypothetically, nodes' graph needs fewer communications).

Is there a direct way to do so?

Title: ADMM for Linear Program over graph


问:
Choosing suitable polynomial degree based on information in advection stencil
I'm working on a finite volume advection scheme for unstructured meshes which uses a multidimensional polynomial weighted least squares fit for interpolating from cell centres onto faces.

In 2D, the advection scheme typically uses a 4x3 stencil of cell values in the domain interior to fit a polynomial that is cubic in x and quadratic in y.  The weights resulting from the SVD are expensive to compute but they can be precomputed during initialisation.

There are two cases where which wish to omit some higher order terms:
1. Near boundaries (which may be irregular in shape), the stencil is necessarily smaller
2. If some cell values are very closely spaced on certain meshes, we do not have enough information to sensibly fit high order terms

I'm looking for a technique for dealing with these cases.  Can the SVD help me determine which terms to omit?  Or are there other techniques I could consider?


问:
Solving system of equations with zeros on diagonal
I am using finite elements, Newton Raphson technique and MUMPS direct solver, to solve for $P$ in this equation:

$$\frac{\delta K}{\delta x} \frac{\delta P}{\delta x}= \frac{\delta K B}{\delta x} + A$$

$A$ and $B$ are constant parameters. $K$ is a parameter varying over space. $P$ is the solution variable.

$K$ is distributed like this:

Thus, $K$ becomes zero! So I'll have zeros on the diagonal of my matrix and the solver will eventually print out that the matrix is singular.

Are there any techniques to still solve this problem?


问:
Solving compressible inviscid Euler equations with shockwaves in polar coordinates
For the past several weeks I was attempting to adapt Lax-Wendroff or some similar scheme for polar coordinates. The process was complicated due to me being unable to find step-by-step derivations of used schemes, so i had to guess the process for cartesian coordinates and to repeat the same for polar ones.

I managed to do it for MacCormack's scheme, and it exploded even on somewhat smooth, shockwave-less test.

Due to geometry of my system and because of other equations that govern it (I'm trying to add hydrodynamics into 2d plasm discharge model), I was trying to stick to polar coordinates. Hovewer, now I'm quite certain that I won't be able to find ready-for-use scheme for euler equations in polar coordinates. To derive one on my own seems implausible, too.

Thus my questions arise.

Are there ready-to-use robust solvers for inviscid compressible Euler equations with sharp gradients in polar coordinates?

If there are none, how do I attempt to derive one? Which scheme should I try to derive?

Is it even feasible to try to stick to polar coordinates? Or should I solve Euler equations in cartesian coordinates, and then interpolate values for the rest of my model? If so, should I use uniform grid with cropped cells on the border, or should I use some finite-volume method?

If I should move away from polar coordinates, which scheme should I use? There are quite a lot of schemes, there are different types of error correction mechanisms, etc - but can someone with experience suggest something?


问:
Symplectic integration of PDE
I consider ordinary wave equation
$$ u_{tt} - u_{xx} = 0 $$
with initial conditions
$$u(x, 0) = \exp (-2 x^2) \\
u_t(x, 0) = 0 $$

To solve this problem I approximate $u_{xx}$ with 4-th order difference scheme and find time evolution with symplectic 4-th order method.

I use following parameters:
$$x_{min} = -500 \\
x_{max} = 500 \\
dx = 0.5 \\
dt = 0.1 $$

The problem is that, for example, in the point $x=0$ this algorithm perfectly works for $t < 1000$ (note that this value equals $x_{max} - x_{min} = 2 x_{max}$ But then evolution looks quite strange. Here is a picture of this.

I tried to fix boundary value manually (at every step set $u(x_{min}, t) = u(x_{max}, t) = 0, \ u_t (x_{min}, t) = u_t (x_{max}, t) = 0 $ )
and also tried to use "exponential suppression" to avoid radiation from boundaries.
The $u_{xx}$ approximation with finite difference is used only for points which are more than 2 steps distant from boundaries to avoid problems.

It seems that the problem lies in the presence of boundaries, because such symplectic algorithm works fine for harmonic oscillator (where there are no boundaries) for any time value.

UPD: The stencil I use for partial derivative:
$$u_{xx} (x) = \frac{-u(x-2) + 16 u(x-1) - 30 u(x) + 16 u(x+1) - u(x+2)}{12 (dx)^2} $$

Title: Symplectic integration of PDE


问:
Does there exist a Fourier transform algorithm for perturbed data?
Assuming I have a length-$n$ real vector $x$ and have already computed its Fourier transform $\hat x$ (in time $O(n\log n)$), I would like to compute the Fourier transform of $y = x + \delta x$, where $\delta x$ is a small perturbation.

Does there exist a method for computing the Fourier transform of $y$ in less than $O(n\log n)$ time? Perhaps in linear time (e.g., asymptotically constant number of vector arithmetic operations)? In other words, is it possible to take advantage of knowing an "approximate" Fourier transform ($\hat x$) of the input data ($x+\delta x$)?

I have in mind the example where $x$ and $y$ are consecutive timesteps in solving a system of ODEs, where $|y-x|=O(\delta t)$ is small, but there might be other examples I don't know of.


问:
Clustering by shared subsequence
I have a question that relates to the classical "longest common subsequence" problem. I'll give the background to the problem, but you could skip to the formulation below if you like

Let's think of genes on a bacterial chromosome as beads on a string - atomic, unique objects, with an ordering.
I'm deliberately ignoring some things, and want to look at a simplified model. So let's say the chromosome is the interval
[1,1000], and there are 1000 genes, one at each integer in the interval.

Now suppose I have a set of 100 chromosomes, each from a different individual in a population, each with their 1000 genes,
but the ordering is allowed to differ.

In other words, I have 100 different permutations of the integers 1..1000.
Let's look at an example

1 2 3 4 5 6 7 8 9 10

1 8 9 2 3 4 5 6 7 10

1 4 3 2 7 8 9 10 5 6

1 4 3 2 10 8 9 5 6 7

When I look at these, I see two things. First, that the first two both contain
1 2 3 4 5 6 7 10
as a subsequence.
Second, that the final two contain
1 4 3 2 8 9 5 6

I want to construct  an algorithm for splitting a list of permutations into subsets like that,
such that within the subsets, all permutations contain a common subsequence.
A trivial solution is to make each permutation a subset of size 1 of course, so what I really want
is to minimize the number of subsets, and maximize the length of the subsequences.

There is a lot of previous work on finding the longest common subsequence of a fixed set (https://en.wikipedia.org/wiki/Longest_common_subsequence_problem) but I can't find any works addressing how to partition your strings to maximise the length of common subsequence within subsets.

Does anyone know if this is a classical/solved problem also?

PS I posted a (slightly less well formulated) version of this on mathOverflow and the organisers sent me over here.

Title: Clustering by shared subsequence


问:
My discretization of a wave equation in first-order form does not give correct solutions. What should I do?
I haven't much experience with conservation laws, shocks, etc. After reformulating my wave equation to 1st order system (velocity-stress):
$$
\frac{\partial v}{\partial t} + A \frac{\partial v}{\partial x} + B \frac{\partial v}{\partial y} = 0
$$
where: $v=(v_{x},v_{y},\sigma_{x},\sigma_{y},\sigma_{xy}) \\ A = A(\rho, \mu, \lambda)\\ B = B(\rho, \mu, \lambda)$

I put my Dirichlet BC to the first component (some constant vx) and obtained a stiff behavior. I've tried the same models with 2nd order equation and saw beautiful waves propagating. Here I get immediate propagation and loads of instabilities with explicit solvers (FD, RK4). The finer the mesh, the more high-frequency oscillations I get. Maybe I'm used to displacements and don't understand this formulation...

My question is:
Is there any way of solving this kind of problem with ordinary tools or should I bother with Riemann solvers, etc. ? My long-term plan is to extend it to DG. There's plenty of papers on Godunov scheme, etc. but I couldn't quite catch what's going on though...

EDIT:
Yes, it's a linear hyperbolic system. I used the term "stiff" in a more 'colloquial' manner. I can't find a stable explicit timestepping scheme for this system with Dirichlet BC. Even for extremely low dt (1e-20), way below CFL I get immediate waves in the whole domain. On the other hand, for implicit, the solution does not propagate correctly.
For instance, sinusoidal excitation u1 = y = sin(100*t) with 2nd order system looks like this (mesh is obviously too coarse but quite fine):

For 1st order system cosinusoidal excitation v1 = y' = 100*cos(100*t) like this:

The same mesh size, the same timestepping (Crank-Nicolson). I suppose this is not physically correct...

EDIT: The wave does not propagate from the middle too:


问:
Complexity of direct solvers?

Possible Duplicate:

  How to reorder variables to produce a banded matrix of minimum bandwidth?

What is the time and space complexity of direct sparse solvers (e.g., UMFPACK, SUPERLU, PARDISO, etc.) for finite element matrices in 2D and 3D? How does the complexity splits into the creation of the inverse and the forward/backward substitution? Can anybody make some notes on this topic and/or post here some references?

Title: Complexity of direct solvers?


问:
COMSOL - Implementing Perfectly Matched Layers (3D)
If there are any COMSOL gurus here, I want to simulate a 3D RF system infinite in 1 direction. I realize I need to use PML's, but am unsure how to configure them. it seems that simply scaling them by lambda isn't enough - anyone know if I need to configure the material itself in a special way?
Thanks!

EDIT: I want to solve a problem containing a infinite (in +z direction) cylindrical waveguide filled with vacuum with PEC walls. How can I replace (some of) the infinite waveguide with a PML?


问:
Solve for $D$ in $R^{T}DSDR = Id$
Given that $R$ is a rectangular matrix, $D$ is a diagonal, square matrix and $S$ being a square matrix along with the fact that both $D$ and $S$ are invertible. $S$ in this specific case can be assumed to be symmetric, positive-definite.

$R$ is actually the first few columns of the discrete Fourier Transform matrix. Thus, $R^{T}R=Id $. Note $RR^{T}$ need not be $Id$. Also, $R$ has both dimensions of the same order (but with more rows than columns). $R$ also has a full row rank.

I would like to obtain $D$ in terms of the given matrices, which are $R$ and $S$

Title: Solve for $D$ in $R^{T}DSDR = Id$