Help with Runge-Kutta 4

Help Homework
1

Hello, I’m trying to program a subroutine that solves differential equations using the Runge-Kutta 4 method. To try if it’s correct I’m implementing it on a pendulum equation so I can plot it and check if the result is correct. Although i though it was a very simple formula, I am having some struggles because when excecuting I have the following message : Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

This messages somewtimes occur and I tried to search for the memory error but can seem to find it on the code. Can someone give it a go?

Thanks!

PROGRAM PREPA8
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud
	INTEGER:: N, i, iter, nequ
	DOUBLE PRECISION:: yin_0(2), yy1(2)
	EXTERNAL:: EDO, pendolSimple
	COMMON/DADES/V

	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	grav = 3.71d0 !m/s^2
	longitud = 0.45d0 !m

	n=1300 !número de passos
	!Inicialitzem els valors amb les condicions inicials
	yin_0(1) = 0.02d0
	yin_0(2) = 0d0
	dx=0.0003d0

	do i = 1, N
		call RungeKutta4order(nequ, dx, yin_0,yy1, pendolSimple)
		write(*,*) dx*i, yy1
		!canviem les variables pel bucle
		yin_0=yy1
	enddo
	
END PROGRAM PREPA8

	IMPLICIT NONE
	INTEGER:: nequ, i
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), k1(nequ), k2(nequ), k3(nequ), k4(nequ), yaux(nequ)
	DOUBLE PRECISION:: x, dx

	x=0d0 !realment no necessitem x

	call func(nequ, x, yyin, k1) !k1 = f(x0, y0)
	yaux=yyin+(dx+k1/2d0)
	call func(nequ, x+(dx/2d0), yaux, k2)
	yaux=yyin +(dx*k2/2d0)
	call func(nequ, x+(dx/2d0), yaux, k3)
	yaux=yyin+(dx*k3)
	call func(nequ, x+dx, yaux, k4)

	yyout=yyin+(dx/6d0)*(k1+k2*2d0+2d0*k3+k4)

END SUBROUTINE RungeKutta4order

SUBROUTINE pendolSimple(yin, f)
	IMPLICIT NONE
	DOUBLE PRECISION:: yin(2), f(2), grav, longitud, massa
	grav = 3.71d0 !m/s^2
	longitud = 0.45d0 !m
	f(1) = yin(2)
	f(2) = -(grav/longitud) * sin(yin(1))
	return

END SUBROUTINE pendolSimple ````
2

Hard to say since this example isn’t complete.

  • you aren’t defining nequ anywhere
  • set some compiler flags to warn you of uninitialized variables.
  • put all your routines in a module
  • use a runge-kutta library instead of reinventing the wheel
1 Like
3

What do you mean it isn’t complete?
I never use modules and they specifically asked me to do a Runge Kutta subroutine. I inicialized nequ but still getting the error

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0 0x101e3aafe
#1 0x101e39cdd
#2 0x7fff2041dd7c
#3 0x101e28659
#4 0x101e28978
#5 0x101e28adf
zsh: segmentation fault ./a.out

	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud
	INTEGER:: N, i, iter
	INTEGER, PARAMETER:: nequ=2
	DOUBLE PRECISION:: yin_0(nequ), yy1(nequ),
	EXTERNAL:: pendolSimple
	COMMON/DADES/V

	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	grav = 3.71d0 !m/s^2
	longitud = 0.45d0 !m

	n=1300 !número de passos
	!Inicialitzem els valors amb les condicions inicials
	yin_0(1) = 0.02d0
	yin_0(2) = 0d0
	dx=0.0003d0

	do i = 1, N
		call RungeKutta4order(2, dx, yin_0,yy1, pendolSimple)
		write(*,*) dx*i, yy1
		!canviem les variables pel bucle
		yin_0 = yy1
	enddo
	
END PROGRAM PREPA8

SUBROUTINE RungeKutta4order(nequ, dx,yyin,yyout, func)
	IMPLICIT NONE
	INTEGER:: nequ, i
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), k1(nequ), k2(nequ), k3(nequ), k4(nequ), yaux(nequ)
	DOUBLE PRECISION:: x, dx

	x=0d0 !realment no necessitem x

	call func(nequ, x, yyin, k1) !k1 = f(x0, y0)
	yaux=yyin+(dx+k1/2d0)
	call func(nequ, x+(dx/2d0), yaux, k2)
	yaux=yyin +(dx*k2/2d0)
	call func(nequ, x+(dx/2d0), yaux, k3)
	yaux=yyin+(dx*k3)
	call func(nequ, x+dx, yaux, k4)

	yyout=yyin+(dx/6d0)*(k1+k2*2d0+2d0*k3+k4)

