文档章节

RK4

ludlows
 ludlows
发布于 2014/10/07 22:55
字数 475
阅读 14
收藏 0
点赞 0
评论 2
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


© 著作权归作者所有

共有 人打赏支持
ludlows
粉丝 0
博文 15
码字总数 4195
作品 0
海淀
程序员
加载中

评论(2)

ludlows
ludlows

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

. 什么语言 没看懂...
fortran 这语言太老了啊
明天以后
明天以后
. 什么语言 没看懂...
利用RK4更新四元數並轉換成歐拉角及 OpenGL 的 Viewer

開發平台 : 硬體規格 Odroid-U2 作業系統 Ubuntu 13.10 Codename: saucy 下面都是一些基本的程式可示需求修改, 之後我會把算法跟GUI整再一起開發. 四元數更新 : Windows : dev-c++ 專案 : qu...

cctsao1008
2014/11/28
0
2

没有更多内容

加载失败,请刷新页面

加载更多
CORS 跨域实践

本文首发于个人微信公众号《andyqian》,期待你的关注~ 前言 系统通常都是由单体应用逐渐演化而来,演化成为前后端分离的分布式应用。在享受分布式系统带来的诸多好处之时,随之而来的也有不...

andyqian
8分钟前
6
0
开源 java CMS - FreeCMS2.8 会员管理

项目地址:http://www.freeteam.cn/ 会员组管理 会员管理 会员管理 从左侧管理菜单点击会员管理进入。 添加会员 在会员列表下方点击“添加”按钮。 填写相关属性后点击“保存”按钮即可。 编...

freeteam
9分钟前
0
0
bboss升级至 v5.0.6.8版本,改善对Elasticsearch SQL 的支持

v5.0.6.8功能改进如下: (1)持久层支持支持Elasticsearch SQL,使用参考文档:玩转Elasticsearch SQL功能 (2)解决持久层/elasticsearch模板变量解析多层级不起作用问题 (3)完善国际化功能 (4...

linux-tao
11分钟前
0
0
扫码二维码跳转到某个网站

添加maven依赖 <dependency><groupId>com.google.zxing</groupId><artifactId>core</artifactId><version>3.0.0</version></dependency><dependency><groupId>com.goog......

gaomq
17分钟前
0
0
Windows平台下搭建Git服务器的图文教程

Git没有客户端服务器端的概念,但是要共享Git仓库,就需要用到SSH协议(FTP , HTTPS , SFTP等协议也能实现Git共享,此文档不讨论),但是SSH有客户端服务器端,所以在windows下的开发要把自己...

MKChan
23分钟前
0
0
告警系统主脚本&告警系统配置文件&告警系统监控项目

20.20 告警系统主脚本 准备工作 定义监控系统的各个目录,然后再去定义主脚本,因为是分布式的,所以需要每一台机器都需要定义,事先创建好各个脚本和各个目录,随后脚本直接拷贝过去即可,然...

影夜Linux
24分钟前
0
0
谈谈神秘的ES6——(一)初识ECMAScript

谈谈神秘的ES6——(一)初识ECMAScript 在《零基础入门JavaScript》我们就说过,ECMAScript是JavaScript的核心,是JavaScript语法和语义的解释器,同时也是一个标准。而ECMAScript标准其实也...

JandenMa
今天
1
0
第16章 Tomcat配置

16.1 Tomcat介绍 ####Tomcat介绍 LNMP架构针对的开发语言是PHP语言,php 是一门开发web程序非常流行的语言,早些年流行的是asp,在Windows平台上运行的一种编程语言,但安全性差,就网站开发...

Linux学习笔记
今天
1
0
流利阅读笔记29-20180718待学习

高等教育未来成谜,前景到底在哪里? Ray 2018-07-18 1.今日导读 在这个信息爆炸的年代,获取知识是一件越来越容易的事情。人们曾经认为,如此的时代进步会给高等教育带来众多便利。但事实的...

aibinxiao
今天
12
0
OSChina 周三乱弹 —— 你被我从 osc 老婆们名单中踢出了

Osc乱弹歌单(2018)请戳(这里) 【今日歌曲】 @小鱼丁:分享五月天的单曲《后来的我们 (电影《后来的我们》片名曲)》: 《后来的我们 (电影《后来的我们》片名曲)》- 五月天 手机党少年们想...

小小编辑
今天
810
20

没有更多内容

加载失败,请刷新页面

加载更多

下一页

返回顶部
顶部