fork download
copy 
  1. #include <stdio.h>
  2. #define N 4
  3.  
  4. int main() {
  5. double A[N][N] = {
  6. {3, 1.5, -6, 4.8},
  7. {1, 1.5, -2, -2.4},
  8. {0, -1.5, -2, -1},
  9. {2, 4, -1.8, -0.6}
  10. };
  11.  
  12. double b[N] = {1.2, 0.6, -2.4, 0};
  13. double L[N][N] = {0};
  14. double U[N][N];
  15. double y[N] = {0};
  16. double x[N] = {0};
  17. double LU[N][N] = {0};
  18. double Ax[N] = {0};
  19.  
  20. int i, j, k;
  21.  
  22. for(i = 0; i < N; i++) {
  23. for(j = 0; j < N; j++) {
  24. U[i][j] = A[i][j];
  25. }
  26. }
  27.  
  28. for(i = 0; i < N; i++) {
  29. L[i][i] = 1.0;
  30. }
  31.  
  32. for(i = 0; i < N - 1; i++) {
  33. for(j = i + 1; j < N; j++) {
  34.  
  35. double factor = U[j][i] / U[i][i];
  36. L[j][i] = factor;
  37.  
  38. for(k = 0; k < N; k++) {
  39. U[j][k] = U[j][k] - factor * U[i][k];
  40. }
  41. }
  42. }
  43. for(i = 0; i < N; i++) {
  44.  
  45. y[i] = b[i];
  46.  
  47. for(j = 0; j < i; j++) {
  48. y[i] -= L[i][j] * y[j];
  49. }
  50. }
  51.  
  52. for(i = N - 1; i >= 0; i--) {
  53.  
  54. x[i] = y[i];
  55.  
  56. for(j = i + 1; j < N; j++) {
  57. x[i] -= U[i][j] * x[j];
  58. }
  59.  
  60. x[i] /= U[i][i];
  61. }
  62. printf("L行列:\n");
  63.  
  64. for(i = 0; i < N; i++) {
  65. for(j = 0; j < N; j++) {
  66. printf("%8.3f ", L[i][j]);
  67. }
  68. printf("\n");
  69. }
  70. printf("\nU行列:\n");
  71.  
  72. for(i = 0; i < N; i++) {
  73. for(j = 0; j < N; j++) {
  74. printf("%8.3f ", U[i][j]);
  75. }
  76. printf("\n");
  77. }
  78.  
  79. printf("\nyベクトル:\n");
  80.  
  81. for(i = 0; i < N; i++) {
  82. printf("%8.3f\n", y[i]);
  83. }
  84.  
  85.  
  86. printf("\n連立方程式の解 x:\n");
  87.  
  88. for(i = 0; i < N; i++) {
  89. printf("x%d = %8.3f\n", i + 1, x[i]);
  90. }
  91.  
  92. for(i = 0; i < N; i++) {
  93. for(j = 0; j < N; j++) {
  94.  
  95. LU[i][j] = 0;
  96.  
  97. for(k = 0; k < N; k++) {
  98. LU[i][j] += L[i][k] * U[k][j];
  99. }
  100. }
  101. }
  102.  
  103. printf("\nLU行列（LU = A の確認）:\n");
  104.  
  105. for(i = 0; i < N; i++) {
  106. for(j = 0; j < N; j++) {
  107. printf("%8.3f ", LU[i][j]);
  108. }
  109. printf("\n");
  110. }
  111.  
  112. for(i = 0; i < N; i++) {
  113.  
  114. Ax[i] = 0;
  115.  
  116. for(j = 0; j < N; j++) {
  117. Ax[i] += A[i][j] * x[j];
  118. }
  119. }
  120.  
  121. printf("\nAxベクトル（Ax = b の確認）:\n");
  122.  
  123. for(i = 0; i < N; i++) {
  124. printf("%8.3f\n", Ax[i]);
  125. }
  126.  
  127. printf("\n元のbベクトル:\n");
  128.  
  129. for(i = 0; i < N; i++) {
  130. printf("%8.3f\n", b[i]);
  131. }
  132. return 0;
  133. }
  134.  
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
L行列:
   1.000    0.000    0.000    0.000 
   0.333    1.000    0.000    0.000 
   0.000   -1.500    1.000    0.000 
   0.667    3.000   -1.100    1.000 

U行列:
   3.000    1.500   -6.000    4.800 
   0.000    1.000    0.000   -4.000 
   0.000    0.000   -2.000   -7.000 
   0.000    0.000    0.000    0.500 

yベクトル:
   1.200
   0.200
  -2.100
  -3.710

連立方程式の解 x:
x1 =   81.052
x2 =  -29.480
x3 =   27.020
x4 =   -7.420

LU行列（LU = A の確認）:
   3.000    1.500   -6.000    4.800 
   1.000    1.500   -2.000   -2.400 
   0.000   -1.500   -2.000   -1.000 
   2.000    4.000   -1.800   -0.600 

Axベクトル（Ax = b の確認）:
   1.200
   0.600
  -2.400
  -0.000

元のbベクトル:
   1.200
   0.600
  -2.400
   0.000
https://ideone.com/9zoTac
language:
C (gcc 8.3)
created:
1 day 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!