问:
A better way to compute a double integral involving a infinite series?
Let $D_{\nu}(.)$ is the parabolic cylinder function

(http://mathworld.wolfram.com/ParabolicCylinderFunction.html)

And $\Gamma(.)$ is the Gamma function.

Define

$s_y(\mu,\nu,t,z)=2^{\nu}\displaystyle\sum_{k=0}^{\infty}\frac{\Gamma(\nu+k)e^{-(\nu t+z+2kt)^2/(4y)}}{\sqrt{2\pi}y^{1+\mu/2}\Gamma(\nu) k!} D_{\nu+1}(\frac{\nu t+z+2kt}{\sqrt{y}})$

and let

$$es_y(\mu,\nu,t,x,z)=\displaystyle\sum_{k=0}^{\infty}\frac{(-z)^k}{k!}
s_y(\mu+k,\nu+k,t,x+z+kt) $$

and let

$$is_y(\nu,t,r,z,x)=\displaystyle\sum_{l=0}^{\infty}\frac{x^{\nu+2l}}{\Gamma(\nu+l+1)l!}es_y(1+\nu+2l,1+\nu+2l,t,r,z) $$

Assume that $u,w$ are variables, the other variables are know. I would like to compute an integral like this
$$I=\int_{0}^{\infty}\int_0^{\infty}f(u,w) is_a(2\nu,cu,0,w+d,e\sqrt{w}) dudw$$

the integral $I$ can be approximated by truncating the series at large $n$.
When i did so, the result converges very slowly, (not very accurate).

My questions: is there a better way to approximate the above double integral ?
Do we have a program language that can use to find $s_y(),es_y(.),is_y(.)$ exactly ( I mean without truncating the series)?


问:
Hessian eigenvalues in 4D-VAR data assimilation
I am using variational data assimilation (4D-VAR) to estimate emissions of anthropogenic greenhouse gases using a rather complex atmospheric transport model. Hence, the optimal solution to my problem is obtained by minimizing a (non-linear) cost functional $J$.

The model is weakly nonlinear so I can expect the inverse of Hessian of the 4D-VAR cost functional - let's call it $\nabla^2 J^{-1}$ - at the optimal (analysis) point $x_a$ to be a decent approximation to the posterior error covariance matrix $P$. As per standard 4D-VAR theory, the leading eigenpairs of $\nabla^2 J$ (approximately) span the maximum uncertainty space after the assimilation. This all makes sense when we are at a local optimum - i.e., $\nabla^2 J$ is (semi-)positive definite, but what to do if we have not yet reached a local minimum and $\nabla^2 J$ and some of its leading eigenvalues are negative? Is it "safe" to exclude the corresponding eigenpairs from my approximation to the uncertainty space?

ContentLicense: CC BY-SA 3.0


问:
Linear vs Non Linear inverse problems: Does non-linearity help?
This is not a typical question with a deterministic answer. If this is not the right place, feel free to close it.

For the past one year I have been working on various kinds of inverse problem. Most of them are related to parameter estimation, ranging from small linear differential equation to equations of fluid dynamics. And consistently, I observed that it was easier to find a solution to the non linear problem(using same optimization method, let's say steepest descent). By easy I mean it requires lower number of iterations for the same decrease(very subjective). Although, it is nearly impossible to compare different cases exactly, this trend is surprising to me.

The question here is, did anyone else observe similar behaviour? Is there a qualitatively explanation for this behaviour?

This also raises a follow up question of whether we can quantify(from the initial jacobian may be?) the level of difficulty involve in finding a solution?(or estimate the order of magnitude of the number of optimizer steps)


问:
preconditioned Uzawa method with Petsc
I am trying to improve the resolution of a Stokes problem (P2/P1 on unstructured mesh) defined by the matrix $M$:

$M=
\begin{pmatrix}
A_u & 0 & B_u \\
0 & A_v & B_v\\
B_u^T & B_v^T & 0
\end{pmatrix}$

Current Implementation

I am currently solving it using a Uzawa method, with a Krylov method accelerated by a preconditioner. This preconditioner is build by additive composition of the matrices $A_p$ and $A_m$.

I use Petsc to solve each sub-problems, $A_u$, $A_v$, $A_p$ and $A_m$ (in format mpiaij), but I handle the Uzawa method separately.

$A_u$ and $A_v$ are solved using cg, ASM preconditioning and level 2 icc:

```<code>-ksp_type cg -pc_type asm -sub_pc_type icc -sub_pc_factor_levels 2
</code>```

$A_p$ is solved using cg with AMG preconditioning:

```<code>-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg
</code>```

$A_m$ is solved using cg without preconditioning:

```<code>-ksp_type cg -pc_type none
</code>```

Goal

To improve the speed and allow for easier tuning , I would like to have Petcs do the all the Uzawa process. I am able to pass $M$ in mpiaij format to Petsc. But the tool I use can only pass one, and only one matrix to Petsc, but I could concatenate $M$+$A_p$+$A_m$ and give that to Petsc:

$\begin{pmatrix}
M & 0 & 0\\
0 & A_p & 0\\
0 & 0 & A_m
\end{pmatrix}$

How can I tell Petsc to separate $M$ from $A_p$/$A_m$ and to use $A_p$ and $A_m$ as preconditioner matrices for the Shur complement in a Uzawa method?

Title: preconditioned Uzawa method with Petsc


问:
Solving an ODE using shooting method
I have been trying to solve the following nonlinear ordinary differential equation:

$$-\Phi''-\frac{3}{r}\Phi'+\Phi-\frac{3}{2}\Phi^{2}+\frac{\alpha}{2}\Phi^{3}=0$$

with boundary conditions

$$\Phi'(0)=0,\Phi(\infty)=0.$$

My solution is supposed to reproduce the following plots:

Now, to produce the plots given above, I wrote the following Mathematica code:

```<code>α = 0.99;
Φlower = 0;
Φupper = 5;
For[counter = 0, counter <= 198, counter++,
    Φ0 = (Φlower + Φupper)/2;
    r0 = 0.00001;
    Φr0 = Φ0 + (1/16) (r0^2) (2 (Φ0) - 3 (Φ0^2) + α (Φ0^3));
    Φpr0 = (1/8) (r0) (2 (Φ0) - 3 (Φ0^2) + α (Φ0^3));
    diffeq = {-Φ''[r] - (3/r) Φ'[r] + Φ[r] - (3/2) (Φ[r]^2) + (α/2) (Φ[r]^3) == 0, Φ[r0] == Φr0, Φ'[r0] == Φpr0};
    sol = NDSolve[diffeq, Φ, {r, r0, 200}, Method -> "ExplicitRungeKutta"];
    Φtest = Φ[200] /. sol[[1]];
    Φupper = If[(Φtest < 0) ||   (Φtest > 1.2), Φ0, Φupper];
    Φlower = If[(Φtest < 1.2) && (Φtest > 0),   Φ0, Φlower];
]
Plot[Evaluate[{Φ[r]} /. sol[[1]]], {r, 0, 200},
     PlotRange -> All, PlotStyle -> Automatic]
</code>```

In the code, I used Taylor expansion at $r=0$ due to the $-\frac{3}{r}\Phi'$ term. Moreover, I used shooting method and continually bisected an initial interval from $\Phi_{\text{upper}}=5$ to $\Phi_{\text{lower}}=0$ to obtain more and more precise values of $\Phi(0)$.

With the code above, I was able to produce the plots for $\alpha = 0.50, 0.90, 0.95,0.96,0.97$. For example, my plot for $\alpha = 0.50$ is as follows:

However, my plot for $\alpha = 0.99$ does not converge to the required plot:

Can you suggest how I might tackle this problem for $\alpha = 0.99$? Also, is there an explanation for the plots shooting upwards and oscillating after a prolonged asymptotic trend towards the positive $r$-axis?

ContentLicense: CC BY-SA 4.0

Title: Solving an ODE using shooting method


问:
Manipulating/Extracting Data and Developing Methods - Language Choice
As a general programming enthusiast and aspiring Bioinformatician student I have an intermediate understanding of computing (languages) as well as Java, and to a lesser extent C++. Having knowledge in a declarative as SQL and the intent to learn a functional as R; I found myself in need of interpreted/scripting language for quick data management tasks. Such would be used in conjunction with a compiled language for more complex algorithms. I am unbiased as I have yet to commit to mastering a single language; albeit I hope to make that decision soon and a major asset like scripting will play a role in that choice. Leading me to the following:

When paring a compiled and scripting language, which languages
fundamentally compliment the opposing features of the other to the greatest   degree while retaining the ability to communicate by call/implementing their use?

My intent in this case would be able to “plug” one into another to compensate an individual’s shortcomings (if not, at least conceptually) to holistically provide the greatest breath of benefit. Preference would of course be the design and features of the language as in contrast to more creative “hacking” solutions.

(IE Language II has the benefit… but by nature is poor/slow/lacks at… therefore Language A would suit it well.)

Considerations

Memory Management and the Compile-Time versus Execution-Time Tradeoff

Data Types, Generics/Reflection, and Inheritance

Platform/OS Compatibility (Including Source and Code Version)

Ease of use is not so much a consideration as in the case of learning and added utility, but is undeniably a factor in the case of collaboration. I would greatly appreciate any insight as well as the inclusion of anything else that I may have overlooked if applicable. Thanks in advance.


问:
CAS Problem with integrals
I got this problem thrown at me, unfortunately I lack context at the moment but I thought Maple or Mathematica would solve it anyways.

I have a function $f$ over $x$ and $y$ such as this (Maple)

```<code>f := (y-1.298*exp(-1.311*x)-.3435*exp(.3602*x))*(x-3.059865611)*(y+1.875*exp(-.8384*x)+0.7033e-1*exp(.838*x))*(y+.5774*x+3.533)*(.1291*x^2+0.475e-1*y^2+.6927*x+.7080)*(y-.5773*x-3.533)*c;
</code>```

I have boundary condition so that c needs to be determined to satisfy:

$$
\frac{d^2f}{dx^2}+\frac{d^2f}{dy^2} = -2
$$

```<code>solution := solve(diff(f, x, x)+diff(f, y, y) = -2, c)
fsubs := subs(c = solution, f)
</code>```

So far so good, now I want to integrate
$$
\int \left(\frac{df}{dy}-y\right) dx
$$

```<code>integrate(diff(fsubs, y)-y, x)
</code>```

this is where things no longer work. I may be impatient but I thought ten minutes would be enough.

Any ideas?

Title: CAS Problem with integrals


问:
Comparison between of higher order interpolations
A while ago I came up with an algorithm which can be used to numerically solve optimal control problems, which basically came down to discretizing the control input $u(t)$ and interpolating this to get a continuous function. I tried to optimize the cost function for each discrete value by using the golden search method and because of the interpolation this would mean I only have to integrate this over a small section in time. The fact that each discrete value of $u(t)$ only affects a small number of elements also allows to use parallel computing.

For the initial problem, I tried to optimize, I used linear interpolation, since the cost function was only a function of $u(t)$ and $u'(t)$. However recently I tried to expand on this, such that I would also be able to use higher order derivatives of $u(t)$. For linear interpolation you only need the discrete points at the start and end of an element,

$$
f_{i}(\xi) = u_{i} (1 - \xi) + u_{i+1} \xi,
$$

where the $u(t)$ from $t_i$ to $t_{i+1}$ is defined as,

$$
u(t) = f_i\left(\frac{t - t_i}{t_{i+1} - t_i}\right).
$$

In order to ensure that $u'(t)$ is continuous (such that second derivative exists) I had to set the slope at the end of an element equal to the slope at the beginning of the next element. I decided on defining this slope using the central difference method, which meant using two extra points, one before and one after the element,

$$
f'_i(0) = \frac{u_{i+1} - u_{i-1}}{2},
$$

$$
f'_i(1) = \frac{u_{i+2} - u_{i}}{2}.
$$

Using these additional constraints yields the following function for an element,

$$
f_i(\xi) = u_{i-1} \frac{-\xi^3 + 2\xi^2 - \xi}{2} + u_i \frac{3\xi^3 - 5\xi^2 + 2}{2} + u_{i+1} \frac{-3\xi^3 + 4\xi^2 + \xi}{2} + u_{i+2} \frac{\xi^3 - \xi^2}{2}.
$$

Instead of the central difference method I also could have used backward or forward difference method, which would require one less point to be used, but I would think that interpolation ought to be symmetrical. These points can also used in the central difference method for the second derivative. Using this as a constraint for the initial and final second derivative yields the following function for an element,

$$
f_i(\xi) = u_{i-1} \frac{2\xi^5 - 5\xi^4 + 3\xi^3 + \xi^2 - \xi}{2} + u_i \frac{-3\xi^5 + 15\xi^4 - 9\xi^3 - 2 \xi^2 + 2}{2} + u_{i+1} \frac{3\xi^5 - 15\xi^4 + 9\xi^3 + \xi^2 + \xi}{2} + u_{i+2} \frac{-\xi^5 + 5\xi^4 - 3\xi^3}{2}.
$$

However will this function have any disadvantages compared to the previous function? For example would the higher order polynomial make it harder to approximate sections of the cost function integral?


问:
What is the best way to stop an asynchronous-heterogeneous MPI program?
I want to stop immediately an MPI program, that has N processors and each processor runs a different algorithm (with different computational complexities).
Processors share information to each other using an unidirectional token ring topology. As soon as one processor ends it should notify its neighbor and all processors should stop immediately as soon as they receive the interruption.
Here is an sketch of the execution program at each processor:

```<code>1. Processor_i(array of n data) {
2.   while (stopping condition is not satisfied) {
3.     Run algorithm_i(data) for the ith processor
4.     Receive information from its source node (at most read 2 messages)
5.     if (the incomming message indicates interruption)
6.       break;
7.     Send its own information to its destiny node (maximum 2 messages)
8.   }
9.   Save local solution
10.  Send interruption message to its destiny node
11. }
</code>```

Notice that in the step 3, the computational complexity may be, for example, O(n^2), O(n log n) or O(n^3) for processors 1, 2 and 3, respectively.

My problem is that the messages are buffering by MPI and the interruption message is not reading immediately.
I was thinking in writing a file with the stopping instruction, but the filesystem is not shared.

How could I deal with this problem?


问:
Implementation of no-slip boundary conditions in lattice Boltzmann method fluid simulation
My faculty advisor recommended that I take a look at the lattice Boltzmann method as an introduction to scientific computing and potentially an undergraduate honors thesis topic. I cooked up a some Java code (kind of unconventional, but I only know Java and Python right now), and I had a simple 2D simulation. Upon testing a 2D channel flow, I discovered many things wrong with the behavior.

The velocities at the nodes are not stable. In general they all approach zero. By 10,000th time steps, they are on the order of $10^{-15}$

Despite the parabolic velocity profile one would expect, the magnitudes of the x-component of the velocity switches between positive and negative along the x-axis.

The reference I'm using explains the general algorithm clearly, but I have yet to find a clear description of the implementation of the no-slip boundary conditions. This is where I think my problem lies. I'm attempting to implement the on-grid boundary conditions (pgs. 3-4) and simulate the Poiseuille flow described in article.

How exactly are the no-slip boundary conditions initialized? Right, now I'm operating under the assumption that you calculate the equilibrium functions based on initial density and velocity.

There may be an even larger problem with my code if this isn't the problem. The source code is located here if anyone is curious. Keep in mind this was the most intuitive approach to implementing the LBM for me and not necessarily to most efficient.


问:
non convex, non linear optimization involving matrix differential equation solution
I'm trying to develop an inferential procedure for a multivariate dependent  Markov process. Basically, the procedure could be considered as a non linear regression, with a known dependence structure among observation belonging to the same time-point and independent from the others. The non linear trend correspond to a the solution of a matrix differential equation. A more detailed description follows:

Let define a matrix $\Theta$ of parameters $\theta_i$
\begin{array}{ccc}
\theta_1 & \theta_4 & \theta_7 \\
\theta_2 & \theta_5 & \theta_8 \\
\theta_3 & \theta_6 & \theta_9  \end{array}
each entry $\theta_i$ could be unconstrained or linearly constrained (both equality and inequality).
This matrix, multiply by a  given, fixed matrix $V$ and a vector of previous time point $\vec{X}$ observation, governs the derivative of a stochastic process.
As a consequence, in my objective function $\sum(\vec{Y}-f(\Theta,\vec{X}))^2$, I need to calculate the solutions of the matrix differential equation, $f(\Theta,\vec{X})$. To calculate the solutions, I need to calculate a matrix $A$, combining the entries of $\Theta$ according to $V$, calculate eigen decompositions of $A$ and finally the solution $\hat Y=f(\Theta,\vec{X})$ (exponential).

The objective function is non-convex and has many local minima. Runnig some simulations studies, I verified that the global minimum is located at the true parameters values, and the objective is convex in the closed neighbourhood.
Until now, I used a constrained Gauss-Newton method to estimate parameters, and, given good initial values, it works. Relaxing some conditions, the method to calculate initial values is not good enough to guarantee the Gauss-Newton method to converge to global minimum (it get stucks in a closed local minimum). So, I'm now looking for a global optimizer. All my code is in R and RCplex and I have some experience with C++. Is Cplex able to solve non-convex,linearly constrained problem and flexible enough to allow me for solve the matrix differential equation? Are there any alternative program? Any suggestion?
Thanks.


问:
What should be the number of boundary conditions of a PDE
As far as I know, for getting a unique solution to a PDE we should impose some boundary conditions to the PDE. "The number of required auxiliary conditions is determined by the highest order derivative in each independent variable. "

My questions are: In any numerical method, especially in finite element method, is it valid:

1-) to use more boundary conditions than the highest order derivative in each independent variable.

2-) to use less boundary conditions than the highest order derivative in each independent variable.

