 #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 RungeKutta'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 stepsizes 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 RungeKutta'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 RungeKutta 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 stepsize 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 stepsize 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 RungeKutta'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 stepsize 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(yvalor_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(yvalor_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.
Attachments

26.1 KB Views: 1,225