# Starnge problem

I’m writing a rpogram in Fortran for computing Boundary Value Problem ODE.
Here is my code:

module_BVP.f90

``````module module_BVP

public :: Step, Tridiagonal_Solv, ReadFromFile, WriteToFile, channel

contains

!==========  Step =============================================================================================
!           - a,     - b
!     - n.
function Step(n, a, b) result(h)
implicit none

integer,			 intent(in)	:: n
double precision,	 intent(in) :: b, a
real							:: h

h = (b-a) / (n-1)

end function Step

!======================================================================================================================

!============== Subroutine Values =================

subroutine Tridiagonal_Solv(n, a, b, alpha, beta)
implicit none
integer						:: i, j, n
double precision					:: a, b !
real,			  dimension(:)		:: q, w, t, y
real,			  dimension(:)		:: sigma, teta
double precision					:: p, l, z
double precision,		  intent(in) :: alpha, beta
double precision, dimension(n)             :: x, r, vector_B
double precision, dimension(n,n)			 :: matrix_A

x(0) = a
x(n) = b
p = -1 / (Step(n, a, b)*Step(n, a, b))
z = 2 / (Step(n, a, b)*Step(n, a, b)) + 1
l = 2 / Step(n, a, b)

!    x(i)  r(x(i))
do i = 2, n-1
x(i) = a + i*Step(n, a, b)
r(i) = 2*sin(x(i))
end do

!    nxn .
do i = 1, n
do j = 1, n
matrix_A(i,j) = 0
end do
end do

!         (   )
matrix_A(1,1) = -3 / (2*Step(n, a, b))
matrix_A(1,2) = l
matrix_A(1,3) = -1 / (2*Step(n, a, b))

!
do i = 2, n-1
matrix_A(i,i-1) = p
matrix_A(i,i) = z
matrix_A(i,i+1) = p
vector_B(i) = r(i)
end do

!          (   )
matrix_A(n,n-2) = 1 / (2*Step(n, a, b))
matrix_A(n, n-1) = l
matrix_A(n,n) = 3 / (2*Step(n, a, b))

vector_B(1) = alpha
vector_B(n) = beta

y(0) = 0
y(n+1) = 0
sigma(0) = 0
teta(0) = 0
w(1) = 0
t(n) = 0

do i = 2, n-1
w(i) = matrix_A(i,i-1)
end do

do i = 1, n
q(i) = matrix_A(i,i)
end do

do i = 1, n-1
t(i) = matrix_A(i,i+1)
end do

! ,
do i = 1, n
sigma(i) = -t(i)/(q(i) + w(i)*sigma(i-1))
teta(i) = (r(i) - w(i)*teta(i-1))/(q(i) + w(i)*sigma(i-1))
end do

y(n) = teta(n)

! ,
do i = n-1, 1, -1
y(i) = sigma(i)*y(i+1) + teta(i)
end do

end subroutine Tridiagonal_Solv
!==============================================

subroutine ReadFromFile(n, a, b, alpha, beta)
implicit none
integer				 :: n, free, err
double precision			 :: a, b
double precision			 :: alpha, beta

call channel(free)
open(unit = free, file = "Values.dat", status = "old", action = "read")
do
read(free, *, iostat = err) n, a, b, alpha, beta
if(err == 0) then
print *
print *, "     n     a           b            alpha     beta"
print "(i7, f8.2, f16.10, f11.3, f10.3)", n, a, b, alpha, beta
print *
print *, "Data reading complete successfully!"
exit
else if(err < 0) then
print *, "Error in data reading!"
exit
end if
end do
close(free)

!===============================================

!================= Subroutine _Write ========

subroutine WriteToFile(n, x, y)
implicit none
integer						   :: free, n, i
double precision, dimension(n) :: x, y

call channel(free)
open(unit = free, file = "Results.dat", status = "replace", action = "write")
write(free, *) " i      x(i)	     y(i)"
write(free, *)
do i = 1, n
write(free, fmt = "(i3, f13.8, f13.8)"), i, x(i), y(i)
end do
print *
print *, "The results was writen in file 'Results.dat'"
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
``````

