2014計算機実験1_1

Size: px

Start display at page:

Download "2014計算機実験1_1"

たつやなぐも
5 years ago
Views:

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

3 数学モデリングのプロセス問題点の抽出定義仮定数式化万有引力の法則 m すべての物体は引き合う r mm F =G 2 r M モデルの検証モデルによる説明将来予測解釈 F: 万有引力 (kg m s-2) G: 万有引力定数 (m s kg ) 解析数値計算 M: 地球の質量 (kg) により解を得る m: 落下する物質の質量 (kg) r: 2物質間の距離 (m) 画像引用:

4 法則と微分方程式ニュートンの運動方程式物体の加速度 a は物体の質量 m に ma = F Fの力で質量mの物体を押したら t秒後の速度vはどれくらい変わる? 反比例し物体に働く力 F に比例する速度加速度 dx v(t) = dt 2 dv d x a(t) = = 2 dt dt x(t): 位置 m m 2 d x 2 dt F =F 2階常微分方程式

5 ? x(0) mg ma(t) = mg g: t = 0 x = x(0), υ = υ(0) d 2 x t dt = g 2 x v(t) = gt + v(0) x(t) = 1 2 gt2 + v(0)t + x(0) t x(t) ( ) t

6 t y = y(t) yʹ, yʹ,... n f (t,y,y,y,...,y (n) )=0 f (t,y, dy dt,..., dn y dt n )=0 n t y 2

7 dy :, dt = f (t,y) y(t 0 )=y 0 y = y(t) f (t,y)=ky k = 0.5 y 0 = 0.5 dy dt = ky y(t) = y 0 e kt y t

8 ,, dy dx = f (x)g(y) dy dx = f y x =

9 dx dt = f (t, x) = lim t 0 x(t + t) x(t) t t Δt = ti+1 - ti = h f (t i, x i ) x i+1 t i+1 x i t i x i+1 x i + hf(t i, x i ) x(t) = x i, x(t + t) = x i+1 x 0 t i x 1, x 2, x 3,...

10 オイラー法のアルゴリズム区間 [t0, tn] を等間隔 h で分割しステップ ti から ti+1 の x の増加率を h f(ti, xi)で与える x dx dt while (t <= tn){ ti+1 = ti + h, xi+1 = xi + h f(ti, xi) } 解析解 x2 x1 x0 f(t0,x0) t0 f(t1,x1) t1 h 数値解 t2 t

11 x x x2 x1 x0 x2 x1 x0 f(t0,x0) f(t 1,x1) t t t0 t1 t2 t0 t1 t2 h

12 オイラー法の精度 x(t)をt = ti の周りでテイラー展開するとテイラー展開による多項式での近似オイラー法 h2 x(ti+1 ) = x(ti + h) = x(ti ) + hx (ti ) + x (ti ) + 2 x(ti+1 ) = x(ti ) + h f (ti, xi ) 誤差オイラー法はテイラー展開の1次の項までしか考慮しない区間[ t0, t0 + T ]を刻み幅hで分割した時の終端 t0 + T での誤差の累積 2 h 2 T h = x (ti ) h 2 x (ti ) 1回のステップにおける誤差全ステップ数刻み幅 h を 1/10 にすれば累積誤差も約 1/10 になる (ただし計算量は 10倍になる)

13 (, 1 ) #define DT 0.01 /* */ #define STEP_MAX 1000 /* DT*STEP_MAX = 10.0 */ double fn(double, double); /* */ void euler(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; } } void euler(double t, double x, double *x_out, double h){... } double fn(double t, double x){... } dx dt = f (t,x) t = 0 t = DT*STEP_MAX

14 (, 1 ) euler ti x(ti), h ti+1 = ti + h x(ti+1) ti x(ti) x(ti+1) void euler(double t, double x, double *x_out, double h){ double tmp; tmp = x + h*fn(t,x); *x_out = tmp; } x_out double fn(double t, double x){... } tmp ti x(ti) x(ti+1) h

15 Mathematica (1/2) Mathematica, 1) ti x(ti) ti x(ti) 2) GNOME Mathematica $ mathematica 3) Mathematica SetDirectory SetDirectory[ ~/keisanki / ]

16 Mathematica (2/2) 4) ReadList[ file, {type, type}] file type ( {x (ti), y (x(ti))} ) ( ) data data = ReadList[ data_file, {Real, Real}] 5) ListPlot[list] ( ) PlotRange->{0, 10} 0 10 ListPlot[ data, Joined->True, PlotRange->All ]

17 Mathematica Mathematica C fig1=listplot[ data, Joined->True, PlotRange- >All ] fig2=plot[exp[t], {t, 0, 10}, PlotStyle ->Red] fig3=show[fig1, fig2] 1. fig1 2. e t 0 t 10 fig2 3. fig2 fig2 fig3

18 fig3 dx = x(1 x) dt x(0) = x t 2.5 x t h 0.5 h 0.05

19 Mathematica Mathematica C deq = { x [t] == (1 - x[t])x[t], x[0] == 0.01} sol = NDSolve[ deq, x[t], {t, 0, 10}] Plot[ Evaluate[ x[t]/.sol ], {t, 0, 10}] x[t]/.sol/.t->5 1. deq Mathematica =, == 2. deq 0 t 10 sol 3. t /. 4. x(5)

