基礎セミナー課題
080230547 亀谷直樹
1.
ソース
dA/dt=Fin
– A/ι
Euler法で考えると、A(n+1)= A(n)+⊿t (Fin-A(n)/ι)。
dAをdx、Finをｆ(=10.0)、ιをlambda(=1.0)として、計算する。
1: PROGRAM exponential_decay
2: ! t02.f90 --
3: ! solve the differential
equation dx/dt =
f – x/(lambda)
4: IMPLICIT NONE
5: REAL :: lambda = 1.0
6: REAL :: x = 1.0
7: REAL :: t = 0.0
8: REAL :: f = 10.0 ! 定数を追加
9: REAL :: dxdt
10: REAL :: timespan
= 10.0
11: INTEGER, PARAMETER :: n =
10000 ! 値を変えてdtを変化
12: INTEGER :: i
13: REAL :: dt
14: OPEN
(UNIT = 15, FILE = "decay.dat", STATUS =
"REPLACE")
15: 110
FORMAT(1X, 2F10.5)
16: dt = timespan / real(n)
17: DO i = 1, n
18: dxdt = f
- x / lambda
19: x = x + dxdt
* dt
20: t = t + dt
21: WRITE (UNIT = 15, FMT = 110) t , x
22: END
DO
23: CLOSE
(UNIT = 15)
24: END
PROGRAM exponential_decay
2.
実行結果と考察
gnuplotの図を挿入できなかったので、最終計算結果だけをのせる。
厳密解は、ι*Fin =
10.0*1.0 = 10.0。
nの値を変えてdtを変化させた。
|
nの値 |
dt |
計算回数 |
最終結果 |
|
10 |
1 |
10 |
10.00000 |
|
100 |
0.1 |
100 |
9.99976 |
|
1000 |
0.01 |
1000 |
9.99961 |
|
10000 |
0.001 |
10000 |
9.99952 |
t = 100.0に変えて実行してみた。
|
nの値 |
dt |
計算回数 |
最終結果 |
|
100 |
1 |
100 |
10.00000 |
|
1000 |
0.1 |
1000 |
10.00000 |
|
10000 |
0.01 |
10000 |
9.99995 |
|
100000 |
0.001 |
100000 |
9.99952 |
dtが大きいほど、誤差が大きくなるという、おかしな結果になってしまった。
ただし、計算回数が、多くなると、精度はあがるようだ。