问:
How to solve global equation in Comsol only once?
Let's assume that I have a comsol model for solving time-dependent problem. In order to find initial value for one variable I need to solve simple algebraic equation before starting simulation, therefore I added Global ODEs and DAEs module to the model. Now it works good, but it seems like equation from Global ODEs and DAEs module is solved on each time iteration and its obviously redundant calculations, because this equation need to be solved only once.
Questions:

1) Does equation from Global ODEs and DAEs module is solved on each time iterarion or it is only my suggestion? Maybe there is some way how I can check it?

2) If this equation is really solved on each iteration again how can I adjust my solver in order to solve it only once?

ContentLicense: CC BY-SA 3.0


问:
Problem with cell size and boundary conditions in transient cylindrical conduction
I am attempting to model the steady state behavior of a cylinder using the finite volume method (FVM) subjected to a variety of boundary conditions in Matlab. First off, I am treating the cylinder as being axisymmetric so I am only determine the temperature profile in the r-Z plane. I set up a 2D grid so that the r coordinate is on the 'x' axis and the Z coordinate is the 'y' axis. Since the cylinder is axisymmetric then the heat flux at r = 0 is 0 (insulated). I am also applying an insulated boundary condition to the bottom surface at Z=0. The side surface has convection with a constant convection coefficient and air temperature of 310 K. The top of the cylinder, z=Z, there is a heat flux of 1000 W/m^2, convection from the air at 310K, and also surface to ambient radiation with the surroundings being 298 K. The initial temperature of cylinder is 298 K (initial guess value). The problem that I am having is that the temperature distribution matches Comsol's results with a grid of 40x40 and a thermal conductivity of 1. When I increases the number of cells to 50x50, 60x60, etc, or if i change the thermal conductivity to something other than 1 then the temperature profile is still what it should be but all the temperature are significantly lower than expected. My control loop in Matlab is as follows:

Use initial values and problem constants to solve for initial surface temperature of top and side.

Use surfacte temperatures and problem geomety to determine aE, aN, aS, aW, aP, and b for every cell in the grid.

Use line by line TDMA in both directions to determine the temperature profile.

Calculate the error between the new temperature profile and the previous one.

Update the surface temperatures and aP and b because of the nonlinear radiation term that uses previous iterate values.

Perform line by line TDMA again and check for convergence. If the error is not less than 0.001 then repeat the procedure.

For the boundary conditions, I used an energy balance for the top surface and the side surface. The top surface includes surface to ambient radiation and I calculated the surface temperature to be
$$T_{b}=\frac{\alpha G+hT_{f}+\frac{kT_{p}}{\delta z}+\epsilon \sigma(T_{sur}^{4}+3T_{b}^{*4})}{h+4 \epsilon \sigma T_{b}^{*3}+\frac{k}{\delta z}}$$

I also do the same thing with the surface temperature on the side the is experiencing heating through convection. Since the top surface temperature depends on its previous iterate value, I use an initial guess value of 298 for the top surface. I then calculate the initial surface temperature based on that and the equation i provided for the first iteration,

The problem lies in the fact that the surface tempreature depends on both the thermal conductivity and the grid spacing. When I use many cells like 80x80, the top surface temperature for the first TDMA iteration is almost 298 which is hard to believe because its being heated by a 1000 W/m^2 heat flux and its also being heated by the 310 K air temperature. This causes the first iteration to have very little error so convergence only produces a max temperature of like 325 when it should be 450. But if I use a grid of like 40x40 then the surface temperature for the first iteration is much much higher, and convergence causes the final temperature profile to look almost exactly like the expected.

The same exact thing happens with the thermal conductivity, if i change the value then the initial surface temperatures change dramatically which causes the convergence values to be much lower than the expected.

It seems like it all has to do with that initial surface temperature but I am not sure how else to treat it. Any help would be amazing as I have spent many nights til 5am trying to figure out how I could fix this. My contact email is in my profile if you would like to see my matlab code of if you would like to talk to me more about my problem.


问:
assembly matrices in finite element method
I'm trying to construct the right–hand side of my 2D Poisson's equation in Matlab. I used the vertex rule in order to approximate the integral:

```<code>function F1 = rightVer(p,t,f)

%area triangle
it1=t(1,:);
it2=t(2,:);
it3=t(3,:);
area = ((p(1,it2)-p(1,it1)).*(p(2,it3)-p(2,it1))-(p(1,it3)-p(1,it1)).*(p(2,it2)-p(2,it1)))/2;

area = abs(area)/3;
vt1 = area.*f(p(1,it1),p(2,it1));
vt2 = area.*f(p(1,it2),p(2,it2));
vt3 = area.*f(p(1,it3),p(2,it3));

vt = [vt1';vt2';vt3'];
elem = t';
elem(:,4)=[];
N=size(p,2);
F1  =  accumarray(elem(:),vt,[N  1]);
end
</code>```

I also construct it with the pde toolbox with the function assema() but the two f which I obtained are the same except for the last 4 numbers. Can anyone help me?

ContentLicense: CC BY-SA 4.0


问:
Characteristic length of differential element of cylinder surface?
I am trying to find the Nusselt number for a small element of the outside of a cylinder that has a height of $\Delta z$. I found the average Grashof number of a surface as

$$Gr_{L}=\frac{\beta \rho (T_{s}-T_{\infty})L_{c}^{3} }{\nu^{2}}$$

If I wanted to find $Gr_{x}$ of a certain element, would I use a characteristic length of $\Delta z$ or could I find the z-location of the elements centroid and use that as the $L_{c}$. For example, If the cell was located at a z coordinate of 1 m, would that be the characteristic length or would the characteristic length be the height of the differential element?

I was thinking that I could treat each element as a small isothermal cylinder and have the characteristic length of each element be $\Delta z$

Also, does anyone know the a Nusselt number correlation for a vertical cylinder with a Prandtl number of less than one? I have only found one for when it is greater than 1. Thank you


问:
Updating matrix diagonal with Woodbury matrix identity and maintaining numerical accuracy
I have a dense matrix A and its corresponding inverse $A^{-1}$. The Woodbury matrix identity states:

$$ (A + UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1} $$

I wish to perform small updates to the diagonal of $A$ and compute corresponding $A^{-1}$. The updates can actually be described as a constant $k$ multiplying an identity matrix: $kI$ and in my case $U=V^T$ where $V$ is a tall matrix with just ones and zeros that selects the desired diagonal positions of $A$ to modify.

Everything is fine so far, except $k$ needs to be quite large (e.g., $10^6$). I am noticing that for large $k$, changes made are not "reversible" without significant error. What I mean by reversible is performing two update steps, one with $kI$ on the original $A$ and another with $-kI$ on the updated matrix to recover the original $A^{-1}$. If $k$ is small enough, it works fine, but otherwise I run into serious accuracy issues and the resulting matrix does not match the initial one closely. For my target application, I need to be able to apply $k$ to diagonal elements or $-k$ at will, so my question is: Is there a clever way to do this which preserves numerical accuracy even for large values of $k$?

Note: I am using double precision for all values stored and for computation. Furthermore, $k$ is a real scalar, but the matrix $A$ is complex, although that should not really matter in this case.


问:
The traveling salesman problem - Using Space Renormalization

Image attached is where I am at the moment. Blue dots=points/cities, Black x's represent central points in each box that contains at least one city, and pink dots represent the midpoint of these centre points