20 1 (1) (2) (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

21 Q1 A1 fprintf ( ) FILE *fp fp = fopen( data.txt, w ); /* */ fprintf(fp, %f %f n, t, x); fclose(fp); /* */

22 Q2 C sin(t) A2 2 #include <math.h> - lm Q3 Mathematica sin(t) A3 Sin[t] S

23 ( ) 1

24 (, ) #define DIM 2 /* ( ) */ #define DT 0.01 /* */ #define STEP_MAX 1000 /* DT*STEP_MAX = 10.0 */ void derives(double, double[], double[], int); /* */ void euler(double, double[], double[], int, 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; } } void euler(double t, double x, double *x_out, double h){... } double fn(double t, double x){... } t = 0 t = DT*STEP_MAX dx dt = f (t, x)

25 (, ) derives ti x(ti), n derivatives[] ti x(ti) dx/dt or f(t, x) void derives(double t, double x[], double derivatives[],int n){ } derivatives[] n

26 (, ) void euler(double t, double x[], double x_out[], int n, double h){ int i; double x_tmp[dim], derivatives[dim]; derivs(t, x, derivatives, n); /* */ } euler ti x(ti), n, h ti+1 = ti + h x(ti+1) ti x(ti) x(ti+1) n h for(i=0; i<dim; i++) /* Euler */ x_tmp[i] = x[i] + h*derivatives[i]; for(i=0; i<dim; i++) x_out[i] = x_tmp[i]; /* */

27 double fn(double t, double x){ } 1 fn derives f (t, x) void derives(double t, double x[], double derivatives[],int n){ derivatives[] }

28 ...??? 2 double fn1(double t, double x1, double x2){ } double fn2(double t, double x1, double x2){ } f1(t, x1, x2) f2(t, x1, x2)

29 2 ε 1.0, 0.5, t, x1, x2 Mathematica 2. x1, x2 Mathematica x(0), x(1), x(2),... dx 1 dt dx 2 dt = x 1 1 = x 1 3 x3 1 x 2 Van der Pol oscillator equation

30 Mathematica (1/4) Mathematica, t i x1(ti) x2(ti) 1) ) Mathematica % /usr/local/bin/mathematica & 3) Mathematica SetDirectory SetDirectory[ ~/keisanki / ]

31 Mathematica (2/4) 4) {ti, x1(ti), x2(ti)} data data = ReadList[ data_file, {Real, Real, Real}]; 5) - data {{t 0, x 0, y 0 }, {t 1, x 1, y 1 }, {t 2, x 2, y 2 },... } Transpose( ) data datat datat = Transpose[data]; datat {{t 0, t 1, t 2,... }, {x 0, x 1, x 2,... }, {y 0, y 1, y 2,...}}

32 Mathematica (3/4) 4) - datat 1 2 ( t, x1(t) ) datax datax = Transpose[ {datat[[1]], datat[[2]] } ]; ( t, y(t) ) datay ( x(t), y(t) ) dataxy datay = Transpose[ {datat[[1]], datat[[3]] } ]; dataxy = Transpose[ {datat[[2]], datat[[3]] } ];

33 Mathematica (4/4) 5) ListPlot t x(t), y(t) gx = ListPlot[ datax, Joined->True, PlotRange->All] gy = ListPlot[ datay, Joined->True, PlotRange->All] 2 Show Show[gx, gy] x(t), y(t) ListPlot[ dataxy, Joined->True, PlotRange->All]

34 Mathematica deq2 = { x1 [t] == x2[t] x1[t]^2, x2 [t] == 2 x1[t] x2[t], x1[0] == x2[0] == 1} sol2 = NDSolve[ deq2, {x1[t], x2[t]}, {t, 0, 20}] Plot[ Evaluate[ x1[t]/.sol2 ], {t, 0, 20}] Plot[ Evaluate[ x2[t]/.sol2 ], {t, 0, 20}] ParametricPlot[ Evaluate[ {x1[t], x2[t]}/.sol2 ], {t, 0, 20}] 1. deq 2. deq2 0 t 20 sol2 3. t 4. ( x 1 (t), x 2 (t) )

35 LaTeX 1) Mathematica EPS EPS: Encapsulated PostScript Mathematica Edit ---> Save Selection As ---> EPS 2) LaTex LaTeX \usepackage{graphicx} EPS.eps \includegraphics{graph.eps} ( ) 5cm \begin{cetner} \includegraphics[width = 5cm]{graph.eps} \end{center}

36 Tips ( ) 2.5 x dx/dt = x(1-x) ( ) h = 0.5 ( ) x(0) = t

37 / (2002) ( ) / (1994) / (2007)

38 : 1) Burghes, D. and Borrie, M. (1990), 2) (2009) H21 1 3) (2000) Mathematica 2,

OK (S) vncviewer UNIX EDS vncviewer : VNC server: eds.efc.sec.eng.shizuoka.ac.jp:51 OK 2

Mathematica I (2001 5 31, 6 7 ) UNIX EDS vncviewer Internet Exploler http://www.efc.sec.eng.shizuoka.ac.jp/admin/pubsoft/ vncviewer.exe : 1 OK (S) vncviewer UNIX EDS vncviewer : VNC server: eds.efc.sec.eng.shizuoka.ac.jp:51