For example:

1-) In this article and in this article, while the highest order derivative of a spatial variable is in third order authors used two boundary conditions for that variable.

2-) In this article while the highest order derivative of a spatial variable is in third order authors used four boundary conditions for that variable.

I can increase the number of articles like this.


问:
Load balancing/partitioning with unknown weights
For a grid-based numerical simulation, I am looking for a load balancing/partitioning algorithm that not only distributes my grid elements, but also determines (approximates) their respective weights. Does anyone know of existing approaches for this problem or can give me some pointers into which kind of mathematical field I should be looking for a possible solution?

Details

Let's say I have $N_e$ grid elements $e_i$, each with a certain computational weight $w_i$. The $w_i$ are all within an order of magnitude, but I don't know their exact value a priori. The elements are now to be distributed among $N_p$ processes, preserving their order. I can only measure the cumulative weight of all elements on one process but not their individual contributions. However, I can repartition the grid multiple times to get multiple measurements with different decompositions.

What I am searching for is a method that lets me determine an estimate for the $w_i$ with a relatively low number of decompositions $N_d$. The numbers involved  here are in the range of $N_e = \mathcal{O}(10^9)$, $N_p = \mathcal{O}(10^5)$, and $N_d = \mathcal{O}(10^2)$.


问:
BiCopter simulation in Matlab
I derived a bicopter dynamical model with two servos and two BLDC motors. And now am trying to simulate it using Matlab.
As base for simulation I used this paper and this code

Unfortunately, the obtained model always blowing up in case of asymmetrical control input signals. For example, with identical speed on both propellers and 5 degrees tilting angle for one of the thrust vectors, the system is blowing up in less the half second.

Here is simulation plots:
dt = 0.0001, first 1300 points.

This result obtained for constant 5 degree tilting angle for one of the motors and equal motor speeds. Angular acceleration just grows to the infinity in less than 1 second.

Matlab code:

```<code>% Perform a simulation of a copter, from t = 0 through t = 1.
% As an argument, take a controller function. This function must accept
% a struct containing the physical parameters of the system and the current
% gyro readings. The controller may use the strust to store persistent state, and
% return this state as a second output value. If no controller is given,
% a simulation is run with some pre-determined inputs.
% The output of this function is a data struct.

function result = simulate(controller, tstart, tend, dt)
    % Physical constants.
    g = 9.81;
    m = 0.95;

    d = 0.35;
    h = 0.08;

    I_body = diag([5e-3, 5e-3, 10e-3]);
    I_react=1e-3;
    I_gyro=1e-3;

    I = struct('I_body', I_body, 'I_react', I_react, 'I_gyro', I_gyro);

    kd = 0.25;

    kt = 0.00001*9.8;
    kf = 0.0015*9.8;

    % Simulation times, in seconds.
    if nargin < 4
        tstart = 0;
        tend = 1;
        dt = 0.0001;
    end
    ts = tstart:dt:tend;

    % Number of points in the simulation.
    N = numel(ts);

    % Output values, recorded as the simulation runs.
    xout = zeros(3, N);
    xdotout = zeros(3, N);
    thetaout = zeros(3, N);
    thetadotout = zeros(3, N);

    inputout_beta = zeros(2, N);
    inputout_motor_w = zeros(2, N);

    inputout_beta_speed = zeros(2, N);
    inputout_beta_acc = zeros(2, N);

    inputout = struct('inputout_beta', inputout_beta, 'inputout_beta_speed', inputout_beta_speed ,...
        'inputout_beta_acc', inputout_beta_acc,'inputout_motor_w', inputout_motor_w);

    % Struct given to the controller. Controller may store its persistent state in it.
    controller_params = struct('dt', dt, 'I', I, 'kf', kf, 'd', d, 'h', h, 'm', m, 'g', g);

    % Initial system state.
    x = [0; 0; 10];
    xdot = zeros(3, 1);
    theta = zeros(3, 1);

    % If we are running without a controller, do not disturb the system.
    if nargin == 0
        thetadot = zeros(3, 1);
    else
        % With a control, give a random deviation in the angular velocity.
        % Deviation is in degrees/sec.
        deviation = 300;
        thetadot = deg2rad(2 * deviation * rand(3, 1) - deviation);
    end

    ind = 0;
    for t = ts
        ind = ind + 1;

        % Get input from built-in input or controller.
        if nargin == 0
            i = input(t); %built-in generator
        else
            [i, controller_params] = controller(controller_params, thetadot);
        end

                        % Store simulation state for output.
        % don't go through the ground
            if (x(3) < 0)
                x(3) = 0;
                theta(1) = 0; % avoid integration windup
                theta(2) = 0;
            end

        %%%%% Compute forces, torques, and accelerations.
        omega = thetadot2omega(thetadot, theta);
        a = acceleration(i, theta, m, g, kf, kt, kd, xdot);
        omegadot = angular_acceleration(i, omega, I, d, h, kf, kt);

        %%%% Advance system state. Euler method.
        omega = omega + dt * omegadot;
        thetadot = omega2thetadot(omega, theta);

        theta = theta + dt * thetadot;
        xdot = xdot + dt * a;
        x = x + dt * xdot;

        linacc(:, ind) = a;
        angacc(:, ind) = omegadot;a
        xout(:, ind) = x;
        xdotout(:, ind) = xdot;
        thetaout(:, ind) = theta;
        omegaout(:, ind) = omega;
        inputout(:, ind) = i;

    end
    % Put all simulation variables into an output struct.
    result = struct('linear_coord', xout, 'angular_coord', thetaout, 'linear_velocity', xdotout, ...
                    'angular_velocity',omegaout ,'angular_acc', angacc, 'linear_acc', linacc,'t', ts, 'dt', dt, 'input', inputout);

end

% Arbitrary test input.
function in = input(t)
inputout_beta= zeros(2, 1);
inputout_motor_w= zeros(2, 1);

    inputout_beta(1) = 0.15;%-0.1*cos(t);
    inputout_beta(2) = 0;%0.1*cos(t);
    inputout_motor_w(1)= 3000*pi/30;%100*pi/30*cos(t) %rpm to radian per second
    inputout_motor_w(2)= 3000*pi/30;%101*pi/30*cos(t) %rpm to radian per second

    inputout_beta_speed(1)=0;
    inputout_beta_speed(2)=inputout_beta_speed(1);
    inputout_beta_acc = zeros(2, 1);

    in = struct('inputout_beta', inputout_beta, 'inputout_beta_speed', inputout_beta_speed ,...
        'inputout_beta_acc', inputout_beta_acc,'inputout_motor_w', inputout_motor_w);
end

% Compute thrust given current inputs and thrust coefficient.
function F = force(w_motor, beta_motor, kf)
     w_l=w_motor(1);
     w_r=w_motor(2);
     beta_l=beta_motor(1);
     beta_r=beta_motor(2);

    %Right power plant force
    f_r=kf*w_r;
    %Left power plant force
    f_l=kf*w_l;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Left motor force rot mat
    R_Pl_Bm=[cos(beta_l), 0, sin(beta_l);
          0,          1,       0;
         -sin(beta_l),0, cos(beta_l)];
    %Right motor force rot mat
    R_Pr_Bm=[cos(beta_r), 0, sin(beta_r);
          0,          1,       0;
         -sin(beta_r),0, cos(beta_r)] ;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    F_l=[0; 0; f_l];
    F_r=[0; 0; f_r];
    F=R_Pl_Bm*F_l+R_Pr_Bm*F_r;
end

% Compute torques, given current inputs, length, drag coefficient, and thrust coefficient.
function T_Bm = torque(input, d, h, k_f, k_d, I)
     w_l=input.inputout_motor_w(1);
     w_r=input.inputout_motor_w(2);

     beta_l=input.inputout_beta(1);
     beta_r=input.inputout_beta(2);

     beta_l_speed=input.inputout_beta_speed(1);
     beta_r_speed=input.inputout_beta_speed(2);

     beta_l_acc=input.inputout_beta_acc(1);
     beta_r_acc=input.inputout_beta_acc(2);

    %Right power plant torque
    tau_drag_r=k_d*w_r;
    %Left power plant torque
    tau_drag_l=k_d*w_l;

    %Right power plant force
    f_r=k_f*w_r;
    %Left power plant force
    f_l=k_f*w_l;

    F_l=[0; 0; f_l];
    F_r=[0; 0; f_r];

    OlOb=[0;d;h];
    OrOb=[0;-d;h];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Left motor force rot mat
    R_Pl_Bm=[cos(beta_l), 0, sin(beta_l);
              0,          1,       0;
             -sin(beta_l),0, cos(beta_l)];
    %Right motor force rot mat
    R_Pr_Bm=[cos(beta_r), 0, sin(beta_r);
             0,          1,       0;
             -sin(beta_r),0, cos(beta_r)] ;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     I_gyro = I.I_gyro;
     I_react = I.I_react;

     % TORQUES
        %Gyroscopic moments, longitudinal tilting
        Mgyro_Pl_beta_l=[0;-I_gyro*w_l*beta_l_speed;0];
        Mgyro_Pr_beta_r=[0;I_gyro*w_r*beta_r_speed;0];
        %Gyroscopic momentin body frame
        Mgyro_Bm=R_Pl_Bm*Mgyro_Pl_beta_l+R_Pr_Bm*Mgyro_Pr_beta_r;

        %Fan torques due aerodynamic drag of the propeller
        Mdrugr_Pr=[0;0;-tau_drag_r];
        Mdrugl_Pl=[0;0;tau_drag_l];
        %Drag torques in body frame
        Mdrug_Bm=R_Pl_Bm*Mdrugl_Pl+R_Pr_Bm*Mdrugr_Pr;

        %Moment from trust vectoring in body frame
        M_thrv_Bm=cross((R_Pl_Bm*F_l),OlOb)+cross((R_Pr_Bm*F_r),OrOb);

        %Adverse Inertial Reaction moment
        M_react_Pr=[0;-I_react*beta_r_acc;0];
        M_react_Pl=[0;-I_react*beta_l_acc;0];
        %Adverse Inertial Reaction moment, body frame
        M_react_Bm=R_Pr_Bm*M_react_Pr+R_Pl_Bm*M_react_Pl;
        %FULL TORQUE
        %Body frame
        T_Bm=Mgyro_Bm+Mdrug_Bm+M_thrv_Bm+M_react_Bm;
end

% Compute acceleration in inertial reference frame
% Parameters:
%   g: gravity acceleration
%   m: mass of quadcopter
%   kf: thrust coefficient
%   kt: global drag coefficient
function A = acceleration(inputs, angles, m, g, kf, kt, kd, linear_velocity)
    F_g = [0; 0; -g];
    R_Bm_Wi = rotation(angles);
    F_Bm = force(inputs.inputout_motor_w, inputs.inputout_beta, kf);
    %Friction as a force proportional to the linear velocity in each direction.
    F_a = -kd * linear_velocity;
    %mV=R*F
    A = R_Bm_Wi * F_Bm * 1/m + F_g + F_a;
end

% Compute angular acceleration in body frame
% Parameters:
%   I: inertia matrix
function omegad = angular_acceleration(inputs, omega, I, d, h, kf, kt)
    T = torque(inputs, d, h, kf, kt, I);
    omegad = I.I_body\(T - cross(omega, I.I_body * omega));
end

% Convert derivatives of roll, pitch, yaw to omega.
function omega = thetadot2omega(thetadot, angles)
    phi = angles(1);
    theta = angles(2);
    psi = angles(3);
    W = [1,   0,                   -sin(phi);
         0    cos(phi),    cos(phi)*sin(phi);
         0   -sin(phi), cos(theta)/cos(phi)];
    omega = W * thetadot;
end

% Convert omega to roll, pitch, yaw derivatives
function thetadot = omega2thetadot(omega, angles)
    phi = angles(1);
    theta = angles(2);
    psi = angles(3);
    W = [1,     sin(phi)*tan(theta),cos(phi)*tan(theta);
         0      cos(phi),                     -sin(phi);
         0      sin(phi)/cos(theta), cos(phi)/cos(phi)];
    thetadot = W * omega;
end
</code>```

