///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: LinearDEq.java                                     //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of the book, "An Introduction to        //
//     Computational Physics, 2nd Edition," written by Tao Pang and      //
//     published by Cambridge University Press on January 19, 2006.      //
//                                                                       //
// (2) No warranties, express or implied, are made for this program.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// An example of solving the boundary value problem of
// a linear equation via the Runge-Kutta method.

import java.lang.*;
public class LinearDEq {
  static final int n = 100, m = 5;
  public static void main(String argv[]) {
    double y1[] = new double [n+1];
    double y2[] = new double [n+1];
    double y[] = new double [2];
    double h = 1.0/n;

 // Find the 1st solution via Runge-Kutta method
    y[1]  = 1;
    for (int i=0; i<n; ++i) {
      double x = h*i;
      y = rungeKutta(y, x, h);
      y1[i+1] = y[0];
    }

 // Find the 2nd solution via Runge-Kutta method
    y[0] = 0;
    y[1] = 2;
    for (int i=0; i<n; ++i) {
      double x = h*i;
      y = rungeKutta(y, x, h);
      y2[i+1] = y[0];
    }

 // Superpose two solutions found
    double a = (y2[n]-1)/(y2[n]-y1[n]);
    double b = (1-y1[n])/(y2[n]-y1[n]);
    for (int i=0; i<=n; ++i)
      y1[i] = a*y1[i]+b*y2[i];

 // Output the result in every m points
    for (int i=0; i<=n; i+=m)
      System.out.println(y1[i]);
  }

// Method to complete one Runge-Kutta step.

  public static double[] rungeKutta(double y[],
    double t, double dt) {
    int k = y.length;
    double k1[] = new double[k];
    double k2[] = new double[k];
    double k3[] = new double[k];
    double k4[] = new double[k];
    k1 = g(y, t);
    for (int i=0; i<k; ++i) k2[i] = y[i]+k1[i]/2;
    k2 = g(k2, t+dt/2);
    for (int i=0; i<k; ++i) k3[i] = y[i]+k2[i]/2;
    k3 = g(k3, t+dt/2);
    for (int i=0; i<k; ++i) k4[i] = y[i]+k3[i];
    k4 = g(k4, t+dt);
    for (int i=0; i<k; ++i)
      k1[i] = y[i]+dt*(k1[i]+2*(k2[i]+k3[i])+k4[i])/6;
    return k1;
  }

// Method to provide the generalized velocity vector.

  public static double[] g(double y[], double t) {
    int k = y.length;
    double v[] = new double[k];
    v[0] = y[1];
    v[1] = -Math.PI*Math.PI*(y[0]+1)/4;
    return v;
  }
}