I am doing a short project using Lab Windows to approximate the minimal path between points on a 2D plane. (Graphically only, I don't need an actual path length)

To do this I am using real space re-normalization.

http://scitation.aip.org/docserver/fulltext/aip/journal/cip/10/6/1.168588.pdf?expires=1387106816&id=id&accname=guest&checksum=350319AE8AA74E9002D87940052A0140 - An article by U. Yoshiyuki - "Solving the Traveling Salesman Problem By a Statistical Physics Method" describes exactly what I am trying to do

So far I have:
Plotted random integer valued points on a graph in my user interface

Divided my plotting space evenly into boxes, i.e. a 2x2 graph (Just by plotting lines)

Created a new point in the centre of each box if there was a random point generated in that box

Found the nearest edge point (for the 2x2 graph only)

```<code>     /*
     * Travelling Salesman Problem
     * A Graphical Approximation using Space Re-normalization

     * Assuming there are  3 or more cities
     */

     #include "toolbox.h"
     #include <cvirte.h>
     #include <userint.h>
     #include "TSP.h"

     static int panelHandle;
     void draw_borders(int,int,int,int);/*function prototype*/

     int main (int argc, char *argv[])
     {
if (InitCVIRTE (0, argv, 0) == 0)
    return -1;  /* out of memory */
if ((panelHandle = LoadPanel (0, "TSP.uir", PANEL)) < 0)
    return -1;
DisplayPanel (panelHandle);
RunUserInterface ();
DiscardPanel (panelHandle);
return 0;
     }

     int CVICALLBACK Quit (int panel, int control, int event,
    void *callbackData, int eventData1, int eventData2)
     {
switch (event)
         {
    case EVENT_COMMIT:
        QuitUserInterface (0);
        break;
}
return 0;
     }

     int CVICALLBACK ClearPlot (int panel, int control, int event,
    void *callbackData, int eventData1, int eventData2)
     {
switch (event)
{
    case EVENT_COMMIT:
        DeleteGraphPlot(panelHandle,PANEL_2X2GRAPH,-1,1);
        DeleteGraphPlot(panelHandle,PANEL_4X4GRAPH,-1,1);
        DeleteGraphPlot(panelHandle,PANEL_8X8GRAPH,-1,1);
        DeleteGraphPlot(panelHandle,PANEL_16X16GRAPH,-1,1);
        DeleteGraphPlot(panelHandle,PANEL_GRAPH_PATH,-1,1);
        break;
}
return 0;
     }

     int CVICALLBACK Go (int panel, int control, int event,
    void *callbackData, int eventData1, int eventData2)
     {
switch (event)
{
    case EVENT_COMMIT:

        int x[1024],y[1024]; /* City co-ordinates */
        int i,j,k,num=10; /* num is the number of cities to be plotted */
        int xi=0,xf=1024,yi=0,yf=1024; /* boundary points om graph */
        int flag=0;

        /* Structures to store central points for each box containing a city*/
        typedef struct{
            int x,y;
        } Segment1;
        Segment1 Seg1[100][100];

        typedef struct{
            int x,y;
        } Segment2;
        Segment2 Seg2[100][100];

        typedef struct{
            int x,y;
        } Segment3;
        Segment3 Seg3[100][100];

        typedef struct{
            int x,y;
        } Segment4;
        Segment4 Seg4[100][100];

        /*2x2*/
        int quad1[2][2]; /*count for cities*/
        int grid1x2d[2][2],grid1y2d[2][2];/*grid containing x and y coordinate of central point locations in each box which contains a city*/
        int grid1xx[4], grid1yy[4];/*same grid stretched to 1d for plotting*/
        int mid1x2d, mid1y2d;

        /*4x4*/
        int quad2[4][4];
        int grid2x2d[4][4],grid2y2d[4][4];
        int grid2xx[16], grid2yy[16];

        /*8x8*/
        int quad3[8][8];
        int grid3x2d[8][8],grid3y2d[8][8];
        int grid3xx[64], grid3yy[64];

        /*16x16*/
        int quad4[16][16];
        int grid4x2d[16][16],grid4y2d[16][16];
        int grid4xx[256], grid4yy[256];

        draw_borders(xi,xf,yi,yf);/*Draw borders for clarity*/

        /* Plot Cities */
        for(i=0;i<num;i++){

            /*represent each city as a random point*/
            x[i] = (int)Random(0,1024);
            y[i] = (int)Random(0,1024);
            PlotPoint(panelHandle,PANEL_2X2GRAPH,x[i],y[i],VAL_CROSS,VAL_BLUE); /*plot city on each graph*/
            PlotPoint(panelHandle,PANEL_4X4GRAPH,x[i],y[i],VAL_CROSS,VAL_BLUE);
            PlotPoint(panelHandle,PANEL_8X8GRAPH,x[i],y[i],VAL_CROSS,VAL_BLUE);
            PlotPoint(panelHandle,PANEL_16X16GRAPH,x[i],y[i],VAL_CROSS,VAL_BLUE);
            PlotPoint(panelHandle,PANEL_GRAPH_PATH,x[i],y[i],VAL_CROSS,VAL_BLUE);

            /*Print out city coordinates*/
            char resultx[100];
            sprintf(resultx,"x co-ordinate of city %d is %d\n",i+1,x[i]);
            SetCtrlVal(panelHandle,PANEL_CITYCOORDINATES,resultx);

            char resulty[100];
            sprintf(resulty,"y co-ordinate of city %d is %d\n",i+1,y[i]);
            SetCtrlVal(panelHandle,PANEL_CITYCOORDINATES,resulty);
        }

        SetCtrlVal(panelHandle,PANEL_MESSAGE,"Cities Plotted");

        ////////////*2x2 Graph*////////////////////////////////////////////////////////////////////////////
        /* First find centre point of each box which contains a city
         * Second find midpoint of centre point locations of each box that contains at least one city
         */

        /* Set count in all 4 boxes to be zero*/
        for (i=0;i<2;i++){ /*x direction*/
            for (j=0; j<2; j++) /*y direction*/
                quad1[i][j]=0;
        }

        /* For every city, add 1 to count */
        for (i=0;i<num;i++)
            quad1[x[i]/512][y[i]/512]++;

        /* Set up a grid containing the x,y coordinates of central point locations of each box that contains at least one city*/
        for (i=0;i<2;i++){
            for (j=0; j<2; j++){
                if (quad1[i][j]>0){
                    grid1x2d[i][j]=i*512+256;  grid1y2d[i][j] =j*512+256;

                    /* Assign these values to a structure to call on them later when finding midpoint*/
                    Seg1[i][j].x = grid1x2d[i][j]; Seg1[i][j].y = grid1y2d[i][j];
                }
                else{
                    /* Set up a grid containing dummy x,y coordinates if there are no cities in a particular box */
                    grid1x2d[i][j]=-1;  grid1y2d[i][j] =-1;

                    Seg1[i][j].x = grid1x2d[i][j]; Seg1[i][j].y = grid1y2d[i][j];
                }
            }
        }

        /* Find edge points which are found if representative points are connected*/
        /* If the 2X2 Graph has an empty box then there will be an edgepoint in centre of graph*/
        for(i=0;i<2;i++){
            for(j=0;j<2;j++){
                if (quad1[i][j]<=0){
                                      PlotPoint(panelHandle,PANEL_2X2GRAPH,512,512,VAL_BOLD_X,VAL_DK_GREEN);
                    PlotPoint(panelHandle,PANEL_4X4GRAPH,512,512,VAL_BOLD_X,VAL_MAGENTA);
                    PlotPoint(panelHandle,PANEL_8X8GRAPH,512,512,VAL_BOLD_X,VAL_MAGENTA);
                    PlotPoint(panelHandle,PANEL_16X16GRAPH,512,512,VAL_BOLD_X,VAL_MAGENTA);
                }
            }
        }

        for(i=0;i<2;i++){
            j=0;
            if( quad1[i][j]>0 && quad1[i][j+1]>0){
                mid1x2d= (Seg1[i][j].x + Seg1[i][j+1].x) /2 ;
                mid1y2d= (Seg1[i][j].y + Seg1[i][j+1].y) /2 ;
                PlotPoint(panelHandle,PANEL_2X2GRAPH,mid1x2d,mid1y2d,VAL_BOLD_X,VAL_DK_GREEN);
                PlotPoint(panelHandle,PANEL_4X4GRAPH,mid1x2d,mid1y2d,VAL_BOLD_X,VAL_MAGENTA);
                PlotPoint(panelHandle,PANEL_8X8GRAPH,mid1x2d,mid1y2d,VAL_BOLD_X,VAL_MAGENTA);
                PlotPoint(panelHandle,PANEL_16X16GRAPH,mid1x2d,mid1y2d,VAL_BOLD_X,VAL_MAGENTA);
            }
        }
        for(j=0;j<2;j++){
            i=0;
            if( quad1[i][j]>0 && quad1[i+1][j]>0){
                mid1x2d= (Seg1[i][j].x + Seg1[i+1][j].x) /2 ;
                mid1y2d= (Seg1[i][j].y + Seg1[i+1][j].y) /2 ;
                PlotPoint(panelHandle,PANEL_2X2GRAPH,mid1x2d,mid1y2d,VAL_BOLD_X,VAL_DK_GREEN);
                PlotPoint(panelHandle,PANEL_4X4GRAPH,mid1x2d,mid1y2d,VAL_BOLD_X,VAL_MAGENTA);
                PlotPoint(panelHandle,PANEL_8X8GRAPH,mid1x2d,mid1y2d,VAL_BOLD_X,VAL_MAGENTA);
                PlotPoint(panelHandle,PANEL_16X16GRAPH,mid1x2d,mid1y2d,VAL_BOLD_X,VAL_MAGENTA);
            }
        }

        /*stretch into 1d array for plotting*/
        for (k=0;k<4;k++){
            grid1xx[k]=grid1x2d[k%2][(int)k/2];
        }
        for (k=0;k<4;k++){
            grid1yy[k]=grid1y2d[k%2][(int)k/2];
        }

        PlotXY (panelHandle, PANEL_2X2GRAPH, grid1xx , grid1yy ,4, VAL_INTEGER, VAL_INTEGER, VAL_SCATTER,VAL_BOLD_X, VAL_THIN_LINE, 1, VAL_BLACK);
        ////////////////////////////////////////////////////////////////////////////////////////////////////

        ////////////*Optimal Paths*////////////////////////////////////////////////////////////////////////////

        ////////////////////////////////////////////////////////////////////////////////////////////////////

        ////////////*4x4 Graph*////////////////////////////////////////////////////////////////////////////
        /* Set count in all 4 boxes to be zero*/
        for (i=0;i<4;i++){ /*x direction*/
            for (j=0; j<4; j++) /*y direction*/
                quad2[i][j]=0;
        }

        /* For every city, add 1 to count */
        for (i=0;i<num;i++)
            quad2[x[i]/256][y[i]/256]++;

        /* Set up a grid containing the x,y coordinates of central point locations of each box that contains at least one city*/
        for (i=0;i<4;i++){
            for (j=0; j<4; j++){
                if (quad2[i][j]>0){
                    grid2x2d[i][j]=i*256+128;  grid2y2d[i][j] =j*256+128;

                    /* Assign these values to a structure to call on them later when finding midpoint*/
                    Seg2[i][j].x = grid2x2d[i][j]; Seg2[i][j].y = grid2y2d[i][j];
                }
                else{
                    /* Set up a grid containing dummy x,y coordinates if there are no cities in a particular box */
                    grid2x2d[i][j]=-1;  grid2y2d[i][j] =-1;

                    Seg2[i][j].x = grid2x2d[i][j]; Seg2[i][j].y = grid2y2d[i][j];
                }
            }
        }

        /*stretch into 1d array for plotting*/
        for (k=0;k<16;k++)
            grid2xx[k]=grid2x2d[k%4][(int)k/4];

        for (k=0;k<16;k++)
            grid2yy[k]=grid2y2d[k%4][(int)k/4];

        PlotXY (panelHandle, PANEL_4X4GRAPH, grid2xx , grid2yy ,16, VAL_INTEGER, VAL_INTEGER, VAL_SCATTER, VAL_BOLD_X, VAL_THIN_LINE, 1, VAL_BLACK);
        ////////////////////////////////////////////////////////////////////////////////////////////////////

        ////////////*8x8 Graph*////////////////////////////////////////////////////////////////////////////
        /* Set count in all 4 boxes to be zero*/
        for (i=0;i<8;i++){ /*x direction*/
            for (j=0; j<8; j++) /*y direction*/
                quad3[i][j]=0;
        }

        /* For every city, add 1 to count */
        for (i=0;i<num;i++)
            quad3[x[i]/128][y[i]/128]++;

        /* Set up a grid containing the x,y coordinates of central point locations of each box that contains at least one city*/
        for (i=0;i<8;i++){
            for (j=0; j<8; j++){
                if (quad3[i][j]>0){
                    grid3x2d[i][j]=i*128+64;  grid3y2d[i][j] =j*128+64;

                    /* Assign these values to a structure to call on them later when finding midpoint*/
                    Seg3[i][j].x = grid3x2d[i][j]; Seg3[i][j].y = grid3y2d[i][j];
                }
                else{
                    /* Set up a grid containing dummy x,y coordinates if there are no cities in a particular box */
                    grid3x2d[i][j]=-1;  grid3y2d[i][j] =-1;

                    Seg3[i][j].x = grid3x2d[i][j]; Seg3[i][j].y = grid3y2d[i][j];
                }
            }
        }

        /*stretch into 1d array for plotting*/
        for (k=0;k<64;k++)
            grid3xx[k]=grid3x2d[k%8][(int)k/8];

        for (k=0;k<64;k++)
            grid3yy[k]=grid3y2d[k%8][(int)k/8];

        PlotXY (panelHandle, PANEL_8X8GRAPH, grid3xx , grid3yy ,64, VAL_INTEGER, VAL_INTEGER, VAL_SCATTER, VAL_BOLD_X, VAL_THIN_LINE, 1, VAL_BLACK);
        ////////////////////////////////////////////////////////////////////////////////////////////////////

        ////////////*16x16 Graph*////////////////////////////////////////////////////////////////////////////
        /* Set count in all 4 boxes to be zero*/
        for (i=0;i<16;i++){ /*x direction*/
            for (j=0; j<16; j++) /*y direction*/
                quad4[i][j]=0;
        }

        /* For every city, add 1 to count */
        for (i=0;i<num;i++)
            quad4[x[i]/64][y[i]/64]++;

        /* Set up a grid containing the x,y coordinates of central point locations of each box that contains at least one city*/
        for (i=0;i<16;i++){
            for (j=0; j<16; j++){
                if (quad4[i][j]>0){
                    grid4x2d[i][j]=i*64+32;  grid4y2d[i][j] =j*64+32;

                    /* Assign these values to a structure to call on them later when finding midpoint*/
                    Seg4[i][j].x = grid4x2d[i][j]; Seg4[i][j].y = grid4y2d[i][j];
                }
                else{
                    /* Set up a grid containing dummy x,y coordinates if there are no cities in a particular box */
                    grid4x2d[i][j]=-1;  grid4y2d[i][j] =-1;

                    Seg4[i][j].x = grid4x2d[i][j]; Seg4[i][j].y = grid4y2d[i][j];
                }
            }
        }

        /*stretch into 1d array for plotting*/
        for (k=0;k<256;k++)
            grid4xx[k]=grid4x2d[k%16][(int)k/16];

        for (k=0;k<256;k++)
            grid4yy[k]=grid4y2d[k%16][(int)k/16];

        PlotXY (panelHandle, PANEL_16X16GRAPH, grid4xx , grid4yy ,256, VAL_INTEGER, VAL_INTEGER, VAL_SCATTER,VAL_BOLD_X, VAL_THIN_LINE, 1, VAL_BLACK);
        ////////////////////////////////////////////////////////////////////////////////////////////////////

        break;
}
return 0;
</code>```

}

