forums.silverfrost.com

[ognik]

Hi - student trying to write a 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 next apply the 4-step formula that I get obviously wrong values for Y(n+1) (using debug). I've looked at this every which way until my head hurts, can't see what I'm doing wrong. Help would be nice javascript:emoticon('')
[code:1:f5c5c2e3a5]
program AB4step
!Summary: Apply Adams-Bashforth 4 step method to dy/dx = -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, 3

[DanRRight]

Typically students like you who do not care what they post do not get much attention on any forums

[ognik]

? What have I done wrong?

[mecej4]

[ognik]

My apologies and thanks for the advice. I did review the code before submitting - and it was fine in review, however after submitting it lost the last few lines, I will remember the length limit.

Actually my program turned out to be correct - I was simply using too large a steps-size, which appears to cause A-B 4th to oscillate. But thanks for the advice - I am learning Fortan and clearly using arrays is better than my clumsy names .