# Fortran 90. Euler and Runge-Kutta's method

• Fortran
• fluidistic
In summary: That's the whole point of Runge-Kutta, so I'm not sure I'm helping much. : )In summary, the conversation discusses the creation and testing of two programs that numerically solve an ODE using Euler and Runge-Kutta's methods. The third program combines these two methods but produces errors. The exact solution and error calculations are also mentioned. The individual programs show relatively accurate results, but the combined program does not work as expected. Another iterative method is suggested for further testing.
fluidistic
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:
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"

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"
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 $y(t)$ in $y'(t)=-y(t)+\sin (2 \pi t)$ with the initial condition $y(0)=1$ 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"
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)

call euler(y_0,a,b,h)
close(1)
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

• comparison.jpg
26.1 KB · Views: 1,369
What does the graph of the non-combined routines look like?

rcgldr said:
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

• view.jpg
19.2 KB · Views: 1,056
• weird.jpg
18.1 KB · Views: 1,013
Here with h=0.2.

#### Attachments

• smallh.jpg
20.3 KB · Views: 1,074
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:
rcgldr said:
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 :/

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

• hmm.jpg
24.1 KB · Views: 1,023
Last edited:
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:
rcgldr said:
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

• okay.jpg
24.9 KB · Views: 995
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

• strange.jpg
25.3 KB · Views: 955
fluidistic said:
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.

rcgldr said:
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

• clueless.jpg
23.9 KB · Views: 1,285
• precision.jpg
13.2 KB · Views: 898
Last edited:
(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

• num_int21.jpg
24.8 KB · Views: 933
• num_int22.jpg
20.4 KB · Views: 950
Last edited:
rcgldr said:
(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.
I'd never have guessed what the error was... thanks once again.

#### Attachments

• trapeziumandHeun.jpg
18.7 KB · Views: 707
• rungekutta.jpg
14.5 KB · Views: 1,032
• bigEuler.jpg
19.2 KB · Views: 977
rcgldr said:
You're calculating the exact solution after you've done an integration step, which corresponds to (t+h) ... not (t)...

fluidistic said:
Runge-Kutta's method is by far the most accurate.
I got the same results using a spreadsheet.

fluidistic said:
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:
Ah thanks rcgldr! This makes sense.

## 1. What is Fortran 90 and how is it different from other versions of Fortran?

Fortran 90 is a programming language commonly used in scientific and engineering applications. It is an updated version of Fortran 77 and includes new features such as modules, pointers, and dynamic memory allocation. These additions make Fortran 90 more powerful and easier to use compared to its predecessors.

## 2. What is Euler's method and how is it used in Fortran 90?

Euler's method is a numerical method used for solving ordinary differential equations. It involves approximating the solution by taking small steps based on the derivative of the function. In Fortran 90, Euler's method can be implemented using loops and arrays to perform the necessary calculations.

## 3. How does Runge-Kutta's method compare to Euler's method in Fortran 90?

Runge-Kutta's method is a higher-order numerical method for solving ordinary differential equations. It is more accurate than Euler's method as it takes into account multiple derivatives at each step. In Fortran 90, Runge-Kutta's method can also be implemented using loops and arrays, but it requires more computations compared to Euler's method.

## 4. Can Fortran 90 be used for other scientific applications besides solving differential equations?

Yes, Fortran 90 can be used for various scientific applications such as data analysis, simulations, and statistical analysis. Its efficient and fast processing makes it a popular choice for scientific programming.

## 5. Are there any limitations or challenges when using Fortran 90 for scientific calculations?

One limitation of Fortran 90 is that it is not as user-friendly as other programming languages, and it may require a steep learning curve for beginners. Additionally, it may not be suitable for applications that require complex data structures or graphics. However, these challenges can be overcome with practice and familiarity with the language.

Replies
26
Views
5K
Replies
5
Views
1K
Replies
8
Views
1K
Replies
4
Views
1K
Replies
4
Views
2K
Replies
8
Views
1K
Replies
4
Views
3K
Replies
2
Views
5K
Replies
4
Views
2K
Replies
7
Views
2K