# Fortran 90. Euler and Runge-Kutta's method

1. Jun 22, 2011

### 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:
Code (Text):
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"

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
Code (Text):
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"
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:
Code (Text):
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"

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)

call euler(y_0,a,b,h)
close(1)
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.

#### Attached Files:

• ###### comparison.jpg
File size:
26.1 KB
Views:
533
2. Jun 22, 2011

### rcgldr

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

3. Jun 22, 2011

### fluidistic

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.

#### Attached Files:

File size:
19.2 KB
Views:
341
• ###### weird.jpg
File size:
18.1 KB
Views:
322
4. Jun 22, 2011

### fluidistic

Here with h=0.2.

#### Attached Files:

• ###### smallh.jpg
File size:
20.3 KB
Views:
325
5. Jun 22, 2011

### 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

Last edited: Jun 23, 2011
6. Jun 23, 2011

### fluidistic

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

7. Jun 23, 2011

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

#### Attached Files:

• ###### hmm.jpg
File size:
24.1 KB
Views:
302
Last edited: Jun 23, 2011
8. Jun 23, 2011

### rcgldr

Did you try the predictor-corrector method?

Euler_trapezoidal_example.htm

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

Last edited: Jun 23, 2011
9. Jun 23, 2011

### fluidistic

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. :/

#### Attached Files:

• ###### okay.jpg
File size:
24.9 KB
Views:
267
10. Jun 23, 2011

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

#### Attached Files:

• ###### strange.jpg
File size:
25.3 KB
Views:
238
11. Jun 23, 2011

### rcgldr

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.

12. Jun 23, 2011

### fluidistic

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.

#### Attached Files:

File size:
23.9 KB
Views:
311
• ###### precision.jpg
File size:
13.2 KB
Views:
243
Last edited: Jun 23, 2011
13. Jun 24, 2011

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

#### Attached Files:

File size:
25.3 KB
Views:
254
• ###### num_int22.jpg
File size:
21.5 KB
Views:
261
Last edited: Jun 24, 2011
14. Jun 24, 2011

### fluidistic

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
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.
I'd never have guessed what the error was... thanks once again.

#### Attached Files:

File size:
18.7 KB
Views:
297
File size:
14.5 KB
Views:
299
• ###### bigEuler.jpg
File size:
19.2 KB
Views:
287
15. Jun 24, 2011

### rcgldr

I got the same results using a spreadsheet.

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.

Last edited: Jun 24, 2011
16. Jun 24, 2011

### fluidistic

Ah thanks rcgldr! This makes sense.