BVP.f90

``````program BVP
use module_BVP

implicit none

integer				     :: n
double precision			     :: a, b
double precision			     :: alpha, beta
double precision, dimension(n) :: x, y

call ReadFromFile(n, a, b, alpha, beta)
call Tridiagonal_Solv(n, a, b, alpha, beta)
call WriteToFile(n, x, y)

end program BVP
``````

There is no error when I compile it. But when run it in the Results.dat file there is no results writen. Only zeros in all columns. Why is that happen.

Hi ferry2,

You need to pass in your main program’s “x” and “y” variables to Tridiagonal_Solv. Otherwise, the “x” and “y” in Tridiagonal_Solv are local to the subroutine.

Though, I’m a bit surprised the code made it to writing out a file. It should get a segmentation violation in Tridiagonal_Solv since several of your arrays are assumed-shape but not used as arguments. Also, since x and y’s size wont be known until the file is read in, they should be allocatable arrays. Finally, be sure to declare the arrays with a “0” lower bound. In Fortran, the default is 1.

• Mat

Now my code is:

module_BVP.f90

``````module module_BVP

public :: Step, Tridiagonal_Solv, ReadFromFile, WriteToFile, channel

contains

!========== Функция Step =============================================================================================
! Изчисляваме големината на стъпката по зададени ляв край на интервала - a, десен край на интервала - b
! и брой на итерациите - n.
function Step(n, a, b) result(h)
implicit none

integer,			 intent(in)	:: n
double precision,	 intent(in) :: b, a
real							:: h

h = (b-a) / (n-1)

end function Step

!======================================================================================================================

!============== Subroutine Values =====================================================================================
! Подпрограма, в която изчисляваме всички стойности на x_i
subroutine Tridiagonal_Solv(n, a, b, alpha, beta, x, y)
implicit none
integer										  :: i, j, n ! i и j са броячи за циклите, а n е размерността на матриците и векторите
double precision							  :: a, b ! а и b са съответно ляв и десен край на интервала
real, allocatable, dimension(:)				  :: q, w, t ! q, w, t - Три едномерни масива за коефициентите на системата;

real, allocatable, dimension(:)				  :: sigma, teta ! Eдномерни масиви за векторите от прогонъчни коефициенти
double precision							  :: p, l, z ! Променливи, в които се съхраняват изчислени изрази
double precision							  :: alpha, beta ! alpha и beta са стойностите на граничните условия
double precision, allocatable, dimension(:)	  :: x, r, vector_B, y ! х е вектор за променливата "х", r е вектор за функцията r(x)
!y - Eдномерен масив за вектора от неизвестни
double precision, allocatable, dimension(:,:) :: matrix_A ! A e тридиагонална матрица, в която

allocate(q(n), w(n), t(n), y(0:n+1), sigma(0:n), teta(0:n), x(0:n), r(n), vector_B(n), matrix_A(n,n))

x(0) = a
x(n) = b
p = -1 / (Step(n, a, b)*Step(n, a, b))
z = 2 / (Step(n, a, b)*Step(n, a, b)) + 1
l = 2 / Step(n, a, b)

! Пресмятаме стойностите на x(i) и r(x(i))
do i = 2, n-1
x(i) = a + i*Step(n, a, b)
r(i) = 2*sin(x(i))
end do

! Създаваме си нулева nxn матрица.
do i = 1, n
do j = 1, n
matrix_A(i,j) = 0
end do
end do

! Записваме елементите стоящи на първия ред на матрицата (от първото гранично условие)
matrix_A(1,1) = -3 / (2*Step(n, a, b))
matrix_A(1,2) = l
matrix_A(1,3) = -1 / (2*Step(n, a, b))

! Записваме всички останали елементи
do i = 2, n-1
matrix_A(i,i-1) = p
matrix_A(i,i) = z
matrix_A(i,i+1) = p
vector_B(i) = r(i)
end do

