Try   HackMD

Solving Linear Algebra Problem

Linear equation

Linear equations try to solve following kind of equation
\begin{equation} \textbf{A}\textbf{x} = \textbf{b} \tag{1} \end{equation}
where \(\textbf{A}\) is a NxN matrix, \(\textbf{x}\) is a Nx1 and \(\textbf{b}\) is a Nx1 column. There are two ways to solve it by using linear package, LAPACK. First is using the subroutine \(dgesv\) which solve \(\textbf{x}\) immediately. Second method is that we can find the inverse matrix of \(\textbf{A}\) , then multiplied by \(\textbf{b}\). Finding inverse matrix in LAPACK require two subroutine, \(dgetrf\) and \(dgetri\).

Example 1

In the first example we first solve \(\textbf{x}\) and then verified the result by multiply \(\textbf{x}\) by the matrix \(\textbf{A}\). I then solve this equation for \(10^6\) times to compare its run time to another method. Finally I computer the error of solution by calculate the difference between \(\textbf{A}\textbf{x}\) and
\(\textbf{b}\) of two method.

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

As you can see both run time and accuracy are superior when using subroutine \(dgesv\).

Example 2

However, when the matrix \(\textbf{A}\) is a symmetric matrix the later method has a higher accuracy despite a slower speed.

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Code





































































program main
    ! Example program that solves the matrix
    ! equation A * x = B using LAPACK.
    implicit none

    double precision :: a(5, 5)  ! Matrix A.
    double precision :: b(5)     ! Matrix B.
    ! Pivot indices (list of swap operations).
    double precision :: pivot(5) 
    integer :: rc       ! Return code.
    double precision :: aa(5, 5), bb(5), aaa(5, 5)
    double precision :: work2(5)
    double precision :: start, finish
    integer :: i

    a = reshape([4.d0 , 1.d0 , 1.d0 , 0.d0 , 1.d0,  &
                 -1.d0, 3.d0 , 1.d0 , 1.d0 , 0.d0,  &
                 2.d0 , 1.d0 , 5.d0 , -1.d0, -1.d0, &
                 -1.d0, -1.d0, -1.d0, 4.d0 , 0.d0,  &
                 0.d0 , 2.d0 , -1.d0, 1.d0 , 4.d0], &
                 [5, 5])
    b = [6.d0, 6.d0, 6.d0, 6.d0, 6.d0]

    ! method one
    call cpu_time(start)
    do i = 1, 1000000
        aa = a
        bb = b
        call dgesv(5, 1, aa, 5, pivot, bb, 5, rc)
    end do
    call cpu_time(finish)

    if (rc == 0) then
        write(*, '(a)') repeat('=', 50)
        print *, 'Using dgesv:'
        print '(a, 5(f0.4, ", "))', 'x = ', bb
    else
        print '(a, i0)', 'Error: ', rc
    end if
    print '(a, 5(f0.4, ", "))', 'Ax = ', matmul(a, bb)
    write(*, *) 'cost(run 10^6 times)', finish-start, 's'
    write(*, *) 'accuracy: ', sum(matmul(a, bb)-b)

    ! method two
    write(*, '(a)') repeat('=', 50)
    print *, 'Using dgetrf+dgetri:'

    call cpu_time(start)
    do i = 1, 1000000
        aaa = a
        call dgetrf(5, 5, aaa, 5, pivot, rc)
        call dgetri(5, aaa, 5, pivot, work2, 5, rc)
    end do
    call cpu_time(finish)

    if (rc == 0) then
        print '(a, 5(f0.4, ", "))', 'x = ', matmul(aaa, b)
        print '(a, 5(f0.4, ", "))', 'Ax = ', matmul(a, matmul(aaa, b))
    else
        print '(a, i0)', 'Error: ', rc
    end if
    write(*, *) 'cost(run 10^6 times)', finish-start, 's'
    write(*, *) 'accuracy: ', sum(matmul(a, matmul(aaa, b))-b)
    write(*, '(a)') repeat('=', 50)
end program main

Eigen-system

