I've successfully (well I think so) made 2 different programs that can numerically solve an ODE using Euler and Runge-Kutta's methods.(adsbygoogle = window.adsbygoogle || []).push({});

Here they are:

andCode (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 programI 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.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 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:

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

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.

**Physics Forums - The Fusion of Science and Community**

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

Loading...

Similar Threads - Fortran Euler Runge | Date |
---|---|

How to buy Fortran Programming Software Commersial License | Jan 30, 2018 |

Euler method in Fortran | Nov 12, 2013 |

Fortran program for oscillator using Euler method | Oct 22, 2013 |

Euler Method in Fortran - HELP | Oct 30, 2008 |

Fortran 90, Euler's method help | May 4, 2008 |

**Physics Forums - The Fusion of Science and Community**