问:
minimization a quadratic form with linear constraints (Prospective method SMIC 74 method)
my problem is a ecuation i don't understand, I have no idea how to solve that, is a classic minimisation programme of quadratic form with linear constraints, here is the ecuation:

$$\sum_{i,j}^n[P(i/j)P(j) - \sum_{k=1}^rt(ijk)\prod_k]^2 + \sum_{i,j}^n[P(i/ \bar j)P(\bar j) - \sum_{k=1}^rs(ijk)\prod_k]^2$$

subject to the contrainsts:

$$\sum_{k=1}^r \prod_k =1, \prod_k \ge 0 , \text{ For all k}$$

and a example:

$\text{minimize }[(0.30 -  \prod_1, -  \prod_5)^2  + (0.15 -  \prod_3, -  \prod_7)^2
 + (0.30 -  \prod_1, -  \prod_3)^2  + (0.04 -  \prod_5, -  \prod_7)^2
 + (0.24 -  \prod_1, -  \prod_5)^2  + (0.04 -  \prod_2, -  \prod_6)^2
 + (0.24 -  \prod_1, -  \prod_2)^2 + (0.24 -  \prod_5, -  \prod_6)^2
 + (0.27 -  \prod_1, -  \prod_3)^2 + (0.35 -  \prod_2, -  \prod_4)^2
 + (0.35 -  \prod_1, -  \prod_2)^2 + (0.20 -  \prod_3, -  \prod_4)^2]$

All these equations I got from this paper called "SMIC 74—A method for constructing and ranking scenarios", i must applicate this in a programming languaje (preferably web), can help me with this forecasting method SMIC?


问:
Some questions about fractional derivative
I have read some papers about time-fractional PDEs solved by finite element method. The time-fractional derivative is the Caputo derivative defined by
$$
\frac{\partial^{\alpha}u}{\partial t^{\alpha}}=\frac{1}{\Gamma(1-\alpha)}\int_{0}^t\frac{\partial u}{\partial s}\frac{ds}{(t-s)^{\alpha}},~~~~0<\alpha<1.
$$

Now, we consider the following time-fractional diffusion equation
\begin{align}
&\frac{\partial^{\alpha}u}{\partial t^{\alpha}}-\Delta u=f,~~in~~\Omega,\\
&u=0,~~~~~~~~~on~~~\partial\Omega,\\
&u(0)=u_0,
\end{align}
where $\Omega$ is a convex and bounded domain.

1.What does mean the singularity of the time-fractional derivative?

2.In what cases, the singularity of the time-fractional derivative can lead to low smoothness of solution $u$.

3.What does mean the regularity of a function is not smooth about time?

Could anyone give me some explanations?

ContentLicense: CC BY-SA 3.0


问:
Produce equation of curve, given some coordinates

Possible Duplicate:

  Get equation for a curve which intersects x at seemingly randomly distributed points?

I tried asking something similar to this before, but I guess I didn't explain very well. Hopefully this will be more clear.

I need to get an equation which describes a curve with certain points.

For example, I have a list of 5-10 different points along the x-axis where a curve either touches or crosses the x-axis.

Now, how can I get an equation to describe any curve (at least one of the possible curves) which satisfies this?

I'm only concerned with x values within the range of the given values. So, if the curve interacts with the x-axis at other points outside this range, that would be ok too, it doesn't matter....

So to put it another way,

if y = 0; and 0 < x < 100; also x = (list of given points eg. 0,4,17,19,56,67,100)

It doesn't matter what the other x values are when y ≠ 0. And it doesn't matter what other x values might be outside the given range (in this case when x > 100).

Is there any website or free software which can take the coordinates and spit out an equation? Or is there a method you can explain, which I can use to write my own program for doing so?

To be more clear, I'm given a list of points all where y=0, How can I get an equation such that, someone else could look at the equation alone, without also seeing the points I used to produce the equation, but then be able to use that equation to get those same points.... to tell me where the line crosses x within the range I specify?

Also,
If multiple possibilities exist, the shorter length of equation the better.


问:
Spurious oscillations in diffusion-reaction problems with finite volume
I have successfully solved the multi-species diffusion-reaction equation
\begin{equation}
\frac{\partial c_i}{\partial t} = \nabla \cdot (d_i(x)\nabla c_i) + s_i(x,t), \quad \quad (1)
\end{equation}
with dicontinuous source term
\begin{equation}
 s(x,t) =  \left\{ \begin{array}{rcl} s_1(x,t) & \text{for} & 0<x\le x' , \\ s_2(x,t) &\text{for} & x'<x\le 1. \end{array}\right.
\end{equation}
In fact, the variable diffusion coefficient also takes slightly different forms in these two region. The equation is discretized with the conservative central scheme
\begin{equation}
c_j'(t) = \frac{1}{h^2}\biggl( d(x_{j+1/2})(c_{j+1}(t)-c_j(t))- d(x_{j-1/2})(c_{j}(t)-c_{j-1}(t))\biggr) + s(x_j,t).
\end{equation}
I used SUNDIALS package for time integration.

Now I'm trying to solve the problem in a  concentrated solution
\begin{equation}
\frac{\partial c_i}{\partial t} = \nabla \cdot \sum_k  D_{ik}(x)\nabla c_k + s(x,t),
\end{equation}
where $D_{ik}$ may be negative. In dilute solution approximation (solvent concentration $c_0 \rightarrow \infty$), we have $D_{ik} \rightarrow d_{i}\delta_{ik}$. However, the solution exhibit spurious oscillations.

The blue line is sloved with $c_0=10^{20}$ or equation (1); purple and green lines with $c_0=10^{10}$. The difference between purple and green line is the average scheme used for the diffusion coefficient at the cell interface, the former is with simple average and the latter with harmonic average. As shown in the figures below, there are oscillations near the boundaries and the point $x'$. (I should mention that the vary first osscilations began at the left boundary before they built up at the right boundary and the point $x'$).

It is odd that even with the oscillations, the scheme with harmonic average never fails to converge at each timestep and follows closely the result of analytical dilution scheme (equation (1)); where as the simple average scheme fails at the one tenth of the time domain. At much later time the solutions look like

Are there discretizatons within the finite volume scheme that avoid such oscillations?


问:
How to impose Neumann boundary conditions in finite volume problems?
I'm trying to better understand finite volume methods and have started coding up a basic script to solve the diffusion equation $$u_t = u_{xx}$$ which has the finite volume form: $$\frac{\bar{u}^{n+1}_i\ - \bar{u}^n_i }{\Delta t} = \hat{F}_{i+\frac{1}{2}}^n - \hat{F}_{i-\frac{1}{2}}^n$$ where $\bar{u_i}$ is the average value of $u$ over cell $i$ and $\hat{F}$ is the numerical approximation to the flux at a cell interface.

I am using an ENO approximation in order to reproduce the point values $u$ and the derivative values $u_x$ as well. I'm aware of the ENO method to reproduce the numerical flux at each cell, however, I'm confused on how I can implement no flux boundary conditions using this method. Instead, here's what I've tried:

Using my array of known averages, $\bar{u}_i^n$, construct two coefficient matrices. One that approximates the value of the flux, $u_x$, at the left cell interface (which I'll call $L$) and one approximates it at the right cell interface (which I'll call $R$), giving me an equation that looks like $$\bar{u}^{n+1}_i = \bar{u}^n_i - \frac{\Delta t}{\Delta x}(R-L)\bar{u}^n_i$$. This allows me to just make the coefficients of the first row of $L$ and the last row of $R$ all zeros in order to implement no flux boundary conditions. However, when doing it this way, I seem to get something that is obviously non-conservative after I advance one time step:

This is clearly due to my implementation of the flux conditions, as the interior nodes all give correct values within the prescribed error.

Does anyone know how to solve this problem?


问:
How To Calculate Theoretical CPU FLOPS?
I actually find the formulae for peak theoretical performance:

Node performance in GFlops = (CPU speed in GHz) x (number of CPU cores) x (CPU instruction per cycle) x (number of CPUs per node)

CPU Speed and CPU cores are easy. but how can i know the cpu operations per cycle?
For example, my qualcomm processor specs are as following:

CPU Clock Speed Up to 1.2 GHz
CPU Cores:Quad-core CPU /4x ARM Cortex A7
CPU Bit Architecture:32-bit
Can someone please help me with this.

Title: How To Calculate Theoretical CPU FLOPS?


问:
Numerical diffusion in during advection of a free surface in an FE context
I am currently working on a project where a two-phase flow is considered. The phases are described using a level set approach and a signed distance function from the interface between the phases where positive values are liquid and negative values are gas.

I consider a Stoke's flow which gives me a velocity field $v$ and a pressure field $p$. Then I use $v$ to advect the level set using the advection equation

$\frac{\partial \phi}{\partial t} + v \cdot\nabla\phi=0$.

I am new to this type of problem but I think I have some basic understanding on the numerical (diffusion) problems arising from this equation due to discretization and I know that I need some kind of stabilization. My research concerns the modelling of two-phase flow in a porous material and the results needs to be quite accurate, hence the popular SUPG scheme does not quite cut it.

My question is simply, what schemes (for Finite Elements) exists, that gives a better solution than the SUPG scheme? And furthermore, is the any good literature on the topic that discuss, and preferably, compares any stabilization schemes?


问:
How to implement chaotic sender and receiver with ordinary differential equations?
I am referring this  paper, and trying to implement chaotic sender and receiver, to decode message, as given in section $V^{th}$ Chaotic signal masking.
The process that I want to implement is figure 6 of the same paper, where we have Lorenz system given by equations:
:  \begin{align}
\frac{\mathrm{d}x}{\mathrm{d}t} &= \sigma (y - x), \\
\frac{\mathrm{d}y}{\mathrm{d}t} &= x (\rho - z) - y, \\
\frac{\mathrm{d}z}{\mathrm{d}t} &= x y - \beta z.
\end{align}  [from Wikipedia]and $\sigma=16; \rho=45.6; \beta=8$ for real applications the parameters have been scaled.

My MATLAB code is given here:

```<code>clear all;
close all;
clc;
%Solution for the Lorenz equations in the time interval [0,100] with initial conditions [1,1,1].
sigma=16;
beta=4;
rho=45.6;
f = @(t,a) [-sigma*a(1)/10 + sigma*a(2)/10; rho*a(1)/10 - a(2)/10 - (a(1)*a(3))/200; -beta*a(3)/20 + (a(1)*a(2))/100];
%'f' is the set of differential equations and 'a' is an array containing values of x,y, and z variables.
%'t' is the time variable
[t,a] = ode45(f,[0 10],[1 1 1]);%'ode45' uses adaptive Runge-Kutta method of 4th and 5th order to solve differential equations
% plot3(a(:,1),a(:,2),a(:,3))
m_t= rand(1,size(a,1))<=0.5;
s_t=a(:,1)+ m_t';
f = @(t,a_r) [-sigma*a_r(1)/10 + sigma*a(1)/10; rho*a_r(1)/10 - a(1)/10 - (a_r(1)*a(1))/200; -beta*a(1)/20 + (a_r(1)*a(1))/100];
[t,a_r] = ode45(f,[0 10],[1 1 1]);
m_t_est=s_t(1:size(a_r,1),:)-a_r(:,1);
l=m_t(:,1:size(a_r,1))'-m_t_est;
plot(l)
</code>```

part of this code has been taken from Wikipedia: Lorenz system, and modified suitably, my problem is twofold:

I am not sure as how this ODE is working as integrating from 1 to 10 the size I am getting is like size(a_r(:,1))=61. How to decide the size of given function?

The plot of difference between received waveform and transmitted is like this:

which is quite fine, for some samples but increases dramatically.

Am I doing this the right way or something is really wrong?


问:
How to sort 13 enormous polynomials, each having terms with 5 variables ($1>h_0>h_1>h_2>h_3>h_4$)?
Both algebraically or with software (numerically/computationally, for example) is acceptable. Here is a how one of the polynomials looks like ($L<1$ is a constant).

I have attached a text file with all the 13 polynomials that I need to sort in ascending order - more precisely, I have a target order I need to achieve, so I need to define the ranges of the variables $h_0$ - $h_4$ where that order would be valid. Here is one representative polynomial sample:

$$
7 h_0 h_1^2 h_2^2 h_3 L^3 +
29 h_0 h_1^2 h_2^2 h_3 h_4^2 L^4 +
3 h_0 h_2^2 h_4 L^2 +
5 h_0 h_2^2 h_4^3 L^3 +
\ldots
$$

Each polynomial has approximately 200 terms(!).


