ルンゲクッタ法の説明と実装例[C言語]

3日かけてルンゲクッタ法を理解する所までにきました．

とは言うものの，この記事では1階常微分方程式についての4次ルンゲクッタについて説明しています．

前回，前々回の記事である

・オイラー法と修正オイラー法の説明と実装例
・ニュートン・コーツの公式とシンプソンの公式の導出と実装

を読むと今回の記事で説明するルンゲクッタ法の導出部分が理解しやすいかと思います．

1.(4次)ルンゲクッタ法の簡単な説明

詳しい説明はしたところであまり興味がない人が多い気がするので，詳しい説明は最後にすっ飛ばして，実装するにおいて最低限の部分だけ示しておきます．

次の初期値がわかっている1階常微分方程式を解く事を考えます．
$$\begin{eqnarray*} f(x_0) &=& y_0 \\ \frac{df(x)}{dx} &=& g(x,f(x)) \end{eqnarray*}$$ この時，4次ルンゲクッタ法では次のようにして計算します． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta x$$ ただしk1,...,k4は $$\begin{eqnarray*} k_1 &=& g(x_k,f(x_k)) \\ k_2 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_1 \frac{\Delta x}{2}) \\ k_3 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_2 \frac{\Delta x}{2}) \\ k_4 &=& g(x_{k+1},f(x_k)+k_3 \Delta x) \end{eqnarray*}$$ として計算しています．
この辺の導出は一番最後でしているので気になる人は見てってくださいな．

2.実装例

何を実装しようか悩んで適当にWikipediaを漁っていたらロジスティック方程式というモデルがあったので，これのシミュレーションをしてみたいと思います．

Wikipediaによれば，生物の個体数がどのように変化していくかを説明するモデルの一つであるとのこと．

求める常微分方程式はこれ． $$\frac{dN(t)}{dt} = r \frac{K-N(t)}{K} N(t)$$ rは一匹の個体の増加率，N(t)は時間tにおける生物の個体数，Kは個体数のしきい値で，最終的に安定するNの値．

例えば水槽なら最大でも20匹ぐらいしか入らない大きさならK=20みたいな値になると思う．

ちなみに厳密解は普通の変数分離法で多少計算が面倒だけど求まる．初期値をN_0とすると次のようになる． $$N(t) = \frac{N_0 K e^{rt}}{K - N_0 + N_0 e^{rt}}$$ 厳密解はさておき，この常微分方程式にルンゲクッタ法を適用する．

まず見通しを良くするために，1で示した式のg(x,f(x))に当たる部分g(t,N(t))を定義しておく． $$\begin{eqnarray*} \frac{dN(t)}{dt} &=& r \frac{K-N(t)}{K} N(t) \\ &=& g(t,N(t)) \end{eqnarray*}$$ そうすると$\Delta t$を微小な値として，$t_k = k \Delta t$とすると，漸化式はこのようになる． $$\int_{t_k}^{t_{k+1}} g(t,N(t)) dt \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta t$$ ただしk1,...,k4は次の式で計算する． $$\begin{eqnarray*} k_1 &=& g(t_k,N(t_k)) \\ k_2 &=& g(t_{k+\frac{1}{2}},N(t_k) + k_1 \frac{\Delta t}{2}) \\ k_3 &=& g(t_{k+\frac{1}{2}},N(t_k) + k_2 \frac{\Delta t}{2}) \\ k_4 &=& g(t_{k+1},N(t_k)+k_3 \Delta t) \end{eqnarray*}$$ それではC言語による実装．関数g,を定義することで，数式とほぼ変わらないプログラムが書けます．実装の係数は適当なのでアレですが．

	/* 4次ルンゲクッタ法によるロジスティック方程式の解の計算 */
	#include <stdio.h>
	#define K 1000.0 // 個体数のしきい値
	#define r 0.9 // 1個体の増加率
	#define N_LOOP 10000 // ループ回数
	// 積分する関数g(t,N)
	double g(double t,double N_t){
	return r*(K-N_t)*N_t/K;
	}
	int main(void){
	int k; // ループカウンタ
	double dt = 1e-3; // 刻み幅Δt
	double N = 1.0; // 個体数t（初期値を代入しておく)
	double k1,k2,k3,k4; // 計算に使う諸変数
	double a,m,b; // t_k,t_k+0.5,t_k+1を格納する
	for( k = 0 ; k < N_LOOP ; k++ ){
	// t_k,t_k+0.5,t_k+1を計算
	a = k*dt;
	m = (k+0.5)*dt;
	b = (k+1.0)*dt;
	// k1,...,k4の値を計算
	k1 = g(a,N);
	k2 = g(m,N+k1*dt*0.5);
	k3 = g(m,N+k2*dt*0.5);
	k4 = g(b,N+k3*dt);
	// 積分して代入
	N = N + (dt/6.0)*(k1+2.0*k2+2.0*k3+k4);
	// 値を表示
	printf("%lf,%e\n",b,N);
	}
	return 0;
	}

