-
Notifications
You must be signed in to change notification settings - Fork 0
/
lorenz_attractor.f90
82 lines (67 loc) · 1.7 KB
/
lorenz_attractor.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
program main
implicit none
real*8::h,y,x,z,K,L,u_limit,M,t
real*8::sigma,rho,beta
real*8::x1,y1,z1
!================================
!constants
sigma=10.0
rho=28.0
beta=8.0D0/3
t=0
h=0.003333
u_limit=200
! initial conditions
x=10.0
y=-0.01
z=9.0
!second point
x1=10.00001
y1=-0.01
z1=9.0
OPEN(20,file='lorenz_attractor.dat')
do while(t<=u_limit)
CALL KLM(x,y,z,f,g,h,K,L,M)
x=x+M
y=y+K
z=z+L
CALL KLM(x1,y1,z1,f,g,h,K,L,M)
x1=x1+M
y1=y1+K
z1=z1+L
t=t+h
write(20,*)x,y,z,x1,y1,z1
enddo
CLOSE(20)
contains
real*8 function f2(x,y,z) !derivative of x
real*8::x,y,z
f2 = sigma*(y-x)
end function f2
real*8 function f(x,y,z) !derivative of y
real*8::x,y,z
f = rho*x-x*z-y
end function f
real*8 function g(x,y,z) !derivative of z
real*8::x,y,z
g = x*y-beta*z
end function g
subroutine KLM(x,y,z,f,g,h,K,L,M) !subroutine for Runge Kutta (order 4)
real*8::x,f,g,y,z,k1,k2,k3,k4,h,K,L,l1,l2,l3,l4,M,m1,m2,m3,m4
m1=h*f2(x,y,z)
k1=h*f(x,y,z)
l1=h*g(x,y,z)
m2=h*f2(x+h/2,y+k1/2,z+l1/2)
k2=h*f(x+h/2,y+k1/2,z+l1/2)
l2=h*g(x+h/2,y+k1/2,z+l1/2)
m3=h*f2(x+h/2,y+k2/2,z+l2/2)
k3=h*f(x+h/2,y+k2/2,z+l2/2)
l3=h*g(x+h/2,y+k2/2,z+l2/2)
m4 = h*f2(x+h,y+k3,z+l3)
k4 = h*f(x+h,y+k3,z+l3)
l4 = h*g(x+h,y+k3,z+l3)
K = (k1+k4)/6 + (k2+k3)/3
L = (l1+l4)/6 + (l2+l3)/3
M = (m1+m4)/6 + (m2+m3)/3
end subroutine KLM
end program main