Fortran 90. Euler and Runge-Kutta's method

  • Fortran
  • Thread starter fluidistic
  • Start date
  • #1
fluidistic
Gold Member
3,744
128
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:
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
Code:
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 [itex]y(t)[/itex] in [itex]y'(t)=-y(t)+\sin (2 \pi t)[/itex] with the initial condition [itex]y(0)=1[/itex] 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:
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.
 

Attachments

Answers and Replies

  • #2
rcgldr
Homework Helper
8,721
543
What does the graph of the non-combined routines look like?
 
  • #3
fluidistic
Gold Member
3,744
128
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.
 

Attachments

  • #4
fluidistic
Gold Member
3,744
128
Here with h=0.2.
 

Attachments

  • #5
rcgldr
Homework Helper
8,721
543
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:
  • #6
fluidistic
Gold Member
3,744
128
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
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
fluidistic
Gold Member
3,744
128
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.
 

Attachments

Last edited:
  • #8
rcgldr
Homework Helper
8,721
543
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:
  • #9
fluidistic
Gold Member
3,744
128
Have you tried the predictor-corrector algorithm yet to see if that helps?

Euler_trapezoidal_example.htm

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

Attachments

  • #10
fluidistic
Gold Member
3,744
128
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.
 

Attachments

  • #11
rcgldr
Homework Helper
8,721
543
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.
 
  • #12
fluidistic
Gold Member
3,744
128
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.
 

Attachments

Last edited:
  • #13
rcgldr
Homework Helper
8,721
543
(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.
 

Attachments

Last edited:
  • #14
fluidistic
Gold Member
3,744
128
(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.
 

Attachments

  • #15
rcgldr
Homework Helper
8,721
543
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.
 
Last edited:
  • #16
fluidistic
Gold Member
3,744
128
Ah thanks rcgldr! This makes sense.
 

Related Threads on Fortran 90. Euler and Runge-Kutta's method

Replies
2
Views
12K
  • Last Post
Replies
19
Views
5K
Replies
5
Views
960
  • Last Post
Replies
2
Views
7K
  • Last Post
Replies
10
Views
17K
Replies
19
Views
13K
Replies
7
Views
1K
Replies
17
Views
19K
  • Last Post
Replies
5
Views
7K
Replies
1
Views
6K
Top