- #1

fluidistic

Gold Member

- 3,934

- 263

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:

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

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.

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

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