fork(1) download 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define EPSILON 1e-5 
  5. #define MAX_ITER 100
  6.  
  7. double f1(double x, double y) { return x * x + y * y - 4; }
  8. double f2(double x, double y) { return pow(x, 4) + x - y * y + 1; }
  9.  
  10. double df1_dx(double x, double y) { return 2 * x; }
  11. double df1_dy(double x, double y) { return 2 * y; }
  12. double df2_dx(double x, double y) { return 4 * x * x * x + 1; }
  13. double df2_dy(double x, double y) { return -2 * y; }
  14.  
  15. void solve(int attempt, double x0, double y0) {
  16.     double x = x0, y = y0;
  17.     int iter = 0;
  18.  
  19.     printf("\n--- 試行 #%d: (初期値 x0=%.2f, y0=%.2f) ---\n", attempt, x0, y0);
  20.     printf("%-5s %-15s %-15s\n", "反復", "x", "y");
  21.  
  22.     while (iter <= MAX_ITER) {
  23.         double F1 = f1(x, y);
  24.         double F2 = f2(x, y);
  25.  
  26.         // 各反復のx, yを表示
  27.         printf("%-5d %-15.10f %-15.10f\n", iter, x, y);
  28.  
  29.         // 収束判定: |f1| < 10^-5 かつ |f2| < 10^-5
  30.         if (fabs(F1) < EPSILON && fabs(F2) < EPSILON) {
  31.             printf(">> 結果: 反復回数 %d 回で収束しました。\n", iter);
  32.             printf(">> 解: (x = %.5f, y = %.5f)\n", x, y);
  33.             return;
  34.         }
  35.  
  36.         double J11 = df1_dx(x, y);
  37.         double J12 = df1_dy(x, y);
  38.         double J21 = df2_dx(x, y);
  39.         double J22 = df2_dy(x, y);
  40.  
  41.         double det = J11 * J22 - J12 * J21;
  42.         if (fabs(det) < 1e-12) {
  43.             printf(">> 結果: 失敗 (ヤコビ行列の行列式が小さすぎます)。\n");
  44.             return;
  45.         }
  46.  
  47.         x += (-F1 * J22 + F2 * J12) / det;
  48.         y += (-J11 * F2 + J21 * F1) / det;
  49.         iter++;
  50.     }
  51.     printf(">> 結果: 失敗 (%d 回の反復内に収束しませんでした)。\n", MAX_ITER);
  52. }
  53.  
  54. int main() {
  55.     // 異なる3組の初期値(x0, y0)で計算すること
  56.     double initial_x[] = {1.2, -1.0, 0.5};
  57.     double initial_y[] = {1.5, 1.0, -1.8};
  58.  
  59.     for (int i = 0; i < 3; i++) {
  60.         solve(i + 1, initial_x[i], initial_y[i]);
  61.     }
  62.  
  63.     return 0;
  64. }
Success #stdin #stdout 0s 5332KB
comments (?)
stdin 
Standard input is empty
stdout 
--- 試行 #1: (初期値 x0=1.20, y0=1.50) ---
反復 x               y              
0     1.2000000000    1.5000000000   
1     1.0338246703    1.7362735971   
2     1.0011115000    1.7317235976   
3     1.0000012343    1.7320504817   
>> 結果: 反復回数 3 回で収束しました。
>> 解: (x = 1.00000, y = 1.73205)

--- 試行 #2: (初期値 x0=-1.00, y0=1.00) ---
反復 x               y              
0     -1.0000000000   1.0000000000   
1     -1.4000000000   1.6000000000   
2     -1.2902943018   1.5334924859   
3     -1.2759112634   1.5402281288   
4     -1.2756822609   1.5403359482   
>> 結果: 反復回数 4 回で収束しました。
>> 解: (x = -1.27568, y = 1.54034)

--- 試行 #3: (初期値 x0=0.50, y0=-1.80) ---
反復 x               y              
0     0.5000000000    -1.8000000000  
1     1.3750000000    -1.6986111111  
2     1.1035857261    -1.6899229972  
3     1.0098267771    -1.7293333385  
4     1.0000957562    -1.7320249467  
5     1.0000000092    -1.7320508051  
>> 結果: 反復回数 5 回で収束しました。
>> 解: (x = 1.00000, y = -1.73205)
https://ideone.com/mT1sGm
language:
C (gcc 8.3)
created:
12 hours 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!