H25 1 2 1 seto@ics.nara-wu.ac.jp
Euler xi+1 x ( xi f(ti, xi) ) t ti ti+1 h
Euler Euler x xi+1 k 1 = hf (t i, x i ) xi+1 x i+1 = x i + hf t i + h 2, x i + k 1 2 xi f t i + t i+1 2, x i + x i+1 2 t ti ti + h/2 ti+1 h
x(t) t = ti x(t i+1 )=x(t i + h)=x(t i )+hx (t i )+ h2 2 x (t i)+ Taylor O(h 2 ) ( 2 )
Runge-Kutta ( ) Euler 4 k 1 = f (t i, x i ) k 2 = f (t i + h 2, x i + h 2 k 1) k 3 = f (t i + h 2, x i + h 2 k 2) k 4 = f (t i + h, x i + hk 3 ) Taylor O(h 4 ) ( 4 ) x i+1 = x i + h k 1 6 + k 2 3 + k 3 3 + k 4 6
Runge-Kutta ( ) 4 x k 1 = f (t i, x i ) 1 k 2 = f (t i + h 2, x i + h 2 k 1) 2 xi+1 4 k 3 = f (t i + h 2, x i + h 2 k 2) 3 3 k 4 = f (t i + h, x i + hk 3 ) 4 xi 1 2 x i+1 = x i + h k 1 6 + k 2 3 + k 3 3 + k 4 6 ti ti + h/2 ti+1 t h
Runge-Kutta ( ) 4 x k 1 = f (t i, x i ) 1 xi+1 k 2 = f (t i + h 2, x i + h 2 k 1) 2 k 3 = f (t i + h 2, x i + h 2 k 2) 3 k 4 = f (t i + h, x i + hk 3 ) 4 xi 1 2 x i+1 = x i + h k 1 6 + k 2 3 + k 3 3 + k 4 6 ti ti + h/2 ti+1 t h
Runge-Kutta ( ) 4 x k 1 = f (t i, x i ) 1 xi+1 k 2 = f (t i + h 2, x i + h 2 k 1) 2 k 3 = f (t i + h 2, x i + h 2 k 2) 3 3 k 4 = f (t i + h, x i + hk 3 ) 4 xi 1 2 x i+1 = x i + h k 1 6 + k 2 3 + k 3 3 + k 4 6 ti ti + h/2 ti+1 t h
Runge-Kutta ( ) 4 x k 1 = f (t i, x i ) 1 k 2 = f (t i + h 2, x i + h 2 k 1) 2 xi+1 4 k 3 = f (t i + h 2, x i + h 2 k 2) 3 3 k 4 = f (t i + h, x i + hk 3 ) 4 xi 1 2 x i+1 = x i + h k 1 6 + k 2 3 + k 3 3 + k 4 6 ti ti + h/2 ti+1 t h
(, 1 ) #define DT 0.01 /* */ #define STEP_MAX 1000 /* DT*STEP_MAX = 10.0 */ double fn(double, double); /* */ void runge_kutta(double, double, double*, double); /* */ main(){ long step; double t, x, x_next; dx dt = f (t,x) } x=0.1; /* */ for(i=0; i<step_max; i++){ t = i*dt; euler(t, x, &x_next, DT); x = x_next; } t = 0 t = DT*STEP_MAX void runge_kutta(double t, double x, double *x_out, double h){... } double fn(double t, double x){... }
(, 1 ) runge_kutta ti x(ti), h ti+1 = ti + h x(ti+1) ti x(ti) x(ti+1) void runge_kutta(double t, double x, double *x_out, double h){ h k 1 = f (t i, x i ) k 2 = f (t i + h 2, x i + h 2 k 1) k 3 = f (t i + h 2, x i + h 2 k 2) k 4 = f (t i + h, x i + hk 3 ) k 1 = hf(t i, x i ) k 2 = hf(t i + h 2, x i + k 1 2 ) k 3 = hf(t i + h 2, x i + k 2 2 ) k 4 = hf(t i + h, x i + k 3 ) } x i+1 = x i + h k 1 6 + k 2 3 + k 3 3 + k 4 6 x i+1 = x i + k 1 6 + k 2 3 + k 3 3 + k 4 6
H22 1 1 3 (1) (2) 1. (x(0) = 1) h = 0.1 h = 0.05 ( ) h 2. Mathematica (1) dx dt = x 0 t 1 (2) dx dt = x sin (t) 0 t 50
4 (1)-(3) h = 0.05 Mathematica (1) dx dt = t 1 x 2, x(0) x(0) = 0.05 = 0 0 t 5 (2) t dx dt = x + t 2 + x 2, x(0) = 1 1 1 t 10 (3) dx dt = x 1 + t + cos t, x(0) = 1 0 t 100
(1) dx dt = t 1 x 2, x(0) = 0.05 1: dy dx = f (x)g(y) (2) t dx dt 1: = x + t 2 + x 2, x(0) = 1 2: dx dy dx = f y x x 2 + a 2 = sinh 1 x a + C (3) dx dt sinh -1 x ( ) Mathematica ArcSinh[t] = x 1 + t + cos t, x(0) = 1 1: 1 x + P(t)x = Q(t) 2:? µ = e P(t)dt 3: d (1 + t)x =? dt
1
1 k 1 = f(t, x) k 2 = f(t + h 2, x + h 2 k 1) k 3 = f(t + h 2, x + h 2 k 2) k 4 = f(t + h, x + hk 3 ) x(t + h) = x(t) + h k 1 6 + k 2 3 + k 3 3 + k 4 6
(, ) #define DIM 2 /* ( ) */ #define DT 0.01 /* */ #define STEP_MAX 1000 /* DT*STEP_MAX = 10.0 */ void derives(double, double[], double[], int); /* */ dx dt = f (t, x) void runge_kutta(double, double, double*, double); /* */ main(){ long step; double t, x, x_next; } x=0.1; /* */ for(i=0; i<step_max; i++){ t = i*dt; euler(t, x, &x_next, DT); x = x_next; } t = 0 t = DT*STEP_MAX void runge_kutta(double t, double x, double *x_out, double h){... } double fn(double t, double x){... }
(, ) derives t x(ti) derivatives[] ti x(ti) dx/dt or f(t, x) void derives(double t, double x[], double derivatives[]){ f(t, x) derivatives[] derivative[0] = a*(x[1]- x[0]); }
(, ) runge t x(ti), h x(ti+1) x_out[] void derives(double t, double x[], double x_out[], double h){ } k 1 = f(t, x) ti k 2 = f(t + h 2, x + h 2 k 1) k 3 = f(t + h 2, x + h 2 k 2) k 4 = f(t + h, x + hk 3 ) x(t + h) = x(t) + h k 1 6 + k 2 3 + k 3 3 + k 4 6 x(ti) x(ti+1) k 1 ~k 4 derives[] x_out[] h
-
Lotka-Volterra dx dt = ax dy dt = bxy cy + dxy a: b: c: d:
5 a = 1, b = 0.03, c = 1, d = 0.025 x(0), y(0)
1963 E. N. Lorenz dx 1 dt dx 2 dt dx 3 dt = a (x 2 x 1 ) = bx 1 x 2 x 1 x 3 = x 1 x 2 cx 3 a: b: c:
6 (1) a = 10, b = 28, c = 8.0/3 3 (a) (x1, x2, x3) = (0, 2, 0) (b) (x1, x2, x3) = (0, 1.99, 0) (2)(1) (a) (b) x1 x1? t
Mathematica (1/4) Mathematica, ti x1(ti) x2(ti) x3(ti) 1) 0.000 1.00 2.00 3.00 0.001 1.02 1.97 3.03 0.002 1.05 1.94 3.08... 2) GNOME Mathematica $ mathematica 3) Mathematica SetDirectory SetDirectory[ ~/keisanki-1-2010/ ]
Mathematica (2/4) 4) {ti, x1(ti), x2(ti), x3(ti)} data data = ReadList[ data_file, {Real, Real, Real, Real}]; 5) - data {{t0, x1(t0), x2(t0), x3(t0)}, {t1, x1(t1), x2(t1), x3(t1)},... } Transpose( ) data datat datat = Transpose[data]; datat {{t0, t1, t2,... }, {x1(t0), x1(t1), x1(t2),... }, {x2(t0), x2(t1), x2(t2),...}, {x3(t0), x3(t1), x3(t2),... }}
Mathematica (3/4) 6) - datat 1 ( t, x1(t) ), ( t, x2(t) ), ( t, x3(t) ) datax = Transpose[ {datat[[1]], datat[[2]] } ]; datay = Transpose[ {datat[[1]], datat[[3]] } ]; dataz = Transpose[ {datat[[1]], datat[[4]] } ]; 3 ( x1(t), x2(t), x3(t) ) dataxyz dataxyz = Transpose[ {datat[[2]], datat[[3]], datat[[4]] } ];
Mathematica (4/4) 7) ListPlot t x1(t), x2(t), x3(t) gx = ListPlot[ datax, Joined->True, PlotRange->All] gy = ListPlot[ datay, Joined->True, PlotRange->All] gz = ListPlot[ dataz, Joined->True, PlotRange->All] 3 Show Show[gx, gy, gz] 3 x1(t), x2(t), x3(t) ListPointPlot3D[dataXYZ, BoxRatios->{1,1,1}, ViewPoint->{0.810, -2.475, 2.160}]
: 1) (2009) H21 1 2) Hirsch, M.W., Smale S., Devaney, R.L. (2007) 2 - -,