Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Euler method in Fortran

  1. Nov 12, 2013 #1
    Hi everybody,
    I am programming a new code for a problem.
    The problem is numerically solving the Simple Harmonic Motion using the Euler method. This approach is just an approximate solution and not a exact solution, however when I run the code successfully and plot my data, it comes up as an exact solution and the plot is exactly the same as sinusoidal harmonic motion. (while I am expecting to detect error!)
    I'm really not sure where I am going wrong here. Any help would be appreciated. Thank you!

    Here is my code:

    program SHM_Euler

    implicit none
    real, dimension(100) :: x, v, t
    real, dimension(100) :: a
    integer :: i
    real :: dt !(Step Size)

    x(1) = 0.01745
    t(1) = 0.00
    dt = 0.1

    write (*,*) ' t x(t) v(t) '

    do i=1, 100
    t(i+1) = t(i) + dt

    v(i) = -0.0546*sin(3.1305*t(i)) !Velocity
    a(i)=-0.1710*cos(3.1305*t(i)) !Acceleration

    x(i+1) = x(i) + v(i)*dt
    v(i+1) = v(i) + a(i)*dt

    write(6,'(3(e12.5,3x))') t(i+1), x(i+1), v(i+1)

    t(i) = t(i+1)
    x(i) = x(i+1)

    end do

    end program SHM_Euler
  2. jcsd
  3. Nov 12, 2013 #2


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    It would help if you described the details of the SHM problem.

    Where did the expressions for v(i) and a(i) come from?

    I have a problem with the indexing of your t, x, v, and a variables. Using the i+1 index means that when i = 100, you are calculating variables of index = 101, when the maximum dimension is only 100. This is careless programming.

    It's not necessarily outside the realm of possibility that a numerical method like Euler could give exact solutions to a certain differential equation.
  4. Nov 12, 2013 #3


    User Avatar

    Staff: Mentor

    It seems to me that you are assuming a sinusoidal solution here. You should start with the differential equation for SHM which gives you

    a(i) = -c*x(i)

    where c is a constant that characterizes the system, e.g. for a mass on a spring c = k/m where k is the spring constant and m is the mass. Then find v(i+1) and x(i+1) as you've done.

    Also, these statements should not be there:

    They replace the current values of t, x, and v with the new ones that you have just calculated. The net effect is to shift each data point "down" one position in the arrays, wiping out the initial values that you set at the beginning of the program.

    And as SteamKing noted, don't try to go past the ends of your arrays!
    Last edited: Nov 12, 2013
  5. Nov 13, 2013 #4
    Simple Harmonic Motion (SHM) is a sinusoidal motion in time and demonstrates a single resonant frequency. In this equation position is x=A Cos(wt) which A is an amplitude and w is a frequency. Similarly, we can derive velocity and acceleration by differentiating from x (v=dx/dt & a=dv/dt). I put an amount for amplitude and frequency.

    You are right. I have changed index to real, dimension(0:100) :: x, v, t and t(i) = t(i-1) + dt and it was run successful but I still have the same problem.
  6. Nov 13, 2013 #5
    The replace statement is a part of Euler method that substitute an old value to a new one.
  7. Nov 13, 2013 #6


    User Avatar

    Staff: Mentor

    You need to "replace old values with new ones" if you are storing them (one after the other) in the same variables, without using arrays, so you can advance from one step to the next. However, you're using arrays which contain the results for all steps, and you advance from one step to the next by incrementing the array index.

    It's the difference between doing something like this for a trivial constant-velocity problem (without using an array):

    Code (Text):

    x = 0
    t = 0
    dt = 0.1
    v = 5

    do i = 1, 100
        write (*, *) t, x
        t_next = t + dt
        x_next = x + v * dt
        t = t_next
        x = x_next
    end do
    versus this (using an array):

    Code (Text):

    x(1) = 0
    t(1) = 0
    dt = 0.1
    v = 5

    do i = 1, 100
        write (*, *) t(i), x(i)
        t(i+1) = t(i) + dt
        x(i+1) = x(i) + v * dt
    end do
    The practical difference between the two versions is that at the end of the first one, you have only the final values of x and t available for future calculations; but at the end of the second one, you have all the values for the intermediate steps available in the arrays.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook