|
forums.silverfrost.com Welcome to the Silverfrost forums
|
View previous topic :: View next topic |
Author |
Message |
ognik
Joined: 04 Feb 2015 Posts: 19
|
Posted: Mon Feb 23, 2015 2:05 am Post subject: Adams-Bashforth 4 step method diverging |
|
|
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 |
|
Back to top |
|
|
DanRRight
Joined: 10 Mar 2008 Posts: 2816 Location: South Pole, Antarctica
|
Posted: Mon Feb 23, 2015 6:15 am Post subject: |
|
|
Typically students like you who do not care what they post do not get much attention on any forums |
|
Back to top |
|
|
ognik
Joined: 04 Feb 2015 Posts: 19
|
Posted: Mon Feb 23, 2015 11:35 pm Post subject: |
|
|
? What have I done wrong? |
|
Back to top |
|
|
mecej4
Joined: 31 Oct 2006 Posts: 1886
|
Posted: Tue Feb 24, 2015 1:51 pm Post subject: Re: |
|
|
ognik wrote: | ? What have I done wrong? |
Here are a few hints.
Your code is pasted in-line, causing the indentation to get lost, and making the code unpleasant to read. Use code.../code tags around your code fragments (surrounded by brackets).
Your code is longer than the forum limit, so the end part is trimmed (have you noticed this?), making the program incomplete and unusable. Keep them short, or place longer files elsewhere and provide a link.
You state that the result is wrong, but do not state what you obtained, what you think the correct result is, and how close you expect the two to be to each other. Explain your reasoning.
Here is a fragment of code for you to work with. Add declarations before it, complete the AB-4 part, add output statements after the fragment, and try. After you get the code to run and produce satisfactory results, you can consider using x, y and f values at only four x-levels instead of using arrays (these are convenient and make the code easy to read, but could be considered wasteful of memory).
Code: |
!
! Initial condition
x(0)=0.; y(0)=1.; h=0.01
!
! Euler step
!
f(0)=fn(x(0),y(0))
x(1)=h; y(1)=y(0)+h*f(0)
f(1)=fn(x(1),y(1))
!
! AB-2 steps
!
DO i=2,3
x(i)=x(i-1)+h
y(i)=y(i-1)+(3.*f(i-1)-f(i-2))*0.5*h
f(i)=fn(x(i),y(i))
END DO
!
! AB-4 steps
!
DO i=4,N
|
|
|
Back to top |
|
|
ognik
Joined: 04 Feb 2015 Posts: 19
|
Posted: Sun Mar 01, 2015 10:23 pm Post subject: |
|
|
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 . |
|
Back to top |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|