END SUBROUTINE RungeKutta4order

SUBROUTINE pendolSimple(yin, f)
	IMPLICIT NONE
	DOUBLE PRECISION:: yin(2), f(2), grav, longitud, massa
	grav = 3.71d0 !m/s^2
	longitud = 0.45d0 !m
	f(1) = yin(2)
	f(2) = -(grav/longitud) * sin(yin(1))
	return

END SUBROUTINE pendolSimple
4

What strikes me with your code:

Your Subroiutine pendolSimple is defined with two paramters but you call it from your RungeKutta subroutine with four.

Ah, and by the way, if grav is the Earth’s gravity you should check this value. And I doubt you would need double precision for this task, especially when you run it on a 64 bit machine.
Cheers
Norbert

5

Thank you, it was indeed the fact the number of parameters.

1 Like
6

I think an advantage of putting all functions and subroutines into a module is that they will have an explicit interface and the compiler will check that you are calling the function with the correct arguments. What I found weird of your example is that you got this mistake only at run time and the compiler was not able to flag this during compilation.
Using modules makes your code easier to debug and to find mistakes

1 Like
7

I think the code you are pasting here is not the code you are compiling. A compiler would have caught some of these syntax errors (such as the missing subroutine RungeKutta4order declaration in the first one, or the extra comma in the second one).

1 Like
8

Obviously, this is a student’s homework. Which is ok, but apparently you missed at least one line of code while copy/pasting.

The fact they asked you to do a subroutine doesn’t mean said subroutine can’t be in a module.
Also, if you want to learn Fortran, you really need to learn how to use modules, so that old fashion statements (which you use in the code, like common and external) won’t be needed. The code will be much easier to understand (and expand later on.)

Fortran is much easier to learn compared to other languages like C or C++, but you still need to put some effort into learning it…

2 Likes
9

Thank you for the advise, will start to use modules from now on.
I have encountered another problem while creating a subroutine that is supposed to return the excited level of a particle. To do so I created a subroutine integrate that integrates the Schrodinger equation using Runge-Kutta4 but it seems that it’s not updating how it’s suppose to be yyin and yyout at each iteration. Does anyone know what the problem might be?

SUBROUTINE integrate(nequ, E0, a, b, N, yyin, phi)
	IMPLICIT NONE
	DOUBLE PRECISION:: E0, phi, a, b
	DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ), local_yyin(nequ)
	INTEGER:: nequ, i, N
	EXTERNAL:: EDO

	dx=(b-a)/N
	local_yyin=yyin
	DO i = 1, N
		call RungeKutta4order(nequ, dx, local_yyin, yyout, edo)
		print*, dx*i, yyin, yyout
		local_yyin=yyout
	ENDDO

	phi=local_yyin(1)


END SUBROUTINE integrate

SUBROUTINE secant(E1, E2, phi1, phi2, E3)
	IMPLICIT NONE
	DOUBLE PRECISION:: E1, E2, phi1, phi2, E3

	E3=(E1 * phi2 - E2 * phi1) / (phi2 - phi1)
END SUBROUTINE secant

SUBROUTINE tir(nequ, E1, E2, error, a,b,N, yyin)
	IMPLICIT NONE
	DOUBLE PRECISION:: E1, E2, E3, a, b
	DOUBLE PRECISION:: phi1, phi2, phi3, dx, x, error, yyin(nequ), yyout(nequ)
	INTEGER:: i, nequ, N, iter

	i=0 !incialitzo
	iter=500
	!Pas 1 i 2

	do i = 1, iter
		call integrate(nequ, E1, a, b, N,yyin, phi1)
		call integrate(nequ, E2, a, b, N, yyin, phi2)
		print*, phi1, phi2

		call secant(E1, E2, phi1, phi2, E3)
		call integrate(nequ, E3, a, b, N,yyin, phi3)

		if (abs(phi3).LT.error) then
			write(*,*) 'Ha convergit', i, E3
			exit

	    else
	    	E1 = E2
		    E2 = E3
		end if
	end do	


