Fortran 90. Euler and Runge-Kutta's method

[fluidistic]

I've successfully (well I think so) made 2 different programs that can numerically solve an ODE using Euler and Runge-Kutta's methods.
Here they are:
program test
implicit none
real(8)::a,b,h,y_0,t

write(*,*)"Enter the interval a,b, the value of the step-size h and the value of y_0"
read(*,*)a,b,h,y_0

open(unit=1, file="resultadoEuler.dat",status='old')
call euler(y_0,a,b,h)
close(1)

contains
subroutine euler(y_0,a,b,h)
implicit none
real(8), intent(inout)::a,h,b, y_0
real(8):: y,t
t=a
y=y_0
do while(t<=b)
y=y+h*f(y,t)
write(1,*)t, y
t=t+h
end do
end subroutine
function f(y,t)
implicit none
real(8)::f,y,t,pi
pi=acos(-1d0)
f=-y+sin(pi*2d0*t)
end function
end program
and program solvingwithRungeKutta
implicit none
real(8)::a,b,h,y_0,t

write(*,*)"Enter the interval a,b, the value of the step-size h and the value of y_0"
read(*,*)a,b,h,y_0
open(unit=1, file="resultadoRungeKutta.dat",status='old')
call RungeKutta(y_0,a,b,h)
close(1)

contains
subroutine RungeKutta(y_0,a,b,h)
implicit none
real(8), intent(inout)::a,h,b, y_0
real(8):: y,t, F_1, F_2, F_3, F_4
t=a
y=y_0
do while(t<=b)

F_1=h*f(y,t)
F_2=h*f(y+(1d0/2d0)*F_1, t+h*(1d0/2d0))
F_3=h*f(y+(1d0/2d0)*F_2, t+h*(1d0/2d0))
F_4=h*f(y+F_3,t+h)

y=y+(F_1+2d0*F_2+2d0*F_3+F_4)/6d0
write(1,*)t, y
t=t+h
end do
end subroutine
function f(y,t)
implicit none
real(8)::f,y,t,pi
pi=acos(-1d0)
f=-y+sin(pi*2d0*t)
end function
end program

I used these codes to solve for y(t) in y'(t)=-y(t)+\sin (2 \pi t) with the initial condition y(0)=1 in the interval [0,1]. I entered various step-sizes h such as 0.1, 0.01 and 0.001.
I also have the exact solution to this ODE so that I can calculate the error for Euler and Runge-Kutta's methods. And this is where I realize there's an error somewhere. I get an error about the same value for both methods.
This happens when I made another code that combine both codes.
Here it is:
program test2
implicit none
real(8)::a,b,h,y_0,t,y, pi

write(*,*)"Enter the interval a,b, the value of the step-size h and the value of y_0"
read(*,*)a,b,h,y_0


open(unit=32, file="valores_exactos.dat")
pi=acos(-1d0)
t=a
do while (t<=b)
y=(1d0+(2d0*pi)/(1d0+4d0*pi**2d0))*exp(-t)+(sin(2d0*pi*t)-2d0*pi*cos(2d0*pi*t))/(1d0+4d0*pi**2d0)
write(32,*)t, y
t=t+h
end do
close(32)

open(unit=1, file="resultadoEuler.dat",status='old')
call euler(y_0,a,b,h)
close(1)
open(unit=99, file="resultadoRungeKutta.dat",status='old')
call RungeKutta(y_0,a,b,h)
close(99)

contains
subroutine euler(y_0,a,b,h)
implicit none
real(8), intent(inout)::a,h,b, y_0
real(8):: y,t, valor_exacto, pi, error
pi=acos(-1d0)
t=a
y=y_0
open(unit=2, file="error_Euler.dat")
do while(t<=b)
y=y+h*f(y,t)
valor_exacto=(1d0+(2d0*pi)/(1d0+4d0*pi**2d0))*exp(-t)+(sin(2d0*pi*t)-2d0*pi*cos(2d0*pi*t))/(1d0+4d0*pi**2d0)
error=abs(y-valor_exacto)
write(1,*)t, y, valor_exacto
write(2,*)error
t=t+h
end do
close(2)
end subroutine

subroutine RungeKutta(y_0,a,b,h)
implicit none
real(8), intent(inout)::a,h,b, y_0
real(8):: y,t, F_1, F_2, F_3, F_4, valor_perfecto, pi, error_kutta
pi=acos(-1d0)
t=a
y=y_0
open(unit=100, file="error_RungeKutta.dat")
do while(t<=b)
F_1=h*f(y,t)
F_2=h*f(y+(1d0/2d0)*F_1, t+h*(1d0/2d0))
F_3=h*f(y+(1d0/2d0)*F_2, t+h*(1d0/2d0))
F_4=h*f(y+F_3,t+h)
y=y+(F_1+2d0*F_2+2d0*F_3+F_4)/6d0
valor_perfecto=(1d0+(2d0*pi)/(1d0+4d0*pi**2d0))*exp(-t)+(sin(2d0*pi*t)-2d0*pi*cos(2d0*pi*t))/(1d0+4d0*pi**2d0)
error_kutta=abs(y-valor_perfecto)
write(1,*)t, y, valor_perfecto
write(100,*)error_kutta
t=t+h
end do
close(100)
end subroutine
function f(y,t)
implicit none
real(8)::f,y,t,pi
pi=acos(-1d0)
f=-y+sin(pi*2d0*t)
end function
end program

