/*
 * heat1d_i_ng.java -- 陰的スキームで熱方程式を解く
 */

import java.io.*;    // BufferedReader
import corejava.*;   // Format

public class heat1d_i_ng {
  static double pi = Math.PI;

  /* 初期条件を与える関数。複数登録して番号 nfunc で選択する。 */
  static double f(double x, int nfunc) {
    /* f(x)=max(x,1-x) */
    if (nfunc == 1) {
      if (x <= 0.5)
        return x;
      else
        return 1.0 - x;
    }
    /* f(x)=1 */
    else if (nfunc == 2)
      return 1.0;
    else if (nfunc == 3)
      return Math.sin(pi * x);
    /* f(x)= 変な形をしたもの(by M.Sakaue) */
    else if (nfunc == 4) {
      if (x <= 0.1)
        return 5 * x;
      else if (x <= 0.3)
        return -2 * x + 0.7;
      else if (x <= 0.5)
        return 4.5 * x - 1.25;
      else if (x <= 0.7)
        return -x + 1.5;
      else if (x <= 0.9)
        return x + 0.1;
      else
        return -10 * x + 10.0;
    }
    /* やはり非対称な形をしたもの(by M.Sakaue) */
    else if (nfunc == 5) {
      if (x <= 0.2)
        return -x + 0.8;
      else if (x <= 0.3)
        return 4 * x - 0.2;
      else if (x <= 0.4)
        return -4 * x + 2.2;
      else if (x <= 0.7)
        return -x + 1.0;
      else if (x <= 0.8)
        return 4 * x - 2.6;
      else
        return -x + 1.5;
    }
    else
      return 0.0;
  }

  /* 配列 u[] の内容(u[0],u[1],...,u[n])を一行 5 個ずつ表示する。 */
  static void print100(int N, double t, double [] u) {
    int i;
    Format fmt3d = new Format("%3d");
    Format fmt12e = new Format("%12.4e");
    Format.print(System.out, "T= %12.4e\n", t);
    for (i = 0; i < 5; i++)
      System.out.print("  I     u(i)    ");
    System.out.println("");
    for (i = 0; i <= N; i++) {
      System.out.print(fmt3d.form(i) + fmt12e.form(u[i]));
      if (i % 5 == 4)
        System.out.println("");
    }
    System.out.println("");
  }
  
  /* 三重対角行列の LU 分解 (pivoting なし) */
  static 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項方程式を解く */
  static 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];
  }

  /* 以前描いたグラフを消すか？ */
  private static final boolean erase_always = false;
  /* 数値を表示するか? */
  private static final boolean print_numerical_data = true;
  
  public static void main(String args[]) {
    int N, nfunc, i, n, nMax;
    double a, b, theta, h, Tmax, tau, c1, c2, c3, c4, lambda, t;
    double xleft, ybottom, xright, ytop;
    /* グラフを書き換える時間間隔 */
    double dt;
    /* 何ステップおきにグラフを書き換えるか */
    int skip;
    
    /* 区間の左端、右端 */
    a = 0.0;
    b = 1.0;

    try {
      /* 入力 */
      BufferedReader d = new BufferedReader(new InputStreamReader(System.in));
      System.out.print("入力して下さい : nfunc(1..5)=");
      nfunc = Integer.parseInt(d.readLine());
      System.out.print("入力して下さい : θ=");
      theta = Double.valueOf(d.readLine()).doubleValue();
      System.out.print("入力して下さい : N=");
      N = Integer.parseInt(d.readLine());
      System.out.print("入力して下さい : λ=");
      lambda = Double.valueOf(d.readLine()).doubleValue();
      
      /* ベクトルの準備 */
      double [] AL = new double [N];
      double [] AD = new double [N];
      double [] AU = new double [N];
      double [] ff = new double [N];
      double [] u  = new double [N+1];
      
      /* */
      h = (b - a) / N;
      tau = lambda * h * h;
      Format.print(System.out, " 時間の刻み幅τ= %g になりました。\n", tau);

      /* 入力 */
      System.out.print("入力して下さい : 最終時刻 Tmax=");
      Tmax = Double.valueOf(d.readLine()).doubleValue();
      System.out.print("入力して下さい : グラフ書き換え時間間隔 Δt=");
      dt = Double.valueOf(d.readLine()).doubleValue();
      skip = (int) Math.rint(dt / tau);
      
      /* 係数行列の準備 */
      c1 = - theta * lambda;
      c2 = 1.0 + 2.0 * theta * lambda;
      c3 = 1.0 - 2.0 * (1.0 - theta) * lambda;
      c4 = (1.0 - theta) * lambda;
      for (i = 1; i < N; i++) {
        AL[i] = c1;
        AD[i] = c2;
        AU[i] = c1;
      }
      
      /* 初期値 */
      for (i = 0; i <= N; i++)
        u[i] = f(a + i * h, nfunc);
      
      /* 出力（t=0） */
      /* 出力（ここでは画面に数値を表示すること）は FORTRAN と C で
         かなり異なるので、直接の翻訳は出来ない。ごたごたするので、
         print100 という関数にまとめた */
      t = 0.0;
      print100(N, t, u);
      
      /* 係数行列の LU 分解 */
      trilu1(N - 1, AL, AD, AU);
      
      /* 繰り返し計算 */
      nMax = (int) Math.rint(Tmax / tau);
      for (n = 1; n <= nMax; n++) {
        /* 連立1次方程式を作って解く */
        for (i = 1; i < N; i++)
          ff[i] = c3 * u[i] + c4 * (u[i - 1] + u[i + 1]);
        trisol1(N - 1, AL, AD, AU, ff);
        for (i = 1; i < N; i++)
          u[i] = ff[i];
        /* 境界条件 */
        u[0] = 0.0;
        u[N] = 0.0;
        
        /* 時刻 */
        t = n * tau;
        
        if (n % skip == 0) {
          /* 数値データの表示 */
          if (print_numerical_data)
            print100(N, t, u);
        }
      }
    }
    catch(IOException e) {
      System.err.println("IO Error");
      System.exit(1);
    }

  }
}