view raw RK4_KogisticEquation.c hosted with ❤ by GitHub

シミュレーション結果．

このシミュレーションだとほとんど違いがわかりませんね．誤差をプロットするとこうなります．

オーダーが1e-5といったところでしょうか．誤差が振動しているのがわかります．シミュレーションの精度としては十分なのでしょう．

4次のルンゲクッタ法では4次の項まで級数展開が一致しているとかで，誤差は$O(\Delta t^5)$に比例して小さくなるとのことです．

流石にルンゲクッタ法でシミュレーションすれば精度どうなの？とか言われても大丈夫そう……？まだその辺はよく知らないけれど．

3.ルンゲクッタ法の説明

手元の本には書かれていないのではっきりしたことはわからないのですが軽く調べたところによると，修正オイラー法は2次ルンゲクッタ法とも呼ぶらしいです．

なので2次ルンゲクッタについては前回の記事を参照してください．ここでは4次ルンゲクッタについて説明します．

今回はとりあえず1階微分方程式についての4次ルンゲクッタについて説明するので，次の常微分方程式について考えます．
$$f(x_0) = y_0 \\ \frac{df(x)}{dx} = g(x,f(x))$$ ここで　$x_k = x_0 + k \Delta x$　と置くと　$f(x_{k+1})$　と　$f(x_{k})$　の漸化式を得ます． $$\begin{eqnarray*} \int_{x_k}^{x_{k+1}} \frac{df(x)}{dx} dx &=& f(x_{k+1})-f(x_k) \\ &=& \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \\ \\ ∴f(x_{k+1}) &=& f(x_k) + \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \end{eqnarray*}$$

この $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分をどう計算するかによってオイラー法，修正オイラー法，ルンゲクッタ法が分類できます．

オイラー法では矩形積分，修正オイラー法では台形積分，ルンゲクッタ法ではシンプソンの公式を使います．

シンプソンの公式は次のようになります． $$\int_{x_k}^{x_{k+1}}h(x) dx \approx \frac{1}{6}(h(x_k)+4h(x_k+\frac{\Delta x}{2}) + h(x_{k+1}) )\Delta x$$ これを使って積分値を計算すれば良いわけです．$x_{k+\frac{1}{2}} = x_k+\frac{1}{2}\Delta x \\ $ として $$\begin{eqnarray*} g_k &=& g(x_k,f(x_k)) \\ g_{k+\frac{1}{2}} &=& g(x_{k+\frac{1}{2}},f(x_{k+\frac{1}{2}})) \\ g_{k+1} &=& g(x_{k+1},f(x_{k+1})) \\ \end{eqnarray*}$$ という値を考えると， $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ は次のように近似できます． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(g_k + 4g_{k+\frac{1}{2}} + g_{k+1} )\Delta x$$
しかしこの式には未知の値$f(x_{k+\frac{1}{2}})$や$f(x_{k+1})$を使っています．これらの値を矩形積分で求めて計算することもできますが，良い精度が得られません．

