numeric analysis from vb to fortran

[oswaler]

I have been writing in VB for a long time and need to get strong in FORTRAN quick (I've never used FORTRAN). I wrote a program in VB that integrates (using numeric analysis) x^2 from 0 to 1 (which should be 0.3333) and it works fine. As the number of iterations through the loop increases, the result becomes more and more accurate.

I tried to run the same program in FORTRAN with slight syntax changes and I get a very weird result. If I have the loop iterate 1000 times I get a result of 0.332828. If i run the loop 10 000 times I get 0.33293. If i run the loop 1 million times I get 0.3386. If i run the loop 1 billion times I get 3.05E-05. As I increase the number of iterations the accuracy decreases. The strange thing is that the answer is getting bigger for a while and then smaller the more times the loop runs. I should be getting better accuracy with more iterations. The answer seems to hover around an accurate value until I run a lot of iterations and then it plummets.

I apologize for putting a FORTRAN question in a VB forum but the FORTRAN forums I could find weren't very good and I figure most math/engineering people have had a lot of experience with FORTRAN. Thanks for any help you can give.

Here is the vb code and then the FORTRAN code for the same program:

Sub CalcIntegral()


Dim i As Long
Dim dxval As Double, Lower As Double, Upper As Double
Dim dx As Double, x As Double, y As Double

'Set upper and lower limits for integration
Upper = 1
Lower = 0

'dxval represents the number of segments used in the integration
dxval = 1000000000

'dx represents the width of each segment
dx = (Upper - Lower) / dxval
x = Lower

'A For…Next loop calculates the area of each segment and adds up all the areas as a Riemann sum
For i = 1 To dxval

y = y + (x ^ 2 * dx)
x = x + dx

Next i

msgbox y

End Sub

NOW THE FORTRAN VERSION---------------------------------

PROGRAM ERIC
IMPLICIT NONE

real,parameter :: lower=0.0,upper=1.0
real :: dx,x,y,dxval
integer :: i

dxval=1000000000
y=0.0
dx=(upper-lower)/dxval
x=0.0

DO i = 1,dxval, 1
y = y + ((x**2)*dx)
x = x + dx

END DO

write(*,*) y


END PROGRAM

[oswaler]

After posting the last message I tried translating another program from VB to FORTRAN and got the same sort of bizarre result. This program finds the numeric derivative for the function C(x)=200.0 + 0.05 * x + 0.0001 * x^2 where x=500. I used the definition lim as h approaches 0 f'(x)=(f(x+h)-f(x))/h. In vb this works fine and gives the proper answer of 0.15. The smaller my h value, the better the accuracy of the result.

When I put it in FORTRAN I find 2 strange things:

1) If I make the h value less than 1/21846 then I get an answer of 0.0000.

2) For h value of up to around 1/10000 I get an answer of 0.152588 which is close to the actual answer. Once the h value decreases from there (1/11000,1/12000, etc.) the answer gets higher and higher until it maxes out at about .33 or so when the h value is 1/21846.

As in the example in the previous post, the answers are approximately correct until I start asking for more precision, and then the answers diverge wildly from the true answer.

Here's the VB and FORTRAN code for this one. Again, thanks for any help anyone can give.

Sub calcderiv()

Dim origx As Double
Dim h As Double
Dim df As Double
Dim f As Double
Dim x As Double
Dim fprime As Double

x=500

'Define h as a very small value
h = 1 / 2109999999

'Hold original value of x to be used later
origx = x

'set x value to x plus a very small amount
x = x + h

'calculate the function with the slightly increased value of x
df = 200 + 0.05 * x + 0.0001 * x ^ 2

'set x back to original value
x = origx

'calculate the function with the original value of x
f = 200 + 0.05 * x + 0.0001 * x ^ 2

‘find the numeric derivative
fprime = (df - f) / h
MsgBox fprime

End Sub

THE FORTRAN CODE --------------------------------

PROGRAM DERIV
IMPLICIT NONE


REAL :: origx,h,df,f,x,fprime

x = 500.0

h=1/1000000.0
origx = x

x = x + h

df = 200.0 + 0.05 * x + 0.0001 * x**2

x = origx

f = 200.0 + 0.05 * x + 0.0001 * x**2

fprime = (df - f) / h
write(*,*) fprime

END PROGRAM

[oswaler]

In case anyone cares, I found the solution. The x and y variable in the integration program needed to be set to double precision instead of real.

-Eric