Logical error

ferry2 · Jul 11, 2012

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:

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?

melmacianalf · Jul 11, 2012

ferry2: May be you can post a sample 'Input_data.dat' ?

Gordon Shumway

ferry2 · Jul 11, 2012

I attached at the end of the post.

gummibaer · Jul 12, 2012

Well, Ferry, this looks like a problem that would need some effort to work out. I doubt, that anyone will be willing to spend a few hours of time to worm his / her way into your code and find out why you do not receive the results that you expect. You know, nobody but you knows what your code is supposed to do.

My maybe less than 2 cents worth of advice would be, that you yourself debug your code, you are the only living expert on it. If you do not have an online debugger available in your IDE, then printout the values of the variables and constants that contribute to your results. Check which of them shows data as you expect them to be and which do not. For those that do not, printout the values again that affect them and check again which data are reasonable and which not. And so on. And so on. This might be tedious, but this the only relyable way to solve your problem. And it sure will solve your problem if you have enough perseverance. It is not so very uncommon that the programmer finds out after he finished testing his code, that while debugging he did the computing for his first sets of data more or less completely by hand / calculator parallel to that what his computer did.

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Logical error

ferry2

Technical User

melmacianalf

Technical User

ferry2

Technical User

gummibaer

Programmer

Similar threads

Part and Inventory Search

Sponsor