2 ver..2

Transcription

1 Title2 年 度 カオス フラクタル 講 義 ノート Author(s) 井 上, 純 一 Citation Issue Date 2 DOI Doc URLhttp://hdl.handle.net/25/46977 Right Type learningobject Additional There are other files related to th Information above URL. File Information CF2_LectureNote.pdf ( 講 義 ノート 一 括 ダウンロ Instructions for use Hokkaido University Collection of Scholarly and

2 2 ver..2

3 2 3. ver..2,, 29, ver.., 2, 2 2 2,.,,,,,, ver.2.x. 2 OCW, HUSCAP : ( ), ( ),,.,, ( ).,, ( ).,,,, ( V:. HP ), (,, ).,, ( ).,,.,,, /. 23 8

4 2? 6.,, : : / ( ) x n+ = 4x n ( x n ) L 4 (x) : L a 4 (x) :,

5 L a (x) n a , : : Neugart : : : : /

6 : q

7 2,,,,,..?,.,,,,,,..,,,. -, 3,.,,.., x, y.,, x, x.8.6 y x : x y. ( ). y.,,?,,,,,,..2,.4.. 6

8 2, ( 2 ). -exp(-x) y x 2: x y.,,.,, ( 3 )., ( ) ( )., x, y y = ax () ( a > )., 2 3,., 2 y = exp( ax) (2), 3 y = ax exp( x) (3)., x = y = exp( ax) ax + a2 2 x2 + (4) y = ax exp( x) ax ax 2 + (5), x, ()., 7

9 2.4 x*exp(-x) y x 3: x y.,,., x. 2., ((4)(5) 2 )., A( x A ) B( x B ),., y A+B, () y A+B = a(x A + x B ) = ax A + ax B = y A + y B (6), A B., (2) (3)., (2) y A+B = exp[ a(x A + x B )] exp( ax A ) + exp( ax B ) = y A + y B (7), (7),.,, ( c > )/ ( c < ), x ( ) cx., () y = a(cx) = cax = cy (8), c., (2) y = exp[ a(cx)] c{ exp( ax)} = cy (9), (9), c., x, y, (6) (8), 2 (x) ( ).,,, ( ) ( ),.,,. 8

10 2.,,.,, ((2) (3) ) ( ) (), ( )., x,,,..,.., :,.,,. ( ) 4( )., 4( ). 4: ( )..,, (99) p ,,.., K +, Na + Cl,, ( 7 mv )., 4( ) ( 55 mv) V ms, ( )., ms,,

11 2 S w S w N S N 5: ( ). S i (i =,, N) S i = ( ) S i = ( ). S i w is i w S θ S = S =., S i = ( ) S i = ( ) 2, h w S + w 2 S w N S N w S () w (w, w 2,, w N ) () S (S, S 2,, S N ) (2). S h h > θ : S (3) h < θ : S (4)., Θ(x) Θ(x) = { (x ) (x < ) (5), θ = ( w ) S S = Θ(h) = Θ(w S) (6), h S., S i h i = w i S i, S : S (i) = ah i, N S S = a(h + h h N ) = ah + ah ah N = S () + S (2) + + S (N) (7)

12 2, S., (5)(6) (7).,,., N h N i= h i ɛ S.,, S S = a(h + ɛ) = ah + aɛ = S + aɛ S (8), ɛ a aɛ., S = Θ(h + ɛ) (9), h ɛ, S = Θ(h + ɛ) = Θ(h) = S, h ɛ, h > ɛ S = Θ(h + ɛ) = Θ(h) = S, S.,,,.,,,., V: : ( ),,.,.,. 2 :, l m, g. ( ) θ d(l θ) m dt = mg sin θ (2) d 2 θ dt 2 = g sin θ (2) l

13 2., θ, θ = d 2 θ dt 2 = g l θ (22),, θ = e iλt, λ : λ 2 + (g/l) = λ = ±i g/l, α, β θ = α e i g/lt + β e i g/lt g g g g = (α + β) cos t + i(α β) sin l l t = A sin l t + B cos l t (23) (α = (A ib)/2, β = (A + ib)/2)., θ, sin θ θ =, d 2 θ dt 2 = g ) (θ θ3 l 3! + θ5 5!, θ ( θ 3 /3! + θ 5 /5! ).,,., θ = e iλt (24), λ λ 2 = g ) ( θ2 l 3! + θ4 5!, θ (t ) λ, ( θ 2 /3! + θ 4 /5! )., θ = e iλt. (2),,,.,.,.. (24) (25) 3, N. N, N(t)., N, N dn dt = λn (26). λ.,. N, N N dn N = λ t dt (27) 2

14 2, log(n/n ) = λt, N(t) = N e λt (28).,, (, ),,., (26). dn dt = λn µ N 2 = λn( µn) (29) λ, µ >.. µ, N. N, N : λµn 2.,,.,,, ( ),.,, 4.,,.,, ( ) ( ).,, 5. /,,,, / (, 2 ).,,,,.. 3. N,,,. N, : dn dt ( = µλn 2 ) µn (3) 4,,,. 5,,., HUSCAP,. 3

15 2, X = (/µn), dx/dn = /µn 2 dn dt = dx dt ( dx dn ) = µn 2 dx dt (3), X dx dt = λx (32), X = X e λt., N N(t) = µ( X e λt ) (33). t =, N = N, X N, X = (/µn ) N(t) = ) (34) µ (µ N e λt ,.., : ( d ) dt µn : ( = λ ) µn (35) dn dt = λn( µn) (36),., 2,. ( ) 3.2. (35) ( ) ( ) µn t µn t ( = λ ) µn t (37), ( /µn t ) ( /µn ), ( λ) ( ) ( = ( λ) ) ( = ( λ) 2 ) ( = ( λ) t ) µn t µn t µn t 2 µn (38) 6 (θ) (θ ). N N. 4

16 2, N t N t = ( )} (39) µ { ( λ) t µn., (34) (39) 6 9 cont. discrete 8 7 N t 6: ((34), ) ((39), ). µ =, N =, λ =..., µ =, N =, λ =.., (36). N t N t = λn t ( µn t ) (4) N t = ( + λ)n t µλnt 2 (4) 7., λ = a, µ = a/(a ), a >, N t = an t ( N t ) (42),. 7, dn/dt (N t N t t )/ t t =, dn/dt = N t N t. 5

17 2 (42), 2,,,,., ( ). (4), N t N =, λ =., µ = t = t =, (34)(39) ( 6)... :, ( I, II ) C.,.,,,,, ( URL ). 6

18 2, µ =, N =, λ =. µ =, N =, λ =. µ =., N = 2, λ =.5 µ =., N = 2, λ =.5 µ =., N = 2, λ =.5 µ =., N = 2, λ =. N t 7:. (λ, µ, N ) 7,,.,.,., λ =., µ =, N = N t = ( +.)N t.nt 2 (43), N = ( +.)N.N 2 = ( +.). = (44) N 2 = ( +.). = (45),.,,., 5 λ = a, µ = a/(a ), a >, a = ,,, ( (33)(39) ).., 2,., (42). 7

19 2.8.6 N Logistic map discrete-linear con.-linear t 8: λ = a, µ = a/(a ), a = 3.9 (cont-linear) (discrete-linear), (Logistic map).,,. 4, : N t = ( + λ)n t µλnt 2 (46),,.,., (46), N t = (a/µ( a))x t+,., x n : x n+ = ax n ( x n ) (47) (47) x = x., y = f(x) = ax( x) y = x, (x, y) = (x, ), (x, x ) (y = x x = x ), (x, f(x )) (y = x y = f(x) ), (f(x ), f(x )) (y = x x = f( ) ),, x, (47) x, x, x 2,,., a = 2.5, a = 4, a > a. 8

20 x(-x) x x(-x) x.62 path path.2 x n+.6.2 x n x n x n x n x n x n x n 3. x(-x) 3.9 x(-x) x x.8 path.8 path.6.6 x n+ x n x n x n 9:. y = f(x) = ax( x) y = x, (x, y) = (x, ), (x, x ) (y = x x = x ), (x, f(x )) (y = x y = f(x) ), (f(x ), f(x )) (y = x x = f( ) ),,. a = 2.5, 2.9, 3., 3.9. x =..., a = 3. x 2., a, a = 3.9, x. a = ,. 4.2 a = 3.9, x =., x x =. x =.., ( ).., n = 7, 8,,.,, a = 3.9., 9 a = 2.2, 2.9, (, ). 9

21 2.8.6 x n.4.2 x =. x = n : x n+ = 3.9x n ( x n ), x =. x = / ( ) a = 3.9., 2,., 2. x n+ = ax n ( x n ) (48) x n+ = ax n ax 2 n (49) 2 I, 2., (A ) (48),, (49) 3 / ),.,,.,, a = 3.9 = 39/, x =.5 = /2, (48) x = 39 4, x 2 = 52 6, x 3 = , (5),.,, 2

22 2,, 9., n (49)., (49) x = 39 4, x 2 = , (5),, x 2 2,, n (48).,,,,.,..,,.8.8 x.6.4 x n Type- Type Type- Type n.8.8 x.6.4 x n Type- Type n Type- Type-2 : (48)( Type-), (49)( Type-2). x =.5. n 67. n 67. 9, 32 (C int ), z 2 (b 3, b 3,, b ) z = 3 i= b i2 i, (b 3, b 3,, b ) = (,,, ),, z = = 2 32 = ,, (, ). 2

23 2 4.4 a = 3.9, a = 4.,,.,. : x n+ = f(x n ), x k f k (x), f k (x) = x, f k., f 6 (x) = x f 6, f 6 (x) = f 3 (f 3 (x)) = f 3 (x) = x (52), f 6, 3., p n, f n (x) = x x f p (x) = x x., f(x) = 4x( x), f k (x) = x x. a = 4., f(x) = 4x( x) (53) f 2 (x) = 4{4x( x)}( {4x( x)}) = 6x 8x x 3 64x 4 (54) f 3 (x) = (55), 2, x n+ = f(x n ) = 4x n ( x n ) -4 y = f,, f 4 y = x x.,, y = f 2 y = x y = f y = x 2 2 (2 )., f 2 (x) = x 2 x, x : 5x 8x x 3 64x 4 = (56). x =, 3/4, x(x 3/4), x 2 2., n n,., a = 4 x n+ = f(x n ) = 4x n ( x n ),, x n+ = 4x n ( x n ), x n+ = 4x n ( x n ) x n = sin 2 θ n sin 2 θ n+ = 4 sin 2 θ n ( sin 2 θ n ) = 4 sin 2 θ n cos 2 θ n = sin 2 2θ n (57) f k.. 22