END SUBROUTINE tir
SUBROUTINE RungeKutta4order(nequ, dx, yyin, yyout, func)
	IMPLICIT NONE
	INTEGER:: nequ, i
	DOUBLE PRECISION, DIMENSION(nequ):: yyin, yyout, k1, k2, k3, k4, yaux
	DOUBLE PRECISION:: x, dx

	call func(nequ, x, yyin, k1) !k1 = f(x0, y0)

	yaux=yyin+(dx*k1/2d0)
	x=x+0.5d0*dx
	call func(nequ, x, yaux, k2) !k2
	
	yaux=yyin +(dx*k2/2d0)
	x=x+0.5d0*dx
	call func(nequ, x, yaux, k3) !k3
	
	yaux=yyin+(dx*k3)
	x=x+0.5d0*dx
	call func(nequ, x+dx, yaux, k4) !k4

	yyout=yyin+(dx/6d0)*(k1+k2*2d0+2d0*k3+k4)
	
END SUBROUTINE RungeKutta4order

!L'edos en aquest cas és equivalent a la del pèndol simple en aproximacions petites.
!yin=(phi, dphi/dx)
!dyout=(dphi/dx, d^2phi/dx^2)
SUBROUTINE EDO(nequ, x, yin, dyout)
	IMPLICIT NONE
	DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V
	INTEGER:: nequ
	COMMON/ENERGIES/E, V

	dyout(1)= yin(2)
	dyout(2) = 2D0*(V-E)* yin(1) !Correspon a la segona derivada que ens dona l'equació de Schrodinger

END SUBROUTINE EDO

PROGRAM PREPA8
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b
	INTEGER:: N, i, iter
	INTEGER, PARAMETER:: nequ=2
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), phi(200), phi_f(2)
	EXTERNAL:: EDO, pendolSimple
	COMMON/DADES/V, pi

	OPEN(10, file='prepa8.dat')
	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	
	yyin(1)=0.02d0
	yyin(2)=0.0d0
	E1 = 1d0
	E2= 3d0
	error=10-5d0
	call tir(nequ, E1, E2, error, 0d0,1d0,500, yyin)
END PROGRAM PREPA8
10

The use of double precision or not does not depend on the architecture of the machine.

1 Like
11

At least, the variables in the ENERGY common block are not assigned any value.

Please, do yourself a favor: stop using common blocks and put everything in a module.

4 Likes
12

Anyway, even 32 bit precision should be sufficient to check the outcome. But of course ‘DOUBLE PRECISION’ does not do any harm.

13

Normally double precision isn’t needed even for a double pendulum with a friction factor added in the equations - unless someone needs insane precision, I guess.

Note that in this case we have a fixed step Runge-Kutta method - which is fine for educational purposes, but should never be used for real-world problems. Obviously this is student’s homework for some class though, and a fixed step Runge-Kutta is used to introduce the general idea of such methods.

14

I do define the variables in the COMMON at my main program (at the end of the code). I was told to use COMMONS but it’s definitely not necessary in this case so I will delete it. I still have the problem at the subroutine integrate where the values I print during the loop are always the initial values, they don’t get updated even though I wrote yyin=yyout at the end of the loop. I believe this is why the subroutine tir doesn’t work.
I am learning to program and debugging my own codes but it can get very frustrating so help and advice is always great.
Can someone see the problem with the code?? I really would appreciate some help

15

The first thing you should do, and do it religiously in any program, is to always add the intent attribute in all subroutine/function arguments see, e.g., this for a quick description of the attribute.)

For example, in your integrate subroutine, I guess phi is intent(out) - or is it intent(inout)? Do yourself a favor and learn to use intent. Not only it makes your code more readable, but doing that alone may help you find where the problem is. Do it, and I’ll have a better look in the code.

Also, when you ask for help, please try to post a minimal program that reproduces the problem. Posting long code makes people way less willing to help.

2 Likes
16

Whoever told you that in 2023 has a few decade to catch up with Fortran…

You are printing yyin, but your update is yyin_local=yyout. I guess you want to print yyin_local

2 Likes
17 
PROGRAM PREPA8
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b
	INTEGER:: N, i, iter
	INTEGER, PARAMETER:: nequ=2
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), phi(200), phi_f(2)
	EXTERNAL:: EDO, pendolSimple
	COMMON/DADES/V, pi

	OPEN(10, file='prepa8.dat')
	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	
	yyin(1)=0.02d0
	yyin(2)=0.0d0
	E1 = 1d0
	E2= 1.1d0
	error=10-5d0
	call tir(nequ, E1, E2, error, 0d0,1d0,500, yyin)

	CLOSE(10)

END PROGRAM PREPA8

