1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

FORTRN90: Euler Midpoint Method for SHO

  1. Dec 3, 2014 #1
    1. The problem statement, all variables and given/known data
    Write a program to simulate motion of simple harmonic oscillator.
    Initial conditions: Let ω = 1, x(t=0) = 1, v(t=0) = 0.
    Integrate over 30 seconds in intervals of 0.05s.

    2. Relevant equations
    δ2x / δt2 = -ω2x
    As set of 2 coupled ODE's; x' = v, v' = -w2x

    3. The attempt at a solution

    When I printed out the results, i noticed I was getting larger and larger ith numbers of v each time I ran. I am obtaining a graph that represents SHO for the first half a wavelength. Basically, I'm not sure where my error is coming from. I've not included the DISLIN part of the script as this isn;t where the problem is arising from.

    PROGRAM
    eulershm
    IMPLICIT NONE
    INTEGER, PARAMETER
    :: n = 600
    INTEGER :: i = 1
    REAL, DIMENSION(n):: x, v, t
    REAL, PARAMETER :: dt = 0.05, tmax = 30, w = 1
    v = 0
    x = 1
    DO WHILE (t(i).LE.tmax)
    CALL evalF(v(i), v(i+1), x(i), x(i+1), t(i), t(i+1), w, dt)
    i = i+1
    END DO
    END PROGRAM
    eulershm

    SUBROUTINE evalF(v1, v2, x1, x2, t1, t2, w1, tau)
    IMPLICIT NONE
    REAL
    : v1, v2, x1, x2, t1, t2
    REAL :: w1, tau
    x2 = x1 + v1*tau
    v2 = v1 - (w1**2)*tau
    t2 = t1 + tau
    END SUBROUTINE evalF




    Many thanks to whoever can help.
     
  2. jcsd
  3. Dec 3, 2014 #2

    Mark44

    Staff: Mentor

    Fortran is a call-by-reference language, which means that the actual arguments to a subroutine can be changed if the formal parameters are changed in the subroutine definition. What I mean by this is that in your main program, you call evalF with parameters v(i) and v(i + 1) and others. These are what I'm referring to as "actual" arguments. The corresponding formal parameters of this subroutine are v1, v2 and others. Inside your sub you reset x2, v2, and t2, which causes the actual arguments v(i + 1), x(i + 1), and t(i + 1) to change in the main program. I suspect that is what you're seeing.

    Based on the name you chose for your subroutine, evalF, a better choice might be to write it as a function that returns a function value.

    BTW, and this is a minor point, your code wouldn't be considered a "script."
     
  4. Dec 3, 2014 #3
    Hi Mark,
    When I bring in the actual arguments , the (i+1) arguments are empty array spaces, and my method has been to simply to overwrite that empty space with the value obtained from the (i) arguments. Then the do loop will simply make what was the (i+1) the (i) and repeat for the next empty space, terminating when t(i) = tmax.
    Point is, I'm not sure what your identifying as the error source because isn;t this what you just described?

    Do you mean a better choice is to call a function instead of a subroutine?
    And lastly, what defines a "script"?
     
  5. Dec 3, 2014 #4

    Mark44

    Staff: Mentor

    It's been a long while since I wrote any Fortran code, so it's not clear to me what these statements do.
    Code (Text):
    v = 0
    x = 1
    Both v and x are arrays, so I guess the intent is to initialize all elements of the arrays to 0 and 1, respectively. A better way to do things, IMO, is just to initialize the first elements of the two arrays like this:
    Code (Text):
    v(1) = 0
    x(1) = 1
    I don't see that you have initialized your t array anywhere, at least in the code you provided. That will cause problems in this line:
    Code (Text):
    [B]DO WHILE[/B] (t(i)[B].LE.[/B]tmax)
    You should never evaluate a variable before you have given it a value. It will have a value, but its value will be whatever "garbage" happens to be lying around in that memory address.

    Also, and this is something I see a lot in Fortran code, use some white space. The compiler doesn't care at all, but packing everything in with no spaces between things makes it a lot harder for humans to read. For example, the conditional expression in your while statement would be easier to read as
    Code (Text):
    t(i) .LE. tmax
    ; i.e., with spaces before and after the .LE. operator.

    As mentioned already, you don't have any comments, so it's not obvious what your subroutine is supposed to do. Since its job apparently is to fill in three arrays, it shouldn't be a function, contrary to what I said earlier. What I was going by was the name of the routine, evalF, which suggests (inaccurately) that it is evaluating a function. IOW, its name is misleading, which again makes it harder for humans to discern what you are trying to do. To make it easier on us poor humans, give names to your routines that reflect what it is they do.

    A "script" would be a piece of code that is processed by a scripting language such as JavaScript or LaTeX. Scripting languages are generally interpreted languages that are processed on the fly, as opposed to compiled languages in which the code is translated to object code all at once, and then the object code is executed.
     
  6. Dec 3, 2014 #5
    Hi Mark,
    Thanks there is a lot of helpful information there. As is often the case, seems that its the minor errors that are preventing progression. And the non-initialisation of t (as you mentioned) is what is most likely causing the ith values to keep ascending.
    I tried running after going through, and I am reading the following;

    note - I wrote % instead of a hash symbol because my laptop doesn't have one!!

    Program received signal SIGBUS: Access to an undefined portion of a memory object.
    Backtrace for this error: %0x7FB26DEFD7D7
    %0x7FB26DEFDDDE
    %0x7FB26DB54C2F

    /usr/local/bin/gf95link: line 104: 3219 Bus error (core dumped)./$bname

    Bus errors are associated with transfer of data to and from main and subroutine, but I don't see the problem with the original code. Any idea?
     
    Last edited: Dec 3, 2014
  7. Dec 3, 2014 #6

    DrClaude

    User Avatar

    Staff: Mentor

    I don't think the problem comes from the initialization of t. You are not using the Euler midpoint method, but the simpler Euler method, which is unstable.
     
  8. Dec 3, 2014 #7
    Of course, a silly mistake not to mention I am current performing simple before applying midpoint. Do you think then upon advancement to midpoint this problem should be resolved?
     
  9. Dec 3, 2014 #8

    Mark44

    Staff: Mentor

    I suspect that in one call to your subroutine you have passed in an array element that you didn't declare. You have declared your v, x, and t arrays to have 600 elements. I'm guessing that in the call that generated the error above, the index was 601.

    Since you know the number of elements in each array, the most reasonable thing to use would be a counting loop:
    Code (Text):
    DO I = 1, N
    , rather than the while loop with the condition you have.
     
  10. Dec 4, 2014 #9

    DrClaude

    User Avatar

    Staff: Mentor

    I looked in more details, and I found the main error in your program. Have you really implemented this set of equations?

    Thats said, I checked and the Euler method with w = 1 and dt = 0.05 is unstable: the oscillations grow with time. You will still be able to see some oscillations, provided you implement your ODEs correctly.

    Also, you should follow Mark's advice.
     
  11. Dec 4, 2014 #10
    I did consider this previously, went with the while so that the termination parameter could be simply changed but I suppose neither is easier than the other.

    Once again I am at a loss, I thought that the vn+1 and xn+1 equations were sufficient, since these give me values of v and x at all intervals and w is a parameter. Im obviously getting something completely wrong, sorry about this. I think I need to go back to the beginning with it to understand where I am going wrong.
     
  12. Dec 4, 2014 #11

    DrClaude

    User Avatar

    Staff: Mentor

    You should have the program calculate the number of iterations and use that.

    My point is that the program does not calculate the set of first-orer ODEs you wrote. You made a mistake when going from paper to program.
     
  13. Dec 4, 2014 #12
    I still don't understand your implication - are you saying that my program doesn't but should? All I have done is put them in there with the relevant equations because on my paper work, they were relevant in getting me to the v(i+1) and x(i+1) values. Could you go from the basics with what your trying to say please?
     
  14. Dec 4, 2014 #13

    DrClaude

    User Avatar

    Staff: Mentor

    These two lines of code
    do not correspond to these ODEs
    One of the two lines of code is in error.

     
  15. Dec 4, 2014 #14
    Im missing the factor of x in equation v2 aren't I...
     
  16. Dec 4, 2014 #15

    DrClaude

    User Avatar

    Staff: Mentor

    Indeed. :)
     
  17. Dec 4, 2014 #16
    So sorry for taking up your time with such a minor error! Hopefully I won't be back on this thread now. Thanks though DrClaude and Mark!
     
  18. Dec 4, 2014 #17

    DrClaude

    User Avatar

    Staff: Mentor

    Can't talk for everyone, but I spend time on PF to learn and for the pleasure of helping others, be it for small and big things! Don't hesitate to ask more questions.
     
  19. Dec 4, 2014 #18
    Thanks again DrClaude! Good to know!
     
  20. Dec 4, 2014 #19

    Mark44

    Staff: Mentor

    I second what DrClaude said. All of the homework helpers, science advisors, and mentors here are unpaid. Helping others is a big motivation for us.
     
  21. Dec 4, 2014 #20
    It's good to know that there are still vocational physicists rather than careerists.
    Anyhow, I have implemented the advice and it has definitely gotten me a lot further than before. Im still focusing on simple euler, but I have a couple of clear problems still.
    The outputted data, and consequently the graph of v vs. t (as well as with x) is not the same upon each run. And on ANY of the runs, I am not observing anything remotely like simple harmonic motion.
    I'll put here the revised code and let you take another look (sorry I have to type it all out as cant copy and paste over from virtual machine).
    Also, if you could tell me how to simply put both the velocity and position curves on a single graph that would be helpful for future reference.

    PROGRAM eulershm
    USE DISLIN
    IMPLICIT NONE
    INTEGER, PARAMETER ::
    n = 600
    INTEGER:: i
    REAL, DIMENSION (n) :: x, v, t
    x(i) = 1 !Initialization
    v(i) = 0
    t(i) = 0
    DO i = 1,n
    CALL evalS(v(i), v(i+1), x(i), x(i+1), t(i), t(i+1))
    PRINT*, "v(",i,")=", v(i) !Print out values of v (for clarification)
    END DO

    CALL
    METAFL('XWIN')
    CALL DISINI()
    CALL GRAF(0.,30.,0.,1.,-1.,1.,-1.,0.1)
    CALL COLOR('RED')
    CALL CURVE(t,v,i)
    DO i = 1,n
    CALL RLSYMB 93, t(i), v(i))
    END DO

    END PROGRAM eulershm

    SUBROUTINE evalS(v1, v2, x1, x2, t1, t2)
    IMPLICIT NONE
    REAL, INTENT(IN)::v1,x1,t1
    REAL, INTENT(OUT)::v2, x2, t2
    REAL, PARAMETER:: w = 1, dt = 0.05
    v2 = v1 - ((w**2)*x1)*dt
    x2 = x1 + v1*dt
    t2 = t1 + dt
    END SUBROUTINE evalS
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: FORTRN90: Euler Midpoint Method for SHO
  1. Euler's Method (Replies: 0)

Loading...