问:
4th order Symplectic integration of a 2 degree of freedom Hamiltonian system
```<code>#f,g are the dH/dq and -dH/dp respectively

def f(v1,v2):

    alpha = (a*v1 - b*v2)/k
    beta =  (a*v2 - c*v1)/k

    dot1 = (c-a)/(3*k) +  (a*np.exp(alpha)/k - c*np.exp(beta)/k)/(1 + np.exp(alpha) + np.exp(beta))
    dot2 = (b-a)/(3*k) +  (a*np.exp(beta)/k - b*np.exp(alpha)/k)/(1 + np.exp(alpha) + np.exp(beta))
    return dot1, dot2

def g(u1,u2):

    dot1 = -1/3 + np.exp(u1)/(1 + np.exp(u1)+ np.exp(u2))
    dot2 = -1/3 + np.exp(u2)/(1 + np.exp(u1)+ np.exp(u2))
    return -dot1, -dot2

#these are the constants for 4th order symplectic rk method where c1 is for canonical position and c2 for canonical momentum cordinates. (As given by Ruth and Forest)

c1 = np.array((1.0/(2.0*(2.0-2.0**(1.0/3.0))), (1.0-2.0**(1.0/3.0))/(2.0*(2.0-2.0**(1.0/3.0))),(1.0-2.0**(1.0/3.0))/(2.0*(2.0-2.0**(1.0/3.0))), 1.0/(2.0*(2.0-2.0**(1.0/3.0)))))

c2 = np.array((1.0/(2.0-2.0**(1.0/3.0)), -2.0**(1.0/3.0)/(2.0-2.0**(1.0/3.0)),1.0/(2.0-2.0**(1.0/3.0)), 0.0))

#this is the function where the integration is happening and returns the 4 arrays which are time series of the phase-space variables

def sol(x0):

    u1 = np.zeros(n+1)
    u2 = np.zeros(n+1)
    v1 = np.zeros(n+1)
    v2 = np.zeros(n+1)

    u1[0] = x0[0]
    u2[0] = x0[1]
    v1[0] = x0[2]
    v2[0] = x0[3]
    for k in xrange (n):
            u10 = u1[k]
            u20 = u2[k]
            v10 = v1[k]
            v20 = v2[k]
            for i in xrange(len(c2)):
                    U1 = u10 + c1[i]*h*f(v10,v20)[0]
                    U2 = u20 + c1[i]*h*f(v10,v20)[1]
                    V1 = v10 + c2[i]*h*g(U1,U2)[0]
                    V2 = v20 + c2[i]*h*g(U1,U2)[1]
                    u10 = U1
                    u20 = U2
                    v10 = V1
                    v20 = V2

            U1 = u10 + h*c1[-1]*f(v10,v20)[0]
            U2 = u20 + h*c1[-1]*f(v10,v20)[1]
            u10 = U1
            u20 = U2
            u1[k+1] = U1
            u2[k+1] = U2
            v1[k+1] = V1
            v2[k+1] = V2

    return u1,u2,v1,v2
</code>```

When I run the code for some values of parameters (like a,b,c, etc) I get an overflow error in dot2. For other parameters I am getting a phase space which is not bounded, but a Hamiltonian system is supposed to have a bounded trajectory.

Is something wrong with my algorithm? Kindly point the error. I can provide the whole code and the paper, whose results I am trying to reproduce.


问:
Is there an algorithm on python for analyzing swarming behavior or a modeling algorithm to predict swarming behavior?
Is there an algorithm on python to analyze swarming behavior or model swarming behavior for predictions?</p>


问:
Breather solutions of Sine-Gordon Using Finite Differences
I'm attempting to simulate a standing breather of the form
$$ u(x,t)=4\tan^{-1}\left(\sqrt{3}\cos\left(\frac{t}{2}\right)sech\left(\frac{\sqrt{3}x}{2}\right)\right)$$
for the Sine-Gordon equation
$$u_{tt}-u_{xx}+\sin u  =0 $$
on the real line by using centered differences in both time and space. Because of reflections and the use of zero Dirichlet conditions, I get something like this after just a bit of time:

How does one avoid this and make it seem like the simulation is on the whole of $\mathbb{R}$?

Code in MATLAB for reference

```<code>l=40;delta=0.08;T=80;h=0.1;
t=linspace(0,T,T/delta);
x=linspace(-l,l,l/h);
u=zeros(length(x),length(t));
u0=4*atan(sqrt(3)*sech(sqrt(3)*x/2));
ut0=0*x;
for n=2:length(x)-1
        u(n,1)=u0(n);
        u(n,2)=(delta/h)^2/2*(u0(n+1)+u0(n-1))-(delta)^2/2*sin(u0(n))+(1-(delta/h)^2)*u0(n)+delta*ut0(n);
end
for m=2:length(t)-1
    for n=2:length(x)-1
        u(n,m+1)=(delta/h)^2*(u(n+1,m)+u(n-1,m))-delta^2*sin(u(n,m))+2*(1-(delta/h)^2)*u(n,m)-u(n,m-1);
    end
    plot(x,u(:,m+1));
    ylim([-4 4]);
    xlabel('x');
    ylabel('u(x,t)');
    pause(0.001);
end
</code>```


问:
Modelling fracture flow in stress field in COMSOL
I am new to COMSOL and I want to model fluid flow in stress field. As starting point I want a simple 2D box having a fracture from right to left while pressure gradient vertically. I want to know the threshold value of that pressure from top which allows flow of fluid in fracture. Can anyone please suggest me which modules of COMSOL should I be using. I see simple fluid flow in fracture examples but not where impact of pressure controls the flow!</p>


问:
modelling surface coverage of spheres with spheres of varying size distribution
Suppose you have spheres (e.g. cells). You add particles to the suspended cells which are smaller but do have a size distribution rather than a single diameter. The particles stick to the cells. I would like to model how much of the cell surface is covered by particles and how much of the covered fraction is covered by which particle size fraction. I have no idea how to start (probably with a pen and a paper).

Are there any computer packages which can do similar calculations. Is there any Literature I havent found because I do not know what terms to look for? How would you approach such a task ?


问:
How to use "fill_factor" in spilu (scipy)
I am using spilu (incomplete LU decomposition for sparse matrix) in scipy:
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.spilu.html

However there is an option call "fill_factor" which is quite unclear to me.

If fill_factor=100, does it mean the matrix is full; and fill_factor=50 means the matrix is half-full?


问:
Is MATLAB feasible for app deployment versus lower level languages
MATLAB/Octave is a pretty useful language for prototyping scientific computation use cases, however I would like to know if the sheer nature of being a high level language bottlenecks how fast the performance can get. (I'm not talking about writing optimal code like vectorization, but a general impact on CPU Time, like how C++ can be made faster than equivalent Python programs by really large magnitudes).

My use case involves running a lot of statistical computations (like Sample Entropy for a 1D temporal time series) and I don't have equivalent code ready in C right now to test my hypothesis.

As we may agree upon as a rule of thumb, lower level languages outperform higher level languages at the cost of code readability and learning curves, but I had also heard that MATLAB interfaces super fast Fortran libraries for many of its numeric computations so I'm unsure how to proceed without unnecessarily wasting my time.


问:
Spectrum of the Laplace operator
I am studying the discretization of Poisson's equation in $1D$. In Matlab I created different discretization matrices (Laplace operator) according to different sizes of the mesh:

```<code>% Parameters
N = [40, 80, 160, 1280]; % Dimension 1-D grid

% Laplace operators
L1 = -1/((1/N(1))^2) * full( ( spdiags(repmat([1,-2,1], N(1)-1, 1), -1:1, N(1)-1, N(1)-1) ) );
L2 = -1/((1/N(2))^2) * full( ( spdiags(repmat([1,-2,1], N(2)-1, 1), -1:1, N(2)-1, N(2)-1) ) );
L3 = -1/((1/N(3))^2) * full( ( spdiags(repmat([1,-2,1], N(3)-1, 1), -1:1, N(3)-1, N(3)-1) ) );
L4 = -1/((1/N(4))^2) * full( ( spdiags(repmat([1,-2,1], N(4)-1, 1), -1:1, N(4)-1, N(4)-1) ) );
</code>```

I then plotted the eigenvalues of the matrices using

```<code>[V1, D1] = eig(L1);
plot(diag(D1) / norm(diag(D1),2), '--');
</code>```

What I see is that as I increase $N$ (so the step size goes to zero), the difference between the biggest and the smallest eigenvalues decreases. Therefore the condition number is getting better. Is there anything else I can see from such a plot?

Title: Spectrum of the Laplace operator


问:
How do we derive the elemental equation for Discontinuous Galerkin method using Centered Numerical Flux?
If we take a 1st order polynomial approximation in each cell, we can find the (2x2) -mass matrix, differentiation matrix, and flux matrix through the integration of Lagrange polynomials.

However, I am not able to comprehend the simplified equation which we obtain on using centered numerical flux. It seems that we get the result same as up-winding/ down-winding if we ensure continuity of solution.

According to me, if we multiply the element equation by its mass inverse we should be getting the final form, but it doesn't seem to be the case. I would be glad if someone points out my mistake in understanding.

I am attaching the content, which I am referring to: -


问:
IDE for Scientific computing using Fortran
I am looking to learn Fortran for Scientific Computing, which IDE setup for this exists that supports debugging, code completion capabilities and offers granular Running time Analytics. I am interested in a free software that accomplishes the above requirements and is available on Mac.</p>


问:
How to reshape matrix into row-major order for MKL DSS?
I would like to use MKL to solve a sparse linear system. I chose the DSS (Direct Sparse Solver) interface, which implements the following steps:

```<code>//(1).define the non-zero structure of the matrix
dss_define_structure(handle, sym, rowIndex, nRows, nCols, columns, nNonZeros);
//(2).reorder the matrix
dss_reorder(handle, opt, 0);
//(3).factor the matrix
dss_factor_real(handle, type, values);
//(4).get the solution vector
dss_solve_real(handle, opt, rhs, nRhs, solValues);
//(5).deallocate solver storage
dss_delete(handle, opt);
</code>```

According to my test, the DSS uses a column-major ordering. That means

```<code>//column = 3
{ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 }
</code>```

is equivalent to

```<code>{{1,  6,  11},
 {2,  7,  12},
 {3,  8,  13},
 {4,  9,  14}
 {5, 10,  15}}
</code>```

For instance, there is a sparase linear system $\mathbf A_{5\times 5}\mathbf X_{5\times 3}=\mathbf B_{5\times 3}$

where $\mathbf A$ is symmetric sparse array

```<code>//A stored with CSR3 format
rowIndex = { 0, 5, 6, 7, 8, 9 };
columns = { 0, 1, 2, 3, 4, 1, 2, 3, 4 };
values = { 9, 1.5, 6, .75, 3, 0.5, 12, .625, 16 };
//B
rhs[5*3] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15}
//X
solValues[15]
</code>```

By calling the above DSS interface, the solValues is

```<code>{-326.333, 983, 163.417, 398, 61.5,
 -844.667, 2548, 423, 1028, 159,
-1363, 4113, 682.583, 1658, 256.5}
</code>```

Namely

$$
\mathbf X_{5\times 3}=\left(
\begin{array}{ccc}
 -326.333 & -844.667 & -1363 \\
 983 & 2548 & 4113 \\
 163.417 & 423 & 682.583 \\
 398 & 1028 & 1658 \\
 61.5 & 159 & 256.5 \\
\end{array}
\right)
$$

In my application, the matrix is row-major. How to deal with this problem?


问:
Problem in analyzing the program of Gauss Jordan Inverse problem
I had to code a program which calculates Inverse of a matrix by Gauss-Jordan Inverse method , I was trying to analyse the program and then code it myself.

the link

http://hullooo.blogspot.in/2011/02/matrix-inversion-by-gauss-jordan.html

I understand it all except two things -

1) I don't understand how it does the row transformations.

2) While doing the Jordan method why the $j$ loop starts over and ends till last implying we are always pivoting that is moving the largest element up and hence I am unable to see a stable way of doing the left side part of the Augmented matrix to be an Identity matrix.

For reference to Jordan method I am using this link -

https://www.mathsisfun.com/algebra/matrix-inverse-row-operations-gauss-jordan.html

Any help is great!
Am I leaving something or interpreting in a wrong way?


问:
Reduce large sparse linear operators to memory efficient loops?
I'm dealing a lot with large sparse linear operators these days and I'm quite new to them. A lot of the matrices I deal with originate with only a few unique integers, however, there are lots of them. As a result I have tons of memory locked up for the same values and fetching them from RAM seems like quite the bottle neck. It seems like a well designed for loop could replace these matrices for matrix multiplication, but writing custom loops for that seems like a chore. Are there tools or techniques to do this?</p>


问:
Discretization of a multi-function term
I'm trying to do discretization to the following system:

$\frac{{\partial \rho }}{{\partial t}} =  - \frac{{\partial \rho }}{{\partial x}}u - \rho \frac{{\partial u}}{{\partial x}}$

$\frac{{\partial u}}{{\partial t}} =  - u\frac{{\partial u}}{{\partial x}} - \frac{\partial }{{\partial x}}\left[- {\frac{1}{{\sqrt \rho  }}\frac{{{\partial ^2}(\sqrt \rho  )}}{{\partial {x^2}}}} \right]$