Using \(dgeev\) we can solve the eigenvalue and eigenvector of the matrix \(\textbf{A}\). Following print the result of five eigenvalue and five eigenvector. By definition the relation between eigenvalue and eigenvector is given by
\begin{equation} \textbf{A} * v(j) = \lambda(j) * v(j) \tag{2} \end{equation}
So I use this relation to test the accuracy of my result.

Example 1

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Example 2

Same procedure for finding another eigenvalue and eigenvector of second matrix.

Wigner's semi-circle and Girko's Circular Laws in random matrix

Real-Symmetric NxN matrix

I create a symmetric matrix with the size of \(5000 \times 5000\) that fill with uniform distribution in \([-0.5, 0.5]\). Then I solve the eigenvalue of this system, since this matrix is symmetric, all the eigenvalue are real. I plot the probability distribution of eigenvalue, the distribution when the size of the array become infinite is given by
\begin{equation} P(\lambda) = \dfrac{2}{\pi R^2}\sqrt{R^2-x^2} \tag{3} \end{equation}

Code










































































program main
    ! find the eigenvalue and eigenvector of asymmetric matrix A
    ! equation A * x = B using LAPACK.
    implicit none

    integer, parameter :: n=5
    double precision :: a(n, n) 
    integer :: rc, lwork
    double precision :: aa(n, n)
    double precision :: wr(n), wi(n), work(1000), vl(n, n), vr(n, n)
    integer :: i, j
    double complex :: lambda(n), v(n, n)

    a = reshape([4.d0 , 1.d0 , 1.d0 , 0.d0 , 1.d0,  &
                 -1.d0, 3.d0 , 1.d0 , 1.d0 , 0.d0,  &
                 2.d0 , 1.d0 , 5.d0 , -1.d0, -1.d0, &
                 -1.d0, -1.d0, -1.d0, 4.d0 , 0.d0,  &
                 0.d0 , 2.d0 , -1.d0, 1.d0 , 4.d0], &
                 [5, 5])
    aa = a

    lwork = -1
    CALL DGEEV('N', 'V', n, aa, n, wr, wi, vl, n, &
                vr, n, work, lwork, rc )
    lwork = MIN( 1000, INT( work( 1 ) ) )
    CALL DGEEV('N', 'V', n, aa, n, wr, wi, vl, n, &
                vr, n, work, lwork, rc )

    ! write eigenvalue
    write(*, *) 'Eigenvalue'
    do i = 1, n
        write(*, *) i, wr(i), wi(i)
    end do
    write(*, *) repeat('=', 70)

    lambda = 0.d0
    do i = 1, n
        lambda(i) = complex(wr(i), wi(i))
    end do

    ! write eigenvector
    write(*, *) 'Eigenvector'
    do i = 1, n
        j = 1
        do while( j.le.n )
            if( wi(j).eq.0.d0) then
                write(*, 9998, advance='no') vr( i, j )
                v(i, j) = complex(vr(i, j), 0.d0)
                j = j + 1
            else
                write(*, 9999, advance='no') vr( i, j ), vr( i, j+1 )
                write(*,9999,advance='no') vr( i, j ), -vr( i, j+1 )
                v(i, j) = complex(vr(i, j), vr(i, j+1))
                v(i, j+1) = complex(vr(i, j), -vr(i, j+1))
                j = j + 2
            end if
        end do
        write(*,*)
    end do
    write(*, *) repeat('=', 70)

    write(*, *) 'check: A * v(j) - lambda(j) * v(j) ?= 0'
    do j = 1, 5
        print 9999, sum(matmul(a, v(:, j)) - lambda(j) * v(:, j))
    end do
    write(*, *) repeat('=', 70)

9998 FORMAT( 11(:,1X,F6.2) )
9999 FORMAT( 11(:,1X,'(',F6.2,',',F6.2,')') )
end program main

Girko’s circle law

If we use more general asymmetric N-by-N matrix the eigenvalue would be a complex number. I can the plot distribution of real part and imaginary part in the complex plane. I then scale the data by \(\sqrt{12/n}\), so it seats on a unit circle. As you can see by increasing the size of the matrix, the dot are filling more in the unit circle.