24 x x x x 2: y = f,, f 4 y = x. f(x) = 4x( x). f 2, f 4, f, f 3.,, x n+ = 4x n ( x n ) sin 2 θ n+ = sin 2 2θ n (58)., sin 2 θ n+2 sin 2 θ n+2 = sin 2 2θ n+ sin 2 2θ n+ = 4 sin 2 θ n+ cos 2 θ n+ = 4 sin 2 θ n+ ( sin 2 θ n+ ) = 4 sin 2 2θ n cos 2 2θ n = sin 2 4θ n (59) sin 2 θ n+2 = sin 2 4θ n (6)., p, θ n+p sin 2 θ n+p = sin 2 2 p θ n (6) 2.,, m θ n+p = ± 2 p θ n + mπ (62) 2. 23

25 2., p (x n = sin 2 θ n = x n+p = sin 2 θ n+p ), l θ n+p = ± θ n + lπ (63). (62)(63), L l m θ n = Lπ 2 p ± (64). L, L =., (x =, 3/4), p =, θ n = π, π/3, x n+p = sin 2 θ n+p = sin 2 2 p θ n, x n = sin 2 π =, x n = sin 2 π/3 = ( 3/2) 2 = 3/4, OK., 2, p = 2 θ n = π/3, π/5, θ n = π/3, x n = sin 2 π/3 = ( 3/2) 2 = 3/4, 2. θ n = π/5 x n = sin 2 π/5 =.345,, (6) x n+p = sin 2 θ n+p = sin 2 2 p θ n, x n+ = sin 2 2π/5 =.94, x n+3 = sin 2 4π/5 = sin 2 (π π/5) = x n., , 2.,.. 2. : x n+ = 2x n (65). x =., n, x n 3.,,,. : 4/ =.9999, x n+ =.9999x n. 24

26 2 2 4 : x n+ = 2x n = { 2x n ( x n /2) 2 2x n (/2 < x n ) (66)., I #include <stdio.h> #include <math.h> #define x. /* */ #define N main() { FILE *pt; /* */ double x; int i; if((pt = fopen("tent_map.dat","wt"))!=null){ /* */ for(i =, x=x; i <= N; i++){ x =.-fabs( *x); /* */ fprintf(pt,"%d %lf\n",i,x); /* */ } } fclose(pt); /* */ }. n, n 3. Remark T (x) = 2x, x =. = /, x n+ = T (x n ) x = 2 = 4 5 = 5 x 2 = 2 5 = 3 5 = 2 5 x 3 = 4 5 = 5 = 4 5 x 4 = 8 5 = 3 5 = 2 5 x 5 = 4 5 = 4 5 = 4 5 = 4, y = 2x x [, ]. 25

27 2 Tent map y n n 3: x n., 2/5 =.4, 4/5 =.8 2., x =.3 = 3/,, x = 3/5, x 2 = 4/5, x 3 = 2/5, x 4 = 4/5, x 5 = 2/5,, 2/5, 4/5 2., x.γ (γ, 5), 2/5, 4/5 2, x =.γ γ 2 (γ 2, 5), x =.γ γ 2 γ N (γ N, 5) 2 5 N. 4 x =.9 x =. x = x n n 4: x =., x = , x = , x =. 2, N = 27 x = =

28 2.., x =., ,,.5, x n., x = , 2 = a = 4., a = 3.9 : x n+ = ax n ( x n ), : x, x, x 2,,., [, ]., x n [, ], x n., ( ).., ( n = ), x n,.,. for(k=; k<=999;k++){ count[k]=; /* */ } for(i =, x=x; i <= N; i++){ x =.-fabs( *x); for(k=; k<=999; k++){ if((x>=.*k) && (x<.*(k+))){ count[k]++; }else{ count[k]=count[k]; } } } for(k=; k<=999; k++){ 27

29 2 } } fprintf(pt,"%lf %lf\n",.*k,(double)count[k]*/n);, 5., N = x N = x 5: x n..,.,,. 5. B(x). x n+ = B(x n ) = { 2x n ( x n < /2) 2x n (/2 x n ) (67), x < /2 B(x) = T (x), /2 x B(x) = T (x) , 2. x, 2., x, x 2 x = α 2 + α α α n 2 n + (68) x = (.α α 2 α n ) 2 (69) 28

30 2 5.,, x =.5.5 = n + (7), 2 x = (.) 2, x = = n + = 2 = (7) 2, 2 x = (. ) 2., x (68)., x =.29. 2,, 2, (α ),,.,., 2,,, α., α.., 2. x = =.48, x = (.) 2,.48 2 =.836 x = (.) 2., =.672 x = (.) 2,, 2,.672 = , =.244 x = (.) 2, =.4288 x = (.) 2 2.,,, 2 = = (...) 2 (72) (68), x < /2, α =, x x = α α αn 2 n + = (.α 2α 3 α n ) 2 (73)., x x < /2 B(x) = 2x ( α2 B((.α 2 α 3 α n ) 2 ) = B α ) αn 2 n + ( α2 = α ) αn 2 n + = α α α n 2 n + α n+ 2 n + = (.α 2α 3 α n ) 2 (74). x = (.α 2 α 3 α n ) 2, B(x) B 2 (x) B 2 ((.α 2 α 3 α n ) 2 ) = B((.α 2 α 3 α n ) 2 ) ( α2 = 2 + α α n 2 n + α ) n+ 2 n + = α α α n 2 n 2 + α n+ 2 n + = (α 2.α 3 α 4 α n ) 2 (75) 5 ( ) ( ),. 29

31 2., α 2 =, x >, α 2 =,, B((.α 2 α 3 α n ) 2 ) = (.α 3 α 4 α n ) 2., 2 (, )., /2 x, B(x) = 2x, x x = (.α α 2 α n ) 2 (α = ),. ( α B((.α α 2 α n ) 2 ) = B 2 + α α ) n 2 n + ( α = α α ) n 2 n + = α + α α α n 2 n + = (α ) + α α α n 2 n + = (.α 2α 3 α n ) 2 (76)., α =.,, , x n+ = f(x n ) x, x,, x n, [, ] [ɛ, ɛ 2 ] x, x,, x n,.., [, ] 4, I = [, /4), I 2 = [/4, /2), I 3 = [/2, 3/4), I 4 = [3/4, ], /4 = (.) 2, /2 = (.) 2, 3/4 = (.) 2, 2 2 I (.) 2, I 2 (.) 2, I 3 (.) 2, I 4 (.) 2., 4., x = / 2., 2 x = 2 = = (...) 2 (77),,. 2,, I, I 2, I 3, I 4. x = (...) 2 I 3 B(x ) = (...) 2 I 2 B 2 (x ) = (...) 2 I 4 B 3 (x ) = (...) 2 I 3 3

32 2 B 8 (x ) = (...) 2 I., 4,,, (, I = 7, I 2 = 9, I 3 =, I 4 = 6) 4., [, ],,, ( ). 5.3, T (x) B(x) T (x) = B(x).,, /2 x x = (.α α 2 α 3 α n ) 2 T (x), 2 (. ) 2 T ((.α α 2 α 3 α n ) 2 ) = B((.α α 2 α 3 α n ) 2 ) = (. ) 2 (.α 2 α 3 α n ) 2 = (.( α 2 )( α 3 ) ( α n ) ) 2 (.β 2 β 3 β n ) 2 (78)., β k α k (k = 2, 3, ), α k β k, α k β k. (, x < /2, x = (.α 2 α 3 α n ) 2, T (x) = B(x), T ((.α 2 α 3 α n ) 2 ) = (.α 2 α 3 α n ) 2.), /2 x T (x) 2,,., [, ],, [, ],.,, 4,, I (.) 2, I 2 (.) 2, I 3 (.) 2, I 4 (.) 2.,. 5.4 : x n+ = 4x n ( x n )?, 6.,,,, x n =, x n., x n+ = L 4 (x n ) 6 :x n = sin 2 θ n, sin 2 θ n+ = sin 2 2θ n, a = 4 { θ n+ = 6, L a(x) = ax( x). 2θ n ( θ n < π/4) (79) π 2θ n (π/4 θ n π/2) 3

33 N = x N = x 6: :x n+ = 4x n ( x n ) x n. [, ]..,, y n = (2/π)θ n y n, y n ( ).,, (a = 4)., B(x) T (x) L 4 (x).,, L 4 (x) x x = x + δx B(x), T (x). L 4 (x). 5.5, [, ], P (x) P (x) =, P (x)dx = (8)., : x, x 2,, x n,, P (x) = lim n n n δ(x x i ) (8) 7., P (x) =, : L 4 (x).,. i= 7 P (x)., (8) (8). n P (x)dx = lim δ(x x i )dx = lim n n n n n = i= 32

34 2 a = 4, x n = sin 2 θ n, θ n = (π/2)y n, y n, P (y) =., x n = sin 2 θ n θ n dx n dθ n = 2 sin θ n cos θ n = 2 x n ( x n ) (82)., θ n = (π/2)y n y n, dθ n /dy n = (π/2), dx n = 2 x n ( x n ) dθ n = 2 x n ( x n ) π 2 dy n (83) dy n = dx n π x n ( x n ) (84)., P (y)dy x P (y)dy = dy = = dx π x( x) = P (x)dx (85) P (x) = π x( x) (86). 7,., N = (./3.4)*(./sqrt(x*(-x))) numerical analytical x n x x n 7: x n+ = 4x n( x n) x n x n ( ). (86). x n, x n+., : y = 4x( x) (, 5). a = 4., g(x) P (x) lim n n n g(x i ) = i= g(x)p (x)dx (87) 33

35 2 8., a 4 (8). 8 a( 4)., P (x). ( ).,,., 4,,, I I 2 I 3 I 4 I 4,,,.,.,,,.,,. 4 8 P(x) a = 3.58 P(x) a = x x a = P(x) 25 2 P(x) 3 25 a = x x 8: : x n+ = L a (x) = ax( x). (8)., [, ]. 3 8,. 34

36 2 [, ] 8, / 3., / ,. 35

37 2 3 [, ] 8, / 3.,.. #include<stdio.h> #include<math.h> #include<stdlib.h> #define x (./sqrt(3.)) #define N 8 int count[2]; main() { FILE *pt; double x; int i,bin,bin2,bin3,bin4,bin5,bin6,bin7,bin8; bin=bin2=bin3=bin4=bin5=bin6=bin7=bin8=; printf("x=%lf\n",x); if((pt = fopen("test3.dat","wt"))!=null){ fprintf(pt,"%d.",); for(i =,x=x; i < N; i++){ x = 2.*x+.; if(x<.){ x = x; printf("%lf\n",x); count[i]=; fprintf(pt,"%d",); }else if(x>=){ printf("%lf\n",x); x = x-.+.; count[i]=; fprintf(pt,"%d",); } } } fclose(pt); for(i=; i< N; i++){ if((count[i]==) && (count[i+]==) && (count[i+2]==)) bin++; }else if((count[i]==) && (count[i+]==) && (count[i+2]==)){ bin2++; 36

38 2 }else if((count[i]==) && (count[i+]==) && (count[i+2]==)){ bin3++; }else if((count[i]==) && (count[i+]==) && (count[i+2]==)){ bin4++; }else if((count[i]==) && (count[i+]==) && (count[i+2]==)){ bin5++; }else if((count[i]==) && (count[i+]==) && (count[i+2]==)){ bin6++; }else if((count[i]==) && (count[i+]==) && (count[i+2]==)){ bin7++; }else if((count[i]==) && (count[i+]==) && (count[i+2]==)){ bin8++;} } printf("%d %d %d %d %d %d %d %d\n",bin,bin2,bin3,bin4,bin5,bin6,bin7,bin8); } 2 8, 2., , 2.,