ルンゲクッタ法ではこれらの値を二段階の計算で計算します．
$$\begin{eqnarray*} k_1 &=& g(x_k,f(x_k)) \\ k_2 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_1 \frac{\Delta x}{2}) \\ k_3 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_2 \frac{\Delta x}{2}) \\ k_4 &=& g(x_{k+1},f(x_k)+k_3 \Delta x) \end{eqnarray*}$$ $g_{k+\frac{1}{2}}$に当たる部分をk2,k3という形で計算しています．これらを先ほどのシンプソンの公式に代入することで以下のような近似を行うものを4次のルンゲクッタ法と呼ぶそうです． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta x$$

今日はここまで．

ニュートン・コーツの公式とシンプソンの公式の導出と実装[C言語]

ルンゲクッタ法はどうやって導かれたのか謎だったのですが，『情報処理入門コース　数値計算』によれば，シンプソンの公式から来ているとのこと．
なのでシンプソンの公式がどういう経緯でできたのかメモしていきます．

※数式がうまく表示されない時はPCのブラウザから見ることを推奨します

1.ニュートン・コーツの積分公式

まずラグランジュ補間の公式が与えられているものとして使います．f(x)が近似したい式，p(x)が近似式．x_0,...,x_nが関数f(x)上のわかっている点． $$p(x)=\sum_{i=0}^{n}f(x_i)N(i,x) \\ N(i,x) = \prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j}$$ ラグランジュ補間の出処はよくわかんないけど関数N(x)がクロネッカーのデルタ関数のようになっていることから意味をお察しください．とりあえずラグランジュ補間については前書いたのでそっちも見てくださいな．

次に，このラグランジュ補間で求めた近似式を用いて関数f(x)を区間[a,b]で積分する事を考えます． $$a = x_0,x_1,x_2,...,x_n = b$$ という様に区間をn分割し，上記のラグランジュ補間を適用すると次のようになります． $$\begin{eqnarray*} S &=& \int_{a}^{b}f(x)dx \\ & \approx & \int_{a}^{b}p(x)dx \\ &=& \int_{a}^{b} \sum_{i=0}^{n}f(x_i)N(i,x) dx \\ \end{eqnarray*}$$ ここで少しわかりづらいのですが，x_iは定数なのでf(x_i)はxの関数では無いことから， $$\begin{eqnarray*} S & \approx & \int_{a}^{b} \sum_{i=0}^{n}f(x_i)N(i,x) dx \\ &=& \sum_{i=0}^{n}f(x_i) \int_{a}^{b}N(i,x) dx \end{eqnarray*}$$ という式が得られます．ここで関数N(i,x)はx_0,...,x_nが[a,b]をn等分した値，すなわち $$x_k = a + k\Delta x \\ \Delta x = \frac{b-a}{n}$$ と表せることを用いて，次の変数変換をすることを考えます． $$t = \frac{x-a}{\Delta x}\\ dt = \Delta x \, dx$$ ここで，N(i,x)は上記のtを代入して次のように変形したものをM(i,t)とする． $$\begin{eqnarray*} N(i,x) &=& \prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j} \\ &=& \prod_{j=0,j \neq i}^{n} \frac{(a+ t \Delta x) - (a + j \Delta x)}{(a+i \Delta x) - (a + j \Delta x ) } \\ &=& \prod_{j=0,j \neq i}^{n} \frac{t-j}{i-j} \\ &=& M(i,t) \end{eqnarray*}$$ このN(i,x)を変形した関数M(i,t)を用いると先ほどの積分はこうなる． $$\begin{eqnarray*} S & \approx & \sum_{i=0}^{n}f(x_i) \int_{a}^{b}N(i,x) dx \\ &=& \sum_{i=0}^{n}f(x_i) \Delta x \int_{0}^{n}M(i,t) dt \end{eqnarray*}$$ M(i,t)は変数xや，区間[a,b]などに依存しないから， $$C_i = \int_{0}^{n}M(i,t) dt$$ と置いて，先に手計算で計算しておけば使い回しが効いて，最終的に次のようになります． $$S \approx \sum_{i=0}^{n}f(x_i) C_i \Delta x$$ この公式をニュートン・コーツの積分公式と呼ぶとのこと．

