fork download 
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <stdlib.h>
  4.  
  5. #define EPS pow(10.0, -8.0)
  6. #define KMAX 100
  7. #define N 3
  8.  
  9. void newton2(double x, double y, double z);
  10. double *dvector(int i, int j);
  11. void free_dvector(double *a, int i);
  12. double **dmatrix(int nr1, int nr2, int nl1, int nl2);
  13. void free_dmatrix(double **a, int nr1, int nr2, int nl1, int nl2);
  14. double f(double x, double y, double z);
  15. double g(double x, double y, double z);
  16. double h(double x, double y, double z);
  17. double f_x(double x, double y, double z);
  18. double f_y(double x, double y, double z);
  19. double f_z(double x, double y, double z);
  20. double g_x(double x, double y, double z);
  21. double g_y(double x, double y, double z);
  22. double g_z(double x, double y, double z);
  23. double h_x(double x, double y, double z);
  24. double h_y(double x, double y, double z);
  25. double h_z(double x, double y, double z);
  26. double *gauss(double **a, double *b);
  27. double vector_norm1(double *a, int m, int n);
  28.  
  29. int main(void)
  30. {
  31.   double x, y, z;
  32.  
  33.   printf("初期値 x0, y0, z0 を入力してください : ");
  34.   scanf("%lf %lf %lf", &x, &y, &z);
  35.  
  36.   newton2(x, y, z);
  37.  
  38.   return 0;
  39. }
  40.  
  41. void newton2(double x, double y, double z)
  42. {
  43.   int i, k = 0;
  44.   double *xk, *d, **J;
  45.  
  46.   d = dvector(1, N);
  47.   xk = dvector(1, N);
  48.   J = dmatrix(1, N, 1, N);
  49.  
  50.   xk[1] = x; xk[2] = y; xk[3] = z;
  51.  
  52.   do {
  53.     /* 右辺ベクトル */
  54.     d[1] = -f(xk[1], xk[2], xk[3]);
  55.     d[2] = -g(xk[1], xk[2], xk[3]);
  56.     d[3] = -h(xk[1], xk[2], xk[3]);
  57.  
  58.     /* ヤコビ行列 */
  59.     J[1][1] = f_x(xk[1], xk[2], xk[3]);
  60.     J[1][2] = f_y(xk[1], xk[2], xk[3]);
  61.     J[1][3] = f_z(xk[1], xk[2], xk[3]);
  62.     J[2][1] = g_x(xk[1], xk[2], xk[3]);
  63.     J[2][2] = g_y(xk[1], xk[2], xk[3]);
  64.     J[2][3] = g_z(xk[1], xk[2], xk[3]);
  65.     J[3][1] = h_x(xk[1], xk[2], xk[3]);
  66.     J[3][2] = h_y(xk[1], xk[2], xk[3]);
  67.     J[3][3] = h_z(xk[1], xk[2], xk[3]);
  68.  
  69.     /* ガウス法 */
  70.     d = gauss(J, d);
  71.  
  72.     for (i = 1; i <= N; i++)
  73.       xk[i] += d[i];
  74.     k++;
  75.   } while (vector_norm1(d, 1, N) > EPS && k < KMAX);
  76.  
  77.   if (k == KMAX)
  78.     puts("解が見つかりませんでした。");
  79.   else
  80.     printf("解は x = %f, y = %f, z = %f\n", xk[1], xk[2], xk[3]);
  81.  
  82.   free_dmatrix(J, 1, N, 1, N);
  83.   free_dvector(d, 1);
  84.   free_dvector(xk, 1);
  85. }
  86.  
  87. double *dvector(int i, int j)
  88. {
  89.   double *a;
  90.  
  91.   if ((a = (double *)malloc(((j - i + 1) * sizeof(double)))) == NULL) {
  92.     puts("メモリが確保できません (from dvector)");
  93.     exit(1);
  94.   }
  95.  
  96.   return a - i; /* この動作はi > 0 で未定義な気がするが */
  97. }
  98.  
  99. void free_dvector(double *a, int i)
  100. {
  101.   free((void *)(a + i));
  102. }
  103.  
  104. double **dmatrix(int nr1, int nr2, int nl1, int nl2)
  105. {
  106.   int i, nrow, ncol;
  107.   double **a;
  108.  
  109.   nrow = nr2 - nr1 + 1;
  110.   ncol = nl2 - nl1 + 1;
  111.  
  112.   if ((a = (double **)malloc(nrow * sizeof(double *))) == NULL) {
  113.     puts("メモリが確保できません (from dmatrix 行)");
  114.     exit(1);
  115.   }
  116.  
  117.   a -= nr1;
  118.  
  119.   for (i = nr1; i <= nr2; i++)
  120.     if ((a[i] = (double *)malloc(ncol * sizeof(double))) == NULL) {
  121.       puts("メモリが確保できません (from dmatrix 列)");
  122.       exit(1);
  123.     }
  124.   for (i = nr1; i <= nr2; i++)
  125.     a[i] -= nl1;
  126.  
  127.   return a;
  128. }
  129.  
  130. void free_dmatrix(double **a, int nr1, int nr2, int nl1, int nl2)
  131. {
  132.   int i;
  133.  
  134.   for (i = nr1; i <= nr2; i++)
  135.     free((void *)(a[i] + nl1));
  136.   free((void *)(a + nr1));
  137. }
  138.  
  139. /* -y＾2 * z＾2 - y＾2 + 24*y*z - z＾2 -13 = 0 */
  140. double f(double x, double y, double z)
  141. {
  142.   return -y * y * z * z - y * y + 24.0 * y * z - z * z - 13.0;
  143. }
  144.  
  145. /* -x＾2 * z＾2 - x＾2 + 24*x*z - z＾2 -13 = 0 */
  146. double g(double x, double y, double z)
  147. {
  148.   return -x * x * z * z - x * x + 24.0 * x * z - z * z - 13.0;
  149. }
  150.  
  151. /* -x＾2 * y＾2 - x＾2 + 24*x*y - y＾2 -13 = 0 */
  152. double h(double x, double y, double z)
  153. {
  154.   return -x * x * y * y - x * x + 24.0 * x * y - y * y - 13.0;
  155. }
  156.  
  157. /* ここから偏微分 */
  158. double f_x(double x, double y, double z)
  159. {
  160.   return 0.0;
  161. }
  162.  
  163. double f_y(double x, double y, double z)
  164. {
  165.   return -2.0 * y * z * z - 2.0 * y + 24.0 * z;
  166. }
  167.  
  168. double f_z(double x, double y, double z)
  169. {
  170.   return -2.0 * z * y * y + 24.0 * y - 2.0 * z;
  171. }
  172.  
  173. double g_x(double x, double y, double z)
  174. {
  175.   return -2.0 * x * z * z - 2.0 * x + 24.0 * z;
  176. }
  177.  
  178. double g_y(double x, double y, double z)
  179. {
  180.   return 0.0;
  181. }
  182.  
  183. double g_z(double x, double y, double z)
  184. {
  185.   return -2.0 * z * x * x + 24.0 * x - 2.0 * z;
  186. }
  187.  
  188. double h_x(double x, double y, double z)
  189. {
  190.   return -2.0 * x * y * y - 2.0 * x + 24.0 * y;
  191. }
  192.  
  193. double h_y(double x, double y, double z)
  194. {
  195.   return -2.0 * y * x * x + 24.0 * x - 2.0 * y;
  196. }
  197.  
  198. double h_z(double x, double y, double z)
  199. {
  200.   return 0.0;
  201. }
  202.  
  203. /* ガウスの消去法(部分ピボット選択) */
  204. double *gauss(double **a, double *b)
  205. {
  206.   int i, j, k, ip;
  207.   double alpha, tmp;
  208.   double amax, eps = pow(2.0, -50.0);
  209.  
  210.   for (k = 1; k <= N - 1; k++) {
  211.     amax = fabs(a[k][k]); ip = k;
  212.     for (i = k + 1; i <= N; i++)
  213.       if (fabs(a[i][k]) > amax) {
  214.         amax = fabs(a[i][k]);
  215.         ip = i;
  216.       }
  217.  
  218.     if (amax < eps)
  219.       puts("行列は正則ではありません。");
  220.  
  221.     if (ip != k) {
  222.       for (j = k; j <= N; j++)
  223.         tmp = a[k][j], a[k][j] = a[ip][j], a[ip][j] = tmp;
  224.       tmp = b[k], b[k] = b[ip], b[ip] = tmp;
  225.     }
  226.  
  227.     for (i = k + 1; i <= N; i++) {
  228.       alpha = -a[i][k] / a[k][k];
  229.       for (j = k + 1; j <= N; j++)
  230.         a[i][j] += alpha * a[k][j];
  231.       b[i] += alpha * b[k];
  232.     }
  233.   }
  234.  
  235.   b[N] /= a[N][N];
  236.  
  237.   for (k = N - 1; k >= 1; k--) {
  238.     tmp = b[k];
  239.     for (j = k + 1; j <= N; j++)
  240.       tmp -= a[k][j] * b[j];
  241.     b[k] = tmp / a[k][k];
  242.   }
  243.  
  244.   return b;
  245. }
  246.  
  247. /* 1-ノルム */
  248. double vector_norm1(double *a, int m, int n)
  249. {
  250.   int i;
  251.   double norm = 0.0;
  252.  
  253.   for (i = m; i <= n; i++)
  254.     norm += fabs(a[i]);
  255.  
  256.   return norm;
  257. }
  258.  
Not running #stdin #stdout 0s 0KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
https://ideone.com/5quhq
language:
C (gcc 8.3)
created:
12 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!