6.2.2 HeatCanvas_all.java

// HeatCanvas_all.java

import java.awt.*;
import java.awt.event.*;

public class HeatCanvas_all extends Canvas {

  static final boolean DEBUG = false; // true;
  private static final String message = "1D Heat Equation";
  // 問題にとって基本的なパラメーター
  //   N:      区間の分割数
  //   Tmax:   最終時刻
  //   lambda: λ=τ/h^2
  //   A,B: 　 境界条件
  //   theta 　θ法
  //   erabuBC 境界条件のチョイスに用いた変数
  //   erabuF  f(x)のチョイスに用いた変数
  private int N, erabuBC, erabuF, anime;
  private double Tmax, lambda, theta, A, B;
  // 座標系の変換のためのパラメーター
  private int CanvasX = 400, CanvasY = 400;
  private double ratiox, ratioy, X0, Y0;
  //キャンバス
  private HeatCanvas_all c1;
  private Thread th;

  // コンストラクター
  public HeatCanvas_all() {
    super();
    space(-0.1, -0.1, 1.1, 1.1);
  }
  public HeatCanvas_all(int cx, int cy) {
    super();
    CanvasX = cx; CanvasY = cy;
    space(-0.1, -0.1, 1.1, 1.1);
  }

  public void compute(double lambda, int N, 
                        double Tmax, double theta, double A, double B) {
    this.lambda = lambda;
    this.N = N;
    this.Tmax = Tmax;
    this.theta = theta;
    this.A = A;
    this.B = B;
    repaint();
  }

  public void compute(double lambda, int N, 
                        double Tmax, double theta, 
                        double A, double B, int erabuBC) {
    this.lambda = lambda;
    this.N = N;
    this.Tmax = Tmax;
    this.theta = theta;
    this.A = A;
    this.B = B;
    this.erabuBC = erabuBC;
    repaint();
  }

  public void compute(double lambda, int N, 
                        double Tmax, double theta, 
                        double A, double B, int erabuBC, int erabuF) {
    this.lambda = lambda;
    this.N = N;
    this.Tmax = Tmax;
    this.theta = theta;
    this.A = A;
    this.B = B;
    this.erabuBC = erabuBC;
    this.erabuF = erabuF;
    
    repaint();
  }

  public void compute(double lambda, int N, 
                        double Tmax, double theta, 
                        double A, double B, int erabuBC, int erabuF, int anime) 
  {
    this.lambda = lambda;
    this.N = N;
    this.Tmax = Tmax;
    this.theta = theta;
    this.A = A;
    this.B = B;
    this.erabuBC = erabuBC;
    this.erabuF = erabuF;
    this.anime = anime;

    repaint();
  }

  // 初期条件
  private double f(double x, int erabuF) {
    if (erabuF == 0) {
      if (x <= 0.5)
        return x;
      else
        return 1.0 - x;
    }
    else if (erabuF == 1) {
      return 1.0;
    }
    else if (erabuF == 2){
      return Math.sin(2 * x * Math.PI);
    }
    else
      return x;             
  }
  // 座標変換の準備
  private void space(double x0, double y0, double x1, double y1) {
    X0 = x0; Y0 = y0;
    ratiox = CanvasX / (x1 - x0);
    ratioy = CanvasY / (y1 - y0);
  }
  // ユーザー座標 (ワールド座標系) をウィンドウ座標 (デバイス座標系)
  private int wx(double x) {
    return (int)(ratiox * (x - X0));
  }
  // ユーザー座標 (ワールド座標系) をウィンドウ座標 (デバイス座標系)
  private int wy(double y) {
    return CanvasY - (int)(ratioy * (y - Y0));
  }
  // x[], y[] の内容をグラフにする
  private void drawGraph(Graphics g, double x[], double u[]) {
    for (int i = 0; i < N; i++)
      g.drawLine(wx(x[i]), wy(u[i]), wx(x[i + 1]), wy(u[i + 1]));
  }