2.シンプソンの公式

うえで求めたニュートン・コーツの積分公式にn=2を代入したものが世間一般にシンプソンの公式と呼ばれているものとのこと．C_iの値は上式に突っ込んで $$\begin{eqnarray*} C_0 &=& \int_0^2 \frac{t-1}{-1} \frac{t-2}{-2} dt &=& \frac{1}{3} \\ C_1 &=& \int_0^2 \frac{t}{1} \frac{t-2}{-1} dt &=& \frac{4}{3} \\ C_2 &=& \int_0^2 \frac{t}{2} \frac{t-1}{1} dt &=& \frac{1}{3} \\ \end{eqnarray*}$$ となる．この公式のすごいところはこの係数をあらかじめ求めておけば面倒な計算を端折ることができるところだと思います．
この値をニュートン・コーツの積分公式に代入すると， $$\begin{eqnarray*} \int_{a}^{b} f(x) dx &=& \sum_{i=0}^{2}f(x_i) C_i \Delta x\\ &=& (\frac{f(x_0)}{3} +\frac{4f(x_1)}{3} +\frac{f(x_2)}{3}) \Delta x\\ &=& (\frac{f(x_0)}{3} +\frac{4f(x_1)}{3} +\frac{f(x_2)}{3}) {\frac{b-a}{2}}\\ &=& \frac{1}{6}(b-a)(f(x_0)+4f(x_1) +f(x_2)) \end{eqnarray*}$$ となります．

シンプソンの公式を図で表すとこんな感じ．青色の線が二次式で近似した曲線．その曲線の積分を求めてる．

ちなみにn=3の場合などはニュートン・コーツの公式のWikipediaに載ってたりするのでそちらを参照すると良いかと思います．

3.プログラム

抽象的な話はこの辺にして，実際にプログラムを書いてみます．試しに$f(x) = \frac{4}{x^2+1}$を[0,1]区間で積分してみます．

ちなみに$\int_{0}^{1}f(x)dx$の厳密解はπです．

まず $\Delta x = \frac{1}{n}$ とし，$x_k = k \Delta x$ と置くと， $$\begin{eqnarray*} \int_{0}^{1}f(x)dx = \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}}f(x) dx \end{eqnarray*}$$ という微小区間 $[x_k,x_{k+1}]$ の積分となります．この微小な区間に対してシンプソンの公式を適用すると次のようになります． $$\int_{x_k}^{x_{k+1}}f(x) dx \approx \frac{1}{6}(f(x_k)+4f(x_k+\frac{\Delta x}{2}) + f(x_{k+1}) )\Delta x$$ 実際はf(x_k)の足し算がかぶっているのでもう少しスマートにできるのですが，今回は愚直にそのままこれを実装します．

	/* シンプソンの公式による1/(x^2+1)の計算 */
	#include <stdio.h>
	// 被積分関数
	double f(double x){
	return 4.0/(x*x+1.0);
	}
	int main(void){
	int k; // ループカウンタ
	int N = 1000; // 分割数
	double dx = 1.0/N; // 刻み幅Δx
	double S = 0.0; // 和を格納する変数
	double a; // x_kの値を格納する
	double b; // x_{k+1}の値を格納する
	double m; // a,bの中間の値を格納する
	for( k = 0 ; k < N ; k++ ){
	// 近似を行う3点を計算
	a = k * dx;
	m = a + dx / 2.0;
	b = (k+1) * dx;
	// f_k + g(x_k,f(x_k))をfに代入
	S += (1.0/6.0)*(f(a)+4.0*f(m)+f(b))*dx;
	// 値の表示
	printf("%lf,%e\n",b,S);
	}
	return 0;
	}

view raw SimpsonsRule.c hosted with ❤ by GitHub

実行結果

結果が3.141593なので円周率が概ねちゃんと求まっているので大丈夫そうです．

オイラー法と修正オイラー法の説明と実装例[C言語]

