ランダムウォークの境界条件・偏微分方程式の数値計算

B L06(2018-05-22 Tue) : Time-stamp: 2018-05-22 Tue 21:53 JST hig,, 2, multiply transf http://hig3.net L06 B(2018) 1 / 38

L05-Q1 Quiz : 1 M λ 1 = 1 u 1 ( ). M u 1 = u 1, u 1 = ( 3 4 ) s (s 0)., u 1 = 1 7 ( 3 4 ). 2 p(1) = M p(0) = 1 3 ( 1 2 ). 3 M ( ) λ 1 λ 2, u 1, u 2, λ 1 = 1, λ 2 = 1 6, u 1 = ( 3 4 ) s, u 2 = ( 1 1 ) s (s 0). λ 1 = 1,. u 1, u 2 s = 1 ( p(t) ). L06 B(2018) 2 / 38

p(t) =M t p(0) =(UDU 1 ) t p(0) = ( ) ( ) ( ) λ t 1 u 1 u 1 0 2 0 λ t u 1 u 2 p(0) 2 U 1 p(0), ( a b ). ( ) p(t) = u 1 λ t 1 u 2λ t 2 ( a b ) =a u 1 λ t 1 + b u 2 λ t 2. p(0), p(0) = a u 1 + b u 2 a, b a = 1 7, b = 4 7., p(t) = 1 7 ( 3 4 ) + 4 ( 1 ) 7 1 ( 1 6 )t. L06 B(2018) 3 / 38

t + p(t), u 1 1 7 (= ).. 1 6 < 1, p(t) 1 7 ( 3 4 ) (t + ) p(x,t) 1 0.75 4/7 0.5 3/7 0.25 0 0 1 2 3 4 t p(0,t) p(1,t) 4 p(t) = 1 7 ( 3 4 ) + 1 ( 1 ) 14 1 ( 1 6 )t L06 B(2018) 4 / 38

p( ) = 1 7 ( 3 4 ) ( ), p(t) p( ) = ( 4 1 ) 7 1 1 6 t = 4 2 7 1 6 t log p(t) p( ) =t log 1 6 + log 4 2 7, t log, log 2. log p(x,t)-p( ) 0-1 -2-3 -4-5 -6-7 -8 log p(t)-p( ) 0 1 2 3 4 t L06 B(2018) 5 / 38

3, ( 2 ) 3, t. L06-Q2 Quiz : ( 1 λ = 1 1 11 ) 3 2, p(t) = 1 3 ( 11 1 ) 1 6 L06-Q3 Quiz : ( 1 ) ( 0 ( 3 1 4 )t + 1 1 ) 2 6 ( 1 1 4 )t 1 L06 B(2018) 6 / 38

1, E[(X(t) + 1) 2 ] = p(t) = 1 7 ( 3 4 ) + 4 7 1 (x + 1) 2 p(x, t) x=0 ( 1 1 ) ( 1 6 )t. =(0 + 1) 2 ( 3 7 + 4 7 ( 1 6 )t ) + (1 + 1) 2 ( 4 7 + 4 7 ( 1 6 )t ) = 19 7 12 7 ( 1 6 )t 2,, t + u 1 = 1 7 ( 3 4 ) (0 + 1)2 3 7 + (1 + 1)2 4 7. P (X(t) > 0) = E[1 [X>0] (X)] =(0 1) p(t) =p(1, t) = 4 7 4 7 ( 1 6 )t L06 B(2018) 7 / 38

L05-Q4 Quiz : ( λ = 1 ( ), 10 ) ( 01 ) s + k (s 0 k 0). 1, 0 1 ( 1 ( ) 10 ) (, u 1 =, u 2 = 1 01 ) 2, s u 1 + (1 s) u 2 (0 s 1) 0 1. ( λ = 1 3, 0 ) +1 s(s 0)., ( ( 0 ) ), u 3 =. p(0) +1 1 1 p(t) = a u 1 1 t + b u 2 1 t + c u 3 ( 1 3 )t. L06 B(2018) 8 / 38

( 1 p(0) = 1 11 ) 2 a, b, c, 0 p(t) = 1 2 u 1 + 1 2 u 2 + 1 4 u 3 ( 1 3 )t ( p( ) = 1 21 ) 4. ( 1 2 p(0) = 1 11 ) 3 a, b, c, 1 p(t) = 1 3 u 1 + 2 3 u 2 + 0 u 3 ( 1 3 )t ( p( ) = 1 11 ) 3. 1 3 {0} {1, 2}.,. L06 B(2018) 9 / 38

5 +1 : 6 p(x, t) L06 B(2018) 10 / 38

C I x = 0,..., m 1 m. p(t), p(x, t) 1 double p [ m ] = { 1. 0, 0. 0,...., 0. 0 } ; /. m / 2 / {p ( 0, t ), p ( 1, t ), p ( 2, t ),..., p (m 1, t )} / M = ( p 00 p 01 p 10 p 11 ) = ( 0.1 0.3 0.9 0.7 ) 1 double M[ ] [ 2 ] = { { 0. 1, 0. 3 }, 2 { 0. 9, 0. 7 } } ; / 2 / {{p 00, p 01 }, {p 10, p 11 }} q = M p q x = y M xy p y. 1 p [ ] p ( x, 0 ) ; 2 p ; 3 f o r ( t ){ 4 pn=m p ; / / 5 p=pn ; 6 p ; 7 } L06 B(2018) 11 / 38

