掃き出しプログラム

#include <stdio.h>
#include <stdlib.h>
#include <math.h>  

main()
{
int i,j,k,a,b;
double m,n,h,h2,o;
static double L[700][701];

printf("分割数n=");
scanf("%lf",&n);

h = 1 / n;
h2 = h*h-2;
a = n+1;
b = n+2;

printf("n=%lf, h=%lf, h^2-2=%lf, a=%d, b=%d\n\n",n,h,h2,a,b);

/*行列入力*/
L[0][0] = 1;
L[n][n] = 1;

for(i=1; i<a-1; i++){
L[i][i] = h2;
if(i != 1){
L[i][i-1] = 1;
}
if(i != a-2){
L[i][i+1] = 1;
}
}

for(i=1; i<n; i++){
L[i][a] = -1*h*h*h*i;
}

/*掃き出し実行*/
for(k=0; k<a; k++){
for(i=0; i<a; i++){
if(i != k){
m = L[i][k];
for(j=0; j<b; j++){
L[i][j] = L[i][j] - ((m / L[k][k]) * L[k][j]);
}
}
}

m = L[k][k];
for(j=k; j<b; j++){
L[k][j] = L[k][j] / m;
}
}

/* 出力 */
for(i=0;i<a;i++){
m = (i*h);
if( m==.2 || m==.4 || m==.6 || m==.8 || (.599999<m && m<.600001)) printf("u%lf = %lf \n",m,L[i][a]);
}

return 0;
}

/*
10
u0.200000 = 0.036132
u0.400000 = 0.062842
u0.600000 = 0.071084
u0.800000 = 0.052550

20
u0.200000 = 0.036106
u0.400000 = 0.062798
u0.600000 = 0.071035
u0.800000 = 0.052514

40
u0.200000 = 0.036100
u0.400000 = 0.062787
u0.600000 = 0.071022
u0.800000 = 0.052505

80
u0.200000 = 0.036098
u0.400000 = 0.062784
u0.600000 = 0.071019
u0.800000 = 0.052503

160
u0.200000 = 0.036098
u0.400000 = 0.062783
u0.600000 = 0.071019
u0.800000 = 0.052503

320
u0.200000 = 0.036098
u0.400000 = 0.062783
u0.600000 = 0.071018
u0.800000 = 0.052503

640
u0.200000 = 0.036098
u0.400000 = 0.062783
u0.600000 = 0.071018
u0.800000 = 0.052502

参考文献：　http://ufcpp.net/study/linear/determinant.html#sweeping
*/