FORTRN90: Euler Midpoint Method for SHO

[SalfordPhysics]

Homework Statement

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.

Homework Equations

δ²x / δt² = -ω²x
As set of 2 coupled ODE's; x' = v, v' = -w²x

The Attempt at a Solution

[/B]
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.

 

[Mark44]

Homework Statement

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.

Homework Equations

δ²x / δt² = -ω²x
As set of 2 coupled ODE's; x' = v, v' = -w²x

The Attempt at a Solution

[/B]
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.

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."

 

[SalfordPhysics]

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."

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"?

 

[Mark44]

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"?

It's been a long while since I wrote any Fortran code, so it's not clear to me what these statements do.

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:

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:

[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

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.

 

[SalfordPhysics]

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?

 

[DrClaude]

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.

 

[SalfordPhysics]

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.

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?

 

[Mark44]

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?

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:

DO I = 1, N

, rather than the while loop with the condition you have.

 

[DrClaude]

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

As set of 2 coupled ODE's; x' = v, v' = -w²x

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.

 

[SalfordPhysics]

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:

DO I = 1, N

, rather than the while loop with the condition you have.

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.

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.

Once again I am at a loss, I thought that the v_n+1 and x_n+1 equations were sufficient, since these give me values of v and x at all intervals and w is a parameter. I am 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.

 

[DrClaude]

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.

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

Once again I am at a loss, I thought that the v_n+1 and x_n+1 equations were sufficient, since these give me values of v and x at all intervals and w is a parameter.

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.

 

[SalfordPhysics]

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.

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?

 

[DrClaude]

These two lines of code

x2 = x1 + v1*tau
v2 = v1 - (w1**2)*tau

do not correspond to these ODEs

As set of 2 coupled ODE's; x' = v, v' = -w²x

One of the two lines of code is in error.

 

[SalfordPhysics]

Im missing the factor of x in equation v2 aren't I...

 

[DrClaude]

Im missing the factor of x in equation v2 aren't I...

Indeed. :)

 

[SalfordPhysics]

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!

 

[DrClaude]

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!

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.

 

[SalfordPhysics]

Thanks again DrClaude! Good to know!

 

[Mark44]

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.

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.

 

[SalfordPhysics]

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. I am 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 can't 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

 

[Mark44]

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. I am 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 can't 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

I see a few problems, some of them possibly due to typing transcription errors.
1. Your initializations aren't right.

x(i) = 1                                        !Initialization
v(i) = 0
t(i) = 0

Since i is uninitialized, the three statements above are setting a random element of each array to 1 or 0.
The above should be

x(1) = 1                                        !Initialization
v(1) = 0
t(1) = 0

You should never use a variable before it's been initialized, other than to set it to a value.
2. Your DO loop runs from i = 1 to 600, inclusive, and your arrays are declared to have elements with indexes from 1 through 600. In your last call to evalS, you are passing v(600), v(601), x(600), x(601), t(600), and t(601). Your code then attempts to set values for memory that's outside the bounds of the v, x, and t arrays. Not good.
3. This call is missing its parentheses: CALL RLSYMB 93, t(i), v(i))
4. Add some spaces to make this easier to read: CALL GRAF(0.,30.,0.,1.,-1.,1.,-1.,0.1)
Here it is with a space after each comma: CALL GRAF(0., 30., 0., 1.,-1., 1.,-1., 0.1)

 

[SalfordPhysics]

I see a few problems, some of them possibly due to typing transcription errors.
1. Your initializations aren't right.

x(i) = 1                                        !Initialization
v(i) = 0
t(i) = 0

Since i is uninitialized, the three statements above are setting a random element of each array to 1 or 0.
The above should be

x(1) = 1                                        !Initialization
v(1) = 0
t(1) = 0

You should never use a variable before it's been initialized, other than to set it to a value.
2. Your DO loop runs from i = 1 to 600, inclusive, and your arrays are declared to have elements with indexes from 1 through 600. In your last call to evalS, you are passing v(600), v(601), x(600), x(601), t(600), and t(601). Your code then attempts to set values for memory that's outside the bounds of the v, x, and t arrays. Not good.
3. This call is missing its parentheses: CALL RLSYMB 93, t(i), v(i))
4. Add some spaces to make this easier to read: CALL GRAF(0.,30.,0.,1.,-1.,1.,-1.,0.1)
Here it is with a space after each comma: CALL GRAF(0., 30., 0., 1.,-1., 1.,-1., 0.1)

Thanks Mark
1) I had previously declared i = 1, but removed this and forgot to change my initialisation.
2) This is the thing that has caught me out now numerous times. Should I set then the array size (n) as 601, and then run from i = 1 to (n-1)?
3) Thats a mistake from re-typing it all out, sorry!

Having made these changes, finally obtained a function close to simple harmonic oscillation that increases in amplitude over time (as consequential error of the program).
I have uploaded a screenshot so you can admire your great contributions!

 

[Mark44]

2) This is the thing that has caught me out now numerous times. Should I set then the array size (n) as 601, and then run from i = 1 to (n-1)?

I would leave the array size alone, but run the loop from 1 to N - 1. On the final iteration (when i = 599), you'll be reading v(599) and setting v(600). Same for the other two arrays.

 

[SalfordPhysics]

Thanks Mark I would go with that but I need to simulate up to t=30 seconds, and since i=1 is for t=0, i=2 is first point making i=601 the 600th interval (the 30second point).

 

[Mark44]

A better solution, IMO, is to declare the array so that it starts with index 0.

REAL, DIMENSION (0: n) :: x, v, t

With n as a parameter (600), the array above will have 601 elements. By default, arrays in Fortran start with index 1, but you can have them start at other indexes, as in the example above.

 

[DrClaude]

Edit: I was reply to an old post without noticing there were newer ones.