Mercurial > op > rk_sakamoti

comparison Lorenz2Langford.diff @ 6:d90925ba642e default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
add a patch file
author "uncorrelated zombie" <uncorrelated@yahoo.co.jp>
date Tue, 23 Jul 2024 02:36:26 +0900
parents
children
comparison
equal deleted inserted replaced
5:0ca8347e4c47 6:d90925ba642e
1 diff -r 925250ea21a4 rk_Lorenz.R
2 --- a/rk_Lorenz.R Sun Jul 21 06:20:46 2024 +0900
3 +++ b/rk_Lorenz.R Tue Jul 23 02:35:22 2024 +0900
4 @@ -1,25 +1,25 @@
5 dll <- paste("rksub", .Platform$dynlib.ext, sep="")
6 dyn.load(dll)
7
8 -max_i <- as.integer(5000)
9 +max_i <- as.integer(5000*10)
10
11 r <- .Fortran("rk_sakamoti",
12 - as.double(c(0.5, 0.5, 0.5)),
13 - max_i,
14 + as.double(c(1.0, 1.0, 1.0)),
15 + max_i,
16 matrix(0, max_i, 4))
17 dyn.unload(dll)
18
19 m <- r[[3]]
20 colnames(m) <- c("t", "x", "y", "z")
21
22 -png(filename = "Lorenz system.png",
23 +png(filename = "Langford.png",
24 width = 800, height = 600, bg="white", type="cairo")
25 library(scatterplot3d)
26 scatterplot3d(m[,2], m[,3], m[,4],
27 highlight.3d = FALSE,
28 col.axis = "blue", col.grid = "lightblue",
29 type = "l", color = "purple", lwd = 2,
30 - main = "Lorenz system",
31 + main = "Langford Attractor",
32 xlab = colnames(m)[2],
33 ylab = colnames(m)[3],
34 zlab = colnames(m)[4],
35 @@ -27,3 +27,16 @@
36 cex.lab = 2)
37 dev.off()
38
39 +subplot <- function(m, n){
40 + m2 <- m[n, ]
41 + scatterplot3d(m2[,2], m2[,3], m2[,4],
42 + highlight.3d = FALSE,
43 + col.axis = "blue", col.grid = "lightblue",
44 + type = "l", color = "purple", lwd = 2,
45 + main = "Langford Attractor",
46 + xlab = colnames(m)[2],
47 + ylab = colnames(m)[3],
48 + zlab = colnames(m)[4],
49 + cex.main = 3,
50 + cex.lab = 2)
51 +}
52 diff -r 925250ea21a4 rkmod.f90
53 --- a/rkmod.f90 Sun Jul 21 06:20:46 2024 +0900
54 +++ b/rkmod.f90 Tue Jul 23 02:35:22 2024 +0900
55 @@ -10,10 +10,15 @@
56 real(kind=dpkind) :: time
57 real(kind=dpkind),dimension(:) :: wk
58 real(kind=dpkind),dimension(size(wk)) :: rhside
59 - real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0
60 - rhside(1)=a*(wk(2)-wk(1))
61 - rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3)
62 - rhside(3)=wk(1)*wk(2)-c*wk(3)
63 +! real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0
64 +! rhside(1)=a*(wk(2)-wk(1))
65 +! rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3)
66 +! rhside(3)=wk(1)*wk(2)-c*wk(3)
67 + real(kind=dpkind),parameter :: a = 0.6
68 + rhside(1) = (wk(3) - 0.7)*wk(1) - 3.5*wk(2)
69 + rhside(2) = 3.5*wk(1) + (wk(3) - 0.7)*wk(2)
70 + rhside(3) = a + wk(3) - wk(3)**3 / 3 - (wk(1)**2 + wk(2)**2)*(1 + 0.25*wk(3))
71 +
72 end function
73 !Rから呼び出すサブルーチン
74 subroutine calc(x, max_i, r)