with the boundary conditions:

$\rho (x,t = 0),\,\,\,\,u(x,t = 0)$

For the first equation everything is good:

$\frac{{\rho _i^{n + 1} - \rho _i^n}}{{\Delta t}} =  - \frac{{\rho _{i + 1}^n - \rho _i^n}}{{\Delta x}}u_i^n - \rho _i^n\frac{{u_{i + 1}^n - u_i^n}}{{\Delta x}}$

And I can explicitly solve for $\rho _i^{n + 1}$.

But for the second equation. i'm struggling with the term:

$\frac{\partial }{{\partial x}}\left[ { - \frac{1}{{\sqrt \rho  }}\frac{{{\partial ^2}}}{{\partial {x^2}}}\left( {\sqrt \rho  } \right)} \right]$

So my questions are as follows:

Should I use the chain rule in a straight forward manner like this:

$\frac{\partial }{{\partial x}}\left[ { - \frac{1}{{\sqrt \rho  }}\frac{{{\partial ^2}}}{{\partial {x^2}}}\left( {\sqrt \rho  } \right)} \right] = \frac{1}{\rho }\left[ { - \frac{1}{{2{\rho ^2}}}{{\left( {\frac{{\partial \rho }}{{\partial x}}} \right)}^3} + \frac{1}{\rho }\left( {\frac{{\partial \rho }}{{\partial x}}} \right)\left( {\frac{{{\partial ^2}\rho }}{{\partial {x^2}}}} \right) - \frac{1}{{2\rho }}\left( {\frac{{{\partial ^3}\rho }}{{\partial {x^3}}}} \right)} \right]$

and afterwards discretize the derivatives using forward Euler like in the first equation ?

If I'm using the approach above, how should I calculate the third derivative ?

Should I use some other clever way to discretize the term ? maybe calculating the square root of $\rho$ before hand ? I don't know if it's more efficient that way

Title: Discretization of a multi-function term


问:
Rhie--Chow interpolation on PDE level
The Rhie--Chow [1] interpolation seems to be a standard tool in the Finite-Volume discretization of incompressible flows.

It is commonly defined on the discrete level [2].

In the lecture notes [3] -- Chapter 6, page 56 [2] -- it says that

We find that the Rhie-Chow interpolation is the same as
  adding a pressure term, which is proportional to a third [spatial] derivative. [...] In the continuity equation [...] the term is proportional to a fourth-order derivative term [...]

Is there any work on how this third or fourth order term actually looks like? Or is there an approximation of the NSE on the PDE level such that a, say, 1st order FVM discretization would directly lead to the Rhie--Chow interpolated scheme?

Title: Rhie--Chow interpolation on PDE level


问:
Eigenvalues using QR iteration
I'm trying find the eigenvalues of a matrix A using QR iteration with Householder. I used this code which I found from Cornell University that decomposes QR with Householder.

```<code> A=  [3.56829648331616, -2.19026841424281,0.557474633283884,-0.126072015319561;
    -3.25393787896303, 4.81285873937612,-2.79568990543933,0.632240894103258
    1.14081451177546,   -4.13904528762030,  6.19190034744431,-2.98172499459682
    -0.639981731971678, 2.32194922534595,   -7.91889162828103,6.11802605281534]

[m,n] = size(A);
Q = eye(m); % Orthogonal transform so far
R = A; % Transformed matrix so far

for j = 1:n
% -- Find H = I-tau*w*w’ to put zeros below R(j,j)
normx = norm(R(j:end,j));
s = -sign(R(j,j));
u1 = R(j,j) - s*normx;
w = R(j:end,j)/u1;
w(1) = 1;
tau = -s*u1/normx;
% -- R := HR, Q := QH
R(j:end,:) = R(j:end,:)-(tau*w)*(w'*R(j:end,:));
Q(:,j:end) = Q(:,j:end)-(Q(:,j:end)*w)*(tau*w)';
end
</code>```

Now that I have the correct  R and Q verified through MATLAB, I should be able to find the eigenvalues by multiplying R and Q.

Let $E=RQ$,

where the diagonal elements of E must contain the eigenvalues and is symmetric.

However, upon doing so its not giving me the correct eigenvalues.

Do I have to write another iteration to find the eigenvalues with Q and R?

Title: Eigenvalues using QR iteration


问:
How to prove the relationship $\|I_h u\|_{0,\infty}\leq C$ for $u\in H^2(\Omega)$?
If $I_h$ is the bilinear interpolation operator, i.e.
$I_h: H^2(\Omega)\rightarrow V_h$, where $V_h=\{v_h\in L^2(\Omega)|v_h|_K\in \mathrm{span}\{1,x,y,xy\},~~K
\in\mathcal{T}_h\}$

i know that $\|I_hu\|_{0,\infty}\leq C(u)$ for $u\in H^2(\Omega)$. However, i dont know how to prove it. Can any one give me some suggestions?

In addition, if $I_h$ is a nonconforming interpolation operator(such as Crouzeix-Raviart element, the degree of freedom is $\int_l uds$). More precisely,
$I_h: H^1(\Omega)\rightarrow V_h$, where $V_h=\{v_h\in L^2(\Omega)|v_h|_K\in \mathrm{span}\{1,x,y\},~~K
\in\mathcal{T}_h\}$. Does the above relationship still hold? If so, how to prove it?

If the $\|\nabla I_hu\|_{0,\infty}\leq C(u)$ holds, $u$ should belong to $H^{?}(\Omega)$?

I'm very glad to receive some suggestions.


