数値解析
2023年度前期 第14週 [7月20日]
静岡大学大学院
工学研究科機械工学専攻
ロボット・計測情報分野
創造科学技術大学院
情報科学専攻三浦 憲二郎講義アウトライン [7月20日]•常微分方程式•1段階法と多段階法•オイラー法•ホイン法•ルンゲ・クッタ法

微分方程式の重要性：復習•自然現象や工学モデルの記述•微分方程式が用いられる．
•常微分方程式と偏微分方程式
理科系の数学入門 (斎藤秀司・三松佳彦・戸瀬信之編)
第4巻 常微分方程式 高崎金久著，日本評論社

常微分方程式の起源は物理学，特に力学の研究にある．動力学の基本法則は 17世
紀にニュートンによって見出された．それは物体の加速度と物体に働く力の比例関係を表すものであり，物体の座標を時間の関数と見なすとき，２階までの導関数を含む常微分方程式にになる．さらに時代が下って， 19世紀に電気・磁気・熱現象が研究対象になると，時間と空間の両方に関する物理量の変化の法則を書き下すために偏微分方
程式も用いられるようになった．
今日では，理工系の学問を学ぶ人にとって（程度の差はあっても）微分方程式に関する知識は必要不可欠である．さらに，近年は経済学においても微分方程式の重要性が高まっている．もちろん，純粋に数学的な観点から見ても，微分方程式は重要かつ興味深いものであり，解析学の中で中心的な位置を占めている．

•簡単な問題を除き、解析解を求めるのは困難->
数値解法が必要不可欠

1段法と多段法 p.154•f(x,y)は a≦ｘ≦b, y∊ℜで定義された2変数関数とする．1階の常微分方程式y’=f(x,y) を考える．このとき，任意の x∊ [a,b]に対して y’=f(x,y(x))となるようなC1級の関数y=y(x)を常微分方程式の一般解または解という．y0が与えられたとき，y(a)=y0 を満たすような解を求める問題を常微分方程式の初期値問題という．

•数値解法
•１段階法
•多段階法

1段階法と多段階法 p.154•区間[a,b]を等間隔 a=x0<x1<x2<…<xn=b, h=(b-a)/n に分割
φ(x,y) は与えられた関数，関数値y(xi)の数値解をYiとする．
•１段階法

•多段階法

定数αj(j=0,1,…,k-1)

右辺にYi+1が含まれている公式：陰公式
いない公式：陽公式

1段階法と多段階法 p.155•１段階法•次数pの精度を持つ，またはp次の精度をもつという．
•局所離散化誤差•多段階法•次数pの精度を持つ，またはp次の精度をもつという．
•局所離散化誤差

オイラー法 p.155•[a,b]を n 等分し, 各分点 xi を
h : 刻み幅
yi = y(xi) とすると，テイラーの定理より•オイラー法 （1次の精度を持つ1段階法）

オイラー法：プログラム p.157#include <stdio.h>

実行結果
分割数を入力してください--- > 10x=0.1 y=1.1
x=0.2 y=1.22
x=0.3 y=1.362
x=0.4 y=1.5282
x=0.5 y=1.72102
x=0.6 y=1.943122
x=0.7 y=2.197434
x=0.8 y=2.487178
x=0.9 y=2.815895
x=1 y=3.187485/* 関数の定義 */double func(double x, double y);/* オイラー法 */void euler(double x, double y, double a, double b, int n,double (*f)(double, double) );int main(void){
int n;
printf("分割数を入力してください--->");
scanf("%d",&n);euler( 0.0, 1.0, 0.0, 1.0, n, func );return 0;}/* オイラー法 */void euler(double x, double y, double a, double b, int n,double (*f)(double, double) ){double h;int i;}

h = (b-a)/n; /* 刻み幅 */for ( i = 0 ; i < n ; i++) {y = y + h * (*f)(x,y);x += h;printf("x=%f ¥t y=%f ¥n", x, y );}/* 関数の定義 */double func(double x, double y){return( x + y );
}

3.5
3
2.5
2
系列1

1.5
1
0.5
0

0

0.2

0.4

0.6

0.8

11.2