Newton-Euler equations:

$$
\ddot{x}  m=[K_f  w_l  c(\beta_l) + K_f  w_r  c(\beta_r)]  [s(\phi)  s(\psi) + c(\phi)  c(\psi)  s(\theta)] + c(\psi)  c(\theta)  [K_f  w_l  s(\beta_l) + K_f  w_r  s(\beta_r)]
$$

$$
\ddot{y}  m =c(\theta)  s(\psi)  [K_f  w_l  s(\beta_l) + K_f  w_r  s(\beta_r)] - [K_f  w_l  c(\beta_l) + K_f  w_r  c(\beta_r)]  [c(\psi)  s(\phi) - c(\phi)  s(\psi)  s(\theta)]
$$

$$
\ddot{z}  m = c(\phi)  c(\theta)  (K_f  w_l  c(\beta_l) + K_f  w_r  c(\beta_r)) - s(\theta)  (K_f  w_l  s(\beta_l) + K_f  w_r  s(\beta_r)) - g  m
$$

$$
I_{xx}  \ddot{\phi}+\dot{\theta}  \dot{\psi}  (I_{zz}-I_{yy} )=K_t  w_l  s(\beta_l) - K_t  w_r  s(\beta_r) + K_f  d  w_l  c(\beta_l) - K_f  d  w_r  c(\beta_r)
$$

$$
I_{yy}  \ddot{\theta}+\dot{\psi}  \dot{\phi}  (I_{xx}-I_{zz} )=I_g  \dot{\beta_r}  w_r - I_r  \ddot{\beta_r} - I_g  \dot{\beta_l}  w_l - I_r  \ddot{\beta_l} - K_f  h  w_l  s(\beta_l) - K_f  h  w_r  s(\beta_r)
$$

$$
I_{zz}  \ddot{\psi}+\dot{\theta}  \dot{\phi}  (I_{zz}-I_{yy} )=K_t  w_l  c(\beta_l) - K_t  w_r  c(\beta_r) - K_f  d  w_l  s(\beta_l) + K_f  d  w_r  s(\beta_r)
$$

Where:
$d,h$ - distances to CG

$I$- Moments of inertia

$k_f,k_t$ - experimentally obtained coefficients

Problem is independent from simulation step used. I also tried RK4 realization with the same result.
How can I correct this behavior?
What is the cause of the singularity?

Title: BiCopter simulation in Matlab


问:
Jacobi method converging then diverging
I am working to solve Poisson's equation in 2D axisymmetric cylindrical coordinates using the Jacobi method. The $L^2$ norm decreases from $\sim 10^3$ on the first iteration (I have a really bad guess) to $\sim 0.2$ very slowly. Then, the $L^2$ norm begins to increase over many iterations.

My final matrix is weakly diagonally dominate, except for the 2nd order Neumann condition at $r = 0$.

Can I make a small tweek to make this work, is it numeric or do I need a new algorithm?

My geometry is parallel plates with sharp points at $r = 0$ on both plates.

My boundary conditions are
$$\left. \frac{\partial V}{\partial r} \right|_{r=0} = 0$$

Although I would like my second radial BC to be
$$\left. \frac{\partial V}{\partial r} \right|_{r=\infty} = 0$$
I settled for
$$\left. \frac{\partial V}{\partial r} \right|_{r=a} = 0$$

Then Dirichlet conditions at the upper and lower boundaries
$$V(r, L(r) ) = V_0$$
$$V(r, U(r) ) = V_L$$

where
$$L(r) = \begin{cases}
 & 0 \text{ if } r \geq R_L \\
 & H_L (1 - \frac{r}{R_L} ) \text{ if } r \leq R_L
\end{cases}$$

and

$$U(r) = \begin{cases}
 & H \text{ if } r \geq R_U \\
 & H + H_U (\frac{r}{R_U} - 1 ) \text{ if } r \leq R_U
\end{cases}$$

Title: Jacobi method converging then diverging


问:
What is the most complete open and linked dataset to use for taxonomy?
I want to ask questions like 'what family is this species in, and which other species are in this family?'.

I am aware of Wikispecies/Wikidata but are there any other good open and linked databases/APIs that would be worth looking at?


问:
Discontinuity at Interface
The equation at the left of the interface is

\begin{equation}
\displaystyle\frac{\partial C_i}{\partial t} = D_i \nabla^2 C_i - z_i \frac{D_i}{RT}F \nabla \cdot (C_i \nabla \phi_2)
\end{equation}

The equation at the right of the interface is

\begin{equation}
\displaystyle\frac{\partial \epsilon C_i}{\partial t} = \nabla \cdot \biggl(\epsilon D_i \nabla C_i + z_i \frac{ D_i}{RT}F\epsilon C_i \nabla \phi_2 \biggr) + S(x)
\end{equation}

with an extra source term accounted for pore wall reaction and precipitation (also an extra dependent variable $\epsilon$ for porosity).

At $t=1e-5$ the concentration for $C_1$ looks like

At $t=1e-2$ the concentration for $C_1$ looks like

But since the value of the source term is negative, I suppose it should be a upward parabola at the interface.
From time $t=1e-3$ to $t=1e-2$, the distribution for $\phi_2$ changes dramatically at the left of the interface though.

Some concentrations on the other hand are smooth across the domain though, but $C_1$ and $C_8$, which have the two largest initial condition values, are discontinuous like the picture shown at the interface.

