A..1 Newton-glsc.c

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

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

typedef double *vector;

vector new_vector(int n)
{
  return malloc(sizeof(double) * n);
}

double dotprod(int n, vector x, vector y)
{
  int i;
  double s;
  s=0;
  for (i=0; i< n; i++)
    s += x[i] * y[i];
  return s;
}

#define G_DOUBLE
#include <glsc.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(void)
{
  int N, i, k;
  vector al,ad,au,akl,akd,aku;
  vector U, x;
  double h, h2, du, H;
  double win_width, win_height, w_margin, h_margin;

  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;
  }

  /* GLSC初期化 */
  win_width = 150.0; win_height = 150.0; w_margin = 5.0; h_margin = 5.0;
  g_init("NewtonMeta", win_width  + 2 * w_margin, win_height + 2 * h_margin);
  /* 画面とメタファイルの両方に記録する */
  g_device(G_BOTH);
  /* 座標系の定義: [-0.2,1.2]×[-2.0, 20.0] という閉領域を表示する */
  g_def_scale(0,
              -0.2, 1.2, -2.0, 20.0,
              w_margin, h_margin, win_width, win_height);
  /* 線を二種類用意する */
  g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
  g_def_line(1, G_BLACK, 0, G_LINE_DOTS);
  /* 表示するための文字列の属性を定義する */
  g_def_text(0, G_BLACK, 3);
  /* 定義したものを選択する */
  g_sel_scale(0); g_sel_line(0); g_sel_text(0);

  /* タイトル表示 */
  g_text(40.0, 20.0, "-u''=u^2 in (0,1), u(0)=u(1)=0");
  /* 点線 */
  g_sel_line(1);
  /* 座標軸 */
  g_move(-0.2, 0.0); g_plot(1.2, 0.0); g_move(0.0, -2.0); g_plot(0.0, 20.0);
  /* 高さ 10 */
  g_move(-0.2, 10.0); g_plot(1.2, 10.0);
  for (k = 1; k < 100; k++) {
    /* A U^k の計算 */
    mul_mv(N - 1, x+1, al+1, ad+1, au+1, U+1);
    /* A_k=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
    g_sel_line(0);
    g_move(0.0, U[0]);
    for (i = 1; i <= N; i++)
      g_plot(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);
  }
  g_sleep(-1.0);
  g_term();
  return 0;
}