学科的に流体シミュレーションとか多分そのうちやると思うのですが，ルンゲクッタ法を理解していなかったのでとりあえず常微分方程式の数値計算の基本であるオイラー法と修正オイラー法について記事を書いてみます．

1.オイラー法

オイラー法は精度は悪いですが，非常に簡単で計算が高速です．研究などでは恐らく使われることはありませんが，ゲームプログラミング等をしているとシューティングの弾の運動などは大抵オイラー法が使われています．
今回は説明のための例として，次のような常微分方程式があったとします．今回は初期値f(x_0)がわかっているものとします． $$f(x_0) = y_0 \\ \frac{df(x)}{dx} = g(x,f(x))$$ ここで　$x_k = x_0 + k \Delta x$　と置くと　$f(x_{k})$　の値は次のようにして求められることがわかります． $$\begin{eqnarray*} f(x_{k}) &=& f(x_0) + \int_{x_0}^{x_k} g(x,f(x)) dx \\ &=& f(x_0) + \sum_{i=0}^{k-1} \int_{x_{i}}^{x_{i+1}} g(x,f(x)) dx \end{eqnarray*}$$ ここで　$f(x_{k+1})$　と　$f(x_{k})$　の漸化式を考えます． $$\begin{eqnarray*} \int_{x_k}^{x_{k+1}} \frac{df(x)}{dx} dx &=& f(x_{k+1})-f(x_k) \\ &=& \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \\ \\ ∴f(x_{k+1}) &=& f(x_k) + \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \end{eqnarray*}$$ この $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分がキモです．この部分をどうするかで精度が大きく変わってきます．

オイラー法ではg(x,f(x))が区間$[x_k,x_{k+1}]$でほとんど変化しないものとして，次のように近似しています． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx g(x_k,f(x_k)) \Delta x$$ これを漸化式に代入したものがオイラー法です． $$f(x_{k+1}) \approx f(x_k) + g(x_k,f(x_k)) \Delta x$$

2.修正オイラー法

オイラー法は非常に簡単なのですが精度がよろしくないです．なのでもう少し精度を上げる必要があります．

修正オイラー法はその改良版の一つです．もっとも，これもあまり制度が良いというわけではないですが．

先ほど， $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分がキモと言いました．オイラー法では関数g(x,f(x))があまり変化しないものとして計算していました．イメージとしてはこんな感じです．

矩形近似とか長方形近似とか呼ばれているタイプの積分方法です．見るからに精度がよろしくなさそうなのがわかります．

そこで台形公式というものを使って修正オイラー法というものが提案されました．イメージとしてはこんな感じです．

名前の通り台形で近似するというものです．数式で表すと次のようになります．台形の面積は((高さ)×（上辺）+(下辺))÷2という小学生の時に習った公式を用いて，

$$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx (g(x_k,f(x_k)) + g(x_{k+1},f(x_{k+1}))) \frac{\Delta x}{2}$$ という近似式で計算します．しかしこれを漸化式に代入すると $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},f(x_{k+1}))) \frac{\Delta x}{2}$$ となり，$f(x_{k+1})$ を求めるために $f(x_{k+1})$ を使っているのでこれでは使えません．なので右辺の $f(x_{k+1})$ だけは，オイラー法で使った矩形近似を用いて次のように求めておきます． $$\widehat{f(x_{k+1})} = f(x_k) + g(x_k,f(x_k)) \Delta x$$ この　$\widehat{f(x_{k+1})}$　を用いて漸化式に代入したものが修正オイラー法です． $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},\widehat{f(x_{k+1})})) \frac{\Delta x}{2}$$

3.利用例

ここまで抽象的すぎたのでここで一度，次の1階微分方程式を数値的に解いてみたいと思います． $$\frac{df(x)}{dx} - 2f(x) = x$$ まぁ適当に作った式なので意味は無いです．ちなみに厳密解はCを定数として $$f(x) = -\frac{1}{4}(2x+1) + Ce^{2x}$$ と表されます．ちなみに定数Cは関数f(x)の初期値f(0) = aとした時，$C = a + \frac{1}{4}$になります．
今回は初期値 f(0) = a = 1 として数値的に求めてみます．