ホイン法 p.158•式(8.2）より
2変数のテイラーの定理より

1変数のテイラーの定理よりホイン法 p.158•次数2の1段階法

ホイン法 （Heun method）•ホイン法アルゴリズム

ホイン法：プログラム p.159#include <stdio.h>

実行結果
分割数を入力してください--- > 10
x=0.1
y=1.11x=0.2
y=1.24205
x=0.3
y=1.398465
x=0.4
y=1.581804
x=0.5
y=1.794894
x=0.6
y=2.040857
x=0.7
y=2.323147
x=0.8
y=2.645578
x=0.9
y=3.012364
x=1
y=3.428162/* 関数の定義 */double func(double x, double y);/* ホイン法 */void heun(double x, double y, double a, double b, int n,double (*f)(double, double) );int main(void){
int n;
printf("分割数を入力してください--->");
scanf("%d",&n);heun( 0.0, 1.0, 0.0, 1.0, n, func );return 0;}/* ホイン法 */void heun(double x, double y, double a, double b, int n,double (*f)(double, double) ){double k1, k2, h;
int i;
h = (b-a)/n;for ( i = 0 ; i < n ; i++){k1 = f(x,y); k2 = f(x+h, y+h*k1);y = y + h/2.0 * ( k1 + k2 );x += h;printf("x=%f ¥t y=%f ¥n", x, y );}}/* 関数の定義 */double func(double x, double y){return( x + y );
}

4
3.5
32.5
2

系列1

1.5
1
0.5
00

0.2

0.4

0.6

0.8

11.2

ルンゲ・クッタ法 p.160•1段階法でもっともよく使われる． 次数4のルンゲ・クッタ法(Runge-Kutta)ルンゲ・クッタ法：プログラム p.161#include <stdio.h>#include <stdlib.h>

実行結果
分割数を入力してください--- > 10
x=0
y=1
x=0.1
y=1.110342
x=0.2
y=1.242805
x=0.3
y=1.399717
x=0.4
y=1.583648
x=0.5
y=1.797441
x=0.6
y=2.044236
x=0.7
y=2.327503
x=0.8
y=2.651079
x=0.9
y=3.019203
x=1
y=3.436559

double *dvector(long i, long j); /* ベクトル領域の確保 */void free_dvector(double *a, long i); /* 領域の解放 */double func(double x, double y); /* 関数の定義 *//* ルンゲ・クッタ法 */double *rk4( double y0, double *y, double a, double b, int n,double (*f)(double, double) );int main(void){
double *y, h, a=0.0, b=1.0, y0=1.0 ;int i, n;printf("分割数を入力してください--->");scanf("%d",&n);y = dvector( 0, n );/* 領域の確保 */y = rk4( y0, y, a, b, n, func ); /* ルンゲ・クッタ法 *//* 結果の表示 */h = (b-a)/n ;/* 刻み幅 */for ( i = 0 ; i <= n ; i++) {printf("x=%f ¥t y=%f ¥n", a+i*h, y[i] );}free_dvector( y, 0 ); /* 領域の解放 */return 0;}/* ルンゲ・クッタ法 */double *rk4( double y0, double *y, double a, double b, int n,double (*f)(double, double) ){double k1, k2, k3, k4, h, x;
int i;
h = (b-a)/n;/* 初期値の設定 */y[0] = y0; x = a;/* ルンゲ・クッタ法 */for ( i = 0 ; i < n ; i++) {k1 = f(x,y[i]); k2 = f(x+h/2.0, y[i]+h*k1/2.0);
k3 = f(x+h/2.0, y[i]+h*k2/2.0);
k4 = f(x+h, y[i]+h*k3);
y[i+1] = y[i] + h/6.0 * ( k1 + 2.0*k2 + 2.0*k3 + k4 );x += h;}}return y;

4
3.5
32.5
2

系列1

1.5
10.5
00

0.2

0.4

0.6

0.811.2

ルンゲ・クッタ型•微分積分学入門第二課，一松 信，近代科学社 p.176•注意：教科書 p.162 のαi -> βi，ki -> ki, λi -> ωi,β0 -> α11, β1 ->α21, γ1 -> α22, β2 -> α31, γ2 -> α32, δ2 -> α33•微分方程式の初期値問題 y’=f(x,y), y(x0)=y0•ynがわかったとき，xn+1=xn+h での yn+1の値を求める．
•k1, k2, … , kp を補助変数としてもし条件 i≦j のとき αij=0 が成り立てば，

ルンゲ・クッタ型公式の導出•微分積分学入門第二課，一松 信，近代科学社 p.182•f(x,y)が十分に滑らかとすると，f’(x,y)の高階導関数は，合成関数の微分の公式により

fx, fy はx, y に関する偏微分ルンゲ・クッタ型公式の導出•基準点(x,y)=(xn,yn)においてテイラー展開して3次の項までとる．

•h, h2, h3, h4 の項からなる．kj に関する式を代入し，全体としてh5以上の項になる部分を捨てて整理する．

ルンゲ・クッタ型公式の導出•∑ωi ki を真の解のテイラー展開の1次以上４次以下の項

ルンゲ・クッタ型公式の導出

ルンゲ・クッタ型公式の導出
段 : fを計算する回数
位 : 一致する次数
1段1位オイラー法

2段2位ホイン法

4段4位ルンゲ・クッタ法

まとめ•常微分方程式•1段階法と多段階法•オイラー法•ホイン方•ルンゲ・クッタ法

T.Kanai, U.Tokyo