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: 落下する物質の質量 (kg) 3-2 -1 r: 2物質間の距離 (m) 画像引用: http://www.matumura-ringo.com/info.html

法則と微分方程式 ニュートンの運動方程式 物体の加速度 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階常微分方程式

? 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

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

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 2 1 0 1 2 2 1 0 1 2 t

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

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,...

オイラー法のアルゴリズム 区間 [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

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

オイラー法の精度 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倍になる)

(, 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

(, 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

Mathematica (1/2) Mathematica, 1) ti x(ti) 0.000 1.00 0.001 1.02 0.002 1.05... ti x(ti) 2) GNOME Mathematica $ mathematica 3) Mathematica SetDirectory SetDirectory[ ~/keisanki-1-2010/ ]

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 ]

Mathematica Mathematica C 1. 2. 3. 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

fig3 dx = x(1 x) dt x(0) = 2 2.5 x 2.0 1.5 1.0 0.5 0 2 4 6 8 10 t 2.5 x 2.0 1.5 1.0 0.5 0 2 4 6 8 10 t h 0.5 h 0.05

Mathematica Mathematica C 1. 2. 3. 4. 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)

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

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

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

( ) 1

(, ) #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)

(, ) 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

(, ) 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]; /* */

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

...??? 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)

2 ε 1.0, 0.5, 0.1 1-2. 1. 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

Mathematica (1/4) Mathematica, t i x1(ti) x2(ti) 1) 0.000 1.00 2.00 0.001 1.02 1.97 0.002 1.05 1.94... 2) Mathematica % /usr/local/bin/mathematica & 3) Mathematica SetDirectory SetDirectory[ ~/keisanki-1-2010/ ]

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

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

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]

Mathematica 1. 2. 3. 4. 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) )

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}

Tips ( ) 2.5 x 2.0 1.5 1.0 1 dx/dt = x(1-x) ( ) h = 0.5 ( ) x(0) = 2 0.5 0 2 4 6 8 10 t

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

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