FORTRAN: second-order ODE with Euler Method

[Cecily]

Homework Statement

Dear all, please help. I have tried this question and came up with strange numbers, my fortran is definitely not correct. Please help!

When the effect of the air resistance is taken into account, the equation of motion for a particle of mass
m falling vertically in a gravitational field is given by
m*y''= -m*g - b*v (*)
where g is the gravitational acceleration and v = dy/dt is the velocity. The coefficient b specifies the strength of the interaction between the air and the particle. The initial conditions are y(0) = 1000 m and v(0) = 0 ms^-2. Take m = 1 kg and g = 9.8 ms^-2 in the following.Write a Fortran program that uses the Euler method to solve the initial value problem.
Consider the case without air resistance b = 0. Plot y(t) against t until y = 0 (i.e., the particle has
reached the ground) for two different time steps Δt = 1 s and 0.01 s in one figure. Plot also the
exact solution in the same figure for comparison.

Homework Equations

The Attempt at a Solution

 

[Cecily]

here is my definitely wrong program...

PROGRAM lab5a
IMPLICIT NONE

INTEGER :: i
REAL :: dt, y_0, v_0, f
REAL, DIMENSION(200) :: y, v, t
REAL, PARAMETER :: epsilon = 1D-6

y(1) = y_0
v(1) = v_0

dt = 1.D0
DO i=2
y(i) = y(1) + v(i-1)*dt
v(i) = v(i-1) + dt*f(y(i-1), v(i-1))
t(i) = t(i-1) + dt

WRITE (10, *) t(i), y(i)

IF (y(i) < epsilon) exit

END DO

dt = 1D-2

DO i=1
y(i) = y(1) + v(i-1)*dt
v(i) = v(i-1) + dt*f(y(i-1), v(i-1))
t(i) = t(i-1) + dt

WRITE (11, *) t(i), y(i)

IF (y(i) < epsilon) exit

END DO

END PROGRAM lab5a

!============================
REAL FUNCTION f(y, v)
IMPLICIT NONE
REAL :: y, v
REAL :: m=1, g=9.8, b=0

f = -g - b / m * v

END FUNCTION

 

[Mark44]

here is my definitely wrong program...

PROGRAM lab5a
IMPLICIT NONE

INTEGER :: i
REAL :: dt, y_0, v_0, f
REAL, DIMENSION(200) :: y, v, t
REAL, PARAMETER :: epsilon = 1D-6

y(1) = y_0
v(1) = v_0

dt = 1.D0
DO i=2
y(i) = y(1) + v(i-1)*dt
v(i) = v(i-1) + dt*f(y(i-1), v(i-1))
t(i) = t(i-1) + dt

WRITE (10, *) t(i), y(i)

IF (y(i) < epsilon) exit

END DO

dt = 1D-2

DO i=1
y(i) = y(1) + v(i-1)*dt
v(i) = v(i-1) + dt*f(y(i-1), v(i-1))
t(i) = t(i-1) + dt

WRITE (11, *) t(i), y(i)

IF (y(i) < epsilon) exit

END DO

END PROGRAM lab5a

!============================
REAL FUNCTION f(y, v)
IMPLICIT NONE
REAL :: y, v
REAL :: m=1, g=9.8, b=0

f = -g - b / m * v

END FUNCTION

Your DO loops are not right. For the first one, you have
DO i=1
.
.
.
END DO

You need to indicate the final value of your loop control variable, like this:
DO i=1 , 10
.
.
.
END DO

Your other DO loop has the same problem.