! Subroutine for numerical integration using the Runge-Kutta 4th order method
! Purpose: Numerically integrates a system of ordinary differential equations (ODEs)
! Inputs:
!   nequ: Number of equations in the system
!   E0: Initial value for a parameter (input)
!   a, b: Integration limits (input)
!   N: Number of integration steps (input)
!   yyin: Initial values of the dependent variables (input)
! Outputs:
!   phi: Final value of the first dependent variable after integration (output)
SUBROUTINE integrate(nequ, E0, a, b, N, yyin, phi)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: nequ, N
  DOUBLE PRECISION, INTENT(IN) :: E0, a, b
  DOUBLE PRECISION, INTENT(IN) :: yyin(nequ)
  DOUBLE PRECISION, INTENT(OUT) :: phi
  DOUBLE PRECISION :: dx, yyout(nequ)
  INTEGER :: i
  EXTERNAL :: EDO, RungeKutta4order  ! External function for ODE and Runge-Kutta method

  ! Calculate the step size for integration
  dx = (b - a) / N

  ! Perform numerical integration using Runge-Kutta 4th order method
  DO i = 1, N
    CALL RungeKutta4order(nequ, dx, yyin, yyout, edo)
    WRITE(*,*) dx * i, yyout
    yyin = yyout
  ENDDO

  ! Set the final value of the first dependent variable
  phi = yyin(1)
END SUBROUTINE integrate

! Subroutine for the secant method to find the root of a function
! Purpose: Implements the secant method
! Inputs:
!   E1, E2: Initial guesses for the root (input)
!   phi1, phi2: Function values corresponding to the initial guesses (input)
! Outputs:
!   E3: Improved estimate for the root (output)
SUBROUTINE secant(E1, E2, phi1, phi2, E3)
  IMPLICIT NONE
  DOUBLE PRECISION, INTENT(IN) :: E1, E2, phi1, phi2
  DOUBLE PRECISION, INTENT(OUT) :: E3

  ! Apply the secant method formula to improve the root estimate
  E3 = (E1 * phi2 - E2 * phi1) / (phi2 - phi1)
END SUBROUTINE secant

! Subroutine for solving a nonlinear equation using the secant method and numerical integration
! Purpose: Solves a nonlinear equation iteratively
! Inputs:
!   nequ: Number of equations in the system (input)
!   E1, E2: Initial guesses for the root (input)
!   error: Tolerance for convergence (input)
!   a, b: Integration limits (input)
!   N: Number of integration steps (input)
!   yyin: Initial values of the dependent variables (input)
SUBROUTINE tir(nequ, E1, E2, error, a, b, N, yyin)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: nequ, N
  DOUBLE PRECISION, INTENT(IN) :: E1, E2, error, a, b
  DOUBLE PRECISION, INTENT(INOUT) :: yyin(nequ)  ! Changed to INOUT for modification inside the subroutine
  INTEGER :: i, iter

  ! Initialize counters
  i = 0
  iter = 500

  ! Iteratively apply the secant method along with numerical integration
  DO i = 1, iter
    CALL integrate(nequ, E1, a, b, N, yyin, phi1)
    CALL integrate(nequ, E2, a, b, N, yyin, phi2)

    ! Display intermediate results
    PRINT*, 'Iteration:', i, 'E1:', E1, 'E2:', E2, 'phi1:', phi1, 'phi2:', phi2

    ! Apply the secant method to get an improved root estimate
    CALL secant(E1, E2, phi1, phi2, E3)
    CALL integrate(nequ, E3, a, b, N, yyin, phi3)

    ! Check for convergence
    IF (ABS(phi3) < error) THEN
      WRITE(*,*) 'Converged after', i, 'iterations. Root:', E3
      EXIT
    ELSE
      ! Update guesses for the root
      E1 = E2
      E2 = E3
    ENDIF
  ENDDO
END SUBROUTINE tir

These are the subroutines I’m using alongside the Runge-Kutta and Schrodinger diferential equation(EDO). I added a print at the integrate subroutine to check on the values. I’m getting all the time the same values (the intial ones) yyout even though I defined the update yyin=yyout. I think Runge-Kutta is okay because I tried it with a pendulum equation and it worked just fine. Does anyone see the problem

1 Like
18

Sorry, but the code you have posted in your last message won’t compile at all, so it can’t be the code you are running. If you want some help try posting the exact code you are running.

1 Like
19

