B..2.2 C言語のプログラム例

/*
 * poisson1d.c --- 1次元Poisson方程式
 *
 *   - u''(x)=f(x)      (0<x<1)
 *   u(0)=α
 *   u'(1)=β
 *
 * を差分法で解くプログラム。cglsc コマンドでコンパイル＆実行出来る。
 *
 * u(x)=(x-1/3)^2 の場合、α=1/9, β=4/3  --- 2次関数なので厳密に解ける
 */

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

#define MAXN (100);

// 三重対角行列の Gauss の消去法による LU 分解
void trilu1(int, double *, double *, double *);
// LU分解を用いて三項方程式を解く
void trisol1(int, double *, double *, double *, double *);

// Poisson方程式の右辺
double f(double x)
{
  return - 2.0;
}

// 厳密解
double solution(double x)
{
  return (x - 1.0 / 3) * (x - 1.0 / 3); 
}

int main(void)
{
  double alpha, beta;
  int i, n;
  double *x, *al, *ad, *au, *u;
  double h, h2;
  alpha = 1.0 / 9.0;
  beta = 4.0 / 3.0;

  printf("n="); scanf("%d", &n);

  x = malloc(sizeof(double) * (n + 1));
  al = malloc(sizeof(double) * (n + 1));
  ad = malloc(sizeof(double) * (n + 1));
  au = malloc(sizeof(double) * (n + 1));
  u = malloc(sizeof(double) * (n + 1));
  if (ad == NULL || al == NULL || au == NULL || u == NULL) {
    fprintf(stderr, "メモリの割り当てに失敗しました。\n");
    exit(1);
  }
  h = 1.0 / n;
  for (i = 0; i <= n; i++)
    x[i] = i * h;
  h2 = h * h;
  for (i = 1; i <= n; i++) {
    ad[i] = 2;
    au[i] = al[i] = -1;
    u[i] = h2 * f(i * h);
  }
  /* 第N成分は 1/2 倍する */
  ad[n] = 1;
  u[n] *= 0.5;
  /* 非同次境界条件 */
  u[1] += alpha;
  u[n] += beta * h;
  /* 連立1次方程式を解く */
  trilu1(n, al, ad, au);
  trisol1(n, al, ad, au, u);

  /* グラフの準備 */
  u[0] = alpha;
  g_init("POISSON1D", 170.0, 170.0);
  g_device(G_BOTH);
  g_def_scale(0,
	      -0.2, 1.2, -0.2, 1.2,
	      10.0, 10.0, 150.0, 150.0);
  g_sel_scale(0);
  g_move(-0.2, 0.0); g_plot(1.2, 0.0);
  g_move(0.0, -0.2); g_plot(0.0, 1.2);

  /* 厳密解を描く */
  g_move(x[0], solution(x[0]));
  for (i = 1; i <= n; i++)
    g_plot(x[i], solution(x[i]));
  printf("クリックして下さい\n");
  g_sleep(-1.0);
  /* 差分解を赤く描く */
  g_line_color(G_RED);
  g_move(x[0], u[0]);
  for (i = 1; i <= n; i++)
    g_plot(x[i], u[i]);

  printf("x=1 での厳密解 %f, 差分解 %f, 誤差 %e\n",
          solution(1.0), u[n], fabs(solution(1.0) - u[n]));
  g_sleep(-1.0);
  g_term();
  
  return 0;
}

/* 3項方程式 (係数行列が三重対角である連立1次方程式のこと) Ax=b を解く
 *
 *   入力
 *     n: 未知数の個数
 *     al,ad,au: 連立1次方程式の係数行列
 *       (al: 対角線の下側 i.e. 下三角部分 (lower part)
 *        ad: 対角線       i.e. 対角部分 (diagonal part)
 *        au: 対角線の上側 i.e. 上三角部分 (upper part)
 *        つまり
 *
 *         ad[1] au[1]   0   .................. 0
 *         al[2] ad[2] au[2]   0  ............. 0
 *           0   al[3] ad[3] au[3]  0 ......... 0
 *                         ....................
 *                         al[n-1] ad[n-1] au[n-1]
 *                           0     al[n]   ad[n]

 *        al[i] = A_{i,i-1}, ad[i] = A_{i,i}, au[i] = A_{i,i+1},
 *        al[1], au[n] は意味がない)
 *
 *     b: 連立1次方程式の右辺の既知ベクトル
 *        (添字は 1 から。i.e. b[1],b[2],...,b[n] にデータが入っている。)
 *   出力
 *     al,ad,au: 入力した係数行列を LU 分解したもの
 *     b: 連立1次方程式の解
 *   能書き
 *     一度 call すると係数行列を LU 分解したものが返されるので、
 *     以後は同じ係数行列に関する連立1次方程式を解くために、
 *     関数 trisol1() が使える。
 *   注意
 *     Gauss の消去法を用いているが、ピボットの選択等はしていな
 *     いので、ピボットの選択をしていないので、係数行列が正定値である
 *     などの適切な条件がない場合は結果が保証できない。
 */

void trid1(int n, double *al, double *ad, double *au, double *b)
{
  trilu1(n,al,ad,au);
  trisol1(n,al,ad,au,b);
}

/* 三重対角行列の LU 分解 (pivoting なし) */
void trilu1(int n, double *al, double *ad, double *au)
{
  int i;
  /* 前進消去 (forward elimination) */
  for (i = 1; i < n; i++) {
    al[i + 1] /= ad[i];
    ad[i + 1] -= au[i] * al[i + 1];
  }
}

/* LU 分解済みの三重対角行列を係数に持つ3項方程式を解く */
void trisol1(int n, double *al, double *ad, double *au, double *b)
{
  int i;
  /* 前進消去 (forward elimination) */
  for (i = 1; i < n; i++) b[i + 1] -= b[i] * al[i + 1];
  /* 後退代入 (backward substitution) */
  b[n] /= ad[n];
  for (i = n - 1; i >= 1; i--) b[i] = (b[i] - au[i] * b[i + 1]) / ad[i];
}