Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran 90. Euler and Runge-Kutta's method

  1. Jun 22, 2011 #1

    fluidistic

    User Avatar
    Gold Member

    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"
    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 (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"
    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 (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"
    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.
     

    Attached Files:

  2. jcsd
  3. Jun 22, 2011 #2

    rcgldr

    User Avatar
    Homework Helper

    What does the graph of the non-combined routines look like?
     
  4. Jun 22, 2011 #3

    fluidistic

    User Avatar
    Gold Member

    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:

  5. Jun 22, 2011 #4

    fluidistic

    User Avatar
    Gold Member

    Here with h=0.2.
     

    Attached Files:

  6. Jun 22, 2011 #5

    rcgldr

    User Avatar
    Homework Helper

    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
  7. Jun 23, 2011 #6

    fluidistic

    User Avatar
    Gold Member

    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 :/
     
  8. Jun 23, 2011 #7

    fluidistic

    User Avatar
    Gold Member

    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
      hmm.jpg
      File size:
      24.1 KB
      Views:
      278
    Last edited: Jun 23, 2011
  9. Jun 23, 2011 #8

    rcgldr

    User Avatar
    Homework Helper

    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
  10. Jun 23, 2011 #9

    fluidistic

    User Avatar
    Gold Member

    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
      okay.jpg
      File size:
      24.9 KB
      Views:
      240
  11. Jun 23, 2011 #10

    fluidistic

    User Avatar
    Gold Member

    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:

  12. Jun 23, 2011 #11

    rcgldr

    User Avatar
    Homework Helper

    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.
     
  13. Jun 23, 2011 #12

    fluidistic

    User Avatar
    Gold Member

    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:

    Last edited: Jun 23, 2011
  14. Jun 24, 2011 #13

    rcgldr

    User Avatar
    Homework Helper

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

    Last edited: Jun 24, 2011
  15. Jun 24, 2011 #14

    fluidistic

    User Avatar
    Gold Member

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

    Attached Files:

  16. Jun 24, 2011 #15

    rcgldr

    User Avatar
    Homework Helper

    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
  17. Jun 24, 2011 #16

    fluidistic

    User Avatar
    Gold Member

    Ah thanks rcgldr! This makes sense.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Fortran 90. Euler and Runge-Kutta's method
  1. Euler to Runge-Kutta 4 (Replies: 19)

Loading...