问:
Growing error from a smooth initial condition for Fisher KPP equation
I'm studying the Fisker-KPP equation on the line (and in $]0, 100[$ numerically):
$$
\partial_t u = \Delta_{xx} u + u(1-u)
$$
I notice a behavior I don't understand with a smooth initial condition $u_0$ that has the following form:
$$
u_0(x) =
\left\{
\begin{aligned}
&1 \quad \mbox{if} \quad |x-50| < 10 \\
&\exp{\left(\frac{1}{3^2} - \frac{1}{||x-50|-13|^2}\right)}
\left(1 - \exp{\left( -\frac{1}{||x-50|-10|^2}\right)}\right) \quad \mbox{if} \quad 10 < |x-50| < 13 \\
&0 \quad \mbox{if} \quad |x-50| > 13
\end{aligned}
\right.
$$
I'm using as a numerical scheme the Strang splitting, and a code taken from here:

```<code>"""
solve a scalar diffusion-reaction equation:

 phi_t = kappa phi_{xx} + (1/tau) R(phi)

using operator splitting, with implicit diffusion

M. Zingale
"""

#from __future__ import print_function

import numpy as np
from scipy import linalg
from scipy.integrate import ode
#import sys
import matplotlib.pyplot as plt

def frhs(t, phi, tau):
    """ reaction ODE righthand side """
    return 0.25*phi*(1.0 - phi)/tau

def jac(t, phi):
    return None

def react(gr, phi, tau, dt):
    """ react phi through timestep dt """

    phinew = gr.scratch_array()

    for i in range(gr.ilo, gr.ihi+1):
        r = ode(frhs,jac).set_integrator("vode", method="adams",
                                         with_jacobian=False)
        r.set_initial_value(phi[i], 0.0).set_f_params(tau)
        r.integrate(r.t+dt)
        phinew[i] = r.y[0]

    return phinew

def diffuse(gr, phi, kappa, dt):
    """ diffuse phi implicitly (C-N) through timestep dt """

    phinew = gr.scratch_array()

    alpha = kappa*dt/gr.dx**2

    # create the RHS of the matrix
    R = phi[gr.ilo:gr.ihi+1] + \
        0.5*alpha*(    phi[gr.ilo-1:gr.ihi] -
                   2.0*phi[gr.ilo  :gr.ihi+1] +
                       phi[gr.ilo+1:gr.ihi+2])

    # create the diagonal, d+1 and d-1 parts of the matrix
    d = (1.0 + alpha)*np.ones(gr.nx)
    u = -0.5*alpha*np.ones(gr.nx)
    u[0] = 0.0

    l = -0.5*alpha*np.ones(gr.nx)
    l[gr.nx-1] = 0.0

    # set the boundary conditions by changing the matrix elements

    # homogeneous neumann
    d[0] = 1.0 + 0.5*alpha
    d[gr.nx-1] = 1.0 + 0.5*alpha

    # dirichlet
    #d[0] = 1.0 + 1.5*alpha
    #R[0] += alpha*0.0

    #d[gr.nx-1] = 1.0 + 1.5*alpha
    #R[gr.nx-1] += alpha*0.0

    # solve
    A = np.matrix([u,d,l])
    phinew[gr.ilo:gr.ihi+1] = linalg.solve_banded((1,1), A, R)

    return phinew

def est_dt(gr, kappa, tau):
    """ estimate the timestep """

    # use the proported flame speed
    s = np.sqrt(kappa/tau)
    dt = gr.dx/s
    return dt

class Grid(object):

    def __init__(self, nx, ng=1, xmin=0.0, xmax=1.0, vars=None):
        """ grid class initialization """

        self.nx = nx
        self.ng = ng

        self.xmin = xmin
        self.xmax = xmax

        self.dx = (xmax - xmin)/nx
        self.x = (np.arange(nx+2*ng) + 0.5 - ng)*self.dx + xmin

        self.ilo = ng
        self.ihi = ng+nx-1

        self.data = {}

        for v in vars:
            self.data[v] = np.zeros((2*ng+nx), dtype=np.float64)

    def fillBC(self, var):

        if not var in self.data.keys():
            sys.exit("invalid variable")

        vp = self.data[var]

        # Neumann BCs
        vp[0:self.ilo+1] = vp[self.ilo]
        vp[self.ihi+1:] = vp[self.ihi]

    def scratch_array(self):
        return np.zeros((2*self.ng+self.nx), dtype=np.float64)

    def initialize(self):
        """ initial condition """

        phi = self.data["phi"]
        length1 = 10.
        length2 = 13.
        epsilon   = length2 - length1
        phi[:] = np.maximum( \
                 np.exp( 1./epsilon**2 - 1./(np.abs(np.abs(self.x-50.)-length2))**2) * \
                 ( 1. - np.exp( -1./(np.abs(np.abs(self.x-50.)-length1))**2)) * \
                 ( self.x >  50.-length2 ) * \
                 ( self.x <  50.+length2 ) \
                 , \
                 ( self.x >= 50.-length1 ) * \
                 ( self.x <= 50.+length1 ) \
                 )

def interpolate(x, phi, phipt):
    """ find the x position corresponding to phipt """

    idx = (np.where(phi >= 0.5))[0][0]
    xs   = np.array([x[idx-1],   x[idx],   x[idx+1]])
    phis = np.array([phi[idx-1], phi[idx], phi[idx+1]])

    xpos = 0.0

    for m in range(len(phis)):
        # create Lagrange basis polynomial for point m
        l = None
        n = 0
        for n in range(len(phis)):
            if n == m:
                continue

            if l == None:
                l = (phipt - phis[n])/(phis[m] - phis[n])
            else:
                l *= (phipt - phis[n])/(phis[m] - phis[n])

        xpos += xs[m]*l

    return xpos

def evolve(nx, kappa, tau, tmax, dovis=1, return_initial=0):
    """
    the main evolution loop.  Evolve

     phi_t = kappa phi_{xx} + (1/tau) R(phi)

    from t = 0 to tmax
    """

    # create the grid
    gr = Grid(nx, ng=1, xmin = 0.0, xmax=100.0,
              vars=["phi", "phi1", "phi2"])

    # pointers to the data at various stages
    phi  = gr.data["phi"]
    phi1 = gr.data["phi1"]
    phi2 = gr.data["phi2"]

    # initialize
    gr.initialize()

    phi_init = phi.copy()

    # runtime plotting
    if dovis == 1: plt.ion()

    t = 0.0
    while t < tmax:

        dt = est_dt(gr, kappa, tau)

        if t + dt > tmax:
            dt = tmax - t

        # react for dt/2
        phi1[:] = react(gr, phi, tau, dt/2)
        gr.fillBC("phi1")

        # diffuse for dt
        phi2[:] = diffuse(gr, phi1, kappa, dt)
        gr.fillBC("phi2")

        # react for dt/2 -- this is the updated solution
        phi[:] = react(gr, phi2, tau, dt/2)
        gr.fillBC("phi")

        t += dt

        if dovis == 1:
            plt.clf()
            plt.plot(gr.x, phi)
            plt.grid()
            plt.xlim(gr.xmin,gr.xmax)
            plt.ylim(0.0,1.0)
            plt.title("Reaction-Diffusion, <span class="math-container">$t = {:3.2f}$</span>".format(t))
            plt.draw()
            plt.pause(0.1)

    if return_initial == 1:
        return phi, gr.x, phi_init
    else:
        return phi, gr.x

kappa = 1.0
tau = 0.25
nx = 256
tmax1 = 1.0
phi1, x1 = evolve(nx, kappa, tau, tmax1)
</code>```

As far as I can tell, the initial condition being of class $\cal{C}^{\infty}$, and $1$ being stable, the solution should remain $1$ where the initial condition is $1$. But this is not what I observe.

Is this a numerical artifact?

ContentLicense: CC BY-SA 4.0


问:
Tridiagonal Solver in Python
I have a code I'm working on that involves solving a 1D Schrodinger equation using a Crank-Nicolson time step.  The code is written in NumPy/SciPy, and I was doing a bit of profiling and discovered that the bulk of the simulation time is taken up with calls to the linear solver (spsolve in this case).  I'm using a simple discretization such that my matrices are symmetric and tridiagonal, and I was wondering if anyone had tried doing things like implementing the Thomas algorithm directly and if this was able to beat spsolve.</p>

Title: Tridiagonal Solver in Python


问:
Convergence rate of Picard iterations
Given a first order ODE $y'(x)=f(x,y)$ with the initial condition $y(x_0)=x_0$ such that it satisfies Picard thoerem of existence and uniquness, one can compute the solution by Picard iterations :

$$ y_0 (x) \equiv y_0 \, ,$$ and
$$ y_{n+1}(x) = y_0 + \int\limits_{x_0}^x f(t,y_n (t)) \, dt \, .$$

Computationally speaking,
how fast is this method? Do we have a theorem about its convergence rate? Do we have particulary "bad" examples for its convergence?

Generally speaking, I know this method is not the usual go-to for numerical solutions of ODE's, but why?

Title: Convergence rate of Picard iterations


问:
Quadrature in finite element methods | How should I compute integrals involving the solution of the last time step?
Let $\Delta\subseteq\mathbb R^2$ denote the triangle spanned by $(0,0)$, $(1,0)$ and $(0,1)$ and $$\mathbb P_r(\Delta):=\left\{p:\Delta\to\mathbb R\mid p(x)=\sum_{|\alpha|\le r}\lambda_\alpha x^\alpha\text{ for all }x\in\Delta\text{ for some }(\lambda_\alpha)_{|\alpha|\le r}\subseteq\mathbb R\right\}$$ for $r\in\mathbb N_0$.

I'm numerically solving a PDE on a rectangle $\Lambda=(0,a)\times(0,b)$ triangulated using quadratic Lagrange finite elements as depicted in the following figure:

In the assembly of the linear system, I need to compute integrals of the form $$\int_\Delta(u^0\circ F_{\tilde\Delta})\psi\;,\tag1$$ where $u^0$ is the solution of the previous time step, $F_{\tilde\Delta}$ is the transformation from a finite element $\tilde\Delta$ to $\Delta$ and $\psi\in\mathbb P_2(\Delta)$.

How should I compute these integrals?

By definition of the finte element space, $\left.u^0\right|_{\tilde\Delta}\in\mathbb P_2(\tilde\Delta)$. So, the integrand in $(1)$ belongs to $\mathbb P_4(\Delta)$.

That's why I guess I should use a quadrature scheme which exactly integrations $\mathbb P_4(\Delta)$-functions. However, since I only now the values of $u^0\circ F_{\tilde\Delta}$ at $a^0:=(0,0)$, $a^1:=(1,0)$, $a^2:=(0,1)$, $a^3:=(1/2,1/2)$, $a^4:=(0,1/2)$ and $a^5:=(1/2,0)$, the Gaussian quadrature ansatz $$\int_\Delta f\approx\frac12\sum_{i=0}^5w_if(a^i)\tag2$$ is somehow limited. In fact, with these $a^i$ it's only possible to exactly integrate $\mathbb P_2(\Delta)$-functions (and we have $w_0=w_1=w_2=0$, $w_3=w_4=w_5=1/3$).

So, is there no possibility to do better?

Please suppose that we have a fixed mesh (i.e. no adaptive refinement) we need to deal with.


问:
Search direction for CG method
I am studying optimization methods and I was able to understand and derive the search direction

$$ p_k = r_{k-1} + \beta p_{k-1} $$

for Conjugate Gradient Method, with

$$ \beta = -\frac{p_{k-1}^TAr_{k-1}}{p_{k-1}^TAp_{k-1}}. $$

In this expression $p_k$ is the search direction and $r_k$ the residual. Nevertheless, in a paper [1] I found the following expression ($d$ is the direction, $g$ is the gradient and $y_{k-1} = (g_k-g_{k-1}))$:

$$ d_1 = -Hg_1,\\ d_k = -Hg_k + \frac{y_{k-1}^THg_k}{y_{k-1}^Td_{k-1}}d_{k-1} .$$

I am not able to understand if this is the same as $p_k$. In the paper the author writes that the basic CG is with $H$ being the identity matrix.

Is $\frac{y_{k-1}^THg_k}{y_{k-1}^Td_{k-1}} = \beta$?

[1] A Relationship between the BFGS and Conjugate Gradient Algorithms and Its Implications for New Algorithms, L. Nazareth.

Title: Search direction for CG method


问:
Extrapolation after successive finite element refinement
I wish to compute a certain function $\lambda$ (in my case an eigenvalue of the Laplacian) to a certain accuracy, which I wish to guarantee, if possible. I started with a finite element method. I know that in 2D there are way better methods, but my goal is to work on general surfaces, so that's why I work with finite elements for now. Let $h$ be the size of the elements.

I start from some $h_0$ and keep refining using the midpoint refinement procedure (replace each triangle with four triangle constructed using midpoints). This procedure gives me an approximation of $\lambda$ depending on $h$: $\lambda_h,\lambda_{h/2},...,\lambda_{h/2^n}...$ Due to the size of the computations involved I can only do $9$ or $10$ refinements which seem to get me $4$ digits of precision.

After reading a few chapters from the extremely useful book SIAM 100 digit challenge I decided to try applying a basic extrapolation technique to the sequence of approximations $\lambda_h,\lambda_{2h},...,\lambda_{h/2^n}$. To my surprise the extrapolation seems to double the number of correct digits.

I searched a bit and did not found references justifying such a behavior. So my question is:

Do you know any papers dealing with this kind of approximations:

compute some finite element approximations for successive mesh refinements

apply an extrapolation procedure

guarantee the extra precision of the extrapolated result ?


问:
How to avoid the Broyden's jacobian approximation becoming poorer with the number of iterations?
I have to solve many times a nonlinear system of the form
$$f(x) = b^{(n)}$$
inside a loop.
The function $f$ is expensive to compute and I do not have its jacobian, so I have tried the good Broyden's method. As the initial guesses of the solution, $x_0^{(n+1)}$, and the jacobian, $J_0^{(n+1)}$, at the (n+1)th iteration of the loop, I am using the solution of the previous one, $x^{(n)}$, and the approximation to the jacobian at $x^{(n)}$ provided by Broyden's method, $B^{(n)}$.

I have observed that, along the simulation, the Broyden's method needs more and more iterations to find the solutions. I guess this happens because the approximation $B^{(n)}$ to the jacobian that the method uses becomes poorer and poorer (sometimes it even becomes singular). Is there any way to avoid this?

PD: In this link http://www.math.hkbu.edu.hk/~zeng/Teaching/math3620/Broyden.pdf it is said that Powell ("A Hybrid Method for Nonlinear Equations") presented a modification to Broyden's jacobian update so that $B$ converges to the jacobian as $x$ converges to the solution. I cannot find any paper in which this is explained. Does anyone know where to find this modification explained?


问:
boundary condition in 2-D planar finite difference problem
I'm working on a 2-D planar finite difference code. My differencing scheme at the boundary nodes involves introducing ghost nodes in the computation.
My code also involves a multi-dimensional derivative:
$$
\frac{d^2u}{dxdy} \approx \frac{u_{i+1,j+1} - u_{i+1,j-1} - u_{i-1, j+1} + u_{i-1,j-1}}{4\Delta x\Delta y}
$$

This scheme becomes an issue at a corner node. Let me illustrate:
$$
g1\  \ g2 \ \ g3 \\
g4\ \  d1 \ \  d2 \\
g5\ \  d3 \ \ d4 \\
$$

the g's represent the ghost codes and the d's represent the real domain nodes that I'm solving for. x-direction is horizontal; y-direction is vertical. Dirichlet boundary conditions are specified for both the top and left boundaries. Let's say $u=u_1$ at the top boundary and $u=u_2$ at the left boundary. What should be the boundary condition at ghost node g1, which lies along both the x and y boundaries?


问:
scalable parallel mesh/amr on unstructured grid
I am trying to code a scalable parallel AMR for unstructured grid. There seems to be three approaches for this

a) Store some global grid info on each processor and partition with parmetis (The easiest approach I am pursuing). I have never used ParMetis before (just metis) so I do not know how one would be able to read a given mesh file in parallel for partitioning, and then later on store the file in a format readable in parallel.

b) Store global grid on rank=0, use Metis and then scatter. I don't like this approach because it is more involved than (a), and you still need one node that is able to store the whole grid anyway.

c) Use an approach like p4est, which only works for structured quads/hexes. I am not sure how much global information is stored on each processor but it seems pretty scalable. Can the techniques used for the mesh encoding in p4est be used for unstructured grids as well. It seems to me even if one stores only one integer, say to indicate to which processor an element belongs, sooner or later scalability will be limited due to not being able to store that information one a node. What I think is done in p4est is to partition a coarse grid (able to be stored on each node), and do refinement after partitioning, but this approach can not be used for unstructured AMR, right ?

Thanks for any input.


问:
Kahan Summation for Three-Term Recurrences
Kahan summation applies to summation problems, but not to three-term recurrence relations. However, a three-term recurrence shares many of the features of a summation-albeit with a rescaling step at each iteration. Hence there seems to be hope for a Kahan summation for three-term recurrences. Has this idea been worked out in the literature or is the idea flawed in some way?</p>


问:
Estimate $L_2$ norm of a elliptic problem with unknown exact solution on finite element method
I have the elliptic problem $$-\Delta u = 1,\,\,\Omega\subset\mathbb{R}^2$$
with $u=0$ on $\partial\Omega,$ with $\Omega=[-1,1]^2\backslash([0,1]\times[-1,0])$ and I want to estimate the $L_2$ error of linear finite element method on triangular mesh, if $u_h$ is the finite element solution and $u$ is the exact solution (unknown),
$$\|u-u_h\|^2_{L_{2}(\Omega)}=\int_{\Omega}|u-u_h|^2\,dx.$$

Ι though to solve this problem with quadratic finite elements and set $u=\tilde{u}_h,$ where $\tilde{u}_h$ is the solution of quadratic finite elements.


问:
How to perform Taylor series expansion consisting of cell averaged derivatives in a computational element?
I have encountered the term cell "Cell averaged derivatives" and "Taylor series expansion using cell averaged derivatives about centroid" in the context of polynomial representation in any cell (computational element with some higher order interpolation polynomial nodes).

Here is the link to the below paper:
http://dept.ku.edu/~cfdku/papers/AIAA-2009-605.pdf

I have also read the reference of this reference paper
http://people.math.gatech.edu/~yingjie/publications/centr_DG_ovlpcell_new_iii.pdf
Section 4.1

I tried to find alternate references and easier sources of explanation, but everywhere it has been mentioned that the coefficients of the above Taylor series expansion can be found easily.

However, I am still unable to grasp the procedure properly. It would be nice if someone could help me understand which seems to be trivial for everyone.


问:
Pseudoinverse of a large sparse matrix in r
This question was moved from Cross-Validated: https://stats.stackexchange.com/questions/274042/pseudoinverse-of-large-sparse-matrix-in-r

I am trying to calculate the pseudoinverse of a large sparse matrix in R using the singular value decomposition. The matrix is roughly 240,000 x 240,000, and I have it stored as type dgCMatrix. I have tried using pinv, ginv, and other standard pseudoinverse functions, but they error out due to memory constraints. I then tried to opt for the sparse matrix svd provided by the package irlba, which I was then going to use to compute the pseudoinverse using the standard formula after converting all outputs to sparse matrices. My code is here:

```<code>lim = 40
digits = 4
SVD =irlba(L,lim)
tU = round(SVD$u,digits)
nonZeroU = which(abs(U)>0,arr.ind = T)
sparseU = sparseMatrix(i=nonZeroU[,2],j=nonZeroU[,1],x = U[nonZeroU])
V = round(SVD$v,digits)
nonZeroV = which(abs(V)>0,arr.ind = T)
sparseV = sparseMatrix(i=nonZeroV[,1],j=nonZeroV[,2],x = U[nonZeroV])
D = as(Diagonal(x=1/SVD$d),"sparseMatrx")
pL =D%*%sparseU
pL = sparseV%*%pL
</code>```

I am able to get to the last line without an issue, but then I get an error due to memory constraints that says

