/*
 * 各種１段階法による微分方程式初期値問題の数値解析
 * 基本プログラム
 * コンパイル：cc diffeq.c -o diffeq -lm
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*
 * 問題：初期値問題 dy/dx = f(x,y) = y，y(0) = 1 を解け．
 */

// 初期値 y_0, Y_0
#define	Y0 	1.0

// 区間 [x_0, x_N]
#define X0	0.0
#define	XN	1.0

// 刻み幅 h
#define H	1.0
// #define H	1.0e-1
// #define H	1.0e-2
// #define H	1.0e-3

// 既知関数 f(x,y)
double f(double x, double y)
{
	return (y);
}

// 厳密解：y(x) = exp(x)
double exact_func(double x)
{
	return (exp(x));
}

// オイラー法の漸化式：Y_n+1 = Y_n + h f(x_n, Y_n)
// （１次の１段階法，１次のルンゲ-クッタ法）
double euler(double xn, double yn)
{
	return (yn + H*f(xn, yn));
}

/* 課題：実装せよ．
// ホイン法の漸化式：Y_n+1 = Y_n + ...
// （２次の１段階法，２次のルンゲ-クッタ法）
double heun(double xn, double yn)
{
	double	f1, f2;

	f1 = ...
	f2 = ...
	return (yn + (f1+f2)*H/2.0);
}
*/

/* 課題：実装せよ．
// ルンゲ-クッタ法の漸化式：Y_n+1 = Y_n + ...
// （４次の１段階法）
double rk4(double xn, double yn)
{
	...
}
*/

// 任意の１段階法による解析
// name：解析手法の名前の文字列
// approx_func：近似解の漸化式への関数ポインタ
void solver(char *name, double (*approx_func)(double x, double y))
{
	double	e;	// 絶対誤差ε
	double	yx;	// 厳密解 y
	double	ya;	// 近似解 Y_n, Y_n+1
	double	xn;	// x_n = x_0 + n h
	int	n = 0;

	double	delta = H/10.0;	// 微小値

	ya = Y0;
	printf("%s：\n", name);
	printf("n	x_n		近似解Y_n	厳密解y(x_n)	絶対誤差ε\n");

/* ［知識］
 * 理論的には一致するハズの２つの数値が，
 * 計算上の誤差によって，現実には一致しない場合がある．
 * 具体的に，for (...; xn <= XN; ...) では，
 * xn の最終値が XN に達しない場合がある．
 * 次の for 文では，微小値 delta の導入によって，
 * この誤差の影響を回避した．
 */
	for (xn = X0; xn <= XN + delta; xn += H) {
		yx = exact_func(xn);		// y(x_n)
		e = ya - yx;			// ε = Y_n - y(x_n)
		printf("%d	%f	%f	%f	%e\n", n, xn, ya, yx, e);
		ya = approx_func(xn, ya);	// Y_n+1 = Y_n + ...
		n++;
	}
	printf("\n");
}
/* ［知識］
 * 「関数ポインタ」とは，関数名を記憶するためのポインタ（変数）である．
 * （一方，「ポインタ関数」は，ポインタ（アドレス）を返す関数．互いに別物．）
 * 例えば今，この関数が solver("...", euler）として呼び出されている場合，
 * 関数ポインタは approx_func = euler; となり，
 * 上のコード ya = approx_func(xn, ya); は，
 * 実質的に，ya = euler(xn, ya); と同じことになる．
 */

int main()
{
	printf("初期値問題：dy/dx = y, y(0) = 1\n");
	printf("厳密解：y = e^x\n");
	printf("刻み幅：%e\n\n", H);

	solver("オイラー法", euler);
//	solver("ホイン法", heun);
//	solver("ルンゲ-クッタ法", rk4);
	return (0);
}