How to fix the Euler method error?

[koushan]

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
v(1)=0

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)
v(i)=v(i+1)

end do

end program SHM_Euler

 

[SteamKing]

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.

 

[jtbell]

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

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:

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

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!

 

[koushan]

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.

 

[koushan]

The replace statement is a part of Euler method that substitute an old value to a new one.

 

[jtbell]

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):

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):

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.