39 ,, L a (x), L 2 a(x) = L a (L a (x)),, L n a(x) y = x, [ ] [ ] [ ].,, a = 2 : L 2 (x) = 2x( x), L 2 2(x), L 2 3(x), L 3 4(x) x 9., a = 4.8 n = 2 x.8 n = 4 x x x.8 n = x.8 n = 3 x x x 9: y = L 2,, L4 2 y = x. L 2(x) = 2x( x). 2 p.4 7, a = 4 L n 4 (x)., n L n a(x) ( ),, w a (n). p.4 7, a = 4 w 4 (n) n, w 4 (n) = 2 n., a = 2 9, n w 2 (n) = 2., 38

40 2. dl n 2 (x)/dx = x = /2. : d df(x) {f(g(x))} = dx dg(x) dg(x) dx (88), dl n 2 (x) dx n 2 (x) n 2 2 (x) = dl n 2 (x) dl n 2 (x) dl dl2 n 2 (x) dl dl n 3 2 (x) dl2 2(x) dl 2 (x) dl 2(x) dx = 2( 2L n 2 (x)) 2( 2L n 2 2 (x)) 2( 2L 2 (x)) 2( 2x) = (89)., L 2 (x) = L 2 2 = = L2 n (x) = /2 x x = /2, L n 2 (x) = n x = /2 9, L n 2 n w 2 (n) = , : h a. log w a (n) h a = lim n n (9), a = 2, a = 4 log 2 h 2 = lim n n = (9) log 2 n n log 2 h 4 = lim = lim = log 2 (92) n n n n,.,.,,,., n. 7 2, x n+ = f(x n ) x x x ( ) x x = tan θ = f (x ) (93) 9., L (n) 2 (/2) = /2 n. L 2(x) = /2 4x 2 4x + = (2x ) 2 = n =. L n 2 (/2) = /2, Ln+ 2 (/2) = L 2 (L n 2 (/2)) = L 2(/2) = /2 n +., L (n) 2 (/2) = /2 n. x = /2 L 2 (x) = /2 4x 2 4x + = (2x ) 2 = x = /2., L 2 2 = L 2(L 2 (x)) = /2 L 2 (x) L 2 (x) = /2, L 2 (x) = /2 x /2,, L 2 2 (x) = /2 x /2.. 39