```<code>     void draw_borders(int xi,int xf,int yi,int yf){ /*function definition*/

        /*Draw borders for clarity*/

        /*2x2*/
        PlotLine(panelHandle,PANEL_2X2GRAPH,xi,(yf+yi)/2,xf,(yf+yi)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_2X2GRAPH,(xf+xi)/2,yi,(xf+xi)/2,yf,VAL_RED);

        /*4x4*/
        PlotLine(panelHandle,PANEL_4X4GRAPH,xi,(yf+yi)/2,xf,(yf+yi)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_4X4GRAPH,(xf+xi)/2,yi,(xf+xi)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_4X4GRAPH,(xi+(xf+xi)/2)/2,yi,(xi+(xf+xi)/2)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_4X4GRAPH,(xf+(xf+xi)/2)/2,yi,(xf+(xf+xi)/2)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_4X4GRAPH,xi,(yf+(yf+yi)/2)/2,xf,(yf+(yf+yi)/2)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_4X4GRAPH,xi,(yi+(yf+yi)/2)/2,xf,(yi+(yf+yi)/2)/2,VAL_RED);

        /*8x8*/
        PlotLine(panelHandle,PANEL_8X8GRAPH,xi,(yf+yi)/2,xf,(yf+yi)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,(xf+xi)/2,yi,(xf+xi)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,(xi+(xf+xi)/2)/2,yi,(xi+(xf+xi)/2)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,(xf+(xf+xi)/2)/2,yi,(xf+(xf+xi)/2)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,xi,(yf+(yf+yi)/2)/2,xf,(yf+(yf+yi)/2)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,xi,(yi+(yf+yi)/2)/2,xf,(yi+(yf+yi)/2)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,(xi+(xi+(xf+xi)/2)/2)/2,yi,(xi+(xi+(xf+xi)/2)/2)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,(xf+(xi+(xf+xi)/2)/2)/2,yi,(xf+(xi+(xf+xi)/2)/2)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,(xf+(xf+(xf+xi)/2)/2)/2,yi,(xf+(xf+(xf+xi)/2)/2)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,(xi+(xf+(xf+xi)/2)/2)/2,yi,(xi+(xf+(xf+xi)/2)/2)/2,yf,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,xi,(yi+(yi+(yf+yi)/2)/2)/2,xf,(yi+(yi+(yf+yi)/2)/2)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,xi,(yi+(yf+(yf+yi)/2)/2)/2,xf,(yi+(yf+(yf+yi)/2)/2)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,xi,(yf+(yf+(yf+yi)/2)/2)/2,xf,(yf+(yf+(yf+yi)/2)/2)/2,VAL_RED);
        PlotLine(panelHandle,PANEL_8X8GRAPH,xi,(yf+(yi+(yf+yi)/2)/2)/2,xf,(yf+(yi+(yf+yi)/2)/2)/2,VAL_RED);

        /*16x16*/
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,(yf+yi)/2,xf,(yf+yi)/2,VAL_RED);//y=midline
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xf+xi)/2,yi,(xf+xi)/2,yf,VAL_RED);//x=midline
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xi+(xf+xi)/2)/2,yi,(xi+(xf+xi)/2)/2,yf,VAL_RED);//x=1/4 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xf+(xf+xi)/2)/2,yi,(xf+(xf+xi)/2)/2,yf,VAL_RED);//x=3/4 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,(yf+(yf+yi)/2)/2,xf,(yf+(yf+yi)/2)/2,VAL_RED);//y=3/4 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,(yi+(yf+yi)/2)/2,xf,(yi+(yf+yi)/2)/2,VAL_RED);//y=1/4 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xi+(xi+(xf+xi)/2)/2)/2,yi,(xi+(xi+(xf+xi)/2)/2)/2,yf,VAL_RED);//x=1/8 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xf+(xi+(xf+xi)/2)/2)/2,yi,(xf+(xi+(xf+xi)/2)/2)/2,yf,VAL_RED);//x=5/8 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xf+(xf+(xf+xi)/2)/2)/2,yi,(xf+(xf+(xf+xi)/2)/2)/2,yf,VAL_RED);//x=7/8 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xi+(xf+(xf+xi)/2)/2)/2,yi,(xi+(xf+(xf+xi)/2)/2)/2,yf,VAL_RED);//x=3/8 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,(yi+(yi+(yf+yi)/2)/2)/2,xf,(yi+(yi+(yf+yi)/2)/2)/2,VAL_RED);//y=1/8 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,(yi+(yf+(yf+yi)/2)/2)/2,xf,(yi+(yf+(yf+yi)/2)/2)/2,VAL_RED);//y=3/8 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,(yf+(yf+(yf+yi)/2)/2)/2,xf,(yf+(yf+(yf+yi)/2)/2)/2,VAL_RED);//y=7/8 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,(yf+(yi+(yf+yi)/2)/2)/2,xf,(yf+(yi+(yf+yi)/2)/2)/2,VAL_RED);//y=3/8 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi+((xi+(xi+(xf+xi)/2)/2)/2)/2,yi,xi+((xi+(xi+(xf+xi)/2)/2)/2)/2,yf,VAL_RED);//x=1/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi+((xi+(xf+(xf+xi)/2)/2)/2)/2,yi,xi+((xi+(xf+(xf+xi)/2)/2)/2)/2,yf,VAL_RED);//x=3/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi+((xf+(xi+(xf+xi)/2)/2)/2)/2,yi,xi+((xf+(xi+(xf+xi)/2)/2)/2)/2,yf,VAL_RED);//x=5/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi+((xf+(xf+(xf+xi)/2)/2)/2)/2,yi,xi+((xf+(xf+(xf+xi)/2)/2)/2)/2,yf,VAL_RED);//x=7/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xi+xf)/2+((xi+(xi+(xf+xi)/2)/2)/2)/2,yi,(xi+xf)/2+((xi+(xi+(xf+xi)/2)/2)/2)/2,yf,VAL_RED);//x=9/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xi+xf)/2+((xi+(xf+(xf+xi)/2)/2)/2)/2,yi,(xi+xf)/2+((xi+(xf+(xf+xi)/2)/2)/2)/2,yf,VAL_RED);//x=11/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xi+xf)/2+((xf+(xi+(xf+xi)/2)/2)/2)/2,yi,(xi+xf)/2+((xf+(xi+(xf+xi)/2)/2)/2)/2,yf,VAL_RED);//x=13/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,(xi+xf)/2+((xf+(xf+(xf+xi)/2)/2)/2)/2,yi,(xi+xf)/2+((xf+(xf+(xf+xi)/2)/2)/2)/2,yf,VAL_RED);//x=15/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,yi+((yi+(yi+(yf+yi)/2)/2)/2)/2,xf,yi+((yi+(yi+(yf+yi)/2)/2)/2)/2,VAL_RED);//y=1/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,yi+((yi+(yf+(yf+yi)/2)/2)/2)/2,xf,yi+((yi+(yf+(yf+yi)/2)/2)/2)/2,VAL_RED);//y=3/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,yi+((yf+(yf+(yf+yi)/2)/2)/2)/2,xf,yi+((yf+(yf+(yf+yi)/2)/2)/2)/2,VAL_RED);//y=7/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,yi+((yf+(yi+(yf+yi)/2)/2)/2)/2,xf,yi+((yf+(yi+(yf+yi)/2)/2)/2)/2,VAL_RED);//y=5/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,((yi+yf)/2)+((yf+(yi+(yf+yi)/2)/2)/2)/2,xf,((yi+yf)/2)+((yf+(yi+(yf+yi)/2)/2)/2)/2,VAL_RED);//y=13/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,((yi+yf)/2)+((yf+(yf+(yf+yi)/2)/2)/2)/2,xf,((yi+yf)/2)+((yf+(yf+(yf+yi)/2)/2)/2)/2,VAL_RED);//y=15/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,((yi+yf)/2)+((yi+(yf+(yf+yi)/2)/2)/2)/2,xf,((yi+yf)/2)+((yi+(yf+(yf+yi)/2)/2)/2)/2,VAL_RED);//y=11/16 line
        PlotLine(panelHandle,PANEL_16X16GRAPH,xi,((yi+yf)/2)+((yi+(yi+(yf+yi)/2)/2)/2)/2,xf,((yi+yf)/2)+((yi+(yi+(yf+yi)/2)/2)/2)/2,VAL_RED);//y=9/16 line
</code>```

}

Now I need to code up the following:

Save the optimal paths for any 2x2 box as a look up table

Draw the optimal path for each graph, and save where it hits the edge of its box. This edge point must then be used in the next step (when I make the graph 8x8, then 16x16, etc. until I have found minimum path to a good approximation)

And how can I pass through the graph in an ordered way to do this, and how can I save the optimal paths?

Thanks


问:
how to find outward normal for robin codition

I wrote the code to fix this problem; now  I should validate it using a test function to see the error in my code. I can't figure out how starting from u can derive the function g on the edge, in particular I don't manage to compute the outward normal.

This is my code to calculate g: https://dpaste.de/hzYm (p,t,e are output of initmesh)


问:
how to visualize velocity from Lagrangian particle tracking method
I have few particle's position and velocity information as a function of time.
Particle is following a trajectory path, while doing so, one particle may disappear, other particle may appear.

Particle position can be represented as a point and velocity as a vector at that point. If we think of particle position, it is nothing but a point in the flow path or pipe or channel
(I lack that information or I dont know about it), so I want to reverse engineer point position as kind of geometry position and map or interpolate velocity contour on those points, and see how flow is changing with time.

I dont want to know how to visualize using Matlab or Python.

is there any example or link in Matlab or Python that has reversed engineer Lagrangian particle tracking method.


问:
Quadratic programming problem involving permutation matrices
Does anyone know a good algorithm for quickly finding an approximate solution to the following problem?

Given two square matrices $A$ and $B$, minimize $\| P A P^\top - B \|$ over all permutation matrices $P$.

I have heard that there are several types of algorithms for these kinds of problems, like iterative improvement, simulated annealing, tabu search, genetic algorithms, evolution strategies, ant algorithms, and scatter search. I am looking for existing software.


问:
Computing linear combinations of sines and cosines (phasors)
I have a finite series that looks like this:

$f(t) = \sum^n_{i=0} A_i cos(\Theta_i + \omega_i t) + B_i sin(\Theta_i + \omega_i t)$

That is, a finite series of pairs of phasors.

What's the state of the art for numerically calculating the value of the series at a given $t$?  I'm interested in numerically finding the root of this series, so it's important to minimize numerical inaccuracies, particularly around the zeroes.

Just naively performing the cosine and sine calls produces a lot of numerical noise, which causes my root finding to fail.  I've built up a solution that works using Kahan summation and various trig identities (see this), but it's really, really slow.

The N I have is not very large, so an answer doesn't have to necessarily scale well to more terms.

If anyone knows any links or papers that are relevant I'd appreciate it.


问:
How to reduce RAM requirement in PETSC and SLEPCreading large binary matrices
I have two 50000 x 50000 binary matrices A and B for solving A*x=lambda*B*x eigenvalue problem. These matrices are sparse. I am trying to solve using PETSC and SLEPC. My memory requirement shoots off like more than 200 GB of RAM !

I used mpirun in 16 cores and 96GB RAM with swap of 130GB.

Is there a way to solve this problem.? Am I doing something wrong ?
Kindly let me know.


问:
How do I perform chebyshev interpolation from a to b with custom angle range?
Typically Chebyshev interpolation from $-1$ to $1$ with angle from $0$ to $\pi$:

$\xi_j=\cos \left ({\pi j \over N}\right )$

$x_j=(1+\xi_j) * {L \over 2}$

$w$:

$w_0=\pi/(2N)$

$w_{1,...,N-1}=\pi/(N)$

$w_N=\pi/(2N)$

I appear to have a need in the same with angle from $0$ to $\pi/2$ only, while solving a collocation grid problem where mesh is very fine inside of the computational domain and gets coarser near computational domain boundary. Literature review does not appear to be fruitful.


问:
Immersed boundary method
I'm trying to use immersed boundary method for the 3D flow problem (Navier-Stokes equations), but I'm maybe misunderstood something in this method. Main principles I took from this book. I use the continuous forcing approach and the quasi-rigid immersed boundary, i.e. the boundary force computed from this equation , where k is the stiffness, x ref is the initial node position.

