1 // 菊地文雄, 山本昌宏「微分方程式と計算機演習」, 山海堂 (1991), p. 58
2 /* ** FINITE DIFFERENCE METHOD **
3 ** FOR -D2U/Dx2=F **
4 N : NUMBER OF DIVISIONS
5 H : MESH SIZE
6 A,B,C : COEFFICIENT MATRIX
7 F : FREE TERM & SOLUTION VECTOR
8 */
9
10 #include <stdio.h>
11 #include <stdlib.h>
12 #define NDIM (100)
13
14 void trid(double a[], double b[], double c[], double f[], int n);
15
16 int main(void)
17 {
18 int n, i;
19 double h, h2;
20 double a[NDIM], b[NDIM], c[NDIM], f[NDIM], x[NDIM];
21 printf("INPUT : N=\n");
22 scanf("%d", &n);
23 if (n <= 1) exit(0);
24 h = 1.0 / n;
25 h2 = h * h;
26 for (i=0; i<n; i++) {
27 a[i] = -1.0 / h2;
28 b[i] = 2.0 / h2;
29 c[i] = -1.0 / h2;
30 f[i] = 1.0;
31 }
32 b[n-1] = 1.0 / h2;
33 f[n-1] = 0;
34 trid(a,b,c,f,n);
35 printf("MESH SIZE=%12.4e\n", h);
36 printf("** NODAL VALUES OF SOLUTION **\n");
37 for (i = 0; i < 5; i++)
38 printf(" I U(I) ");
39 printf("\n");
40 for (i=0; i<n; i++) {
41 printf("%3d%12.4e ", i+1, f[i]);
42 if ((i+1) % 5 == 0 || (i+1) == n)
43 printf("\n");
44 }
45 }
46
47 void trid(double a[], double b[], double c[], double f[], int n)
48 {
49 int i;
50 double aa;
51 // FORWARD ELIMITNATION
52 for (i = 0; i < n-1; i++) {
53 aa = a[i+1] / b[i];
54 b[i+1] -= aa * c[i];
55 f[i+1] -= aa * f[i];
56 }
57 // BACKWARD SUBSTITUTION
58 f[n-1] /= b[n-1];
59 for (i=n-2; i>= 0; i--)
60 f[i] = (f[i] - c[i] * f[i+1]) / b[i];
61 }
桂田 祐史
2018-06-08