+1 : 5 +1 : 6 p(x, t) L06 B(2018) 12 / 38

+1 : X(2) :, ( ) t = 2 x = 3 t = 20 x < 0 X(T ) : T ( ) x(t ) : T ( 1 ) (x(0), x(1), x(2),..., x(t )) : (path) ( ) L06 B(2018) 13 / 38

+1 : ϕ(x(t)) C B(2018)L02 ϕ(x(t )) = 1 N N ϕ(x(t ) (n) ) n=1 1 / 1 / 2 f o r ( n ){ 3 / 2 / 4 f o r ( t ){ 5 / 3 / 6 x=x+getrandom ( g e t u n i f o r m ( ) ) ; 7 / 4 / 8 } 9 / 5 / 10 } 11 / 6 / sum1=0, sum1+=phi(x) ( ϕ(x) = x 3 ), printf( %f,(double)sum1/n)? L06 B(2018) 14 / 38

+1 : C ϕ(x) = 1 [ ] (x) phi? sum1. L06 B(2018) 15 / 38

+1 : ϕ(x(0), X(1),..., X(T ))! : X(0), X(1),..., X(T ) 10. 1 i n t path [TMAX] ; 2 / 1 / 3 f o r ( n ){ 4 / 2 / 5 f o r ( t ){ 6 / 3 / 7 x=x+getrandom ( g e t u n i f o r m ( ) ) ; 8 path [ t ]=x ; / 4 / 9 } 10 / 5 / 11 } 12 / 6 / 13 14 i n t p h i ( i n t path [ ], i n t tend ){ 15 / path [ 0 ],.., path [ tend ] / 16 r e t u r n 0 ; / or 1 / 17 } sum1=0, sum1+=phi(path,t), printf( %f,(double)sum1/n)? phi 1,0. sum count. L06 B(2018) 16 / 38

5 +1 : 6 p(x, t) L06 B(2018) 17 / 38

