Linear shooting method to solve a B.V.P

Question

I have wrote a code to approximate the solution of a boundary value problem:

  x'' = p(t)x'(t)+q(t)x(t)+r(t)
  x(b) = beta in [a,b] 

by using Runge-Kutta method of order N = 4.

Do you have any idea to make it short? I don't know how to call a subroutine as an argument in another subroutine. I had to copy the rk4 subroutine for each equation of U and V.

program test_linsht
! 
! Linear Shooting method
! To approximate the solution of boundary value problem 
!           x'' = p(t)x'(t)+q(t)x(t)+r(t)
!'          x(b) = beta in [a,b] 
! and using Runge-Kutta method of order N = 4
! 
!
!
implicit none
integer,parameter :: n = 2
integer,parameter :: nstep = 20
real(8) :: U(nstep)             ! U contains the solutions at each step
real(8) :: V(nstep),out(nstep)  ! V just like U for another equation 
real(8) :: a,b,h,alph,beta,w
real(8) :: x(0:n)
integer :: i
a = 0.d0
b = 4.0
alph = 125.d-2
beta  = -95.d-2
h = (b-a)/nstep
! x = t0 x0 x'0
x = (/0.d0, alph, 0.d0/)

call rks4U(n,h,x,nstep,U)
a = 0.d0
b = 4.0
x = (/0.d0, 0.d0, 1.d0 /)
h = (b-a)/nstep

call rks4V(n,h,x,nstep,V)

! calculate the solution to the boundary value broblem
w = (beta - U(nstep))/V(nstep)
write (*,*) w,U(nstep),V(nstep)
!w = 0.485884
out = U + w * V
do i = 1,nstep
    write(*,100) i,i*h,U(i), w*V(i), out(i)
enddo
100 format (1x,I5,f10.2,f15.6,f15.6,f15.6)
end program

subroutine xprsysU(n,x,f)
! x prime of system of equations

    real(8), dimension(0:n) :: x,f
    integer :: n
    ! time is introduced as a new differential equation
    f(0) = 1.d0
    f(1) = x(2)
    f(2) = 2.*x(0)/(1+x(0)*x(0))*x(2)- 2./(1.+x(0)*x(0))*x(1) + 1.0
end subroutine

subroutine xprsysV(n,x,f)
! x prime of system of equations

    real(8), dimension(0:n) :: x,f
    integer :: n
    ! time is introduced as a new differential equation
    f(0) = 1.d0
    f(1) = x(2)
    f(2) = 2.*x(0)/(1+x(0)*x(0))*x(2)- 2./(1.+x(0)*x(0))*x(1) 
end subroutine



subroutine rks4U(n,h,x,nstep,xout)

implicit none
      real(8) :: x(0:n)
      real(8) :: y(0:n), f(0:n,4)  
      real(8) :: xout(nstep) 
      integer :: i,k,n,nstep
      real(8) :: h
      !write(*,200) "k","t","x","y"
      !write(*,100) 0,x
      f = 0.
out:  do k = 1,nstep
          call xprsysU(n,x,f(0,1))
in1:      do i = 0,n
              y(i) = x(i) + 0.5*h*f(i,1)
              !print *, f(i,1)
          end do in1          
          call xprsysU(n,y,f(0:n,2))
in2:      do i = 0,n
              y(i) = x(i) + 0.5*h*f(i,2)      
          end do in2
          call xprsysU(n,y,f(0:n,3))    
in3:      do i = 0,n
              y(i) = x(i) + h*f(i,3)       
          end do in3
          call xprsysU(n,y,f(0:n,4))    
in4:      do i = 0,n
              x(i) = x(i) + (h/6.0)* (f(i,1) + 2.0*(f(i,2) + f(i,3)) + f(i,4)) 
          end do in4
          xout(k) = x(1)          
      end do out
      !100 format (1x,I5,f10.2,f15.8,f15.8)
      !200 format (1x,A4,2A10,A16)
end subroutine rks4U 


subroutine rks4V(n,h,x,nstep,xout)