41 2 x = f (x ) x (94)., f (x ) = (df(x)/dx) x=x. x n+ θ x x x x+ x x n 2: x = x x = x.,, x x 2 ( ), x 2 = f (x ) x = f (x ) f (x ) x (95)., n x n = f (x n ) f (x n 2 ) f (x ) f (x ) x (96)., n, x f (x n ) f (x n 2 ) f (x ) f (x ) x ( ).,, ( ).,, f (x i ) δ (δ > ) x n δ n x = e n log δ x (n ) (97),., log δ. λ = ( ) n log xn x } {{ } = n log f (x n ) f (x n 2 ) f (x ) f (x ) = n { } n log f (x n ) + + log f (x ) = lim log f (x i ) (98) n n i= 4

42 2 λ,., x n / x x n x = e nλ (99), λ >, x n / x n e nλ. λ <, x,, 2. λ.,,,., : x n+ = 2x n = { 2x n ( x n /2) 2 2x n (/2 < x n ) (), x f (x) = T (x) = 2, n λ = lim log T n log 2 (x i ) = lim = log 2 > () n n n n i=,., : x n+ = B(x n ) = { 2x n ( x n < /2) 2x n (/2 x n ) (2), f (x) = T (x) = 2, λ = log 2, L 4 (x) a = 4 L 4 (x) = 4x( x), L 4(x) = 4( 2x), x., (98) 4

43 2, P (x),. n λ = lim log f (x i ) = n n i= dxp (x) log L 4(x) = log 4( 2x) dx π x( x) (3). x. x = /2, x = sin 2 θ 2 /2 log 4( 2x) dx π x( x) π/4 log(4 cos 2θ) = 2 2 sin θ cos θ π sin θ cos θ = log π π/4 log(cos 2θ)dθ log 4 + I (4). 2,,. I = 2 π π/2 log(cos θ)dθ = 2 π π/2 log(sin θ)dθ (5), : θ π/2 θ. 2 3 I = 2 2 π = 2 2 π = 2 2 π = 2 2 π { π/2 log(cos θ)dθ + π/2 π/2 π/2 = 2 log π log(cos θ sin θ)dθ ( ) sin 2θ log 2 π/2 {log(sin 2θ) log 2} π log(sin θ)dθ } log(sin θ)dθ = 2 log 2 + I 2 (6). 2θ θ., I, I = log 2., L 4 (x),,. λ = log 4 log 2 = log 2 (7) 7..4 : L a 4 (x).,,. 42

44 2, (98) 2., L a 4 (x) L a(x) = a 2x, n λ a = lim log(a 2x i ) (8) n n i=., x, x,., n, n. 2., n =, a.5 λ -.5 numerical log(2) a 2: L a (x). n =, a... ( 5, 6 )., a, a,., 2 [.82 :.85] 22( ),, 22( )., ( ) :,,., x ( x ), P (x) = δ(x x ), : δ(x x )f(x)dx = f(x ) λ = δ(x x ) log a 2x dx = log a 2x = log a 2(a ) a (9) 2 x, x,.,, ( ) (98).. 43

45 λ - λ a a 22: 2.., ax ( x ) = x, x = (a )/a. 23( )., (9).,, T : (x, x 2,, x T ), P (x) = T T δ(x x i ) () i=, λ = T T i= δ(x x i ) log a 2x dx = T T log a 2x i () i=.,, x (x, x 2,, x T ), T. 23( ) [, 4] a λ - -5 log(x*abs(-2*(x-)/x)) numerical λ a a 23: ( ). [, 4] a., a,., 44

46 2. 4 a,. T a (x) = a( 2x ), x n+ = T a (x n ) (2) a /2 a, a, 2. 45

47 カオス フラクタル 講義ノート 2 年度版 担当: 大学院情報科学研究科 井上 純一 課題 4 の解答例 まずは前回学んだリアプノフ指数の計算例を見ておこう. テント写像はパラメータ a を導入するこ とで次のように一般化することができる. Ta (x) = a( 2x ), xn+ = Ta (xn ) (3) リアプノフ指数の定義から, (3) 式に基づいて生成される軌道の各点での Ta (x) が求まればリア プノフ指数が計算できる. これは明らかに x の値に依らず 2a であるから数値計算するまでもなく n n log 2a log Ta (xi ) = lim = log 2a = lim n n n n i= λa (4) である. ということは, a を /2 a の範囲で変化させると必ずリアプノフ指数は正の値にな り, パラメータ a で一般化されたテント写像から生成される軌道はカオスとなる. また, a を変化し つつ, 写像からの軌道の取りうる値を重ね打ちすると興味深い現象が見られる. その結果をリアプ ノフ指数の a 依存性とともに図 24 に載せる. generalized tent map log(2*x) x a.8.9 図 24: 一般化されたテント写像の分岐図と対応するリアプノフ指数. 参考までにロジスティック写像 La6=4 (x) のリアプノフ指数を計算するためのプログラムの一部を 以下に載せておこう2. /* 写像が陽に与えられた場合のリアプノフ指数の数値計算 */ for(a=.; a <= 4.; a = a +.){ for(i =,sum=.,x=x; i <= N; i++){ // x = a*(.-fabs(-.999*x)); /* テント写像ではこちらを使う */ x= a*x*(-x); sum = sum + log(a*fabs(.-2.*x)); /* 写像の出力結果 x に対して和を計算 */ } 2 ここではファイルのオープン/クローズを使わないが, ファイル名を指定してリダイレクションすればよい. ここは 46 ページ目

48 2 printf("%lf %lf\n",a,(double)sum/n); } 8 L a (x), a,,,. a,, a.. 8. L a (x) a 25: a. 25 a. a : x n+ = ax n ( x n ), x n 22., 25. /* */ for(a=.; a <= 4.; a = a +.){ for(i =,sum=.,x=x; i <= N; i++){ // x = a*(.-fabs(-.999*x)); /* */ 22,,., ( ). x,. 47

49 2 } x= a*x*(-x); if(i >= Nskip){ printf("%lf %lf",a,x); } } 8.., a a = a = 3 2., x a, L a (x) = x x, x =,, x = (a )/a , a = 3 x = (a )/a x n+ = L a (x n ).8.6 map -(/a) x a 26: a x = (a )/a., a > 3 x = (a )/a, 2., x = (a )/a a > 3,.,, {dl a (x)/dx x=x }, { dl a (x) dx x=x } = { a } 2(a ) = 2 a > (5) a,, ( a)(3 a) >, a > 3., x = (a )/a a > 3 ( < a < 3 )., L a (x) = x x =,.,,., L a (x) x = ( ) {dl a (x)/dx x= } = a >, a > x =., a x = ( ). 48

50 a > 3 2, 2. L 2 a(x) = L a (L a (x)) = x. x, ( ax + a ) (x =, (a )/a ) L 2 a(x) x = x( ax + a ){a 2 x 2 a(a + )x + a + } = (6), a 2 x 2 a(a + )x + a + = x (+) 2 = a + + (a + )(a 3), x ( ) 2 2a = a + (a + )(a 3) 2a 2 2, a., x (±) 2 (7).8 map (a++sqrt((a+)*(a-3)))/(2*a) (a+-sqrt((a+)*(a-3)))/(2*a).6 x a 27: 2 a x (+) 2, x( ) 2. L 2 a(x) x (±) dl 2 a(x) dx = dl2 a(x) dl a (x) dl a(x) = a( 2L a (x)) a( 2x) (8) dx, x (±) 2., x (+) 2, L a (x (+) 2 ) = ax(+) 2 ( x(+) 2 ) = a + + (a + )(a 3) 2 = a + (a + )(a 3) 2a a (a + )(a 3) 2a,, a( 2L (+) a (x 2 )) = (a + )(a 3). = x ( ) 2 (9) a( 2x (+) 2 ) = (a + )(a 3) (2) 23 2, L a(x) L 2 a(x). 49

51 2, dl 2 a(x) dx x=x (+) 2 = ( (a + )(a 3) )( (a + )(a 3) + ) = (a 2 2a 4) (2)., a 2 2a 4 <,, (a + )(a + 6 )(a 3)(a 6) <., a 3 > a 6 <,, 3 < a < + 6 x (+) x ( ) 2,, L2 a(x ( ) 2 ) = x(+) 2 a( 2L 2 a(x ( ) 2 )) = (a + )(a 3) (22) a( 2x ( ) 2 ) = + (a + )(a 3) (23), (8) dl 2 a(x) dx x ( ) 2 = (a 2 2a 4) = dl2 a(x) dx x (+) 2 (24)., x ( ) 2 x(+) 2 3 < a < n 25, x (±) 2 a = + 6, a > , L 4 a(x) = x x =, (a )/a, x (±) 2, dl 4 a(x)/dx,. 25,, a,.,,., L 4 a(x) = x a(a(a(ax( x))( (ax( x))))( (a(ax( x))( (ax( x)))))) ( a(a(ax( x))( (ax( x))))( (a(ax( x))( (ax( x)))))) x = (25), x 6. x =, (a )/a, x (±) 2, x(x (a )/a)(x x (+) 2 )(x x( ) 2 ), (25), 2.,., 2 I ( ) : (a + )(a + 6 ). 5

52 a, a,., a.., x = x = (a )/a : a =, x = (a )/a, x = x (±) 2 : a = 3, x = x (±) 2 4 ( ) : a 2 = + 6., a 2 a, a a,, a 2 a = = =.225 r (26) a a 3 2. : (a 3 a 2 )/(a 2 a ) a 3, ( ).,, : (a n+2 a n+ )/(a n+ a n ) (26) r.,., a n+2 a n+ = r(a n+ a n ) (27) a 2 a = r(a a ) a 3 a 2 = r(a 2 a ) a a = r(a a 2 ). a = a, a a a = r(a a ) (28) a = a ra r = = 3.58 (29). a = ,. 4 8( ). 8..5, a > a 28., (4 3 )., 3, 25,.,, ( ),, II. 25, II ( ). 5

53 カオス フラクタル 講義ノート 2 年度版 担当: 大学院情報科学研究科 井上 純一 x x a a 図 28: ロジスティック写像のパラメータ a に対する軌道の分岐図. 図 25 の一部を拡大したもの. 8.2 参考: ニュートン法の復習 f (x) = a2 x2 a(a + )x + a + = の手の方程式 (x の多項式 = ) の数値解を求めたい場合, ニュートン法と呼ばれる手法を用いるとよい. この方法は方程式 f (x) = の解の候補 x = x で曲線 y = f (x) に接する接線の方程式: y = f (x ) + f (x )(x x ) (3) の x 軸との交点 : x x f (x ) f (x ) (3) を解候補 x の改良版とするもので, これを推し進めて n-番目の解候補と n + -番目の解候補の間 に成り立つ漸化式 xn+ = xn f (xn ) f (xn ) (32) を反復的に解き, この収束点を f (x) = の解とするものである. 注: 2 分法 というやり方でも解を求めることができる. 詳しくは を参照してみるとよい. レポート課題 5 2 次方程式: f (x) a2 x2 a(a + )x + a + = をニュートン法により求めるプログラムを作成せよ. このとき x の初期値 x を [, ] の一様乱数により与え, 初期値をいくつか変えることで全ての実解が 自動的に求まるようなプログラムを組むとなお良い. ここは 52 ページ目

54 2 a 3 < a < + 6,, a. x, dl 2 a(x)/dx,

55 2 5. /* */ #include <stdio.h> #include <math.h> double func(x,a) /* f(x)=. f(x) */ double x; double a; { return (pow(a*x,2)-a*(a+)*x+a+); } double Dfunc(x,a) /* f(x).. */ double x; double a; { return (2*a*a*x-a*(a+)); } main() { FILE *fpr; int i,imax=; double x,y,a; if((fpr = fopen("lower.dat", "wt"))!=null){ for(a=3.; a<=.+sqrt(6.); a=a+.){ /* a */ for(i =, x=.; i <= imax; i++){ x = x - func(x,a)/dfunc(x,a); if(fabs(y-x)<.e-5){ fprintf(fpr,"%lf %lf\n",a,x); break; }else{ y = x; } }} } fclose(fpr); }, a 2, 2 x,, x 54

56 2., a, a x, a x a a x,,. a x. for(a=3.,x=.; a<=.+sqrt(6.); a=a+.){ /* a */ for(i = ; i <= imax; i++){ x upper branch lower branch (a++sqrt((a+)*(a-3)))/(2*a) (a+-sqrt((a+)*(a-3)))/(2*a) a 29: 2. 9, : x n+ = ax n ( x n ) (33),., ( ), ( ). ( 5 ), 2, ( x ), ( a ),., ( ), ( ).,,,, 55

57 2. ( ),,, 26., 6,,, ( ).,. 9. :, ( ),., /,,., ( ),,,.,,,,,. (, ),, ( / ),,,.,, ( ).,, ( ), ( ),,,,.,,,. 958 Alban William Housego Phillips.,,.,,,,,, ,.,, Neugart (24), 2 26, ( ). 56

58 inflation (%) unemployment (%) 3:..,. 9.2 Neugart Neugart.,,.,. 9.2., t 27 U t, i.,, o t., t., t + U t+ ( ) t U t U t+ = U t + i( U t ) o t U t (34)., t + U t+, U t t : o t U t, : i( U t ). Neugart i,. : x t+ = ax t ( x t ) (35) 27 t. 57

59 s in Japan U**(-.263) s in Japan U**(-.7667) π.5 π U 99s in Japan U**(-.359) U 2s in Japan U**( )-.6 π.2 π U U 3:. a., i,., t o t. t. o t = J s + J c,t U t + d( U t ) (36) J c,t = γ(m π t ) (37), (36) J s + J c,t. 2 J s,., 2 J c,t J s,, Neugart (37) m π t. γ < γ <.,, J s + J c,t,,. ( t ), 2 d( U t ),, U t U t + d( U t )., d On the job searching., (36) t, 58

60 2., π e,t, w b,t, ( ) w p.,. Neugart, t +, t π t, a ( a ),. π e,t+ = aπ t + ( a)π e,t (38), π t,, w p, t w b,t.. w p = ( µ)y (39), y, µ.,, ( 32 ).,, /µ.6.5 /µ demand demand 32: /µ (demand).,..,,, (y),,.., y =., 28,.., t w b,t = ( b)u t (4)., b < b <, 29,,,,,, y = 3., t 28, ( ),, ( ). 29,. 3,,. 59

61 2. π t = δ ( π e,t + w ) b,t w p = ( π e,t + µ ( b)u ) t w p δ µ (4),, t w t w b,t w p w p (42)., w t w p, t w b,t w p., (4), t ( w t > ), t, π e,t,, ( w t < ),,.,,,,,,., (4) δ, (4) w p = w b,t, ( )., π t U t., (34) (36),(37) (4) (38) J s + γ(m π t ) U t+ = U t + i( U t ) U t U t + d( U t ) π t+ = ( µ δ µ + aπ t + ( a) ( b δ µ ( δπ t µ ( b)u t µ ( U t + i( U t ) U t J s + γ(m π t ) U t + d( U t ) )) )) (43) (44)., U t π t,., (43)(44) U = U t+ = U t, π = π t+ = π t U µ m(δ )( µ) = (45) b π = m (46)., π, m,, π = m (, (46) ), (43), U = U t+ = U t, π = π t+ = π t = m J s U + i( U ) U U + d( U ) = U (47) J s = i( U )(U + d( U )) U (48) J s., (48) (43)(44). 6

62 (43)(44) i,. 33 i =.8 U t π t U t.6 π t t t 33: i =.8 U( ) ( )., d =., b =.5, µ =.4, γ =.5, δ = 2, a =.5, m =.3. t = 2.,,,., i =. i =.85 U t, π t,. 34., i, d =., b =.5, µ =.4, γ =.5, δ = 2, a =.5, m =.3., i.399, i,., m d, U.8 π i i 34: i U( ) ( ). i =.8., d =., b =.5, µ =.4, γ =.5, δ = 2, a =.5, m =.3.,, U π. 6

63 i( ) U π , U π π U 35: U π ( )., d =., b =.5, µ =.4, γ =.5, δ = 2, a =.5, i =.8, m =.3.., 36,,, t, t +.,., Neugart, U t+.8.6 π t U t π t : i =.8 U( ) ( ) t, t+., d =., b =.5, µ =.4, γ =.5, δ = 2, a =.5, m =.3. 62

64 U t, π t,,, ( ) P (U), P (π). 37., Neugart P(U).4.3 P(π) U π :, P (U)( ), P (π)( ). d =., b =.5, µ =.4, γ =.5, δ = 2, a =.5, m =.3.,,,.,,,. ( ):,., Neugart,. HUSCAP. ( ) (, ),., HP: ( ): (6/7) ( ).,,.. 63

65 2 :., : dn dt = λn( µn) (49),.,,,,,.,,,.,,, (, ).,, ( ).,,,.. dx dt dy dt dz dt = y z (5) = x + ay (5) = b + z(x c) (52). a, b, c,,... : dy dt = f(t, y) (53)., f. h t n = t + nh, n =, 2, 3, (54) t (t ), t n+ = t + (n + )h, t n+ t n = h ( t n+ = t n + h)., h y(t n+ ) = y(t n + h) = y(t n ) + dy h + O(h 2 ) dt tn,y n = y(t n ) + f(t n, y n )h + O(h 2 ) (55) 64

66 2., y(t n ) = y n, y(t n+ ) = y n+ y n+ = y n + f(t n, y n )h + O(h 2 ) (56) t n = t + nh (57). (56), n (57) t n y n., h 2 (O(h 2 ) )..2 : dy/dt = f(t, y) (t n, y n ), (t n, y n ) (t n+, y n+ ) (t n+/2, y n+/2 ),., (54) t n+/2 = t + (n + /2)h y n+ = y(t n+ ) = y(t n+/2 + h/2) t n+ = t n+/2 + h/2, t n = t n+/2 h/2 (58) = y n+/2 + f(t n+/2, y n+/2 ) y n = y(t n ) = y(t n+/2 h/2) = y n+/2 f(t n+/2, y n+/2 ) ( ) h + ( ) f h (t n+/2, y n+/2 ) + O(h 3 ) 2 (59) ( ) h + ( ) f h (t n+/2, y n+/2 ) + O(h 3 )(6) 2., f df/dt., (59)(6) y n+ y n = f(t n+/2, y n+/2 )h + O(h 3 ) (6) y n+/2 = y(t n+/2 ) = y(t n + h/2) = y(t n ) + h 2 f(t n, y n ) = y n + k 2 (62), (53) y n+ = y n + k 2 + O(h 3 ) (63) k = hf(t n, y n ) (64) k 2 = hf(t n + h/2, y n + k /2) (65) t n = t + nh (66) 65

67 2, O(h 2 ) O(h 3 ) 3. (63)(64)(65)(66) 2.,,., O(h 5 ), 4, ( 4 ). y n+ = y n + (k + 2k 2 + 2k 3 + k 4 ) + O(h 5 ) 6 (67) k = hf(t n, y n ) (68) k 2 = hf(t n + h/2, y n + k /2) (69) k 3 = hf(t n + h/2, y n + k 2 /2) (7) k 4 = hf(t n + h, y n + k 3 ) (7) t n = t + nh (72)..2. :,., 4 (67)-(72) : dn dt = λn( µn) (73). 38., 32.,. #include <stdio.h> #include <math.h> #define h. #define mu. #define lambda.5 double func(y) double y; /* */ {return (lambda*y*(.-mu*y));} /* f */ main() { FILE *fpr; int i,imax=; double y,k,k2,k3,k4,k; 3 h. h =., O(h 2 ) O(h 3 ). 32 ( y(t) y n ). 66

68 N = 2, λ =.5, µ =. /(.-(.-.5)*exp(-.5*x)) N =, λ =., µ = /(-(-.)*exp(-.*x)) N = 2, λ =.5, µ =. /(.-(.-.5)*exp(-.5*x)) N t 38: 4. if((fpr = fopen("rk4.dat", "wt"))!=null){ for(i =,y=2.; i <= imax; i++){ k = h*func(y); k2 = h*func(y+.5*k); k3 = h*func(y+.5*k2); k4 = h*func(y+k3); k = (k+2.*k2+2.*k3+k4)/6.; y = y + k; fprintf(fpr,"%lf %lf\n",i*h,y); } } fclose(fpr); }.2.2 2: 3 (5)(5)(52).. /* */ #include <stdio.h> #include <math.h> #include <stdlib.h> #define h. #define a.5 #define b.5 67

69 2 #define c 5.7 double func(x,y,z) /* */ double x; double y; double z; {return (-(y+z));} double func2(x,y,z) /* 2 */ double x; double y; double z; {return (x+a*y);} double func3(x,y,z) /* 3 */ double x; double y; double z; {return (b+z*(x-c));} main() { FILE *fpr; int i,imax=5; double x,kx,k2x,k3x,k4x,kx,y,ky,k2y,k3y,k4y,ky,z,kz,k2z,k3z,k4z,kz; if((fpr = fopen("rossler.dat", "wt"))!=null){ for(i =,x=-.,y=-.,z=.; i <= imax; i++){ kx = h*func(x,y,z); ky = h*func2(x,y,z); kz = h*func3(x,y,z); k2x = h*func(x+.5*kx,y+.5*ky,z+.5*kz); k2y = h*func2(x+.5*kx,y+.5*ky,z+.5*kz); k2z = h*func2(x+.5*kx,y+.5*ky,z+.5*kz); k3x = h*func(x+.5*k2x,y+.5*k2y,z+.5*k2z); k3y = h*func2(x+.5*k2x,y+.5*k2y,z+.5*k2z); k3z = h*func3(x+.5*k2x,y+.5*k2y,z+.5*k2z); k4x = h*func(x+k3x,y+k3y,z+k3z); k4y = h*func2(x+k3x,y+k3y,z+k3z); k4z = h*func3(x+k3x,y+k3y,z+k3z); kx = (kx+2.*k2x+2.*k3x+k4x)/6.; ky = (ky+2.*k2y+2.*k3y+k4y)/6.; kz = (kz+2.*k2z+2.*k3z+k4z)/6.; x = x + kx; y = y + ky; z = z + kz; fprintf(fpr,"%lf %lf %lf %lf\n",i*h,x,y,z); 68

70 2 } } } fclose(fpr);,, 3 x, y, z , x, y, z. 2 x(t) y(t) 6 z(t) t t 39: 4. x(t), y(t) z(t). 39, xyz 3. 2 x, 3 y, 4 z., I gnuplot gnuplot gnuplot> set parametric gnuplot> splot lossler.dat u 2:3:4 w line ( gnuplot> gnuplot ) , 4 3 x-y, x-z 4. [ ]. 4, gnuplot. 69

71 2 a = b =.5, c = 5.7 z y x 5 2 4:. a = b =.5, c = 5.7. x 2 : dx/dt = v d 2 x dt 2 ɛ( x2 ) dx dt + x = (74) dv dt dx dt = ɛ( x 2 )v x (75) = v (76) x, v 34,. ( ), ɛ, x, v (gnuplot ). ( ) : (6/2) ( ) 3 WING,, 34. 7

72 y z x x 4: 4 3 x-y, x-z.,. 7

73 2 6 dv dt dx dt = ɛ( x 2 )v x (77) = v (78)., ɛ, ( ). x, v ( ) (77) v = dx/dt, v dv dt, E = x dx dt (79) ( ) d v 2 dt 2 + x2 = (8) 2 v 2 + x 2 = ( 2E) 2 (8), 2E = v 2 + x2 (x, v t = ). ɛ = (77)(78) 4, (x, v) ε = : Harmonic 3 2 v x 42: ɛ =. 2E = v 2 + x2 = = = 3.6. gnuplot. #include <stdio.h> #include <math.h> #include <stdlib.h> #define h. 72

74 2 #define epsilon. double func(x,y) double x; double y; {return (y);} double func2(x,y) double x; double y; {return (epsilon*(.-x*x)*y-x);} main() { FILE *fpr; int i,imax=5; double x,kx,k2x,k3x,k4x,kx,y,ky,k2y,k3y,k4y,ky,z,kz,k2z,k3z,k4z,kz; if((fpr = fopen("test.dat", "wt"))!=null){ for(i =,x=.,y=3.; i <= imax; i++){ kx = h*func(x,y); ky = h*func2(x,y); k2x = h*func(x+.5*kx,y+.5*ky); k2y = h*func2(x+.5*kx,y+.5*ky); k3x = h*func(x+.5*k2x,y+.5*k2y); k3y = h*func2(x+.5*k2x,y+.5*k2y); k4x = h*func(x+k3x,y+k3y); k4y = h*func2(x+k3x,y+k3y); kx = (kx+2.*k2x+2.*k3x+k4x)/6.; ky = (ky+2.*k2y+2.*k3y+k4y)/6.; x = x + kx; y = y + ky; fprintf(fpr,"%lf %lf %lf\n",i*h,x,y); } } fclose(fpr); } ɛ. 43 (, ) ɛ =., ɛ =.5,.5., ɛ =.5 x, v.,,.,.,,. 73

75 カオス フラクタル 講義ノート 2 年度版 4 担当: 大学院情報科学研究科 井上 純一 3 ε =. ε = v v x x ε = x, ε =.5 v, ε = v x t 2 図 43: 有限の の場合のアトラクタの様子と =.5 の場合の時間変化 (右).. レスラー方程式とその振る舞い 前回みたように, レスラー方程式: dx dt dy dt dz dt = y z (82) = x + ay (83) = b + z(x c) (84) をルンゲ クッタ法を用いて数値的に解き, その軌道を 3 次元に描画すると図 44 のようになる. こ の図を描画するのに, 方程式の初期値は原点近傍 (x, y, z) = (.,.,.) に選んである. この 図をみるからに複雑なアトラクタであることがわかるが, この振る舞いを方程式の形からある程度 説明することができる. まず, この図 44 から, 軌道がスタートしてからしばらくは x-y 平面内に停滞していることがわか る. そこで, 方程式 (82)(83)(84) において, x-y 平面内近傍での初期の振舞いがどうであるかを 調べるために, z ' とおいてみると, (82)(83) 式から dx dt dy dt = y (85) = x + ay (86) が得られる. (85) 式の両辺をを時間で微分し, それに (86) 式を代入すると, 次の x に関する 2 階 ここは 74 ページ目

76 2 a=.398, b=2, c=4 z x y 44:. a =.398, b = 2, c = 4.. d 2 x dt 2 adx dt + x = (87), x = e pt, p 2 ap + =, 2 p = a ± a (88), < a < 2 4 a 2 /2 ω, (87) x = e at/2 { αe iωt + βe iωt} = e at/2 {A cos ωt + B sin ωt} (89)., y (86) y = dx dt = eat/2 {( ) ( ) } aa ab 2 + Bω cos ωt + 2 Aω sin ωt (9) (A, B )., x-y e at/2 ( )., z,, 44,., 75

77 2 z,, (84)., x c, dz/dt = : z = b/(x c) >.,, x, x > c dz/dt = b + z(x c) > (b > ( 44 b = 2)), z ( b ). 45 z x y 2 45: 44, z. x > c z., x z., 45, z, x-y.., z (82) dx/dt,, x., z x > c, x < c, dz/dt = b + z(x c) < z,, x-y (, ),. 44 3,,., gnuplot gnuplot> set view 6,5 76

78 2 (6,5. ).,,. 3, 2., x-y x θ y cos θ x sin θ =, z ( 44 ) y cos θ x sin θ = z a=.398, b=2, c=4 6 5 z x y θ = θ = 6 θ = 2 θ = 8 θ = 24 θ = 3 z r 46: y cos θ x sin θ = z.., θ θ = θ = 3 o,,.,., (x t, x t+ ), 35., ( )., : 35. dx dt dy dt = ax + ay (9) = cx y xz (92) 77

79 2 dz dt = bz + xy (93) a, b, c a=, b=8/3, c=28 z x y 47: (9)(92)(93). a =, b = 8/3, c = a =, b = 8/3, c = 3, z = 2, 3, 4, (, ). 78

80 2 8,. a =, b = 8/3, c = 3 z = 2, 3, 4 48., 2 5 z = 2 z = 3 z = 4 y x 48: a =, b = 8/3, c = 3 z = 2, 3, 4. x = x + kx; y = y + ky; z = z + kz; arg = fabs(z-2); if((i>=iskip) && (arg <5.e-6)){ fprintf(fpr,"%lf %lf\n",x,y); }, z = 2 (x, y). 2, x-y 3 2.,. 79

81 カオス フラクタル 講義ノート 2 年度版 2. 担当: 大学院情報科学研究科 井上 純一 埋め込み次元とアトラクタの再構成 レスラー方程式: dx dt dy dt dz dt = y z (94) = x + ay (95) = b + z(x c) (96) を数値的に解いた後, 例えば x(t) のみに着目し, 適当な時間遅れ τ を導入し, 次のような m 次元ベ クトルを作る. (x(t), x(t + τ ),, x(t + (m )τ )) (97) そこで, この m 次元ベクトルを逐次プロットすることで, 前回みた奇妙なアトラクタが再構成できるか を調べてみよう. 図 49 に次元を m = 2 と選んだ場合の (x(t), x(t+τ )), (y(t), y(t+τ )), (z(t), z(t+τ )) をぞれぞれプロットした. この図より, x のみ, y, z のみの軌道を適切な時間遅れ τ, 埋め込み次元 6 3 τ=2 5 τ= y(t+τ) x(t+τ) x(t) y(t) τ= 5 z(t+τ) z(t) 図 49: 再構成. 埋め込み次元を m = 2 に選んでいる. m で切り出すと, x, y, z を 3 次元にプロットした場合とほぼ同じアトラクタを再構成することがで きる. この例では m = 2, τ =, 2 に選んでいるが, 与えられたデータが方程式からの時系列ではな く性質が不明なものである場合, これら 2 つのパラメータを適切に決定しなければならない. 時間 遅れ τ は x(t) と x(t + τ ) の自己相関関数, あるいは, 相互情報量を評価することで見積もることが できるが, ここではそれには触れない. 次元 m に関しては次に述べる 相関次元 がその指針を与 ここは 8 ページ目

82 2.,.,. 2.2,, 2.,.., ɛ ( )., n(ɛ). n(ɛ),. D = lim ɛ log n(ɛ) log ɛ (98), l D., l ɛ l/2ɛ, n(ɛ) = l/2ɛ. log(l/2ɛ) D 2 ( l ) = lim = lim ɛ log ɛ ɛ,. { log(l/2) log ɛ } + log ɛ log ɛ = (99), l n(ɛ) = l 2 /πɛ 2 log(l 2 /πɛ 2 ) D 2 ( l ) = lim = lim ɛ log ɛ ɛ { log(l 2 /π) log ɛ } + 2 log ɛ log ɛ = 2(2),.,.. 2,, 2, 3,,. [ ]. 2.3, {x(t i )}: i =, 2,, N τ, m : r i (x(t i ), x(t i + τ),, x(t i + (m )τ)), i =,, N (2)., r i ɛ, r j : j =,, N p i p i = N N Θ(ɛ r i r j ) (22) j= 8

83 2., Θ(x) ( ), x, x <., r i r j r i r j = x(t i ) x(t j ) + x(t i + τ) x(t j + τ) + + x(t i + (m )τ) x(t j + (m )τ) (23). N i= D 2 = lim p2 i ɛ log ɛ (24). N i= p 2 i = N 2 N i= j= (24) C(ɛ). N Θ(ɛ r i r j ) C(ɛ) (25) D 2 = lim ɛ log C(ɛ) log ɛ (26), ɛ, D 2 = log C(ɛ)/log ɛ, C(ɛ) = ɛ D2 (27). (25) ɛ C(ɛ), D 2. x(t) τ =. 2 = 2 m =, 2, 3, 4, 5 (x(t i ), x(t i + τ),, x(t i + (m )τ)) : i =,, N, 5., m, - log C(ε) m = m = 2 m = 3 m = 4 m = log ε 5:, x(t) τ = 2 m =, 2, 3, 4, 5 (x(t i ), x(t i + τ),, x(t i + (m )τ)) : i =,, N. N = 3. m = 2 D 2 2, m. 82

84 2 D 2 = 2.,, /. m = 5. /********************************************************************/ /* */ /********************************************************************/ #include <stdio.h> #include <math.h> #define h. #define a.398 #define b 2 #define c 4 #define imax 5 #define N 3 double g[n]; double g2[n]; double g3[n]; double g4[n]; double g5[n]; double xx[imax]; double func(x,y,z) double x; double y; double z; {return (-(y+z));} double func2(x,y,z) double x; double y; double z; {return (x+a*y);} double func3(x,y,z) double x; double y; double z; {return (b+z*(x-c));} double STEP(x) double x; { if(x>=.){ return (.); 83

85 2 }else{ return (.); } } double DISTANCE(i,j) int i; int j; { return (fabs(g[i]-g[j])+fabs(g2[i]-g2[j])+fabs(g3[i]-g3[j]) +fabs(g4[i]-g4[j])+fabs(g5[i]-g5[j])); } double COUNT(epsilon,n) double epsilon; int n; { int i,j; double sum; for(i=,sum=.; i<n; i++){ for(j=; j<n; j++){ sum = sum + STEP(epsilon-DISTANCE(i,j)); } } return ((double)sum/(n*n)); } main() { FILE *fpr; int i,k; double x,kx,k2x,k3x,k4x,kx,y,ky,k2y,k3y,k4y,ky,z,kz,k2z,k3z,k4z,kz; double arg,arg2,epsilon; if((fpr = fopen("m5.dat", "wt"))!=null){ for(i =,x=.,y=.,z=.; i < imax; i++){ kx = h*func(x,y,z); ky = h*func2(x,y,z); kz = h*func3(x,y,z); k2x = h*func(x+.5*kx,y+.5*ky,z+.5*kz); k2y = h*func2(x+.5*kx,y+.5*ky,z+.5*kz); k2z = h*func3(x+.5*kx,y+.5*ky,z+.5*kz); k3x = h*func(x+.5*k2x,y+.5*k2y,z+.5*k2z); k3y = h*func2(x+.5*k2x,y+.5*k2y,z+.5*k2z); k3z = h*func3(x+.5*k2x,y+.5*k2y,z+.5*k2z); k4x = h*func(x+k3x,y+k3y,z+k3z); k4y = h*func2(x+k3x,y+k3y,z+k3z); 84

86 2 k4z = h*func3(x+k3x,y+k3y,z+k3z); kx = (kx+2.*k2x+2.*k3x+k4x)/6.; ky = (ky+2.*k2y+2.*k3y+k4y)/6.; kz = (kz+2.*k2z+2.*k3z+k4z)/6.; x = x + kx; xx[i]=x; y = y + ky; z = z + kz; } for(k=; k< N; k++){ g[k] = xx[k]; g2[k] = xx[k+2]; g3[k] = xx[k+4]; g4[k] = xx[k+6]; g5[k] = xx[k+8]; } for(epsilon=.; epsilon < ; epsilon=epsilon+.){ arg = COUNT(epsilon,N); arg2 = epsilon; fprintf(fpr,"%lf %lf\n", log(arg2), log(arg)); } } fclose(fpr); } 2.4 q. D q = q lim ɛ log n(ɛ) i= pq i log ɛ (28), n(ɛ) ɛ ( ), N = n(ɛ)., q = 2 N log i= D 2 = lim p2 i ɛ log ɛ (29),., q = n(ɛ) log D = lim ɛ log ɛ i= = lim ɛ log n(ɛ) log ɛ (2),., q =. q /,, / q q 85