Linear boundary-value problem ODE

The last question are trying to solve a boundary-value ODE, using Finite-Difference method we can easily transform this ODE into a problem of linear algebra in which I can use previous result to solve it. I divide the interval into different number, and we can see that the method converge to the exact solution very fast. It only require few point to converges to the exact solution to this ODE.








































































































































































program main
    ! find the eigenvalue and eigenvector of asymmetric matrix A
    ! equation A * x = B using LAPACK.
    implicit none

    integer, parameter :: n=5000
    integer, parameter :: lmax = 100000
    double precision :: a(n, n) 
    integer :: rc, lwork
    double precision :: aa(n, n)
    double precision :: wr(n), wi(n), work(lmax), vl(n, n), vr(n, n)
    integer :: i, j
    double complex :: lambda(n)
    double precision :: rand

    call random_seed()
    do i = 1, n
        do j = 1, i
            call random_number(rand)
            ! symmetric matrix
            a(i, j) = rand - 0.5d0
            a(j, i) = rand - 0.5d0
        end do
    end do
    aa = a

    lwork = -1
    CALL DGEEV('N', 'N', n, aa, n, wr, wi, vl, n, &
                vr, n, work, lwork, rc )
    lwork = MIN( lmax, INT( work( 1 ) ) )
    CALL DGEEV('N', 'N', n, aa, n, wr, wi, vl, n, &
                vr, n, work, lwork, rc )
!
    ! write eigenvalue
    ! write(*, *) 'Eigenvalue'
    ! do i = 1, n
        ! write(*, *) i, wr(i), wi(i)
    ! end do
    ! write(*, *) repeat('=', 70)

    open(66, file='a.dat', status='unknown')
    lambda = 0.d0
    do i = 1, n
        lambda(i) = complex(wr(i), wi(i))
        write(66, *) wr(i), wi(i)
    end do

    ! write eigenvector
    ! write(*, *) 'Eigenvector'
    ! do i = 1, n
        ! j = 1
        ! do while( j.le.n )
            ! if( wi(j).eq.0.d0) then
                ! write(*, 9998, advance='no') vr( i, j )
                ! v(i, j) = complex(vr(i, j), 0.d0)
                ! j = j + 1
            ! else
                ! write(*, 9999, advance='no') vr( i, j ), vr( i, j+1 )
                ! write(*,9999,advance='no') vr( i, j ), -vr( i, j+1 )
                ! v(i, j) = complex(vr(i, j), vr(i, j+1))
                ! v(i, j+1) = complex(vr(i, j), -vr(i, j+1))
                ! j = j + 2
            ! end if
        ! end do
        ! write(*,*)
    ! end do
    ! write(*, *) repeat('=', 70)
!
    ! write(*, *) 'check: A * v(j) - lambda(j) * v(j) ?= 0'
    ! do j = 1, 5
        ! print 9999, sum(matmul(a, v(:, j)) - lambda(j) * v(:, j))
    ! end do
    ! write(*, *) repeat('=', 70)