```<code>Error in sparseV %*% pL :
  Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 105
</code>```

Of course I could piece together the pseudoinverse entry by entry using a for loop and vector multiplications, but I would like to be able to calculate it using a simple function that takes advantage of the sparsity of the resultant pseudoinverse matrix. Due to the nature of the original matrix (it is a graph laplacian), I know that the pseudoinverse should also be a sparse matrix. Any help would be greatly appreciated!


问:
Outflow boundary condition - second derivative of velocity
Consider a fluid flow simulation in a pipe. At the outflow, instead of explicitly imposing a boundary condition, I linearly extrapolate information from the interior (for velocity components). This will translate to $$\frac{\partial ^2 u}{\partial x^2} =\frac{\partial ^2 v}{\partial x^2}=\frac{\partial ^2 w}{\partial x^2}= 0$$ at the boundary($x$ being axial direction). This kind of numerical condition has been used successfully throughout the literature (see page 262 of this paper).

What I am interested is - what does this condition physically mean?

Any discussion would be greatly appreciated.


问:
Solve scalar quadratic equation
I am implementing a Trust-Region method using Dogleg for the search direction defined as $\tau p^U$ for $0 \leq \tau \leq 1$ and $p^U + (\tau -1)(p^B-p^U)$ if $1 \leq \tau \leq 2$. To compute $\tau$ we have to solve the following equation:

$\| p^U + (\tau - 1)(p^B-p^U) \|^2 = \Delta ^2.$

I tried to expand it in order to arrive at an equation of the form $ax^2+bx+c=0$. I obtain:

$$ (\tau -1)^2\|p^B-p^U\|^2 + 2(\tau -1)(p^U,p^B-p^U) + \|p^U\|^2 - \Delta ^2 = 0 .$$

I do not not if it is correct or not and how to proceed. Any suggestion?

Title: Solve scalar quadratic equation


问:
Matlab CODE of ADER scheme
I read  ADER WENO Finite Volume scheme for hyperbolic conservation laws with source term. I want to implement this method in Matlab coding. I got some problems in coding during to computation of
numerical flux. Please help me the Matlab coding
$$u_t+f(u)_x=S(u,x,t),\; u(0,x)=u_0(x),\; x\in \mathbb{R}, \;t>0,\;u(t,-\infty)=u(t,\infty)=0.$$</p>

Title: Matlab CODE of ADER scheme


问:
MATLAB: solving multiple ODE systems in parallel
I have a system of parameterized ODEs that I would like to solve using MATLAB and its ode45 solver, and was wondering if it is possible to perform such a task in parallel. For example, I have a simple model of a rocket falling in gravity, where I would like to try different parameter values of tburnstart and tburnend, which tell the program when to turn a rocket's engines on and off:

```<code>function xdot = equations_of_motion_1D_landing(t, x, G, g, M, Isp, T, tburnstart, tburnend)

    if t < tburnstart
        xdot = [
                x(3); %xdot1: x-velocity
                x(4); %xdot2: y-velocity
                -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2)); %xdot3: x-acceleration
                -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2)); %xdot4: y-acceleration
                0 %xdot5: engine mass flow rate
                ];
    else if  t >= tburnstart && t <= tburnend
        xdot = [
                x(3); %xdot1: x-velocity
                x(4); %xdot2: y-velocity
                -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2)); %xdot3: x-acceleration
                -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2))+(T/x(5)); %xdot4: y-acceleration
                -T/(g*Isp) %xdot5: engine mass flow rate
                ];
    else
        xdot = [
                x(3); %xdot1: x-velocity
                x(4); %xdot2: y-velocity
                -(G*M*x(1))/((x(1)^2+x(2)^2)^(3/2)); %xdot3: x-acceleration
                -(G*M*x(2))/((x(1)^2+x(2)^2)^(3/2)); %xdot4: y-acceleration
                0 %xdot5: engine mass flow rate
               ];
    end
end
</code>```

Currently I test different values for tburnstart and tburnend using the following code in my main file:

```<code>%Setting up burn start and end times
tmin = 0;
tmax = 200;
tburnstart = (tmax-tmin).*rand(pop_size,1)+tmin; %burn start times
tburnend = (tmax-tburnstart).*rand(pop_size,1) + tburnstart + 10; %burn end times (must be larger than start times) (+10 makes sure tspan has more than 1 element)
pop_init = [tburnstart, tburnend]; %creating intial population

for i = 1:length(pop_init)
    tburnstart = pop_init(i,1);
    tburnend = pop_init(i,2);
    %tspan = [0:tburnend];
    tspan = 0:0.1:tburnend;
    [t,x] = ode45(@(t,x) equations_of_motion_1D_landing(t, x, G, g, M, Isp, T, tburnstart, tburnend), tspan, xinit);
end
</code>```

However, using a for loop to evaluate ode45 is extremely slow, and I was wondering if there is a way to evaluate all parameters simultaneously to save on computing time.


问:
Stability of nonlinear partial differential equation
I want to find an expression for the stability of the nonlinear Poisson equation. I know about von Neumann stability analysis which applies to linear equations as far as I know. Any suggestion how to analyse the stability?

The  equation I am trying to solve is:

$$-\frac{\partial ^{2}\phi}{\partial \tau^{2}}+\alpha\frac{1}{r'^{2}}\frac{\partial}{\partial r'}\left(r'^{2}\frac{\partial\phi}{\partial r'}\right) = \left(\frac{\gamma}{\phi^{n+1}}+\beta\right)\,. $$

And my discretisation looks like:

$$ a\phi_{i-2,j}-2a\phi_{i-1,j}+b\phi_{i,j+1}+ c\phi_{i,j-1}+d= e \phi_{i,j}+\frac{h}{\phi_{i,j}^{n+1}} $$

$a, b, c, d, e, h$ are constants, and $i$ and $j$ indices stand for time and radius.


问:
Discontinuos Galerkin Method - inhomogeneous advection problem
I'm currently trying to get into this topic.  I've learned that the basic scheme for the advection problem ($D_{x}u+a*D_{t}u=0$) can be solved in a scheme like $$ M^{k}\frac{d}{dt}u^{k}_{h}-(S^{k})^{T}au^{k}_{h}=-(au_{h})^{*}\psi(x^{k}_{r})+(au_{h})^{*}\psi(x^{k}_{l})$$Here, we have $$ M^{k}_{ij}=(\psi_{i},\psi_{j}),  S^{k}_{ij}=(\psi_{i},\frac{d\psi_{j}}{dx})$$ $\psi_{i}$ are the basis functions.
 I think I understand this scheme but I'm not quite how to deal with the case of an inhomogeneous PDE.
 How has the scheme to be modified if I've got $f(x,t)$ instead of 0 above?


问:
Composite Laminate Mass Matrix Help
So, as the title of the post implies, I'm making an FEM code for a composite laminate plate using the Mindlin plate theory. And, as I am trying to get the mode shapes for the plate, I needed to make a mass matrix. Following some other examples and some past code someone gave me, I managed to make one.

But the problem is, when I change the number of elements, like for example, from 4x4 to 5x5 elements, some of the mode shapes, the 4th one specifically, flips upside down. And that is the case when I change to number of elements. Sometimes the 1st and 2nd mode shapes are flipped. Sometimes all of them are. But I still can't find what is wrong.

What could be causing this? Is there some checklist of simple things I should check to correct my mass matrix? If there is an issue you think that might be a problem, please be as specific as possible in coding terms. Thanks in advance!

Title: Composite Laminate Mass Matrix Help


问:
Nonlinear 2D thermoconductivity equation(numerical solution)
I have to write a solver for 2D equation:
$$\partial_t u = u^2(\partial_x ^2 u + \partial_y ^2 u)$$
I try to use explicit method:
$$\partial_t u = \frac{u_{i,j}^{k+1} - u_{i,j}^k}{\tau}$$ and $$\partial_x ^2 u = \frac{u_{i-1,j}^k - 2u_{i,j}^k + u_{i+1,j}^k}{h^2}$$ and tha same for $\partial_y ^2 u$ with points $j-1, j, j+1$. Also I use: $$ u_{i,j}^k = \frac{u_{i+1,j}^k + u_{i,j}^k}{2}$$ for the coefficient before $\partial_x ^2$ and the same for $y$ with points $j+1, j$.
I wrote the following code in Mathematica:

```<code>nn = 6;
mm = 200;
h = 1./nn;
tau = 1./mm;
v0[x_, y_] := 1.;

Table[{x[i] = i*h, y[j] = j*h}, {i, 0, nn}, {j, 0, nn}];
Table[t[i] = i*t, {i, 0, mm}];
Table[u[i, j, 0] = v0[x[i], y[j]], {i, 0, nn}, {j, 0, nn}];
Table[u[0, j, k] = 0, {j, 0, nn - 1}, {k, 0, mm}];
Table[u[nn, j, k] = 0, {j, 0, nn - 1}, {k, 0, mm}];
Table[u[i, 0, k] = 0, {i, 0, nn - 1}, {k, 0, mm}];
Table[u[i, nn, k] = 0, {i, 0, nn - 1}, {k, 0, mm}];

Do[
  u[i, j, k + 1] =
   u[i, j, k] +
    tau*(((u[i - 1, j, k] - 2*u[i, j, k] + u[i + 1, j, k])/
          h^2)*(( u[i, j, k] + u[i - 1, j, k])/
           2)^2 + ((u[i, j - 1, k] - 2*u[i, j, k] + u[i, j + 1, k])/
          h^2)*((u[i, j, k] + u[i, j - 1, k])/2)^2), {k, 0, mm}, {i,
   1, nn - 1}, {j, 1, nn - 1}];

g[k_] := Table[u[i, j, k], {i, 0, nn}, {j, 0, nn}];
p = Table[
  ListPlot3D[g[k], PlotRange -> {{0, nn}, {0, nn}, {0, 6}}], {k, 0,
   mm, 20}]
</code>```

I use the simpliest conditions:
$u(0,0,t) = u(10,10,t) = 0$ and $u(x,y,0) = 1$


问:
Mathematical physics applications in present-day image processing
During the past few years several important areas of image processing and image classification or generation became dominated by convolutional neural networks.

I'm interested if there are any methods coming from mathematical physics context (like methods designed for solving ill-posed problems, spectral analysis, image deblurring and deringing), that outperform neural network-based approaches in 2017 for some common computer vision or image processing problem. Or methods that don't have any neural network-based rivals. Maybe in the field of biomedical image analysis (just a guess)?

The deeper and more specific the answer is, the better.


问:
How to do numerical computation of $L^p$ norm of a $p$ dimensional trigonometric polynomial
Id like to know methods for  numerical computation of $L^2$ norm of a two dimensional trigonometric polynomial.

I have the coefficients. If I want to compute the L^1 norm, I can do so by sampling in time at a rate greater than twice of the order of the trig polynomial (Shannon sampling theorem) and then doing summation of samples which can approximate the L^1 norm integral quite well. But for L^2 norm I need to square in time domain and cannot use this procedure as squaring is non linear?

PS : I dont want to use Parseval's theorem, as it wont help if I want to compute L^3 norm of a 3 dimensional trigonometric polynomial. I want a general procedure for $L^p$ norm of a $p$ dimensional trig polynomial.


问:
How do we implement Parameter free generalised Moment limiter in 1D Case in Discontinuous Galerkin methods?
I am referring to this paper:-

"A Parameter-Free Generalized Moment Limiter for High-
Order Methods on Unstructured Grids " by Michael Yang and Z.J. Wang.

http://dept.ku.edu/~cfdku/papers/AIAA-2009-605.pdf

They have implemented a general limiter which doesn't require any user input. In the process, they have represented the solution in the cell in terms of the average of derivatives (Eqn. 28 of paper).

Can someone please explain how was this representation attained. The understanding of this is a quite high priority for me. Please help.


问:
How to implement outgoing wave boundary condition
I am solving the one dimensional wave equation:

$0=\Box\phi = -\partial_t^2\phi + \partial_r^2\phi ,$

using a Crank-Nicolson finite difference scheme, in the domain $r\in[0,R]$. First, I define $\xi\equiv\partial_r\phi$ and $\Pi\equiv\partial_t\phi$, so that the equations of motion can be written as follows:

$
0 = \partial_t\xi - \partial_r\Pi, \\
0 = \partial_t\Pi - \partial_r\xi,
$

Specifically, the stencils I'm using for the interior grid points are:

$
(\partial_t u)^{n+1/2}_j = \frac{u^{n+1}_j - u^n_j}{\Delta t} + \mathcal{O}\left((\Delta t)^2\right) \\
(\partial_r u)^{n+1/2}_j = \frac{1}{4\Delta r}\left(u^{n+1}_{j+1} - u^{n+1}_{j-1} + u^n_{j+1} - u^n_{j-1}\right) + \mathcal{O}\left((\Delta t)^2, (\Delta r)^2\right)
$

