## RK4 原

ludlows

``````MODULE Global_data
!Symbolic names for kind types of single- and double-precision reals:
INTEGER, PARAMETER :: SP = KIND(1.0)
INTEGER, PARAMETER :: DP = KIND(1.0D0)
!Frequently used mathematical constants (with precision to spare):
REAL(DP), PARAMETER :: Pi=3.141592653589793238462643383279502884197_dp
! order of the problem
INTEGER, PARAMETER :: no_of_equations=2
! parameters for the equation
REAL(DP) :: y_initial(no_of_equations)
REAL(DP) :: x_start,x_finish,h,kappa
END MODULE Global_data

MODULE RungeKutta_method
! here we just need dp, kappa, h and no_of_equations
! if we set wp=>sp we can change to single precision numbers
USE global_data, ONLY : wp=>dp,h,kappa,no_of_equations
IMPLICIT NONE

CONTAINS

FUNCTION func(x,y)
!    we assign y to have 'no_of_equations's elements throughout to keep the program consistent
REAL(wp) , INTENT(IN) :: x,y(no_of_equations)
REAL(wp) :: func(no_of_equations)
!   The ODE here is y'' + kappa y' + x y = 0
!   so that the set of first ODEs is
!    y_1' = y_2
!    y_2' = - kappa y_2 - x y_1
! where we write y_1 = y and y_2 = y'
func(1) = y(2)
func(2) = -kappa*y(2) - x*y(1)

END FUNCTION func

FUNCTION rungekutta4(x,y)

REAL(wp) , INTENT(IN) :: x,y(no_of_equations)
REAL(wp),DIMENSION(no_of_equations) :: rungekutta4,k1,k2,k3,k4
! Note here that k1,k2,... etc are two element arrays.
k1 = h*func(x,y)
k2 = h*func(x+0.5_wp*h,y+0.5_wp*k1)
k3 = h*func(x+0.5_wp*h,y+0.5_wp*k2)
k4 = h*func(x+h,y+k3)
rungekutta4 = y + (k1+2._wp*(k2+k3) + k4)/6._wp

END FUNCTION rungekutta4

END MODULE RungeKutta_method

PROGRAM euler
! we need all the variables from global data
USE global_data, ONLY : wp=>dp,kappa,h,x_start,x_finish,y_initial,no_of_equations
! RungeKutta_method gives us our rungekutta method function
USE RungeKutta_method
IMPLICIT NONE

INTEGER :: i,no_of_steps
REAL(wp) :: x,y(no_of_equations)

! Setup initial conditions
! set x_0 = 0
x_start = 0._wp
! set x_n = 1
x_finish = 1._wp
! Set y(x=0) = 0
y_initial(1)=0._wp
! Set y'(x=0) = 1
y_initial(2)=1._wp
! Set the number of steps
no_of_steps = 10
! Set the parameter kappa from the ODE
kappa = 5._wp
! Set the step size
h = (x_finish-x_start)/dble(no_of_steps)

! Open a file to write to and output some info to screen
WRITE(6,*)' Run with ::',no_of_steps,' steps.'
OPEN(unit=10,file="runge_kutta.dat")

! Set the start values of x and y
x = x_start
y = y_initial
WRITE(6,'(3(E15.8,1X))')x,y
WRITE(10,'(3(E15.8,1X))')x,y

! Loop over the required no. of steps to reach x_finish
DO i=1,no_of_steps
! update the values of x and y
y = rungekutta4(x,y)
x = x + h

! pick out only 10 value to print at.
if(mod(i,no_of_steps/10)==0)THEN
! print to screen and file
WRITE(6,'(3(E15.8,1X))')x,y
WRITE(10,'(3(E15.8,1X))')x,y
END if

END DO

CLOSE(unit=10)

END PROGRAM euler``````

### 评论(2)

#### 引用来自“明天以后”的评论

. 什么语言 没看懂...
fortran 这语言太老了啊
. 什么语言 没看懂...

cctsao1008
2014/11/28
0
2

IDE 插件新版本发布，开发效率 “biu” 起来了

28分钟前
6
0

56分钟前
3
0
Kubernetes Client-go Informer 源码分析

4
0

11. 重置密码 密钥和密码都支持远程登陆， 二选一 两个都可以登陆， 密钥相对于密码来说，相对安全一点 本地登陆无法是用密钥 修改密码 root 用户 passwd root 修改普通用户 passwd usernam...

miko0089

6
0
bash特性

1.支持别名 alias 2.命令替换 \$(COMMANS) 或者 `COMMAND` 3. bash支持的引号： `` :命令替换 ""：弱引用，可以实现变量替换 ''：强引用，不实现变量替换 4.文件名通配 globbing：（man 7 glo...

3
0