After $t=1e-2$ the RADAU5 solver converges very slowly (I guess it won't converge to a meaningful values anyway).
I'm using cell-centered finite volume for it is easier to implement the Neumann boundary conditions. The discontinuous point is at the first cell right at the left of the interface. I'm using $N=100$; for $N=400$ the discontinuous point even show earlier.

What could be the problems? And how should I resolve it?

Thanks!

Note:

It's assumed that concentration and potential are continuous at the interface, which in fact they are, but in a cell-centered grid, I'm thinking I can't just naively do nothing at the interface.

Title: Discontinuity at Interface


问:
Oscillating convergence plot in my Resilient BackPropagation (RPROP) implementation
I have implemented in matlab a neural network that uses rprop's algorithm to update its weights.

Strangely the error on the training set does not converge to a local minimum, but oscillates.

Here is the graphics plot of the error function on trainingset:

Here is the rprop's algorithm:

function [net] = resilientBackPropagation(net, DW, DB, ETA_PLUS, ETA_MINUS)

for i=1:length(net.W)
    productDW = net.DW{i}.*DW{i};
    productDB = net.DB{i}.*DB{i};

```<code>indDW_gt_0 = find(productDW > 0);
indDB_gt_0 = find(productDB > 0);
indDW_lt_0 = find(productDW < 0);
indDB_lt_0 = find(productDB < 0);
indDW_eq_0 = find(productDW == 0);
indDB_eq_0 = find(productDB == 0);

net.deltaMarginW{i}(indDW_gt_0) = min(ETA_PLUS.*net.deltaMarginW{i}(indDW_gt_0), 50);
net.deltaMarginB{i}(indDB_gt_0) = min(ETA_PLUS.*net.deltaMarginB{i}(indDB_gt_0), 50);

   net.deltaMarginW{i}(indDW_lt_0) = max(ETA_MINUS.*net.deltaMarginW{i}(indDW_lt_0), exp(-6));
net.deltaMarginB{i}(indDB_lt_0) = max(ETA_MINUS.*net.deltaMarginB{i}(indDB_lt_0), exp(-6));

    net.W{i}(indDW_gt_0) = net.W{i}(indDW_gt_0)-(sign(DW{i}(indDW_gt_0).*net.deltaMarginW{i}(indDW_gt_0)));
net.B{i}(indDB_gt_0) = net.B{i}(indDB_gt_0)-(sign(DB{i}(indDB_gt_0).*net.deltaMarginB{i}(indDB_gt_0)));

  DW{i}(indDW_lt_0) = 0;
DB{i}(indDB_lt_0) = 0;

     net.W{i}(indDW_eq_0) = net.W{i}(indDW_eq_0)-(sign(DW{i}(indDW_eq_0).*net.deltaMarginW{i}(indDW_eq_0)));
net.B{i}(indDB_eq_0) = net.B{i}(indDB_eq_0)-(sign(DB{i}(indDB_eq_0).*net.deltaMarginB{i}(indDB_eq_0)));

Totale nella Rete Neurale.
net.DW{i} = DW{i};
net.DB{i} = DB{i};
</code>```

end

What is the problem?
thanks


问:
Oscillating convergence in my Resilient BackPropagation (RPROP) implementation
I have implemented in matlab a neural network that uses rprop's algorithm to update its weights.

Strangely the error on the training set does not converge to a local minimum, but oscillates.

Here is the graphics plot of the error function on trainingset:

Here is the rprop's algorithm:

```<code>function [net] = resilientBackPropagation(net, DW, DB, ETA_PLUS, ETA_MINUS)

for i=1:length(net.W) % Ciclo sul Numero di Strati della Rete

    % Calcolo il Prodotto tra le Derivate dei pesi e dei bias di E e E^-1 dell'i-esimo Strato
    productDW = net.DW{i}.*DW{i};
    productDB = net.DB{i}.*DB{i};

    % Prelevo gli Indici dei pesi e dei bias dell'i-esimo Strato con Prodotto Positivo, Negativo e Nullo
    indDW_gt_0 = find(productDW > 0);
    indDB_gt_0 = find(productDB > 0);
    indDW_lt_0 = find(productDW < 0);
    indDB_lt_0 = find(productDB < 0);
    indDW_eq_0 = find(productDW == 0);
    indDB_eq_0 = find(productDB == 0);

    % Calcolo i Margini dei pesi e dei bias con Prodotto Positivo
    net.deltaMarginW{i}(indDW_gt_0) = min(ETA_PLUS.*net.deltaMarginW{i}(indDW_gt_0), 50);
    net.deltaMarginB{i}(indDB_gt_0) = min(ETA_PLUS.*net.deltaMarginB{i}(indDB_gt_0), 50);

    % Calcolo i Margini dei pesi e dei bias con Prodotto Negativo
    net.deltaMarginW{i}(indDW_lt_0) = max(ETA_MINUS.*net.deltaMarginW{i}(indDW_lt_0), exp(-6));
    net.deltaMarginB{i}(indDB_lt_0) = max(ETA_MINUS.*net.deltaMarginB{i}(indDB_lt_0), exp(-6));

    % Aggiornamento dei pesi e dei bias della Rete con Prodotto Positivo
    net.W{i}(indDW_gt_0) = net.W{i}(indDW_gt_0)-(sign(DW{i}(indDW_gt_0).*net.deltaMarginW{i}(indDW_gt_0)));
    net.B{i}(indDB_gt_0) = net.B{i}(indDB_gt_0)-(sign(DB{i}(indDB_gt_0).*net.deltaMarginB{i}(indDB_gt_0)));

    % Aggiornamento dei pesi e dei bias della Rete con Prodotto Negativo
    DW{i}(indDW_lt_0) = 0;
    DB{i}(indDB_lt_0) = 0;

     % Aggiornamento dei pesi e dei bias della Rete con Prodotto Nullo
    net.W{i}(indDW_eq_0) = net.W{i}(indDW_eq_0)-(sign(DW{i}(indDW_eq_0).*net.deltaMarginW{i}(indDW_eq_0)));
    net.B{i}(indDB_eq_0) = net.B{i}(indDB_eq_0)-(sign(DB{i}(indDB_eq_0).*net.deltaMarginB{i}(indDB_eq_0)));

    % Memorizzazione delle Derivate dei pesi e dei bias Aggiornati dell'Errore Totale nella Rete Neurale.
    net.DW{i} = DW{i};
    net.DB{i} = DB{i};

end
</code>```

What is the problem?

ps:the dataset used is banana_dataset and in the plot the x-axis represents the epoch of training and the y-axis the error.

ETA_MINUS = 0.5;
ETA_PLUS = 1.2


问:
Adaptive Finite Element Method - Laplace
I'm currently attempting to turn my code for solving the laplace equation using finite element approximations into an adaptive one using the dual weighted residual as my error estimator: i.e. my indicator for element i: $η_i=(f[i]−∑_0^N A_{ij}U_h[j])∗z[i]$. i.e. I set up f, A and z over a very fine mesh and solve U_h over the normal mesh then interpolate it so i can perform the above operation.

I am solving for the dual solution z relatively simply; i.e. B(v,z) = Q(v) where B is the stiffness matrix and Q is a vector of 0s with 1 at the quantity of interest.

However the errors do not behave how they should; as an example i got the mesh to just increase from 10 in uniform multiples of 10 until 150. However these are the error estimates (as solved above):

```<code>0.0335462065038
0.0335462065038
0.0335462065038
0.0335462065038
0.0335462065038
0.0333608591916
0.0335462065038
0.0334614504804
0.0333158830798
0.0335462065038
0.0334317537862
0.0333344104201
0.0335441544821
0.0334327061457
</code>```

They seem to drop then grow, drop then grow repeatedly and i have no idea why, does anybody have any clues what could be causing this?

The following is my code

```<code>def rhs(x0):
    rhsfunction = np.empty(len(x0))
    for i in range(len(x0)):
        rhsfunction[i] = math.sin(math.pi*x0[i])
    return rhsfunction

def Poisson_Stiffness(x0):
    """Finds the Poisson equation stiffness matrix with any non uniform mesh x0"""

    x0 = np.array(x0)
    N = len(x0) - 1 # The amount of elements; x0, x1, ..., xN

    h = x0[1:] - x0[:-1]

    a = np.zeros(N+1)
    a[0] = 1 #BOUNDARY CONDITIONS
    a[1:-1] = 1/h[1:] + 1/h[:-1]
    a[-1] = 1/h[-1]
    a[N] = 1 #BOUNDARY CONDITIONS

    b = -1/h
    b[0] = 0 #BOUNDARY CONDITIONS

    c = -1/h
    c[N-1] = 0 #BOUNDARY CONDITIONS: DIRICHLET

    data = [a.tolist(), b.tolist(), c.tolist()]
    Positions = [0, 1, -1]
    Stiffness_Matrix = diags(data, Positions, (N+1,N+1))

return Stiffness_Matrix

def NodalQuadrature(x0):
    """Finds the Nodal Quadrature Approximation of sin(pi x)"""

    x0 = np.array(x0)
    h = x0[1:] - x0[:-1]
    N = len(x0) - 1

    approx = np.zeros(len(x0))
    approx[0] = 0 #BOUNDARY CONDITIONS

    for i in range(1,N):
        approx[i] = math.sin(math.pi*x0[i])
        approx[i] = (approx[i]*h[i-1] + approx[i]*h[i])/2

    approx[N] = 0 #BOUNDARY CONDITIONS

return approx

def Solver(x0):

    Stiff_Matrix = Poisson_Stiffness(x0)

    NodalApproximation = NodalQuadrature(x0)
    NodalApproximation[0] = 0

    U = scipy.sparse.linalg.spsolve(Stiff_Matrix, NodalApproximation)

return U

def Dualsolution(rich_mesh,qoi): #BOUNDARY CONDITIONS?
    """Find Z from stiffness matrix Z = K^-1 Q over richer mesh"""

    K = Poisson_Stiffness(rich_mesh)
    Q = np.zeros(len(rich_mesh))
    if qoi in rich_mesh:

        qoi_rich_node = rich_mesh.tolist().index(qoi)
        Q[qoi_rich_node] = 1.0
    else:

        nearest = find_nearest(rich_mesh,qoi)
        if nearest < qoi:

            qoi_lower_node = rich_mesh.tolist().index(nearest)
            qoi_upper_node = qoi_lower_node+1

        else:

            qoi_upper_node = rich_mesh.tolist().index(nearest)
            qoi_lower_node = qoi_upper_node-1

        Q[qoi_upper_node] = (qoi - rich_mesh[qoi_lower_node])/(rich_mesh[qoi_upper_node]-rich_mesh[qoi_lower_node])
        Q[qoi_lower_node] = (rich_mesh[qoi_upper_node] - qoi)/(rich_mesh[qoi_upper_node]-rich_mesh[qoi_lower_node])
    Z = scipy.sparse.linalg.spsolve(K,Q)
    return Z

def Indicators(richx0,basex0,qoi):
"""interpolate U, eta_i = r_i z_i = (f_i - sum(AijUj)Zi"""
    U = Solver(basex0)
    U_inter = interp1d(basex0,U)
    U_rich = U_inter(richx0)
    f = rhs(richx0)
    A = Poisson_Stiffness(richx0).tocsr()
    Z = Dualsolution(richx0,qoi)

    eta = np.empty(len(richx0))
    for i in range(len(eta)):
        eta[i] = f[i]
        for j in range(len(richx0)):
            eta[i] = eta[i] - A[i,j]*U_rich[j]
        eta[i] = eta[i] * Z[i]
</code>```

I understand theres alot of code here but i have no idea where else to turn. I'm almost certain my functions for solving the finite element method are correct (i.e. the Stiffness matrix, nodal quadrature and the solver). Any help would be appreciated alot.
Thanks, James

EDIT, this is my plot of the error indicator which must be incorrect but i cannot figure out what is wrong with it.


问:
Molecular dynamics: Lennard-Jones forces in 2D system
As I am new to all these stuff I would like to verify whether my procedure is correct.

I need to calculate the Lennard-Jones potential, potential energy of the system and the forces.

The procedure is as follows:

1) Read the input data: side lenght of the square -a, number of particles -N, the parameters $\sigma$, $\epsilon$

2) For $i=1,n$ insert randomly n points in x and y axis.

2a) For $j=i+1,n+1$ calculate the difference to square between
particles in x and y axis.

2b) If $r_{ij} \leq 2^{1/6}\sigma$ go back to 2a) step

2c) End of the both loops

3) For $i=1,n$

3a) For $j=i+1,+1$ calculate the lennard jones potential $$u_{ij}=4\epsilon[(\sigma/r_{ij})^{12}-(\sigma/r_{ij})^6]$$ where $r_{ij}$ is being calculated in 2a) step

3b) calculate the potential energy of the system
$$U_N=\sum_{i=1}^N\sum_{j=i+1}^N u_{ij}(r_{ij})$$

3c) Calculate the force $$F=\sum_{j=1,j\neq i}^NF_{ij}$$ where $$F_{ij}=(-\nabla_ju_{ij}(r_{ij}))$$

3d) End of the both loops

I would really grateful for checking, whether it is even correct.
I would also appreciate any books connected to this problem - I am currently using Allend Tildesley, but it is rather difficult for me at start.


问:
Difference between fast and normal Givens Rotations?
would someone be so kind as to explain me the difference between the ordinary givens-rotation and the fast givens-rotation?

I know that the fast givens Rotation reduces the Count of operations to approximately the half. How is this possible?

Thanks in advance for any hints :-)


问:
How to implement conjugate gradient method to minimize this nonlinear action?
Given a 2D stochastic differential equation:
\begin{align}
\dot{x}_{i}=f_{i}(\textbf{x})+g_{ij}\xi_{j}(t),
\end{align}
where $i=2$, $g_{ij}g_{jk}=2\epsilon\delta_{ik}$ and $f(\textbf{x})=-\nabla\phi(\textbf{x})+l(\textbf{x})$ with $\nabla\phi(\textbf{x})\cdot l(\textbf{x})=0$, we have its action functional for trajectories:
\begin{align}
S_{T}(x(t))&=\int_{0}^{T}dt\frac{1}{4\epsilon}(\dot{x}_{i}-f_{i})^{2}.
\end{align}

Let us now consider a 2D potential function with double well in $x$-direction: $\phi=\frac{1}{2}(y-1)^{2}-\frac{1}{2}(x-1)^{2}+\frac{1}{4}(x-1)^{4}$, and a given $l(\textbf{x})$.
I want to minimize the corresponding discretized action for the set of trajectories connecting the two lowest energy points:
\begin{align}
S_{T}(x(t))=\Delta t\sum_{k=1}^{M-1}\sum_{i=1}^{2}\Big[\frac{x_{i}^{k+1}-x_{i}^{k}}{\Delta t}-\frac{1}{2}(f_{i}^{k+1}+f_{i}^{k})\Big]^{2},
\end{align}
where the trajectory $x(t)$ is divided into $M-1$ time segments.

How to implement the conjugate gradient method to minimize this nonlinear action?


问:
should boundary conditions be effecting moving mesh results?
I have a question on the use of moving mesh to solve the inviscid euler equations. I have solved the following equations:
$$\frac{\partial}{\partial t}\left[\begin{array}{c}
\rho\\
\rho u
\end{array}\right]=\frac{\partial}{\partial x}\left[\begin{array}{c}
\rho(u-w)\\
\rho u(u-w)+p
\end{array}\right],$$
where p is a linear function of $p=f(\rho)$ and $w$ is the mesh velocity.

I have used the integral form of FVM, treated the right hand side in its conservative form and have evaluated the results at the cell centres. At the left BC I have wall conditions, so $u_L = w_L = 0$ and at the right hand BC, $u_R = w_R = f(t)$. I have used linear interpolation for the other boundary terms. I am moving the mesh using $dx/dt = w$.

I know that moving the mesh should not effect the solution, however, my results seem to be effected since for instance for the $\rho u (u-w)$ term of the momentum equation, on the right hand side, when $w=0$, I have $\rho_R u_R u_R$ but when $w\neq 0$ I will have $\rho_R u_R (0-0) = 0$.
Can someone please let me know if I am doing the right thing here, and if it is possible to avoid this issue somehow? Thank you.


问:
Variational Monte Carlo to calculate local energy of hydrogen like ions in python
I'm writing up a code to calculate the local energy for electrons in hydrogen like ions for a given wavefunction. My code is giving me weird results, which leads me to believe something is wrong.

Firstly, am I going about accepting/rejecting moves the right way? And am I doing anything weird in the accumulation phase which gives me the wrong answer?

Basically, I know from mathematica that my local energy for the trial wavefunction I'm using should be -0.48, however I keep getting a whole splash of results, with little consistency (eg 0.6, -0.38, 0.7 etc.). Am I doing anything weird?

EDIT: sorry about the hard to read code, I'm new to this coding thing and am not very sure I'm doing the right thing. Essentially, what I am trying to do is sample my local energy with the probability distribution of my trial wavefunction. I'm not sure if I'm accumulating the sampled energy correctly.

for reference here is my code

