.
Ｃプログラミング (2014)
.
第 11 回
– 宿題２の解説，演習２ –
松田七美男
2014 年 12 月 18 日
. 宿題 2：概要
.
強制振動の過渡応答
.
一次元の強制振動の初期値問題
1
1
x
¨ = −x − x˙ + cos t, x(0) = 0, x(0)
˙
= v0
4
2
を，4 次のルンゲ-クッタ法で数値解析し変位 x(t) の絶対値 |x| の
最大値を数値的に求めなさい．
.
. 初期速度 v0 を「学籍番号÷ 10」，としなさい．
例）1Y13B078 ⇒ v0 = 78/10 = 7.8
1
. 刻み幅 h は始め 0.01 とし，|x| が最大となる時間 t1 を 0.01
の精度で見つける．続いて t1 − 0.01 から，刻み幅を細かく
h = 1 × 10−6 に変えて再計算して |x| が最大になる時間を
1 × 10−6 の精度で見つける．
2
. 解析する時間の最大値 Tmax は，30 程度としなさい．
3
. 宿題 2：振動の様子
v0 = 1, 2, 5, 10 とした場合の振動の様子を図に示します．変位の
絶対値が最も大きくなるのは t = 2, 5 付近の初めの 2 つピークの
どちらかであることが一目瞭然です．
10
v0 = 1
2
5
10
8
6
x(t)
4
2
0
-2
-4
-6
-8
0
5
10
15
t
.
20
25
30
. 宿題 2：振動の様子
v0 = 1, 2, 5, 10 とした場合の振動の様子を図に示します．変位の
絶対値が最も大きくなるのは t = 2, 5 付近の初めの 2 つピークの
どちらかであることが一目瞭然です．
10
v0 = 1
2
5
10
8
6
x(t)
4
2
0
-2
-4
-6
-8
0
5
10
15
t
.
20
25
30
. 解答例
#include
#include
#include
#include
#define
#define
#define
#define
<stdio.h>
<stdlib.h>
<math.h>
"rk4fix.h"
Tmax
H
HH
V0
while (t < Tmax){
xtmp = r[0];
vtmp = r[3];
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
t1 = t;
xm1 = xtmp;
vm1 = vtmp;
}
}
30
0.01
1e-6
7.8
vec6 mov(double t, double r[6])
{
vec6 ret;
ret.f[0]
ret.f[1]
ret.f[2]
ret.f[3]
ret.f[4]
ret.f[5]
=
=
=
=
=
=
h = HH;
tm = t1;
t = t1 - H;
r[0] = xm1;
r[3] = vm1;
r[3];
0;
0;
-r[0] - 0.25*r[3] + cos(0.5*t);
0;
0;
while (t <= t1 + 2*H){
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
tm = t;
}
}
printf("xmax = %.8f at tm = %.7f\n", xmax, tm);
return ret;
}
int main(int argc, char **argv)
{
double t = 0, h = H, t1 = 0, tm = 0,
xtmp, xm1, vtmp, vm1, xmax = 0,
r[6] = {0, 0, 0, V0, 0, 0};
return 0;
.
}
. 解答例
#include
#include
#include
#include
#define
#define
#define
#define
<stdio.h>
<stdlib.h>
<math.h>
"rk4fix.h"
Tmax
H
HH
V0
30
0.01
1e-6
7.8
while (t < Tmax){
xtmp = r[0];
vtmp = r[3];
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
t1 = t;
xm1 = xtmp;
vm1 = vtmp;
}
}
定数の定義
vec6 mov(double t, double r[6])
{
vec6 ret;
ret.f[0]
ret.f[1]
ret.f[2]
ret.f[3]
ret.f[4]
ret.f[5]
=
=
=
=
=
=
h = HH;
tm = t1;
t = t1 - H;
r[0] = xm1;
r[3] = vm1;
r[3];
0;
0;
-r[0] - 0.25*r[3] + cos(0.5*t);
0;
0;
while (t <= t1 + 2*H){
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
tm = t;
}
}
printf("xmax = %.8f at tm = %.7f\n", xmax, tm);
return ret;
}
int main(int argc, char **argv)
{
double t = 0, h = H, t1 = 0, tm = 0,
xtmp, xm1, vtmp, vm1, xmax = 0,
r[6] = {0, 0, 0, V0, 0, 0};
return 0;
.
}
. 解答例
#include
#include
#include
#include
#define
#define
#define
#define
<stdio.h>
<stdlib.h>
<math.h>
"rk4fix.h"
Tmax
H
HH
V0
30
0.01
1e-6
7.8
while (t < Tmax){
xtmp = r[0];
vtmp = r[3];
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
t1 = t;
xm1 = xtmp;
vm1 = vtmp;
}
}
定数の定義
vec6 mov(double t, double r[6])
{
vec6 ret;
ret.f[0]
ret.f[1]
ret.f[2]
ret.f[3]
ret.f[4]
ret.f[5]
=
=
=
=
=
=
h = HH;
tm = t1;
t = t1 - H;
r[0] = xm1;
r[3] = vm1;
r[3];
微分方程式の記述
0;
0;
-r[0] - 0.25*r[3] + cos(0.5*t);
0;
0;
while (t <= t1 + 2*H){
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
tm = t;
}
}
printf("xmax = %.8f at tm = %.7f\n", xmax, tm);
return ret;
}
int main(int argc, char **argv)
{
double t = 0, h = H, t1 = 0, tm = 0,
xtmp, xm1, vtmp, vm1, xmax = 0,
r[6] = {0, 0, 0, V0, 0, 0};
return 0;
.
}
. 解答例
#include
#include
#include
#include
#define
#define
#define
#define
<stdio.h>
<stdlib.h>
<math.h>
"rk4fix.h"
Tmax
H
HH
V0
30
0.01
1e-6
7.8
while (t < Tmax){
xtmp = r[0];
前の値の保存
vtmp = r[3];
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
t1 = t;
xm1 = xtmp;
vm1 = vtmp;
}
}
定数の定義
vec6 mov(double t, double r[6])
{
vec6 ret;
ret.f[0]
ret.f[1]
ret.f[2]
ret.f[3]
ret.f[4]
ret.f[5]
=
=
=
=
=
=
h = HH;
tm = t1;
t = t1 - H;
r[0] = xm1;
r[3] = vm1;
r[3];
微分方程式の記述
0;
0;
-r[0] - 0.25*r[3] + cos(0.5*t);
0;
0;
while (t <= t1 + 2*H){
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
tm = t;
}
}
printf("xmax = %.8f at tm = %.7f\n", xmax, tm);
return ret;
}
int main(int argc, char **argv)
{
double t = 0, h = H, t1 = 0, tm = 0,
xtmp, xm1, vtmp, vm1, xmax = 0,
r[6] = {0, 0, 0, V0, 0, 0};
return 0;
.
}
. 解答例
#include
#include
#include
#include
#define
#define
#define
#define
<stdio.h>
<stdlib.h>
<math.h>
"rk4fix.h"
Tmax
H
HH
V0
30
0.01
1e-6
7.8
while (t < Tmax){
xtmp = r[0];
前の値の保存
vtmp = r[3];
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
t1 = t;
xm1 = xtmp; 最大値の探索と
vm1 = vtmp; 相応する値の保存
}
}
定数の定義
vec6 mov(double t, double r[6])
{
vec6 ret;
ret.f[0]
ret.f[1]
ret.f[2]
ret.f[3]
ret.f[4]
ret.f[5]
=
=
=
=
=
=
h = HH;
tm = t1;
t = t1 - H;
r[0] = xm1;
r[3] = vm1;
r[3];
微分方程式の記述
0;
0;
-r[0] - 0.25*r[3] + cos(0.5*t);
0;
0;
while (t <= t1 + 2*H){
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
tm = t;
}
}
printf("xmax = %.8f at tm = %.7f\n", xmax, tm);
return ret;
}
int main(int argc, char **argv)
{
double t = 0, h = H, t1 = 0, tm = 0,
xtmp, xm1, vtmp, vm1, xmax = 0,
r[6] = {0, 0, 0, V0, 0, 0};
return 0;
.
}
. 解答例
#include
#include
#include
#include
#define
#define
#define
#define
<stdio.h>
<stdlib.h>
<math.h>
"rk4fix.h"
Tmax
H
HH
V0
30
0.01
1e-6
7.8
while (t < Tmax){
xtmp = r[0];
前の値の保存
vtmp = r[3];
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
t1 = t;
xm1 = xtmp; 最大値の探索と
vm1 = vtmp; 相応する値の保存
}
}
定数の定義
vec6 mov(double t, double r[6])
{
vec6 ret;
ret.f[0]
ret.f[1]
ret.f[2]
ret.f[3]
ret.f[4]
ret.f[5]
=
=
=
=
=
=
h = HH;
tm = t1;
t = t1 - H;
r[0] = xm1;
r[3] = vm1;
r[3];
微分方程式の記述
0;
0;
-r[0] - 0.25*r[3] + cos(0.5*t);
0;
0;
while (t <= t1 + 2*H){
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
tm = t;
}
}
printf("xmax = %.8f at tm = %.7f\n", xmax, tm);
return ret;
}
int main(int argc, char **argv)
{
double t = 0, h = H, t1 = 0, tm = 0,
xtmp, xm1, vtmp, vm1, xmax = 0,
r[6] = {0, 0, 0, V0, 0, 0};
return 0;
.
時間を戻して値設定
}
. 解答例
#include
#include
#include
#include
#define
#define
#define
#define
<stdio.h>
<stdlib.h>
<math.h>
"rk4fix.h"
Tmax
H
HH
V0
30
0.01
1e-6
7.8
while (t < Tmax){
xtmp = r[0];
前の値の保存
vtmp = r[3];
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
t1 = t;
xm1 = xtmp; 最大値の探索と
vm1 = vtmp; 相応する値の保存
}
}
定数の定義
vec6 mov(double t, double r[6])
{
vec6 ret;
ret.f[0]
ret.f[1]
ret.f[2]
ret.f[3]
ret.f[4]
ret.f[5]
=
=
=
=
=
=
h = HH;
tm = t1;
t = t1 - H;
r[0] = xm1;
r[3] = vm1;
r[3];
微分方程式の記述
0;
0;
-r[0] - 0.25*r[3] + cos(0.5*t);
0;
0;
t1
時間を戻して値設定
H
H
while (t <= t1 + 2*H){ 2H の範囲にある
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
tm = t;
}
}
printf("xmax = %.8f at tm = %.7f\n", xmax, tm);
return ret;
}
int main(int argc, char **argv)
{
double t = 0, h = H, t1 = 0, tm = 0,
xtmp, xm1, vtmp, vm1, xmax = 0,
r[6] = {0, 0, 0, V0, 0, 0};
return 0;
.
max
}
. 解答例
#include
#include
#include
#include
#define
#define
#define
#define
<stdio.h>
<stdlib.h>
<math.h>
"rk4fix.h"
Tmax
H
HH
V0
30
0.01
1e-6
7.8
while (t < Tmax){
xtmp = r[0];
前の値の保存
vtmp = r[3];
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
t1 = t;
xm1 = xtmp; 最大値の探索と
vm1 = vtmp; 相応する値の保存
}
}
定数の定義
vec6 mov(double t, double r[6])
{
vec6 ret;
ret.f[0]
ret.f[1]
ret.f[2]
ret.f[3]
ret.f[4]
ret.f[5]
=
=
=
=
=
=
h = HH;
tm = t1;
t = t1 - H;
r[0] = xm1;
r[3] = vm1;
r[3];
微分方程式の記述
0;
0;
-r[0] - 0.25*r[3] + cos(0.5*t);
0;
0;
t1
時間を戻して値設定
H
H
while (t <= t1 + 2*H){ 2H の範囲にある
rk4fixv6(mov, t, r, h);
t += h;
if (fabs(r[0]) > xmax) {
xmax = fabs(r[0]);
tm = t;
10−6 の精度で最大値の探索
}
}
printf("xmax = %.8f at tm = %.7f\n", xmax, tm);
return ret;
}
int main(int argc, char **argv)
{
double t = 0, h = H, t1 = 0, tm = 0,
xtmp, xm1, vtmp, vm1, xmax = 0,
r[6] = {0, 0, 0, V0, 0, 0};
return 0;
.
max
}
. 演習 2
.
平方根関数の自作
.
√
平方根 x を計算する関数 double mysqrt( double x) を作成しなさい．
ただし，x の値に応じて以下の場合分けを行って算出するものと
します．
.
. 0.01 ≤ x ≤ 100 の値に対しては，方程式 f (ξ) = ξ 2 − x = 0 の解を
ニュートン法 (試行回数 8) で求めなさい．これは，次の漸化式を
計算することになります．
(
)
1
x
ξ0 = x, ξk+1 =
ξk +
2
ξk
1
√
. 上記範囲外については x × (10)2n = 10n √x と変換することに
より算出しなさい．
.3 引数 x の平方根を，自作 mysqrt() と数学関数ライブラリの sqrt()
について比較するプログラムを作成し，自分の学籍番号を X とし
て，X, X × 10−4 , X × 105 に対して実行しなさい．
2
. 演習 2 提出に関する注意
Waseda-net Web メールサービスを用いて提
出しなさい．
宛先は [email protected]
件名には，学生番号（英数半角），名前，演習 2 をこの
順に記述しなさい．
例）1Y13B999，△○ ×□，演習 2
〆切は，12 月 18 日 10:30
送信控えの保存を忘れずに
ソース（inv mat.c と lineq invmat.c）と実行
結果の両方が必要．添付不可．
ソースには空白を挿入したり，字下げをほど
こすなど読みやすくなる工夫をしなさい．

Download Report