H

Size: px

Start display at page:

Download "H"

そうさだい
4 years ago
Views:

1 H seto@ics.nara-wu.ac.jp

2 Euler xi+1 x ( xi f(ti, xi) ) t ti ti+1 h

3 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

4 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 )

5 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 k k k 4 6

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 k k k 4 6 ti ti + h/2 ti+1 t h

7 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 k k k 4 6 ti ti + h/2 ti+1 t h

8 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 k k k 4 6 ti ti + h/2 ti+1 t h

9 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 k k k 4 6 ti ti + h/2 ti+1 t h

10 (, 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){... }

11 (, 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 k k k 4 6 x i+1 = x i + k k k k 4 6

12 H (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

13 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) = t 10 (3) dx dt = x 1 + t + cos t, x(0) = 1 0 t 100

14 (1) dx dt = t 1 x 2, x(0) = : 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

15 1

16 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 k k k 4 6

17 (, ) #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){... }

18 (, ) 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]); }

19 (, ) 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 k k k 4 6 x(ti) x(ti+1) k 1 ~k 4 derives[] x_out[] h

20 -

21 Lotka-Volterra dx dt = ax dy dt = bxy cy + dxy a: b: c: d:

22 5 a = 1, b = 0.03, c = 1, d = x(0), y(0)

23 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:

24 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

25 Mathematica (1/4) Mathematica, ti x1(ti) x2(ti) x3(ti) 1) ) GNOME Mathematica $ mathematica 3) Mathematica SetDirectory SetDirectory[ ~/keisanki / ]

26 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),... }}

27 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]] } ];

28 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.160}]

29 : 1) (2009) H21 1 2) Hirsch, M.W., Smale S., Devaney, R.L. (2007) 2 - -,

2014計算機実験1_1

2014計算機実験1_1 H26 1 1 1 seto@ics.nara-wu.ac.jp 数学モデリングのプロセス問題点の抽出定義仮定数式化万有引力の法則 m すべての物体は引き合う r mm F =G 2 r M モデルの検証モデルによる説明将来予測解釈 F: 万有引力 (kg m s-2) G: 万有引力定数 (m s kg ) 解析数値計算 M: 地球の質量 (kg) により解を得る m: 落下する物質の質量