@4ulalia I would suggest that you put the main program and all functions and subroutines that you are using for this project in a single post here, adding also the instructions for compilation. Make sure that the code that you post here compiles and gives the same results. Then we can help :slight_smile:

20 

PROGRAM PREPA8
	IMPLICIT NONE
	DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b
	INTEGER:: N, i, iter
	INTEGER, PARAMETER:: nequ=2
	DOUBLE PRECISION:: yyin(nequ), yyout(nequ), phi(200), phi_f(2)
	EXTERNAL:: EDO, pendolSimple
	COMMON/DADES/V, pi

	OPEN(10, file='prepa8.dat')
	!Dades ctts::
	V=-2.4d0
	pi=4d0*atan(1d0)
	
	yyin(1)=0.02d0
	yyin(2)=0.0d0
	E1 = 1d0
	E2= 1.1d0
	error=10-5d0
	call tir(nequ, E1, E2, error, 0d0,1d0,500, yyin)

CLOSE(10)

END PROGRAM PREPA8

!SUBRUTINE RUNGE-KUTTA D'ORDRE 4
SUBROUTINE RungeKutta4order(nequ, dx, yyin, yyout, func)
	IMPLICIT NONE
	INTEGER:: nequ, i
	DOUBLE PRECISION, DIMENSION(nequ):: yyin, yyout, k1, k2, k3, k4, yaux
	DOUBLE PRECISION:: x, dx

	call func(nequ, x, yyin, k1) !k1 = f(x0, y0)

	yaux=yyin+(dx*k1/2d0)
	x=x+0.5d0*dx
	call func(nequ, x, yaux, k2) !k2
	
	yaux=yyin +(dx*k2/2d0)
	x=x+0.5d0*dx
	call func(nequ, x, yaux, k3) !k3
	
	yaux=yyin+(dx*k3)
	x=x+0.5d0*dx
	call func(nequ, x+dx, yaux, k4) !k4

	yyout=yyin+(dx/6d0)*(k1+k2*2d0+2d0*k3+k4)
	
END SUBROUTINE RungeKutta4order

!L'edos en aquest cas és equivalent a la del pèndol simple en aproximacions petites.
!yin=(phi, dphi/dx)
!dyout=(dphi/dx, d^2phi/dx^2)
SUBROUTINE EDO(nequ, x, yin, dyout)
	IMPLICIT NONE
	DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V
	INTEGER:: nequ
	COMMON/ENERGIES/E, V

	dyout(1)= yin(2)
	dyout(2) = 2D0*(V-E)* yin(1) !Correspon a la segona derivada que ens dona l'equació de Schrodinger

END SUBROUTINE EDO

UBROUTINE integrate(nequ, E0, a, b, N, yyin, phi)
	IMPLICIT NONE
	DOUBLE PRECISION:: E0, phi, a, b
	DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ)
	INTEGER:: nequ, i, N
	EXTERNAL:: EDO

	dx=(b-a)/N
	DO i = 1, N
		call RungeKutta4order(nequ, dx, yyin, yyout, edo)
		write(*,*) dx*i, yyout
		yyin=yyout
	ENDDO

	phi=yyin(1)


END SUBROUTINE integrate

SUBROUTINE secant(E1, E2, phi1, phi2, E3)
	IMPLICIT NONE
	DOUBLE PRECISION:: E1, E2, phi1, phi2, E3

	E3=(E1 * phi2 - E2 * phi1) / (phi2 - phi1)
END SUBROUTINE secant

SUBROUTINE tir(nequ, E1, E2, error, a,b,N, yyin)
	IMPLICIT NONE
	DOUBLE PRECISION:: E1, E2, E3, a, b
	DOUBLE PRECISION:: phi1, phi2, phi3, dx, x, error, yyin(nequ), yyout(nequ)
	INTEGER:: i, nequ, N, iter

	i=0 !incialitzo
	iter=500
	!Pas 1 i 2

	do i = 1, iter
		call integrate(nequ, E1, a, b, N,yyin, phi1)
		call integrate(nequ, E2, a, b, N, yyin, phi2)
		print*, phi1, phi2

		call secant(E1, E2, phi1, phi2, E3)
		call integrate(nequ, E3, a, b, N,yyin, phi3)

		if (abs(phi3).LT.error) then
			write(*,*) 'Ha convergit', i, E3
			exit

	    else
	    	E1 = E2
		    E2 = E3
		end if
	end do	


END SUBROUTINE tir

Here is my original code, it does compile. Again I have the same error I mentioned.