implicit none
      real(8) :: x(0:n)
      real(8) :: y(0:n), f(0:n,4)  
      real(8) :: xout(nstep) 
      integer :: i,k,n,nstep
      real(8) :: h
      !write(*,200) "k","t","x","y"
      !write(*,100) 0,x
      f = 0.
out:  do k = 1,nstep
          call xprsysV(n,x,f(0,1))
in1:      do i = 0,n
              y(i) = x(i) + 0.5*h*f(i,1)
              !print *, f(i,1)
          end do in1          
          call xprsysV(n,y,f(0:n,2))
in2:      do i = 0,n
              y(i) = x(i) + 0.5*h*f(i,2)      
          end do in2
          call xprsysV(n,y,f(0:n,3))    
in3:      do i = 0,n
              y(i) = x(i) + h*f(i,3)       
          end do in3
          call xprsysV(n,y,f(0:n,4))    
in4:      do i = 0,n
              x(i) = x(i) + (h/6.0)* (f(i,1) + 2.0*(f(i,2) + f(i,3)) + f(i,4)) 
          end do in4
          xout(k) = x(1)          
      end do out
      !100 format (1x,I5,f10.2,f15.8,f15.8)
      !200 format (1x,A4,2A10,A16)
end subroutine rks4V 

Answer 1

To pass a function or subroutine to another, you need to declare it with an interface block. The compiler also needs to be aware of the function/subroutine you are calling, which means either using a module to hold the RK4 subroutines or, with small enough code, a contains block.

I'm going to go with the latter since it's easier here due to code size, but the strategy really isn't that different with using a module. I'm also going to snip out most of the code, as I don't think much is needed to show how to do it.

program test_linsht
   implicit none
   ...
   call rk4(n, h, x, nstep, xprsysU, U)
   call rk4(n, h, x, nstep, xprsysV, V)
   ...
contains

subroutine xprsysU(n,x,f)
! x prime of system of equations
   real(8), dimension(0:n) :: x,f
   integer :: n
   ! time is introduced as a new differential equation
   f(0) = 1.d0
   f(1) = x(2)
   f(2) = 2.*x(0)/(1+x(0)*x(0))*x(2)- 2./(1.+x(0)*x(0))*x(1) + 1.0
end subroutine

subroutine rk4(n, h, x, nstep, sub, xout)
   real(8) :: x(0:n)
   real(8) :: y(0:n), f(0:n,4)
   interface
      subroutine sub(n, x, f)
         real(8), dimension(0:n) :: x,f
         integer :: n
      end subroutine sub
   end interface
   ...
   do k=1,nstep
      call sub(n, x, f)
   end do
   ...
end subroutine
end program test_linsht

Note that the interface block has the same structure as the declaration part of the subroutine, so it should be obvious that it is necessary for the inputs of the different functions to be the same, otherwise we'd have a compiler error.

Other odds and ends:

• Fortran is 1-indexed, so it's more natural in the language to use dimension(1:n) (or even dimension(n), which is the same thing as the other), instead of defining it as dimension(0:n). Obviously it can handle an index of 0, but it's just not typical of Fortran programs.
• I don't think the labels in# and out on the do loops are needed in this case. Typically, I'd only recommend them for the cases where an inner loop requires you to increment the outer loop, but that doesn't happen here, it's just clutter.
• It's possibly due to copy-and-paste, but indentation matters for reading purposes. You've got some parts indented and others not, some indented further than others.
• Fortran allows array math, use it. You can turn those do loops in the RK4 routine into y(:) = x(:) + h*f(:,1) etc.
• It is not generally recommended to use real(8) syntax, it's possibly not what you think it is. The Fortran 90 standard introduced kind notation to give you portable definition of arbitrary precision numbers.
• With a contains statement, the parameters are global, so you would not need to pass them to the subroutines. With a module, you would have to declare those parameters in the module as global variables.
• I'd recommend declaring variables as intent(in) and intent(out) whenever and wherever possible, as it helps keep you straight with regards to variables usage (and possibly help optimize the code).