Tags:
1. Feb 23, 2015

### 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)

#### Attached Files:

• ###### AB-4step-9.zip
File size:
1.2 KB
Views:
68
2. Feb 23, 2015

### 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).

3. Feb 23, 2015

### ognik

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

4. Feb 23, 2015

### ognik

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)
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  '
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)'
end do
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

5. Feb 24, 2015

### 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.

6. Feb 24, 2015

### ognik

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