The Crank-Nicolson scheme is implicit but I am solving the implicit equations using Gauss-Seidel iteration to solve for $\xi^{n+1}_j$ and $\Pi^{n+1}_j$. I'm having no trouble evolving my initial data (a Gaussian profile for the $\phi$ field, moving the the left) in the interior grid points. I'm only running into trouble at the boundaries of my simulation.

I would like to implement outgoing wave boundary conditions at the boundaries of my simulation, i.e. at $r=0$

$
0 = \partial_t\xi - \partial_r\xi, \\
0 = \partial_t\Pi - \partial_r\Pi,
$

and at $r=R$

$
0 = \partial_t\xi + \partial_r\xi \\
0 = \partial_t\Pi + \partial_r\Pi
$

For some reason I am having a hard time numerically implementing these boundary conditions. I have tried several approaches (including first and second order upwind schemes) to implement the boundary conditions but when I run my code either I get what appears to be a numerical instability or the the wave solution 'bounces back' into my domain; i.e. it doesn't exit my domain contrary to the boundary condition!

I've been having difficulty finding a good reference/discussion on outgoing wave boundary conditions. Could someone provide some insight on how to implement outgoing wave boundary conditions using an implicit method? Thank you in advance for the help.


问:
Using Newton-Raphson method to solve the hydrostatic equation
I'm trying to use newton-raphson method for nonlinear systems of equations as described in 'Numerical recipies' book in chapter 9.6 to solve the hydrostatic equation for a polytropic star.

For each iteration, I change the radius vector as described in the book like x→x+Δx(Δ can be calculated for each iteration) but I can't ignore the fact that it can do shell crossing - meaning the i shell moves to a higher shell then i+1i+1, this cause a negative volume. I know there is a method to avoid shell crossing by multiplying delta with a factor but I don't know how to calculate it. There is no description in the book about shell crossing.

let me be more specific:

I want to solve $F_i(x_1,....x_N)=0$ set of equations,( when i goes from 1 to N).
the method is to expand $F$ to a taylor series:
$\vec{F}(\vec{x}+\vec{\delta x})=\vec{F}+J\vec{\delta x} +O(x^2)$
, where $J$ is the partial deriviatives of $F$. by setting $\vec{F}(\vec{x}+\vec{\delta x})=0 $  we can get N sets of equations : $J \cdot \vec{\delta x}=-F$ that can be solved by LU decomposition.

the method starts from initial point $x_{old}$. for each step we calculate $\delta x$ and add it to $x_{old}$ , meaning : $x_{new}=x_{old}+\delta x$

if the newton step $\delta x$ minimizes $f=0.5\cdot F\cdot F$ we can say this step is acceptible. By using line searches and backtracking we can guaranty to minimize $f$. Be aware that every solution to $F_i(x_1,....x_N)=0$  is a minimum of $f$ but not the oposite, so we must check if $f$ goes to zero.

My set of equations are:
$F_i=\frac{dP_i/dm_i}{Gm_i/4\pi r_i^4}+1=0$

where $\delta P_i=P_{i-1}-P_i$ , $\delta m_i =0.5\cdot (m_{i+1}-m_{i-1})$ , and from the polytropic equetion: $P_i=K(\frac{m_i-m_{i-1}}{4\pi / 3 (r_i^3-r_{i-1}^3)})^\gamma$

Where do i stuck? in the first iteration. and it's suppose to be bad. the problem is shell crossing. $\delta x$ changes $r$ so much that it changes the order of the cells. meaning, $r_i$ the radius of the star (from zero) is suddenly bigger then $r_{i+1}$.  I know there is a method to correct this by multiplying $\delta$ by some factor. How to calculate it? Is there any other methods to correct it?

Do anyone know how to help?


问:
solving numerically a 2D integral by using simps and quad combined
I have the following function

$$f(x,y) =\left(\frac{1}{\exp(x-E_f)+1}-\frac{1}{\exp(x-y-E_f)+1} \right)\frac{1}{\sqrt{4-xy}}$$

I want to integrate it using Python 3. The domain of the $y$ variable is from sampled data hence I cannot used something like quad. The $x$ variable comes from analytical formula so I can generate whatever value I want for it. My question is how can use simps and quad combined to integrate this formula.


问:
How is the Gastner-Newman equation implemented to create value-by-area cartograms?
There is a paper called "Density-equalizing map projections: Diffusion-based algorithm and applications" by Michael T. Gastner and M. E. J. Newman, which explains their algorithm (which is based in diffusion equations) for generating value-by-area cartograms.

While it explains the theoretical side of the mathematics involved with their algorithm, it doesn't explain how they actually implemented it. I tried to piece it together by looking at the source code from cart, but I don't have the programming knowledge (it's written in c, which I don't know) required to understand it.

If anyone has at least a decent understanding of it and can explain the steps needed to create a cartogram using their algorithm, that would be greatly appreciated. Otherwise, if you have other helpful resources on the topic, those would be good too.


问:
Constrained optimization: Stationary point vs. Nash point
1s question: definition of stationary point for constrained optimization

As far as I know, a stationary point of a constrained optimization problem is a stationary point of the Lagrangian (that has to be differentiable).

Now consider the problem
\begin{equation}
\min_{x\in X} F(x)
\end{equation}
where $F:\mathbb{R}^p\to \mathbb{R}$ is differentiable but possibly non-convex, $X\subset \mathbb{R}^p$ is a closed convex set.

It is easy to obtain the following first-order optimality condition, where $x^*$ is a local minimizer:
$$\nabla F(x^*)^\top(x-x^*) \ge 0\quad\forall x\in X. \qquad (*)$$
In several references that I read (such as this one), $x^*$ is called a stationary point if and only if it satisfies the above condition.

My question is: Is this definition of stationary point widely accepted? Is there any reliable reference on that? And, is there a similar definition if $f$ is non-differentiable?

2nd question: Stationary point vs. Nash (equilibrium) point

Consider now that $F(x)$ takes the form
$$F(x_1,...,x_n)= f(x_1,...,x_n) + \sum_{i=1}^n r_i(x_i)$$
where $f:\mathbb{R}^p\to \mathbb{R}$ is a differentiable multi-convex function, $r_i:\mathbb{R}^{p_i}\to \mathbb{R}$ are extended-value convex functions. Also, $X$ is a closed multi-convex set. (Roughly speaking, the problem of minimizing over one block while the others are fixed is a convex problem.)

This problem is considered in this paper.

Nash equilibrium (eq. (2.3) in the paper, reformulated):
$$x_i^* = \arg\min_{x_i} F(x_1^*,.\ldots.,x_{i-1}^*,x_i,x_{i+1}^*,\ldots,x_n^*) \qquad (2.3),$$
which is equivalent to (eq. (2.4) in the paper):
$$\left(\nabla_{x_i}f(x^*) + p_i^*\right)^\top (x_i - x_i^*) \ge 0 \quad\forall x_i \qquad (2.4),$$
where $p_i^*\in\partial r_i(x_i^*)$.

In Remark 2.2 (right after equation (2.4)), the authors stated that:

In general, the condition (2.4) is weaker than the first-order
  optimality condition. For our problem, a critical point must be a
  Nash point, but a Nash point is not necessarily a critical point.

I am not sure how the authors defined 'critical point' (= 'stationary point') in this case, because if it's the same as in my first question, then I think for this problem, a Nash point must be a stationary point, given that $r_i$ are differentiable.

Indeed, assume that $x^*$ is a Nash point, we have
$$\nabla F(x^*)^\top (x - x^*) = \sum_{i=1}^n \left(\nabla_{x_i}f(x^*) + \nabla r_i(x_i^*)\right)^\top (x_i - x_i^*).$$
Since each term in the sum is non-negative according to $(2.4)$, we must have $\nabla F(x^*)^\top (x - x^*) \ge 0$, i.e., $x^*$ is a stationary point according to the definition $(*)$.

What do you think?

Thanks a lot for your discussion.


问:
Can similarity criterion help reducing calculations in Finite Element Analysis (FEA) problems?
Can similarity criterion help reducing calculations in Finite Element Analysis (FEA) problems?

I'm using NX (a thermal and CFD simulation software) to simulate the temperature and fluid field of a LED application.

It takes a long time each time i change something.
I'm wondering if I can reduce the calculations by scaling down the object I'm calculating.

For example, if cooling solution A is best for an application of size 10. Is the cooling solution A also the best if I scale down the application to size 5?

If this is true, then i can base my calculations on size 5 and reduce drastically the time of calculation.

I believe this is not so simple. For example heat pipes behave differently with different lengths.

But maybe with some corrective parameters, this can be done.
Any advices ?

(I realize that this might not be the place to ask this question, please also advice where i can post this question)

OwnerDisplayName: osager


问:
How can I evaluate the accuracy of my n-body simulation?
I am making an n-body simulation in python. There are many different methods to numerically solve the system of differential equations governing the gravitational interactions between the $n$ particles; the challenge of n-body simulations is to find an integration scheme which is both accurate and fast.

For the moment, I am working with n small (say $n \le 5$). I began with the most naive method of Euler integration, which is first order. Even with two bodies, this method can produce pretty bad errors if the bodies come very close (since then the force between them becomes very large, so their accelerations become very large). Of course, one can always reduce the time-step to get more accuracy, but at the expense of computation time.

Next I would like to implement a higher order method which has better stability properties etc like leapfrog integration or Runge-Kutta.

I expect these methods to give much better results, but the problem is I am unsure how to measure quantitatively the accuracy of my simulation. What metric should I use?

One of the issues is that I can't tell if the simulation is giving realistic results, because I haven't been using the real value of the gravitational constant $G$ and none of the distances/times/masses in my simulation have units.

The reason why I haven't been paying attention to units is because, say if I wanted to do a scale sun-earth-moon simulation; since the sun-earth distance is about 370 times the earth-moon distance, it's not possible for everything to fit inside the display window and also be able to distinguish the earth and the moon.

The other problem is using the real value of G and real masses etc leads to computations involving big numbers. But, keeping everything else fixed, choosing lower values of $G$ reduces the errors (because the forces are then smaller). So I have no real way of telling if the simulation is accurate.

So how can I rigorously test the accuracy of the simulation, and what can I do concerning units?


问:
Discrete operator textbooks
This will be a vague question.

When I was writing a finite element matrix assembly routine, a colleague noticed that I had a bug in my code because the sparsity pattern of the one of the blocks didn't look correct to him.

He explained to me that the block operator is a Laplacian operator, which is the discrete analogue of the continuous second derivative operator, hence there should be no coupling between the off diagonal block terms of the matrix and therefore that I had a bug in my code, since the sparsity pattern did not match. I was blown away by the alternative way of viewing the problem, and I want to learn more.

He didn't know of any reference materials that teach things like that; he just claimed to have learned it by experience. Does anyone know of a textbook or  resource on discrete operators which can help train my brain to think about numerics in this sort of sense?

I come from a mainly computer science background, and have more coding intuition than this sort of mathematical intuition.

Title: Discrete operator textbooks


问:
Optimal synthesis of crank–rocker linkages
I need help programming a script, in MATLAB, that for assigned values compute the optimal transmission angle of a crank-rocker mechanism.

The parameters are:

$L$: crank's angle of oscillation, for example from -20° to 150°

$m$: rocker's angle of oscillation, it's a fixed value.

$t$:  $\tan(L/2)$

$u$: $\tan((L-m)/2)$

$v$: $\tan(m/2)$

There is a condition on $L$, $L<(90+(m/2))$

The equation is a cubic:

$$Q^3 + 2Q^2 - t^2 Q- \frac{t^2(1+t^2)}{u^2}$$

where $Q=t^2/x^2$.

The roots of $Q$ must be real and positive.
$x$ is the variable and it should satisfy the following inequalities:
$x^2 \geq 1$ and $x^2 \leq u^2 t^2$.

Now the lengths of the 4-bar are
\begin{align}
a &= \sqrt{\frac{v^2}{1 + v^2}}\cr
b &= \sqrt{\frac{x^2 v^2}{1 + v^2}}\cr
c &= \sqrt{\frac{x^2 + t^2}{1 + t^2}}\cr
d &= \sqrt{\frac{u^2 + x^2}{1 + u^2}}
\end{align}

I solve the equation for $Q$ then calculate $x$. Now the transmission angle (max and min) depends on $a(x)$, $b(x)$, $c(x)$ and $d(x)$. I have a different transmission angle (max and min) for all the different values of crank's angle of oscillation. For example as, Istated previously from $-20$ to $(90+(m/2))$ (the lower bound is discretionary)
$$\begin{aligned}
 \mathit{angtrasm}_{\min}&= \arccos\left({\frac{b^2 + c^2-(d-a)^2}{2bc}}\right)
\\
\mathit{angtrasm}_{\max} &= \arccos\left({\frac{b^2 + c^2-(d + a)^2}{2bc}}\right)
\end{aligned}
$$

