View difference between Paste ID: bLmP57mW and 4ZXMYJKJ

-subroutine LU(n, m, A, P, Y)
+module global_vars
-    integer, intent(in) :: n, m
+    integer :: N = 100, l = 2
-    integer :: i, j, k, pivot, temp
+    real :: R = 10.0, V_0 = 2, a = 5.0, lk = 0.1, rk = 10.0, dk = 0.05, k
-    double precision :: pivotValue
+    save
-    integer, dimension(:), intent(inout) :: P
+end module
-    double precision, dimension(:), allocatable :: x
+
-    double precision, dimension(:, :), intent(inout) :: A, Y
+function V(x)
-    allocate(x(n))
+    use global_vars
-    do i = 1, n
+    real :: V
-        !searching for max
+    real, intent(in) :: x
-        pivotValue = 0.0
+    if (x < a) then
-        pivot = -1
+        V = -V_0
-        do k = i, n
+    else
-            if (abs(A(k, i)) > pivotValue) then
+        V = 0.0
-                pivotValue = abs(A(k, i))
+
-                pivot = k
+    return
-            end if
+end function
-        end do
+
-        !is matrix invertible?
+recursive function jl(x, l) result(ans)
-        if (abs(pivotValue) < 1.0/10000) then
+    real, intent(in) :: x
-            print *, "matrix is invertiblle"
+    real :: ans
-        end if
+    integer, intent(in) :: l
-        !row swap
+    if (l == 0) then
-        !call dswap(1, P(i), 0, P(pivot), 0)
+        ans = sin(x) / x
-        temp = P(i)
+    else if (l == 1) then
-        P(i) = P(pivot)
+        ans = sin(x) / x**2 - cos(x) / x
-        P(pivot) = temp
+    else
-        call dswap(n, A(i,1:n), 1, A(pivot, 1:n), 1)
+        ans = (2 * (l - 1) + 1) * jl(x, l - 1) / x  - jl(x, l - 2)
-        do j = i + 1, n
+
-            A(j, i) = A(j, i) / A(i, i)
+end function
-            do k = i + 1, n
+
-                A(j, k) = A(j, k) - A(j, i) * A(i, k)
+function jl_hat(x, l)
-            end do
+    real, intent(in) :: x
-        end do
+    integer, intent(in) :: l
     real :: jl, jl_hat
-    !A = L + U - I
+    jl_hat = jl(x, l) * x
-    !Gauss
+    return
-    do k = 1, m
+end function
-        !forward
+
-        do i = 1, n
+recursive function nl(x, l) result(ans)
-            x(i) = Y(P(i), k)
+    real, intent(in) :: x
-            x(i) = x(i) - ddot(i-1, A(i,1:i-1), 1, x(1:i-1), 1)
+    real :: ans
-            !do j = 1, i - 1
+    integer, intent(in) :: l
-            !    x(i) = x(i) - A(i, j) * x(j)
+    if (l == 0) then
-            !end do
+        ans = -cos(x) / x
-        end do
+    else if (l == 1) then
-        !backward
+        ans = -cos(x) / x**2 - sin(x) / x
-        do i = n, 1, -1
+    else
-            if (i < n) then
+        ans = (2 * (l - 1) + 1) * nl(x, l - 1) / x  - nl(x, l - 2)
-                x(i) = x(i) - ddot(n-i, x(i+1:n), 1, A(i,i+1:n), 1)
+
-            end if
+end function
-            !do j = i + 1, n
+
-            !    x(i) = x(i) - A(i, j) * x(j)
+function nl_hat(x, l)
-            !end do
+    real, intent(in) :: x
-            x(i) = x(i) / A(i, i)
+    integer, intent(in) :: l
-        end do
+    real :: nl, nl_hat
-        !answer
+    nl_hat = nl(x, l) * x
-        do i = 1, n
+end function
-            Y(i, k) = x(i)
+
-        end do
+function dfac(x)
     integer, intent(in) :: x
-    deallocate(x)
+    integer :: i, dfac
     dfac = mod(x, 2)
-
+    i = dfac + 2
-program lu_decomposition
+    do while (i <= x)
         dfac = dfac * i
-    integer :: n, m, i, inp, outp
+        i = i + 2
-    integer, dimension(:), allocatable :: P, sP
+
-    double precision, dimension(:,:), allocatable :: A, sA, tempY, Y, sY
+    return
-    interface
+end function
-        subroutine LU(n, m, A, P, Y)
+
-            integer, intent(in) :: n, m
+function F(x, delta)
-            integer :: i, j, k, pivot, temp
+    use global_vars
-            double precision :: pivotValue
+    real :: nl, jl, jl_hat, nl_hat, F, V
-            integer, dimension(:), intent(inout) :: P
+    real, intent(in) :: x, delta
-            double precision, dimension(:), allocatable :: x
+    F = -V(x) * (cos(delta) * jl_hat(k * x, l) - sin(delta) * nl_hat(k * x, l))**2 / k
-            double precision, dimension(:, :), intent(inout) :: A, Y
+    return
-        end subroutine LU
+end function
-    end interface
+
-    !input
+subroutine RungeKutta(h, N, Y)
-    inp = 42
+
-    outp = 43
+    integer :: i, N
-    open(unit = inp, file = "input.txt")
+    real :: k, F, h, x, S(4), Y(N)
-    read(inp, *) n, m
+    do i = 1, N - 1
-    allocate(P(n))
+        x = h * i
-    do i = 1, n
+        S(1) =       F(x          , Y(i)                 )
-        P(i) = i
+        S(2) = 2.0 * F(x + h / 2.0, Y(i) + h * S(1) / 2.0)
         S(3) = 2.0 * F(x + h / 2.0, Y(i) + h * S(2) / 2.0)
-    allocate(sP(n))
+        S(4) =       F(x + h      , Y(i) + h * S(3)      )
-    sP = P
+        Y(i + 1) = Y(i) + h * sum(S) / 6.0
-    allocate(A(n, n))
+
-    allocate(sA(n, n))
+
-    allocate(Y(n, m))
+
-    allocate(sY(n, m))
+program phase
-    allocate(tempY(m, n))
+    use global_vars
-    read (inp, *) A
+
-    A = transpose(A)
+    integer :: i, out = 42, dfac
-    sA = A
+    real :: V, h
-    read (inp, *) tempY
+    real, allocatable :: DELTA(:)
-    Y = transpose(tempY)
+    allocate(DELTA(N))
-    sY = Y
+    h = R / N
-    close(inp)
+    k = lk
-    open(unit = outp, file = "output.txt")
+    open(unit = out, file = "output.txt")
-    !LU decomposition
+    write(out, *) "k               delta"
-    call LU(n, m, A, P, Y)
+    do while (k <= rk)
-    !LU output
+        DELTA(1) = -V(0.0) * (k * h)**(2 * l + 1) * h**2 / ((2 * l + 3) * dfac(2*l + 1)**2)
-    write(outp, *) "A"
+        call RungeKutta(h, N, DELTA)
-    do i = 1, n
+        write(out, *) k, DELTA(N)
-        write(outp, *) A(i, :)
+        k = k + dk
     end do
-    write(outp, *) "P"
+    close(out)
-    write(outp, *) P
+