```<code>import math
import random

def psi(r): #my trial wavefunction
    return (math.exp(-1.2*r))

#hamiltonian
def Ham(r): #the hamiltonian operator
    return (1.44*psi(r)+(2/r)*(-1.2*psi(r)))
def Top(r): #numerator of the El expression
return (-0.5*psi(r)*Ham(r)*r*r)
def But(r): #denomenator of the El expression
    return (psi(r)*psi(r)*r*r)
def Elocal(r): #defines the local energy
    return (Top(r)/But(r))

def r(x,y,z):
    return math.sqrt(math.pow(x,2)+math.pow(y,2)+math.pow(z,2))

x0 = random.uniform(0, 50)
y0 = random.uniform(0, 50)
z0 = random.uniform(0, 50)

#local energy
El = 0

#my movement weight
dr = 0.5

#begin equilibriate
n = 0
for n in range (1000):

    dx = random.uniform(-1,1)*dr
    dy = random.uniform(-1,1)*dr
    dz = random.uniform(-1,1)*dr

    x = x0+dx
    y = y0+dy
    z = z0+dz

    dpsi = psi(r(x,y,z))
    opsi = psi(r(x0,y0,z0))

    W = (dpsi*dpsi)/(opsi*opsi)

    if W >= random.random():
        x0 = x
        y0 = y
        z0 = z
    else:
        x0 = x0
        y0 = y0
        z0 = z0
    n = n+1
# begin accumulate
w = 0
N = 0
for N in range (1000):

    dx = random.uniform(-1,1)*dr
    dy = random.uniform(-1,1)*dr
    dz = random.uniform(-1,1)*dr

    x = x0+dx
    y = y0+dy
    z = z0+dz

    dpsi = psi(r(x,y,z))
    opsi = psi(r(x0,y0,z0))

    W = (dpsi*dpsi)/(opsi*opsi)

    if W >= random.random():
        x0 = x
        y0 = y
        z0 = z
        El = El+Elocal(r(x,y,z))
        w = w+1
    else:
        x0 = x0
        y0 = y0
        z0 = z0
        El = El
    N = n+1
print El/w
</code>```


问:
Area of convex n-dimensional polytope
I am looking for an efficient algorithm to calculate the surface area of an irregular N-dimensional polytope.

I have a description of this polytope both as coordinates of the vertices as as linear inequalities.

Title: Area of convex n-dimensional polytope


问:
what do positive real parts of eigenvalues mean?
I am solving a 1D advection problem of the the form $$d{Q}/dt = [A]{Q}$$ where {Q} is the vector of unknowns and [A] is the matrix of coefficients of spatial discretisation. I have worked out the eigenvalues of [A] to get an idea about the stability of the semi-discretised system and I am getting eigenvalues with positive real parts plus imaginary parts. I know positive real parts are a sign of instability, but my results show the smooth travelling of a wave. I wanted to know what conclusions I can draw from this?
Thanks! </p>


问:
localizing rings in SAGE
I am a bit of SAGE newb, so I apologize if this is a simplistic question.

I want to use SAGE to compute some invariants of elliptic curves defined over rings which aren't the integers $\mathbb{Z}$ or the rationals $\mathbb{Q}$. More specifically, I am interested looking at elliptic curves defined over the ring
$$
\mathbb{Z}[1/6][c_4,c_6]
$$
or
$$
\mathbb{Z}_{(p)}[c_4,c_6]
$$
where $p$ is some prime. Can anyone tell me how produce a localized ring such as $\mathbb{Z}[1/6]$ or $\mathbb{Z}_{(p)}$ in SAGE? Or know where in the manual I could find such information?

Thanks!

Title: localizing rings in SAGE


问:
Calculating tangential electric field intensity on the boundary
I have an object with surface charge $\sigma$ for which I want to caclulate tangential electric field which would correspond to mathematical formulation:
$$
\vec{E}(\vec{x}) = \nabla_x \int \frac{\sigma_y}{|\vec{x} - \vec{y}|} dS_y
$$
where afterwards normal component $\vec{n} (\vec{E} \cdot \vec{n})$ is substracted.

For now I have tried to put this surface charge on surface mesh nodes and then transforming the integral to summations (skipping singularity) assuming that integral value on triangle is average of it on nodes. In this way my calculation totally fails.

As looking why it does not work I started to suspect that higly singular kernel in the integral. To overcome it I look on some different ways:

Interpolating the charge density piecwise linearly and then doing the integration

Push the charges inside the suface about a size of mesh triangle

What other options are available and which do you recommend? Also would it help if I calculate potential either of theese methods and afterwards taking derivative of it?


问:
comparison of stability of two non-linear methods
I have solved a numerical problem using two different sets of non-linear governing equations. I want to get an understanding of the stability of the methods relative to each other. To do so, I solving the problem with each method whilst keeping the time stepping schemes, time step and mesh sizes the same. I then run the code, and observe the largest time step using which stable results are produced or convergence is obtained. I wanted to know whether this makes sense. I know this process does not get definite answer about stability of the methods.

Also what I have found that for this particular problem, stable results are still obtained even with CFL numbers greater than 1 (very large time steps). Do the large time steps I try still need to satisfy CFL?


问:
1-D finite differences with piecewise linear solution
I was having some lectures and I didn't quite understand the following: let's say you have a grid like
$$G=\{ x \in \mathbb{R} : x = x_j = hj, \  j = 0,1,...,n,\ h=1/n\}.$$
And you write in difference form the following 1D diffusion  equation (comma before index indicates derivative along the index):
$$-(au_{,1})_{,1}=s, \qquad x \in \Omega := (0,1).$$

Giving:
$$ \{ -(a_{j-1} + a_j)u_{j-1} + (a_{j-1}+2a_j+a_{j+1})u_j-(a_j+a_{j+1}  )u_{j+1}\}/2h^2 = s_j$$
Now let's say you have a discontinuity of $a$ along some interface $\Gamma$. In our 1D problem, let:
$$a(x) = \epsilon, \ \ 0<x\le x^*, \ \ a(x) = 1, \ \ x^*<x<1.$$

For the boundary conditions $u(0) = 0$ and $u(1) = 1$. This jump can be written:
$$\epsilon \,\underset{x \uparrow x^*}{\lim u_{,1}} = \underset{x \downarrow x^*}{\lim u_{,1}}. $$

Here is what I don't understand: the lecturer claims that by proposing a linear piecewise solution for the diffusion equation it is found that:
$$u= cx, \ \  0\le x< x^*, \ \ u=\epsilon c x + 1 -\epsilon c, \ \ x^*\le x\le1.$$
Where:
$$c= 1/(x^*-\epsilon x^* + \epsilon).$$

This is done with the difference equation, assuming $x_j < x^* \le x_{j+1}$, giving:
$$u_j= \alpha j, \ \  0\le j< k, \ \ u=\beta j + 1 -\beta n, \ \ k+1\le j\le n.$$
Where,
$$ \beta =\epsilon \alpha, \ \ \alpha=\left( \epsilon\frac{1-\epsilon}{1+\epsilon} + \epsilon(n-k)+k\right)^{-1}.$$
Hence:
$$u_k = x_k/(\epsilon h(1-\epsilon)/(1+\epsilon)+(1-\epsilon)x_k+\epsilon).$$

How are those Piecewise Linear Solution built?

Any input will be much appreciated.


问:
How to classify chaotic systems from a stability perspective
I am wondering what chaotic systems are from the perspective of numerical analysis. I am talking about 'deterministic chaos' such as for instance the 'logistic map' exhibits it. That is, the solution after many steps in the iterative solution is highly dependent on the initial conditions. It is clear that it will also be highly sensitive to numerical errors. But are chaotic systems

well posed (Hadamard),

well conditioned,

always unstable in implementations,

susceptible to regularization?


问:
Time dependent self-consistent equations
I am facing the following problem. I need to solve numerically a set of coupled equations </p>

<p>$$i\frac{d}{dt}f_{n}^{(i)}(t) = \left[U\cdot n(n-1) + \mu\cdot n\right]f_{n}^{(i)}(t) - \sqrt{n+1}\Phi_i^{*}\ f_{n+1}^{(i)}-\sqrt{n}\Phi_{i}\ f_{n-1}^{(i)}$$
where $U,\mu$ are just constants and
$$\Phi_{i} = \sum\limits_{n=1}^{N}\left( f_{n-1}^{(i+1)}(t)\right)^{*}\ f_{n}^{(i+1)}(t)\sqrt{n}$$
has to be determined self-consistently. These are coupled differential equations where $i = 1,2,\ldots, M$ and $n = 0,1,,\ldots, N$. What kind of numerical methods can be used in order to solve this problem efficiently? To be clear, I am looking for time evolution of each $f_{n}^{(i)}$ (complex numbers).</p>


问:
Is there a term for Goodhart's Law in the context of optimization?
Let's say I'm optimizing something. To pick an arbitrary example, let's say I'm choosing the shape of some part to maximize strength-to-weight ratio. So I get some FEM software, parametrize the shape, and run gradient ascent or whatever to find the optimum shape.

When I do this procedure, I would generally expect that the calculated optimal strength-to-weight ratio is an overestimate of the actual strength-to-weight ratio for that shape---especially if I'm running the FEM software at low accuracy. Basically, the optimization algorithm will seek out and exploit any inaccuracies in the software to create impossibly high performance metrics.

In the context of economics, this basic insight is famously known as Goodhart's Law. In the context of optimization, how do people refer to this? Is there a term for it? Is there a paper or textbook that I can cite that discusses it?


问:
Maxwellian distribution of velocities with Shake algorithm present
I am writing a code to perform hybrid monte carlo molecular dynamics. To do this, I need to have a code to initialize the velocities of all particles according to a maxwell distribution. The code is as follows

```<code>   subroutine thermalize
        implicit none
        double precision :: MassInv
        double precision :: RATIO
        integer :: i

        !--- temp is temperature set point value
        STDtemp = sqrt(temp)
        do i = 1 ,np
           MassInv = 1.0d0/(mass(position(i)%type))
           MassInv = sqrt(MassInv*kboltz*kcal2amu)

  !---Gaussian() is a functions that returns gaussian random variable n(0,1)
           v(i)%x = Gaussian()*STDtemp*MassInv*ratio
           v(i)%y = Gaussian()*STDtemp*MassInv*ratio
           v(i)%z = Gaussian()*STDtemp*MassInv*ratio

        enddo
      end subroutine thermalize
</code>```

To measure the temperature

```<code>T = 2.0*KE/(kboltz*(3*np - NCONS)
</code>```

where KE is the kinetic energy, np is the number of atoms, kbotlz is boltzmanns constant. I was trying to run some tests with temp=380, and found that my thermalize function was giving me initial kinetic energies that translated to temperatures around T=550 K. When I turned shake off, and hence ncons = 0, the temperature produced by thermalize was 379. So I think I need to modify the thermalize command to account for the constrained DOF's

To try to account for the shake, the ratio of the KE with shake and without shake present for a fixed temperature is

```<code>1/2kT*(3np-ncons)   = 1/2kT*(3np)
</code>```

dividing and rearranging

```<code>KE_shake/KE_noshake = (3np-ncons)/np
</code>```

Thus, I am thinking about multiplying by the sqrt of that ratio.

```<code>v(i)%x = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np)
v(i)%y = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np)
v(i)%z = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np)
</code>```

This gives me the correct temperature (~380). However, when I thermalize my particles, and run an NVE trajectory of a previously equilibrated fluid (temperature should fluctuate around 380), the temperature drops to about 300. I have no idea if this is statistically mechanically correct.

Any thoughts?


问:
ARPACK- Reverse Communication Interface
I am trying to use ARPACK to diagonalize a sparse Hermitian matrix.

My program is executing without error but returning incorrect eigenvalues.  I have tried to trace the problem back.

The work is done by two functions:

znaupd and a function I have defined av to perform matrix vector multiplication.

znaupd takes a vector as input and passes a resultant vector to av.  There is some starting default vector that ARPACK chooses.

Looking at the output of my program, it seems like the default vector is essentially being chosen to be the zero vector, and znaupd just passes this same zero vector to av every time.

Title: ARPACK- Reverse Communication Interface


问:
Find constrained vectors maximizing angles between them - methods?
This is related to a question I had asked earlier, with the distinction that earlier I did not have a non-linear objective functional to minimize. The problem is reproduced below with added information. I used MATLAB to solve the problem posed in that question.

Why I need to minimize the objective functional given below: As it turns out, the constraints I had put on the vectors were not enough to guarantee that the vectors be "as different" from each other as possible. Some of the vectors were completely identical. Therefore, I though it might help to maximize the angles between the vectors, kind of driving them towards orthogonality (right?).

