Lab1-2

import com.quantifind.charts.Highcharts._

trait InitialVals {
  def x0: Double = 1
  def x1: Double = 2
  def y0: Double = 1
  def z0: Double = 0.5
  def makeXList(n: Int): List[Double] = {
    val h = (x1 - x0) / n
    val l = List(x0)
    def process(i: Int, acc: List[Double]): List[Double] = {
      if (i > 1) process(i - 1, acc :+ acc(n-i) + h) else acc
    }
    process(n, l)
  }
}

object DefEuler extends InitialVals {
  def perform(n: Int): (List[Double], List[Double]) = {
    val l     = (List[Double](), List[Double]())
    val h     = (x1 - x0) / n
    val xi    = x0 + (_: Int) * h
    val evalY = (y: Double, z: Double) => y + h * z
    val evalZ = (i: Int, y: Double, z: Double) =>
      z + h * ((1 / (4 * scala.math.pow(y, 3))) - (y / (4 * scala.math.pow(xi(i), 2))))

    def process(i: Int, y: Double, z: Double, acc: (List[Double], List[Double])): (List[Double], List[Double]) = {
      println(y + " " + z)
      val nextY = evalY(y, z)
      val nextZ = evalZ(i, y, z)
      if (i < n) process(i + 1, nextY, nextZ, (acc._1 :+ nextY, acc._2 :+  nextZ)) else acc
    }
    process(0, y0, z0, l)
  }
}

object Euler extends InitialVals {
  def perform(n: Int): (List[Double], List[Double]) = {
    val l     = (List[Double](1.0), List[Double](0.5))
    val h     = (x1 - x0) / n
    val xi    = x0 + (_: Int) * h

    def process(i: Int,
                y: Double,
                z: Double,
                acc: (List[Double], List[Double]))
              : (List[Double], List[Double]) = {
      def newtone(y: Double, z: Double, x: Double, count: Int): (Double, Double) = {
        println(count)

        val f1 = y - acc._1(i) - h * z
        val f2 = z - acc._2(i) - h * (1 / (4 * scala.math.pow(y, 3)) - (y / (4 * scala.math.pow(x, 2))))

        val df1dy = 1
        val df2dy = (h / 4) * ((3 / scala.math.pow(y, 4) + (1 / scala.math.pow(x, 2))))
        val df1dz = -h
        val df2dz = 1

        val detJ = df1dy * df2dz - df1dz * df2dy
        val detA1 = f1 * df2dz - df1dz * f2
        val detA2 = df1dy * f2 - f1 * df2dy

        val nextY = y - detA1 / detJ
        val nextZ = z - detA2 / detJ

        def newtoneCheck(a: Double, nextA: Double): Boolean = {
          val eps = 0.00000001
          scala.math.abs((nextA - a) / a ) < eps
        }

        if (newtoneCheck(y, nextY) && newtoneCheck(z, nextZ)) (nextY, nextZ) else newtone(nextY, nextZ, x, count + 1)
      }

      val nextY = newtone(y, z, xi(i), 0)._1
      val nextZ = newtone(y, z, xi(i), 0)._2

      if (i < n - 1) process(i + 1, nextY, nextZ, (acc._1 :+ nextY, acc._2 :+ nextZ)) else acc
    }
    process(0, y0, z0, l)
  }
}

