２０２３年度前期 数値解析 第１４回講義ノート
２０２３年７月２０日 三浦憲二郎
スライド１ 講義アウトライン•

常微分方程式•

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.155[a,b]を n 等分し, 各分点 xi をh: 「」幅「yi = y(xi) とすると，•」をもつという．

」誤差

スライド５ オイラー法•」を持つ，または p 次の「オイラー法 （「」の定理より」次の精度を持つ「」段階法）

スライド６ オイラー法：プログラム p.157#include <stdio.h>/* 関数の定義 */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;

実行結果
分割数を入力してください--- > 10
x=0.1

y=1.1

x=0.2

y=1.22x=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

}/* オイラー法 */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 );}スライド 7 ホイン法 p.158
次数「」の「」段階法

ホイン法 （Heun method）

ホイン法アルゴリズム

スライド８ ホイン法：プログラム p.159#include <stdio.h>/* 関数の定義 */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;}

実行結果
分割数を入力してください--- > 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

/* ホイン法 */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 );}スライド９ ルンゲ・クッタ法•p.1601 段階法でもっともよく使われる． 次数「」の「スライド１０ ルンゲ・クッタ法：プログラム
p.161#include <stdio.h>#include <stdlib.h>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++) {」法(Runge-Kutta)

実行結果
分割数を入力してください--- > 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

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;}
スライド１１ ルンゲ・クッタ型公式の導出•

微分積分学入門第二課，一松 信，近代科学社 p.182•

f(x,y)が十分に滑らかとすると，f’(x,y)の高階導関数は，合成関数の微分の公式により

fx, fy は x, y に関する偏微分•∑ωi ki を真の解のテイラー展開の 1 次以上４次以下の項

段 : f を計算する回数
位 : 一致する次数
1 段１位 オイラー法
２段１位 ホイン法

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