Problem: You are given a set of 16 vectors, $\{\mathbf{a}^i\}_{i=1}^{16}$, and matrices $D_1$ and $D_2$ (both of dimension $[10\times6]$. Let the following be true:

$\sum_{i=1}^{16} \mathbf{a}^i = \mathbf{1} \in \mathbb{R}^6\;,$

$\mathbf{a}^i \geq \mathbf{0}\;,\forall i\;,$ (i.e. all elements positive),

$(D_1 + D_2)\mathbf{1} = \mathbf{1} \in \mathbb{R}^{10}\;.$

Then, find vectors $\{\mathbf{b}^i\}_{i=1}^{16}$ and $\{\mathbf{c}^i\}_{i=1}^{16}$ that minimize the functional:
$$
\min \sum_{i,j=1}^{16} \biggl\vert \dfrac{\pi}{2} - \cos^{-1}\biggl(\dfrac{\mathbf{d}^i\cdot \mathbf{d}^j}{||\mathbf{d}^i||||\mathbf{d}^j||}\biggr)\biggr\vert\;,
$$
where $\mathbf{d}^i := \left(\begin{array}{c}
\mathbf{b}^i\\
\mathbf{c}^i
\end{array}\right)$, and such that they satisfy the following constraints:

$\mathbf{c}^i - D_1 \mathbf{b}^i = D_2 \mathbf{a}^i\;,\quad \forall i\in \{1,2,\dots,16\}\;,$

$\sum_{i=1}^{24} \mathbf{c}^i = \mathbf{1} \in \mathbb{R}^{10}\;,$

$\sum_{i=1}^{24} \mathbf{b}^i = \mathbf{1} \in \mathbb{R}^6\;,$

$\mathbf{c}^i \geq \mathbf{0}\;,\forall i\;,$ (i.e. all elements of each vector positive),

$\mathbf{b}^i \geq \mathbf{0}\;,\forall i\;,$ (i.e. all elements of each vector positive).

Questions

Could anyone please give me any pointers on how I could go about tackling this problem? I have no idea about any methods I could use to solve this problem computationally.

Any comments on the objective functional I have defined? I was also thinking of $
\min \sum_{i,j=1}^{16} \biggl(\dfrac{\pi}{2} - \cos^{-1}\biggl(\dfrac{\mathbf{d}^i\cdot \mathbf{d}^j}{||\mathbf{d}^i||||\mathbf{d}^j||}\biggr)\biggr)^2\;,
$
but thought that they would be equivalent. Any thoughts? Could I do better than these two?

Is there a way to devise a functional to drive them towards linear independence?

Also, this is NOT homework.

Thanks a lot.


问:
What are the career options for computational scientists other than being in software Industry?
I would like to know the career options for a computational scientist other than getting into software industry. I wish to use the my skills of Numerical Linear Algebra,Parallel Programming, Scientific Visualization,Topological Data Analysis in some Research Lab of any manufacturing company or Scientific company.</p>


问:
Effect of Periodic Boundary condition in Pair Distribution function
When we simulate a bulk molecular system using molecular dynamics in (say) NVE ensemble (microcanonical), we apply some periodic boundary conditions (PBC) to avoid surface effects. What effect does this have on quantities like Pair Distribution function, since it so happens that part of a molecule is near one edge of the box and the other is in a different edge of the box. What effect does this have on PDF, MSD and other such thermodynamic quantities ?

PS : I am currently simulating SPC/E water molecules using cubic box, in some frames I find water molecules which are split apart, but connected by PBC.


问:
FEniCS: both normal and shear stress boundary conditions for elasticity?
I would like to have both the normal (xx) component and shear (xy) component of a 2D (stress) tensor defined on a boundary (y=const, for instance) for an equation which is of the type

$$
\nabla \cdot \sigma = 0
$$

The normal value of $\sigma_{xx}=P$ at $x=0,1$ can be accounted for via a Neumann boundary condition, in the term obtained by integration by parts in weak form:

$$
\int_\Omega  \sigma : \nabla \tau
=  \int_{\partial\Omega}  P \cdot \tau \, d s .
$$

$\tau$ is a vector test function for the trial function $\vec{u}$, which physically is the displacement.
However, I am having difficulty fixing the shear stress at the boundaries $x=0,1$ of the unit square, as a simple example of my goal. How can one best accomplish this in FEniCS?

My question is similar to this one, but as far as I can tell, they are distinct.

A minimal running code in FEniCS (in python) which gives incorrect results for $\sigma$ is below. Thank you!

```<code>from dolfin import *

# Mesh Params ################
E, nu = 1.0, 0.4 #material params
mesh = UnitSquareMesh(100, 100)

print 'Creating Lagrangian function spaces...'
Vv = VectorFunctionSpace(mesh, "Lagrange", 1)
Vf = FunctionSpace(mesh, "Lagrange", 1)
Vt = TensorFunctionSpace(mesh, "Lagrange", 1)

# P is stress vector of stress values normal and tangent to boundary
P = Expression(('(x[0]-xc)','(x[0]-xc)'),xc=0.5)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-7
        return on_boundary and abs(x[0]) < tol

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-7
        return on_boundary and abs(x[0] - 1) < tol

left_boundary = LeftBoundary()
right_boundary = RightBoundary()

boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_parts.set_all(0)
right_boundary.mark(boundary_parts, 1)
left_boundary.mark(boundary_parts, 2)
ds = Measure("ds")[boundary_parts]

##############################
# Define variational problem #
##############################
print 'Defining test and trial functions...'
v = TestFunction(Vf)
tau = TestFunction(Vv)
phi = TrialFunction(Vf)
u = TrialFunction(Vv)

# functions for variational problem
d = u.geometric_dimension()
epsilon = sym(nabla_grad(u))
# PLANE STRESS
sigma = E/(1+nu)* (epsilon) + E*nu/((1-nu**2)) * tr(epsilon)*Identity(d)

# Define normal and tangent directions on mesh
norm = FacetNormal(mesh)
tang = as_vector([norm[1], -norm[0]])

tau_n = dot(tau,norm)
tau_t = dot(tau,tang)
P_n = dot(P,norm)
P_t = dot(P,tang)

#################
# SOLVE
#################
au = inner(sigma, sym(nabla_grad(tau)))*dx
Lu = dot(P_n,tau_n)*ds(1) + dot(P_n,tau_n)*ds(2) + dot(P_t,tau_t)*ds(1) + dot(P_t,tau_t)*ds(2)
u = Function(Vv)
solve(au == Lu, u)

## View the results
sigma = E/(1+nu)* (sym(nabla_grad(u))) + E*nu/((1-nu**2)) * tr(sym(nabla_grad(u)))*Identity(d)
s_V = project( sigma, Vt)
sxx = project(s_V[0,0], Vf);
sxy = project(s_V[0,1], Vf);
syx = project(s_V[1,0], Vf);
syy = project(s_V[1,1], Vf);
plot(sxx, interactive=True, rescale=False,scalarbar=True)
plot(syy, interactive=True, rescale=False,scalarbar=True)
plot(sxy, interactive=True, rescale=False,scalarbar=True)
</code>```


问:
What strategies / decompositions would be useful to solve the following linear system repeatedly if I only care about time to solution?
Let $A\in \mathbb{R}^{n\times n}$ be symmetric positive semidefinite, and $B\in \mathbb{R}^{n\times n}$ be symmetric positive definite.  Suppose $B$ is block diagonal so it is easy to invert.  (We have access to the entries of $B^{-1}$ in the computer).  Let $\omega \in \mathbb{R}$.  Suppose I want to solve the following symmetric and positive definite system
$$
((A-\omega^2 B)B^{-1}(A-\omega^2 B)+\omega^2 B)x = b,
$$

for a range of values of $\omega$.  Note that the RHS, $b$, may depend on $\omega$.

I am interested in strategies to compute $x(\omega)$ repeatedly.  As an idealization, suppose I only care about correctness and the time it takes to produce a solution from a given $\omega$.  I mean, any time spent performing decompositions and preprocessing will be considered free.

What decompositions / strategies would be useful?

Edit

I have been playing around with this today.  This is my proposed solution.

We use a combination of a Cholesky decomposition and an eigen decomposition.  Let $C\in \mathbb{R}^{n\times n}$ be a lower triangular matrix such that $CC^\mathrm{T} = B$.  Then

$$
\begin{align*}
((A&-\omega^2 B)B^{-1}(A-\omega^2 B)+\omega^2 B)\\
&= \cdots \text{just algebra}\\
&= C\left((C^{-1} A C^{-\mathrm{T}} - \omega^2 I)(C^{-1} A C^{-\mathrm{T}} - \omega^2 I)+ \omega^2 I) \right)C^{\mathrm{T}}
\end{align*}
$$

Now we take an eigen decomposition of $C^{-1} A C^{-\mathrm{T}}$, i.e.

$C^{-1}A C^{-\mathrm{T}} = Q\Lambda Q^\mathrm{T}$ for $Q\in \mathrm{R}^{n\times n}$, orthonormal, $\Lambda\in \mathbb{R}^{n\times n}$, diagonal.

Then
\begin{align*}
 C\big((C^{-1} &A C^{-\mathrm{T}} - \omega^2 I)(C^{-1} A C^{-\mathrm{T}} - \omega^2 I)+ \omega^2 I) \big)C^{\mathrm{T}} \\
&= \cdots \text{(just algebra)}\\
&= CQ \big((\Lambda - \omega^2 I)(\Lambda- \omega^2I)+ \omega^2 I\big) Q^\mathrm{T} C^\mathrm{T}
\end{align*}

Now, this system can be trivially inverted, and none of the decompositions depends on $\omega$.

The real question is if any of this can be used to design a preconditioner for a system like this that is too large to be factored in this way.  Any ideas on this would be much appreciated.


问:
Limitations of particle in cell method for high density plasma
Are there any limitations of particle in cell (PIC) method for high density plasma? To be more specific, is modeling of a narrow channel of high density plasma possible or are there any limitations connected with PIC approximation?</p>


问:
Accelerated convergence for Sparse NMF
In the Non-Negative Matrix factorization (NMF), you basically compute an approximation of a given matrix $V \in \mathbb{R}_{+}^{n \times m}$ into matrices $W$ and $H$ such that $W \in \mathbb{R}_{+}^{n \times r}$ and $H \in \mathbb{R}_{+}^{r \times m}$. Or in other words, $V \approx WH$. For the Sparse NMF with sparsity parameter $\lambda\ge 0 $ and sparsity function $g$, you minimize a cost function of this form:
$$F(W, H) = \frac{1}{2}\|V - \bar W H\|^2 + \lambda\sum_{ij}g\left(H_{ij}\right) \qquad\text{where } \bar W =\frac W {\|W\|}$$

And $\|\ldots\|$ can be any differentiable norm. For $g(H_{ij}) =H_{ij}$, the steps for updating $H$ and $W$ are suggested here as follows  for :

\begin{align}H_{ij} \longleftarrow H_{ij}&\odot\frac{V_i^T \bar W_j}{R_i^T\bar W_j + \lambda}\\
W_j \longleftarrow W_j&\odot\frac{\sum_i H_{ij}\left[V_i + \left(R_i^T\bar W_j\right)\bar W_j\right]}{\sum_i H_{ij}\left[R_i + \left(R_i^T\bar W_j\right)\bar W_j\right]}\end{align}

Where $R = WH$ and the matrix division and the $\odot$ multiplication  are done
element-wise.

In this paper, it is suggested to accelerate the convergence as follows:

To speed up the convergence, the multiplicative term can be exponentiated with an acceleration parameter $\delta$ larger than one. For each iteration the costfunction is evaluated and if the cost is smaller than in the previous iteration, $\delta$ is increased and if it is larger, then $\delta$ is decreased.

This is done in the $H$ update step. Following the notations of the updates above, this will be as follows:
$$H_{ij} \longleftarrow H_{ij}\odot\left(\frac{V_i^T \bar W_j}{R_i^T\bar W_j + \lambda}\right)^\delta$$

Question:

How can one guarantee that with this acceleration we'll have fewer computational steps until convergence than without. Or how to ensure a minimum number of increasing and decreasing  of $\delta$ until convergence ?

Title: Accelerated convergence for Sparse NMF


问:
Constructing 2 fold oversampled cosine basis in MATLAB
So I'm trying to construct a 2 fold oversampled cosine basis in MATLAB. I know how to construct the basis as a square matrix using the following command:

```<code>dct(eye(n,n))
</code>```

where dct is the discrete cosine transform function and $n\times n$ is the size of my basis but how would I do this if I wanted my basis to be oversampled? Specifically for it to be $n\times 2n$, i.e. a fat matrix such that the columns lose the property of linear independence?

Any help would be greatly appreciated. Thanks in advance!


问:
Cholesky factorization of a block matrix
Let a Matrix $\mathbf{Y=[A\;b;b^T\;}d]$, where $d$ is a scalar, and $\mathbf{b}$ is a vector.
How to find the cholesky factorization of $\mathbf{Y}$ if i know the cholesky factorization of $\mathbf{A}$.?</p>


