"MathJax" と "prism" のテスト

次の4つの偏微分方程式をまとめて「マクスウェル方程式」と呼びます． \begin{align} \bm{\nabla} \cdot \bm{E}(\bm{r}, t) = \cfrac{\rho(\bm{r}, t)}{\varepsilon_0} \end{align} \begin{align} \bm{\nabla} \cdot \bm{B}(\bm{r}, t) = 0 \end{align} \begin{align} \bm{\nabla} \times \bm{H}(\bm{r}, t) = \bm{i}_e(\bm{r}, t) + \cfrac{\partial \bm{D}(\bm{r}, t)}{\partial t} \end{align} \begin{align} \bm{\nabla} \times \bm{E}(\bm{r}, t) = - \cfrac{\partial \bm{B}(\bm{r}, t)}{\partial t} \end{align}

 

 

for(int i=0; i

 
行列の固有値を「QR法」によって求めるプログラムのC言語ソース・コードを以下に示します．
void QR(int n, double mat_A[n][n], double A[n][n])
{
    //=====================================================
    //initialize
    double** Q = (double**)malloc(sizeof(double*) * n);
    double* Q_row = (double*)malloc(sizeof(double) * n * n);
    for(int i=0; i<n; i++)
        Q[i] = Q_row + n * i;
    
    double epsilon = pow(10, -8);
    int m = n;
    double* row_temp = (double*)malloc(sizeof(double) * n);
    double origin, r, s, c, temp;

    //=====================================================
    //Householder transformation
    Householder(n, mat_A, A);

    //=====================================================
    //QR decomposition
    while(m>1)
    {
        //-------------------------------------------------
        //convergence check
        if( fabs(A[m-1][m-2]) < epsilon )
        {
            m -= 1;
            continue;
        }

        //-------------------------------------------------
        //origin shift
        origin = A[n-1][n-1];
        if( m==n )
            origin = 0;
        for(int i=0; i<m; i++)
            A[i][i] -= origin;

        //-------------------------------------------------
        //initialize "Q"
        for(int i=0; i<m; i++)
        {
            for(int j=0; j<m; j++)     
                Q[i][j] = 0;
            Q[i][i] = 1;
        }

        //-------------------------------------------------
        //calculate "R" and "Q"
        for(int i=0; i<m-1; i++) //"i" is row index: 0 to m-2
        {
            //---------------------------------------------
            //calculate "sin" and "cos" to delete "A[i+1][i]"
            r = sqrt( A[i][i]*A[i][i] + A[i+1][i]*A[i+1][i] );
            if(r==0)
            {
                s = 0;
                c = 0;
            }
            else
            {
                s = A[i+1][i] / r;
                c = A[i][i] / r;
            }

            //---------------------------------------------
            //update "R" (= "A")
            A[i][i] = r;
            A[i+1][i] = 0;
            for(int j=i+1; j<m; j++)
            {
                temp = A[i][j] * c + A[i+1][j] * s;
                A[i+1][j] = -A[i][j] * s + A[i+1][j] * c;
                A[i][j] = temp;
            }

            //---------------------------------------------
            //update "Q"
            //A = (P_n-1 P_n-2 ... P_1)^-1 R = QR
            //Q^T = P_n-1 P_n-2 ... P_1
            for(int j=0; j<m; j++)
            {
                temp = Q[j][i] * c + Q[j][i+1] * s;
                Q[j][i+1] = -Q[j][i] * s + Q[j][i+1] * c;
                Q[j][i] = temp;
            }
        }

        //-------------------------------------------------
        //calculate "RQ" (RQ = AQ, save to "A")
        for(int i=0; i<m; i++) //"i" is row index
        {
            for(int j=i; j<m; j++) //"j" is column index
                row_temp[j] = A[i][j]; //the i-th row of "A"
            
            for(int j=0; j<m; j++) //"j" is column index
            {
                temp = 0;
                for(int k=i; k<m; k++) //"k" is product loop index
                    temp += row_temp[k] * Q[k][j];
                A[i][j] = temp;
            }
        }

        //-------------------------------------------------
        //post-processing of origin shift
        for(int i=0; i<m; i++)
            A[i][i] += origin;
    }

    free(Q_row);
    free(Q);
    return;
}