87 2 ( ), : n(ɛ) i= p i = { n(ɛ) lim q q log n(ɛ) } log p q i= i = lim pq i log p n(ɛ) i q n(ɛ) = p i log p i (2) i= i= pq i i= (q ) lim q q = (22) n(ɛ) lim D i= q = D = lim p i log p i q ɛ log ɛ (23). ( )., D q D 2, D, D.,,. 9 ( a =, b = 8/3, c = 3), m = 5, τ = 2. [ ]. (7/5) [ ]. 86

88 カオス フラクタル 講義ノート 2 年度版 3 担当: 大学院情報科学研究科 井上 純一 フラクタル [カオス編] は前回で終了したので, ここからはフラクタルについて学んでいく. 3. 自己相似性 自然界には縮尺 (スケール) を変えても次々と同じ構造 形が現れるようなものが多い. 例えば, 図 5 の上 2 枚の写真はいずれもブロッコリを撮影したものであるが, それら 2 つを並べた図 5(下) の 枚を見るまで, 同じ大きさのブロッコリに見える. この例からわかるように, ものによっては 図 5: 自己相似性を持つブロッコリ. 対象への視点をズームイン/アウトしても次々と同じ形状が現れ, 我々が日頃その典型的 (特徴的) な大きさを知っている 鉛筆 や タバコ などの第三の物体をそれらの傍らに置かないかぎり, その大きさについて知ることはできない (図では背景にあるノートの罫線が薄く見えるが, これに よってかろうじて 2 つの大小を比較することができる). このような性質を自己相似性と呼ぶ. この 性質は, ここから学習していく フラクタル と呼ばれる複雑図形の特徴の一つとなっている. このブロッコリの例以外にも雲や木, 海岸線など実に様々な対象がフラクタルの性質を持ってい る. また, このようなフラクタル図形の持つ性質を有効に使ったアルゴリズムを構築したり, ある いはコンピュータ グラフィクス (CG) に用いたりなど工学分野への 応用 も数多くなされてい る. ここからは, そのようなフラクタルについて基本的な事柄を学習していくが, [カオス編] 同様, この講義ではフラクタルの応用面/より進んだ内容について紹介/学習する余裕はないので, それら について興味のある者はこの講義を受講後に大学院講義 混沌系工学特論 などを別途聴講すると よいであろう. ここは 87 ページ目

