I wrote a fortran program which numerically solve the non-linear problem -f(x)+sin(f(x))-gamma=0, x is defined in [-l, l] with conditions f'(-l)=he, f'(l)=he, where gamma and he are constants.
This is the subprograms:
This is the main program code:
If I choose n<300 then I receive very big numbers in the results and should not be. Also subprogram Solver not calculates correctly phi. That what I receive in results is the Initial Approximation. I think I have a logical error somewhere, but I can't find it. Could someone help?
This is the subprograms:
Code:
module module_BVP
public :: Step, Tridiagonal_Solv_Mod, grid, Init_Approx, Calculate_del, RQ, Solver, phi_dif, ReadFromFile, WriteToFile, channel
contains
!========== function Step =============================================================================================
function Step(n, a, b) result(h)
implicit none
integer n
double precision b, a
double precision h
common /DIM/ D(1)
integer D
common /PAR/ P(4)
double precision P
n = D(1)
a = P(1)
b = P(2)
h = (b - a) / (n - 1)
end function Step
!######################################################################################################################
!============== Subroutine Tridiagonal_Solv ===========================================================================
subroutine Tridiagonal_Solv_Mod(n, a, b, he, gamma, x, y)
implicit none
integer n
integer i
integer label
double precision a, b, he, gamma
double precision h
double precision s(n), t(n)
double precision down(n), main(n), up(n)
double precision r(n), q(n)
double precision x(n)
double precision y(n)
double precision phi0(n)
common /DIM/ D(1)
integer D
common /PAR/P(4)
double precision P ! Параметри на задачата
common /NI/ NI(2)
integer NI
double precision h32
double precision h12
double precision h2
double precision h22
n = D(1)
a = P(1)
b = P(2)
he = P(3)
gamma = P(4)
label = NI(2)
h = Step(n, a, b)
h32 = 3.0d0/(2.0d0*h)
h2 = h**2
h12 = 1.0d0/(h2)
h22 = 2.0d0/h2
call grid(n, x)
call Init_Approx(n, x, label, phi0)
call RQ(n, phi0, r, q)
down(1) = -h32
do i = 2, n - 1
down(i) = -h12
end do
down(n) = 1.0d0/(2.0d0*h)
main(1) = 2.0d0/h
do i = 2, n - 1
main(i) = q(i) + h22
end do
main(n) = 2.0d0/h
up(1) = -1.0d0/(2.0d0*h)
do i = 2, n - 1
up(i) = -h12
end do
up(n) = h32
s(1) = -(main(1)*up(2) - main(2)*up(1))/(down(1)*up(2) - down(2)*up(1))
t(1) = (up(2)*r(1) - up(1)*r(2))/(down(1)*up(2) - down(2)*up(1))
s(2) = -(up(1)*down(2) - up(2)*down(1))/(main(1)*down(2) - main(2)*down(1))
t(2) = (down(2)*r(1) - down(1)*r(2))/(main(1)*down(2) - main(2)*down(1))
do i = 3, n - 2
s(i) = -up(i)/(main(i) + down(i)*s(i - 1))
t(i) = (r(i) - down(i)*t(i - 1))/(main(i) + down(i)*s(i - 1))
end do
y(n) = ((r(n - 1) - t(n - 2)*down(n - 1))*(s(n - 2)*down(n) + main(n)) - (r(n) - t(n - 2)*down(n))*(s(n - 2)*down(n - 1) + main(n - 1)))/(up(n - 1)*(down(n)*s(n - 2) + main(n)) - up(n)*(down(n - 1)*s(n - 2) + main(n - 1)))
y(n - 1) = (up(n - 1)*(r(n) - down(n)*t(n - 2)) - up(n)*(r(n - 1) - down(n - 1)*t(n - 2)))/(up(n - 1)*(down(n)*s(n - 2) + main(n)) - up(n)*(down(n - 1)*s(n - 2) + main(n - 1)))
do i = n - 2, 1, -1
y(i) = s(i)*y(i + 1) + t(i)
end do
end subroutine Tridiagonal_Solv_Mod
!#####################################################################################################################
!================== Subroutine grid ==================================================================================
subroutine grid(n, x)
implicit none
integer i
common /DIM/ D(1)
integer D
integer n
common /PAR/ P(4)
double precision P
double precision a
double precision b
double precision x(n)
double precision h
n = D(1)
a = P(1)
b = P(2)
h = Step(n, a, b)
do i = 1, n - 1
x(i) = a + (i - 1)*h
end do
x(n) = b
end subroutine grid
!#####################################################################################################################
!================== Subroutine Init_Approx ===========================================================================
subroutine Init_Approx(n, x, label, phi0)
implicit none
integer n
integer i
double precision x(n)
double precision phi0(n)
double precision, parameter :: pi = 3.141592653589793d0
integer label
common /DIM/ D(1)
integer D
common /NI/ NI(2)
integer NI
n = D(1)
label = NI(2)
call grid(n, x)
select case(label)
case(1)
do i = 1, n
phi0(i) = 4.0d0*atan(exp(x(i)))
end do
case(2)
do i = 1, n
phi0(i) = 4.0d0*atan(exp(-x(i)))
end do
case(3)
do i = 1, n
phi0(i) = pi
end do
case(4)
do i = 1, n
phi0(i) = 0
end do
end select
end subroutine Init_Approx
!#####################################################################################################################
!================= Subroutine Calculate_del ==========================================================================
subroutine Calculate_del(n, a, b, phi, del)
implicit none
integer i
double precision del
double precision h
common /DIM/ D(1)
integer D
integer n
double precision phi(n)
double precision F(n)
double precision x(n)
common /PAR/P(4)
double precision P
double precision a
double precision b
double precision he
double precision gamma
double precision h2
n = D(1)
a = P(1)
b = P(2)
he = P(3)
gamma = P(4)
h = Step(n, a, b)
h2 = h**2
del = 0.0d0
call grid(n, x)
do i = 1, n
F(i) = 0
end do
F(1) = (-3.0d0*phi(1) + 4.0d0*phi(2) - phi(3))/2.0d0*h - he
del = F(1)**2
do i = 2, n - 1
F(i) = -(phi(i - 1) - 2*phi(i) + phi(i + 1))/(h2) + dsin(phi(i)) - gamma
del = del + F(i)**2
end do
F(n) = (phi(n - 2) - 4*phi(n - 1) + 3*phi(n))/2.0d0*h - he
del = del + F(n)**2
end subroutine Calculate_del
!#####################################################################################################################
!================== Subroutine RQ ================================================================================
subroutine RQ(n, y, r, q)
implicit none
integer i
common /DIM/ D(1)
integer D
integer n
double precision y(n)
double precision r(n), q(n)
double precision h
common /PAR/ P(4)
double precision P
double precision a
double precision b
double precision he
double precision gamma
double precision h2
n = D(1)
a = P(1)
b = P(2)
he = P(3)
gamma = P(4)
h = Step(n, a, b)
h2 = h**2
r(1) = -(-3.0d0*y(1) + 4.0d0*y(2) - y(3))/(2.0d0*h) + he
q(1) = dcos(y(1))
do i = 2, n - 1
r(i) = (y(n - 1) - 2.0d0*y(i) + y(i + 1))/(h2) - dsin(y(i)) + gamma
q(i) = dcos(y(i))
end do
r(n) = -(y(n - 2) - 4.0d0*y(n - 1) + 3.0d0*y(n))/(2.0d0*h) + he
q(n) = dcos(y(n))
end subroutine RQ
!#####################################################################################################################
!================== Subroutine Solver ================================================================================
subroutine Solver(n, x, y, a, b, he, gamma, phi0, phi, flag)
implicit none
integer n
integer i, k
integer label
double precision x(n)
double precision y(n)
double precision tau
double precision EPSABS
double precision tau_min
integer MITER
integer flag
double precision del0, del1
double precision error
double precision phi0(n), phi(n)
double precision a, b, he, gamma
common /PAR/P(4)
double precision P
common /NI/ NI(2)
integer NI
common /NUMER_DP/ NDP(5)
double precision NDP
common /DIM/ D(1)
integer D
n = D(1)
tau_min = NDP(2)
MITER = NI(1)
EPSABS = NDP(1)
flag = 0
a = P(1)
b = P(2)
he = P(3)
gamma = P(4)
label = NI(2)
call Init_Approx(n, x, label, phi0)
do i = 1, n
phi(i) = phi0(i)
end do
k = 0
call Calculate_del(n, a, b, phi0, del0)
error = sqrt(del0)
do while(error .gt. EPSABS)
call Tridiagonal_Solv_Mod(n, a, b, he, gamma, x, y)
do i = 1, n
phi0(i) = phi(i) + y(i)
end do
call Calculate_del(n, a, b, phi0, del1)
tau = del0 / (del0 + del1)
if(tau .lt. tau_min) then
tau = tau_min
else if(tau .gt. 0.9d0) then
tau = 1.0d0
end if
do i = 1, n
phi(i) = phi(i) + tau*y(i)
end do
k = k + 1
call Calculate_del(n, a, b, phi, del0)
error = sqrt(del0)
if(k .eq. MITER) then
flag = 1
return
end if
end do
return
end subroutine Solver
!#####################################################################################################################
!================= Subroutine phi_dif ================================================================================
subroutine phi_dif(n, phi, d_phi)
implicit none
integer i
common /DIM/ D(1)
integer D
integer n
double precision phi(n)
double precision d_phi(n)
common /PAR/ P(4)
double precision P
double precision h
double precision a
double precision b
double precision he
n = D(1)
a = P(1)
b = P(2)
he = P(3)
h = Step(n, a, b)
d_phi(1) = (-3*phi(1) + 4*phi(2) - phi(3))/(2*h)
do i = 2, n - 1
d_phi(i) = (phi(i + 1) - phi(i - 1))/(2*h)
end do
d_phi(n) = (phi(n - 2) - 4*phi(n - 1) + 3*phi(n))/(2*h)
end subroutine phi_dif
!#####################################################################################################################
!================= Subroutine ReadFromFile ===========================================================================
subroutine ReadFromFile(n, free)
implicit none
integer n, free, err, label
logical FExist
common /DIM/ D(1)
integer D
common /PAR/ P(4)
double precision P
common /NI/ NI(2)
integer NI
double precision tau_min
double precision EPSABS
common /NUMER_DP/ NDP(5)
double precision NDP
inquire(file = 'Input_data.dat', exist = FExist)
if(.not. FExist) then
write(*, *) '* Error: Data file not found...'
stop
end if
call channel(free)
open(unit=free, file = "Input_data.dat", status = "old", action = "read")
write(*, *) "Read data file..."
read(free, *) D(1)
read(free, *) P(1)
read(free, *) P(2)
read(free, *) P(3)
read(free, *) P(4)
read(free, *) NI(1)
read(free, *) NI(2)
read(free, *) NDP(1)
read(free, *) NDP(2)
close(free)
n = D(1)
label = NI(2)
end subroutine ReadFromFile
!#####################################################################################################################
!================= Subroutine WriteToFile ============================================================================
subroutine WriteToFile(n, a, b, x, phi, d_phi)
implicit none
integer free, n, i
double precision x(n)
double precision phi(n)
double precision d_phi(n)
double precision a, b
common /PAR/ P(4)
double precision P
a = P(1)
b = P(2)
call channel(free)
open(unit=free, file="Results.dat", status="replace", action="write")
write(free, *) "i x(i) phi(x(i)) phi'(x(i))"
write(free, *)
do i = 1, n
write(free, fmt = "(i5, f15.10, f15.10, f15.10)"), i, x(i), phi(i), d_phi(i)
end do
print *
print *, "The results was writen in file 'Results.dat'"
print *
close(free)
end subroutine WriteToFile
!#####################################################################################################################
!================== CHANNEL ==========================================================================================
SUBROUTINE channel (NFREE)
LOGICAL FISOPEN
INTEGER NFREE
FISOPEN = .TRUE.
NFREE = 9
DO WHILE (FISOPEN)
NFREE = NFREE + 1
INQUIRE (NFREE, OPENED = FISOPEN)
END DO
RETURN
END SUBROUTINE channel
!#####################################################################################################################
end module module_BVP
This is the main program code:
Code:
program BVP
use module_BVP
implicit none
integer nn
parameter (nn = 10000)
integer n
integer i
integer free
double precision a
double precision b
double precision he
double precision gamma
double precision x(nn)
double precision y(nn)
double precision phi0(nn)
double precision phi(nn)
double precision d_phi(nn)
integer flag
integer label
common /PAR/P(4)
double precision P
common /NI/ NI(1)
integer NI
common /NUMER_DP/ NDP(5)
double precision NDP
common /DIM/ D(1)
integer D
n = D(1)
a = P(1)
b = P(2)
he = P(3)
gamma = P(4)
flag = 0
call ReadFromFile(n, free)
call Init_Approx(n, x, label, phi0)
call Solver(n, x, y, a, b, he, gamma, phi0, phi, flag)
call phi_dif(n, phi, d_phi)
call grid(n, x)
call WriteToFile(n, a, b, x, phi, d_phi)
end program BVP
If I choose n<300 then I receive very big numbers in the results and should not be. Also subprogram Solver not calculates correctly phi. That what I receive in results is the Initial Approximation. I think I have a logical error somewhere, but I can't find it. Could someone help?