fork download 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. // 円周率の定義
  5. #define PI 3.14159265358979323846
  6.  
  7. // 対象の関数 f(x) = 1 / (2 + cos(pi * x))
  8. double f(double x) {
  9.     return 1.0 / (2.0 + cos(PI * x));
  10. }
  11.  
  12. // ラグランジュ基底多項式 l_j(x) の計算
  13. double lagrange_basis(int j, double x, const double nodes[], int n) {
  14.     double l = 1.0;
  15.     for (int k = 0; k <= n; k++) {
  16.         if (k != j) {
  17.             l *= (x - nodes[k]) / (nodes[j] - nodes[k]);
  18.         }
  19.     }
  20.     return l;
  21. }
  22.  
  23. // 補間多項式 P_n(x) の値の計算
  24. double interpolant(double x, const double nodes[], const double values[], int n) {
  25.     double p = 0.0;
  26.     for (int j = 0; j <= n; j++) {
  27.         p += values[j] * lagrange_basis(j, x, nodes, n);
  28.     }
  29.     return p;
  30. }
  31.  
  32. int main(void) {
  33.     int i, j, k;
  34.     int n = 7; // 7次のチェビシェフ補間（節点数 8個）
  35.     double nodes[8];
  36.     double values[8];
  37.  
  38.     // 0. 積分の真値を高精度（十分な分割数による数値積分）で近似
  39.     int ref_steps = 100000;
  40.     double a = -1.0, b = 1.0;
  41.     double h_ref = (b - a) / ref_steps;
  42.     double true_integral = 0.5 * (f(a) + f(b));
  43.     for (i = 1; i < ref_steps; i++) {
  44.         true_integral += f(a + i * h_ref);
  45.     }
  46.     true_integral *= h_ref;
  47.     printf("積分の真値 (近似): %.12f\n\n", true_integral);
  48.  
  49.     // ==========================================
  50.     // (1) 7次のチェビシェフ補間の節点と値の計算
  51.     // ==========================================
  52.     for (k = 0; k <= n; k++) {
  53.         // チェビシェフ節点: x_k = cos((2k+1)pi / (2(n+1)))
  54.         nodes[k] = cos((2.0 * k + 1.0) * PI / (2.0 * (n + 1.0)));
  55.         values[k] = f(nodes[k]);
  56.     }
  57.  
  58.     printf("--- (1) 7次チェビシェフ補間の節点 ---\n");
  59.     for (k = 0; k <= n; k++) {
  60.         printf("x_%d = %10.6f, f(x_%d) = %10.6f\n", k, nodes[k], k, values[k]);
  61.     }
  62.     printf("\n");
  63.  
  64.     // ==========================================
  65.     // (2) 補間多項式の [-1, 1] での積分と真値との誤差
  66.     // ==========================================
  67.     // 補間多項式を十分に細かい分割（例: 1000分割）で数値積分する
  68.     int interp_steps = 1000;
  69.     double h_interp = (b - a) / interp_steps;
  70.     double poly_integral = 0.5 * (interpolant(a, nodes, values, n) + interpolant(b, nodes, values, n));
  71.     for (i = 1; i < interp_steps; i++) {
  72.         double xi = a + i * h_interp;
  73.         poly_integral += interpolant(xi, nodes, values, n);
  74.     }
  75.     poly_integral *= h_interp;
  76.  
  77.     double error_interp = fabs(poly_integral - true_integral);
  78.  
  79.     printf("--- (2) 補間多項式の積分と誤差 ---\n");
  80.     printf("補間多項式の積分値: %.12f\n", poly_integral);
  81.     printf("真値との誤差      : %.2e\n\n", error_interp);
  82.  
  83.     // ==========================================
  84.     // (3) 7分割・8分割の台形則による数値積分
  85.     // ==========================================
  86.     // 7分割の台形則
  87.     int n7 = 7;
  88.     double h7 = (b - a) / n7;
  89.     double trap_7 = 0.5 * (f(-1.0) + f(1.0));
  90.     for (i = 1; i < n7; i++) {
  91.         trap_7 += f(-1.0 + i * h7);
  92.     }
  93.     trap_7 *= h7;
  94.     double err_trap_7 = fabs(trap_7 - true_integral);
  95.  
  96.     // 8分割の台形則
  97.     int n8 = 8;
  98.     double h8 = (b - a) / n8;
  99.     double trap_8 = 0.5 * (f(-1.0) + f(1.0));
  100.     for (i = 1; i < n8; i++) {
  101.         trap_8 += f(-1.0 + i * h8);
  102.     }
  103.     trap_8 *= h8;
  104.     double err_trap_8 = fabs(trap_8 - true_integral);
  105.  
  106.     printf("--- (3) 台形則による数値積分と誤差 ---\n");
  107.     printf("7分割の台形則 積分値: %.12f, 誤差: %.2e\n", trap_7, err_trap_7);
  108.     printf("8分割の台形則 積分値: %.12f, 誤差: %.2e\n", trap_8, err_trap_8);
  109.  
  110.     return 0;
  111. }
Success #stdin #stdout 0.01s 5320KB
comments (?)
stdin 
Standard input is empty
stdout 
積分の真値 (近似): 1.154700538379

--- (1) 7次チェビシェフ補間の節点 ---
x_0 =   0.980785, f(x_0) =   0.998182
x_1 =   0.831470, f(x_1) =   0.879572
x_2 =   0.555570, f(x_2) =   0.547553
x_3 =   0.195090, f(x_3) =   0.354863
x_4 =  -0.195090, f(x_4) =   0.354863
x_5 =  -0.555570, f(x_5) =   0.547553
x_6 =  -0.831470, f(x_6) =   0.879572
x_7 =  -0.980785, f(x_7) =   0.998182

--- (2) 補間多項式の積分と誤差 ---
補間多項式の積分値: 1.154838461837
真値との誤差      : 1.38e-04

--- (3) 台形則による数値積分と誤差 ---
7分割の台形則 積分値: 1.154929577465, 誤差: 2.29e-04
8分割の台形則 積分値: 1.154761904762, 誤差: 6.14e-05
https://ideone.com/mLivZX
language:
C (gcc 8.3)
created:
2 days 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!