  // 三重対角行列のLU分解
  private void trilu(int N, double [] AL, double [] AD, double [] AU) {
    for (int i = 1; i < N; i++) {
      AL[i+1] /= AD[i];
      AD[i+1] -= AU[i] * AL[i+1];
    } 
  }
  // 三項方程式を解く
  private void trisol(int N, double [] AL, double [] AD, 
  double [] AU, double [] F) {
    int i;
    for (i = 1; i < N; i++)
      F[i+1] -= F[i] * AL[i+1]; 
    F[N] /= AD[N];   
    for (i = N - 1; i >= 1; i--)
      F[i] = (F[i] - AU[i] * F[i+1]) / AD[i];
  }   
  
  public void paint(Graphics g) {
    int i;
    double h = 1.0 / N;
    double tau = lambda * h * h;
    double dt = 0.01;
    int Jmax = (int) (Tmax / tau);
    int skip = (int) (dt / tau + 0.5);
    double c1 = -theta * lambda;
    double c2 = 1.0 + 2.0 * theta * lambda;
    double c3 = 1.0 - 2.0 * (1.0 - theta) * lambda;
    double c4 = (1.0 - theta) * lambda;
    double AL[] = new double[N+1];
    double AD[] = new double[N+1];
    double AU[] = new double[N+1];
    double ff[] = new double[N+1];
    double u[] = new double[N+1];
    double x[] = new double[N+1];
    for (i = 1; i < N; i++) {
      AL[i] = c1; AD[i] = c2; AU[i] = c1;
    }
    // BC (Boundary Condition) を選択した時に変更
    if (erabuBC == 1) {
      AD[1] = 1.0 + theta * lambda;
      AD[N-1] = 1.0 + theta * lambda;
    }
    else if(erabuBC == 2) {
      AD[1] = 1.0 + theta * lambda;
    }
    else if(erabuBC == 3) {
      AD[N-1] = 1.0 + theta * lambda;
    }
    
    for (i = 0; i <= N; i++) 
      x[i] = i * h;

    // タイトルを表示する
    g.setColor(Color.red);
    g.setFont(new Font("Helvetica", Font.BOLD, 24));
    g.drawString(message, 40, 50);

    // 初期値を計算する
    for (i = 0; i <= N; i++)
      u[i] = f(i * h, erabuF);

    // 初期値のグラフを描く
    g.setColor(Color.black);
    drawGraph(g, x, u);

    trilu(N - 1, AL, AD, AU);

    // 時間に関するループ
    for (int j = 1; j <= Jmax; j++) {
      ff[1] = c3 * u[1] + c4 * (u[0] + u[2])+ A * theta * lambda;
      ff[N-1] = c3 * u[N-1] + c4 * (u[N-2] + u[N]) + B * theta * lambda;

      // BC (Boundary Condition) を選択した時に変更
      if (erabuBC == 1) {
        ff[1] = c3 * u[1] + c4 * (u[0] + u[2])- A * theta * lambda * h;
        ff[N-1] = c3 * u[N-1] + c4 * (u[N-2] + u[N]) + B * theta * lambda * h;
      }
      else if (erabuBC == 2) {
        ff[1] = c3 * u[1] + c4 * (u[0] + u[2])- A * theta * lambda * h;
      }
      else if (erabuBC == 3) {
        ff[N-1] = c3 * u[N-1] + c4 * (u[N-2] + u[N]) + B * theta * lambda * h;
      }
      for (i = 2; i < N - 1; i++) 
          ff[i] = c3 * u[i] + c4 * (u[i-1] + u[i+1]);

      trisol(N - 1, AL, AD, AU, ff);

      for (i = 1; i < N; i++)
        u[i] = ff[i];

      u[0] = A; u[N] = B;

      // BC (Boundary Condition) を選択した時に変更
      if (erabuBC == 1) {
        u[0] = u[1] - A * h; u[N] = u[N-1] + B * h;
      }
      else if (erabuBC == 2) {
        u[0] = u[1] - A * h;
      }
      else if (erabuBC == 3) {
        u[N] = u[N-1] + B * h;
      }

      if (DEBUG)
        for (i = 0; i <= N; i++)
          System.out.println("u[" + i + "]=" + u[i]);
      // 適当な間隔 (dt=0.01) でグラフを描く
      if (j % skip == 0)
        drawGraph(g, x, u);
      try {
        th.sleep(10 * anime);
      } catch(InterruptedException e){}
    }
  }
}