object Scala0101 { def main(args:Array[String]) { println(3 + 5) println(3 - 5) println(3 * 5) println(Math.pow(3, 5)) println(5 / 3) println(5.0 / 3) println(5 / 3.0) println(5 % 3) printf("%d\n", 3 * 5) print(3 * 5 + "\n") println("%3d".format(3 * 5)) println("%23.20f".format(5 / 3.0)) } } /* Z:\>scala Scala0101.scala 8 -2 15 243.0 1 1.6666666666666667 1.6666666666666667 2 15 15 15 1.66666666666666670000 */
object Scala0102 { def main(args:Array[String]) { val i = 3 * 5; printf("3 * 5 = %d\n", i); } }
object Scala0103 { def main(args:Array[String]) { for (i <- 1 to 9) { print(i + ", ") } println() } }
object Scala0104 { def main(args:Array[String]) { for (i <- 1 to 9) { if (i % 3 == 0) { print(i + ", ") } } println() for (i <- 1 to 9 if i % 3 == 0) { print(i + ", ") } println() } }
object Scala0105 { def main(args:Array[String]) { var sm = 0 for (i <- 1 to 99) { if (i % 3 == 0) { sm += i } } println(sm) } }
(1 to 9)
(1 to 9).filter(n => n % 3 == 0)
(1 to 99).filter(n => n % 3 == 0).reduce((x, n) => x + n) (1 to 99).filter(n => n % 3 == 0).sum
// 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数 def sn(a:Int, lim:Int) = { val n = lim / a // 項数:n = 上限:lim / 公差:a val l = n * a // 末項:l = 項数:n * 公差:a (a + l) * n / 2 // 総和:sn = (初項:a + 末項:l) * 項数:n / 2 } // 3 の倍数の合計 sn(3,999)
val n = 10000
n * (n + 1) / 2
val n = 5000
n * (n + 1)
val n = 5000 (Math.pow(n, 2)).toInt
val n = 1000
n * (n + 1) * (2 * n + 1) / 6
val n = 100 (Math.pow(n, 2) * Math.pow((n + 1), 2) / 4).toInt
val n = 10 val a = 2 val r = 3 a * (Math.pow(r, n).toInt - 1) / (r -1)
(0 to 9) (0 to 9).map(n => n * 3 + 5) (0l to 9l).map(n => n * 3 + 5).product
// 等差数列の積 def prod(m: Long, d: Int, n: Int): Long = { n match { case 0 => 1 case _ => m * prod(m + d, d, n - 1) } } // 初項 5, 公差 3, 項数 10 の数列の積 prod(5, 3, 10)
// 階乗を求める関数 def Fact(n: Int): Int = { n match { case 0 => 1 case _ => n * Fact(n - 1) } } // 10の階乗 Fact(10)
// 下降階乗冪 def FallingFact(x: Int, n: Int): Int = { n match { case 1 => x case _ => x * FallingFact(x - 1, n - 1) } } // 10 から 6 までの 総乗 FallingFact(10, 5)
// 上昇階乗冪 def RisingFact(x: Int, n: Int): Int = { n match { case 1 => x case _ => x * RisingFact(x + 1, n - 1) } } // 10 から 14 までの 総乗 RisingFact(10, 5)
// 階乗 def Fact(n: Int): Int = { n match { case 0 => 1 case _ => n * Fact(n - 1) } } // 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) val n = 10 val r = 5 Fact(n) / Fact(n - r) // 下降階乗冪 def FallingFact(x: Int, n: Int): Int = { n match { case 1 => x case _ => x * FallingFact(x - 1, n - 1) } } // 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) val n = 10 val r = 5 FallingFact(n, r)
// 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) val n = 10 val r = 5 Math.pow(n, r)
// 組合せ def Comb(n: Int, r: Int): Int = { (n, r) match { case (_, 0) => 1 case (_, 1) => n case (_, _) if n == r => 1 case (_, _) => Comb(n - 1, r - 1) + Comb(n - 1, r) } } // 異なる 10 個のものから 5 個取ってできる組合せの総数 val n = 10 val r = 5 Comb(n, r)
// 組合せ def Comb(n: Int, r: Int): Int = { (n, r) match { case (_, 0) => 1 case (_, 1) => n case (_, _) if n == r => 1 case (_, _) => Comb(n - 1, r - 1) + Comb(n - 1, r) } } // 異なる 10 個のものから重複を許して 5 個とる組合せの総数 val n = 10 val r = 5 Comb(n + r - 1, r)
// 自作の正弦関数 def mySin(x:Double, n:Int, nega:Boolean, numerator:Double, denominator:Double, y:Double):Double = { val m = 2 * n val denom = denominator * (m + 1) * m val num = numerator * x * x val a = num / denom // 十分な精度になったら処理を抜ける if (a <= 0.00000000001) y else y + mySin(x, n + 1, !nega, num, denom, if (nega) a else -a) } (0 to 24). map(_ * 15). // 0° 〜 360° のうち 30°, 45° の倍数 filter(n => (n % 30 == 0) || (n % 45 == 0)). // ラジアン map(degree => (degree, Math.toRadians(degree))). // 標準の正弦関数 と 自作の正弦関数 map(t => (t._1, Math.sin(t._2), mySin(t._2, 1, false, t._2, 1.0, t._2))). // 標準関数との差異 foreach {t => t match { case (degree, d1, d2) => println("%3d : %13.10f - %13.10f = %13.10f".format(degree, d1, d2, d1 - d2)) } } (0 to 24). map(_ * 15). filter(n => (n % 30 == 0) || (n % 45 == 0)). foreach { degree => val radian = Math.toRadians(degree) // 自作の正弦関数 val d1 = mySin(radian, 1, false, radian, 1.0, radian) // 標準の正弦関数 val d2 = Math.sin(radian) // 標準関数との差異 println("%3d : %13.10f - %13.10f = %13.10f".format(degree, d1, d2, d1 - d2)); }
// 自作の余弦関数 def myCos(x:Double, n:Int, nega:Boolean, numerator:Double, denominator:Double, y:Double):Double = { val m = 2 * n val denom = denominator * m * (m - 1) val num = numerator * x * x val a = num / denom // 十分な精度になったら処理を抜ける if (a <= 0.00000000001) y else y + myCos(x, n + 1, !nega, num, denom, if (nega) a else -a) } (0 to 24). map(_ * 15). filter(n => (n % 30 == 0) || (n % 45 == 0)). foreach { degree => val radian = Math.toRadians(degree) // 自作の余弦関数 val d1 = myCos(radian, 1, false, 1.0, 1.0, 1.0) // 標準の余弦関数 val d2 = Math.cos(radian) // 標準関数との差異 System.out.println("%3d : %13.10f - %13.10f = %13.10f".format(degree, d1, d2, d1 - d2)); }
// 自作の正接関数 def myTan(x:Double, x2:Double, n:Int, t:Double):Double = { val denom = x2 / (n - t) val num = n - 2 if (num <= 1) x / (1 - denom) else myTan(x, x2, num, denom) } (0 to 12). map(_ * 15). filter(n => n % 180 != 0). map(_ - 90). foreach { degree => val radian = Math.toRadians(degree) val x2 = radian * radian // 自作の正接関数 val d1 = myTan(radian, x2, 15, 0.0) // 15:必要な精度が得られる十分大きな奇数 // 標準の正接関数 val d2 = Math.tan(radian) // 標準関数との差異 System.out.println("%3d : %13.10f - %13.10f = %13.10f".format(degree, d1, d2, d1 - d2)); }
// 自作の指数関数 def myExp(x:Double, n:Int, numerator:Double, denominator:Double, y:Double):Double = { val denom = denominator * n val num = numerator * x val a = num / denom // 十分な精度になったら処理を抜ける if (Math.abs(a) <= 0.00000000001) y else y + myExp(x, n + 1, num, denom, a) } (0 to 20). map(n => (n - 10) / 4.0). foreach { x => // 標準の指数関数 val d1 = Math.exp(x) // 自作の指数関数 val d2 = myExp(x, 1, 1.0, 1.0, 1.0) // 標準関数との差異 System.out.println("%5.2f : %13.10f - %13.10f = %13.10f".format(x, d1, d2, d1 - d2)); }
// 自作の指数関数 def myExp(x:Double, x2:Double, n:Int, t:Double):Double = { val denom = x2 / (n + t) val num = n - 4 if (num < 6) 1 + ((2 * x) / (2 - x + denom)) else myExp(x, x2, num, denom) } (0 to 20). map(n => (n - 10) / 4.0). foreach { x => // 標準の指数関数 val d1 = Math.exp(x) // 自作の指数関数 val x2 = x * x val d2 = myExp(x, x2, 30, 0.0) // 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる // 標準関数との差異 System.out.println("%5.2f : %13.10f - %13.10f = %13.10f".format(x, d1, d2, d1 - d2)) }
// 自作の対数関数 def myLog(x2:Double, numerator:Double, denominator:Double, y:Double):Double = { val denom = denominator + 2 val num = numerator * x2 * x2 val a = num / denom // 十分な精度になったら処理を抜ける if (Math.abs(a) <= 0.00000000001) y else y + myLog(x2, num, denom, a) } (1 to 20). map(_ / 5.0). foreach { x => // 標準の対数関数 val d1 = Math.log(x) // 自作の対数関数 val x2 = (x - 1) / (x + 1) val d2 = 2 * myLog(x2, x2, 1.0, x2) // 標準関数との差異 System.out.println("%5.2f : %13.10f - %13.10f = %13.10f".format(x, d1, d2, d1 - d2)) }
// 自作の対数関数 def myLog(x:Double, n:Int, t:Double):Double = { var n2 = n var x2 = x if (n > 3) { if (n % 2 == 0) n2 = 2 x2 = x * (n / 2) } val t2 = x2 / (n2 + t) if (n <= 2) x / (1 + t2) else myLog(x, n - 1, t2) } (1 to 20). map(_ / 5.0). foreach { x => // 標準の対数関数 val d1 = Math.log(x) // 自作の対数関数 val d2 = myLog(x - 1, 27, 0.0) // 27:必要な精度が得られる十分大きな奇数 // 標準関数との差異 System.out.println("%5.2f : %13.10f - %13.10f = %13.10f".format(x, d1, d2, d1 - d2)) }
// 自作の双曲線正弦関数 def mySinh(x:Double, n:Int, numerator:Double, denominator:Double, y:Double):Double = { val m = 2 * n val denom = denominator * (m + 1) * m val num = numerator * x * x val a = num / denom // 十分な精度になったら処理を抜ける if (Math.abs(a) <= 0.00000000001) y else y + mySinh(x, n + 1, num, denom, a) } (0 to 20). map(_ - 10). foreach { x => // 自作の双曲線正弦関数 val d1 = mySinh(x, 1, x, 1.0, x) // 標準の双曲線正弦関数 val d2 = Math.sinh(x) // 標準関数との差異 System.out.println("%3d : %17.10f - %17.10f = %13.10f".format(x, d1, d2, d1 - d2)); }
// 自作の双曲線余弦関数 def myCosh(x:Double, n:Int, numerator:Double, denominator:Double, y:Double):Double = { val m = 2 * n val denom = denominator * m * (m - 1) val num = numerator * x * x val a = num / denom // 十分な精度になったら処理を抜ける if (Math.abs(a) <= 0.00000000001) y else y + myCosh(x, n + 1, num, denom, a) } (0 to 20). map(_ - 10). foreach { x => // 自作の双曲線余弦関数 val d1 = myCosh(x, 1, 1.0, 1.0, 1.0) // 標準の双曲線余弦関数 val d2 = Math.cosh(x) // 標準関数との差異 System.out.println("%3d : %17.10f - %17.10f = %13.10f".format(x, d1, d2, d1 - d2)); }
object Scala0601 { def main(args: Array[String]) { val a:Double = 0 val b:Double = 1 // 台形則で積分 var n:Int = 2 for (j <- 1 to 10) { var h:Double = (b - a) / n var s:Double = 0 var x:Double = a for (i <- 1 to (n - 1)) { x += h s += f(x) } s = h * ((f(a) + f(b)) / 2 + s) n *= 2 // 結果を π と比較 println("%3d : %13.10f, %13.10f".format(j, s, s - Math.PI)) } } def f(x:Double):Double = { 4 / (1 + x * x) } }
object Scala0602 { def main(args: Array[String]) { val a:Double = 0 val b:Double = 1 // 中点則で積分 var n:Int = 2 for (j <- 1 to 10) { var h:Double = (b - a) / n var s:Double = 0 var x:Double = a + (h / 2) for (i <- 1 to n) { s += f(x) x += h } s *= h n *= 2 // 結果を π と比較 println("%3d : %13.10f, %13.10f".format(j, s, s - Math.PI)) } } def f(x:Double):Double = { 4 / (1 + x * x) } }
object Scala0603 { def main(args: Array[String]) { val a:Double = 0 val b:Double = 1 // Simpson則で積分 var n:Int = 2 for (j <- 1 to 5) { var h:Double = (b - a) / n var s2:Double = 0 var s4:Double = 0 var x:Double = a + h for (i <- 1 to (n / 2)) { s4 += f(x) x += h s2 += f(x) x += h } s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 var s:Double = (s2 + s4) * h / 3 n *= 2 // 結果を π と比較 println("%3d : %13.10f, %13.10f".format(j, s, s - Math.PI)) } } def f(x:Double):Double = { 4 / (1 + x * x) } }
object Scala0604 { def main(args: Array[String]) { val a:Double = 0 val b:Double = 1 val t = Array.ofDim[Double](7, 7) // 台形則で積分 var n:Int = 2 for (i <- 1 to 6) { var h:Double = (b - a) / n var s:Double = 0 var x:Double = a for (j <- 1 to (n - 1)) { x += h s += f(x) } // 結果を保存 t(i)(1) = h * ((f(a) + f(b)) / 2 + s) n *= 2 } // Richardsonの補外法 n = 4 for (j <- 2 to 6) { for (i <- j to 6) { t(i)(j) = t(i)(j - 1) + (t(i)(j - 1) - t(i - 1)(j - 1)) / (n - 1) if (i == j) { // 結果を π と比較 println("%3d : %13.10f, %13.10f".format(j, t(i)(j), t(i)(j) - Math.PI)) } } n *= 4 } } def f(x:Double):Double = { 4 / (1 + x * x) } }
object Scala0701 { // データ点の数 - 1 val N = 6 def main(args: Array[String]) { // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット val x = (0 to N).map(_ * 1.5 - 4.5) val y = x.map(f) // 0.5刻みで 与えられていない値を補間 val d1 = (0 to 18).map(_ * 0.5 - 4.5) val d2 = d1.map(f) val d3 = d1.map(lagrange(_, x, y)) (d1 zip d2 zip d3).foreach { case ((d1, d2), d3) => println("%5.2f\t%8.5f\t%8.5f\t%8.5f".format(d1, d2, d3, d2 - d3)) } } // 元の関数 def f(x:Double) = { x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // Lagrange (ラグランジュ) 補間 def lagrange(d:Double, x:IndexedSeq[Double], y:IndexedSeq[Double]) = { var sum_list = List[Double](0) for (i <- 0 to N) { var prod_list = List(y(i)) for (j <- 0 to N) { if (i != j) { prod_list = ((d - x(j)) / (x(i) - x(j)))::prod_list } } sum_list = prod_list.product::sum_list } sum_list.sum } }
object Scala0702 { // データ点の数 - 1 val N = 6 def main(args: Array[String]) { // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット val x = (0 to N).map(_ * 1.5 - 4.5) val y = x.map(f) // 0.5刻みで 与えられていない値を補間 val d1 = (0 to 18).map(_ * 0.5 - 4.5) val d2 = d1.map(f) val d3 = d1.map(neville(_, x, y)) (d1 zip d2 zip d3).foreach { case ((d1, d2), d3) => println("%5.2f\t%8.5f\t%8.5f\t%8.5f".format(d1, d2, d3, d2 - d3)) } } // 元の関数 def f(x:Double) = { x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // Neville (ネヴィル) 補間 def neville(d:Double, x:IndexedSeq[Double], y:IndexedSeq[Double]) = { val w = Array.ofDim[Double](N + 1, N + 1) for (i <- 0 to N) w(0)(i) = y(i) for (j <- 1 to N) { for (i <- 0 to N - j) w(j)(i) = w(j-1)(i+1) + (w(j-1)(i+1) - w(j-1)(i)) * (d - x(i+j)) / (x(i+j) - x(i)) } w(N)(0) } }
object Scala0703 { // データ点の数 - 1 val N = 6 def main(args: Array[String]) { // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット val x = (0 to N).map(_ * 1.5 - 4.5) val y = x.map(f) // 差分商の表を作る val d = Array.ofDim[Double](N + 1, N + 1) for (j <- 0 to N) d(0)(j) = y(j) for (i <- 1 to N) { for (j <- 0 to N - i) d(i)(j) = (d(i-1)(j+1) - d(i-1)(j)) / (x(j+i) - x(j)) } // n階差分商 val a = (0 to N).map(d(_)(0)) // 0.5刻みで 与えられていない値を補間 val d1 = (0 to 18).map(_ * 0.5 - 4.5) val d2 = d1.map(f) val d3 = d1.map(newton(_, x, a)) (d1 zip d2 zip d3).foreach { case ((d1, d2), d3) => println("%5.2f\t%8.5f\t%8.5f\t%8.5f".format(d1, d2, d3, d2 - d3)) } } // 元の関数 def f(x:Double) = { x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // Newton (ニュートン) 補間 def newton(d:Double, x:IndexedSeq[Double], a:IndexedSeq[Double]) = { var sum_list = List(a(0)) for (i <- 1 to N) { var prod_list = List(a(i)) for (j <- 0 to i - 1) { prod_list = (d - x(j))::prod_list } sum_list = (prod_list.product)::sum_list } sum_list.sum } }
object Scala0704 { // データ点の数 - 1 val N = 6 val Nx2 = 13 def main(args: Array[String]) { // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット val x = (0 to N).map(_ * 1.5 - 4.5) val y = x.map(f) val yd = x.map(fd) // 差分商の表を作る val z = (0 to Nx2).map(_ / 2).map(x(_)) val d = Array.ofDim[Double](Nx2 + 1, Nx2 + 1) for (i <- 0 to Nx2) d(0)(i) = y(i / 2) for (i <- 1 to Nx2) { for (j <- 0 to Nx2 - i) if (i == 1 && j % 2 == 0) d(i)(j) = yd(j / 2) else d(i)(j) = (d(i-1)(j+1) - d(i-1)(j)) / (z(j+i) - z(j)) } // n階差分商 val a = (0 to Nx2).map(d(_)(0)) // 0.5刻みで 与えられていない値を補間 val d1 = (0 to 18).map(_ * 0.5 - 4.5) val d2 = d1.map(f) val d3 = d1.map(hermite(_, z, a)) (d1 zip d2 zip d3).foreach { case ((d1, d2), d3) => println("%5.2f\t%8.5f\t%8.5f\t%8.5f".format(d1, d2, d3, d2 - d3)) } } // 元の関数 def f(x:Double) = { x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // 導関数 def fd(x:Double) = { 1 - Math.pow(x,2) / 2 + Math.pow(x,4) / (4 * 3 * 2) } // Hermite (エルミート) 補間 def hermite(d:Double, z:IndexedSeq[Double], a:IndexedSeq[Double]) = { var sum_list = List(a(0)) for (i <- 1 to Nx2) { var prod_list = List(a(i)) for (j <- 0 to i - 1) { prod_list = (d - z(j))::prod_list } sum_list = (prod_list.product)::sum_list } sum_list.sum } }
object Scala0705 { // データ点の数 - 1 val N = 6 def main(args: Array[String]) { // 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット val x = (0 to N).map(_ * 1.5 - 4.5) val y = x.map(f) // 3項方程式の係数の表を作る val a = Array.ofDim[Double](N) val b = Array.ofDim[Double](N) val c = Array.ofDim[Double](N) val d = Array.ofDim[Double](N) for (i <- 1 to N - 1) { a(i) = x(i) - x(i-1) b(i) = 2.0 * (x(i+1) - x(i-1)) c(i) = x(i+1) - x(i) d(i) = 6.0 * ((y(i+1) - y(i)) / (x(i+1) - x(i)) - (y(i) - y(i-1)) / (x(i) - x(i-1))) } // 3項方程式を解く (ト−マス法) val g = Array.ofDim[Double](N) val s = Array.ofDim[Double](N) g(1) = b(1) s(1) = d(1) for (i <- 2 to N - 1) { g(i) = b(i) - a(i) * c(i-1) / g(i-1) s(i) = d(i) - a(i) * s(i-1) / g(i-1) } val z = Array.ofDim[Double](N + 1) z(0) = 0 z(N) = 0 z(N-1) = s(N-1) / g(N-1) for (i <- N - 2 to 1 by -1) z(i) = (s(i) - c(i) * z(i+1)) / g(i) // 0.5刻みで 与えられていない値を補間 val d1 = (0 to 18).map(_ * 0.5 - 4.5) val d2 = d1.map(f) val d3 = d1.map(spline(_, x, y, z)) (d1 zip d2 zip d3).foreach { case ((d1, d2), d3) => println("%5.2f\t%8.5f\t%8.5f\t%8.5f".format(d1, d2, d3, d2 - d3)) } } // 元の関数 def f(x:Double) = { x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2) } // Spline (スプライン) 補間 def spline(d:Double, x:IndexedSeq[Double], y:IndexedSeq[Double], z:IndexedSeq[Double]) = { // 補間関数値がどの区間にあるか var k = -1 for (i <- N to 1 by -1) { if (d <= x(i)) k = i - 1 } if (k < 0) k = N val d1 = x(k+1) - d val d2 = d - x(k) val d3 = x(k+1) - x(k) (z(k) * Math.pow(d1,3) + z(k+1) * Math.pow(d2,3)) / (6.0 * d3) + (y(k) / d3 - z(k) * d3 / 6.0) * d1 + (y(k+1) / d3 - z(k+1) * d3 / 6.0) * d2 } }
object Scala0801 { // 重力加速度 val g = -9.8 // 空気抵抗係数 val k = -0.01 // 時間間隔(秒) val h = 0.01 def main(args: Array[String]) { // 角度 val degree = 45 val radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 val v = 250 * 1000 / 3600 // 水平方向の速度 val vx = v * Math.cos(radian) // 鉛直方向の速度 val vy = v * Math.sin(radian) // 位置 val x = 0.0 val y = 0.0 // Euler法 euler(1, vx, vy, x, y) } def euler(i:Int, vx:Double, vy:Double, x:Double, y:Double):Unit = { // 経過秒数 val t = i * h // 位置 val wx = x + h * vx val wy = y + h * vy println("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f".format(t, vx, vy, wx, wy)) // 速度 val wvx = vx + h * fx(vx, vy) val wvy = vy + h * fy(vx, vy) if (wy >= 0.0) euler(i+1, wvx, wvy, wx, wy) else () } // 空気抵抗による水平方向の減速分 def fx(vx:Double, vy:Double) = { k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 def fy(vx:Double, vy:Double) = { g + (k * Math.sqrt(vx * vx + vy * vy) * vy) } }
object Scala0802 { // 重力加速度 val g = -9.8 // 空気抵抗係数 val k = -0.01 // 時間間隔(秒) val h = 0.01 def main(args: Array[String]) { // 角度 val degree = 45 val radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 val v = 250 * 1000 / 3600 // 水平方向の速度 val vx = v * Math.cos(radian) // 鉛直方向の速度 val vy = v * Math.sin(radian) // 位置 val x = 0.0 val y = 0.0 // Heun法 heun(1, vx, vy, x, y) } // Heun法 def heun(i:Int, vx:Double, vy:Double, x:Double, y:Double):Unit = { // 経過秒数 val t = i * h // 位置・速度 val wx2 = x + h * vx val wy2 = y + h * vy val wvx2 = vx + h * fx(vx, vy) val wvy2 = vy + h * fy(vx, vy) val wx = x + h * ( vx + wvx2 ) / 2 val wy = y + h * ( vy + wvy2 ) / 2 val wvx = vx + h * (fx(vx, vy) + fx(wvx2, wvy2)) / 2 val wvy = vy + h * (fy(vx, vy) + fy(wvx2, wvy2)) / 2 println("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f".format(t, wvx, wvy, wx, wy)) if (wy >= 0.0) heun(i+1, wvx, wvy, wx, wy) else () } // 空気抵抗による水平方向の減速分 def fx(vx:Double, vy:Double) = { k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 def fy(vx:Double, vy:Double) = { g + (k * Math.sqrt(vx * vx + vy * vy) * vy) } }
object Scala0803 { // 重力加速度 val g = -9.8 // 空気抵抗係数 val k = -0.01 // 時間間隔(秒) val h = 0.01 def main(args: Array[String]) { // 角度 val degree = 45 val radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 val v = 250 * 1000 / 3600 // 水平方向の速度 val vx = v * Math.cos(radian) // 鉛直方向の速度 val vy = v * Math.sin(radian) // 位置 val x = 0.0 val y = 0.0 // 中点法 midpoint(1, vx, vy, x, y) } def midpoint(i:Int, vx:Double, vy:Double, x:Double, y:Double):Unit = { // 経過秒数 val t = i * h // 位置・速度 val wvx1 = h * fx(vx, vy) val wvy1 = h * fy(vx, vy) val wvx2 = vx + wvx1 / 2 val wvy2 = vy + wvy1 / 2 val wvx = vx + h * fx(wvx2, wvy2) val wvy = vy + h * fy(wvx2, wvy2) val wx = x + h * wvx2 val wy = y + h * wvy2 println("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f".format(t, wvx, wvy, wx, wy)) if (wy >= 0.0) midpoint(i+1, wvx, wvy, wx, wy) else () } // 空気抵抗による水平方向の減速分 def fx(vx:Double, vy:Double) = { k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 def fy(vx:Double, vy:Double) = { g + (k * Math.sqrt(vx * vx + vy * vy) * vy) } }
object Scala0804 { // 重力加速度 val g = -9.8 // 空気抵抗係数 val k = -0.01 // 時間間隔(秒) val h = 0.01 def main(args: Array[String]) { // 角度 val degree = 45 val radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 val v = 250 * 1000 / 3600 // 水平方向の速度 val vx = v * Math.cos(radian) // 鉛直方向の速度 val vy = v * Math.sin(radian) // 位置 val x = 0.0 val y = 0.0 // Runge-Kutta法 rungekutta(1, vx, vy, x, y) } def rungekutta(i:Int, vx:Double, vy:Double, x:Double, y:Double):Unit = { // 経過秒数 val t = i * h // 位置・速度 val wx1 = h * vx val wy1 = h * vy val wvx1 = h * fx(vx, vy) val wvy1 = h * fy(vx, vy) val wvx5 = vx + wvx1 / 2 val wvy5 = vy + wvy1 / 2 val wx2 = h * wvx5 val wy2 = h * wvy5 val wvx2 = h * fx(wvx5, wvy5) val wvy2 = h * fy(wvx5, wvy5) val wvx6 = vx + wvx2 / 2 val wvy6 = vy + wvy2 / 2 val wx3 = h * wvx6 val wy3 = h * wvy6 val wvx3 = h * fx(wvx6, wvy6) val wvy3 = h * fy(wvx6, wvy6) val wvx7 = vx + wvx3 val wvy7 = vy + wvy3 val wx4 = h * wvx7 val wy4 = h * wvy7 val wvx4 = h * fx(wvx7, wvy7) val wvy4 = h * fy(wvx7, wvy7) val wx = x + ( wx1 + wx2 * 2 + wx3 * 2 + wx4) / 6 val wy = y + ( wy1 + wy2 * 2 + wy3 * 2 + wy4) / 6 val wvx = vx + (wvx1 + wvx2 * 2 + wvx3 * 2 + wvx4) / 6 val wvy = vy + (wvy1 + wvy2 * 2 + wvy3 * 2 + wvy4) / 6 println("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f".format(t, wvx, wvy, wx, wy)) if (wy >= 0.0) rungekutta(i+1, wvx, wvy, wx, wy) else () } // 空気抵抗による水平方向の減速分 def fx(vx:Double, vy:Double) = { k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 def fy(vx:Double, vy:Double) = { g + (k * Math.sqrt(vx * vx + vy * vy) * vy) } }
object Scala0804 { // 重力加速度 val g = -9.8 // 空気抵抗係数 val k = -0.01 // 時間間隔(秒) val h = 0.01 def main(args: Array[String]) { // 角度 val degree = 45 val radian = degree * Math.PI / 180.0 // 初速 250 km/h -> 秒速に変換 val v = 250 * 1000 / 3600 // 水平方向の速度 val vx = v * Math.cos(radian) // 鉛直方向の速度 val vy = v * Math.sin(radian) // 位置 val x = 0.0 val y = 0.0 // Runge-Kutta-Gill法 rungekuttagill(1, vx, vy, x, y) } def rungekuttagill(i:Int, vx:Double, vy:Double, x:Double, y:Double):Unit = { // 経過秒数 val t = i * h // 位置・速度 val wx1 = h * vx val wy1 = h * vy val wvx1 = h * fx(vx, vy) val wvy1 = h * fy(vx, vy) val wvx5 = vx + wvx1 / 2 val wvy5 = vy + wvy1 / 2 val wx2 = h * wvx5 val wy2 = h * wvy5 val wvx2 = h * fx(wvx5, wvy5) val wvy2 = h * fy(wvx5, wvy5) val wvx6 = vx + wvx1 * ((Math.sqrt(2.0) - 1) / 2) + wvx2 * (1 - 1 / Math.sqrt(2.0)) val wvy6 = vy + wvy1 * ((Math.sqrt(2.0) - 1) / 2) + wvy2 * (1 - 1 / Math.sqrt(2.0)) val wx3 = h * wvx6 val wy3 = h * wvy6 val wvx3 = h * fx(wvx6, wvy6) val wvy3 = h * fy(wvx6, wvy6) val wvx7 = vx - wvx2 / Math.sqrt(2.0) + wvx3 * (1 + 1 / Math.sqrt(2.0)) val wvy7 = vy - wvy2 / Math.sqrt(2.0) + wvy3 * (1 + 1 / Math.sqrt(2.0)) val wx4 = h * wvx7 val wy4 = h * wvy7 val wvx4 = h * fx(wvx7, wvy7) val wvy4 = h * fy(wvx7, wvy7) val wx = x + ( wx1 + wx2 * (2 - Math.sqrt(2.0)) + wx3 * (2 + Math.sqrt(2.0)) + wx4) / 6 val wy = y + ( wy1 + wy2 * (2 - Math.sqrt(2.0)) + wy3 * (2 + Math.sqrt(2.0)) + wy4) / 6 val wvx = vx + (wvx1 + wvx2 * (2 - Math.sqrt(2.0)) + wvx3 * (2 + Math.sqrt(2.0)) + wvx4) / 6 val wvy = vy + (wvy1 + wvy2 * (2 - Math.sqrt(2.0)) + wvy3 * (2 + Math.sqrt(2.0)) + wvy4) / 6 println("%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f".format(t, wvx, wvy, wx, wy)) if (wy >= 0.0) rungekuttagill(i+1, wvx, wvy, wx, wy) else () } // 空気抵抗による水平方向の減速分 def fx(vx:Double, vy:Double) = { k * Math.sqrt(vx * vx + vy * vy) * vx } // 重力と空気抵抗による鉛直方向の減速分 def fy(vx:Double, vy:Double) = { g + (k * Math.sqrt(vx * vx + vy * vy) * vy) } }
object Scala0901 { def main(args: Array[String]) { val a = 1.0 val b = 2.0 println("%12.10f".format(bisection(a, b))) } def bisection(a:Double, b:Double):Double = { // 区間 (a, b) の中点 c = (a + b) / 2 val c = (a + b) / 2 println("%12.10f\t%13.10f".format(c, c - Math.sqrt(2))) val fc = f(c) if (Math.abs(fc) < 0.0000000001) c else { if (fc < 0) { // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 bisection(c, b) } else { // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 bisection(a, c) } } } def f(x:Double) = { x * x - 2 } }
object Scala0902 { def main(args: Array[String]) { val a = 1.0 val b = 2.0 println("%12.10f".format(falseposition(a, b))) } def falseposition(a:Double, b:Double):Double = { // 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 val c = (a * f(b) - b * f(a)) / (f(b) - f(a)) println("%12.10f\t%13.10f".format(c, c - Math.sqrt(2))) val fc = f(c) if (Math.abs(fc) < 0.0000000001) c else { if (fc < 0) { // f(c) < 0 であれば, 解は区間 (c, b) の中に存在 falseposition(c, b) } else { // f(c) > 0 であれば, 解は区間 (a, c) の中に存在 falseposition(a, c) } } } def f(x:Double) = { x * x - 2 } }
object Scala0903 { def main(args: Array[String]) { val x = 1.0 println("%12.10f".format(iterative(x))) } def iterative(x0:Double):Double = { val x1 = g(x0) println("%12.10f\t%13.10f".format(x1, x1 - Math.sqrt(2))) if (Math.abs(x1 - x0) < 0.0000000001) x1 else iterative(x1) } def g(x:Double) = { (x / 2) + (1 / x) } }
object Scala0904 { def main(args: Array[String]) { val x = 2.0 println("%12.10f".format(newton(x))) } def newton(x0:Double):Double = { val x1 = x0 - (f0(x0) / f1(x0)) println("%12.10f\t%13.10f".format(x1, x1 - Math.sqrt(2))) if (Math.abs(x1 - x0) < 0.0000000001) x1 else newton(x1) } def f0(x:Double) = { x * x - 2 } def f1(x:Double) = { 2 * x } }
object Scala0905 { def main(args: Array[String]) { val x = 2.0 println("%12.10f".format(bailey(x))) } def bailey(x0:Double):Double = { val x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0))))) println("%12.10f\t%13.10f".format(x1, x1 - Math.sqrt(2))) if (Math.abs(x1 - x0) < 0.0000000001) x1 else bailey(x1) } def f0(x:Double) = { x * x - 2 } def f1(x:Double) = { 2 * x } def f2(x:Double) = { 2 } }
object Scala0906 { def main(args: Array[String]) { val x0 = 1.0 val x1 = 2.0 println("%12.10f".format(secant(x0, x1))) } def secant(x0:Double, x1:Double):Double = { val x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) println("%12.10f\t%13.10f".format(x2, x2 - Math.sqrt(2))) if (Math.abs(x2 - x1) < 0.0000000001) x2 else secant(x1, x2) } def f(x:Double) = { x * x - 2 } }
object Scala1001 { val N = 3 def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(9,2,1,1),Array(2,8,-2,1),Array(-1,-2,7,-2),Array(1,-1,-2,6)) var b:Array[Double] = Array(20,16,8,17) var c:Array[Double] = Array(0,0,0,0) // ヤコビの反復法 jacobi(a,b,c) println("X") disp_vector(c) } // ヤコビの反復法 def jacobi(a:Array[Array[Double]], b:Array[Double], x0:Array[Double]) = { var finish:Boolean = false while (!finish) { var x1:Array[Double] = Array(0,0,0,0) finish = true for (i <- 0 to N) { x1(i) = 0 for (j <- 0 to N) if (j != i) x1(i) += a(i)(j) * x0(j) x1(i) = (b(i) - x1(i)) / a(i)(i) if (Math.abs(x1(i) - x0(i)) > 0.0000000001) finish = false } for (i <- 0 to N) x0(i) = x1(i) if (!finish) disp_vector(x0) } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } }
object Scala1002 { val N = 3 def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(9,2,1,1),Array(2,8,-2,1),Array(-1,-2,7,-2),Array(1,-1,-2,6)) var b:Array[Double] = Array(20,16,8,17) var c:Array[Double] = Array(0,0,0,0) // ガウス・ザイデル法 gauss(a,b,c) println("X") disp_vector(c) } // ガウス・ザイデル法 def gauss(a:Array[Array[Double]], b:Array[Double], x0:Array[Double]) = { var finish:Boolean = false while (!finish) { finish = true for (i <- 0 to N) { var x1:Double = 0 for (j <- 0 to N) if (j != i) x1 += a(i)(j) * x0(j) x1 = (b(i) - x1) / a(i)(i) if (Math.abs(x1 - x0(i)) > 0.0000000001) finish = false x0(i) = x1 } if (!finish) disp_vector(x0) } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } }
object Scala1003 { val N = 3 def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1)) var b:Array[Double] = Array(8,17,20,16) // ピボット選択 pivoting(a,b) println("pivoting") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("forward elimination") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 後退代入 backward_substitution(a,b) println("X") disp_vector(b) } // 前進消去 def forward_elimination(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N - 1) { for (row <- pivot + 1 to N) { var s:Double = a(row)(pivot) / a(pivot)(pivot) for (col <- pivot to N) a(row)(col) -= a(pivot)(col) * s b(row) -= b(pivot) * s } } } // 後退代入 def backward_substitution(a:Array[Array[Double]], b:Array[Double]) { for (row <- (0 to N).reverse) { for (col <- (row + 1 to N).reverse) b(row) -= a(row)(col) * b(col) b(row) /= a(row)(row) } } // ピボット選択 def pivoting(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) { // 各列で 一番値が大きい行を 探す var max_row:Integer = pivot var max_val:Double = 0 for (row <- pivot to N) { if (Math.abs(a(row)(pivot)) > max_val) { // 一番値が大きい行 max_val = Math.abs(a(row)(pivot)) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp:Double = 0 for (col <- 0 to N) { tmp = a(max_row)(col) a(max_row)(col) = a(pivot)(col) a(pivot)(col) = tmp } tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp } } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 2次元配列を表示 def disp_matrix(matrix:Array[Array[Double]]) = { matrix.foreach { row => row.foreach { col => print("%14.10f\t".format(col)) } println() } } }
object Scala1004 { val N = 3 def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1)) var b:Array[Double] = Array(8,17,20,16) // ピボット選択 pivoting(a,b) println("pivoting") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("forward elimination") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 後退代入 backward_substitution(a,b) println("X") disp_vector(b) } // 前進消去 def forward_elimination(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) { for (row <- 0 to N) { if (row != pivot) { var s:Double = a(row)(pivot) / a(pivot)(pivot) for (col <- pivot to N) a(row)(col) -= a(pivot)(col) * s b(row) -= b(pivot) * s } } } } // 後退代入 def backward_substitution(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) b(pivot) /= a(pivot)(pivot) } // ピボット選択 def pivoting(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) { // 各列で 一番値が大きい行を 探す var max_row:Integer = pivot var max_val:Double = 0 for (row <- pivot to N) { if (Math.abs(a(row)(pivot)) > max_val) { // 一番値が大きい行 max_val = Math.abs(a(row)(pivot)) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp:Double = 0 for (col <- 0 to N) { tmp = a(max_row)(col) a(max_row)(col) = a(pivot)(col) a(pivot)(col) = tmp } tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp } } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 2次元配列を表示 def disp_matrix(matrix:Array[Array[Double]]) = { matrix.foreach { row => row.foreach { col => print("%14.10f\t".format(col)) } println() } } }
object Scala1005 { val N = 3 def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1)) var b:Array[Double] = Array(8,17,20,16) // ピボット選択 pivoting(a,b) println("pivoting") println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("LU") disp_matrix(a) // Ly=b から y を求める (前進代入) var y:Array[Double] = Array(0,0,0,0) forward_substitution(a,b,y) println("Y") disp_vector(y) // Ux=y から x を求める (後退代入) var x:Array[Double] = Array(0,0,0,0) backward_substitution(a,y,x) println("X") disp_vector(x) } // 前進消去 def forward_elimination(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N - 1) { for (row <- pivot + 1 to N) { var s:Double = a(row)(pivot) / a(pivot)(pivot) for (col <- pivot to N) a(row)(col) -= a(pivot)(col) * s // これが 上三角行列 a(row)(pivot) = s // これが 下三角行列 // b(row) -= b(pivot) * s // この値は変更しない } } } // 前進代入 def forward_substitution(a:Array[Array[Double]], b:Array[Double], y:Array[Double]) { for (row <- 0 to N) { for (col <- 0 to row - 1) b(row) -= a(row)(col) * y(col) y(row) = b(row) } } // 後退代入 def backward_substitution(a:Array[Array[Double]], y:Array[Double], x:Array[Double]) { for (row <- (0 to N).reverse) { for (col <- (row + 1 to N).reverse) y(row) -= a(row)(col) * x(col) x(row) = y(row) / a(row)(row) } } // ピボット選択 def pivoting(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) { // 各列で 一番値が大きい行を 探す var max_row:Integer = pivot var max_val:Double = 0 for (row <- pivot to N) { if (Math.abs(a(row)(pivot)) > max_val) { // 一番値が大きい行 max_val = Math.abs(a(row)(pivot)) max_row = row } } // 一番値が大きい行と入れ替え if (max_row != pivot) { var tmp:Double = 0 for (col <- 0 to N) { tmp = a(max_row)(col) a(max_row)(col) = a(pivot)(col) a(pivot)(col) = tmp } tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp } } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 2次元配列を表示 def disp_matrix(matrix:Array[Array[Double]]) = { matrix.foreach { row => row.foreach { col => print("%14.10f\t".format(col)) } println() } } }
object Scala1006 { val N = 3 def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5,2,3,4), Array(2,10,6,7), Array(3,6,15,9), Array(4,7,9,20)) var b:Array[Double] = Array(34,68,96,125) println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("LL^T") disp_matrix(a) // Ly=b から y を求める (前進代入) var y:Array[Double] = Array(0,0,0,0) forward_substitution(a,b,y) println("Y") disp_vector(y) // L^Tx=y から x を求める (後退代入) var x:Array[Double] = Array(0,0,0,0) backward_substitution(a,y,x) println("X") disp_vector(x) } // 前進消去 def forward_elimination(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) { var s:Double = 0; for (col <- 0 to pivot - 1) s += a(pivot)(col) * a(pivot)(col) // ここで根号の中が負の値になると計算できない! a(pivot)(pivot) = Math.sqrt(a(pivot)(pivot) - s) for (row <- pivot + 1 to N) { s = 0 for (col <- 0 to pivot - 1) s += a(row)(col) * a(pivot)(col) a(row)(pivot) = (a(row)(pivot) - s) / a(pivot)(pivot) a(pivot)(row) = a(row)(pivot) } } } // 前進代入 def forward_substitution(a:Array[Array[Double]], b:Array[Double], y:Array[Double]) { for (row <- 0 to N) { for (col <- 0 to row - 1) b(row) -= a(row)(col) * y(col) y(row) = b(row) / a(row)(row) } } // 後退代入 def backward_substitution(a:Array[Array[Double]], y:Array[Double], x:Array[Double]) { for (row <- (0 to N).reverse) { for (col <- (row + 1 to N).reverse) y(row) -= a(row)(col) * x(col) x(row) = y(row) / a(row)(row) } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 2次元配列を表示 def disp_matrix(matrix:Array[Array[Double]]) = { matrix.foreach { row => row.foreach { col => print("%14.10f\t".format(col)) } println() } } }
object Scala1007 { val N = 3 def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5,2,3,4), Array(2,10,6,7), Array(3,6,15,9), Array(4,7,9,20)) var b:Array[Double] = Array(34,68,96,125) println("A") disp_matrix(a) println("B") disp_vector(b) println() // 前進消去 forward_elimination(a,b) println("LDL^T") disp_matrix(a) // Ly=b から y を求める (前進代入) var y:Array[Double] = Array(0,0,0,0) forward_substitution(a,b,y) println("Y") disp_vector(y) // DL^Tx=y から x を求める (後退代入) var x:Array[Double] = Array(0,0,0,0) backward_substitution(a,y,x) println("X") disp_vector(x) } // 前進消去 def forward_elimination(a:Array[Array[Double]], b:Array[Double]) { for (pivot <- 0 to N) { var s:Double = 0 // pivot < k の場合 for (col <- 0 to pivot - 1) { s = a(pivot)(col) for (k <- 0 to col - 1) s -= a(pivot)(k) * a(col)(k) * a(k)(k) a(pivot)(col) = s / a(col)(col) a(col)(pivot) = a(pivot)(col) } // pivot == k の場合 s = a(pivot)(pivot) for (k <- 0 to pivot - 1) s -= a(pivot)(k) * a(pivot)(k) * a(k)(k) a(pivot)(pivot) = s } } // 前進代入 def forward_substitution(a:Array[Array[Double]], b:Array[Double], y:Array[Double]) { for (row <- 0 to N) { for (col <- 0 to row - 1) b(row) -= a(row)(col) * y(col) y(row) = b(row) } } // 後退代入 def backward_substitution(a:Array[Array[Double]], y:Array[Double], x:Array[Double]) { for (row <- (0 to N).reverse) { for (col <- (row + 1 to N).reverse) y(row) -= a(row)(col) * a(row)(row) * x(col) x(row) = y(row) / a(row)(row) } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 2次元配列を表示 def disp_matrix(matrix:Array[Array[Double]]) = { matrix.foreach { row => row.foreach { col => print("%14.10f\t".format(col)) } println() } } }
object Scala1101 { val N = 3 // ベキ乗法で最大固有値を求める def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5.0, 4.0, 1.0, 1.0) ,Array(4.0, 5.0, 1.0, 1.0) ,Array(1.0, 1.0, 4.0, 2.0) ,Array(1.0, 1.0, 2.0, 4.0)) var x:Array[Double] = Array(1.0 ,0.0 ,0.0 ,0.0) // ベキ乗法 val (lambda, x1) = power(a, x) println() println("eigenvalue") println("%14.10f".format(lambda)) println("eigenvector") disp_vector(x1) } // ベキ乗法 def power(a:Array[Array[Double]], x:Array[Double]): (Double, Array[Double]) = { var x0:Array[Double] = normarize(x) var x1:Array[Double] = Array(0.0 ,0.0 ,0.0 ,0.0) var lambda:Double = 0.0 var e:Double = 0.0 var k = 1 do { e = x0.sum // 1次元配列を表示 print("%3d\t".format(k)) disp_vector(x0) k += 1 // 行列の積 x1 = A × x0 x1 = Array(0.0 ,0.0 ,0.0 ,0.0) for (i <- 0 to N) for (j <- 0 to N) x1(i) += a(i)(j) * x0(j) // 内積 var (t1, t2) = (x0 zip x1).map(t => t match {case (f, s) => (s * s, s * f)}).unzip // 固有値 lambda = t1.sum / t2.sum // 正規化 (ベクトル x1 の長さを1にする) x1 = normarize(x1) x0 = x1.clone } while (k <= 100 && Math.abs(e - x1.sum) >= 0.00000000001) return (lambda, x0) } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 正規化 (ベクトルの長さを1にする) def normarize(x:Array[Double]):Array[Double] = { var s:Double = x.map(n => n * n).sum s = Math.sqrt(s) x.map(n => n / s) } }
object Scala1102 { val N = 3 // 逆ベキ乗法で最小固有値を求める def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5.0, 4.0, 1.0, 1.0) ,Array(4.0, 5.0, 1.0, 1.0) ,Array(1.0, 1.0, 4.0, 2.0) ,Array(1.0, 1.0, 2.0, 4.0)) var x:Array[Double] = Array(1.0 ,0.0 ,0.0 ,0.0) // LU分解 forward_elimination(a) // 逆ベキ乗法 val (lambda, x1) = inverse(a, x) println() println("eigenvalue") println("%14.10f".format(lambda)) println("eigenvector") disp_vector(x1) } // 逆ベキ乗法 def inverse(a:Array[Array[Double]], x:Array[Double]): (Double, Array[Double]) = { var x0:Array[Double] = normarize(x) var x1:Array[Double] = Array(0.0 ,0.0 ,0.0 ,0.0) var lambda:Double = 0.0 var e:Double = 0.0 var k = 1 do { e = x0.sum // 1次元配列を表示 print("%3d\t".format(k)) disp_vector(x0) k += 1 // Ly = b から y を求める (前進代入) var b = x0.clone var y:Array[Double] = Array(0.0 ,0.0 ,0.0 ,0.0) forward_substitution(a,y,b) // Ux = y から x を求める (後退代入) x1 = Array(0.0 ,0.0 ,0.0 ,0.0) backward_substitution(a,x1,y) // 内積 var (t1, t2) = (x0 zip x1).map(t => t match {case (f, s) => (s * s, s * f)}).unzip // 固有値 lambda = t2.sum / t1.sum // 正規化 (ベクトル x1 の長さを1にする) x1 = normarize(x1) x0 = x1.clone } while (k <= 100 && Math.abs(e - x1.sum) >= 0.00000000001) return (lambda, x0) } // LU分解 def forward_elimination(a:Array[Array[Double]]) { for (pivot <- 0 to N - 1) { for (row <- pivot + 1 to N) { var s:Double = a(row)(pivot) / a(pivot)(pivot) for (col <- pivot to N) a(row)(col) -= a(pivot)(col) * s // これが 上三角行列 a(row)(pivot) = s // これが 下三角行列 } } } // 前進代入 def forward_substitution(a:Array[Array[Double]], y:Array[Double], b:Array[Double]) { for (row <- 0 to N) { for (col <- 0 to row - 1) b(row) -= a(row)(col) * y(col) y(row) = b(row) } } // 後退代入 def backward_substitution(a:Array[Array[Double]], x:Array[Double], y:Array[Double]) { for (row <- (0 to N).reverse) { for (col <- (row + 1 to N).reverse) y(row) -= a(row)(col) * x(col) x(row) = y(row) / a(row)(row) } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.foreach { col => print("%14.10f\t".format(col)) } println() } // 正規化 (ベクトルの長さを1にする) def normarize(x:Array[Double]):Array[Double] = { var s:Double = x.map(n => n * n).sum s = Math.sqrt(s) x.map(n => n / s) } }
object Scala1103 { val N = 3 // LR分解で固有値を求める def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5.0, 4.0, 1.0, 1.0) ,Array(4.0, 5.0, 1.0, 1.0) ,Array(1.0, 1.0, 4.0, 2.0) ,Array(1.0, 1.0, 2.0, 4.0)) var e:Double = 0.0 var k = 1 do { // LU分解 var (l, u) = decomp(a) // 行列の積 a = multiply(u, l) // 対角要素を表示 print("%3d\t".format(k)) disp_eigenvalue(a) k += 1 // 収束判定 e = (for (i <- 0 to N) yield { for (j <- 0 to (i - 1)) yield { Math.abs(a(i)(j)) } }).flatten.sum } while (k <= 200 && e >= 0.00000000001) println() println("eigenvalue") disp_eigenvalue(a) } // LU分解 def decomp(a:Array[Array[Double]]): (Array[Array[Double]], Array[Array[Double]]) = { var l = Array.ofDim[Double](N + 1, N + 1) var u = Array.ofDim[Double](N + 1, N + 1) l(0)(0) = 1.0 for (j <- 0 to N) u(0)(j) = a(0)(j) for (i <- 1 to N) { u(i)(0) = 0.0 l(0)(i) = 0.0 l(i)(0) = a(i)(0) / u(0)(0) } for (i <- 1 to N) { l(i)(i) = 1.0 u(i)(i) = a(i)(i) - (for (k <- 0 to i) yield (l(i)(k) * u(k)(i))).sum for (j <- i + 1 to N) { l(j)(i) = (a(j)(i) - (for (k <- 0 to i) yield (l(j)(k) * u(k)(i))).sum) / u(i)(i) u(i)(j) = a(i)(j) - (for (k <- 0 to i) yield (l(i)(k) * u(k)(j))).sum } } return (l, u) } // 行列の積 def multiply(a:Array[Array[Double]], b:Array[Array[Double]]): Array[Array[Double]] = { (for (i <- 0 to N) yield { (for (j <- 0 to N) yield { (for (k <- 0 to N) yield a(i)(k) * b(k)(j) ).sum }).toArray }).toArray } // 対角要素を表示 def disp_eigenvalue(matrix:Array[Array[Double]]) { (for (i <- 0 to N) yield (matrix(i)(i))).map(col => print("%14.10f\t".format(col))) println() } }
object Scala1104 { val N = 3 // QR分解で固有値を求める def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5.0, 4.0, 1.0, 1.0) ,Array(4.0, 5.0, 1.0, 1.0) ,Array(1.0, 1.0, 4.0, 2.0) ,Array(1.0, 1.0, 2.0, 4.0)) var e:Double = 0.0 var k = 1 do { // QR分解 var (q, r) = decomp(a) // 行列の積 a = multiply(r, q) // 対角要素を表示 print("%3d\t".format(k)) disp_eigenvalue(a) k += 1 // 収束判定 e = (for (i <- 0 to N) yield { for (j <- 0 to (i - 1)) yield { Math.abs(a(i)(j)) } }).flatten.sum } while (k <= 200 && e >= 0.00000000001) println() println("eigenvalue") disp_eigenvalue(a) } // QR分解 def decomp(a:Array[Array[Double]]): (Array[Array[Double]], Array[Array[Double]]) = { var x = new Array[Double](N + 1) var q = Array.ofDim[Double](N + 1, N + 1) var r = Array.ofDim[Double](N + 1, N + 1) for (k <- 0 to N) { for (i <- 0 to N) x(i) = a(i)(k) for (j <- 0 to k - 1) { r(j)(k) = (for (i <- 0 to N) yield (a(i)(k) * q(i)(j))).sum r(k)(j) = 0.0 for (i <- 0 to N) x(i) -= r(j)(k) * q(i)(j) } r(k)(k) = Math.sqrt(x.map(n => n * n).sum) for (i <- 0 to N) q(i)(k) = x(i) / r(k)(k) } return (q, r) } // 行列の積 def multiply(a:Array[Array[Double]], b:Array[Array[Double]]): Array[Array[Double]] = { (for (i <- 0 to N) yield { (for (j <- 0 to N) yield { (for (k <- 0 to N) yield a(i)(k) * b(k)(j) ).sum }).toArray }).toArray } // 対角要素を表示 def disp_eigenvalue(matrix:Array[Array[Double]]) { (for (i <- 0 to N) yield (matrix(i)(i))).map(col => print("%14.10f\t".format(col))) println() } }
object Scala1105 { val N = 3 // ヤコビ法で固有値を求める def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5.0, 4.0, 1.0, 1.0) ,Array(4.0, 5.0, 1.0, 1.0) ,Array(1.0, 1.0, 4.0, 2.0) ,Array(1.0, 1.0, 2.0, 4.0)) var v:Array[Array[Double]] = Array(Array(1.0, 0.0, 0.0, 0.0) ,Array(0.0, 1.0, 0.0, 0.0) ,Array(0.0, 0.0, 1.0, 0.0) ,Array(0.0, 0.0, 0.0, 1.0)) // ヤコビ法 val (a1, v1) = jacobi(a, v) println() println("eigenvalue") disp_eigenvalue(a1) println() println("eigenvector") disp_eigenvector(v1) } // ヤコビ法 def jacobi(a: => Array[Array[Double]], v: => Array[Array[Double]]): (Array[Array[Double]], Array[Array[Double]]) = { var max_val = 0.0 var k = 1 do { // 最大値を探す var p = 0 var q = 0 max_val = 0.0 for (i <- 0 to N) { for (j <- i + 1 to N) { if (max_val < Math.abs(a(i)(j))) { max_val = Math.abs(a(i)(j)) p = i q = j } } } // θ を求める var t = 0.0 if (Math.abs(a(p)(p) - a(q)(q)) < 0.00000000001) { // a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0 if (a(p)(p) < 0) t = -t } else { // a_{pp} ≠ a_{qq} のとき t = Math.atan(2.0 * a(p)(q) / (a(p)(p) - a(q)(q))) / 2.0 } // θ を使って 行列 U を作成し、A = U^t × A × U val c = Math.cos(t) val s = Math.sin(t) // U^t × A var t1 = 0.0 var t2 = 0.0 for (i <- 0 to N) { t1 = a(p)(i) * c + a(q)(i) * s t2 = -a(p)(i) * s + a(q)(i) * c a(p)(i) = t1 a(q)(i) = t2 // 固有ベクトル t1 = v(p)(i) * c + v(q)(i) * s t2 = -v(p)(i) * s + v(q)(i) * c v(p)(i) = t1 v(q)(i) = t2 } // A × U for (i <- 0 to N) { t1 = a(i)(p) * c + a(i)(q) * s t2 = -a(i)(p) * s + a(i)(q) * c a(i)(p) = t1 a(i)(q) = t2 } // 対角要素を表示 print("%3d\t".format(k)) disp_eigenvalue(a) k += 1 } while (k <= 100 && max_val >= 0.00000000001) return (a, v) } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.map( col => print("%14.10f\t".format(col))) println() } // 正規化 (ベクトルの長さを1にする) def normarize(x:Array[Double]):Array[Double] = { x.map(n => n / Math.sqrt(x.map(n => n * n).sum)) } // 対角要素を表示 def disp_eigenvalue(matrix:Array[Array[Double]]) { (for (i <- 0 to N) yield (matrix(i)(i))).map(col => print("%14.10f\t".format(col))) println() } // 固有ベクトルを表示 def disp_eigenvector(matrix:Array[Array[Double]]) { for (i <- 0 to N) { disp_vector(normarize((for (j <- 0 to N) yield (matrix(i)(j))).toArray)) } } }
object Scala1106 { val N = 3 // ハウスホルダー変換とQR分解で固有値を求める def main(args: Array[String]) { var a:Array[Array[Double]] = Array(Array(5.0, 4.0, 1.0, 1.0) ,Array(4.0, 5.0, 1.0, 1.0) ,Array(1.0, 1.0, 4.0, 2.0) ,Array(1.0, 1.0, 2.0, 4.0)) // ハウスホルダー変換 var (d, e) = tridiagonalize(a) // QR分解 decomp(a, d, e) println() println("eigenvalue") disp_vector(d) println() println("eigenvector") disp_matrix(a) } // ハウスホルダー変換 def tridiagonalize(a: => Array[Array[Double]]): (Array[Double], Array[Double]) = { var v = new Array[Double](N + 1) var d = new Array[Double](N + 1) var e = new Array[Double](N + 1) for (k <- 0 to N - 2) { d(k) = a(k)(k) var w = (for (i <- 0 to N) yield (if (i < k + 1) 0.0 else a(i)(k)) ).toArray var s = Math.sqrt(w.map(n => n * n).sum) if (w(k + 1) < 0) s = -s if (Math.abs(s) < 0.00000000001) e(k + 1) = 0.0 else { e(k + 1) = -s w(k + 1) += s s = 1 / Math.sqrt(w(k + 1) * s) for (i <- k + 1 to N) w(i) *= s for (i <- k + 1 to N) { v(i) = (for (j <- k + 1 to N) yield (if (j <= i) a(i)(j) * w(j) else a(j)(i) * w(j)) ).sum } s = (for (i <- k + 1 to N) yield (w(i) * v(i)) ).sum / 2.0 for (i <- k + 1 to N) v(i) -= s * w(i) for (i <- k + 1 to N) for (j <- k + 1 to i) a(i)(j) -= w(i) * v(j) + w(j) * v(i) // 次の行は固有ベクトルを求めないなら不要 for (i <- k + 1 to N) a(i)(k) = w(i) } } d(N - 1) = a(N - 1)(N - 1) d(N ) = a(N )(N ) e(0) = 0.0 e(N ) = a(N )(N - 1) // 次の行は固有ベクトルを求めないなら不要 for (k <- (N to 0 by -1)) { var w = (for (i <- 0 to N) yield (if (i < k + 1) 0.0 else a(i)(k)) ).toArray for (i <- k + 1 to N) { v(i) = (for (j <- k + 1 to N) yield (a(i)(j) * w(j)) ).sum } for (i <- k + 1 to N) for (j <- k + 1 to N) a(i)(j) -= v(i) * w(j) for (i <- 0 to N) a(i)(k) = 0.0 a(k)(k) = 1.0 } return (d, e) } // QR分解 def decomp(a: => Array[Array[Double]], d: => Array[Double], e: => Array[Double]) { e(0) = 1.0 var h = N while (Math.abs(e(h)) < 0.00000000001) h -= 1 while (h > 0) { e(0) = 0.0 var l = h - 1 while (Math.abs(e(l)) >= 0.00000000001) l -= 1 var j = 1 do { var w = (d(h - 1) - d(h)) / 2.0 var s = Math.sqrt(w * w + e(h) * e(h)) if (w < 0.0) s = -s var x = d(l) - d(h) + e(h) * e(h) / (w + s) var y = e(l + 1) var z = 0.0 var t = 0.0 var u = 0.0 for (k <- l to h - 1) { if (Math.abs(x) >= Math.abs(y)) { t = -y / x u = 1 / Math.sqrt(t * t + 1.0) s = t * u } else { t = -x / y s = 1 / Math.sqrt(t * t + 1.0) * (if (t < 0) -1.0 else 1.0) u = t * s } w = d(k) - d(k + 1) t = (w * s + 2 * u * e(k + 1)) * s d(k ) = d(k ) - t d(k + 1) = d(k + 1) + t e(k ) = u * e(k) - s * z e(k + 1) = e(k + 1) * (u * u - s * s) + w * s * u // 次の行は固有ベクトルを求めないなら不要 for (i <- 0 to N) { x = a(k )(i) y = a(k + 1)(i) a(k )(i) = u * x - s * y a(k + 1)(i) = s * x + u * y } if (k < N - 2) { x = e(k + 1) y = -s * e(k + 2) z = y e(k + 2) = u * e(k + 2) } } print("%3d\t".format(j)) j += 1 disp_vector(d) } while (j <= 100 && Math.abs(e(h)) >= 0.00000000001) e(0) = 1.0 while (Math.abs(e(h)) < 0.00000000001) h -= 1 } // 次の行は固有ベクトルを求めないなら不要 for (k <- 0 to N - 1) { var l = k for (i <- k + 1 to N) if (d(i) > d(l)) l = i var t = d(k) d(k) = d(l) d(l) = t for (i <- 0 to N) { t = a(k)(i) a(k)(i) = a(l)(i) a(l)(i) = t } } } // 1次元配列を表示 def disp_vector(row:Array[Double]) = { row.map( col => print("%14.10f\t".format(col))) println() } // 2次元配列を表示 def disp_matrix(matrix:Array[Array[Double]]) { for (row <- 0 to N) { for (col <- 0 to N) print("%14.10f\t".format(matrix(row)(col))) println() } } }