89 2 3.2, [ ], r [, ], r, 2. r <. x n+ =.5 x n (24) y n+ =.6 y n (25). r <.2 x n+ =.5 x n (26) y n+ =.5 y n +. (27).2 r <.4 x n+ =.46 x n.32 y n (28) y n+ =.39 x n +.38 y n +.6 (29).4 r <.6 x n+ =.47 x n.5 y n (22) y n+ =.7 x n +.42 y n +. (22).6 r <.8 x n+ =.43 x n +.28 y n (222) y n+ =.25 x n +.45 y n +. (223).8 r x n+ =.42 x n +.26 y n (224) y n+ =.35 x n +.3 y n +.7 (225)., / (24)(25) (26)(27), /5 (28)(29), (22)(22), (222)(223), (224)(225). 2,, 5,, gnuplot ,., r,. 36 Step=, gnuplot with dots. 88

90 カオス フラクタル 講義ノート 2 年度版 担当: 大学院情報科学研究科 井上 純一 Step = Step = Step = Step = 図 52: 計算機を用いて描かれた木. ブロッコリに見えなくもない. #define max #define SEED 23 /* 乱数の種 */ main() { int i; double x,y,xn,yn,r; double ran3(long*); long idum=(seed); FILE *pt; xn =.5; /* 初期値 */ yn =.; /* 初期値 */ if((pt=fopen("tree.dat","wt"))!= NULL){ for(i=; i <= max; i++){ r = ran3(&idum); /* 乱数を入れる */ if(r <.){ xn =.5*x; yn =.6*y; }else if((r >.) && (r <.2)){ xn =.5*x; ここは 89 ページ目

91 2 yn = -.5*y +.; }else if((r >.2) && (r <.4)){ xn =.46*x -.32*y; yn =.39*x +.38*y +.6; }else if((r >.4) && (r <.6)){ xn =.47*x -.5*y; yn =.7*x +.42*y +.; }else if((r >.6) && (r <.8)){ xn =.43*x +.28*y; yn = -.25*x +.45*y +.; }else{ xn =.42*x +.26*y; yn = -.35*x +.3*y +.7; } fprintf(pt, "%f %f\n", xn, yn); x = xn; y = yn; } } fclose(pt); } 3.2.2, 53.,, :. 9

92 2, ( ) , :. S, S, S 2. L, L, L

93 2 S, S, S 2. L, L, L 2, l, L = 3l., L = 3l + (l/2) 3 = ( + /2)3l = (3/2)L., 2 L 2 = (3/2)L = (3/2) 2 L, n L n = (3/2) n L., ( ) n L n 3 = (226) L 2,.,, S = ( 3/4)l 2, S = ( 3/4)(l/2) 3 3 = ( 3/4)l 2 (3/4) = (3/4)S, n S n = (3/4) n S., ( ) n S n 3 = (227) S 4.,,. ( ), /, ( ).. t = t = t = 55:.. 55 ( ),. 3.3,. 92

94 2 3.3., 3,., /3 3,..,. 56: ,.,,,., /3,, /3, 2/3, /9, /3 n., 3. [ ], 2, 3. x, 3 x x = a 3 + a a a n 3 n + (228) (.a a 2 a n...) 3, a i {,, 2} 3., /3 /3 = n + = { 3 + } = 2 3 = /3 (229) 3, /3 = (.2 ) 3, /9 /3 /9 = n + = { 3 + } = 2 9 = /9 (23) 3 93

