£Ã¥×¥í¥°¥é¥ß¥ó¥°(2018) - Âè11²ó – ½ÉÂê£²¤Î²òÀâ¡¤±é½¬£² –

(1)

.

Ｃプログラミング

(2018)

第

11 回

–

宿題２の解説，演習２

–

松田七美男

2018

年

12 月

13 日

(2)

宿題

2

：概要（１）

.

空気の粘性抵抗を受ける質点の運動

.

一様な重力加速度

g

の下，速さ

v

に比例する空気の粘性抵抗が作用

する場合の運動方程式を整理すると以下のような独立した連立常微

分方程式が得られます．

dvx

dt

=

−bvx,

dvy

dt

=

−g − bvy

(1)

ここに，

b

は定数です．原点から初速度

v

0 ，水平面から傾角

θ

で投

げ上げられた質点の運動は以下のように求まります

x(t) =

v

0 cos θ

b

(

1 − e

−bt

)

(2)

y(t) =

1 b

(

v

0 sin θ +

g

b

) (

1 − e

−bt

)

−

g

b

t

(3)

後藤憲一，他：

「

基礎

_{力学演習」例題}

11 (

共立出版

) p.14

(3)

宿題

2

：概要（２）

.

1 y

方向の微分方程式を

4 次のルンゲ・クッタ法で数値計

算し，質点が水平面に落下する時間

t

m

を求めなさい．具

体的には，投げ上げて（

y > 0

）から初めて

y < 0

となる

時間

t

m1

を求めればよい．刻み幅は

h = 0.001

程度とし

なさい．

.

2 t

_m

は式

(3)

の左辺を

0 とした方程式の解であり，区間

[ t

m1

− h, t

m1

]

に真の値があると考えられる．２分法によ

る数値解

t

m2

を求めなさい．

.

3 前項の

t

_m2

の前後の値を式

(3)

に代入して，

y(t

_m2

)

が最

も

0 に近いことを確かめなさい．

.

4 到達距離

x(t

_m

)

を算出しなさい．

(4)

宿題

2

：要件

I

4 次ルンゲ・クッタ法については，独自のコードを作成し

ないこと

．すなわち，講義で用いた

rk4fix.o

に含まれ

る関数

rk4fixv6

を呼び出すこととし，

main

関数と連立

常微分方程式の定義関数及び

x(t), y(t)

を表す関数を作成

しなさい．

I

諸定数は，

g = 9.8 m

· s

−2

, b = 1.0 m

−1

, v0

= 9.8 m

· s

−1

とし，傾角

θ

を「学籍番号

mod 60 +20

」度（即ち，

20 から

79 までの整数値）としなさい．

(5)

宿題

2

：実行画面例

θ = 26

度のときの実行画面例を示します．刻み幅

h

を

0.001 ，

2 分法の打ち切り誤差を

1.0 × 10

−13

にしています．

.

$ ./hw2018A-2

Th = 26

tm = 0.7780000000000

t = 0.7770948429848, y = 7.238654e-13

t = 0.7770948429849, y = 3.916867e-13

t = 0.7770948429850, y = 5.950795e-14

t = 0.7770948429851, y = -2.726708e-13

t = 0.7770948429852, y = -6.039613e-13

xm = 4.7587109588997

(6)

解答例

#include "rk4fix.h" #define B 1.0 #define G 9.8 #define V0 9.8 #define EPS 1E-13

double f(double t, double g, double b, double v) {

return (v + g / b) / b *

(1 - exp(-b * t)) - g / b * t; }

vec6 mov(double t, double r[6]) { vec6 ret = { {0} }; ret.f[1] = r[4]; ret.f[4] = -G - B * r[4]; return ret; }

int main(int argc, char **argv) {

int deg = 66, i;

double t = 0, r[6] = { 0 }, h = 0.001; double tL, tU, dt, tm, vx, vy, Th; if (argc > 1) deg = atof(argv[1]); deg = (deg % 60) + 20; printf("Th = %d\n", deg); Th = M_PI / 180 * deg; vy = V0 * sin(Th); vx = V0 * cos(Th); r[4] = vy; do { rk4fixv6(mov, t, r, h); t += h; } while (r[1] >= 0.0); printf("tm = %.13f\n", t); tm = t; tL = tm - h; tU = tm; do { dt = tU - tL; t = (tU + tL) / 2;

if (f(t, G, B, vy) * f(tU, G, B, vy) > 0) tU = t; else tL = t;

} while (dt > EPS); tm = t;

for (i = -2; i < 3; i++) { t = tm + i * EPS;

printf("t = %.13f, y = %.6e\n", t, f(t, G, B, vy)); } printf("xm = %.13f\n", f(tm, 0, B, vx)); return 0; }