object RungeFourth extends InitialVals {
  def perform(n: Int): (List[Double], List[Double]) = {
    val l     = (List[Double](), List[Double]())
    val h     = (x1 - x0) / n

    //def k1(x: Double, y: Double, f: (Double, Double) => Double): Double = f(x, y)
    //def k2(x: Double, y: Double, f: (Double, Double) => Double): Double = f(x + h / 2, y + h * k1(x, y, f) / 2)
    //def k3(x: Double, y: Double, f: (Double, Double) => Double): Double = f(x + h / 2, y + h * k2(x, y, f) / 2)
    //def k4(x: Double, y: Double, f: (Double, Double) => Double): Double = f(x + h, y + h * k3(x, y, f))
    def k1(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x, y)
    def k2(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x + h / 4, y + k1(x, y, f) / 4)
    def k3(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x + h / 2, y + k2(x, y, f) / 2)
    def k4(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x + h, y + k1(x, y, f) - 2 * k2(x, y, f) + k3(x, y, f))

    def fi(fPrevious: Double, x: Double, y: Double, f: (Double, Double) => Double): Double = {
      //fPrevious + ((h / 6) * (k1(x, y, f) + 2 * k2(x, y, f) + 2 * k3(x, y, f) + k4(x, y, f)))
      fPrevious + ((1.0 / 6) * (k1(x, y, f) + 4 * k3(x, y, f) + k4(x, y, f)))
    }

    def process(i: Int, y: Double, z: Double, acc: (List[Double], List[Double])): (List[Double], List[Double]) = {
      val xi    = x0 + i * h
      val evalY = (y: Double, z: Double) => z
      val evalZ = (y: Double, z: Double) => ((1 / (4 * scala.math.pow(y, 3))) - (y / (4 * scala.math.pow(xi, 2))))

      val nextY = fi(y, y, z, evalY)
      val nextZ = fi(z, y, z, evalZ)

      if (i < n) process(i + 1, nextY, nextZ, (acc._1 :+ nextY, acc._2 :+  nextZ)) else acc
    }
    process(0, y0, z0, l)
  }
}

object RungeFifth extends InitialVals {
  def perform(n: Int): (List[Double], List[Double]) = {
    val l     = (List[Double](), List[Double]())
    val h     = (x1 - x0) / n

    def k1(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x, y)
    def k2(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x + h / 2, y + k1(x, y, f) / 2)
    def k3(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x + h / 2, y + k1(x, y, f) / 4 + k2(x, y, f) / 4)
    def k4(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x + h, y - k2(x, y, f) + 2 * k3(x, y, f))
    def k5(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x + (2.0 / 3.0) * h, y + (7.0 / 27.0) * k1(x, y, f) + (10.0 / 27.0) * k2(x, y, f) + (1.0 / 27.0) * k4(x, y, f))
    def k6(x: Double, y: Double, f: (Double, Double) => Double): Double = h * f(x + (1.0 / 5.0) * h, y - (1.0 / 625.0) * (28 * k1(x, y, f) - 125 * k2(x, y, f) + 546 * k3(x, y, f) + 54 * k4(x, y, f) + 378 * k5(x, y, f)))

    def fi(fPrevious: Double, x: Double, y: Double, f: (Double, Double) => Double): Double = {
      fPrevious + ((1.0 / 336) * (14 * k1(x, y, f) + 35 * k4(x, y, f) + 162 * k5(x, y, f) + 125 * k6(x, y, f)))
    }

    def process(i: Int, y: Double, z: Double, acc: (List[Double], List[Double])): (List[Double], List[Double]) = {
      val xi    = x0 + i * h
      val evalY = (y: Double, z: Double) => z
      val evalZ = (y: Double, z: Double) => ((1 / (4 * scala.math.pow(y, 3))) - (y / (4 * scala.math.pow(xi, 2))))

      val nextY = fi(y, y, z, evalY)
      val nextZ = fi(z, y, z, evalZ)

      if (i < n) process(i + 1, nextY, nextZ, (acc._1 :+ nextY, acc._2 :+  nextZ)) else acc
    }
    process(0, y0, z0, l)
  }
}

object Comparator {
  def fetch(l1: List[Double], l2: List[Double], eps: Double): Unit = {
    def process(i: Int, l1: List[Double], l2: List[Double]): Unit = {
      if (l1.length > 0) {
        if (l1.head - l2.head >= eps) println(i + " " + l1.head + " " + l2.head)
        process(i + 1, l1.tail, l2.tail)
      }
    }
    process(0, l1, l2)
  }
}

object Main extends App with InitialVals {
  val num = args(0).toInt
  val eps = args(1).toDouble
  val result0 = RungeFourth.perform(num)
  //val result = DefEuler.perform(num)
  val result = RungeFifth.perform(num)
  val resultX = makeXList(num)
  val resultY = result._1
  val resultY0 = result0._1
  val resultZ0 = result0._2
  val resultZ = result._2

  println("Comparing Y's results of two methods")
  Comparator.fetch(result._1, result0._1, eps)
  println("Comparing Z's results of two methods")
  Comparator.fetch(result._2, result0._2, eps)

  line(resultX, resultY)
  hold
  line(resultX, resultY0)
  unhold
  line(resultX, resultZ)
  hold
  line(resultX, resultZ0)
}