2 1 sim* X(t) = X(t 1) + R(t), P (R(t) = ±1) = 1 2. P (X(3) = 9) = 1. P (X(t) = x) = p(x, t) 2 markov { p(x, t) = 1 2 p(x 1, t 1) + 1 2 p(x + 1, t 1). p(x, 3) = 1(x = 9) 0( ) p(t) = M p(t 1),, < X(t) < +.? 1 2 int x 2 p(x, t) = double p[m], 0 x < M. M,?. 2 p, M? L06 B(2018) 18 / 38

S = {0, 1,..., m 1}. x = 0, x = m 1, ( )?. h = 1/2, m = 6. p(0, t) p(1, t). p(m 2, t) p(m 1, t) p(t) = M p(t 1) 0 h 0 0 0 0 h 0 h 0 0 0 = 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0! p(0, t ) p(1, t ). p(m 2, t ) p(m 1, t ) L06 B(2018) 19 / 38

= 1 x m 2, x = 0, x = m 1 ( ). x = 0 =.. x = 1 x = 0 x = 1 x = 0 x = 1 x = 1 x = 0 x = m 2 ( ). x = 0 x = m 2. x = m 2 x = m 1 X M.,, x = 0, m 1. L06 B(2018) 20 / 38

p M? 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 ( A),p(0, t) =, ( B) 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 ( A), p (0, t) =, ( B) x 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0 h 0 0 0 0 h 0, p(0, t) = p(m 1, t) L06 B(2018) 21 / 38

L06-Q1 Quiz( ) {x} = {0, 1, 2, 3, 4} X(t),. X(t) =X(t 1) + R(t), (t = 1, 2,...) X(0) =2, R(t) (t = 0, 1, 2,...) 1 7 (r = 1) 4 P (R(t) = r) = 7 (r = 0) 2 7 (r = +1) 0 ( )., x = 0, 4.. 1 M. 2. L06 B(2018) 22 / 38

( ),. ( ) ( S = {x},,, ( ) ) ( ) AR(M), ARMA(m), ARIMA(m), ( ) L06 B(2018) 23 / 38

p(x, t) 5 +1 : 6 p(x, t) L06 B(2018) 24 / 38

p(x, t) p(x, t) x, t:, p(x, t),x(t):, t + 1, x ± 1,,, x, t:, p(x, t) X(t):, t + t, x ± x, p(x, t) = 1 2 p(x 1, t 1) + 1 2p(x + 1, t 1) p(x, t + 1) = 1 2 p(x 1, t) + 1 2p(x + 1, t) p(x, t + t) = 1 2 p(x x, t) + 1 2p(x + x, t) : f(x + x) f(x) f (x) x f(x + x) f(x) x f(x) f(x x) x df(x) (x) ( x 0) dx df(x) (x) ( x 0) dx L06 B(2018) 25 / 38

p(x, t) ± t, ± x p(x, t + t) p(x, t) = 1 [(p(x + x, t) p(x, t)) 2 (p(x, t) p(x x, t))] p t (x, t) t =1 2 p t (x, t) t =1 2 ( ) p p (x, t) (x x, t) x x x ( ) p (x, t) x x x x p t (x, t) =1 ( x) 2 2 p (x, t) 2 t x2 p(x, t), u(x, t). 1 ( x) 2 2 t D > 0:. A, p x (x, t). ( x x). L06 B(2018) 26 / 38

p(x, t) (diffusion equation, heat equation) u t (x, t) =D 2 u x 2 (x, t) (x min < x < x max, t > 0) u(x, 0) =x (x min < x < x max ) u(x min, t) =u(x max, t) = 0 (t 0) x,,,, u(x, t) : t, x u:, x, t: L06 B(2018) 27 / 38

p(x, t) (PDE=partial differential equation), u(x 1, x 2,..., x n ), u (t) = 2u(t). x (t) = x(t). x(t): x. u(x, t): u(x) (4 ) http: //www.a.math.ryukoku.ac.jp/~hig/course/compsci2_2013/img/pde-diff.gif L06 B(2018) 28 / 38

p(x, t) 2,, A B L06 B(2018) 29 / 38

p(x, t) D > 0. A u(x, t) = a(2t + x 2 ) + bx + c.,,,. 1 u(x, t) = e x2 2Dt. t, 0, 2πDt Dt N(0, Dt). u(x, t) = e c2 Dt sin(cx). (c R ) L06 B(2018) 30 / 38

p(x, t) L06-Q2 Quiz( ),? 1 B p(x, t) 2 II mx = kx bx. 3 II I x + ax + bx = c. 4-5 A 6 II, L06 B(2018) 31 / 38

p(x, t) L06-Q3 Quiz( ) t, x. ( ). u t (x, t) =2 2 u (x, t) (0 < x < 2π) x2 u(x, 0) = sin(3x) (0 < x < 2π) u(0, t) =u(2π, t) = 0 (t 0) u(x, t) = Ae Bt sin(cx), A, B, C R,. L06 B(2018) 32 / 38

5 +1 : 6 p(x, t) L06 B(2018) 33 / 38

:,. < x < +, p(t) 100. L06-Q4 Quiz( ), S = {x} = {0, 1, 2,..., m 1}. 0 1 0 0. 0 0 1................. M =........ 1 0. 0.. 0 1 1 0 0 0 1 double p [m], q [m] ; p, q, p q = M p 1 i n t m u l t i p l y t r a n s ( double q [ ], double p [ ], i n t m) ;. M 2, M (= O(m 2 ) O(m) ). L06 B(2018) 34 / 38

m 1 i n t m u l t i p l y t r a n s ( double q [ ], double p [ ], i n t m){ 2 i n t x, y ; 3 double M[ ] [ NS ] = { { 0. 0, 1. 0, 0. 0, / /, 0, } ; 4 f o r ( y =0;y<m; y++){ 5 q [ y ]=0; 6 f o r ( x =0;x<m; x++){ 7 q [ y]+=m[ y ] [ x ] p [ x ] ; 8 } 9 } 10 r e t u r n 1 ; 11 } Quiz : 1: 1 i n t m u l t i p l y t r a n s ( double q [ ], double p [ ], i n t m){ 2 i n t x ; 3 f o r ( x =0;x<m 1;x++){ 4 q [ x ]=1.0 p [ x +1] / +0.0 p [ x+2] / ; 5 } 6 q [m 1]=1.0 p [ 0 ] 7 r e t u r n 0 ; 8 } q = M p. m 1 q x = M xy p y = 0 + + 0 + 1 p x+1 + 0 + + 0. y=0 sparse matrix 0. 2,. L06 B(2018) 35 / 38

L06-Q5 Quiz( ), {x} = {0, 1, 2,..., m 1}. 7 2 10 10 0 0 3 5 2 10 10 10 0. 3 5. 0..... M = 10 10... 0...... 2 10 0.... 3 5 2 10 10 10 3 8 0 0 10 10 1 double p [m], q [m] ; p, q, p q = M p 1 i n t m u l t i p l y t r a n s ( double q [ ], double p [ ], i n t m) ;. M 2, M. L06 B(2018) 36 / 38

( )! 2018-06-05 5, 10 ( )+80 ( ), 20..., Visual Studio Excel R. 2018-05-29. (0) (1), (2) p(x, t), (3), (4) n p(x, t), (L03) (L03) M p(t) (L05), (L04,L05), X, M (L06),,, (L06) C. srand rand (L01,L02) (sim11,sim13) (L07) L06 B(2018) 37 / 38

https://learn.math. ryukoku.ac.jp/moodle. : CamScanner on https://download. ios/android moodle.org/mobile https://www.camscanner. com/ /Math 1-614 2018-06-05 ( ) Visual Studio.. L06 B(2018) 38 / 38