Help with Adams-Bashforth program

[ognik]

Hi - trying to write a Fortran95 program to solve a trivial 1st order DE (dy/dx=-xy), the idea being that once it works it can be used for different forms of f(x,y). I use Euler to get the starting point, then A-B's 2 step formula for the next 2 starting points. It is when I apply the A-B 4-step formula that I get obviously wrong values for Y(n+1) (using debug to compare with expected values). I've looked at this every which way until my head hurts, can't see what I'm doing wrong in the program? (code in attached zip file)

 

[Mark44]

It would be helpful if you posted the code right here, rather than posting it in a Zip file that we 1) have to download, and 2) open in some editor. If you have a question about code, it behooves you to make life as easy as possible for those providing help.

If you do this, please enclose your code with BBcode [ code ] and [ /code ] tags (without the extra spaces).

 

[ognik]

Understood thanks. -I seem to loose the code's indenting when pasting, can you advise there a way to keep that?

 

[ognik]

OK - I had 'replace tabs with spaces' set in my editor, code follows:

program AB4step
   !Summary: Apply Adams-Bashforth 4 step method to dy/dx = f(x,y) = -xy given y(0)=1
   implicit none
   integer, parameter  :: dp = selected_real_kind(15, 307)
   real(kind=dp)       :: xN,yN,xNplus1,yNplus1,xNminus1,yNminus1,exact,diff
   real(kind=dp)       :: fN, fNminus1, fNminus2, fNminus3, FUNC
   real, parameter  :: lim=4.0
   real  :: h
   integer         :: iter, numsteps
   integer, parameter  :: n=60     
   character(len=n )  :: stars=REPEAT('*',n)
   !----- header ----   
   print *, " "
   print 10, 'Program AB-4step starts, integrating from 0.0 to ', lim
   print *, 'Enter the number of steps to use (zero or negative to stop)'
   read  (*,*) numSteps  ! ensure integer number os steps, adjust h accordingly
   print *, stars
   print *, " "
   print 20, 'Iter.','X(n+1)','F(n)  ', 'Y(n+1)  ', 'Exact  ','Diff  '
   !----- end header ----   
   do while (numSteps.GT.0)     ! allows trials of diff. values of step size h
   !----- Initialise ----   
     h=lim/numSteps
     yNminus1=1.0  ! y(0)  starting point, n=0
     xNminus1=0.0  ! x=0
     fNminus2=0.0
     fNminus1=FUNC(xNminus1, yNminus1)
     yN=yNminus1+(h*fNminus1)  ! using Eulers 1-step to get 2nd starting value, n=1
     xN=xNminus1+h
   !----- Calculate next 2 points using AB 2step (n = 2,3) --------   
     Do iter=1,2   
       fN=FUNC(xN,yN)
       exact=exp(-(xN**2)/2.0)
       diff=exact-yN
       print 50, iter, xN, fN, yN, exact, diff
       yNplus1=yN+(h*((1.5*fN)-(0.5*fNminus1)))
       xNplus1=xN+h
       yNminus1=yN
       xNminus1=xN
       fNminus3=fNminus2
       fNminus2=fNminus1
       fNminus1=fN
       yN=yNplus1
       xN=xNplus1   
     end do
   !----- Calculate rest of points using AB 4 step (n = 4 onward --------   
     Do iter=iter,numsteps   
       fN=FUNC(xN, yN)
       exact=exp(-(xN**2)/2.0)
       diff=exact-yN
       print 50, iter, xN, fN, yN, exact, diff
       yNplus1=yN+((h/24.0)*((55.0*fN)-(59.0*fNminus1)+(37.0*fNminus2)-(9.0*fNminus3)))
       xNplus1=xN+h
       yNminus1=yN
       xNminus1=xN
       fNminus3=fNminus2
       fNminus2=fNminus1
       fNminus1=fN
       yN=yNplus1
       xN=xNplus1   
     end do
     print *, 'The step size used above was ',h
     print *, 'Enter another number of steps value (zero or negative to stop)'
     read  (*,*) numSteps
   end do
!--- admin ----
   print *, 'User selected end'   
   print *, stars  
10   format (A, F8.3)
20   format (A, 5X, 5(A,8X))
50   format (I2, 5X, F8.4, 2X, 4(F15.9, 3X))
end program AB4step
!---------Functions -----------
function FUNC (x, y)
integer, parameter      :: dp = selected_real_kind(15, 307)
real(kind=dp)           :: x, y, FUNC
   FUNC=-x*y
end function

 

[Mark44]

Here's what I think might be going on. You are using n as if it were a variable, but you have declared it as a parameter, which means that its value doesn't change. In the line that starts with "Do iter = 1, 2" your comment above talks about AB 2step (n = 2, 3). After that loop is finished, iter will be 3, but during the loop iter's value will be 1 in the first iteration and 2 in the second iteration.

In the next loop that starts with "Do iter = iter, numsteps" the comment above says (n = 4 onward), but iter's value starts with 3, not 4.

Also, in the code below, it's possible that you might not be implementing the algorithm correctly.

Do iter=iter,numsteps  
       fN=FUNC(xN, yN)
       exact=exp(-(xN**2)/2.0)
       diff=exact-yN
       print 50, iter, xN, fN, yN, exact, diff
       yNplus1=yN+((h/24.0)*((55.0*fN)-(59.0*fNminus1)+(37.0*fNminus2)-(9.0*fNminus3)))
       xNplus1=xN+h
       yNminus1=yN
       xNminus1=xN
       fNminus3=fNminus2
       fNminus2=fNminus1
       fNminus1=fN
       yN=yNplus1
       xN=xNplus1  
     end do[/quote]
What I would do is verify that the first loop is producing the right numbers, and if so, calculate by hand what the first one or two iterations of the loop above should be producing, and compare these values with what your code is producing.

 

[ognik]

Embarrassed, me. The program is working, but it turns out the algorithm oscillates above a certain step-size. Sorry for any hassle.