上記の微分方程式を変形してみます． $$\begin{eqnarray*} \frac{df(x)}{dx} &=& x + 2f(x) \\ &=& g(x,f(x)) \end{eqnarray*}$$ これでg(x,f(x))が求まりましたので，先ほど示した漸化式に代入してみれば，オイラー法では $$f(x_{k+1}) \approx f(x_k) + (x_k + 2f(x_k)) \Delta x$$ となります.

修正オイラー法では，まず　$\widehat{f(x_{k+1})}$　が $$\begin{eqnarray*} \widehat{f(x_{k+1})} &=& f(x_k) + g(x_k,f(x_k)) \Delta x \\ &=& f(x_k) + (x_k + 2f(x_k)) \Delta x \\ \end{eqnarray*}$$ となることから，次の漸化式に代入してやれば良いです． $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},\widehat{f(x_{k+1})})) \frac{\Delta x}{2} \\ $$ 全部数式を代入すれば,x_kとf(x_k)とΔxだけの式にできますが，数式が長くなるので省略．

もっとも，関数をプログラム内で定義してやればそんなに気にならないと思います．

4.精度の評価

ソースコードは長くなりますので先にオイラー法，修正オイラー法の精度を示しときます．

もっとも，評価と言ってもこの例での真値との差を見るだけですが．以下数値的に求めたものと，真値との差です．
刻み幅Δx=0.001，反復回数1000回です．euler's method 2の方が修正オイラー法です．適切な英語が見つからなかったので……．

1000回繰り返して差が10^-7オーダーですね．とりあえずそれらしい値は取れているのがわかります．

これを見る限り修正オイラー法はオイラー法と比べたら優秀そうではあります．
ただ傾きが指数関数的に増えていきそうなので後日ルンゲクッタ法を実装してそれと比較したいです．

5.ソースコード

C言語による実装．オイラー法はこんな感じです．

	/* Euler法による微分方程式 f' - 2f = x の計算 */
	#include <stdio.h>
	#include <math.h>
	double g(double x,double fx){
	return x + 2.0 * fx;
	}
	int main(void){
	int k; // ループカウンタ
	int N = 10000; // 繰り返し数
	double dx = 1e-3; // 刻み幅Δx
	double f = 1.0; // 関数f(x)の値(初期値)
	double x_0 = 0.0; // 初期値の時のxの値
	double x_k; // x_kの値
	for( k = 1 ; k <= N ; k++ ){
	x_k = x_0 + k * dx; // x_kの値を求める
	f = f + g(x_k,f) * dx; // f_k + g(x_k,f(x_k))をfに代入
	// 値の表示
	printf("%lf,%e\n",x_k,f);
	}
	return 0;
	}

view raw EulersMethod.c hosted with ❤ by GitHub

修正オイラー法はこんな感じです．

	/* 修正Euler法による微分方程式 f' - 2f = x の計算 */
	#include <stdio.h>
	#include <math.h>
	double g(double x,double fx){
	return x + 2.0 * fx;
	}
	int main(void){
	int k; // ループカウンタ
	int N = 10000; // 繰り返し数
	double dx = 1e-3; // 刻み幅Δx
	double f = 1.0; // 関数f(x)の値(初期値)
	double f_bar; // f(x_k+1)を矩形積分でとりあえず求めた値
	double x_0 = 0.0; // 初期値の時のxの値
	double x_k; // x_kの値
	for( k = 1 ; k <= N ; k++ ){
	x_k = x_0 + k * dx; // x_kの値を求める
	// まず矩形積分でf_k+1を見積もって
	f_bar = f + g(x_k,f) * dx;
	// 台形公式で積分する
	f = f + (g(x_k,f) + g(x_k + dx, f_bar)) * dx /2.0;
	// 値の表示
	printf("%lf,%e\n",x_k,f);
	}
	return 0;
	}

view raw EulersMethod2.c hosted with ❤ by GitHub

今日はここまで．