My code is the following

```<code>% Compute transmission angle and lengths of bars
clear all

m = 32;   %rocker’s oscillation
i = 1;
f = 1;

for l = 1: (90 + m/2 );  % crank’s oscillation -pi (there is a problem here i can’t start from -20 for example)
  L(i) = l + 180 ;       %crank’s oscillation
  t = tand(L(i)/2);
  u = tand((L(i)-m)/2);
  v = tand(m/2);
  p = [1 2 -(t^2) -((t^2 + (1 + t^2))/(u^2))];
  if L(i) = = m
    fprintf('Value of crank’s oscillation is not acceptable \n' , L(l));
  else
    Q = roots(p);
    for j = 1:3
      if isreal(Q(j)) && Q(j)>0
        x = abs(t/sqrt(Q(j)));
        if x^2 > = 1 && x^2 < = u^2*t^2
          a = sqrt(v^2/(1 + v^2));
          b = sqrt(x^2*v^2/(1 + v^2));
          c = sqrt((x^2 + t^2)/(1 + t^2));
          d = sqrt((u^2 + x^2)/(1 + u^2));
          angtrasm_min(f) = acosd((b^2 + c^2-(d-a)^2)/(2*b*c));
          angtrasm_max(f) = acosd((b^2 + c^2-(d + a)^2)/(2*b*c));
          scarto_min(f) = abs(90-angtrasm_min(f));
          scarto_max(f) = abs(90-angtrasm_max(f));

          theta(f) = L(i);
          r1(f) = a;
          r2(f) = b;
          r3(f) = c;
          r4(f) = d;
          fprintf('The max transmision angle is %d when theta is %d\n',angtrasm_max(f),L(i));
          fprintf('The max transmision angle is %d when theta is %d\n',angtrasm_min(f),L(i));

          fprintf('.................................\n');
          f = f + 1;
        end
      end
    end
  end
  i = i + 1;
end
</code>```


问:
Use custom defined metric(not matrix) in Javaplex
I'm about to get started trying to use Javaplex for the first time for some application in topological data analysis (TDA) and I want to use a custom defined metric. From looking over the documents, I can do this if I build the distance matrix for all of my points. My question is twofold:

