clear;

// *** 関数の定義を読み出し ***
exec('differential.sci',0);

// *** 計算の設定 ***
xmin = 0;
xmax = 2 * %pi;
x = linspace(xmin,xmax);
dy = cos(x);	// 解析解
N = [0:1:50]';	// 刻み幅

// *** 数値微分 ***
for i = 1:length(N)
  // 刻み幅
  n = N(i);
  dx = 1 / 2 ^ n;
  // 前進差分
  dyf1  = diff_f1(x, dx, sin);
  dyf12 = diff_f1(x, 2 * dx, sin);
  ERRmeanf1(i) = mean(abs(dy - dyf1));
  ERRmaxf1(i)  = max(abs(dy - dyf1));
  DIFmeanf1(i) = mean(abs(dyf12 - dyf1));
  DIFmaxf1(i)  = max(abs(dyf12 - dyf1));
  // 中心差分
  dyf2  = diff_f2(x, dx, sin);
  dyf22 = diff_f2(x, 2 * dx, sin);
  ERRmeanf2(i) = mean(abs(dy - dyf2));
  ERRmaxf2(i)  = max(abs(dy - dyf2));
  DIFmeanf2(i) = mean(abs(dyf22 - dyf2));
  DIFmaxf2(i)  = max(abs(dyf22 - dyf2));
  // 前進差分に対するRomberg1段
  dyf1r  = diff_f1r(x, dx, sin);
  dyf1r2 = diff_f1r(x, 2 * dx, sin);
  ERRmeanf1r(i) = mean(abs(dy - dyf1r));
  ERRmaxf1r(i)  = max(abs(dy - dyf1r));
  DIFmeanf1r(i) = mean(abs(dyf1r2 - dyf1r));
  DIFmaxf1r(i)  = max(abs(dyf1r2 - dyf1r));
  // 中心差分に対するRomberg1段
  dyf2r  = diff_f2r(x, dx, sin);
  dyf2r2 = diff_f2r(x, 2 * dx, sin);
  ERRmeanf2r(i) = mean(abs(dy - dyf2r));
  ERRmaxf2r(i)  = max(abs(dy - dyf2r));
  DIFmeanf2r(i) = mean(abs(dyf2r2 - dyf2r));
  DIFmaxf2r(i)  = max(abs(dyf2r2 - dyf2r));
end

// *** グラフのプロット ***
// *** 誤差のプロット ***
scf(0);
a = gca();
a.data_bounds = [min(N),1E-14; max(N),1];
a.log_flags = "nl";
// 誤差の平均値
plot(N, ERRmeanf1, '-or');	// 前進差分
plot(N, ERRmeanf2, '-om');	// 中心差分
plot(N, ERRmeanf1r, '-ob');	// 前進差分に対するRomberg1段
plot(N, ERRmeanf2r, '-og');	// 中心差分に対するRomberg1段
// 誤差の最大値
plot(N, ERRmaxf1, '-sr');	// 前進差分
plot(N, ERRmaxf2, '-sm');	// 中心差分
plot(N, ERRmaxf1r, '-sb');	// 前進差分に対するRomberg1段
plot(N, ERRmaxf2r, '-sg');	// 中心差分に対するRomberg1段
// *** 刻み幅を変えた際の値の変化 ***
scf(1);
a = gca();
a.data_bounds = [min(N),1E-14; max(N),1];
a.log_flags = "nl";
// 差の平均値
plot(N, DIFmeanf1, '-or');	// 前進差分
plot(N, DIFmeanf2, '-om');	// 中心差分
plot(N, DIFmeanf1r, '-ob');	// 前進差分に対するRomberg1段
plot(N, DIFmeanf2r, '-og');	// 中心差分に対するRomberg1段
// 差の最大値
plot(N, DIFmaxf1, '-sr');	// 前進差分
plot(N, DIFmaxf2, '-sm');	// 中心差分
plot(N, DIFmaxf1r, '-sb');	// 前進差分に対するRomberg1段
plot(N, DIFmaxf2r, '-sg');	// 中心差分に対するRomberg1段