So my 2 first codes seem to work out properly. (I plotted both results along with the exact value of y(t) for 1 thousands points and Runge-Kutta was like approximately 5 times closer to the exact values than Euler's method was, using only my eyes).
However here is the graph I get if I plot the error of both algorithm in code #3. (See attachment).
So it seems my third code doesn't work as my first 2 codes do. I have absolutely NO IDEA what's going on. I took exactly the same subroutines, copied and pasted. I reset every values for both subroutines so that I don't modify any value for the next subroutine.
I'd love some help here, I'm really lost.

[rcgldr]

What does the graph of the non-combined routines look like?

[fluidistic]

What does the graph of the non-combined routines look like?
I don't have the errors in my first codes. However I do have the exact values versus the calculated value for each methods.
I send 2 pictures. One more or less general and the other around where I should have noticed the biggest value of error for Runge-Kutta's method.

This is actually convincing myself that Euler's method doesn't seem that bad at all, even with a small h.
I'll try right now with a small h (around 4 or 5 divisions) and post some pics.

[fluidistic]

Here with h=0.2.

[rcgldr]

I'm not sure what the issue with the combined programs is, but if you want to try another method, the link below is to a description of an iterative (predictor + corrector) method. You could try doing 4 iterations (the initial Euler + 3 trapezoidals, (y[3])) to compare with Runge-Kutta with about the same overhead. With more iterations, this method will converge to specific values for each step, but since it's a linear approximation for each step (trapezoidal), you still need to keep h relatively small.

Euler_trapezoidal_example.htm (http://en.wikipedia.org/wiki/Predictor-corrector_method#Euler_trapezoidal_example)

[fluidistic]

I'm not sure what the issue with the combined programs is, but if you want to try another method, the link below is to a description of an iterative (predictor + corrector) method. You could try doing 4 iterations (the initial Euler + 3 trapezoidals, (y[3])) to compare with Runge-Kutta with about the same overhead. With more iterations, this method will converge to specific values for each step, but since it's a linear approximation for each step (trapezoidal), you still need to keep h relatively small.

Euler_trapezoidal_example.htm (http://en.wikipedia.org/wiki/Predictor-corrector_method#Euler_trapezoidal_example)

Thank you for your help. I'm going to check this out tomorrow, I'm off to bed now (over 2 am).
My last thoughts for today is that not only the program that merges both codes is wrong for Runge Kutta but even my Runge Kutta (2nd code of 1st post) MUST be wrong.
I hate it when I have an error that is so small that I discover it only after looking for hours on what's wrong :/

[fluidistic]

Some news:
I've written a few other codes (also with another methods to solve ODE). I plotted the error for Euler, Heun and Runge-Kutta's methods.
One thing I know for sure now is that there's something wrong in all my codes when it comes to Runge Kutta.
For another ODE than the one here, I've tested Runge-Kutta's method along with Heun and Euler. The interval was [0,10] and the step size was 0.1. Euler was more efficient than Heun by a few. Runge-Kutta's error was almost exactly the same as Heun for all points. It means my Euler's method is more efficient than my Runge-Kutta's one.

If someone could spot my error(s) for Runge-Kutta's code, I'd be glad.
I'll be looking at it over and over again but as of now I'm totally stuck. All "seems" good.

Edit: I'm turning crazy! In fact Heun's method as well as Runge-Kutta's one are supposed to be better than Euler's method. So it means I have errors in both Runge-Kutta's and Heun codes! I've rechecked the algorithm of Runge-Kutta and couldn't spot a single mistake.

[rcgldr]

Did you try the predictor-corrector method?

Euler_trapezoidal_example.htm (http://en.wikipedia.org/wiki/Predictor-corrector_method#Euler_trapezoidal_example)

To convert the euler function to the predictor corrector method, you would add a variable "yg" as the intermediate "guess" value, and then replace the line with y=y+h*f(y,t) with the code fragment shown below starting with yg = ...

real(8):: yg
...
yg = y+h*f(y,t)
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
... repeat the above line more if wanted
y = yg

You could use a variable to hold the expresion t+h since it's used multiple times:

real(8):: yg, th
...
th = t+h
yg = y+h*f(y,t)
yg = y+(h/2d0)*(f(y,t)+f(yg,th))
yg = y+(h/2d0)*(f(y,t)+f(yg,th))
yg = y+(h/2d0)*(f(y,t)+f(yg,th))
... repeat the above line more if wanted
y = yg
...
t = th

[fluidistic]

Have you tried the predictor-corrector algorithm yet to see if that helps?

Euler_trapezoidal_example.htm (http://en.wikipedia.org/wiki/Predictor-corrector_method#Euler_trapezoidal_example)

To conver the euler function to the predictor corrector method, you would add the variable yg as the intermediate "guess" value, and then
replace the line with y=y+h*f(y,t) with the code fragment shown below starting with yg = ...

real(8):: yg
...
yg = y+h*f(y,t)
yg = y+(h/2d0)*(f(y,t)+f(yg,t))
yg = y+(h/2d0)*(f(y,t)+f(yg,t))
yg = y+(h/2d0)*(f(y,t)+f(yg,t))
... repeat the above line more if wanted
y = yg
Just did with your 4 lines for yg.
I took the ODE of the first post.
For h=0.1. Interval [a,b]=[0,1].
See picture attachment. Seems like your algorithm is indeed better than Euler's as it should. The picture represents the error for the 3 methods. So I think Euler and your Euler-trapezium work fine. My Runge-Kutta and my Heun seem not to work well.
I still haven't figured out what's wrong. I've consulted 2 books+several PDF as well as constantly reading my code... I don't see any error. :/

[fluidistic]

I've added my Heun algorithm for the same ODE. Now that's very weird!
The only difference with the previous post is that now I chose h=0.008 which is slightly smaller than 0.01.
My Heun's method is almost like my Runge-Kutta (but give slightly different values), this is so weird...
Have a look at the attached picture. My code is quite long now, I can post it if you like. But there's definitely something wrong with my Runge-Kutta and Heun's codes.

[rcgldr]

Just did with your 4 lines for yg.

yg = y+(h/2d0)*(f(y,t)+f(yg,t))



I corrected my previous post those lines should use f(yg, t+h), instead of f(yg, t):

yg = y+h*f(y,t)
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
... repeat the above line more if wanted
y = yg

Try this change and see if it's any better.

[fluidistic]

I corrected my previous post those lines should use f(yg, t+h), instead of f(yg, t):

yg = y+h*f(y,t)
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
yg = y+(h/2d0)*(f(y,t)+f(yg,t+h))
... repeat the above line more if wanted
y = yg

Try this change and see if it's any better.

Omg... now I'm clueless. Look at what I get: Runge-Kutta, Heun and your Euler/trapezium seem to give almost the same error. While Euler stands alone.
I really have no clue about what's going on.

[rcgldr]

(Note - I updated my response now that I figured out what the problem is)

I see the problem now. You're calculating the exact solution after you've done an integration step, which corresponds to (t+h). In this case, the exact solution needs to also be based on (t+h) and not (t). Move the t=t+h so it's done before you do the exact solution, the code should look like this:

t=t+h
valor_exacto= ...
error=abs(y-valor_exacto)

I used a spreadsheet to test this using the euler + trapezoidal implicit method and these error graphs are the results with step size of .01 and .001. There's still a pattern to the error values, but the maximum error values are relatively small, 0.000085 = 8.5 x 10-5 for step size .01, 0.00000085 = 8.5 x 10-7 for step size .001. The graphs show absolute values, but the actual error values are negative at around .78 to 1.0.

[fluidistic]

(Note - I updated my response now that I figured out what the problem is)

I see the problem now. You're calculating the exact solution after you've done an integration step, which corresponds to (t+h). In this case, the exact solution needs to also be based on (t+h) and not (t). Move the t=t+h so it's done before you do the exact solution, the code should look like this:

t=t+h
valor_exacto= ...
error=abs(y-valor_exacto)

I used a spreadsheet to test this using the euler + trapezoidal implicit method and these error graphs are the results with step size of .01 and .001. There's still a pattern to the error values, but the maximum error values are relatively small, 0.000085 = 8.5 x 10-5 for step size .01, 0.00000085 = 8.5 x 10-7 for step size .001. The graphs show absolute values, but the actual error values are negative at around .78 to 1.0.
Thank you so much. A huge thanks to you!
I modified all codes. Now I see that Euler's method is much less accurate than all the others. Runge-Kutta's method is by far the most accurate. And your Euler/trapezium and my Heun's method are almost equivalent (same order of magnitude error).
I'll post several pictures.
By the way, I still don't understand why what I did was wrong. You said You're calculating the exact solution after you've done an integration step, which corresponds to (t+h). In this case, the exact solution needs to also be based on (t+h) and not (t).
But to me I was calculating the exact solution with exactly the same value for t I used to calculate y. I don't understand why moving the line t=t+h before "valor_exacto" works. :confused:
I'd never have guessed what the error was... thanks once again.

[rcgldr]

You're calculating the exact solution after you've done an integration step, which corresponds to (t+h) ... not (t)...

Runge-Kutta's method is by far the most accurate.I got the same results using a spreadsheet.

I still don't understand why what I did was wrong. You said But to me I was calculating the exact solution with exactly the same value for t I used to calculate y. I don't understand why moving the line t=t+h before "valor_exacto" works.What each integration step calculates is yn+1, using the current yn, tn, and h as input. The initial value y0 corresponds to valor_exacto(t0). Your first integration step calculates y1, which corresponds to valor_exacto(t1), and t1 = t0+h. Note that the Euler-Trapezoidal method and the last term in the Runge-Kutta method also use t+h as input (f(..., t+h)) to calculate those terms.

[fluidistic]

Ah thanks rcgldr! This makes sense.