Can I do it by just defining the metric between two points (it's a simple function to define, just the Euclidean cross-product metric between a Euclidean space and a circle)

and if so,

Would it make sense computationally, or should I just build the distance matrix for each of my datasets? For the context, I have 800 datasets, each with about 40-50 points on a torus in $(x,y,\theta)$-coordinates, so I need to define my torus metric.


问:
Using RK2 Method to solve the simple harmonic oscillator of a horizontal mass on a spring (1D)
Being new to numerical analysis techniques, in particular RK2, I decided the best way to jump in is by using python to solve the well known mass-spring oscillator using RK2 techniques.

My problem is that the step size seems to influence the period of my solution and I do know why and would appreciate it if someone could help me better understand what I am doing wrong. Commenting on the code is welcomed! Thank you.

The problem is as follows (this is not a textbook or homework problem!): A mass is attached to a spring with a linear spring constant. At time $t=0$ the mass is displaced from equilibrium, $x=0$, to the position $x(t=0) = 1 cm$. The mass is released and left to oscillate about equilibrium. Use the method of RK2 to plot the evolution of the masses position as a function of time. The required first order differential equations for the problem are:

$\dot{x} = -\omega x(t)$

$\frac{dx}{dt} = \dot{x}$

Where the RK2 functions I am using are as follows:

$v(t_{n+1}) = v(t_n) -\omega x(t_n)\delta t$

$x(t_{n+1}) = x(t_n) + v(t_n)\delta t $

$ k2_{x} = \delta t (v(t_n) - \omega x(t_n)\frac{\delta t}{2}) $

$ k2_{v} = \delta t (\omega x(t_n) + v(t_n) \frac{\delta t}{2}) $

The last two function are the midpoint values that I think are used for the RK2 approach.

Along with the implemented code :

```<code>import numpy as np
import matplotlib.pyplot as plt

print "RK2 Method for Oscillating Spring"

#Spring Constant
k = 1.0
#Mass
m = 0.02
#Angular frequency squared
w = k/m

period = 2*np.pi*(np.sqrt(1/w))

#Initial Conditions
x_start = 1.0
v_start = 0.0

#Set loop variables
x_position = x_start
velocity = v_start

#Step Size
dt = 0.001

#Arrays to store variable values
x_plot = []
v_plot = []
time = []

#Plot the numerically determined position over
# the time interval 1 to 10 seconds in steps of 0.1
for t in np.arange(1,10,0.1):

    #The two equations below are used in the Euler Method approach
    #Their solution's period, which can be plotted, varies depending on the step size
    # and I know it shouldn't
    #x_position = x_position + dt*velocity
    #velocity = velocity - dt*w*x_position

    #The equations below are used for RK2 method. I don't know why I cant get it to work.
velocity = velocity + dt*(-w*x_position)
    k2_x = dt*(velocity + (-w*x_position*dt/2.0))
    x_position = x_position + k2_x

    k2_v = -dt*w*(x_position + velocity*dt/2.0)
    velocity = velocity + k2_v

    x_plot.append(x_position)
    v_plot.append(velocity)
    time.append(t)

plt.subplot(3,1,1)
plt.plot(time,x_plot)
plt.subplot(3,1,2)
plt.plot(time,v_plot)
plt.subplot(3,1,3)
plt.plot(x_plot,v_plot) # illustrate divergence over number of steps
plt.show()
</code>```

My problem is that I don't understand why varying the step size, dt, that the period of the solution, a sinusoidal function, changes.


问:
Implement Robin boundary condition (finite volume)
I have  a PDE equation with Robin Boundary condition in an annulus system and I should solve it by finite volume method:

\begin{align}
\frac{\partial T_f}{\partial t} - k \left(\frac{\partial^2 T_f}{\partial r^2} + \frac{1}{r} \frac{\partial T_f}{\partial r} \right) = -Q
\end{align}

I do not have any problem in the discretization of inner points by finite volume method,
but in $r_{ID}$ and $r_{OD}$, it has robin boundary condition which I do not know how to discretize it.

\begin{align}
\left.\frac{\partial T_f}{\partial r}\right|_{r=r_{ID}} &= h(T_f-T_w ) \\
\left.\frac{\partial T_f}{\partial r}\right|_{r=r_{OD}} &= h(T_f-T_g )
\end{align}

(h&k are constant)

I would appreciate it if you could help me to discretize this boundary condition.


问:
minimalistic library for classical molecular mechanics force-field in C/C++/Python
I would like some classical molecular mechanics forcefield library or super-simple program which does not have any external dependences, so that I can plug it into my software without worry of "dependency hell". Also I would like to modify it, so I would prefer something with simple and compact source code. It should do just one thing - in summary it should follow unix philosophy

For example there is rdkit but it seems to be too ambitius (lot of functionality I don't need, lot of source code and dependencies to care about). Also many molecular editors like Avogadro and OpenBable have some molecular mechanics implemented ... but it would be much more convenient to have clean simple library than try to extract this functionality out of big software package.

I don't have much requirements on precision of force-field, but I would prefer something which is as general as possible (e.g. UFF) for all elements and exotic bonding topologies although organic chemistry is main target, not something specialized to e.g. proteins like GROMOS, AMBER, CHARM etc.


问:
Quick evaluation of floating point Absolute Error
I need to to find a quick and dirty way to estimate the absolute error introduced by a series of agebraic operations of IEEE single precision floating point numbers, a pessimistic result is ok.

The way I am doing it, is to evaluate the biggest possible value that can be obtained during the chain of operations, find the power of two immediately greater than that number, and multiply such power of two by the number of algebraic operations performed and by the value of the least significant bit in the significand assuming exponent 1 (which is 0.00000011920928955078125)

Does it look like a sensible approach, or am I missing something?


问:
Any method to efficiently compute SVD of a perturbation of matrix $\bf A$ if the SVD of $\bf A$ is already known?
Suppose we know the SVD of matrix $\bf A$, and $\bf B$ is a slight perturbation of $A$ (e.g. $\|{\bf B}-{\bf A}\|_{\text F}$ is relatively small), then is there any method that can efficiently compute the SVD of $\bf B$? That is, can the knowledge of SVD of $\bf A$ be helpful for SVD of $\bf B$?

I searched a little bit and found there are some papers on the bound of perturbation, e.g. Perturbation Theory for the Singular Value
Decomposition, but I currently have no luck in finding a method to compute SVD of perturbation taking advantage of the SVD of the original matrix.

Any help or reference will be very much appreciated!


问:
Backing out a function of parameters from system of nonlinear equations
I have a system of equations that cannot be solved for in closed form:

$F_1(x_1,x_2,\beta)=0 ~\&~ F_2(x_1,x_2,\beta)=0 $

I want to solve for functions $x_1=x_1(\beta) ~\&~ x_2=x_2(\beta)$

My approach so far:

Fix $\beta=\beta_1$

Solve system numerically for $x^1_1,x^1_2$.

Repeat for a different $\beta$ on some pre-defined finite support

Generate a vector of $\mathbf{x_1}=(x^1_1...x^n_1)$ with corresponding  $\mathbf{b}=(\beta_1...\beta_n)$

Approximate the relationship between $x_1$ and $\beta$ using Sieve to get $\hat{x}_1(\beta)$ [repeat for $x_2$]

The problem with this method is that it gets computationally out of hand once I have many equations and parameters due to step 3.

I do know that the system has a solution for any $\beta$, which motivated this method. I also can easily get the Jacobian of $F$ analytically.

Is there a preferable method that is not as computationally burdensome?


问:
FEM/FVM/FD for structural modeling and stability issues due to large structural constants?
I've read that in modeling structures problems, the finite element method (FEM) is typically used. I am unfamiliar with FEM, but I am wondering, in particular, if using FEM, as opposed to finite difference/volume, alleviates some of the numerical stability issues that may arise due to the large structural constants?

I am currently trying to solve the transient linear elastic equations (see https://en.wikipedia.org/wiki/Linear_elasticity#Direct_tensor_form) using FVM. For ease of demonstration, we can consider the 1-D form (just the classical second order hyperbolic wave equation) of this equation in the absence of shear and assuming constant density:

$\frac{\partial^2 u}{\partial t^2} = \frac{2\mu + \lambda}{\rho} * \frac{\partial^2 u}{\partial x^2}$.

Where $u$ is the x-displacement, $\rho$ is the mass density and $\mu$ and $\lambda$ are Lame constants (see http://scienceworld.wolfram.com/physics/LameConstants.html). I solved this by turning the second derivative equation into two first order equations:
Let $$\frac{\partial^2 u}{\partial t^2} = \frac{dz}{dt}$$
where
$$z = \frac{du}{dt}$$
or in vector form:

$$\frac{d\boldsymbol{q}}{t} = \boldsymbol{g}$$
where $\boldsymbol{q} = [u \ \ \ \ \  q]$ and $g = [z \ \  \ \ \ \frac{2\mu+\lambda}{\rho}* \frac{\partial^2u}{\partial x^2}=f $]

In preliminary modeling, I used Explicit Euler and set my constants = 1. The code "appeared" to be stable. Then I set the constants to be representative of steel and some other materials of interest, where the Lame constants fall between O(10^8)-O(10^11). This severely limited my the timestep size that I could use for my code.

Due to this, I decided to try to implement an implicit scheme. To do this, I tried to formulate my problem using the Jacobian matrix as such
$$
\frac{d\boldsymbol{q}}{dt} = A\boldsymbol{q}
$$
where A is Jacobian matrix
$$
A = \frac{d\boldsymbol{g}}{d\boldsymbol{q}} = \begin{bmatrix}
0 & I \\
\frac{df}{du} & 0
\end{bmatrix}
$$

Then I found the eigenvalues of this Jacobian matrix and found eigenvalues with a POSITIVE REAL root, indicating that this problem can't be stable (regardless of the large order of magnitude of the structural constants) at least with most conventional time integration schemes.

Earlier, when I stated that my explicit Euler scheme "appeared" to be stabled. I realized that it was mildly unstable. I had used very small structural constants such that the eigenvalues hugged the imaginary axis and the real roots magnitudes are smaller, but still positive.

Going back to my first question: Does FEM alleviate some of these mentioned issues? Second question: Does anyone have insight into solving such a problem where the jacobian matrix's eigenvalues have positive real roots? Third question: Any recommended schemes for tackling linear elasticity with FV that won't be severely limited by stability with EXPLICIT schemes?


问:
Simulation of a lens, insufficient points
I am simulating the propagation of a light pulse using the equation
$$\frac{\partial}{\partial z}A=\frac{1}{2\cdot k_0}\nabla^2_rA$$
with
$$k_0=\frac{2\pi}{\lambda_0}$$
The propagation with a step size of $dz$ is done using the equation
$$\left(\vec{1}-i\cdot\frac{dz}{2}\nabla^2_r\right)A_{n+1}=\left(\vec{1}+i\cdot\frac{dz}{2}\nabla^2_r\right)A_n$$
with
$$A_0(\vec{r})=\exp\left(-\frac{r^2}{2}\right)$$
Now I add a lens with a focal length of $f$, which is done using the equation
$$A_{0,L}=A_0\exp\left(-i\frac{\left\vert\vec{r}\right\vert^2}{2f}\right)$$
My problem is now: I discretize $A$ along a vector $\vec{r}$ from $0$ to $r_{max}$ with $r_{num}$ elements (radial approach). That works for $F>F_{crit}$. If $F<F_{crit}$, the oscillations of the lens factor are not properly resolved anymore, resulting in artifacts in the simulation.

For my current input values I need at least more than $50000$ points, but even that is not enough due to a small value for $F$.

Using that amount of points increases the amount of computation time and (when displaying the moving pulse) the amount of storage. Is there any way to decrease the amount of points, while still being able to resolve the lens oscillations? I am already using sparse matrices for the calculation, but for storage I have to use a dense matrix.


问:
Implementing Frechet Distance on pandas data frame
I am working on a trajectory analysis project using python and its data science related libraries.

I am planning to implement Frechet Distance algorithm for trajectory analysis, each trajectory has sequence of (x,y) coordinates along with timestamp, speed, dist between consecutive points, etc.

Following is the example df with "voygID" and respective x, y coordinates. (I am not including timestamp, speed, etc columns here, as I assume it wont be required in implementing Frechet distance)

```<code>     voygID         x             y
 0   1             -7.935513     5.103579
 1   1             -7.935781     5.103300
 2   1             -7.936354     5.102726
 3   1             -7.935915     5.102802
 4   2             -7.935306     5.103424
 5   2             -7.945678     5.119876
 6   2             -7.954764     5.128738
 7   3             -7.954888     5.138898
 8   3             -7.955897     5.149808
 9   3             -7.965789     5.156789
 10  3             -7.983457     5.198610
</code>```

In the above data frame (x,y) forms the coordinate data and all the points belonging to the same "voygID" forms one separate trajectory; Therefore in the above df rows(0-3) belonging to voygID = 1 is one trajectory, while points belonging to voygID = 2 is another trajectory and so on. There are thousands of trajectories each with 100s of points.

Now, I want to implement Frechet Distance to compare each trajectory with another and group the similar curves/ trajectories together and form different clusters.

I would kindly request if anyone who has experience in this field, can provide with some useful links for Python implementation of Frechet Distance to do the above task.

Also, I am open to other similarity measure methods like "Dynamic Time Warping" or "Longest common subsequence". Please suggest me if there are preprogrammed python implementations are available.


问:
How to generate a tet10 mesh from a tet4 mesh
I am writing a FEN solver on tet10 elements. I right now have a CFD grid gen tool (Pointwise) which generates tet4 meshes. From this I have the element connectivity list and the node list. How do I go about generating a tet10 mesh from this info.  Is there any particular tool/code to do this. </p>


问:
Implementing Neumann boundary condition for elasticity problem using the finite element method
I am solving the following mixed boundary elasticity problem with the finite element method on the unit circle, using a triangle mesh with 3 nodes on each element. In my problem I know $u_1$  and $t_2$:

With the following weak form:

$$\int_\Omega \sum_{i,j}^2 \sigma_{ij}(u) \epsilon_{ij}(v)= \int_\tau v_2 t_2 $$

I have already implemented it correctly for the pure Dirichlet problem. (so the global stiffness matrix is also good for this problem), but I do something wrong when implementing the Neumann boundary condition and I cannot figure out what. (I still get the right solution for $u_1$ where I have the Dirichlet condition). On the right side of the picture is the fem solution and on the left side is the analytical one:

This is what I do when implementing the Neumann boundary condition:

I go over each triangular element:

I extract the coordinates of each node of the element:
$(x_1,y_1), (x_2,y_2), (x_3, y_3)$
(the elements are numbered in anticlockwise sense)

if $(x_1,y_1)$ and $(x_2,y_2)$ are on the boundary:

I calculate the outward unit normal:
$$ n2=\frac{x_1-x_2}{\sqrt{(y_2-y_1)^2+(x_1-x_2)^2}} $$

I update the force vector F in the place corresponding to $C_e^12$ and $C_e^22$ ($v_2=\sum_{B=1}^3 N^BC_e^B2$) by adding :
$$ 0.5 \sigma_0 n_2 \sqrt{(x_2-x_1)^2+(y_2-y_1)^2}$$

if $(x_1,y_1)$ and $(x_3,y_3)$ are on the boundary:

$n2=\frac{x_3-x_1}{\sqrt{(y_1-y_3)^2 + (x_3-x_1)^2}}$

update with
$ -0.5 \sigma_0 n_2 \sqrt{(x_3-x_1)^2 + (y_3-y_1)^2} $

if $(x_2,y_2)$ and $(x_3,y_3)$ are on the boundary:

$ n2=\frac{x_2-x_3}{\sqrt{(y_3-y_2)^2+(x_2-x_3)^2}} $

update with:
$ 0.5 \sigma_0 n_2 \sqrt{(x_3-x_2)^2 + (y_3-y_2)^2} $


问:
energy computation for BVP with Dirichlet boundary conditions
I am solving quadratic minimization problem
\begin{align}
\min_{x}\ \frac{1}{2} x^T A x -b^T x,
\end{align}

where matrix A results from discretization of Laplacian by FEM method, subjected to Dirichlet BC.

I am having concerns about computing energy.
Should product $ \frac{1}{2} x^T A x -b^T x $ be computed by using matrix A with encoded modification for given BC (zero-ing off diagonals), or should it be done directly on unmodified (pure Neumann) matrix A ?

Is there any special treatment of boundary conditions needed ?


问:
How and when to use IMEX runge-kutta
i have come across method of IMEX Runge-Kutta written by Kennedy and Carpenter in july 2001. The method consist of two kind that is "implicit and explicit"

can i use the ERK@explicit part only to solve the intial value problem, as if i treat the method as another embedded Runge-Kutta?
is it against the algorithm/rule of the method?

Title: How and when to use IMEX runge-kutta


问:
Numerical implementation of the Dirichlet-to-Neumann map
I am solving the Dirichlet problem
$$
  \begin{cases}
   \Delta u = 0, \\
   u|_{\partial D} = f,
  \end{cases}
$$
in a $2d$ domain $D$ using the finite element method. What I want to get is the normal derivative of the solution $u$ on the boundary: $\tfrac{\partial u}{\partial \nu}$. In other words, I'm interested in FEM implementation of the Dirichlet-to-Neumann (Poincare-Steklov) map $\Lambda$. In Wikipedia I've found the following:

When the partial differential equation is discretized, for example by
  finite elements or finite differences, the discretization of the
  Poincaré–Steklov operator is the Schur complement obtained by
  eliminating all degrees of freedom inside the domain.

Hovewer, I could not find an explicit answer to this problem in Google. I'm programming in Matlab and I would extremely appreciate an optimized solution.

Update (30 May 2017). As it was pointed out by Praveen Chandrashekar in comments, I tried to solve this problem by formulating it in a weak form: find $h \in L^2(\partial D)$ such that for any $w \in H^{\frac 1 2}(\partial D)$,
$$
   \int_{\partial D} hw \,dl = \int_D \nabla w \nabla u dS.
$$
As usually, I represent $h = \sum' h_j \psi_j$, $w = \sum w_i \psi_i$, $u = \sum u_l \psi_l$, where $\psi_i$ are the basis elements and $\sum'$ denotes summation over the boundary elements only, and I get
$$
  \sum_i' \sum_j' h_j w_i \int_{\partial D} \psi_i \psi_j dl =  \sum_k \sum_l w_k u_l \int_D \nabla \psi_k \nabla \psi_l dS.
$$
Since $\{w_j\}$ are arbitrary, we get the system
$$
   \sum_j' \int_{\partial D} \psi_i \psi_j dl \cdot h_j = \sum_l\int_D \nabla \psi_k \nabla \psi_l dS \, \cdot u_l.
$$
Introducing the matrices
$$
   B_{ij} = \int_{\partial D} \psi_i \psi_j \, dl, \quad A_{ij} = \int_D \nabla \psi_i \nabla \psi_j dS,
$$
where in the first matrix $i$, $j$ correspond to boundary elements and in the second matrix $i$ corresponds to the boundary elements and $j$ corresponds to all elements, I get the matrix equation
$$
  B h = Au,
$$
for finding $h$.

However, if I implement this approach, the results are far from being precise. For example, consider the following circular domain and boundary values $u_0(x,y) = -y$ (triangulation created using DistMesh):

Then the flux along the boundary must coincide with $u_0(x,y)$ but I get the following noisy picture:

What am I doing wrong?


问:
Python program to simulate the trajectory of a cannonball
I am new into programming. Here, I am trying to simulate the motion of a cannonball until it hits the ground using python programming. I have used both analytical and numerical(Euler method) to plot its trajectory. But, I am not sure if my output is correct for both methods, especially for the analytical solution. For analytical solution, is there any way to show that for t=0, x=1m which is true.
Please look at the code below, and point out my mistakes.</p>

```<code>import matplotlib.pyplot as plt
v=0.0      #Initial 0 velocity
h=0.01     #dt
x=1        #Canonball is at a distance of 1m from the ground initially
t=0.0
m=2.0
g=9.81

ta,xa,xb=[],[],[]   #xa for numerical and xb for analytical solution

while x>=0.0:
    ta.append(t)
    xa.append(x)
    xb.append(-0.5*g*t**2.0)  #Analytical solution
    F=-m*g                       #Numerical solution
    v=v+(F/m)*h
    x=x+v*h
    t=t+h

plt.figure()
plt.plot(ta,xa)
plt.plot(ta,xb,'ro-')
plt.show()
</code>```


问:
Physics of explosions: just vorticity?
eg 1 2. Is it just vorticity? What's actually happening? (Similar: steam engines, volcanoes, clouds).

examples are grid-based, using "vorticity confinement" in Phoenix FD.

EDIT Some techniques: "divergence control" is used to increase the smoke concentration; and fuel particles in the air may combust, further adding smoke at that point (i.e. not at the initial source of the smoke).

A reference: Feldman 2003, Animating Suspended Particle Explosions, ACM Trans. Graph. (Proc. SIGGRAPH) 22 (2003), 708-715

Title: Physics of explosions: just vorticity?


问:
Best way to conduct parallel scaling tests?
I have just finished a new version of a parallel code that I have been working on. I would like to do a strong scaling test, but I am always bothered by how much monotonous work scaling tests seem to require.

Currently, I make folders for each processor count (np1, np2, np4, ...), copy my code into each folder, assign the correct parameters to the PBS file for submission to the cluster, and then submit the jobs one-by-one.

After job completion, I copy the data from each folder onto my machine and peice-by-piece collect the timing data from each folder. Then it goes into MATLAB.

I am curious if any of you have developed more efficient ways to run scaling tests. Do you write scripts to do most of this for you? If so, how do you accomplish that? Bash scripting? Perl, perhaps?

Thanks!


问:
How to get a theoretical background in nonlinear, coupled FEM systems
I'm currently developing simulations for coupled, nonlinear, multi-region systems. Basically, I use the Finite Element Method (FEM) to model each physical quantity in each region. The obtained matrices and the coupling matrices between the quantities are combined into one large system matrix. From this matrix plus the derivatives of the nonlinear functions, a Jacobi matrix is generated. The system is then finally solved by the Newton-Raphson method. This includes the inversion of the Jacobi matrix, which is either done iteratively or directly.

I do not have a specific question, but I observe that I'm running into problems over and over again. Bad convergence, instabilities, oscillations, to name some of them. It's proably because I have almost no theoretical background in the field. FEM itself is not the problem, but the solution of the equations in a coupled and nonlinear scenario.

I'm looking for some good material to get some basic theoretical understanding. What's a good introductory book or lecture about the topic? (I'm a graduated engineer, so at least the basic maths is - hopefully - not the problem.)


问:
Beam propagation over extremely short distances
I like to simulate the free-space propagation of an electric field (beam radius ~ 1-10 cm) over extremely short distances (~ 5-10 cm). The reason is that the field is focused by a lens with a very small focal distance (~ 30 cm). For larger distances, the FFT based method works but it fails for short distances due to aliasing (and I tried extremely high spatial resolutions without success). So I tried to compute the Rayleigh-Sommerfeld diffraction integral directly, but the highly oscillating term in it demands an extremely high spatial resolution as well once the propagation distance becomes too short. The only other method I can think of at the moment is a finite-differences solver for the wave equation itself. But before I try to implement it, I thought about asking if anyone has experience with this kind of problem? And maybe even some advice what I can try to solve it?</p>

OwnerDisplayName: DaPhil