! 9998 FORMAT( 11(:,1X,F6.2) )
! 9999 FORMAT( 11(:,1X,'(',F6.2,',',F6.2,')') )
end program main
\end{minted}
\subsection*{4.f95}
\begin{minted}[frame=lines, linenos, fontsize=\large]
{fortran}
program FD
    ! solve y'' = p(x)y' + q(x)y + r(x)
    ! solve Aw = b
    implicit none
    integer, parameter :: N=101
    double precision, parameter :: pi=4.d0*atan(1.d0)
    double precision, external :: p, q, r
    double precision :: A(N, N), w(N), b(N)
    double precision :: xi, xf, h, wi, wf
    integer :: i
    double precision :: pivot(N), work2(N)
    integer :: rc
    ! if N is 1
    ! double precision :: A(N, N), w(N), b(N)
    
    ! bd condition
    wi = -.3d0
    wf = -.1d0
    xi = 0.d0
    xf = pi / 2.d0
    h = (xf - xi) / (dble(N) + 1.d0)

    ! A matrix
    ! first row
    A(1, 1) = 2.d0 + h**2 * q(xi+h)
    A(1, 2) = -1.d0 + .5d0 * h * p(xi+h)
    do i = 2, N-1
        A(i, i-1) = -1.d0 - .5d0 * h * p(xi+dble(i)*h)
        A(i, i) = 2.d0 + h**2 * q(xi+dble(i)*h)
        A(i, i+1) = -1.d0 + .5d0 * h * p(xi+dble(i)*h)
    end do
    ! last row
    A(N, N-1) = -1.d0 - .5d0 * h * p(xi+dble(N)*h)
    A(N, N) = 2.d0 + h**2 * q(xi+dble(N)*h)

    ! b matrix
    b(1) = -h**2 * r(xi+h) + (1.d0 + .5d0 * h * p(xi+h)) * wi
    do i = 2, N-1
        b(i) = -h**2 * r(xi+dble(i)*h)
    end do
    b(N) = -h**2 * r(xi+dble(N)*h) + (1.d0 + .5d0 * h * p(xi+dble(N)*h)) * wf

    ! find the inverse of A
    call dgetrf(N, N, A, N, pivot, rc)
    call dgetri(N, A, N, pivot, work2, N, rc)

    ! w = A^{-1}b
    w = matmul(A, b)

    ! case: if N = 1
    ! A(1, 1) = 2.d0 + h**2 * q(xi+h)
    ! b(1) = -h**2 * r(xi+h) + (1.d0 + .5d0 * h * p(xi+h)) * wi + (1.d0 - .5d0 * h * p(xi+h)) * wf
    ! w(1) = b(1) / A(1, 1)



    ! output
    open(66, file='d.dat', status='unknown')
    write(66, *) wi
    do i = 1, N
        write(66, *) w(i)
    end do
    write(66, *) wf
    close(66)
end program FD

function p(x)
    implicit none
    double precision :: p, x
    p = 1.d0
end function

function q(x)
    implicit none
    double precision :: q, x
    q = 2.d0
end function

function r(x)
    implicit none
    double precision :: r, x
    r = dcos(x)
end function
Last changed by 
 
Morris Huang
My Personal Website: https://physics-morris.github.io/
0
70

Read more

Ferromagnetic Ising model on a square lattice

Ising model on square lattice are L \times L spin that can take 2 values, either S_i=+1 or S_i=-1. Each spin interact with neighbors and with external field. The Hamiltonian of the system is\begin{equation} H = -\sum_{<ij>} J_{ij}S_iS_j - \sum_i H_iS_i \end{equation}A particular solution for no external field are obtained by Lars Onsager where the second-order phase transition is happened at\dfrac{J}{k_BT_c} = \dfrac{ln(1+\sqrt{2})}{2}. With this in mind, we can use Metropolis algorithm to simulate Ising model with good accuracy. Following four figures show different MC steps, the spin of a 100 \times 100 system with blue represent spin up and red represent spin down. Since the temperature happened around the critical temperature, the magnetization taken very long time to come to equilibrium.

Mar 16, 2023
Non-reversal Random walk (NRRW)

Non-reversal random walk are randomly walk except it cannot reverse it's previous step. Following two figures show different steps of NRRW, compare to previous figure it can walk further since it cannot go reverse its previous step.

Mar 16, 2023
Ideal Random Walk

An ideal 1 dimension random walk with probability of p stepping to right and 1-p to the left are represented by below figure of numbers of particle of first 100 steps. From this we can clearly guess(see) that the probability P(x, N) at N step at location x are higher around the center and lower at the edge, just like a Gaussian distribution.

Mar 16, 2023
Gaussian random variable and CLT

To generate two independent Gaussian random variables x and y we use these equation\begin{equation} x = \sigma \sqrt{-2ln\xi_1}\cos(2\pi \xi_2) \quad,\quad y = \sigma \sqrt{-2ln\xi_1}\sin(2\pi \xi_2) \end{equation}This figure show total number N=10^5 Gaussian random numbers.

Mar 15, 2023
Read more from Morris Huang
Published on HackMD