95 2, /9 = (.2 ) 3.,,, x = /3 + /9 = 4/9 3 4/9 = (.2 ) 3.,,, 3,, 3., [, ] 3,.,, 3., x = a 3 + a a a n 3 n +, a i (23), /3, x /3 x 3 = a a a a n 3 n + (232), x 3, x 3 i a i x/3,, x, x/3.,.. x n+ = { 3x n (x.5) 3x n + 3 (x >.5) (233), x,,,. 94

96 2 (7/5). { 3x n (x.5) x n+ = 3x n + 3 (x >.5) (234)., x =, x x = x =., x = /3 x =, 2 x 2 = 3 3 =, x =., x = /2, x,., x, 2., 57( )., x = /3, /9, x., x = /2, /5, x., x x x = /3 x = /9 x = /2 x = / n x 57:. (234). ( ).., x 57( )..,.,.,.. 95

97 2 4 [ ]. 4. [ ]. x n+ = ax n ( x n ) (235) ( )(235) a x x x 2 x n., a : x, z. z x, y, z z = x + iy, : z n+ = az n ( z n ) x n+ + iy n+ = a(x n + iy n )( x n iy n ) = a(x n x 2 n + yn) 2 + iay n ( 2x n ) (236). a a = a R + ia I,., (236), x n+ = a(x n x 2 n + yn) 2 (237) y n+ = ay n ( 2x n ) (238). (237)(238) x (x, y ), x-y,., 96

98 カオス フラクタル 講義ノート 2 年度版 式は収束せず, x 担当: 大学院情報科学研究科 井上 純一 2 のように発散してしまう. 従って, 全ての可能な初期条件 x2 + y はその初期条件から出発した軌道が [発散してしまうもの] と [有限に留まるもの] との 2 種類に分 類されることがわかる. つまり, 初期条件の作る集合 : {x } は {x } {x x } {x }finite {x x < } でそれぞれが定義される部分集合 {x } と {x }finite に分類される. 連立漸化式 (237)(238) を数値的に解き, 複素平面 x-y 内にこの {x }finite をプロットしてできる 図は, その作り方が上述のように極めて単純なのにも関わらず非常に複雑な形状を持つ. それを図 59 に示そう. この集合 {x }finite のことをジュリア集合と呼ぶ. 図 59 では a = 3.3 に選んだが, こ.25.2 a = 図 59: ジュリア集合. a = 3.3 の場合に得られるジュリア集合は自己相似性を有する. もちろん, a の値をいろいろ変ええれば様々 なジュリア集合が得られる (各自がいろいろ試してみると楽しいと思う). 図 6 にその中から数例 を載せる a = a = a = 図 6: 左から a = 3., 3.5, 及び, a = 3.6 に対するジュリア集合. 4.2 マンデルブロ集合 前節ではロジスティック写像をコントロールするパラメータ a を実数として扱ったが, 先に述べ たようにこれをも複素数に拡張することができる. a = ar + iai として (236) 式に代入し, 実部と ここは 97 ページ目

99 カオス フラクタル 講義ノート 2 年度版 担当: 大学院情報科学研究科 井上 純一 虚部にわけた反復式を書き出してみると xn+ = ar (xn x2n + yn2 ) ai yn ( 2xn ) (239) yn+ = ar yn ( 2xn ) + ai (xn (24) x2n + yn2 ) が得られる. ジュリア集合を求めたときには, a の値を固定し, 反復式 (239)(24) が有限の値に収 まるような初期条件の集合を複素平面内にプロットした. ここでは逆に, 初期条件 {x, y } が一つ 与えられた場合, ar, ai の値を様々変えたときに, 反復式 (239)(24) の解が有限の値を持つような 集合 {ar, ai } を x-y 平面にプロットしたならばどのような図形が得られるのかを考えてみたい. 図 図 6: マンデルブロ集合. x =.5, y = と初期条件を選んである. 6 に x =.5, y = を初期条件に選んだ場合の結果を載せる. このような図形のことをマンデル ブロ集合と呼ぶ. 確率的フラクタル 5 ここからは確率モデルを用いたフラクタルについて見ていく. 5. シェルピンスキー ガスケット 図 62 は前回学んだシェルピンスキー ガスケットであり, これは典型的なフラクタル図形の一 つである38. この図形は次に示すような非常に簡単なルールを用いて生成させることができる. () 図 63 のように頂点を (X, Y ) = (, ) (点 O), (X2, Y2 ) = (2l, ) (点 A), (X3, Y3 ) = (l, 3l) (点 B) とする正三角形を考える. (2) この正三角形内部の任意の点 P() = (x, y ) を選ぶ. (3) 数,,2 をランダムに選び O と P() を結ぶ線分の中点に点 P() = (x, y ) を置く. 38 ガスケット (gasket) とは 詰め物 の意味である. ここは 98 ページ目

100 カオス フラクタル 講義ノート 2 年度版 担当: 大学院情報科学研究科 井上 純一 図 62: シェルピンスキー ガスケット. B P O () P () A 図 63: O が選ばれた場合には P() と O の中点に P() をおく. ここは 99 ページ目

101 2 A P () P () = (x, y ). 2 B P () P () = (x, y ). (4) (3), P () (x, y ) P () (x, y ) P (n) (x n, y n ) 2. x y x y x 2 y 2 x 3 y 3... ( gasket.dat). () (4),. ( ) : gnuplot. 5.2,.,,,,.,, : P () (x, y ) P () (x, y ) P (n) (x n, y n ) 39,.,.,,,,. 64 ( ),., : 2 2.,. 39..

102 2 t = t = t = 64:.. /4 /4 P /4 /4 65: P 4. ( (, )), 2 (, ), 4 ( 65 ), () ( ) = l max l max. (l max ). S(l max /2, l max /2), r max r max = ( 67 ). (2) S r s (= r max + 2.) P(rx, ry). rx = r s cosφ

103 t= : 2. ry = r s sinφ, φ φ [, 2π] ( 68 ). (3) P, 2 2. (4) r (rx) 2 + (ry) 2,. (k) r r k (.), P. (2). (c), r k > r d > r s r d (= r max + 5.), r r d. 69, r, (r r s ), O 2, rx = rx + (r r s )cosθ ry = ry + (r r s )sinθ θ θ [, 2π]. (a), (rx + l max /2, ry + l max /2) 4 (rx + l max /2 +, ry + l max /2), (rx + l max /2, ry + l max /2 + ), (rx + l max /2, ry + l max /2), (rx + l max /2, ry + l max /2 ) (rx + l max /2, ry + l max /2), r max r max = max [r max, ] (rx) 2 + (ry) 2 2

104 2 l max S l max 67:.. P r s S 68: r s.., r max (2). (j) (k),(c),(a) (3). 3. [ ] main., check(), int(), occupy(), jump(), aggregate(), longjump();. (, ) /***********************************************************/ 3

105 2 O 2 r s S r 69: (r r s ) O 2. /* (main ) */ /***********************************************************/ #include<stdlib.h> #include<stdio.h> #include<math.h> #define lmax 62 #define rs (rmax+2.) #define rd (rmax+5.) #define rkill (5*rmax) #define PI /**********************************************/ /* */ /**********************************************/ #define SEED 3 #define SEED2 233 #define SEED3 5 /* (rx,ry) xf[rx][ry]=, xf[rx][ry]=. */ int xf[lmax][lmax]; int rx, ry; double rmax; main() { FILE *ptr; long i,k; rmax=; 4

106 2 init(); /* */ occupy(); /* P */ jump(); /* 4 */ for(i = ; i <= ; i++){ switch(check()) { /* check() (4) k,a,j,c */ case k : occupy(); jump(); break; case a : aggregate(); /* r_max */ if((ptr=fopen("test.dat","at"))!= NULL){ fprintf(ptr, "%d %d \n", rx, ry); /* (rx,ry) */ } fclose(ptr); occupy(); jump(); break; case j : jump(); break; case c : longjump(); /* (r-r_s). (c) */ break; } } } 5.3,.,,?,,, 2, 3.,. 5

107 2 l m(l)., l ( 7 )., l,, m() = 4, m(2) = 9,, m(l) = (l + ) 2, l m(l) = A l 2., A. m(l) = A l 3. 7: ( ).,,, l m(l), l m m(l) = A l d f, d f.,,,, 2, 3,....,., 2,,. d f ( )...,. 6

108 2 4 L l l N, l L = l max ( ), l N(l). log N log l, d f ( ). 9 log(n) log(l) 7

109 2 6 4, [ ] : x n+ = ax 2 n + y n (24) y n+ = bx n (242). () a =.6, b =.2, x n, y n,., (a, b),. (2) x-y, ( 7 )., gnuplot,. (3)., τ, m (τ, m), a =.6, b = y x 7:. 8

110 2 6 : x n+ = 4x n ( x n ) (243),. x = T T x n (244) n=, T,.,,,. P (x) = π x( x) (245),. x = xp (x) = xdx (246) π x( x),, x = x. 6 Remark: : x. x = sin 2 (θ/2), dx = cos(θ/2) sin(θ/2)dθ, x = π π sin 2 (θ/2) sin(θ/2) cos(θ/2) dθ = sin(θ/2) cos(θ/2) 2π π ( cos θ)dθ = 2 (247)., x /2.,, x = x, [, ] / x( x) : I = f(x) x( x) (248),., f(x)., x, x,, x T, T π I = lim f(x i ) (249) T T. (, T ). i= 9

111 2,., I ( ) x =,, ( )., I ( ). ( ). π,,,., /4,., /4 = /4 n (= N) π /4 = 4 = π 4, N, /4., ( ) C.

112 2 () (x, y), x, y [, ], 2. x n+ = 4x n ( x n ), x =. y n+ = 4y n ( y n ), y =. (x, y ), (x, y ),, (x k, y k ),., N N =. (, ) (x k, y k ), /4, n. n/n π/4. pi/4=.... (2) x n. x =. x n for for(i=,i<=n; i++){ x = Logmap(i,x); } Logmap(n,x).. /************************* () **********************************/ #include<stdio.h> // #include<stdlib.h> /* rand() */ #include<math.h> #define N main() { int i,sum; double x,y,p; x=.; /* */ y=.; for(i=,sum=; i<=n; i++){ x=4.*x*(.-x); y=4.*y*(.-y); // x=(double)rand()/rand_max; /* */ // y=(double)rand()/rand_max; if(sqrt((x*x)+(y*y))<.){ sum++;

113 2 }else{ sum= sum; } } p=(double)sum/n; printf("pi/4=%lf\n",p); } /************************* (2) **********************************/ double Logmap(int n, double x) { if(n==){ return (x); }else{ return (4.*Logmap(n-,x)*(.-Logmap(n-,x))); } } (), (x n, y n ) 2. 72, C rand(). N = 7, N = 3 7., /4,, (x n, y n ) = (, ), (, ), (, ) 4., /4,.. π,, π,,, N = 3 7, π 2.882, 3.44.,,., / x( x)., 6., x n y n,.,, y = x,,.. : 8/2( ) 5 ( ), 8-3 4, /4, (x n, y n) = (, ). 2

