D.5 サンプル・プログラム

/*
 * Newton.c
 *  非線形 2 点境界値問題
 *    -u''=u^2 in (0,1), u(0)=u(1)=0
 *  を差分法で離散化して得られる非線形方程式を Newton 法で解く。
 */

#include <stdio.h>
#include <math.h>
#include <matrix.h>
#include "fplot.h"
#include "trid-lu.h"

void mul_mv(int n,
            vector ab,
            vector al, vector ad, vector au,
            vector b)
{
  int i, nm1 = n - 1;
  ab[0]   =                  ad[0] * b[0] + au[0] * b[1];
  for (i = 1; i < nm1; i++)
    ab[i] = al[i] * b[i-1] + ad[i] * b[i] + au[i] * b[i+1];
  ab[nm1]   = al[nm1] * b[nm1-1] + ad[nm1] * b[nm1];
}

double norm(int n, vector x)
{
  return sqrt(dotprod(n, x, x));
}

int main()
{
  int N, i, k;
  vector al,ad,au,akl,akd,aku;
  vector U, x;
  double h, h2, du, H;

  N = 100;
  h = 1.0 / N;
  h2 = h * h;
  al = new_vector(N+1); ad = new_vector(N+1); au = new_vector(N+1);
  akl = new_vector(N+1); akd = new_vector(N+1); aku = new_vector(N+1);
  U = new_vector(N+1);
  x = new_vector(N+1);

  /* 初期値 */
  printf("H (10位でOK)="); scanf("%lg", &H);
  for (i = 0; i <= N; i++)
    U[i] = H;
  U[0] = U[N] = 0.0;
  /* A */
  for (i = 1; i < N; i++) {
    al[i] = - 1.0 / h2; ad[i] = 2.0 / h2; au[i] = - 1.0 / h2;
  }
  openpl(); fspace2(-0.2, -2.0, 1.2, 20.0);

  fmove(0.3, 15.0);
  label("-u''=u^2 in (0,1), u(0)=u(1)=0");
  linemod("dotted");
  fline(-0.2, 0.0, 1.2, 0.0); fline(0.0, -2.0, 0.0, 20.0);
  linemod("solid");

  for (k = 1; k < 100; k++) {
    /* A U^k の計算 */
    mul_mv(N - 1, x+1, al+1, ad+1, au+1, U+1);
    /* F(U^k) の計算 */
    for (i = 1; i < N; i++)
      x[i] -= U[i] * U[i];
    /* F'(U^k) の計算 */
    for (i = 1; i < N; i++) {
      akl[i] = al[i]; akd[i] = ad[i] - 2 * U[i]; aku[i] = au[i];
    }
    /* F'(U^k)^{-1} U(U^k) の計算 */
    trid(N-1, akl+1, akd+1, aku+1, x+1);
    /* */
    du = norm(N-1, x+1);
    printf("du=%g\n", du);
    /* U^{k+1} の計算 */
    for (i = 1; i < N; i++)
      U[i] -= x[i];
    /* */
#ifdef NONE
    for (i = 0; i <= N; i++)
      printf("U[%d]=%g\n", i, U[i]);
#endif
    fmove(0.0, U[0]);
    for (i = 1; i <= N; i++)
      fcont(i * h, U[i]);
    if (du < 1.0e-12)
      break;
  }
  {
    double min, max;
    max = U[0]; min = U[0];
    for (i = 1; i <= N; i++) {
      if (U[i] > max)
        max = U[i];
      else if (U[i] < min)
        min = U[i];
    }
    printf("min=%g, max=%g\n", min, max);
  }
  mkplot("Newton.plot");
  closepl();
  return 0;
}