! Записваме и елементите стоящи на последния ред на матрицата (от второто гранично условие)
matrix_A(n,n-2) = 1 / (2*Step(n, a, b))
matrix_A(n, n-1) = l
matrix_A(n,n) = 3 / (2*Step(n, a, b))

vector_B(1) = alpha
vector_B(n) = beta

y(0) = 0
y(n+1) = 0
sigma(0) = 0
teta(0) = 0
w(1) = 0
t(n) = 0

do i = 2, n-1
w(i) = matrix_A(i,i-1)
end do

do i = 1, n
q(i) = matrix_A(i,i)
end do

do i = 1, n-1
t(i) = matrix_A(i,i+1)
end do

! Цикъл, в който изчисляваме прогонъчните коефициенти
do i = 1, n
sigma(i) = -t(i)/(q(i) + w(i)*sigma(i-1))
teta(i) = (r(i) - w(i)*teta(i-1))/(q(i) + w(i)*sigma(i-1))
end do

y(n) = teta(n)

! Цикъл, в който изчисляваме неизвестните
do i = n-1, 1, -1
y(i) = sigma(i)*y(i+1) + teta(i)
end do

end subroutine Tridiagonal_Solv
!=====================================================================================================================

subroutine ReadFromFile(n, a, b, alpha, beta)
implicit none
integer						 :: n, free, err
double precision			 :: a, b
double precision			 :: alpha, beta

call channel(free)
open(unit = free, file = "Values.dat", status = "old", action = "read")
do
read(free, *, iostat = err) n, a, b, alpha, beta
if(err == 0) then
print *
print *, "     n     a           b            alpha     beta"
print "(i7, f8.2, f16.10, f11.3, f10.3)", n, a, b, alpha, beta
print *
print *, "Data reading complete successfully!"
exit
else if(err < 0) then
print *, "Error in data reading!"
exit
end if
end do
close(free)

!=====================================================================================================================

!================= Subroutine _Write =================================================================================

subroutine WriteToFile(n, x, y)
implicit none
integer										:: free, n, i
double precision, allocatable, dimension(:) :: x, y

allocate(x(n), y(n))

call channel(free)
open(unit = free, file = "Results.dat", status = "replace", action = "write")
write(free, *) " i      x(i)	     y(i)"
write(free, *)
do i = 1, n
write(free, fmt = "(i3, f13.8, f13.8)"), i, x(i), y(i)
end do
print *
print *, "The results was writen in file 'Results.dat'"
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
``````

BVP.f90

``````program BVP
use module_BVP

implicit none

integer										:: n
double precision							:: a, b
double precision							:: alpha, beta
double precision, allocatable, dimension(:) :: x, y

allocate(x(n), y(n))

call ReadFromFile(n, a, b, alpha, beta)
call Tridiagonal_Solv(n, a, b, alpha, beta, x, y)
call WriteToFile(n, x, y)

end program BVP
``````

But I’m geting a very strange result. The data in result file is:

i x(i) y(i)

1 0.05208331 760.09265128
2 0.00005207 760.09240824
3 0.00000000*************
4**************************
5 0.00000000 0.00000000

Hi ferry2,

The stars indicate that your number is too wide for the width provided, however, this is not your main issue.

You want to allocate x and y only once, preferably in your main program. Hence, remove the allocatable clause from x and y’s declaration in WriteToFile and Tridiagonal_Solv as well as removing them from the allocate statements. By allocating x and y in WriteToFile, you’re actually printing our junk values.

Also, since n isn’t known until after you read in your values, you need to wait to allocate x and y until after ReadFromFile is called. Otherwise n’s value is junk.

Since x and y don’t use element zero, you probably want to change “x(0)=a” to “x(1)=a” and remove “y(0)=0”. Also, since the value of Step never changes, why not call it once and store it’s value in a local variable?

• Mat

@mkcolg thankc for the tips. I edit my program and now the problem is that in the “Results.dat” file I’m geting only zeros. Maybe I have a logical error.