114 カオス フラクタル 講義ノート 2 年度版 担当: 大学院情報科学研究科 井上 純一 Logistic map: N= sqrt(-x*x) Random: N= sqrt(-x*x) Logistic map: N=3 sqrt(-x*x).6.8 Random: N=3 sqrt(-x*x) 図 72: 上段は問題に与えられた漸化式を用いた場合の分布図であり, 下段は C 言語で定義されている rand() という疑似 乱数生成関数を用いたときの分布図である. 右側が N = 7, 左側が N = 3 7 と選んである. まで. 全てのコードおよび, 結果のグラフ/図などを全て提出. また, 各課題に関し, 独自の考察な どを含めても良い. [カオス フラクタル] 今後数回の予定: 7 月 9 日: 通常講義. マルチ フラクタル について 7 月 26 日: 最終回. 本講義に関連する研究紹介 ( 時間程度を予定). なお, 後期開講の 情報工学演習 II では, 本講義では扱うことのできなかった, 臨界現象に関す るフラクタル図形とその特徴付けについて学習することになります. 本講義はこの演習とのセット で, カオスとフラクタルについての基本的な知識が得られるように設計しました. ここは 3 ページ目

115 q.,. ɛ p i, D q. D q = q lim ɛ log n(ɛ) i= pq i log ɛ (25).,, ɛ., i p i (ɛ)., q = D q= = lim n(ɛ) log i= ɛ log(ɛ) = lim ɛ log n(ɛ) log(ɛ) (25), ɛ ( 73 ), ɛ, n(ɛ) n(ɛ) ɛ D q= (252), D q= (, )., q =, q = 2. 73: ɛ. 4

116 2 7.2 q (25), q > i, (25)., q,., q,, q. q =, q, ( q ) 4. q (25),.,,, q,,., , [, ].,, 2, p, p., [,.5] 74:. n = 2. 4, ( ). 42,. 5

117 カオス フラクタル 講義ノート 2 年度版 担当: 大学院情報科学研究科 井上 純一 さらに 2 分割し, やはり, 左側に p, 右側に p を割り当て, 回目の分割でできた [, 5, ] の区間 についても同様に右側に p を, 左側に p を割り当てる. このような割り振りでできる 2 分木 の末端にはそれぞれの経路 ( 枝 ) に割り振られた確率の積が与えられる. 従って, 2 回目の分割 では区間 [,.25], [.25,.5], [.5,.75], [.75, ] の 4 区間にそれぞれ p2, p( p), ( p)p, ( p)2 が割り振られることになる (図 74 参照). このプロセスを n 回繰り返すと線分 [, ] は 2n 等分され, その中で確率 pk ( p)n k が割り振 られている小区間は n Ck 個となる. このようにしてできる小区間を横軸, 各小区間に割り当てられ た確率を縦軸にプロットすると図 75 のようになる. この図中の組み込み図には区間 [,.] を拡大 p = e-5.6 6e-5.4 4e-5.2 2e p = 図 75: 二項分岐過程の密度分布. 上図では p =.25 と選んである. 上図の組み込み図は同図のスケールを変えてプロッ トしたもの. 下図はそれぞれ p =.,.4 の場合の二項分岐過程の密度分布. n = 8 と選んでいる. した部分が載せてある. この図からわかるように, この図形は 自己相似性 を持つ. 従って, 具体的に計算機で実行するためのコーディングの際には で 2 分岐の左側を で 右側を表すことに約束すると, n ステップ後の 2 分木の末端 (いわゆる 葉 ) は ( )2, ( )2,, ( )2 等と 2 進数でコードされる. よって, n 個の, の並びのうち k 個の を含む葉には確率 pk ( p)n k を割り当てることにし, この葉の位置する [, ] 区間の値は この葉の 2 進数 (xn xn 2 x )2, (xi = の個数が k 個) を 進数表示したものを 2n でスケー ここは 6 ページ目

118 2., [, ] X. X = n i= x i2 i 2 n (253) (x n x n 2 x ) 2 = ( ) 2 } {{ } k, n k X = n 2 n x i 2 i (254) i= Y = p k ( p) n k (255). (X, Y ) , D q : Z q (ɛ) = n(ɛ) [p i (ɛ)] q (256) i= D q = q lim log Z q (ɛ) ɛ log ɛ (257) 44., (, ). Z q (ɛ) = p n(ɛ) δ(p i (ɛ) p)[p i (ɛ)] q = p i= n(ɛ) δ(p i (ɛ) p)p q (258) i=, n(ɛ) i= δ(p i(ɛ) p), n(ɛ) p, ρ(p), (258). Z q = p ρ(p)p q (259), ρ(p) p., p p k ( p) n k, ρ(p) n C k Z q = = = n nc k [p k ( p) n k ] q k= n nc k p qk ( p) q(n k) k= n nc k (p q ) k {( p) q } n k = [p q + ( p) q ] n (26) k= 43 Y X n C k., n. n = 8,,. 44. q T β = T, i i E i p i = e E i, Z β = i [p i] β = i e βe i,. 7

119 2., ɛ = 2 n, D q = log[p q + ( p) q ] n q log 2 n = log[p q + ( p) q ] q log(/2) (26). p D q p =. p =.25 p = D q.5 D q.5 p =.25 D min D max q q 76: D q. p p =.,.25,.4. p =.25 D q q ±., : D max = log p/ log(/2), D min = log( p)/ log(/2)., q 2 +2, D q,, [D q=+, D q= ]. (26)., p < /2 q +, (26) 2 ( p) q D q=+ q log( p) log( p) = q log(/2) log(/2), q, (26) 2 p q D q=. 76( ). q log p q log(/2) = log p log(/2) (262) (263) q =,, q = +,.,, ,., q, ,,,. 8

120 2 7.4., ɛ i ɛ p i (ɛ) ɛ αi (264) : ɛ., α i,. [α, α + dα] α ρ(α) ɛ f(α) dα (265)., ρ(α), ɛ f(α)., α, f(α). α, Z q = ρ(α) ɛ f(α) [ɛ α ] q dα = ρ(α) ɛ f(α)+qα dα (266)., (266) (256) ρ(α) ɛ f(α) = n(ɛ) δ(α α i ) (267) i=. Z q = n(ɛ) δ(α α i ) ɛ qα dα = ɛ qαi = [p i (ɛ)] q (268) i=. (256).,, Z q α, ɛ ɛ f(α)+qα ( ). α α(q), { f(α) + qα} = (269) α α=α(q) n(ɛ) i= n(ɛ) i= q = f(α) α (27) α=α(q)., f(α) + qα f(α) + qα f(α(q)) + qα(q) + 2! 2 f(α) α 2 (α α(q)) 2 (27) α=α(q) (, 2 f(α)/ α 2 α=α(q) < ), Z q = ρ(α(q))ɛ f(α(q))+qα(q) ɛ 2 f(α) 2 α 2 α=α(q)(α α(q)) 2 dα (272) 9

121 2. ɛ 2 f(α) 2 α 2 α=α(q)(α α(q)) 2 dα = = exp 2 log f(α) ɛ 2 α 2 (α α(q)) 2 dα α=α(q) π (273) 2(log ɛ) 2 f(α) α 2 α=α(q) 46, D q (q )D q = lim ɛ log Z q log ɛ = log ρ(α(q)) 2(log ɛ) 2 f(α) α lim f(α(q)) + qα(q) + lim α=α(q) ɛ log ɛ ɛ log ɛ = f(α(q)) + qα(q) (274). lim ɛ (log ɛ) =, lim ɛ (log ɛ) 3/2 =.,. τ(q) (q )D q (275) τ(q) = f(α(q)) + qα(q) (276) f(α(q)) = qα(q) τ(q) (277) π dτ dq = τ = α + ( τ q + α { q f α ) ( ) α q } α q = α (278)., f(α(q)) = qα(q) τ(q) (279) dτ dq = α (28) 47., 2 f(α)/ α 2 α=α(q) <, f(α)., α (28) τ q (279) f(q(α)) = qα(q) τ(q) α() τ() = τ() = D q= = f() (28) 46 ɛ bα2 dα = e log ɛ bα2 dα = e b(log ɛ)α2 dα = 2 π b log ɛ. 47, (279)(28) f τ., df/dα = y P (y) = αy f = α df dα f = αq f = τ, τ f. 2

122 2., f q = D q=., f(α), (27) q + α = α min, q α = α max., (28), τ(q) αq τ(q + ) α min (+ ) (282) τ(q ) α max ( ) (283), D q = τ(q)/(q ) D q = {α max ( )}/( ) = α max, D q + = {α min (+ )}/(+ ) = α min, α min, α max f. f(α min ) = f(α(q + )) = (+ ) α(q + ) τ(q + ) = (+ ) α min α min (+ ) = (284) f(α max ) = f(α(q )) = ( ) α(q ) τ(q ), (27)(279) = ( ) α max α max ( ) = (285) f(q(α)) = qα(q) τ(q) = f(α) α α=α(q) α(q) τ(q) (286),, : q = f(α) α α=α(q) = (287) : f = α (279) α, ( q)α = τ(q) = ( q)d q α = D q= (288), f f(α()) = α = D q= (289) ,.,. (26), D q =, τ(q) p q + ( p) q q log(/2) (29) τ(q) = (q )D q = pq + ( p) q log(/2) (29) 2

123 2., (28) α = dτ dq = log(/2) [ p q log p + ( p) q ] log( p) p q + ( p) q (292). (292) ξ = ξ = p q p q + ( p) q (293) ( p) q p q + ( p) q (294) α = (ξ log p + ( ξ) log( p)) (295) log(/2),, (279) f = = q log(/2) (ξ log p + ( ξ) log( p)) log(/2) log(pq + ( p) q ) (ξ log ξ + ( ξ) log( p)) (296) log(/2)., α α = f = (ξ log p + ( ξ) log( p)) log(/2) (297) (ξ log ξ + ( ξ) log( ξ)) log(/2) (298)., ξ, p 77( )., D q=, D q=,, α. (293)(294) ξ ξ = ( ) q p (299) p., q =, ξ = /2,, f (298) f = D q= =, α (297) α = log p( p) 2 log(/2) (3)., q =, (299) ξ = p,,, (298)(297) f = α = p log p + ( p) log( p) log(/2) = D q= (3) 22

124 2.8 p =. p =.25 p = p =.25 f = α D q = D q =.6 f f α α.9 p =. p =.25 p = p =. p =.25 p = f α q q 77: f α. p p =.,.25,.4 ( ). p =.25, f D q=, D q=. f q ( ),, α q ( ).. p =.25 77( ).., q. 23

125 2, A4.,.,,., #. : ( ) 5 ( ), 8-3. : 7/26 ( ): 78: ( ). 24