问:
Boundary conditions in the Finite Element Method at only one side of computational domain
I want to solve a Sturm-Liouville problem in 1D, i.e.,
\begin{align}
[p(x)\ u'(x)]'+q(x)u(x) = f(x)
\end{align}
with boundary conditions
\begin{align}
u(0)=a \hspace{1cm} u'(0)=b
\end{align}
How do I implement these conditions?

What I did:

First I multiply both sides by a test function $v\in\lbrace c\in H^1\big\vert v(0)=0, v'(0)=0 \rbrace$ and integrate over the computational domain, i.e. from 0 to 1 in my case.
\begin{align}
-\int\limits_0^1\left(p(x)u'(x)\right)'v(x)dx
+\int\limits_0^1q(x)u(x)v(x)dx
 =
\int\limits_0^1f (x)v(x)dx
\end{align}
 The first term can be transformed by partial integration.
\begin{align}
-p(1)u'(1)v(1) + \underset{=0}{\underbrace{p(0)u'(0)v(0)}} +\int\limits_0^1p(x)u'(x)v'(x)dx
\end{align}
I want to use hermite cubic splines denoted by $\phi_k^l$ with the property
\begin{align}
\left(\phi_k^l\right)^{(m)}(x_j)=\delta_{kj}\delta_{lm}
\end{align}
where $(m)$ denotes the 0th or 1st derivative. Inserting this in the weak formulation above we finally arrive at the matrix equation
\begin{align}
\nonumber
-p(1)u^1_{N +1}&\delta_{0l}\delta_{N+1j}
+a\int\limits_0^1
dx p(x)\left[\phi_0^0\,'(x)\phi_j^l\,'(x) + q(x)\phi_0^0(x)\phi_j^l(x)\right]\\
& +b\int\limits_0^1
dx\, p(x)\left[\phi_0^1\,'(x)\phi_j^l\,'(x) + q(x)\phi_0^1(x)\phi_j^l(x)\right]\\
&+\sum\limits_{i=1}^{N+1}\sum\limits_{k=0}^1u_i^k
\int\limits_0^1
dx\, p(x)\left[\phi_i^k\,'(x)\phi_j^l\,'(x) + q(x)\phi_i^k(x)\phi_j^l(x)\right]\\
=&\int\limits_0^1f (x)\phi_j^l (x)dx \hspace{2cm}j=1,\dots,N+1\ \ l=0,1
\end{align}
Is this derivation correct up to here?

The matrix equation is set up by first mapping every element to the standard intervall $[-1,1]$ and then using Gauss Integration.

Afterwards the global matrix is assembled. Finally boundary conditions have to be implemented according to the above equation. This means we delete the lines and columns corresponding to the coeffiicients of $\phi_0^0$ and $\phi_0^1$. Subtracting the second and third term (the ones with a and b in front) from the left hand side of the above equation on the right hand side and subtraction $p(1)$ at the corresponding position of the matrix we arrive at the final matrix equation. This is visualized in the figure Implementation of boundary conditions.

But using the extremely simple case of $p(x)=q(x)=1$, $f(x)=0$ and using only one domain we can find an analytic solution. The comparison between the analytic and the numerical solution is shown in figure comparision of solutions.
Thus there must be a mistake in the above implementation process (other boundary conditions provided at both ends of the domain give correct solutions.)

For an increasing number of domains the solution does also not converge. Therefore something is obviously wrong but I do not have any idea.

Can you please help me with this?


问:
Solving the Jump function in discontinuous Galerkin method
I am trying to solve the elasticity equation using Crouzeix-Raviart element. Since this can not be solved using continuous FEM I am trying to solve this by discontinuous Galerkin method. In the discontinuous Galerkin formulation the left hand part of equation $a(U,v)$ term has two parts, the first part is the stiffness and the second part provides the stabilization to the equation.

The stabilization term has the jump function $\int[U][v] $, where $[U]=U^+-U^- $ and $[v]=v^+-v^- $. I am finding it difficult to implement this. Could someone please suggest some literature which explains how the jump term is treated in matrix form and is solved.


问:
best package choice for visualizing agent based model simulations in python
This is a visualization/programming design problem. I am working on building a geospatial agent based model, meaning that there are multiple agents interacting within some particular spatial container. The movement of the agents and their interactions are based upon some biased random walks and integro-differential equations. I was trying to figure out which package to use for visualizing this, given some technical constraints. I wrote the model prototype in python, but here is the problem.

I need a way to visualize the interactions between these agents over time.
I looked at matplotlib, but the animation features will run into trouble when I add more than a few agents. I believe on the matplotlib site itself the authors indicate that matplotlib animations are not ready for even medium sized scale.

I want to stay open source--and free. The other choices seem to be Mayavi--which uses VTK. I have used this a bit in the past, but was not sure how well Mayavi can handle visualizing a fair number of agents--say between 100-1000. I also was not sure how actively maintained Mayavi was.

Vispy is another choice. This is also a very new library, so I am not sure how stable it is. It relies on on an interface to WebGL. I was not sure if this was a good choice to use yet, since not so many people have experience with it--not so many blog posts on Vispy yet.

Finally there are also some other hard core math visualization libraries in python like the finite elements solvers: Fenics or others. But I don't have experience with any of these, so not sure if they might work better for this particular problem.

If there are any choices I am missing, please suggest them.

Additional Note:

I have looked at agent based modeling software like Netlogo and Repast, etc. I don't really want to use these because they seem to impose a lot of constraints on the visualization and modeling choices--and I have to learn yet another language :).

UPDATE:

This was a note from the Vispy google group about VTK and mayavi.

```<code>"Mayavi is based on VTK, an old and big library which uses legacy OpenGL
(fixed function pipeline). Mayavi itself is mostly written in C++, and
has trouble attracting contributors. . . . the lead developer has said
its dying a slow death.

Vispy uses modern opengl and is implemented all the way in Python. It
only relies on one of several backends to provide a window to render in.
This makes Vispy more powerful (we can make use of more advanced
shaders) and flexible. Mayavi and VTK have been around for quite a while
though, so it will take some time for Vispy to match the functionality. "
</code>```


问:
Euler's first order Error
I was estimating a falling object's position versus time by using a simple first order step function, where

```<code>for i=2:length(t) % We're using Euler, so we need an initial previous point to start t=t(i);

V=V+dv;

rho_i=    % Calculate our height and then pull rho from data

K=(m/(mg-.5*C_d*A*V^2*rho_i));

 dv=dt/K;

  pos=pos-(V-V_previous)/2*dt;

 V_previous=V;

 velocity_Vector(i)=V;
 pos_Vector(i)=pos;
 end
</code>```

Now I'm wondering, how much will the error be after this thing falls, like many many kilometers?

I know that it's a first order approximation, so there's error...But say I needed to know the bounds on that error.... what would I do to calculate that?

Title: Euler's first order Error


问:
Coupled Diff Equation from Bose Einstein distribution
I am a student doing physics hons and have had very little experience in programming. This semester we are supposed to do a computational project in thermodynamics. I have to solve these two coupled diff eqns:

$$
\begin{aligned}
ω (p,T ) &= \frac{p ^2}{2m} + \frac{2}{N}\sum_q f(p − q)\,n (q) − \frac{1}{N^2}\sum_{s,t} f(s − t)\,n(s)\,n(t) \quad\text{and} \\
n(p) &= \frac{1}{\exp\left[\cfrac{ω (p,T ) − \mu}{kT}\right] − 1}
\end{aligned}
$$

$\omega$ is the energy per boson.

$p$ is the momentum of a boson.

$n(p)$ is the number of bosons in the state with momentum $p$.

$f$ is a function of the form
$$
  f(p)=\frac{1}{2} \left[\epsilon_0- \frac{p^2}{2m}\right]
$$

$\epsilon_0$ is elementary excitation energy at 0 K.

$N$ is total no of bosons.

$T$ is temp and $k$ is a constant.