.

定数の定義

.

微分方程式の記述

.

0 > y を求めるループ

.

二分法による 0 点の検索

.

x(tm

) の表示

(7)

解答例

(1 - exp(-b * t)) - g / b * t; }

int deg = 66, i;

for (i = -2; i < 3; i++) { t = tm + i * EPS;

.

定数の定義

.

微分方程式の記述

.

0 > y を求めるループ

.

二分法による 0 点の検索

.

x(tm

) の表示

(8)

解答例

(1 - exp(-b * t)) - g / b * t; }

int deg = 66, i;

for (i = -2; i < 3; i++) { t = tm + i * EPS;

.

定数の定義

.

微分方程式の記述

.

0 > y を求めるループ

.

二分法による 0 点の検索

.

x(tm

) の表示

(9)

解答例

(1 - exp(-b * t)) - g / b * t; }

int deg = 66, i;

for (i = -2; i < 3; i++) { t = tm + i * EPS;

.

定数の定義

.

微分方程式の記述

.

0 > y を求めるループ

.

二分法による 0 点の検索

.

x(tm

) の表示

(10)

解答例

(1 - exp(-b * t)) - g / b * t; }

int deg = 66, i;

for (i = -2; i < 3; i++) { t = tm + i * EPS;

.

定数の定義

.

微分方程式の記述

.

0 > y を求めるループ

.

二分法による 0 点の検索

.

x(tm

) の表示

(11)

解答例

(1 - exp(-b * t)) - g / b * t; }

int deg = 66, i;

for (i = -2; i < 3; i++) { t = tm + i * EPS;

.

定数の定義

.

微分方程式の記述

.

0 > y を求めるループ

.

二分法による 0 点の検索

.

x(tm

) の表示

(12)

演習

2

：概要

.

sin

関数の自作

.

正弦関数

sin(x)

を計算する関数

double

_mysin(

double

_x)

_{を作成しなさ}

い．ただし，

x

の値に応じて以下の場合分けを行って算出するもの

とします．

.

1 −π < x < π

の値に対しては，以下の級数展開に基づき，項を

k = 16

で打ち切って計算しなさい．

sin x =

∞

∑

k=1

(

−1)

k

−1

x

2k

−1

(2k

− 1)!

= x

∞

∑

k=1

(

−x

2 )

k

−1

(2k

− 1)!

.

2 上記範囲外については変換

sin(x

′

± 2π) = sin(x

′

)

により

−π < x

′

< π

の範囲の

x

′

を求めて算出しなさい．

(13)

演習

2

：要件と実行画面

.

1 Horner

法により，

k = 16

の項から計算を始めるように

組む．

.

2 自分の学籍番号を

X

として，

X

× 10

−2

, X

× 10

−1

, X

に

対して実行しなさい．

X = 62

の場合の実行画面を示します．

.

$ ./ex2018A-2

mysin( 62.00) = -7.391806966492211e-01

sin( 62.00) = -7.391806966492228e-01

mysin(

6.20) = -8.308940281749616e-02

sin(

6.20) = -8.308940281749640e-02

mysin(

0.62) = +5.810351605373051e-01

sin(

0.62) = +5.810351605373051e-01

(14)

演習

2

提出に関する注意

I

ドメイン

.waseda.jp

の送信アドレスから提出しなさい

I

宛先は

matuda-namio@aoni.waseda.jp

I

件名には，学生番号（

英数半角

），名前，演習

2 をこの順

に記述しなさい．

例）

1Y17B999

，△○ ×□，演習

2

I

〆切は，

12 月

13 日

10:30

I

ソース（

ex2018A-2.c

）と実行結果の両方が必要．

添付

不可

．

I

ソースには空白を挿入したり，字下げをほどこすなど読

みやすくなる工夫をしなさい．