微分方程式の解法
とその可視化
Part2
14年11月5日水曜日
常微分方程式の数値解法
前回までの復習
オイラー法とは 微分を１次の差分で近似する．
近似の微分
本当の微分
14年11月5日水曜日
常微分方程式の数値解法
ルンゲ・クッタ法では微分を幾つかの項の和として近似する．
数値解析において常微分方程式の近似解を求める
一連の方法で，以下の方法だけを指す言葉ではない．
14年11月5日水曜日
常微分方程式の数値解法
4段4次ルンゲ・クッタ法
現時点の値からk1を求める．
14年11月5日水曜日
常微分方程式の数値解法
4段4次ルンゲ・クッタ法
k1をdt/2の地点まで延長し，そこから
k2を求める．
14年11月5日水曜日
常微分方程式の数値解法
4段4次ルンゲ・クッタ法
k2を最初の地点に戻し，そこから
k3を求める．
14年11月5日水曜日
常微分方程式の数値解法
4段4次ルンゲ・クッタ法
k3を最初の地点に戻し，そこから
dt離れたk4を求める．
14年11月5日水曜日
常微分方程式の数値解法
4段4次ルンゲ・クッタ法
得られたk1~k4を重み付きで足し合
わせる．
公式完成！
14年11月5日水曜日
常微分方程式の数値解法
オイラー法 V.S. 4段4次ルンゲ・クッタ法
オイラー法
4段4次ルンゲ・クッタ法
直感的でわかりやすい公式
やや複雑な公式
誤差は１次のオーダー
プログラミングも楽
初心者向き
誤差は4次のオーダー
(力のある人は証明してみよ．)
常微分方程式の解法としては
様々な意味で最優秀な公式
余分に４つの記憶領域を使う
14年11月5日水曜日
常微分方程式の数値解法
4段4次ルンゲ・クッタ法の適用例
こんな感じ
問題
4段4次ルンゲ・クッタ法
Exercise!
上記の問題をRK法を用いて解き，
可視化せよ．(ヒント：オイラー法
のプログラムを少し直すだけで完成
できるはず)
14年11月5日水曜日
常微分方程式の数値解法
#include <stdio.h>
#include <math.h>
int
main(void)
{
int i;
double dt = 0.01;
double xN,xN_new;
double k1,k2,k3,k4;
4段4次ルンゲ・クッタ法の適用例
FILE *gp;
gp = popen("gnuplot -persist","w");
fprintf(gp, "plot '-' with lines \n");
xN = 0.01;
for (i = 0; i < 1000; i++)
{
fprintf(gp,"%f %f\n", i * dt, xN);// データの書き込み
k1 = xN;
k2 = xN + dt / 2 * k1;
k3 = xN + dt / 2 * k2;
k4 = xN + dt * k3;
xN_new = xN + dt * (k1 + 2 * k2 + 2 * k3 + k4) /6;//計算
ひみつ
xN = xN_new;
}
fprintf(gp,"e\n");
fflush(gp);
pclose(gp);
return 0;
}
14年11月5日水曜日
常微分方程式の数値解法
Exercise！RK法で計算＆可視化をせよ．また厳密解と比較せよ．
scanf関数
main2.c
main3.c
main4.c
で読み込
めるよう
に．
14年11月5日水曜日
常微分方程式の数値解法
２変数の場合
問題
Exercise!
上記の問題をRK法を用いて解き，可視化せよ．
計算の順番に気をつけること！
14年11月5日水曜日
常微分方程式の数値解法
#include <stdio.h>
#include <math.h>
int
main(void)
{
int i;
double
double
double
double
double
多変数の場合
dt = 0.01;
xN,xN_new;
yN,yN_new;
kx1,kx2,kx3,kx4;
ky1,ky2,ky3,ky4;
xN = 1.0;
yN = 0.0;
for (i = 0; i < 10000; i++)
{
printf("%f %f %f\n",i * dt, xN, yN);
//printf("%f\n", xN * xN + yN * yN);
kx1
kx2
kx3
kx4
=
=
=
=
-yN;
-(yN + dt / 2 * ky1);
-(yN + dt / 2 * ky2);
-(yN + dt * ky3);
ky1
ky2
ky3
ky4
=
=
=
=
xN;
xN + dt / 2 * kx1;
xN + dt / 2 * kx2;
xN + dt * kx3;
xN_new = xN + dt * (kx1 + 2 * kx2 + 2 * kx3 + kx4) / 6;
yN_new = yN + dt * (ky1 + 2 * ky2 + 2 * ky3 + ky4) / 6;
xN = xN_new;
yN = yN_new;
}
return 0;
}
14年11月5日水曜日
main5.c
cf: Euler
cf: Euler
常微分方程式の数値解法
オイラー法 V.S. 4段4次ルンゲ・クッタ法 （誤差対決）
オイラー法が真の解からどれほど離れているかを表す指標
RK法が真の解からどれほど離れているかを表す指標
どちらの方法も の時０に近づくはず！
オイラー法 V.S. 4段4次ルンゲ・クッタ法 （誤差対決）
14年11月5日水曜日
常微分方程式の数値解法
オイラー法 V.S. 4段4次ルンゲ・クッタ法 （誤差対決）
#include <stdio.h>
#include <math.h>
#define
N
(100)
int
main(void)
{
int i;
double T = 3.14159265358979 * 10;
double dt = T / N;
double xN_Euler,xN_Euler_new;
double yN_Euler,yN_Euler_new;
double xN,xN_new,yN,yN_new;
double kx1,kx2,kx3,kx4;
double ky1,ky2,ky3,ky4;
double xs,ys,err_Euler = 0,err_RK = 0;
xN = 1.0; xN_Euler = 1.0;
yN = 0.0; yN_Euler = 0.0;
for (i = 0; i < N; i++)
{
//Exact Solution
xs = cos(i * dt);
ys = sin(i * dt);
//RK
err_RK = err_RK + dt * (fabs(xN - xs) + fabs(yN - ys));//Calc Sa
kx1 = -yN;
ky1 = xN;
kx2 = -(yN + dt / 2 * ky1);
ky2 = xN + dt / 2 * kx1;
kx3 = -(yN + dt / 2 * ky2);
ky3 = xN + dt / 2 * kx2;
kx4 = -(yN + dt * ky3);
ky4 = xN + dt * kx3;
xN_new = xN + dt * (kx1 + 2 * kx2 + 2 * kx3 + kx4) / 6;
yN_new = yN + dt * (ky1 + 2 * ky2 + 2 * ky3 + ky4) / 6;
xN = xN_new;
yN = yN_new;
//Euler
err_Euler = err_Euler + dt * (fabs(xN_Euler - xs) + fabs(yN_Euler - ys));
xN_Euler_new = xN_Euler + dt * (-yN_Euler);
yN_Euler_new = yN_Euler + dt * (xN_Euler);
xN_Euler = xN_Euler_new;
yN_Euler = yN_Euler_new;
}
err_RK = err_RK / T;
err_Euler = err_Euler / T;
printf("err_RK
= %15.15f,dt = % 15.15f\n",err_RK,dt);
printf("err_Euler = %15.15f,dt = % 15.15f\n",err_Euler,dt);
return 0;
}
main6.c
14年11月5日水曜日
厳密解を計算
RK法にて差を計算
オイラー法にて差を計算
常微分方程式の数値解法
オイラー法 V.S. 4段4次ルンゲ・クッタ法 （誤差対決）
Exercise!
main6-1.c
main6.cにおいて
1
155.434715579091460
8720.135117209718373
14年11月5日水曜日
155.434715579091460
常微分方程式の数値解法
オイラー法 V.S. 4段4次ルンゲ・クッタ法 （誤差対決）
Exercise!
main6.cにおいて
とし，誤差がどのようになるか調べ，下記の表を完成させよ
155.434715579091460
x10
x10
x10
14年11月5日水曜日
8720.135117209718373
155.434715579091460
27.716996517146850
0.001606408974007
0.372568841218112
0.000000162268631
0.031937429963379
0.000000000016235
0.003146750675461
0.000000000000019
常微分方程式の数値解法
バネ・ダンパ系の運動方程式
壁
玉
床
14年11月5日水曜日
壁
常微分方程式の数値解法
バネ・ダンパ系の運動方程式
玉
壁
壁
床
玉はバネの伸び に比例して力を受ける．（比例定数は ）
玉はダンパから速度 に比例して抵抗力を受ける．（比例定数は ）
14年11月5日水曜日
常微分方程式の数値解法
バネ・ダンパ系の運動方程式
玉
壁
壁
床
2
d x
dx
+
2
+
dt2
dt
14年11月5日水曜日
2
0x
=0
常微分方程式の数値解法
バネ・ダンパ系の運動方程式
Exercise! 以下のバネ・ダンパ系の運動方程式をRK法を用いて解き，
main7.c
可視化せよ．
2
d x
dx
+
2
+
2
dt
dt
2
0x
=0
← そっと手を離すという条件
← 時刻0でちょっと右にあるという条件
（条件さえ満たせば何でも良いが，例えば
など）
14年11月5日水曜日
レポート
Aコース. 授業で扱ったバネ・ダンパ系の運動方程式をRK法を解く
プログラムにおいて，厳密解と近似解を比較せよ．
ヒント：厳密解は導出することもできるが，図書館（or参考書）な
どで調べることもできる．
O
Bコース. 2重振り子の運動方程式をRK法を
用いて可視化せよ．（初期条件などは適宜決
l0
定すること）また厳密解を明示的に与えるこ
シミュレーション
とができるか考察せよ．
✓000 =
✓100
=
g(2m0 + m1 ) sin ✓0 + gm1 sin(✓0 2✓1 ) + 2m1 sin(✓0 ✓1 ) l0 ✓002 cos(✓0
2l0 (m0 m1 cos2 (✓0 ✓1 ) + m1 )
sin(✓0
✓1 ) (m0 + m1 ) g cos ✓0 + l0 ✓002 + l1 m1 ✓102 cos(✓0
l1 (m0 m1 cos2 (✓0 ✓1 ) + m1 )
✓1 ) + l1 ✓102
l1
✓1 )
m1
14年11月5日水曜日
m0
✓1
✓0
© Copyright 2024 ExpyDoc