Can someone guide me to any simple methods to generate some crude solution to this problem? Based on a paper: "Evaluation of speciﬁc heat for superﬂuid helium between 0 - 2.1 K based on nonlinear theory" By Shosuke Sasaki (arXiv:0807.1361v1 [cond-mat.other] 9 Jul 2008


问:
Application of vector extrapolation methods to convergence to a steady state solution
I'm working on a fluid solver using dual-time stepping and everything works really well, except the convergence in pseudo-time is slow. I'd like to accelerate the convergence. I know multigrid methods are typically used, but I am thinking about exploring another approach that is, hopefully, easier to implement.

I started reading the literature on sequence convergence acceleration. My sequence in this case is the updated solution to the Navier-Stokes equations in pseudo-time and the sequence, if all goes well, converges to a constant field (alternatively, the pseudo-time derivative converges to zero) as an infinite number of steps are taken. My reading in this area has led me to vector extrapolation methods and the related algorithms, but I am unsure of a few points.

For sake of argument, consider the Reduced Rank Extrapolation (RRE) method. The algorithm can be found in Table 3.1 of this report. In summary, you collect a set of solution vectors (solve your equations over a time step and store the updated values). Then, compute the derivative of the solution (or really just store the increment from the fluid solver since that is what we compute anyway) and the resulting vectors become columns in a big matrix. Find the QR factorization of this big matrix and solve a least-squares problem, and use the result to extrapolate your solution vector. Hopefully your convergence accelerates.

First, can these actually help in this case? I think they can, but I'm a bit fuzzy on it all.

Second, can these methods be applied locally or do they need to be applied for all cells in the simulation simultaneously? Everything I have read just shows solution vectors. But I don't know if it's $nvar*ncells$ long, or if I can do the QR factorization in each cell to get an extrapolated vector. There's a big difference in cost between point-wise factorization and total-domain factorization in massively parallel simulations. The linked report says it can be done on subdomains, which seems to imply it is at least a "regional" (as in, not pointwise but at least within a processor) QR factorization.


问:
Estimating the wind speed above the mojave desert
What equation would you all use to estimate the wind speed over the Mojave Desert in December? I'd like to model the speed of the wind up to say 50 thousand feet.

I read up about wind speeds on Wikipedia and they give an analytical formula for wind gradient:

$$\\ v_w(h) = v_{10} \cdot \left(  \frac {h} {h_{10}} \right)^ a$$

I also found this link:
http://digital-desert.com/mojave-preserve/geology/02.html
But it doesn't tell me how to vary altitude.

I would rather be using data that was taken, for example by NOAA, but couldn't find any. Does anybody have a link, algorithm, or recommendation of what I should do?


问:
Probabilistic model to approach problem that is usually dealt with linear programming
I have the following problem. Let's say we have $x_{jk}$ it is an expression value of gene $j$ in a sample $k$. It is the average of expression levels across the cell types $s_{ij}$, weighted by respective proportions $a_{ki}$ ($i = 1 \cdots N$, $N$ is the disease type):

$$
x_{jk} = \sum_{i=1}^{N} a_{ki}s_{ij}
$$

Generally this can be expressed as matrix form

$$
 X = AS
$$

What I want to do is to solve this equation

\begin{align}
&\min_{A}\ || AS- X||^{2}\\
&\text{subject to } \left\{
\begin{array}{c l}
    \sum_{i} a_{ki} = 1\\
    a_{ki} \ge0, \forall i
\end{array}\right.
\end{align}

Typically this is solved by using quadratic programming.  The  number of genes is  large (e.g. ~30K).

I'm wandering if there is any good probabilistic model to deal such kind of problem?


问:
Solving constrained BVP, singular Jacobian
The boundary value problem is

$$
\begin{cases}
\dot{x}_i =
\begin{cases}
(0.5D^{-1}\psi)_i, \text{ if }(0.5D^{-1}\psi)_i \le 0 \\
0 \text{, otherwise}
\end{cases}
\\
\dot{\psi} = 2\Sigma x \\
x(0) = x^0 \\
x(T) = 0
\end{cases}
$$

Here $x^0>0$, $T$ is given, $D>0$ is a diagonal matrix, $\Sigma>0$ is a symmetric matrix.  Actually the constraint is $x(t)\ge 0\forall t$, equivalent to $\dot{x}(t) \le 0$.

For a 2-dimension case, I try to solve it in MATLAB with this code

```<code>function [tint, traj] = liq2mon_v3(x0, D, sigma, T)

  tlow = 0;
  thigh = T;
  L = length(x0); % L=2 in this example
  nsegm = 100; % I know it is for final calc only, not for the mesh
  solinit = bvpinit(linspace(tlow, thigh, nsegm), @guess);
  sol = bvp5c(@bvp4ode, @bvp4bc, solinit);
  tint = linspace(tlow, thigh);
  z = deval(sol, tint);
  traj = z(1:L, :);

  function xpsi = guess(t)
    xpsi = zeros(1, L*2);
    xpsi(1:L) = (x0*(1.0 - 1.0/T*t))';
  end

  function dydt = bvp4ode(t, y)
    xcomp = y(1:L); % L=2 in this example
    psicomp = y(L+1:end);
    u    = D \ psicomp / 2.0;
    dx   = min(u, 0);
    dpsi = 2.0 * sigma * xcomp;
    dydt = [dx' dpsi'];
  end

  function res = bvp4bc(ya, yb)
    xcompa = ya(1:2);
    xcompb = yb(1:2);
    res = [(xcompa - x0)' (xcompb)'];
  end

end
</code>```

and I get an error

```<code>*Unable to solve the collocation equations -- a singular Jacobian  encountered
Error in ==> liq2mon_v3 at 13
sol = bvp5c(@bvp4ode, @bvp4bc, solinit)*
</code>```

However if I change the line

```<code>dx = min(u, 0);
</code>```

to

```<code>u(u>0) = 1e-13;
dx = u;
</code>```

I do get a solution, with $\min(\min(x)) \approx -8\times 10^{-14}$ for the same $x^0$.

Questions:

Is there a way to solve this system in MATLAB (or another app) so that the constraint is not violated?

Maybe this could be solved in Python, even with the current constraint violation?

Why is the Jacobian singular somewhere within MATLAB solver in this case?


问:
Does Conjugate Residual really have convergence properties similar to that of Conjugate Gradient?
I have coded up a toy implementation of Conjugate Residual and have been testing it.  Both wikipedia and the Saad claim that Conjugate Residual and Conjugate Gradient have similar convergence behavior.  For some examples I have been doing I haven't found this to be true.  For a $128 \times 128$ symmetric indefinite matrix, I don't get (machine) convergence at 400 iterations.  This is at odds with my understanding of Conjugate residual being a Krylov space method that should converge in at most 128 iterations.  (The condition number of this matrix is about 120 and every singular value  is repeated since the eigenvalues come in positive and negative pairs.)

I am guessing the statement about similar convergence properties refers to the performance on symmetric positive definite matrices, and numerical experiments are consistent with this.

So my question is, what determines the convergence rate of conjugate residual for a symmetric but indefinite matrix?  Also, is the conjugate residual method supposed to converge in at most $n$ iterations, where $n$ is the size of the system?


问:
Is it possible to run a Solver Foundation solver against a model containing linear and non-linear elements?
This is a follow up question to one I made previously about non-linear equations and ranged real numbers in Solver Foundation.

I acknowledge that where possible, rewriting a problem that is non-linear into linear could avoid some issues, however it is likely for various reasons this won’t be a viable option.

The requirements for the project I am involved with call for a set of constraints to be generated from a data source. The issue is the constraints could be a mix of linear or non-linear problems.

I have tried to come up with a set of example constraints and variables to test against that will hopefully give a clearer picture as to what I am trying to achieve. Does anyone know of a solver capable of handling a model of the kind to be generated by this combination of variables and constraints? We are currently using Solver foundation and have been trying various configurations and solver types.

Example:

Some parameters are created and have been bound to some arbitrary values:

```<code>   P1 = 1
   P2 = 2
   P3 = 3
</code>```

Some decisions with a domain specified as a finite range with a defined step value between each value in the range:

```<code>D1 (1..5) step of 0.1  //(e.g. 1, 1.1, 1.2 … 4.9, 5)
D2 (1..5) step of 1
D3 (10..20) step of 0.1
</code>```

The following example constraints simulate the kinds of constraints likely to be acquired from the data source. Next to each constraint is the expected result of running a solver with a model containing the one constraint, given the values of the above parameters and decision variables. Note that there is no concern with the quality of a solution, simply if there is a solution (feasible) or not (unfeasible).

```<code>C1: P1 = P2 * P3 (infeasible)
C2: P3 = P1 + P2 (feasible)
C3: D1 = P1 * P2 (feasible)
C4: D3 = P1 * P2 (infeasible)
C5: P1 = D1 * D2 (feasible)
C6: P1 = D3 * D3 (infeasible)
C7: D3 = D1 * D2 (feasible)
C8: D2 = D3 * D3 (infeasible)
C9: P3 = P2 / P1 (infeasible)
C10: -P3 = P1 * D1 (infeasible)
</code>```

Any pointers or advice would be most welcome :) moving from a UI background into an optimisation project has been character building to say the least.

Thanks,
d.


问:
How to speed convergence when optimizing a linear objective with nonlinear constraints?
I'm trying to learn how to do optimal rocket trajectory planning. I have a process that works but it converges very slowly; I'm looking for help to understand how to speed that up.

The optimal control problem

This is a simplified version of the linear tangent steering problem:

Imagine a point mass (the "rocket") moving in a 2D plane. No gravity or other external force acts on it. The rocket accelerates itself (up to a maximum), but for now at least it does not change mass. It can instantly rotate to face in any direction.

The goal is to get from a known initial state to a partially-specified final state in minimum time without accelerating above the allowed limit. For my example:

Initial position is at the origin

Initial velocity is zero

Final position is specified in the Y direction (height)

Final position is unspecified in the X direction (down-range)

Final velocity is zero in the Y direction (vertical)

Final velocity is specified in the X direction (horizontal)

Maneuver duration is unspecified (it's the variable being optimized)

Discretizing the problem

I'm breaking the trajectory up into N segments of equal duration. In each segment the rocket has a constant acceleration vector, so the optimal trajectory is approximated by a series of parabolic arcs put end-to-end, with matching positions and velocities at the boundaries between them.

The general approach I'm taking is to start with N very small; say two or three segments. Once I've converged on a solution I can double the number of segments by halving the segment duration and duplicating each segment. Then I iteratively improve that solution until it converges, and repeat.

The problem is simple enough that there are a variety of ways of representing the state and constraints. It's possible to just do the integration in closed form instead of having segment integration constraints, for instance, which gives you a very small number of constraints for the final conditions, plus all the acceleration constraints for each of the segments.

You could opt to connect the segments together with integration constraints instead (specifying that the velocity at the end of a parabolic arc should be consistent with the acceleration over the arc, and that the position should be consistent with the initial and final velocities on that segment). This is more typical and would enable you to go on and add gravity and so forth more readily.

Accelerations could be an explicit part of the state vector, or they could be derived by comparing the initial and final velocity of each segment.

Solving the problem

I'm solving this by taking a candidate solution and alternating two steps:

Move the solution in the direction of the constrained gradient by some amount

Move the solution toward the feasibility boundary if it's violating constraints

To compute the constrained gradient:

$$
\begin{eqnarray*}
\vec{x} &=& \text{state vector} \\
\vec{g} &=& \text{objective direction (the unconstrained objective gradient)} \\
\vec{n}_0, \vec{n}_1, \dotsc &=& \text{constraint directions} \\
c_0, c_1, \dotsc &=& \text{constraint errors} \\
\lambda_0, \lambda_1, \dotsc &=& \text{constraint multipliers} \\
\vec{p} &=& \text{constrained movement direction}
\end{eqnarray*}
$$

The $\vec{g}$ vector is just $-1$ in the duration variable and $0$ everywhere else, since I'm just trying to minimize time. The constraint normals $\vec{n}_i$ and errors $c_i$ are evaluated at the current state $\vec{x}$.

We want the movement direction $\vec{p}$ to be perpendicular to all active constraint normals, which we accomplish by subtracting out the parts of $\vec{g}$ that are not perpendicular to the constraints:

$$
\begin{eqnarray*}
\vec{p} &=& \vec{g} - \lambda_0 \vec{n}_0 - \lambda_1 \vec{n}_1 - \dotsb - \lambda_m \vec{n}_m \\
\vec{p}^\intercal \vec{n}_0 &=& 0 \\
\vec{p}^\intercal \vec{n}_1 &=& 0 \\
\vdots \\
\vec{p}^\intercal \vec{n}_m &=& 0
\end{eqnarray*}
$$

Substituting the equation for $\vec{p}$ into the $m$ equations for perpendicularity lets us solve for the constraint normal multipliers $\lambda_i$, which can then in turn be put into the first equation to yield $\vec{p}$.

Where I'm stuck

The next question is how far to move. $\vec{p}$ can be used as-is, but the constraints for this sort of problem are all essentially spherical, so as we reach the optimum they flatten out and movement slows to a crawl. I've tried scaling $\vec{p}$ up by an arbitrary amount but that doesn't help much, and seems inelegant.

I can think of a couple of possible things to do. I know what the expected rate of objective change will be, by doing a dot product between $\vec{g}$ and $\vec{p}$. If I could get a second derivative for the expected objective change I could guess how long of a step to take to reach the minimum. I'm getting all tangled up trying to figure out how to do that, though.

An even closer fit might be to get an entire acceleration vector for movement in the constrained gradient direction, so that we'd move along a parabola instead of a straight line. This would keep the solution closer to the feasibility surface, so the feasibility correction step would have less to do. I'm not too worried about that, though; the feasibility restoration step works pretty well.

Am I headed in the right direction for solving this? Is there a better way? Are there any good references for this sort of problem?

Thanks in advance.

References

Here are the primary things I've been working from (along with whatever I've been able to find on the Internet):

Practical Methods for Optimal Control and Estimation Using Nonlinear Programming, 2nd ed., by John Betts

Optimal Control Theory: An Introduction, by Donald Kirk


问:
ODE events to switch between 5 equations (friction model)
I am modelling a 1 dof spring-mass-damper system with friction. As first attempt I modelled the friction according to the simple Coulomb model (figure A here http://article.sapub.org/image/10.5923.j.mechanics.20130301.04_006.gif) and it worked fine.
Now I want to add the effect of the static friction when v = 0 m/s (figure C http://article.sapub.org/image/10.5923.j.mechanics.20130301.04_006.gif), but I cannot make it work. In fact I would expect to obtain something similar to the first model, but it is not so.

The equations I modelled are:

for the slipping case ($\dot{x}$ different from 0)

$ m\ddot{x} + c\dot{x} + kx = -m\ddot{y} - \mu_{d} N $ if $\dot{x} > 0$

$ m\ddot{x} + c\dot{x} + kx = -m\ddot{y} + \mu_{d} N $ if $\dot{x} < 0$

for the sticking case ($\dot{x} = 0$) I have 3 equations:

$ m\ddot{x} + c\dot{x} + kx = - m\ddot{y} + kx $ if $ |kx+m\ddot{y}| < \mu_{s} N$

If $ |kx+m\ddot{y}| > \mu_{s} N$ then the equations are

$ m\ddot{x} + c\dot{x} + kx = -m\ddot{y} -\mu_{s} N $ if $ kx+m\ddot{y} > 0$

$ m\ddot{x} + c\dot{x} + kx = -m\ddot{y} +\mu_{s} N $ if $ kx+m\ddot{y} < 0$

where $N$ is the normal force and $-m\ddot{y}$ represents the logarithmic sweep base excitation given as $\ddot{y} = A*\sin{(\phi(t))}$.

So, the starting equation is chosen as follows

```<code>% Initial conditions
q0 = [0;0];

% Choose starting equation
if q0(2) ~= 0
    flag1 = 0;
else
    flag1 = 1;
end
</code>```

Then, to switch from one equation to the other set I need events:

```<code>function [value,isterminal,direction] = events(t,q)
    value = [q(2),q(2),...]; % Find the times when velocity changes sign
    isterminal = [1,1,1]; % Stop the integration when the event functions
                        % are satisfied
    direction = [1,-1,0];
end
</code>```

The first 2 terms in the value look for changes in the velocity direction, but now I need to tell Matlab to choose between one of the 3 equations of the case $\dot{x} = 0$.
The third condition I would add is

```<code>abs(k*q(1)+m*A*sin(2*pi*(60*f0/(R*log(2.))*(2.^(R/60*t)-1))))-mus*N
</code>```

But then if I try to add the condition for the 2 remaining equations I get an error on input matrix dimension here

```<code>[t,q,te,qe,ie] = ode45(odefun,tspan,q0(:,1),options);
</code>```


问:
Eigenvalue-style optimization with quadratic constraints
Suppose $A\in\mathbb{R}^{n\times n}$ is symmetric and positive definite and that we have several symmetric matrices $B_i\in\mathbb{R}^{n\times n}$ that are low-rank and indefinite.  I need an algorithm for solving the following optimization problem:
$$
\begin{align*}
\min_{x\in\mathbb{R}^n}\ & x^\top Ax\\
\textrm{s.t.}\ & \|x\|_2^2=1\\
& x^\top B_i x=0\ \forall i\in\{1,\ldots,k\}.
\end{align*}
$$
Any suggestions?

I am hoping it is somehow equivalent to an eigenvalue problem.  Assorted observations:

If the second constraint is removed, the output is the minimal eigenvector of $A$.

The Lagrange multiplier expression is $Ax=\lambda x + \sum_i \mu_iB_ix$.  So, if we somehow know the dual variables $\{\mu_i\}_{i=1}^k$, then $x$ is an eigenvector of $A-\sum_i\mu_i B_i$ with minimal eigenvalue $\lambda$.  So, we could think of $x$ as a very nonlinear function $x(\mu_1,\ldots,\mu_k)$.

If one of the $B_i$'s is positive definite, then this problem is not satisfiable (implies $x=0$).  But for my specific problem I can prove the feasible region is nonempty.

This can be approached using semidefinite relaxation, but this would be very slow for large $n$.


问:
Applications of algorithm for solving systems of equations with uncertainty
We have been developing algorithms for detecting "robust" zeros of multidimensional functions $f: X\to\Bbb R^n$ where $X$ is an $m$-dimensional domain in $\Bbb R^m$. More precisely, for a given $f$, we can compute an $r>0$ such that every continuous $g$, $\|g-f\|<r$ has a nonempty zero set in $X$ (under some quite weak conditions). This can be used for proving that the initial map $f$ has a zero even in cases when it is given only approximately, such as a list of function values in a multidimensional grid and a Lipschitz constant.

More involved robust properties of the zero set of $f$ could be analysed. For instance we can  find "robust optima" over the zero set with the uncertainty when an additional objective function is given.

Can these algorithms have some natural applications in engineering and/or scientific computations? In some sense, it generalizes solving systems of equations to the case where the input is uncertain.

The mathematics behind the algorithm is nice and nontrivial (and has been the main motivating force for us so far), but the further development makes the best sense if we had specific problems to solve.