module_BVP.f90

``````	module module_BVP

public :: Step, Tridiagonal_Solv, ReadFromFile, WriteToFile, channel

contains

!========== ??????? Step =============================================================================================
! ??????????? ?????????? ?? ???????? ?? ???????? ??? ???? ?? ????????? - a, ????? ???? ?? ????????? - b
! ? ???? ?? ?????????? - n.
function Step(n, a, b) result(h)
implicit none

integer,			 intent(in)	:: n
double precision,	 intent(in) :: b, a
real							:: h

h = (b-a) / (n-1)

end function Step

!======================================================================================================================

!============== Subroutine Values =====================================================================================
! ???????????, ? ????? ??????????? ?????? ????????? ?? x_i
subroutine Tridiagonal_Solv(n, a, b, alpha,beta)
implicit none

integer						   :: n, i
double precision, intent(in)   :: a, b, alpha, beta
real						   :: h ! Ïðîìåíëèâà, â êîÿòî ùå ïàçèì ñòîéíîñòòà íà ñòúïêàòà
double precision, dimension(0:5) :: x, r, u, v, w, z, s, t, y

h = Step(n, a, b)

x(1) = a
r(1) = alpha
do i = 2, n-1
x(i) = a+i*h
r(i) = 2*sin(x(i))
end do
r(n) = beta
x(n) = b

do i = 2, n-1
u(i) = -1/(h*h)
end do

u(n) = 2/h

v(1) = -3/(2*h)
do i = 2, n-1
v(i) = 2/(h*h) + 1
end do
v(n) = 3/(2*h)

w(1) = 2/h
do i = 2, n-1
w(i) = -1/(h*h)
end do

s(0) = 0
t(0) = 0
do i = 1, n
s(i) = -w(i)/(v(i) + u(i)*s(i-1))
t(i) = (r(i) - u(i)*t(i-1))/(v(i) + u(i)*s(i-1))
end do

y(n) = t(n)
do i = n-1, 1, -1
y(i) = s(i)*y(i+1) + t(i)
end do

end subroutine Tridiagonal_Solv
!=====================================================================================================================

subroutine ReadFromFile(n, a, b, alpha, beta)
implicit none
integer						 :: n, free, err
double precision			 :: a, b
double precision			 :: alpha, beta

call channel(free)
open(unit=free,file="Values.dat",status = "old", action = "read")
do
read(free, *, iostat = err) n, a, b, alpha, beta
if(err == 0) then
print *
print *, "     n     a           b            alpha     beta"
print "(i7, f8.2, f16.10, f11.3, f10.3)", n, a, b, alpha, beta
print *
print *, "Data reading complete successfully!"
exit
else if(err < 0) then
print *, "Error in data reading!"
exit
end if
end do
close(free)

!=====================================================================================================================

!================= Subroutine _Write =================================================================================

subroutine WriteToFile(n, x, y)
implicit none
integer										:: free, n, i
double precision, dimension(:) :: x, y

call channel(free)
open(unit=free,file="Results.dat",status="replace",action="write")
write(free, *) " i      x(i)	     y(i)"
write(free, *)
do i = 1, n
write(free, fmt = "(i3, f13.8, f13.8)"), i, x(i), y(i)
end do
print *
print *, "The results was writen in file 'Results.dat'"
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
``````

BVP.f90

``````	program BVP
use module_BVP

implicit none
integer, parameter			   :: nn = 100
integer						   :: n
double precision			   :: a, b
double precision			   :: alpha, beta
double precision, dimension(nn):: x, y

call ReadFromFile(n, a, b, alpha, beta)
call Tridiagonal_Solv(n, a, b, alpha,beta)
call WriteToFile(n, x, y)

end program BVP
``````

Hi ferry2,

I’m assuming your a student, or at least new to programming, so instead of just telling what’s wrong, I wan you to try and figure this out yourself.

Hint: Reread my first post and then look over your main program again.

• Mat