FORTRN90: Euler Midpoint Method for SHO

In summary, the conversation discusses a program to simulate motion of a simple harmonic oscillator with initial conditions of ω = 1, x(t=0) = 1, v(t=0) = 0 and integration over 30 seconds in intervals of 0.05s. The equations used are δ2x / δt2 = -ω2x and a set of 2 coupled ODE's; x' = v, v' = -w2x. The code includes a subroutine called evalF, which may cause errors due to the call-by-reference nature of Fortran. A suggestion is to use a function instead of a subroutine, as well as clarification on what defines a "script."
  • #1
SalfordPhysics
69
1

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


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

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.
 
Physics news on Phys.org
  • #2
SalfordPhysics said:

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


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

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."
 
  • #3
Mark44 said:
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"?
 
  • #4
SalfordPhysics said:
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.
Code:
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:
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:
[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:
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.
 
  • Like
Likes SalfordPhysics
  • #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:
  • #6
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.
 
  • #7
DrClaude said:
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?
 
  • #8
SalfordPhysics said:
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:
Code:
DO I = 1, N
, rather than the while loop with the condition you have.
 
  • Like
Likes SalfordPhysics and DrClaude
  • #9
I looked in more details, and I found the main error in your program. Have you really implemented this set of equations?

SalfordPhysics said:
As set of 2 coupled ODE's; x' = v, v' = -w2x

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.
 
  • Like
Likes SalfordPhysics
  • #10
Mark44 said:
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:
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.

DrClaude said:
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 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. 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.
 
  • #11
SalfordPhysics said:
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.

SalfordPhysics said:
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.
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.
 
  • #12
DrClaude said:
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?
 
  • #13
These two lines of code
SalfordPhysics said:
x2 = x1 + v1*tau
v2 = v1 - (w1**2)*tau
do not correspond to these ODEs
SalfordPhysics said:
As set of 2 coupled ODE's; x' = v, v' = -w2x
One of the two lines of code is in error.

 
  • #14
Im missing the factor of x in equation v2 aren't I...
 
  • #15
SalfordPhysics said:
Im missing the factor of x in equation v2 aren't I...
Indeed. :)
 
  • Like
Likes SalfordPhysics
  • #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!
 
  • #17
SalfordPhysics said:
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.
 
  • #18
Thanks again DrClaude! Good to know!
 
  • #19
DrClaude said:
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.
 
  • #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. 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
 
  • #21
SalfordPhysics said:
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.
Code:
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
Code:
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)
 
  • #22
Mark44 said:
I see a few problems, some of them possibly due to typing transcription errors.
1. Your initializations aren't right.
Code:
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
Code:
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!
 

Attachments

  • eulerSHO.png
    eulerSHO.png
    34 KB · Views: 523
  • #23
SalfordPhysics said:
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.
 
  • Like
Likes SalfordPhysics
  • #24
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).
 
  • #25
A better solution, IMO, is to declare the array so that it starts with index 0.
Code:
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.
 
  • #26
Edit: I was reply to an old post without noticing there were newer ones.
 

Related to FORTRN90: Euler Midpoint Method for SHO

1. What is FORTRN90 and how is it used in the Euler Midpoint Method for SHO?

FORTRAN90 is a high-level programming language commonly used for scientific computing. It is often used to implement numerical methods, such as the Euler Midpoint Method, for solving differential equations. In the Euler Midpoint Method for Simple Harmonic Oscillators (SHO), FORTRAN90 is used to numerically approximate the solution to the SHO differential equation by dividing the time interval into smaller steps and calculating the position and velocity of the oscillator at each step.

2. What is the Euler Midpoint Method and how does it work?

The Euler Midpoint Method is a numerical method used to approximate the solution to a differential equation. It works by dividing the time interval into smaller steps and using the midpoint of each step to calculate the values of the solution. In the case of the SHO, the Euler Midpoint Method calculates the position and velocity of the oscillator at each step, which can then be used to approximate the solution to the SHO differential equation.

3. What are the advantages of using the Euler Midpoint Method for SHO?

There are several advantages to using the Euler Midpoint Method for SHO. Firstly, it is a simple and easy-to-implement method, making it a popular choice for solving differential equations in scientific computing. Additionally, the method is relatively accurate and stable, providing a good approximation of the true solution to the SHO differential equation.

4. What are some limitations of the Euler Midpoint Method for SHO?

Like any numerical method, the Euler Midpoint Method for SHO has its limitations. One of the main limitations is that it is a first-order method, meaning that the error in the approximation will decrease at a rate proportional to the step size. This can result in a less accurate solution compared to higher-order methods. Additionally, the method may not be suitable for all types of differential equations and may require a smaller step size to achieve a desired level of accuracy.

5. How is the accuracy of the Euler Midpoint Method for SHO evaluated?

The accuracy of the Euler Midpoint Method for SHO can be evaluated by comparing the approximate solution to the true solution of the SHO differential equation. This can be done by varying the step size and observing how the error changes. Additionally, the method can be compared to other numerical methods to determine its accuracy and suitability for different types of SHO problems.

Similar threads

  • Engineering and Comp Sci Homework Help
Replies
6
Views
4K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
2K
Replies
10
Views
711
  • Engineering and Comp Sci Homework Help
Replies
4
Views
1K
  • Differential Equations
Replies
1
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
5
Views
1K
  • Calculus and Beyond Homework Help
Replies
1
Views
1K
  • Introductory Physics Homework Help
Replies
17
Views
6K
  • Engineering and Comp Sci Homework Help
Replies
10
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
18
Views
2K
Back
Top