fork(1) download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. // 計算する値の個数
  5. #define N 20
  6.  
  7. // 円周率
  8. #define M_PI acos(-1.0)
  9.  
  10. // 微分係数を計算する点
  11. #define x0 (M_PI / 8)
  12.  
  13. // 真の値 = sqrt(2)
  14. #define d0 sqrt(2)
  15.  
  16. // 関数 f
  17. double f(double x) { return sin(2 * x); }
  18.  
  19. // 相対誤差を計算する関数
  20. double err(double x) { return fabs(x - d0) / d0; }
  21.  
  22. // D_f
  23. double dif_f(double x, double h) { return (f(x + h) - f(x)) / h; }
  24.  
  25. // D_m
  26. double dif_m(double x, double h) { return (f(x + h) - f(x - h)) / (2 * h); }
  27.  
  28. // 結果は引数として受け取った配列 b に格納する
  29. // 配列 a の長さは L で与えられ、結果は b[0] から b[L-2] までに入る
  30. void accel(double a[], int L, int p, double b[]) {
  31.     for (int k = 0; k < L - 1; k++)
  32.         b[k] = (pow(2, p) * a[k + 1] - a[k]) / (pow(2, p) - 1);
  33. }
  34.  
  35. // オーダーを計算
  36. // 危険：算術エラーに対する処理を行わないので、サンドボックスでの実験以外で使ってはならない
  37. void order(double a[], double L) {
  38.     for (int k = 0; k < L - 2; k++)
  39.         printf("%2d: %lf\n", k, -log2((a[k + 2] - a[k + 1]) / (a[k + 1] - a[k])));
  40. }
  41.  
  42. // 計算結果を表示
  43. void print_result(double a[], double h[], int L) {
  44.     for(int k = 0; k < L; k++)
  45.         printf("%2d: %e  %.16lf  %e\n", k, h[k], a[k], err(a[k]));
  46. }
  47.  
  48. int main(void) {
  49.     double h[N], a[4][N];
  50.  
  51.     for(int k = 0; k < N; k++) {
  52.         // 刻み幅
  53.         h[k] = pow(2, -k - 3);
  54.  
  55.         // D_m の計算値
  56.         a[0][k] = dif_m(x0, h[k]);
  57.     }
  58.  
  59.     //order(a[0], N);
  60.     accel(a[0], N, 2, a[1]);
  61.  
  62.     //order(a[1], N - 1);
  63.     accel(a[1], N - 1, 4, a[2]);
  64.  
  65.     print_result(a[2], h, 10);
  66.  
  67.     return 0;
  68. }
Success #stdin #stdout 0s 10320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
 0: 1.250000e-01  1.4142135613039182  7.560223e-10
 1: 6.250000e-02  1.4142135623563752  1.182280e-11
 2: 3.125000e-02  1.4142135623728378  1.819737e-13
 3: 1.562500e-02  1.4142135623730776  1.240373e-14
 4: 7.812500e-03  1.4142135623731134  1.287476e-14
 5: 3.906250e-03  1.4142135623731000  3.454203e-15
 6: 1.953125e-03  1.4142135623730718  1.648597e-14
 7: 9.765625e-04  1.4142135623730308  4.553268e-14
 8: 4.882812e-04  1.4142135623734347  2.400671e-13
 9: 2.441406e-04  1.4142135623730057  6.327473e-14
https://ideone.com/2ciBgi
language:
C (gcc 8.3)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!