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