滑らかなものばかりだったので滑らかさのないものを扱ってみる。 両側を壁で囲まれた一次元の質点の運動について、 壁で瞬間的に衝突、反射するとすれば運動方程式からの数値解析が難しくなるので、 壁がばねで支えられているとして、その復元力をもって有限…

CやJavaに慣れるとFORTRANでの接尾辞なしの浮動小数点型リテラルが単精度であることを忘れてしまう。 前項でのLorenz方程式の初期値問題で与えた初期値 double precision, parameter :: x0(3) = (/1d-1, 1d-1, 1d-1/) ! initial coordinate at t0 が double …

有名どころでLorenz方程式を解いてみる。 方程式のパラメータはLorenzが最初の論文で示したものを使用する。 program ode use runge_kutta_4 implicit none double precision, parameter :: a = 10, b = 28, c = 8d0/3 ! Lorenz's parameters double precisi…

太陽を周回する地球の運動について見てみる。 解くべき運動方程式は時刻tでの軌道面における地球の直交座標を(x(t),y(t))として、 である。太陽は原点から動かないものとする。GMは万有引力定数と太陽質量の積である。 program ode use runge_kutta_4 implic…

RK4モジュールを使って落体の運動方程式。 せっかく計算本体を別モジュールにしたにも関わらず記述量が変わっていないのがなんとも。 program ode use runge_kutta_4 implicit none double precision, parameter :: g = 9.80665 ! gravitational acceleratio…

RK4を使って一階連立微分方程式を解くためのモジュール。 サブルーチンset_funcでのを設定し、 関数nextで時刻tの値から時刻t+hの値を算出する。 FORTRAN 2003以降でないとコンパイルできないと思う。 module runge_kutta_4 implicit none interface functio…

落体の運動方程式をRK4により数値的に解くFORTRANプログラム。 Euler法によるプログラムのうちy(t)とv(t)からy(t+h)とv(t+h)を求める漸化式部分をRK4に変更しただけ。 同じ計算をnextyとnextvで行っているのでおいしくないプログラムになっている。 もっとま…

落体の運動方程式をEuler法により数値的に解くFORTRANプログラム。 program ode implicit none double precision, parameter :: g = 9.80665 ! gravitational acceleration double precision, parameter :: t0 = 0 ! initial t double precision, parameter …

定番のclassical 4th-order Runge-Kutta法、略してRK4を使ってみる。 ...snip double precision function nextx(t, x, h) double precision, intent(in) :: t, x, h double precision :: k1, k2, k3, k4 k1 = f(t, x) k2 = f(t + h / 2, x + h / 2 * k1) k3 …

Heun法を使ってみる。 ...snip double precision function nextx(t, x, h) double precision, intent(in) :: t, x, h double precision :: k1, k2 k1 = f(t, x) k2 = f(t + h, x + h * k1) nextx = x + h / 2 * (k1 + k2) end function nextx ...snip 計算し…

Euler法よりもマシなはずの中点法で計算してみる。 Euler法との違いは漸化式のの定義部分だけである。 ...snip double precision function nextx(t, x, h) double precision, intent(in) :: t, x, h double precision :: k1, k2 k1 = f(t, x) k2 = f(t + h /…

関数についての常微分方程式の初期値問題 の解析解はである。 これをEuler法で数値解析するプログラムをFORTRAN、古流書式でなく90/95で書いてみる。 微分方程式の初期値問題が の形の場合、刻み幅をとし、 と定義すると、Euler法は、 と書けるので、素直に…

の導関数はなので、 x=0のときy'=0となり、ニュートン・ラフソン法においてはx=0は初期値として不適である。 ところが2階中心差分による微分係数ではでであり、 誤って不適切なパラメータが与えられると初期値x=0でも不幸にも正しい解に収束してしまったりす…

はさみうち法でやってみる。 PROGRAM INTERS IMPLICIT NONE DOUBLE PRECISION ERR,HIGH0,LOW0,Y1,Y2,X,Y,HIGH,LOW,MID PARAMETER(ERR=1.0D-10,HIGH0=1.0,LOW0=0.0) Y1(X)=X*TAN(X) Y2(X)=SQRT(1-X**2) HIGH=HIGH0 LOW=LOW0 DO MID=(HIGH+LOW)/2 Y=Y2(MID)-Y1…

導関数を与えて微分係数を算出する方法。 PROGRAM INTERS IMPLICIT NONE DOUBLE PRECISION Y,YD,Y1,Y2,X,XN,E,X0 PARAMETER(E=1.0D-12,X0=0.9) Y1(X)=X*TAN(X) Y2(X)=SQRT(1-X**2) Y(X)=Y1(X)-Y2(X) YD(X)=(TAN(X)+X/COS(X)**2)-(-X/Y2(X)) X=X0 DO XN=X-Y(X…

文関数で定義すると短く済む。 どんどん現代書法的に非推奨になる。いい感じ。 PROGRAM INTERS IMPLICIT NONE DOUBLE PRECISION Y,YD,Y1,Y2,X,XN,E,H,X0 PARAMETER(E=1.0D-12,H=1.0D-8,X0=0.9) Y1(X)=X*TAN(X) Y2(X)=SQRT(1-X**2) Y(X)=Y1(X)-Y2(X) YD(X)=(Y…

古流記法なFORTRANプログラム。 x-y直交座標系において原点を中心とする単位円と曲線のにおける交点を求める。 曲線のy座標は負の値をとらないので、単位円（の半分）はと書ける。 求めるべきはとなるxである。 普通にニュートン・ラフソン法を使う。 微分係…