My trouble is that in the first time layer I have zero boundary force (because I don't know this force) and nonzero fluid velocity. So the fluid shifts the all boundary nodes and I have something like this (this is the 3D flow around a cylinder that was projected on the Oz axis)

Because the all nodes were shifted they all have the same force including the right half circle. So the streamlines twisted through the immersed boundary (a cylinder).

How can I avoid this?

UPD

I have steady flow problem and maybe it's the source of troubles. I thought that the series of calculation of steady problem can approximates the unsteady and allow to use immersed boundary method, but now I'm not sure about this.

UPD
I have reached the more correct flow, when stiffness has been increased. I didn't understand this, because my calculation had not enough time steps.

Title: Immersed boundary method


问:
Is it possible to generalize the two view Sampson error to multiple view cases in computer vision?
In multiple view geometry of computer vision, there is a geometric error called Sampson error which is very useful in the nonlinear estimation of fundamental matrix.

Can this Sampson error be generalize to multiple view cases so as to use it for multiple focal tensor estimation?


问:
Calculating the number of molecules diffusing out of a volume
I have a system of reactions that are governed by differential equations.

They are reacting inside of a volume with known dimensions i.e lbh. I don't have any other information on their position inside that volume. I simply know there are 200 of A, 300 of B etc.

when molecules pass out of this volume they no longer react in the same way and move around in a Brownian fashion.

I want to know how often molecules leave this volume and I want to simulate this system.

I suppose I want a flux, or a probability distribution showing how likely an individual molecule leaves. I want to be able to "Count" the number of molecules left in V at the end of a timestep.

Ways I have thought about solving this:

Fick's first law:
I get a flux with this, but I thought you couldn't apply it in systems where the concentration isn't constant.

Fick's second law:
I don't really know how I would put it into my simulation

Smirnov density:

$$p(t)=\frac{x_0}{\sqrt{4\pi Dt^3}}\exp\left(-\frac{x_0^2}{4Dt}\right)$$

got this from https://physics.stackexchange.com/questions/93498/simulating-diffusion-from-bulk-to-individual-particles haven't really seen any other information on this. Maybe it's called something different?

With this equation there is a problem in that I don't know where the particles are in the volume.

The other approach I have thought of is using gibson-bruck next reaction and treating a transition as a reaction.

Please provide computational details in your answer, I need to be able to implement this.


问:
Conjugate residual/gradient convergence checking in practice
Let's say we want to solve $Ax=b$ ($A$ symmetric positive /semi/definite) with the conjugate residual/gradient method. $A$ comes from FEM where the mesh is being refined. The exact solution is $x_*$ and the approximate solution after $m$ iterations is $x_m$. Say that the initial guess is $x_0=0$.

We have that
$$\varepsilon=\left(\frac{(x_*-x_m, A(x_*-x_m))}{(x_*-x_m, A(x_*-x_m))}\right)^{1/2}=\left(\frac{(A^{1/2}(x_*-x_m), A^{1/2}(x_*-x_m))}{(A^{1/2}(x_*-x_m), A^{1/2}(x_*-x_m))}\right)^{1/2}\le2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^m$$
with $\kappa$ being the condition number of $A$. But how can I calculate or estimate this in practice? I don't have $x_*$, all I have is $x_m$ and thus $r_m=b-Ax_m$. Sure, I can calculate
\begin{align}
\tilde{\varepsilon}=\left(\frac{(r_m,r_m)}{(r_0,r_0)}\right)^{1/2}&=\left(\frac{(b-Ax_m,b-Ax_m)}{(b-Ax_0,b-Ax_0)}\right)^{1/2}=\left(\frac{(b-Ax_m,b-Ax_m)}{(b-Ax_0,b-Ax_0)}\right)^{1/2} \\
&=\left(\frac{(A(x_*-x_m),A(x_*-x_m))}{(A(x_*-x_0),A(x_*-x_0))}\right)^{1/2}=\left(\frac{(A(x_*-x_m),A(x_*-x_m))}{(b,b)}\right)^{1/2}
\end{align}
but then what? For a fix $A$ I think I'll still have the same asymptotic decrease in the $A$-norm when using $\tilde{\varepsilon}$ as stopping criterion, that is $\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^m$, so if I continue until $\tilde{\varepsilon}<\tilde{\varepsilon}_{\mathrm{max}}$ then the number of iterations needed are $m=c\sqrt{\kappa}$. But even so, $c$ depends on $A$! That means my iteration count could behave oddly if I'm using the algorithm for different matrices $A$ (corresponding to a mesh that has been refined a different number of times), does it not? Using $\tilde{\varepsilon}$ might not give the wanted $\sqrt{\kappa}$ growth of the iteration count. (Now, I'm not sure why $\varepsilon$ is a better stopping criterion than $\tilde{\varepsilon}$, but using one for theory and the other for implementation without understanding how they're connected feels very immoral.)


问:
Stability question (finite difference): dealing with corner nodes
Consider one initial boundary value problem for sphere.
$$\frac{\partial u}{\partial t}=\operatorname{div}A\nabla u +f$$
Here is explicit numerical scheme (we consider that it is stable):

$$\frac{u^{n+1}_{i,j,k}-u^{n}_{i,j,k}}{\Delta t}=A_{i,j,k}{\frac{\partial u}{\partial r}}_{i,j,k,+/-}^{n}+B_{i,j,k}{\frac{\partial u}{\partial \varphi}}_{i,j,k,+/-}^{n}+C_{i,j,k}{\frac{\partial u}{\partial \psi}}_{i,j,k,+/-}^{n}+
R_{i,j,k}{\frac{\partial^{2} u}{\partial r^2}}_{i,j,k}^{n}+\Phi_{i,j,k}{\frac{\partial^{2} u}{\partial \varphi^2}}_{i,j,k}^{n}+\Psi_{i,j,k}{\frac{\partial^{2} u}{\partial \psi^2}}_{i,j,k}^{n}+f^{n}_{i,j,k} ;$$
$u(r,\varphi,\psi,t=0)=u_{0}(r,\varphi,\psi)$;$u(r=r_{0},\varphi,\psi,t)=0;$

Under $+/-$ for operator it is understood that for operator was used upwind scheme; the operators are standart.
It is clear that such scheme works for all domain, except points on $r=r_{0}$(actually, we don't need to compute values here, because they are known because of boundary condition) and for point $r=0$ in the center of the ball.

Consider the following method for evaluating $u^{n+1}$ in the center of the ball: we use the definition of divergence

$$\operatorname{div}(A(x,y,z))=\lim\limits_{R\to 0}\frac {\oint_{B_{R}(x,y,z)}{(A,n)ds}}{V_{R}}, $$
therefore,we can write

$$\frac{u^{n+1}_{1,j,k}-u^{n}_{1,j,k}}{\Delta t}=\frac {\oint_{B_{r_2}}{(A\nabla u,\overline{n})ds}}{V_{r_2}} +\overline{f}^{n}_{1,j,k} ,$$
where $B_{r_2}$ -- a ball of radius equal to the value of unit step for $r$, $\overline{f}$ -- mean value of $f$ on the surface of $B_{r_2}$.

Question: will this method be stable for the center of the point?
If not, what stable method of computing the value in the center of ball you can offer?


问:
COMSOL - 2D box - get for every x an average value over y (line average apply for more lines)
I would like to ask you for your help.

I have 2D box. I know how to do an 2D cut line and compute an average value of someting over this line. (in Results -> Data Sets -> Cut Line 2D and then Results->Derived Values->Line Average and choose Cut Line 2D 1).

But I would need to know this average value for many of such parallel Cut Lines. I need to compute a curve of average value of something in y for every x. Have you any idea how to get it? I have tried to produce in Cut Line 2D the Additional parallel lines and then do average on it but the result was one number. I wanted a number for every cut line.

How can I do it? I would like to get an average value of pressure (in y) as a function of coordinate x.


问:
Approximate convolution of independent Beta variates?
Is there a way to approximate the convolution of Beta variables?

Specifically, I am trying to find an approximation to $g(x_0)$:

$$g(x_0) = \int \delta(x_0-\sum_{i=1}^{n} a_i x_i) \prod_{i=1}^{n} f(\alpha_i,\beta_i;x_i)\mathrm{d}x_i$$

where $\delta(x)$ is the Dirac delta function and

$$f(\alpha_i,\beta_i;x_i)=\frac{x_i^{\alpha_i - 1}(1-x_i)^{\beta_i - 1}}{\mathrm{B}(\alpha_i,\beta_i)}$$

for all $0<x<1$, and $f(\alpha_i,\beta_i;x_i)=0$ otherwise.

For example, it is known that the convolution of Gamma functions can be approximated by another Gamma function, whose mean and variance is the sum of the mean and variances of the convoluted Gamma distributions (1). Is there a similar result for the convolution of Beta distributed variables?

(1) Stewart, Trevor, et al. "A simple approximation to the convolution of gamma distributions." (2007).


问:
How to solve nonlinear optimization with constraints that have singular jacobian
I'm solving a nonlinear constrained optimization with constraint of following form:

$$\mathbf{A}^T\mathbf{A}-\mathbf{I}=\mathbf{0}, \mathbf{B}^T\mathbf{B}-\mathbf{I}=\mathbf{0}$$

where $\mathbf{A}$ and $\mathbf{B}$ are matrices that I'm solving for. Usually the initial guess of $\mathbf{A}$ and $\mathbf{B}$ are diagonal matrices with entries of 1 or -1.

I flatted the entries of $\mathbf{A}$ and $\mathbf{B}$ in column order into a vector as the variable $\mathbf{x}$ to solve for. I tried NLopt and Opt++ to solve my problem, both of them failed. While debugging I found that the gradient of constraint is a singular matrix. $\mathbf{A}$ and $\mathbf{B}$ are not necessarilly square matrices, thus neither is gradient of constraint.

My question is: is there any algorithm to solve this nonlinear constrained optimization, which does not require gradient of the constraint?

=======Edit=========

The objective function is:

$$\textrm{off}(\mathbf{A^T\Lambda_XA})+\textrm{off}(\mathbf{B^T\Lambda_YB})+\mu\left\|\mathbf{F^T\Phi{}A}-\mathbf{G^T\Psi{}B}\right\|_{\mathbf{F}}^2$$

where off() means sum of the squared off-diagonal elements, $\Lambda_X$ and $\Lambda_Y$ are diagonal matrices, $\mathbf{F}$ and $\mathbf{G}$, $\mathbf{\Phi}$ and $\mathbf{\Psi}$ are all known matrices.


问:
Minimizing a quadratic form
I would like to minimize the following quadratic form:
$$
f(\mathbf{\theta}) = (\mathbf{y} - \mathbf{\mu}(\mathbf{\theta}))^T \mathbf{\Sigma}({\mathbf{\theta}})^{-1}
(\mathbf{y} - \mathbf{\mu}(\mathbf{\theta})),
$$
wrt $\mathbf{\theta}$. $\mathbf{y}$ is $d$-dimensional, while $\mathbf{\theta}$ is $p$-dimensional. Generally $d \geq p$ and $\mathbf{\Sigma}_{\mathbf{\theta}}$ is positive definite for any $\mathbf{\theta}$.

The partial derivatives:
$$
\frac{\partial \mathbf{\mu}(\mathbf{\theta})_{i}}{\partial \mathbf{\theta}_j}
$$
are relatively easy to estimate. On the other hand the derivatives:
$$
\frac{\mathbf{\Sigma}({\mathbf{\theta}})_{ij}}{\partial \mathbf{\theta}_k}
$$
are nearly impossible to approximate.

What I'm doing at the moment is optimizing only along $\mathbf{\mu}(\mathbf{\theta})$ while treating $\mathbf{\Sigma}({\mathbf{\theta}})$ as if it was fixed. In other words I'm using steepest descent where at each step I move in the direction:
$$
\frac{\partial \mathbf{\mu}(\mathbf{\theta})}{\partial \mathbf{\theta}_j}^T\frac{\partial f(\mathbf{\theta})}{\partial \mathbf{\mu}(\mathbf{\theta})} =
\frac{\partial \mathbf{\mu}(\mathbf{\theta})}{\partial \mathbf{\theta}_j}^T\mathbf{\Sigma}({\mathbf{\theta}})^{-1}(\mathbf{y} - \mathbf{\mu}(\mathbf{\theta})).
$$

My question is: is there any guarantee doing this procedure converges to the minimum of $f(\mathbf{\theta})$? Thanks!

Additional Information

I have two possible situation:

Finite sample: it doesn't exist a $\mathbf{\theta}$ such that $\mathbf{y} = \mathbf{\mu}(\mathbf{\theta})$.

In the limit: there is a unique $\mathbf{\theta_0}$ such that $\mathbf{y} = \mathbf{\mu}(\mathbf{\theta_0})$.

Title: Minimizing a quadratic form


问:
When a flow changes its character from time-periodic to chaotic?
In this simple benchmark describing the flow around a disc in the plane with stationary bc, for low Reynolds numbers the solution is stationary and for higher RN it is periodic in time. Is there a theory saying, by which RN periodicity disappears and the solution(s) become(s) completely chaotic (if it happens at all)?</p>


问:
Transport Equation in a Tube: Source Term on Boundary
I'm modeling mass transport in a flow reactor. The flow reactor is a tube, which allows me to use cylindrical symmetry in solving the Convection-Diffusion-Reaction (CDR) Equation, which governs the mass transport.

The equation that I'm solving is:

$D ( \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} ) - v(r) \frac{\partial u}{\partial z} = S(r)$

Here, r and z are the coordinates in the respective directions, u is the property I want to solve for (the mass concentration), v is the velocity field, S is a given source and D is the diffusion coefficient.

The Boundary conditions are:

$\frac{\partial c}{\partial r} = 0, $ at $r=0$ and at $r=R$ (where R is the length of the rectangle in the r direction).

The system is sketched in the figure below:

My approach is to use the Method of Lines, i.e. to solve the system in the r direction for each step in the z direction.

I discretised the above equation as follows:

$\frac{du}{dz} = \frac{D}{v_i}(\frac{u_{i+1} - 2u_i + u{i-1}}{2h}+\frac{1}{r_i}\frac{u_{i+1}-u_{i-1}}{2h}) - S_i$

For both boundary conditions I used the central difference approximation, i.e. at r = 0 and r = R:

$\frac{u_{i+1} - u_{i-1}}{2h} = 0$

I am using a grid with two ghost points, as shown in the figure below:

Now, the code works ok, if I do not put the source on any of the boundaries, see the first plot below. However, in the actual experiment the source is in the centre of the flow reactor, i.e. the source is located on the boundary r = 0. When I put the source there, I get quite large oscillations (see the second plot below).

I think that I've either implemented the wrong boundary condition or that I'm using the wrong grid for this kind of computation.

My thoughts on removing those oscillations are:
 1. use a staggered grid, i.e. do not actually put a source on the boundary
 2. use a mixed boundary condition

What would you recommend? Am I missing something else?

The MatLab code is also attached below.

```<code>clear all
close all

%radius of flow reactor in mm
R = 13

%step-size in r direction
hr = 0.2
r = -hr : hr : R + hr

%parameters
D = 1.0 %diffusion
v = 1./(13.1 - r)
v = ones(1,length(r))

%build diagonals a, b, c of matrix A
%
%      |b c 0 0 0|
%      |a b c 0 0|
%  A = |0 a b c 0|
%      |0 0 a b c|
%      |0 0 0 a b|
%
% with a = alpha - beta, b = -2alpha, c = alpha + beta

alpha = ones(length(r)-2,1) ./ v(2:end-1)' .* (D/hr^2)
beta = ones(length(r)-2,1) ./ (v(2:end-1).* r(2:end-1))' .* (D/(2*hr))

diag_a = alpha - beta
diag_b = -2*alpha
diag_c = alpha + beta

A= spdiags([diag_a diag_b diag_c],[-1 0 1],length(r)-2,length(r)-2)

%implement boundary conditions du/dr = 0 at r = 0
A(1,1) = -4*alpha(1)
A(2,1) =  4*alpha(1)

%implement boundary conditions du/dr = 0 at r = R
A(end-1,end)= A(end-1,end)*2

ht = 0.00001
u = zeros(length(r)-2,1);

source = zeros(length(r)-2,1);
source(1:length(r)/2-2) = 1.

n_steps = 500

clear ustored

ustored = []

%step forward in time with explicit Euler
for i = 1:n_steps

    ustored = [ustored u];
    u_next = u(:,1) + ht*(A*u(:,1)+source(:,1));
    u = u_next;

end

ustored = [ustored u]

time = 0:ht:(n_steps)*ht;

mesh(time,r(2:end-1),ustored)
</code>```


问:
Two-chordless cycle extraction from a failed comparability graph recognition
I have implemented a comparability graph recognition algorithm from  M.C. Golumbic's "Algorithmic graph theory and perfect graphs" book. It is hinted in Fekete, Schepers, and van der Veen's "Exact Algorithm for Higher-Dimensional Orthogonal Packing" paper that the algorithm can be modified such that given the negative recognition result you can obtain a two-chordless cycle. But it is not explained how it is done and my knowledge on graph theory isn't anything special.

My current intuition is that if the recognition algorithm fails at level $k$ then a cycle must have been created where a transitive orientation is not possible. And this would mean the set of edges at the point of failure at level $k$ would induce a subgraph that contains the cycle I'm looking for, and then I just need to prune away the edges that have no involvment with the cycle created.

My questions: Is my intuition correct and if not, what does Fekete, Schepers, and van der Veen's "simple to modify" hint mean?

I'm doing this for efficiency, I want to obtain two results with one algorithmic search.

import networkx as nx

def comparability(G):
  classification = dict()
  k = 0
  def explore(i, j):
    # k and classification from outer function
    adjacent_I = set(G.neighbors(i))
    adjacent_J = set(G.neighbors(j))
    classtest = lambda e: e in classification and abs(classification[e]) != k
    for m in adjacent_I:
      if m not in adjacent_J or classtest((j, m)):
        if (i, m) not in classification:
          classification[(i, m)] = k
          classification[(m, i)] = -k
          if not explore(i, m):
            return False
        elif (i, m) in classification and classification[(i, m)] == -k:
          return False
    for m in adjacent_J:
      if m not in adjacent_I or classtest((i, m)):
        if (m, j) not in classification:
          classification[(m, j)] = k
          classification[(j, m)] = -k
          if not explore(m, j):
            return False
        elif (m, j) in classification and classification[(m, j)] == -k:
          return False
    return True
  for i, j in G.edges():
    if (i, j) not in classification:
      k += 1
      classification[(i, j)] = k
      classification[(j, i)] = -k
      if not explore(i, j):
        return False, classification
  return True, classification

G = nx.Graph()
G.add_nodes_from([1,2,3,4])
G.add_edges_from([(1,2),(1,3),(1,4),(3,4)])
print comparability(G) # True

G = nx.Graph()
G.add_nodes_from([1,2,3,4,5,6])
G.add_edges_from([(1,2),(2,3),(2,4),(3,4),(3,5),(4,6)])
print comparability(G) # False
```


问:
Finite volume method
I have question connected with finite volume method. Consider equation

$$\frac{\partial u}{\partial t}=\operatorname{div}A\nabla u +f . \quad x\in \overline{B}_{1} (0)\subset \mathbb{R}^3 -\text{unit ball with its boundary}$$ I solved IBVP for this equation with a help of finite volume method, for grid I took uniform spherical coordinates grid,.
(picture lower is related to topic)

I have a problem with finding flows for the faces of cell: how to approximate the value
$$\oint_{\partial V_{i,j,k}}(A\nabla u, n)dS $$
(where $A$ is a variable matrix)? When $A$ is an identity matrix everything is clear. Could you help me with this issue or to advise literature where this issue is represented?

Title: Finite volume method


问:
Auto labeling algorithm
I have a set of points (2D space), and for every point there's a label (like city names on a map).
I want to find a real-time algorithm that allows labels to avoid overlapping, moving them from their original position if necessary.

I've heard about simulated annealing algorithm, but I can't find a good source to learn how to do this. Do you have some idea where I can find bibliography and practical examples?

Thanks in advance for your replies.

Title: Auto labeling algorithm


问:
How to do something in Modelica every solver iteration?
How is it possible to do something once every solver iteration?

The expressions when change(time) or when time > pre(time) are not permitted.

Best regards,
Koldo

Edit: The final goal is to do a cosimulation from Modelica with other tool that does not have FMI and where the Modelica side signals the iteration delta time.
- Approach 1: Is there any way to fire a 'when' under an algorithm section just before or after an iteration step has been successfully finished? - Approach 2: I wonder if there is something similar to fmiEventUpdate(). It appears in the automatically generated .mo file when doing FMI, and "is called after a time, state or step event occurred ... the function returns once a new consistent state has been found".


问:
Are the eigenvalues of the product matrix of two real symmetric square matrices also real values?
Suppose $A,B \in \mathbb{R}^{n\times n}; A=A^T, B=B^T$, let $C = AB, D =BA$,

If we have all the real eigenvalues of $A$ and $B$, e.g. the eigenvalue decomposition of them:

$A=P\Lambda_1 P^T$,
$B=Q\Lambda_2 Q^T$; both $P$ and $Q$ are orthogonal matrices, and both $\Lambda_1$ and $ \Lambda_2$ are diagnoal matrices;

can we also obtain the eigenvalues of $C$ and $D$ based on these information without solving eigenvalue problems of $C$ and $D$?

Are the eigenvalues of $C$ and $D$ also real?


问:
Matrix completion when the eigenvectors are a tensor product?
Suppose we have incomplete observations of the square matrix $X$. Most matrix completion algorithms assume the matrix is low-rank. What if instead we assume the matrix of eigenvectors is a tensor product: $V = A\otimes B$, $X=VDV^{-1}$? This still dramatically reduces the space of possible matrices, so it may be possible to exactly recover a matrix of this form. How would we solve this algorithmically?</p>


问:
Integrating highly oscillatory functions
I have a logarithmic grid, upon which i have two functions that are similar to this one (this is only the last 100 points):

These are essentially very similar to a Sin function at this point.  I need to take the integral of two of these functions multiplied together, along with the position:

$$
\int_0^R f_1(r)\; r\; f_2(r)\; dr
$$

Unfortunately, this gives me functions similar to the following:

I need to integrate these with high accuracy, but I'm not sure how to do it.  There are a few options that I know of, but I expect there are more that I'm not sure of:

I can interpolate the original functions (is cubic interpolation a good interpolation strategy for sinusoids and products of sinusoids?), then multiply the interpolated functions together at specific points and use something from this scicomp post

I can switch the order of the interpolation and multiplication, multiplying the values I have, then taking the interpolation.

I can do a simpler Simpsons rule approximation, though I'm not sure how to do such a thing on a non-uniform grid.


问:
Sequential approach to solving coupled PDEs
I'm dealing with a coupled system of three transient, non-linear convection-diffusion equations. Let's just say to simplify the problem that they take the following form:
$$
-\nabla\cdot(D_{1}(u_{2},u_{3})\nabla u_{1}) = -\nabla\cdot\mathbf{\chi}(u_{2},u_{3})
$$
$$
\frac{\partial u_{2}}{\partial t} + \nabla\cdot f_{1}(u_{1},u_{2},u_{3}) - \nabla\cdot(D_{2}(u_{2},u_{3})u_{2}) = 0
$$
$$
\frac{\partial u_{3}}{\partial t} + \nabla\cdot f_{2}(u_{1},u_{2},u_{3}) - \nabla\cdot(D_{3}(u_{2},u_{3})u_{3}) = 0
$$

Assume a general set of boundary conditions for each equation:
$$
u_{i} = g_{i,D}, \;\;\;\text{on}\;\Gamma_{D_{i}}
$$
$$
D_{i}\nabla u_{i} = g_{i,N}, \;\;\;\text{on}\;\Gamma_{N_{i}}
$$

My spatial discretization is via DGFEM and the time discretization is just Backward Euler. Anyway, I'll use the superscript $(k)$ to denote the solution at the $k$th time step. This is how I am decoupling the equations:

1) Solve for $u_{1}^{(k+1)}$ using information from the previous time step for $u_{2}$ and $u_{3}$. This effectively linearizes the equation.

$$
-\nabla\cdot(D_{1}(u_{2}^{(k)},u_{3}^{(k)})\nabla u_{1}^{(k+1)}) = -\nabla\cdot\chi(u_{2}^{(k)},u_{3}^{(k)})
$$

2) Use the new $u_{1}^{(k+1)}$ and information from the previous time step for $u_{3}$ to solve for $u_{2}^{(k+1)}$:

$$
\frac{\partial u_{2}}{\partial t} + \nabla\cdot f_{1}(u_{1}^{(k+1)},u_{2}^{(k+1)},u_{3}^{(k)}) - \nabla\cdot(D_{2}(u_{2}^{(k+1)},u_{3}^{(k)})\nabla u_{2}^{(k+1)}) = 0
$$

3) Solve for $u_{3}^{(k+1)}$:

$$
\frac{\partial u_{3}}{\partial t} + \nabla\cdot f_{2}(u_{1}^{(k+1)},u_{2}^{(k+1)},u_{3}^{(k+1)}) - \nabla\cdot(D_{3}(u_{2}^{(k+1)},u_{3}^{(k+1)})\nabla u_{3}^{(k+1)}) = 0
$$

So the first equation is linear and the last two are non-linear. I am using Newton's method to handle the non-linearity. However, there are some bugs in my code and I'm a bit stuck trying to interpret the results.

I know this is long but I'm trying to be as detailed as possible here. Just a couple of questions:

1) I test the simple case where $u_{1} = u_{2} = u_{3} = 5+t$ and I get about machine precision error using linear basis functions. Newton's method is also converging quadratically. I'm using a timestep approximately equal to $h^{p+1}$. I can also increase it to $100h^{p+1}$ which surprisingly doesn't seem to affect the convergence much. Using this time step, I have:

```<code>/* sample Newton iterations */
******** Current Time: 7.812500 ********
u2 Iteration: 1, Error: 4.155186e-01
u2 Iteration: 2, Error: 2.480219e-02
u2 Iteration: 3, Error: 9.169256e-05
u2 Iteration: 4, Error: 1.250692e-09
u2 Iteration: 5, Error: 8.042476e-16
Iterations for u2 convergence: 5

u3 Iteration: 1, Error: 4.155186e-01
u3 Iteration: 2, Error: 2.480219e-02
u3 Iteration: 3, Error: 9.169256e-05
u3 Iteration: 4, Error: 1.250692e-09
u3 Iteration: 5, Error: 7.664044e-16
Iterations for u3 convergence: 5

/* error results */
*********** L2 Error ***********
u1 error in L2-norm:
4.53162e-12
4.64396e-12
6.15316e-12
1.07398e-11

u2 error in L2-norm:
4.12483e-13
7.30247e-13
1.39962e-12
2.76142e-12

u3 error in L2-norm:
4.09102e-13
7.29944e-13
1.39935e-12
2.76121e-12

*********** H1 Error ***********
u1 error in H1-norm:
1.05844e-11
2.03320e-11
3.84458e-11
7.47941e-11

u2 error in H1-norm:
6.40033e-13
1.33164e-12
2.71254e-12
5.47814e-12

u3 error in H1-norm:
6.34670e-13
1.33089e-12
2.71211e-12
5.47862e-12
</code>```

Now the bizarre part is that if I take $u_{1}=5+t$, $u_{2}=6+t$, and $u_{3}=7+t$, my results are awful. I'm thinking maybe the first case was the exception rather than the norm since it's the only example I've tested so far that seems to work. Anyway, the weird part about this one is that I'm getting about machine precision error for $u_{2}$ while the others are flatlining around 1.e0. I'm not sure how this could possible be since the solution to $u_{2}$ depends on the solutions to $u_{1}$ and $u_{3}$, which are incorrect if my code is to be believed. I observe quadratic convergence for $u_{2}$, but obviously not the others...

The time step is around $h^{p+1}$ here, but I observe similar results for smaller time steps.

```<code>******** Current Time: 0.074219 ********
u2 Iteration: 1, Error: 3.908065e-03
u2 Iteration: 2, Error: 2.219552e-06
u2 Iteration: 3, Error: 2.889154e-12
Iterations for u2 convergence: 3

u3 Iteration: 1, Error: 1.373826e-01
u3 Iteration: 2, Error: 1.564705e-02
u3 Iteration: 3, Error: 2.492157e-03
u3 Iteration: 4, Error: 7.872758e-04
u3 Iteration: 5, Error: 1.311116e-04
u3 Iteration: 6, Error: 2.997904e-05
u3 Iteration: 7, Error: 9.158730e-06
u3 Iteration: 8, Error: 1.519301e-06
u3 Iteration: 9, Error: 3.555316e-07
u3 Iteration: 10, Error: 1.067623e-07
u3 Iteration: 11, Error: 1.755777e-08
u3 Iteration: 12, Error: 4.195674e-09
u3 Iteration: 13, Error: 1.242190e-09
u3 Iteration: 14, Error: 2.024999e-10
u3 Iteration: 15, Error: 4.936768e-11
u3 Iteration: 16, Error: 1.443518e-11
u3 Iteration: 17, Error: 2.333785e-12
Iterations for u3 convergence: 17

*********** L2 Error ***********
u1 error in L2-norm:
3.02951e+00
5.43305e+00
5.96871e+00
4.25056e+00

u2 error in L2-norm:
1.45749e-13
2.71039e-13
3.69646e-13
4.95591e-13

u3 error in L2-norm:
6.29385e-01
8.30180e-01
1.11972e+00
1.13907e+00

*********** H1 Error ***********
u1 error in H1-norm:
7.86736e+00
1.81562e+01
3.67674e+01
3.97038e+01

u2 error in H1-norm:
2.33954e-13
2.00202e-12
4.95659e-12
1.03835e-11

u3 error in H1-norm:
2.00438e+00
3.91391e+00
6.24097e+00
7.06979e+00
</code>```

I think maybe the most telling part is the convergence of the Newton iterations. Namely, what does it mean if I'm not observing quadratic convergence? My initial guess is just the solution at the previous time step, and since the time step is quite small, I'm pretty sure the initial guess is not the issue. I know that non-quadratic convergence can imply that the Jacobian is not being computed properly, but the first example is so similar to the second that I can't imagine that's the case...anyone have any insights here?

2) I'm a bit iffy on the spatial discretizaton of the convection term:
$$
\int_{\Omega}\nabla\cdot f_{1}(u_{1},u_{2},u_{3})z = -\int_{\Omega}\nabla z\cdot f_{1}(u_{1},u_{2},u_{3}) + \int_{\Gamma}z(f_{1}\cdot\mathbf{n})
$$

(via Green's identity). The DG discretization is
$$
-\sum_{E\in\mathcal{M}_{h}}\int_{E}\nabla z\cdot f_{1} + \sum_{e\in\Gamma_{\text{int}}}\int_{e}\left\{\left\{f_{1}\cdot\mathbf{n}\right\}\right\}[[z]] + \sum_{e\in\Gamma_{\text{out}}}\int_{e}(f_{1}\cdot\mathbf{n})z + \sum_{e\in\Gamma_{\text{in}}}\int_{e}(f_{1}\cdot\mathbf{n})z
$$

where {{ }} is the average and [[ ]] is the jump. In all of the literature I've come across, the inflow integral is moved to the right-hand side. So if we're considering the second equation (solving for $u_{2}$), is this correct?

$$
\sum_{e\in\Gamma_{\text{in}}}\int_{e}(f_{1}(u_{1},g_{2,D},u_{3})\cdot\mathbf{n})z
$$

So when computing the Jacobian for the second equation, does this inflow term not contribute to the Jacobian at all since $u_{2,D}$ is known?

I'm able to solve a single transient, non-linear convection diffusion equation using Newton's method, so I'm not entirely sure what it is that I'm missing.


问:
Efficient way to do congruent transformation using matrix inverse?
I know a square self-adjoint matrix $S_{vv}$ and I want to find:

$S_{rr} = HS_{vv}H^{\dagger}$

where $\dagger$ denotes conjugate transpose.

I do not know $H$ but I do know $H^{-1}$.

What is the most efficient way to solve something of this form?

I have considered rewriting it as:

$S_{rr} = H(HS_{vv})^{\dagger}$

Finding the LU decomposition of $H^{-1}$, and then doing two lots of forward and back substitutions to solve for $S_{rr}$.  But it feels like there might be a better way?

Any help much appreciated.


问:
Methods for integrating black box functions on a non-uniform grid
If i have some function expressed as points on a non-uniform grid (I'm specifically interested in logarithmic grids, but general results are also interesting), and I want to integrate it, I believe there are two general possibilities:

Interpolate, and use some method like Gauss Quadrature, or Clenshaw-Curtis Quadrature.  The downside is that you now have stacked errors from the interpolation and the integration method.  The upside is that these quadrature methods are known to be very accurate.

Use some integration method like the trapezoidal rule, or Simpson's rule, which use the points that you have available to you.

Is there a general rule of thumb for which method is better?  Does it depend on some basic property of the function (for example: periodic)?


问:
Is it possible to construct such a symmetric matrix with desired eigenvalues?
Suppose a real, dense and asymmetric square matrix $A\in\mathbb{R}^{n\times n}$,
all its eigenvalues $\lambda_i \in \mathbb R$

Is it possible to construct a symmetric matrix $B\in\mathbb{R}^{n\times n}$, the eigenvalues of which are $\lambda_i^2$ to compute absolute values of $\lambda_i$ without the need of solving eigenvalues of $A$?


问:
Fast Eigenvalue and SVD Solver for Structured Matrices
I am looking for a fast Eigenvalue and SVD solver for small dense structured matrices (Hankel and Toeplitz). I have searched for efficient implementations in libraries like MKL but I am not able to find anything specific to structured matrices. I found a set of interesting papers on the subject:

Luk and Qiao, "A Fast Eigenvalue algorithm for Hankel matrices" which uses Lancoz tridiagonalization based on FFT and a QR like diagonalization of the tridiagonal matrix for computing eigenvalues with a complexity of $O(n^2\log(n))$.

Luk and Qiao, "A fast symmetric SVD algorithm for square Hankel matrices"

Do any of the libraries out there or any implementation you know of use this as I would like to reuse as much as possible? Are there any algorithms which achieve bounds better than this for structured matrices?


问:
Markowitz Pivoting to reduce size of a dense integer system
I am dealing with a large sparse integer matrix that I need to find the nullspace of. I've seen Markowitz Pivoting come up in several places discussing similar problems such as here:

http://www.mapleprimes.com/posts/41191-Solving-Sparse-Linear-Systems-In-Maple

And in the MAGMA documentation under the sparse null space algorithm. How can one use this technique (or similarly perhaps, the Structured Gaussian Elimination mentioned above it) to reduce a sparse matrix to a dense, small matrix? I understand that you can eliminate variables with regular Gaussian Elimination, but how does one retrieve those after solving the new system?


问:
PETSc: KSP gives strange discontinuity at the edge between MPI processes
I'm using PETSc's SNES linesearch (NEWTONLS) to solve a system of nonlinear PDEs. SNES+KSP successfully solved a test problem of the same type before without any trouble.

However, with my current problem, if I use the KSP to invert the Jacobian, the intermediate solution develops a "ridge" at the boundary between two MPI processes. The image below shows the solution after one SNES step running with -ksp_gmres_modifiedgramschmidt -pc_type asm -sub_pc_type lu -sub_pc_factor_shift_type NONZERO.

The problem seems to go away if I execute the exact same binary, but using the MUMPS direct solver with -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps, as shown below.

In both cases above SNES decided to start by taking a full Newton step. However, in the KSP run the ridge continues to grow and the step size gets significantly suppressed over the next few steps before failing completely. In the MUMPS run, SNES converged straightforwardly after taking 5 full steps.

Does anyone have any idea what is happening here?


问:
Numerical integral with a weakly singular kernel with a satisfactory precision
I am working on a numerical method for time fractional PDE. One problem is that I must compute a numerical integral of the following form:

$$
\begin{equation}
\int_0^{t_m} (t_m-s)^{-\beta}f(s)ds
\end{equation}
$$

where $t_m$ and $0<\beta<1$ are constants. Noting that the integral has a weakly singular kernel, so classical quadrature, such as Gaussian quadrature, will invalid for the singularity of $(t_m-s)^{-\beta}$ at $s=t_m$. So my question is how to compute this kind of integral numerically with a given precision (because I need a high precision when computing the convergence ratio)?

My try:

I have deduced the Gaussian quadrature with the kernel of $(t_m-s)^{-\beta}$ which have a highest algebraic accuracy. But the formula is too complex so that when I compute the coefficients and nodes MATLAB loss precision which result in a bad precision.

Another method I have tried is to compute this by using integral by parts. Thus, the kernel will be $(t_m-s)^{1-\beta}/(1-\beta)$, which is without singularity. This method is valid when $f(s)$ have a good regularity, say in $C^1$. But when $f(s)$ does not have a continuous derivative, such as $|s|$ or Heaviside function, this method is also invalid.

What is more, $(t_m-s)^{-\beta}$ is only one kind of weakly singular kernel. How about other situations? Is there any general method to conquer this problem?


问:
Locally conservative method for differential generalized eigenvalue problem
I have to approximate the smallest eigenvalue of the following generalized eigenvalue problem

$$
- \nabla \cdot  D(x) \nabla p(x) + \alpha(x) p(x) = \lambda \beta(x) p(x)
$$

over a domain like follows (many more subdomains, this is only an approximation to illustrate the problem)

The eigenvalues are all positive, so to approximate the smallest one I invert the left side of the equation after the spatial discretization, and use Arnoldi method. The functions $D(x)$, $\alpha(x)$ and $\beta(x)$ are piecewise constant over the whole domain, being probably different over different subdomains.

I already have a code to approximate the smallest eigenvalue and associated eigenvector with a Continuous Galerkin FEM.

I was wondering which locally conservative method can I use to solve the eigenvalue problem (i.e., I want the solution to satisfy the balance equation over each subdomain).

What I know:

OPTION 1: I know that using discontinuous Galerkin FEM the balance equation will be satisfied over each subdomain, but it usually deal with the mixed formulation of the problem

$$
- \nabla \cdot  q(x) + \alpha(x) p(x) = \lambda \beta(x) p(x) \\
q(x) - D(x) \nabla p(x) = 0
$$

and it introduces a lot of degrees of freedom to the system.

OPTION 2: I have read that it can also be done using a Hybridizable discontinuous Galerkin, because the we can apply a static condensation technique and the degrees of freedom related to $q$ are the only ones that we have to solve in a coupled way. It works fine for source problems, but I think that we can not apply this technique with an eigenvalue problem.

edit: (why static condensation will not work for eigenproblems?)
If dealing with a source problem, we form the linear system for the hybridizable discontinuous Galerkin as follows:

$$
\left[\begin{array}{c c}
A & B \\
C & D
\end{array}\right]
\left[\begin{array}{c}
p \\
q
\end{array}\right]
=
\left[\begin{array}{c}
r \\
s
\end{array}\right]
$$

We know that matrix $A$ is block diagonal because the degrees of freedom are completelly decoupled (from element to element), so we perform the static condensation to solve the following problem

$$
(D-CA^{-1}B)q = s-CA^{-1} r
$$

and then we recover $p$ with

$$
Ap = r-Bq
$$

When we form the algebraic eigenvalue problem with discontinuous galerking we get

$$
\left[\begin{array}{c c}
A & B \\
C & D
\end{array}\right]
\left[\begin{array}{c}
p \\
q
\end{array}\right]
=
\lambda \left[\begin{array}{c c}
E & 0 \\
0 & 0
\end{array}\right]
\left[\begin{array}{c}
p \\
q
\end{array}\right]
$$

Then, as in a linear system, we try to isolate $p$ in the first equation, and we obtain

$$
p = (A-\lambda E)^{-1}B q
$$

and after substitution in the second equation we get

$$
(D + C(A-\lambda E)^{-1}B ) q = 0
$$

Matrices $A$ and $E$ are block diagonal, so easy to invert, but the eigenvalue $\lambda$ now is inside the matrix to invert, so the eigenvalue problem has been converted into a nonlinear eigenvalue problem. Different linear approximations of  $(A-\lambda E)^{-1}$ (taylor expansion on $\lambda = 0$ or $\lambda = \lambda_0$, for instance) will solve the problem of the nonlinear eigenvalue problem, but then the algorithm will converge much more slower and succesive corrections will be necessaries. Thus static condensation is not a good idea here.

OPTION 3: There are some papers that talk about locally conservative continuos Galerkin methods: The Continuous Galerkin Method Is Locally Conservative, and Locally Conservative Fluxes for the Continuous Galerkin Method. The first one seems to perform some extra calculations to approximate the fluxes, that for my is equivalent to the mixed formulation. And the second one says that is posible to deal with the primal formulation, but also develop the article about a formulation that includes the fluxes $q(x)$.

QUESTION: it is no clear for me that locally conservation can be obtained with a primal formultaion (like in the first equation I posted), moreover taking into accout that we are solving an eigenvalue problem. Can it be done?

Any idea or advice about literature on the topic would be greatly appreciated. I am mainly interested in FEM methods, but any idea is welcome.


问:
Dominant contributions of a quadratic form
Let $\Sigma$ be a covariance matrix (e.g. symmetric positive definite). For arbitrary vectors $\epsilon$, I need to compute $\chi^2 \equiv \epsilon^\top\Sigma^{-1}\epsilon$, which I do using a Cholesky decomposition.

If $\Sigma=I$, then $\chi^2$ would be the 2-norm of $\epsilon$. I would then say that each coefficient of $\epsilon$ had the same weight in $\chi^2$. Similarly, if $\Sigma$ is a diagonal matrix with positive diagonal elements $(\lambda_1,\ldots,\lambda_N)$, $\epsilon_i$ would be weighted by $\lambda_i^{-1}$ :

$$ \chi^2 = \sum_{i=1}^N \frac{\epsilon_i^2}{\lambda_i} $$

It is then trivial to conclude that, in general (i.e. irrespective of the value of $\epsilon$), the coefficients of $\epsilon$ that contribute most to $\chi^2$ are those where $\lambda_i$ is biggest.

Unfortunately, $\Sigma$ is not diagonal. We can of course diagonalize it: $PDP^\top \equiv \Sigma$ ($D$ is diagonal and $P$ is the matrix whose columns are the normalized eigenvectors of $\Sigma$; $P$ is unique if we order the eigenvalues from largest to smallest in $D$). Diagonalization allows us to isolate independent contributions to $\chi^2$, and we get

$$\chi^2 = \sum_{i=1}^N \frac{(p_i^\top\epsilon)^2}{\lambda_i} \qquad \text{where} \quad \Sigma p_i \equiv \lambda_i p_i$$

What previously was a scalar weight now becomes a vector weight, and is $p_i/\sqrt\lambda_i$

I am only interested in dominant contributions to $\chi^2$, i.e. where $\lambda_i$ is small and $p_{ij}$ is large.

Analytically, what is a proper way to determine dominant contributions to $\chi^2 \equiv \epsilon^\top\Sigma^{-1}\epsilon$? Is it possible to tell that some region of $\epsilon$ will be weighted more or less than another one, accounting for the covariances?

Numerically, the $1/\lambda$ dependence of the discussion suggests that the smaller the eigenvalue, the larger the contribution. That is very unfortunate: my $\Sigma$ matrix is very ill-conditioned, and computation of the smallest eigenvalues is probably just gibberish. Is there a way to do implement this discussion either robustly or by avoiding completely the calculation of eigenvalues?


问:
Numerical Quadrature of Oscillating Integral With Non-oscillating part
As you will know there are different numerical integrals (I believe Levin's method is the most popular one) for the numerical quadrature of oscillating integrands which may roughly speaking be written as a product of oscillating part and non-oscillating amplitude $exp(i \omega t) A(t)$ (the form can of course be slightly more general). Integrands of the form $exp(i \omega t) A(t)+B(t)$ may be treated using a method for oscillatory and one for non-oscillatory for the first and second part, respectively. I wonder if there are routines also for the efficient quadrature of integrands which cannot be analytically decomposed as $exp(i \omega t) A(t)+B(t)$ but which nevertheless feature a rapidly oscillatory behavior (you can assume fixed frequency) . I.e. they behave like $exp(i \omega t) A(t)+B(t)$ with some unknown $A$ and $B$. It appears to me that I could in principle do better than using a conventional quadrature for non-oscillating integrands as I could try to find $A$ and $B$ using smooth (and successively refining) fits and then run a quadrature for oscillating integrals in order to obtain a numerical estimate for the integral. I wonder if there are any known efficient methods for such problems.  </p>


问:
Algorithms to compute largest gap between smallest nonzero eigenvalues of sparse symmetric matrix
I am looking mainly for c/c++ implementations but also for theoretical algorithms to compute gaps between smallest positive eigenvalues of symmetric, singular matrix or real numbers. To be precise, I want to find index n of the eigenvalue such that between $\lambda_n$ and $\lambda_{n+1}$ (both are nonzero) there is the largest gap (I assume that eigenvalues are sorted). The matrix I am dealing with is huge, so I am allowed to use only sparse representation. I cannot store all the zero entries due to memory constraints.
Can you recommend any c/c++ software which can compute it to me and uses sparse matrix representations? I have tried to install fortran ARPACK, but with no success.
Maybe you know some papers with algorithms solving this problem. Thanks in advance!</p>


问:
Shift-Invert in Anasazi/Belos using Tpetra Sparse Matrices
I am currently trying to port my code over to Trilinos because the problems that I am working on are too big for LAPACK/ARPACK. Specifically I am computing the generalized eigenvalues/eigenvectors:

$$\mathbf{A}x = \lambda\mathbf{B}x$$

where both $\mathbf{A}$ and $\mathbf{B}$ are large and sparse. I am mainly interested in targeting specific interior eigenvalues using a shift-invert method:

$$\mathbf{A}x - \sigma\mathbf{B}x = \lambda\mathbf{B}x - \sigma\mathbf{B}x$$
$$\mathbf{A}x - \sigma\mathbf{B}x = (\lambda - \sigma)\mathbf{B}x$$
$$\frac{1}{\lambda - \sigma}x = (\mathbf{A}x - \sigma\mathbf{B})^{-1}\mathbf{B}x$$

This maps eigenvalues that are close to $\sigma$ towards $\infty$ so that they are easy to find by iterative methods. To avoid computing the inverse in the previous expression we instead solve the linear system at every iteration:

$$w = \mathbf{B}v_i$$
$$(\mathbf{A}x - \sigma\mathbf{B})v_{i+1} = w$$

Previously I was solving this system using LU decomposition but as the matrices become large, this method is unwieldy. Instead I am trying to use Belos to solve this linear system for the Anasazi eigensolver package. An important point to note is that my matrices are complex-valued and so Epetra is a less than ideal solution, I would rather get Tpetra working.

I've been looking at this Epetra example which defines an Anasazi::EpetraGenOp generalized operator to solve the linear system. This generalized operator basically computes $\mathbf{A}^{-1}\mathbf{B}v_i$ for a shift of $\sigma = 0$.

Anasazi::EpetraGenOp is defined in the Epetra adapter for Anasazi but the Tpetra adapter does not seem to have an equivalent generalized operator.

My question:

What is the correct way to solve a shift-inverted generalized eigenvalue problem using Tpetra/Belos/Anasazi?


问:
Derivation of a Higher Order Compact Alternating Direction Implicit Method
I dont understand how this Higher Order Compact ADI scheme, which is fourth order in time and space, for the wave equation is derived:
I go through the following

Using Tylor's expansion $u(t+h,x,y)$ and $u(t-h,x,y)$ about the point (t,x,y) we have:\
(1)
\begin{equation}
\label{311}
u(t+k,x,y)-2u(t,x,y)+u(t-k,x,y)=k^2\left(\frac{\partial^2 u(t,x,y)}{\partial t^2}\right)+\frac{k^4}{12}\left(\frac{\partial^4 u(t,x,y)}{\partial t^4}\right)+\cdots
\end{equation}
If $u$ is a solution of $u_{tt}=u_{xx}+u_{yy}$ then we have the following:\
(2)

\begin{equation}
\label{312}
\frac{\partial^2 }{\partial t^2}=\frac{\partial^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}
\end{equation}

 And \
 (3)
\begin{equation}
\label{313}
\frac{\partial^4 }{\partial t^4}=\frac{\partial^4 }{\partial x^4}+\frac{\partial^4 }{\partial y^4}
\end{equation} (3)
Substituting equations (2) and (3)  in to equation (1) and $u(t+k,x,y)=u^{n+1}_{i,j}$,$u(t-k,x,y)=u^{n-1}_{i,j}$ and $u(t,x,y)=u^{n}_{i,j}$ shows:\
(4)
\begin{equation}
\label{314}
u^{n+1}_{i,j}-2u^{n}_{i,j}+u^{n-1}_{i,j}=k^2\left(\frac{\partial^2u(t,x,y) }{\partial x^2}+\frac{\partial^2 u(t,x,y)}{\partial y^2}\right)+\frac{k^4}{12}\left(\frac{\partial^4 u(t,x,y)}{\partial x^4}+\frac{\partial^4 u(t,x,y)}{\partial y^4}\right)+\cdots
\end{equation}
Again using Taylor's expansion up to the fourth order difference for the following:\
(5)
\begin{equation}
\label{315}
\frac{\partial^2 u}{\partial x^2}=\frac{1}{h^2}\left(\delta^2_x-\frac{1}{12}\delta_x^4\right)u
\end{equation}
And \
 (6)
\begin{equation}
\label{316}
\frac{\partial^2 u}{\partial y^2}=\frac{1}{h^2}\left(\delta^2_y-\frac{1}{12}\delta_y^4\right)u
\end{equation}

And we also have :
\(7)
\begin{equation}
\label{317}
\frac{\partial^4 u}{\partial x^4}=\frac{\partial^2}{\partial x^2}\left( \frac{\partial^2 u}{\partial x^2}\right)=\frac{1}{h^4}\delta^4_xu
\end{equation}
And \(8)
\begin{equation}
\label{318}
\frac{\partial^4 u}{\partial y^4}=\frac{\partial^2}{\partial y^2}\left( \frac{\partial^2 u}{\partial y^2}\right)=\frac{1}{h^4}\delta^4_yu
\end{equation}
Substituting equations (5)-(8) in to equation (4) shows to the fourth order accuracy:\(9)
\begin{equation}
\label{319}
u^{n+1}_{i,j}-2u^{n}_{i,j}+u^{n-1}_{i,j}=\left(r^2\left( 1+\frac{r^2-1}{12}\delta_x^2\right)\delta^2_x +r^2\left( 1+\frac{r^2-1}{12}\delta_y^2\right)\delta^2_y\right) u_{i,j}^n
\end{equation}
The next step may be applying some operator,but I dont know how to apply and to reach the following ?

$$(1+\frac{1-r^2}{12}\delta_x^2)u^*_{ij}=[2+\frac{1+5r^2}{6}\delta_x^2-\frac{2(1+r^2)}{1-r^2}\delta_y^2]u_{ij}^n-[1+\frac{1-r^2}{12}\delta_x^2]u_{ij}^{n-1}$$
$$(1+\frac{1-r^2}{12}\delta_y^2)u^{n+1}_{ij}=u^*_{ij}+\frac{1-r^2}{12}\delta_y^2[\frac{2(1+10r^2+r^4)}{(1-r^2)^2}u_{ij}^n-u_{ij}^{n-1}]$$

where $h_x=h_y=h$


问:
Solving ODE boundary problem with additional conditions
I need to solve two point bondary problem for ODE. The solution of the problem itself is very easy. The real issue is that when solving arising non linear system of equations it always converges to one solution. However there exist two (physically justified) solutions. To illustrate this problem in detail consider the following case:

The discrete form of ODE is:
$$
y_{i+1}=y_i+\frac{\Delta x}{2} (a\cdot y_i^2+a\cdot y_{i+1}^2)
$$
I want to find $a$ and $y_i$ values. As one can see in each node $x_i$ there are two possible solutions of this equation. Sometimes I want to pick the ones that are placed on the left of the minimum of the parabol, sometimes those on the right. It depends on the value of some additional function:

$$f(x_{i+1},y_{i+1})$$

that is computed in every step.

How do I formulate a problem? Should I look into optimization techniques?


问:
Tracking for two meshes
I'm dealing with physical simulation (position based dynamic). Now I'm trying to realize tracking for two meshes. In order to explain what does it mean, let's assume that we have two similar meshes: one with fine resolution and other with coarse. A need to make in the way that mesh with fine resolution followed the coarse mesh moving. But it should not be exact following, instead it should follow only in general. Let's consider this task in terms of mathematics:
 First we need to build Laplace-Beltrami discrete operator. I use next definition
$$L=V^{-1/2}CV^{-1/2},$$
where $V$ is a diagonal matrix of Voronoi areas, $C$ is a sparse symmetric matrix of cotangent weights.
After that I perform an eigenvalue decomposition
$$L=QAQ^{T}$$
The next step is to build the matrix of harmonic functions $T$;

On Internet I have found that it could be built in the next way
$$T=(V^{-1/2}Q)^{T}$$
So now I have $n*n$ matrix T, where n - is a number of vertices's in the fine mesh.
Now, if have T matrix I can build the next constraint $$T*p_f==TBp_c$$

$p_f$ is a vector of vertices's in the fine mesh

$p_c$ is a vector of vertices in the coarse mesh

$B$ such matrix that $Bp_c=p_f$.

If this constraint is satisfied fine mesh exactly follows for the coarse mesh.
Now, if we discard from matrix $T$ those eigenvectors, which correspond to the biggest eigenvalues we obtain the $k*n$ matrix $T$, and it means that we discarded information about high frequent harmonics in the mesh spectral representation. As result mesh with fine resolution will follow for the coarse only in general but such small details like wrinkle will be different.

In my implementation it works only if I use full $T$ matrix, but when I discard at least one eigenvector, behavior becomes unpredictable.

Question:
Is there any fundamental errors in this approach? I'm confident in all, except of definition for $T$ matrix.

Title: Tracking for two meshes


问:
Spatio-temopral wavelet analysis
Am quite new to wavelet analysis and would like some help. I am performing a spatio-temporal analysis of monthly gridded rainfall data. With PCA, I can reduce the dimension of the rainfall data into a few leading modes, yielding EOFs (which explain spatial variability) and principal components (explaining temporal variability).

I would like to perform a similar analysis with wavelets using the Matlab Wavelet Toolbox. As of now, I am able to decompose a 2D data (spatial decomposition in my case) but unable to take into account the temporal variability in the data.

My first course of action has been to first compress the data with PCA and then perform wavelet decomposition of the leading modes in both the spatial (EOFs) and temporal (PCs) domain.

I am wondering if this is the right way to perform such an analysis and would like suggestions as to how to proceed.

Thanks a lot.

Title: Spatio-temopral wavelet analysis


问:
BFGS Fails to Converge
The model I'm working on is a multinomial logit choice model. It's a very specific dataset so other existing MNLogit libraries don't fit with my data.

So basically, it's a very complex function which takes 11 parameters and returns a loglikelihood value. Then I need to find the optimal parameter values that can minimize the loglikelihood using scipy.optimize.minimize.

Here are the problems that I encounter with different methods:

'Nelder-Mead’: it works well, and always give me the correct answer. However, it's EXTREMELY slow. For another function with a more complicated setup, it takes 15 hours to get to the optimal point. At the same time, the same function takes only 1 hour on Matlab using fminunc (which uses BFGS by default)

‘BFGS’: This is the method used by Matlab. It works well for any simply functions. However, for the function that I have, it always fails to converge and returns 'Desired error not necessarily achieved due to precision loss.’. I've spent lots of time playing around with the options but still failed to work.

'Powell': It quickly converges successfully but returns a wrong answer. The code is printed below (x0 is the correct answer, Nelder-Mead works for whatever initial value), and you can get the data here: https://www.dropbox.com/s/aap2dhor5jyxy94/data.csv

Thanks!

```<code>import pandas as pd
import numpy as np
from scipy.optimize import minimize

# https://www.dropbox.com/s/aap2dhor5jyxy94/data.csv
df = pd.read_csv('data.csv', index_col=0)
dfhh = df.hh
B = df.ix[:,'b0':'b4'].values # NT*5
P = df.ix[:,'p1':'p4'].values # NT*4
F = df.ix[:,'f1':'f4'].values # NT*4
SDV = df.ix[:,'lagb1':'lagb4'].values

def Li(x):
    b1 = x[0] # coeff on prices
    b2 = x[1] # coeff on features
    a = x[2:7] # take first 4 values as alpha
    E = np.exp(a + b1*P + b2*F) # (1*4) + (NT*4) + (NT*4) build matrix (NT*J) for each exp()
    E = np.insert(E, 0, 1, axis=1) # (NT*5)
    denom = E.sum(1)
    return -np.log((B * E).sum(1) / denom).sum()

x0 = np.array([-32.31028223, 0.23965953, 0.84739154, 0.25418215,-3.38757007,-0.38036966])
np.random.seed(0)
x0 = x0 + np.random.rand(6)

 minL = minimize(Li, x0, method='Nelder-Mead',options={'xtol': 1e-8, 'disp': True})
# minL = minimize(Li, x0, method='BFGS')
# minL = minimize(Li, x0, method='Powell', options={'xtol': 1e-12, 'ftol': 1e-12})
print minL
</code>```

Update: 03/07/14 Simpler Version of the Code
Now Powell works well with very small tolerance, however the speed of Powell is slower than Nelder-Mead in this case. BFGS still fails to work.

Title: BFGS Fails to Converge


问:
Self-consistent numerical solution of a set of equations
I am trying to solve an assignment on solving the Bogoliubov de Gennes equations self-consistently in Matlab. BdG equations in 1-Dimension are as follows:-

$$\left(\begin{array}{cc}
-\frac{\hbar^{2}}{2m}\frac{\delta^{2}}{\delta z^{2}}-\mu+V\left(z\right) & \triangle(z)\\
\triangle(z) & \frac{\hbar^{2}}{2m}\frac{\delta^{2}}{\delta z^{2}}+\mu-V(z)
\end{array}\right)\left(\begin{array}{c}
u_{n}(z)\\
v_{n}(z)
\end{array}\right)= \epsilon_{n}\left(\begin{array}{c}
u_{n}(z)\\
v_{n}(z)
\end{array}\right)$$
along with the equations for gap function  $\triangle(z)$ and number density $n(z)$.
   $$\triangle(z)=U\sum_{n}\left(1-2f_{n,}\right)u_{n}(z)v_{n}^{\star}(z)$$
and
$$n(z)=2\sum_{n}|{u_{n}(z)}|^{2}f_{n}+|{v_{n}(z)}|^{2}\left(1-f_{n}\right).$$

For the case of solving the BdG equations in Fourier space in Matlab for the case of a periodic potential and periodic gap function (assumed), we can take $$u_{n}(z)=\sum_{k}\exp\left[ikz\right]U_{n,k}, $$
$$\triangle(z)=\sum_{K}\exp(iKz)T_{K},$$
and
$$ V(z)=\sum_{K}\exp(iKz)P_{K} $$
where the sum is over the reciprocal lattice vectors $K$ leading to the number equation
$$N=2\sum_{n,k}\left[f_{n}|{U_{n,k}}|^{2}+\left(1-f_{n}\right)|{V_{n,k}}|^{2}\right]$$
with $ f_{n}$ as the Fermi distribution function $f_{n}=\frac{1}{\exp(\beta(\epsilon-\mu))+1}$. Solving the set of equations self-consistently for a fixed $N$, I am trying to get a value of chemical potential from the number equation each time after solving the eigenvector components $U_{n,k}$ and $V_{n,k}$, but due to the form of the exponentials in the number equation and sum over large number of them, I am unable to get a correct value of chemical potential out of them using Matlab routines as the root of the equation to put it back into the equations for eigenvector components.

In most cases, I get random values of chemical potential since the equation is more or less insoluble.  How can I avoid this error ? Is there a better way to numerically solve the BdG equations self-consistently ? I also want to do this assignment in real space avoiding finite size effects but started with the Fourier space case to avoid errors associated with discretizing the differential. Please guide and ask for any details you might need.

Following is my MATLAB code to solve the equations in real space but the code does not work as fsolve does not find the mu value.

http://postimg.org/image/v7amx5vd9/full/

http://postimg.org/image/do9pyyk7z/full/


问:
Rank one constraint in optimization, which is nonconvex, but will we find the global maximum?
The objective function is $L = D(B)$, here $D$ is a function of $B$ and $D$ is a convex function. We have a prior information of $B$: it is rank one and we can write $B$ into $B = pq^T$. So we can change the objective function into $L=D(p, q)$. Obviously, for $p$ and $q$, this problem is not a convex problem anymore and we can have many sets of $p$ and $q$. My question is that even though we have a lot of local minimum regarding $p$ and $q$, will we have a global minimum $B$ even though there are many $p$ and $q$? Thanks so much for your answering.   </p>


问:
Modal analysis of structure with aerodynamic damping
I'm using modal decomposition to predict the steady state response of a beam structure to harmonic loading. The structure itself is very lightly damped, but we know from experiments that the aerodynamic drag force as it moves through the air is very significant.

I have the mass and stiffness matrices for the beam ($\boldsymbol{[M]}$ and $\boldsymbol{[K]}$) and so I am able to obtain the eigenvectors and eigenvalues ($\boldsymbol{[V]}$ and $\boldsymbol{[D]}$). The eigenvalues yields a vector containing the natural frequencies of the system (the square root of the diagonal entries of $\boldsymbol{[D]}$, I'll call it $\boldsymbol{\omega}$.

The system is excited by a force vector in the global reference frame $\boldsymbol{F_0}$ at an excitation frequency $\omega_F$.

There is a rotation matrix which rotates the beams from the local to the global reference frame $\boldsymbol{[R]}$.

I have successfully used the procedure to calculate the steady state response without aerodynamic damping:

Get eigenvalues and eigenvectors $\boldsymbol{[V]}$ and $\boldsymbol{[D]}$

Get generalised force vector $\boldsymbol{Q}=\boldsymbol{[V]}^T \boldsymbol{F_0}$

Calculate modal participation factor vector $\boldsymbol{q}$. The individual elements are calculated as follows:

$$q_i=\frac{Q_i}{{\omega_i}^2}\frac{1}{\sqrt{[(1-(\omega_F/\omega_i)^2)^2+(2\zeta_i\omega_F/\omega_i)^2]}}$$

Calculate response from $\boldsymbol{x}=\boldsymbol{[V]q}$

Of course, this drastically overpredicts the amplitude of the oscillations because there is no aerodynamic damping. I've found a few examples for single d.o.f. systems in which the energy lost from the system is calculated (by integrating the drag force x the velocity). The integration results in the following equation for energy loss:
$$\Delta E=\frac{8}{3}\frac{1}{2}\rho A C_D \omega^2 X^3$$
This can then be converted to an equivalent damping coefficient as follows:
$$C_{eq}=\frac{\Delta E}{\pi \omega X^2}$$
Which can be converted to a damping coefficient by doing:
$$\zeta=\frac{C_{eq}}{2 \omega_n}$$
I've modelled a simple spring-mass system with air resistance and it works fine (I have to iterate to find the final position because the amplitude calculation requires knowledge of the damping which requires knowledge of the amplitude so there is a circular reference, but that's no biggy!)

Could someone please explain the equivalent for a multi-d.o.f. system? My initial forays have been unsuccessful, with much higher drag coefficients than would really occur.


问:
What type of boundary condition(s) do I need to define for a process with diffusion of tracer concentration?
I want to simulate the diffusion of initial dye concentration in a small estuary with a 2D depth averaged model. Would a radiation boundary condition for concentration of tracers alone be enough or do I need to define any additional conditions at the open boundary? </p>


问:
Discretizing boundary conditions for vortex methods
I am working on a fluid simulation using vortex methods. For this I must compute the vortex sheet on my boundaries given as:
$$
\gamma(\mathbf{x}) - \frac{1}{\pi}\int_S \frac{\partial}{\partial\mathbf{n}}(\text{log}|\mathbf{x} - \mathbf{x'}|)\gamma(\mathbf{x'})d\mathbf{x'} = -2\mathbf{u}\cdot \mathbf{s}
$$
I have two questions on how to compute the vortex sheet $\gamma$:

The paper of Cooper and Barba state this can be solved using radial basis functions.
I think the hardest part to compute is
$$
\oint \frac{\partial}{\partial n}[\text{log}|\mathbf{x}-\mathbf{x'}|]\phi(|\mathbf{x'} - \mathbf{x_i}) d\mathbf{x'}
$$
What would the discrete version of this be? So how do I compute this term?

If I would use the panel method, how does a discretisation look like?


问:
Turbulence: Suggest research papers similar to this experiment
I am trying to find the effect on turbulence characteristics when flow passes through a converging section. For this experiment I will use a PIV system to gather data about the flow.
However, I could not find any specific literature on flow through symmetric converging channel. (The converging channel I would be using is type of a converging hopper with a circular hole exit.) I have attached a 2d drawing of this

Please advise if you know any similar papers relating to this.


问:
BLACS context value and multiple MPI communicators
I am reposting here a question previously asked in stackoverflow: as suggested in a reply to my question, scicomp could be a better place to obtain some useful comments/suggestions for my problem.

I am trying to make some tests with BLACS/Scalapack (C interface, Intel MKL version) making use multiple MPI communicators: particularly, what I am trying to obtain is a set of BLACS contexts/grids corresponding (one-to-one) to a set of (disjoint) MPI communicators. While I have no problems in working with a single "global" communicator (MPI_COMM_WORLD) I have some difficulties with the multiple communicator case. I hope that you could give me some suggestions for the following problem.

I am a bit confused about the BLACS context variable value update after a call to Cblacs_gridinit: let's suppose to start with a "global" context corresponding to the MPI_COMM_WORLD communicator. I can obtain the corresponding BLACS context with the call:

```<code>MPI_Comm globalCommunicator(MPI_COMM_WORLD);

MKL_INT globalContext(Csys2blacs_handle(globalCommunicator));
</code>```

and create a grid on it with

```<code>Cblacs_gridinit(&globalContext,
                &c_blacsGridOrdering,
                MKL_INT i_nTaskRow,
                MKL_INT i_nTaskCol);
</code>```

with, for example, char c_blacsGridOrdering('R').

The globalContext value in this case is 0.

At some point in my code, which is making use of 6 MPI tasks,
an MPI communicator (localCommunicator) corresponding to the group of 4 tasks with id [0;3] is being created: at this point I would like to create a new BLACS context (localContext) for this "local" communicator and a local grid on it. I can do that with the code

```<code>MKL_INT localContext(Csys2blacs_handle(localCommunicator));

Cblacs_gridinit(&localContext,
                &c_blacsGridOrdering,
                MKL_INT i_nTaskRowLocal,
                MKL_INT i_nTaskColLocal);
</code>```

where the above function call sequence is performed only by the tasks included in the local communicator.

The localContext value, after the call to Csys2blacs_handle, is equal to 1 (for each task of the local communicator) but it is modified and set to 0 by the subsequent call to Cblacs_gridinit.

Obviously this causes some problems in the rest of my code since, for example, if I try to retrieve the MPI communicator corresponding to localContext with a call to Cblacs2sys_handle I get a six task communicator corresponding to the initial globalCommunicator.

Most likely I am making some very stupid error (missing call to some BLACS functions??) or it is not clear to me the interplay between multiple MPI communicators and BLACS context/grids but I can not find what is wrong with my code.

Do you have some suggestions concerning the above problem? Many thanks for your help!

/////////
UPDATE 1
/////////

I have an update for my question which could be useful to find a solution or at least an explanation for the observed problem: the described problem is no more present if I initialize the first grid (the one making use of the global context) on a task grid whose size is such that all available MPI tasks are included, for example with MKL_INT i_nTaskRow(1) and MKL_INT i_nTaskCol(MPI_COMM_WORLD size). Is this behavior expected from BLACS? Many thanks again for your support!