Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help with Adams-Bashforth program

  1. Feb 23, 2015 #1
    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)
     

    Attached Files:

  2. jcsd
  3. Feb 23, 2015 #2

    Mark44

    Staff: Mentor

    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).
     
  4. Feb 23, 2015 #3
    Understood thanks. -I seem to loose the code's indenting when pasting, can you advise there a way to keep that?
     
  5. Feb 23, 2015 #4
    OK - I had 'replace tabs with spaces' set in my editor, code follows:
    Code (Text):

    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
     
     
  6. Feb 24, 2015 #5

    Mark44

    Staff: Mentor

    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.
    Code (Text):
    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.
     
  7. Feb 24, 2015 #6
    Embarrassed, me. The program is working, but it turns out the algorithm oscillates above a certain step-size. Sorry for any hassle.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Help with Adams-Bashforth program
  1. Programming Help (Replies: 4)

  2. Programming help (Replies: 2)

  3. Programing help (Replies: 8)

  4. Help with programming (Replies: 5)

Loading...