WScript.Echo 3 + 5 WScript.Echo 3 - 5 WScript.Echo 3 * 5 WScript.Echo 3 ^ 5 WScript.Echo 5 / 3 WScript.Echo 5 \ 3 WScript.Echo 5 Mod 3 WScript.StdOut.Write 3 * 5 & vbNewLine WScript.StdOut.WriteLine 3 * 5
Option Explicit Dim i: i = 3 * 5 WScript.Echo "3 * 5 = " & i
Option Explicit Dim i For i = 1 To 9 WScript.StdOut.Write i & ", " Next WScript.StdOut.WriteLine
Option Explicit Dim i For i = 1 To 9 If i Mod 3 = 0 Then WScript.StdOut.Write i & ", " End If Next WScript.StdOut.WriteLine
Option Explicit Dim sum: sum = 0 Dim i For i = 1 To 99 If i Mod 3 = 0 Then sum = sum + i End If Next WScript.Echo sum
Option Explicit ' 初項:a, 公差:a で 上限:lim の数列の総和を返す関数 Private Function sn(a, lim) Dim n: n = lim \ a ' 項数:n = 上限:lim / 公差:a Dim l: l = n * a ' 末項:l = 項数:n * 公差:a sn = (a + l) * n \ 2 ' 総和:sn = (初項:a + 末項:l) * 項数:n / 2 End Function ' 3 の倍数の合計 WScript.Echo sn(3, 999)
Option Explicit ' 10000 までの 自然数の和 ' 項目数 n = 10000 Dim n: n = 10000 WScript.Echo n * (n + 1) \ 2
Option Explicit ' 10000 までの 偶数の和 ' 項目数 n = 5000 Dim n: n = 10000 \ 2 WScript.Echo n * (n + 1)
Option Explicit ' 10000 までの 奇数の和 ' 項目数 n = 5000 Dim n: n = 10000 \ 2 WScript.Echo n ^ 2
Option Explicit ' 1000 までの 自然数の2乗の和 Dim n: n = 1000 WScript.Echo n * (n + 1) * (2 * n + 1) \ 6
Option Explicit ' 100 までの 自然数の3乗の和 Dim n: n = 100 WScript.Echo (n ^ 2) * ((n + 1) ^ 2) \ 4
Option Explicit ' 初項 2, 公比 3, 項数 10 の等比数列の和 Dim n: n = 10 Dim a: a = 2 Dim r: r = 3 WScript.Echo (a * ((r ^ n) - 1)) \ (r - 1)
Option Explicit Dim a: a = 5 '初項 5 Dim d: d = 3 '公差 3 Dim n: n = 10 '項数 10 Dim p: p = 1 '積 Dim m Dim i For i = 1 To n m = a + (d * (i - 1)) WScript.StdOut.Write m & " * " & p & " = " p = p * m WScript.Echo p Next
Option Explicit '初項 5, 公差 3, 項数 10 の数列の積を表示する WScript.Echo product(5, 3, 10) Private Function product(m, d, n) If n = 0 Then product = 1 Else product = m * product(m + d, d, n - 1) End If End Function
'階乗を求める関数 Private Function Fact(n) If n <= 1 Then Fact = 1 Else Fact = n * Fact(n - 1) End If End Function '10の階乗 WScript.Echo(Fact(10)) WScript.Echo(10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1)
'下降階乗冪 Private Function FallingFact(x, n) If n <= 1 Then FallingFact = x Else FallingFact = x * FallingFact(x - 1, n - 1) End If End Function '10 から 6 までの 総乗 WScript.Echo(FallingFact(10, 5)) WScript.Echo(10 * 9 * 8 * 7 * 6)
'上昇階乗冪 Private Function RisingFact(x, n) If n <= 1 Then RisingFact = x Else RisingFact = x * RisingFact(x + 1, n - 1) End If End Function '10 から 14 までの 総乗 WScript.Echo(RisingFact(10, 5)) WScript.Echo(10 * 11 * 12 * 13 * 14)
'階乗 Private Function Fact(n) If n <= 1 Then Fact = 1 Else Fact = n * Fact(n - 1) End If End Function '下降階乗冪 Private Function FallingFact(x, n) If n <= 1 Then FallingFact = x Else FallingFact = x * FallingFact(x - 1, n - 1) End If End Function '順列 (異なる 10 個のものから 5 個取ってできる順列の総数) Dim n: n = 10 Dim r: r = 5 WScript.Echo(Fact(n) / Fact(n - r)) WScript.Echo(FallingFact(n, r))
'重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) Dim n: n = 10 Dim r: r = 5 WScript.Echo(n ^ r)
'組合せ Private Function Comb(n, r) If (r = 0) Or (r = n) Then Comb = 1 ElseIf r = 1 Then Comb = n Else Comb = Comb(n - 1, r - 1) + Comb(n - 1, r) End If End Function '組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数) Dim n: n = 10 Dim r: r = 5 WScript.Echo(Comb(n, r))
'組合せ Private Function Comb(n, r) If (r = 0) Or (r = n) Then Comb = 1 ElseIf r = 1 Then Comb = n Else Comb = Comb(n - 1, r - 1) + Comb(n - 1, r) End If End Function '重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数) Dim n: n = 10 Dim r: r = 5 WScript.Echo(Comb(n + r - 1, r))
Option Explicit Const PI = 3.14159265359 Dim degree For degree = 0 To 360 Step 15 If (degree Mod 30 = 0 Or degree Mod 45 = 0) Then Dim radian: radian = degree * PI / 180.0 '自作の正弦関数 Dim d1: d1 = mySin(radian, 1, False, radian, 1.0, radian) '標準の正弦関数 Dim d2: d2 = Sin(radian) '標準関数との差異 WScript.StdOut.Write Right(Space(3) & degree, 3) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1, 10, -1, 0, 0), 13) & " - " WScript.StdOut.Write Right(Space(13) & FormatNumber(d2, 10, -1, 0, 0), 13) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine End If Next '自作の正弦関数 Private Function mySin(ByVal x, ByVal n, ByVal nega, ByVal numerator, ByVal denominator, ByVal y) Dim m: m = 2 * n denominator = denominator * (m + 1) * m numerator = numerator * x * x Dim a: a = numerator / denominator '十分な精度になったら処理を抜ける If (a <= 0.00000000001) Then mySin = y Else If Not nega Then a = -a mySin = y + mySin(x, n + 1, Not nega, numerator, denominator, a) End If End Function
Option Explicit Const PI = 3.14159265359 Dim degree For degree = 0 To 360 Step 15 If (degree Mod 30 = 0 Or degree Mod 45 = 0) Then Dim radian: radian = degree * PI / 180.0 '自作の余弦関数 Dim d1: d1 = myCos(radian, 1, False, 1.0, 1.0, 1.0) '標準の余弦関数 Dim d2: d2 = Cos(radian) '標準関数との差異 WScript.StdOut.Write Right(Space(3) & degree, 3) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1, 10, -1, 0, 0), 13) & " - " WScript.StdOut.Write Right(Space(13) & FormatNumber(d2, 10, -1, 0, 0), 13) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine End If Next '自作の余弦関数 Private Function myCos(ByVal x, ByVal n, ByVal nega, ByVal numerator, ByVal denominator, ByVal y) Dim m: m = 2 * n denominator = denominator * m * (m - 1) numerator = numerator * x * x Dim a: a = numerator / denominator '十分な精度になったら処理を抜ける If (a <= 0.00000000001) Then myCos = y Else If Not nega Then a = -a myCos = y + myCos(x, n + 1, Not nega, numerator, denominator, a) End If End Function
Option Explicit Const PI = 3.14159265359 Dim i For i = 0 To 180 Step 15 If (i Mod 180 <> 0) Then Dim degree: degree = i - 90 Dim radian: radian = degree * PI / 180.0 Dim x2: x2 = radian * radian '自作の正接関数 Dim d1: d1 = myTan(radian, x2, 15, 0.0) '15:必要な精度が得られる十分大きな奇数 '標準の正接関数 Dim d2: d2 = Tan(radian) '標準関数との差異 WScript.StdOut.Write Right(Space(3) & degree, 3) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1, 10, -1, 0, 0), 13) & " - " WScript.StdOut.Write Right(Space(13) & FormatNumber(d2, 10, -1, 0, 0), 13) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine End If Next '自作の正接関数 Private Function myTan(ByVal x, ByVal x2, ByVal n, ByVal t) t = x2 / (n - t) n = n - 2 If (n <= 1) Then myTan = x / (1 - t) Else myTan = myTan(x, x2, n, t) End If End Function
Option Explicit Dim i For i = 0 To 20 Dim x: x = (i - 10) / 4.0 '標準の指数関数 Dim d1: d1 = Exp(x) '自作の指数関数 Dim d2: d2 = myExp(x, 1, 1.0, 1.0, 1.0) '標準関数との差異 WScript.StdOut.Write Right(Space(5) & FormatNumber(x, 2, -1, 0, 0), 5) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1, 10, -1, 0, 0), 13) & " - " WScript.StdOut.Write Right(Space(13) & FormatNumber(d2, 10, -1, 0, 0), 13) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine Next '自作の指数関数 Private Function myExp(ByVal x, ByVal n, ByVal numerator, ByVal denominator, ByVal y) denominator = denominator * n numerator = numerator * x Dim a: a = numerator / denominator '十分な精度になったら処理を抜ける If (Abs(a) <= 0.00000000001) Then myExp = y Else myExp = y + myExp(x, n + 1, numerator, denominator, a) End If End Function
Option Explicit Dim i For i = 0 To 20 Dim x: x = (i - 10) / 4.0 '標準の指数関数 Dim d1: d1 = Exp(x) '自作の指数関数 Dim x2: x2 = x * x Dim d2: d2 = myExp(x, x2, 30, 0.0) '30:必要な精度が得られるよう, 6から始めて4ずつ増加させる '標準関数との差異 WScript.StdOut.Write Right(Space(5) & FormatNumber(x, 2, -1, 0, 0), 5) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1, 10, -1, 0, 0), 13) & " - " WScript.StdOut.Write Right(Space(13) & FormatNumber(d2, 10, -1, 0, 0), 13) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine Next '自作の指数関数 Private Function myExp(ByVal x, ByVal x2, ByVal n, ByVal t) t = x2 / (n + t) n = n - 4 If (n < 6) Then myExp = 1 + ((2 * x) / (2 - x + t)) Else myExp = myExp(x, x2, n, t) End If End Function
Option Explicit Dim i For i = 1 To 20 Dim x: x = i / 5.0 '標準の対数関数 Dim d1: d1 = Log(x) '自作の対数関数 Dim x2: x2 = (x - 1) / (x + 1) Dim d2: d2 = 2 * myLog(x2, x2, 1.0, x2) '標準関数との差異 WScript.StdOut.Write Right(Space(5) & FormatNumber(x, 2, -1, 0, 0), 5) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1, 10, -1, 0, 0), 13) & " - " WScript.StdOut.Write Right(Space(13) & FormatNumber(d2, 10, -1, 0, 0), 13) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine Next '自作の対数関数 Private Function myLog(ByVal x2, ByVal numerator, ByVal denominator, ByVal y) denominator = denominator + 2 numerator = numerator * x2 * x2 Dim a: a = numerator / denominator '十分な精度になったら処理を抜ける If (Abs(a) <= 0.00000000001) Then myLog = y Else myLog = y + myLog(x2, numerator, denominator, a) End If End Function
Option Explicit Dim i For i = 1 To 20 Dim x: x = i / 5.0 '標準の対数関数 Dim d1: d1 = Log(x) '自作の対数関数 Dim d2: d2 = myLog(x - 1, 27, 0.0) '27:必要な精度が得られる十分大きな奇数 '標準関数との差異 WScript.StdOut.Write Right(Space(5) & FormatNumber(x, 2, -1, 0, 0), 5) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1, 10, -1, 0, 0), 13) & " - " WScript.StdOut.Write Right(Space(13) & FormatNumber(d2, 10, -1, 0, 0), 13) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine Next '自作の対数関数 Private Function myLog(ByVal x, ByVal n, ByVal t) Dim n2: n2 = n Dim x2: x2 = x If (n > 3) Then If (n Mod 2 = 0) Then n2 = 2 End If x2 = x * (n \ 2) End If t = x2 / (n2 + t) If (n <= 2) Then myLog = x / (1 + t) Else myLog = myLog(x, n - 1, t) End If End Function
Option Explicit Dim i For i = 0 To 20 Dim x: x = i - 10 '自作の双曲線正弦関数 Dim d1: d1 = mySinh(x, 1, x, 1.0, x) '標準の双曲線正弦関数はない Dim d2: d2 = (Exp(x) - Exp(-x)) / 2 '標準関数との差異 WScript.StdOut.Write Right(Space(3) & x, 3) & " : " WScript.StdOut.Write Right(Space(17) & FormatNumber(d1, 10, -1, 0, 0), 17) & " - " WScript.StdOut.Write Right(Space(17) & FormatNumber(d2, 10, -1, 0, 0), 17) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine Next '自作の双曲線正弦関数 Private Function mySinh(ByVal x, ByVal n, ByVal numerator, ByVal denominator, ByVal y) Dim m: m = 2 * n denominator = denominator * (m + 1) * m numerator = numerator * x * x Dim a: a = numerator / denominator '十分な精度になったら処理を抜ける If (Abs(a) <= 0.00000000001) Then mySinh = y Else mySinh = y + mySinh(x, n + 1, numerator, denominator, a) End If End Function
Option Explicit Dim i For i = 0 To 20 Dim x: x = i - 10 '自作の双曲線余弦関数 Dim d1: d1 = myCosh(x, 1, 1.0, 1.0, 1.0) '標準の双曲線余弦関数はない Dim d2: d2 = (Exp(x) + Exp(-x)) / 2 '標準関数との差異 WScript.StdOut.Write Right(Space(3) & x, 3) & " : " WScript.StdOut.Write Right(Space(17) & FormatNumber(d1, 10, -1, 0, 0), 17) & " - " WScript.StdOut.Write Right(Space(17) & FormatNumber(d2, 10, -1, 0, 0), 17) & " = " WScript.StdOut.Write Right(Space(13) & FormatNumber(d1 - d2, 10, -1, 0, 0), 13) & vbNewLine Next '自作の双曲線余弦関数 Private Function myCosh(ByVal x, ByVal n, ByVal numerator, ByVal denominator, ByVal y) Dim m: m = 2 * n denominator = denominator * m * (m - 1) numerator = numerator * x * x Dim a: a = numerator / denominator '十分な精度になったら処理を抜ける If (Abs(a) <= 0.00000000001) Then myCosh = y Else myCosh = y + myCosh(x, n + 1, numerator, denominator, a) End If End Function
Option Explicit Const PI = 3.14159265359 Const a = 0 Const b = 1 '台形則で積分 Dim n: n = 2 Dim i, j For j = 1 To 10 Dim h: h = (b - a) / n Dim s: s = 0 Dim x: x = a For i = 1 To n - 1 x = x + h s = s + f(x) Next s = h * ((f(a) + f(b)) / 2 + s) n = n * 2 '結果を π と比較 WScript.StdOut.Write Right(Space(2) & j, 2) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(s, 10, -1, 0, 0), 13) & ", " WScript.StdOut.Write Right(Space(13) & FormatNumber(s - PI, 10, -1, 0, 0), 13) & vbNewLine Next Private Function f(x) f = 4 / (1 + x * x) End Function
Option Explicit Const PI = 3.14159265359 Const a = 0 Const b = 1 '中点則で積分 Dim n: n = 2 Dim i, j For j = 1 To 10 Dim h: h = (b - a) / n Dim s: s = 0 Dim x: x = a + (h / 2) For i = 1 To n s = s + f(x) x = x + h Next s = h * s n = n * 2 '結果を π と比較 WScript.StdOut.Write Right(Space(2) & j, 2) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(s, 10, -1, 0, 0), 13) & ", " WScript.StdOut.Write Right(Space(13) & FormatNumber(s - PI, 10, -1, 0, 0), 13) & vbNewLine Next Private Function f(x) f = 4 / (1 + x * x) End Function
Option Explicit Const PI = 3.14159265359 Const a = 0 Const b = 1 'Simpson則で積分 Dim n: n = 2 Dim i, j For j = 1 To 5 Dim h: h = (b - a) / n Dim s2: s2 = 0 Dim s4: s4 = 0 Dim x: x = a + h For i = 1 To n \ 2 s4 = s4 + f(x) x = x + h s2 = s2 + f(x) x = x + h Next s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 = s4 * 4 Dim s: s = (s2 + s4) * h / 3 n = n * 2 '結果を π と比較 WScript.StdOut.Write Right(Space(2) & j, 2) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(s, 10, -1, 0, 0), 13) & ", " WScript.StdOut.Write Right(Space(13) & FormatNumber(s - PI, 10, -1, 0, 0), 13) & vbNewLine Next Private Function f(x) f = 4 / (1 + x * x) End Function
Option Explicit Const PI = 3.14159265359 Const a = 0 Const b = 1 Dim t(6, 6) '台形則で積分 Dim n: n = 2 Dim i, j For i = 1 To 6 Dim h: h = (b - a) / n Dim s: s = 0 Dim x: x = a For j = 1 To n - 1 x = x + h s = s + f(x) Next '結果を保存 t(i,1) = h * ((f(a) + f(b)) / 2 + s) n = n * 2 Next '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 Then '結果を π と比較 WScript.StdOut.Write Right(Space(2) & j, 2) & " : " WScript.StdOut.Write Right(Space(13) & FormatNumber(t(i,j), 10, -1, 0, 0), 13) & ", " WScript.StdOut.Write Right(Space(13) & FormatNumber(t(i,j) - PI, 10, -1, 0, 0), 13) & vbNewLine End If Next n = n * 4 Next Private Function f(x) f = 4 / (1 + x * x) End Function
Option Explicit 'データ点の数 - 1 Private Const N = 6 Dim x(): ReDim x(N) Dim y(): ReDim y(N) '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット Dim i For i = 0 To N Dim d: d = i * 1.5 - 4.5 x(i) = d y(i) = f(d) Next '0.5刻みで 与えられていない値を補間 For i = 0 To 18 d = i * 0.5 - 4.5 Dim d1: d1 = f(d) Dim d2: d2 = lagrange(d, x, y) '元の関数と比較 WScript.StdOut.Write Right(Space(5) & FormatNumber(d, 2, -1, 0, 0), 5) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d1, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d2, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d1 - d2, 5, -1, 0, 0), 8) Next '元の関数 Private Function f(ByVal x) f = x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function 'Lagrange (ラグランジュ) 補間 Private Function lagrange(ByVal d, ByVal x(), ByVal y()) Dim sum: sum = 0 Dim i, j For i = 0 To N Dim prod: prod = y(i) For j = 0 To N If j <> i Then prod = prod * (d - x(j)) / (x(i) - x(j)) End If Next sum = sum + prod Next lagrange = sum End Function
Option Explicit 'データ点の数 - 1 Private Const N = 6 Dim x(): ReDim x(N) Dim y(): ReDim y(N) '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット Dim i For i = 0 To N Dim d: d = i * 1.5 - 4.5 x(i) = d y(i) = f(d) Next '0.5刻みで 与えられていない値を補間 For i = 0 To 18 d = i * 0.5 - 4.5 Dim d1: d1 = f(d) Dim d2: d2 = neville(d, x, y) '元の関数と比較 WScript.StdOut.Write Right(Space(5) & FormatNumber(d, 2, -1, 0, 0), 5) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d1, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d2, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d1 - d2, 5, -1, 0, 0), 8) Next '元の関数 Private Function f(ByVal x) f = x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function 'Neville (ネヴィル) 補間 Private Function neville(ByVal d, ByVal x(), ByVal y()) Dim w(): ReDim w(N, N) Dim i For i = 0 To N w(0,i) = y(i) Next Dim j 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)) Next Next neville = w(N,0) End Function
Option Explicit 'データ点の数 - 1 Private Const N = 6 Dim x(): ReDim x(N) Dim y(): ReDim y(N) '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット Dim i For i = 0 To N Dim d1: d1 = i * 1.5 - 4.5 x(i) = d1 y(i) = f(d1) Next '差分商の表を作る Dim d(): ReDim d(N, N) Dim j For j = 0 To N d(0,j) = y(j) Next 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)) Next Next 'n階差分商 Dim a(): ReDim a(N) For j = 0 To N a(j) = d(j,0) Next '0.5刻みで 与えられていない値を補間 For i = 0 To 18 d1 = i * 0.5 - 4.5 Dim d2: d2 = f(d1) Dim d3: d3 = newton(d1, x, a) '元の関数と比較 WScript.StdOut.Write Right(Space(5) & FormatNumber(d1, 2, -1, 0, 0), 5) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d2, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d3, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d2 - d3, 5, -1, 0, 0), 8) Next '元の関数 Private Function f(ByVal x) f = x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function 'Newton (ニュートン) 補間 Private Function newton(ByVal d, ByVal x(), ByVal a()) Dim sum: sum = a(0) Dim i, j For i = 1 To N Dim prod: prod = a(i) For j = 0 To (i - 1) If j <> i Then prod = prod * (d - x(j)) End If Next sum = sum + prod Next newton = sum End Function
Option Explicit 'データ点の数 - 1 Private Const N = 6 Private Const Nx2 = 13 Dim x(): ReDim x(N) Dim y(): ReDim y(N) Dim yd(): ReDim yd(N) '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット Dim i For i = 0 To N Dim d1: d1 = i * 1.5 - 4.5 x(i) = d1 y(i) = f(d1) yd(i) = fd(d1) Next '差分商の表を作る Dim z(): ReDim z(Nx2) Dim d(): ReDim d(Nx2, Nx2) For i = 0 To Nx2 Dim j: j = i \ 2 z(i) = x(j) d(0,i) = y(j) Next For i = 1 To Nx2 For j = 0 To (Nx2 - i) If i = 1 And j Mod 2 = 0 Then d(i,j) = yd(j \ 2) Else d(i,j) = (d(i-1,j+1) - d(i-1,j)) / (z(j+i) - z(j)) End If Next Next 'n階差分商 Dim a(): ReDim a(Nx2) For j = 0 To Nx2 a(j) = d(j,0) Next '0.5刻みで 与えられていない値を補間 For i = 0 To 18 d1 = i * 0.5 - 4.5 Dim d2: d2 = f(d1) Dim d3: d3 = hermite(d1, z, a) '元の関数と比較 WScript.StdOut.Write Right(Space(5) & FormatNumber(d1, 2, -1, 0, 0), 5) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d2, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d3, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d2 - d3, 5, -1, 0, 0), 8) Next '元の関数 Private Function f(ByVal x) f = x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function '導関数 Private Function fd(ByVal x) fd = 1 - (x ^ 2) / 2 + (x ^ 4) / (4 * 3 * 2) End Function 'Hermite (エルミート) 補間 Private Function hermite(ByVal d, ByVal z(), ByVal a()) Dim sum: sum = a(0) Dim i, j For i = 1 To Nx2 Dim prod: prod = a(i) For j = 0 To (i - 1) If j <> i Then prod = prod * (d - z(j)) End If Next sum = sum + prod Next hermite = sum End Function
Option Explicit 'データ点の数 - 1 Private Const N = 6 Dim x(): ReDim x(N) Dim y(): ReDim y(N) '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット Dim i For i = 0 To N Dim d1: d1 = i * 1.5 - 4.5 x(i) = d1 y(i) = f(d1) Next '3項方程式の係数の表を作る Dim a(): ReDim a(N) Dim b(): ReDim b(N) Dim c(): ReDim c(N) Dim d(): ReDim d(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))) Next '3項方程式を解く (ト−マス法) Dim g(): ReDim g(N) Dim s(): ReDim s(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) Next Dim z(): ReDim z(N) z(0) = 0 z(N) = 0 z(N-1) = s(N-1) / g(N-1) For i = N - 2 To 1 Step -1 z(i) = (s(i) - c(i) * z(i+1)) / g(i) Next '0.5刻みで 与えられていない値を補間 For i = 0 To 18 d1 = i * 0.5 - 4.5 Dim d2: d2 = f(d1) Dim d3: d3 = spline(d1, x, y, z) '元の関数と比較 WScript.StdOut.Write Right(Space(5) & FormatNumber(d1, 2, -1, 0, 0), 5) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d2, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(d3, 5, -1, 0, 0), 8) & vbTab WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d2 - d3, 5, -1, 0, 0), 8) Next '元の関数 Private Function f(ByVal x) f = x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function 'Spline (スプライン) 補間 Private Function spline(ByVal d, ByVal x(), ByVal y(), ByVal z()) '補間関数値がどの区間にあるか Dim k: k = -1 Dim i For i = 1 To N If d <= x(i) Then k = i - 1 Exit For End If Next If k < 0 Then k = N Dim d1: d1 = x(k+1) - d Dim d2: d2 = d - x(k) Dim d3: d3 = x(k+1) - x(k) spline = (z(k) * (d1 ^ 3) + z(k+1) * (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 End Function
Option Explicit Private Const PI = 3.14159265359 '重力加速度 Private Const g = -9.8 '空気抵抗係数 Private Const k = -0.01 '時間間隔(秒) Private Const h = 0.01 '角度 Private Const degree = 45 Private radian: radian = degree * PI / 180.0 '初速 250 km/h -> 秒速に変換 Private v: v = 250 * 1000 \ 3600 '水平方向の速度 Private vx(): ReDim vx(1) vx(0) = v * Cos(radian) '鉛直方向の速度 Private vy(): ReDim vy(1) vy(0) = v * Sin(radian) '経過秒数 Private t: t = 0.0 '位置 Private x: x = 0.0 Private y: y = 0.0 '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx, ByVal vy) fx = k * Sqr(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx, ByVal vy) fy = g + (k * Sqr(vx * vx + vy * vy) * vy) End Function 'Euler法 Dim i: i = 1 Do While (y >= 0.0) '経過秒数 t = i * h '位置 x = x + h * vx(0) y = y + h * vy(0) WScript.StdOut.Write Right(Space(4) & FormatNumber(t, 2, -1, 0, 0), 4) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(vx(0), 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(vy(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(x, 5, -1, 0, 0), 9) & vbTab WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(y, 5, -1, 0, 0), 8) '速度 vx(1) = vx(0) + h * fx(vx(0), vy(0)) vy(1) = vy(0) + h * fy(vx(0), vy(0)) vx(0) = vx(1) vy(0) = vy(1) i = i + 1 Loop
Option Explicit Private Const PI = 3.14159265359 '重力加速度 Private Const g = -9.8 '空気抵抗係数 Private Const k = -0.01 '時間間隔(秒) Private Const h = 0.01 '角度 Private Const degree = 45 Private radian: radian = degree * PI / 180.0 '初速 250 km/h -> 秒速に変換 Private v: v = 250 * 1000 \ 3600 '水平方向の速度 Private vx(): ReDim vx(2) vx(0) = v * Cos(radian) '鉛直方向の速度 Private vy(): ReDim vy(2) vy(0) = v * Sin(radian) '経過秒数 Private t: t = 0.0 '位置 Private x(): ReDim x(2) x(0) = 0.0 Private y(): ReDim y(2) y(0) = 0.0 '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx, ByVal vy) fx = k * Sqr(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx, ByVal vy) fy = g + (k * Sqr(vx * vx + vy * vy) * vy) End Function 'Heun法 Dim i: i = 1 Do While (y(0) >= 0.0) '経過秒数 t = i * h '位置・速度 x(1) = x(0) + h * vx(0) y(1) = y(0) + h * vy(0) vx(1) = vx(0) + h * fx(vx(0), vy(0)) vy(1) = vy(0) + h * fy(vx(0), vy(0)) x(2) = x(0) + h * ( vx(0) + vx(1) ) / 2 y(2) = y(0) + h * ( vy(0) + vy(1) ) / 2 vx(2) = vx(0) + h * (fx(vx(0), vy(0)) + fx(vx(1), vy(1))) / 2 vy(2) = vy(0) + h * (fy(vx(0), vy(0)) + fy(vx(1), vy(1))) / 2 x(0) = x(2) y(0) = y(2) vx(0) = vx(2) vy(0) = vy(2) WScript.StdOut.Write Right(Space(4) & FormatNumber(t, 2, -1, 0, 0), 4) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(vx(0), 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(vy(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(x(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(y(0), 5, -1, 0, 0), 8) i = i + 1 Loop
Option Explicit Private Const PI = 3.14159265359 '重力加速度 Private Const g = -9.8 '空気抵抗係数 Private Const k = -0.01 '時間間隔(秒) Private Const h = 0.01 '角度 Private Const degree = 45 Private radian: radian = degree * PI / 180.0 '初速 250 km/h -> 秒速に変換 Private v: v = 250 * 1000 \ 3600 '水平方向の速度 Private vx(): ReDim vx(1) vx(0) = v * Cos(radian) '鉛直方向の速度 Private vy(): ReDim vy(1) vy(0) = v * Sin(radian) '経過秒数 Private t: t = 0.0 '位置 Private x(): ReDim x(1) x(0) = 0.0 Private y(): ReDim y(1) y(0) = 0.0 '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx, ByVal vy) fx = k * Sqr(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx, ByVal vy) fy = g + (k * Sqr(vx * vx + vy * vy) * vy) End Function '中点法 Dim i: i = 1 Do While (y(0) >= 0.0) '経過秒数 t = i * h '位置・速度 vx(1) = h * fx(vx(0), vy(0)) vy(1) = h * fy(vx(0), vy(0)) Dim wx: wx = vx(0) + vx(1) / 2.0 Dim wy: wy = vy(0) + vy(1) / 2.0 vx(0) = vx(0) + h * fx(wx, wy) vy(0) = vy(0) + h * fy(wx, wy) x(0) = x(0) + h * wx y(0) = y(0) + h * wy WScript.StdOut.Write Right(Space(4) & FormatNumber(t, 2, -1, 0, 0), 4) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(vx(0), 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(vy(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(x(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(y(0), 5, -1, 0, 0), 8) i = i + 1 Loop
Option Explicit Private Const PI = 3.14159265359 '重力加速度 Private Const g = -9.8 '空気抵抗係数 Private Const k = -0.01 '時間間隔(秒) Private Const h = 0.01 '角度 Private Const degree = 45 Private radian: radian = degree * PI / 180.0 '初速 250 km/h -> 秒速に変換 Private v: v = 250 * 1000 \ 3600 '水平方向の速度 Private vx(): ReDim vx(5) vx(0) = v * Cos(radian) '鉛直方向の速度 Private vy(): ReDim vy(5) vy(0) = v * Sin(radian) '経過秒数 Private t: t = 0.0 '位置 Private x(): ReDim x(5) x(0) = 0.0 Private y(): ReDim y(5) y(0) = 0.0 '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx, ByVal vy) fx = k * Sqr(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx, ByVal vy) fy = g + (k * Sqr(vx * vx + vy * vy) * vy) End Function 'Runge-Kutta法 Dim i: i = 1 Do While (y(0) >= 0.0) '経過秒数 t = i * h '位置・速度 x(1) = h * vx(0) y(1) = h * vy(0) vx(1) = h * fx(vx(0), vy(0)) vy(1) = h * fy(vx(0), vy(0)) Dim wx: wx = vx(0) + vx(1) / 2.0 Dim wy: wy = vy(0) + vy(1) / 2.0 x(2) = h * wx y(2) = h * wy vx(2) = h * fx(wx, wy) vy(2) = h * fy(wx, wy) wx = vx(0) + vx(2) / 2 wy = vy(0) + vy(2) / 2 x(3) = h * wx y(3) = h * wy vx(3) = h * fx(wx, wy) vy(3) = h * fy(wx, wy) wx = vx(0) + vx(3) wy = vy(0) + vy(3) x(4) = h * wx y(4) = h * wy vx(4) = h * fx(wx, wy) vy(4) = h * fy(wx, wy) x(0) = x(0) + ( x(1) + x(2) * 2 + x(3) * 2 + x(4)) / 6 y(0) = y(0) + ( y(1) + y(2) * 2 + y(3) * 2 + y(4)) / 6 vx(0) = vx(0) + (vx(1) + vx(2) * 2 + vx(3) * 2 + vx(4)) / 6 vy(0) = vy(0) + (vy(1) + vy(2) * 2 + vy(3) * 2 + vy(4)) / 6 WScript.StdOut.Write Right(Space(4) & FormatNumber(t, 2, -1, 0, 0), 4) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(vx(0), 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(vy(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(x(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.WriteLine Right(Space(9) & FormatNumber(y(0), 5, -1, 0, 0), 9) i = i + 1 Loop
Option Explicit Private Const PI = 3.14159265359 '重力加速度 Private Const g = -9.8 '空気抵抗係数 Private Const k = -0.01 '時間間隔(秒) Private Const h = 0.01 '角度 Private Const degree = 45 Private radian: radian = degree * PI / 180.0 '初速 250 km/h -> 秒速に変換 Private v: v = 250 * 1000 \ 3600 '水平方向の速度 Private vx(): ReDim vx(5) vx(0) = v * Cos(radian) '鉛直方向の速度 Private vy(): ReDim vy(5) vy(0) = v * Sin(radian) '経過秒数 Private t: t = 0.0 '位置 Private x(): ReDim x(5) x(0) = 0.0 Private y(): ReDim y(5) y(0) = 0.0 '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx, ByVal vy) fx = k * Sqr(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx, ByVal vy) fy = g + (k * Sqr(vx * vx + vy * vy) * vy) End Function 'Runge-Kutta-Gill法 Dim i: i = 1 Do While (y(0) >= 0.0) '経過秒数 t = i * h '位置・速度 x(1) = h * vx(0) y(1) = h * vy(0) vx(1) = h * fx(vx(0), vy(0)) vy(1) = h * fy(vx(0), vy(0)) Dim wx: wx = vx(0) + vx(1) / 2.0 Dim wy: wy = vy(0) + vy(1) / 2.0 x(2) = h * wx y(2) = h * wy vx(2) = h * fx(wx, wy) vy(2) = h * fy(wx, wy) wx = vx(0) + vx(1) * ((Sqr(2.0) - 1) / 2) + vx(2) * (1 - 1 / Sqr(2.0)) wy = vy(0) + vy(1) * ((Sqr(2.0) - 1) / 2) + vy(2) * (1 - 1 / Sqr(2.0)) x(3) = h * wx y(3) = h * wy vx(3) = h * fx(wx, wy) vy(3) = h * fy(wx, wy) wx = vx(0) - vx(2) / Sqr(2.0) + vx(3) * (1 + 1 / Sqr(2.0)) wy = vy(0) - vy(2) / Sqr(2.0) + vy(3) * (1 + 1 / Sqr(2.0)) x(4) = h * wx y(4) = h * wy vx(4) = h * fx(wx, wy) vy(4) = h * fy(wx, wy) x(0) = x(0) + ( x(1) + x(2) * (2 - Sqr(2.0)) + x(3) * (2 + Sqr(2.0)) + x(4)) / 6 y(0) = y(0) + ( y(1) + y(2) * (2 - Sqr(2.0)) + y(3) * (2 + Sqr(2.0)) + y(4)) / 6 vx(0) = vx(0) + (vx(1) + vx(2) * (2 - Sqr(2.0)) + vx(3) * (2 + Sqr(2.0)) + vx(4)) / 6 vy(0) = vy(0) + (vy(1) + vy(2) * (2 - Sqr(2.0)) + vy(3) * (2 + Sqr(2.0)) + vy(4)) / 6 WScript.StdOut.Write Right(Space(4) & FormatNumber(t, 2, -1, 0, 0), 4) & vbTab WScript.StdOut.Write Right(Space(8) & FormatNumber(vx(0), 5, -1, 0, 0), 8) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(vy(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.Write Right(Space(9) & FormatNumber(x(0), 5, -1, 0, 0), 9) & vbTab WScript.StdOut.WriteLine Right(Space(9) & FormatNumber(y(0), 5, -1, 0, 0), 9) i = i + 1 Loop
Option Explicit Dim a: a = 1.0 Dim b: b = 2.0 WScript.StdOut.Write Right(Space(12) & FormatNumber(bisection(a, b), 10, -1, 0, 0), 12) & vbNewLine Private Function bisection(ByVal a, ByVal b) Dim c Do While(True) '区間 (a, b) の中点 c = (a + b) / 2 c = (a + b) / 2 WScript.StdOut.Write Right(Space(12) & FormatNumber(c, 10, -1, 0, 0), 12) & vbTab WScript.StdOut.Write Right(Space(12) & FormatNumber(c - Sqr(2), 10, -1, 0, 0), 12) & vbNewLine Dim fc: fc = f(c) If Abs(fc) < 0.0000000001 Then Exit Do If fc < 0 Then 'f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c Else 'f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c End If Loop bisection = c End Function Private Function f(ByVal x) f = x * x - 2.0 End Function
Option Explicit Dim a: a = 1.0 Dim b: b = 2.0 WScript.StdOut.Write Right(Space(12) & FormatNumber(falseposition(a, b), 10, -1, 0, 0), 12) & vbNewLine Private Function falseposition(ByVal a, ByVal b) Dim c Do While(True) '点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c = (a * f(b) - b * f(a)) / (f(b) - f(a)) WScript.StdOut.Write Right(Space(12) & FormatNumber(c, 10, -1, 0, 0), 12) & vbTab WScript.StdOut.Write Right(Space(12) & FormatNumber(c - Sqr(2), 10, -1, 0, 0), 12) & vbNewLine Dim fc: fc = f(c) If Abs(fc) < 0.0000000001 Then Exit Do If fc < 0 Then 'f(c) < 0 であれば, 解は区間 (c, b) の中に存在 a = c Else 'f(c) > 0 であれば, 解は区間 (a, c) の中に存在 b = c End If Loop falseposition = c End Function Private Function f(ByVal x) f = x * x - 2.0 End Function
Option Explicit Dim x: x = 1.0 WScript.StdOut.Write Right(Space(12) & FormatNumber(iterative(x), 10, -1, 0, 0), 12) & vbNewLine Private Function iterative(ByVal x0) Dim x1 Do While(True) x1 = g(x0) WScript.StdOut.Write Right(Space(12) & FormatNumber(x1, 10, -1, 0, 0), 12) & vbTab WScript.StdOut.Write Right(Space(12) & FormatNumber(x1 - Sqr(2), 10, -1, 0, 0), 12) & vbNewLine If Abs(x1 - x0) < 0.0000000001 Then Exit Do x0 = x1 Loop iterative = x1 End Function Private Function g(ByVal x) g = (x / 2) + (1 / x) End Function
Option Explicit Dim x: x = 2.0 WScript.StdOut.Write Right(Space(12) & FormatNumber(newton(x), 10, -1, 0, 0), 12) & vbNewLine Private Function newton(ByVal x0) Dim x1 Do While(True) x1 = x0 - (f0(x0) / f1(x0)) WScript.StdOut.Write Right(Space(12) & FormatNumber(x1, 10, -1, 0, 0), 12) & vbTab WScript.StdOut.Write Right(Space(12) & FormatNumber(x1 - Sqr(2), 10, -1, 0, 0), 12) & vbNewLine If Abs(x1 - x0) < 0.0000000001 Then Exit Do x0 = x1 Loop newton = x1 End Function Private Function f0(ByVal x) f0 = x * x - 2 End Function Private Function f1(ByVal x) f1 = 2 * x End Function
Option Explicit Dim x: x = 2.0 WScript.StdOut.Write Right(Space(12) & FormatNumber(bailey(x), 10, -1, 0, 0), 12) & vbNewLine Private Function bailey(ByVal x0) Dim x1 Do While(True) x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0))))) WScript.StdOut.Write Right(Space(12) & FormatNumber(x1, 10, -1, 0, 0), 12) & vbTab WScript.StdOut.Write Right(Space(12) & FormatNumber(x1 - Sqr(2), 10, -1, 0, 0), 12) & vbNewLine If Abs(x1 - x0) < 0.0000000001 Then Exit Do x0 = x1 Loop bailey = x1 End Function Private Function f0(ByVal x) f0 = x * x - 2 End Function Private Function f1(ByVal x) f1 = 2 * x End Function Private Function f2(ByVal x) f2 = 2 End Function
Option Explicit Dim x0: x0 = 1.0 Dim x1: x1 = 2.0 WScript.StdOut.Write Right(Space(12) & FormatNumber(secant(x0, x1), 10, -1, 0, 0), 12) & vbNewLine Private Function secant(ByVal x0, ByVal x1) Dim x2 Do While(True) x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) WScript.StdOut.Write Right(Space(12) & FormatNumber(x2, 10, -1, 0, 0), 12) & vbTab WScript.StdOut.Write Right(Space(12) & FormatNumber(x2 - Sqr(2), 10, -1, 0, 0), 12) & vbNewLine If Abs(x2 - x1) < 0.0000000001 Then Exit Do x0 = x1 x1 = x2 Loop secant = x2 End Function Private Function f(ByVal x) f = x * x - 2 End Function
Option Explicit Private Const N = 4 Private a: a = Array(Array(9,2,1,1),Array(2,8,-2,1),Array(-1,-2,7,-2),Array(1,-1,-2,6)) Private b: b = Array(20,16,8,17) Private c: c = Array(0,0,0,0) 'ヤコビの反復法 jacobi a,b,c WScript.StdOut.WriteLine "X" disp_vector c 'ヤコビの反復法 Private Sub jacobi(ByVal a, ByVal b, ByRef x0) Do While(True) Dim x1: x1 = Array(0,0,0,0) Dim finish: finish = True Dim i, j For i = 0 To (N - 1) x1(i) = 0 For j = 0 To (N - 1) If i <> j Then x1(i) = x1(i) + a(i)(j) * x0(j) End If Next x1(i) = (b(i) - x1(i)) / a(i)(i) If (Abs(x1(i) - x0(i)) > 0.0000000001) Then finish = False End If Next For i = 0 To (N - 1) x0(i) = x1(i) Next If (finish) Then Exit Sub End If disp_vector x0 Loop End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub
Option Explicit Private Const N = 4 Private a: a = Array(Array(9,2,1,1),Array(2,8,-2,1),Array(-1,-2,7,-2),Array(1,-1,-2,6)) Private b: b = Array(20,16,8,17) Private c: c = Array(0,0,0,0) 'ガウス・ザイデル法 gauss a,b,c WScript.StdOut.WriteLine "X" disp_vector c 'ガウス・ザイデル法 Private Sub gauss(ByVal a, ByVal b, ByRef x0) Do While(True) Dim finish: finish = True Dim i, j For i = 0 To (N - 1) Dim x1: x1 = 0 For j = 0 To (N - 1) If i <> j Then x1 = x1 + a(i)(j) * x0(j) End If Next x1 = (b(i) - x1) / a(i)(i) If (Abs(x1 - x0(i)) > 0.0000000001) Then finish = False End If x0(i) = x1 Next If (finish) Then Exit Sub End If disp_vector x0 Loop End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub
Option Explicit Private Const N = 4 Private a: a = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1)) Private b: b = Array(8,17,20,16) 'ピボット選択 pivoting a,b WScript.StdOut.WriteLine "pivoting" WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '前進消去 forward_elimination a,b WScript.StdOut.WriteLine "forward elimination" WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '後退代入 backward_substitution a,b WScript.StdOut.WriteLine "X" disp_vector b '前進消去 Private Sub forward_elimination(ByRef a, ByRef b) Dim pivot, row, col For pivot = 0 To (N - 2) For row = (pivot + 1) To (N - 1) Dim s: s = a(row)(pivot) / a(pivot)(pivot) For col = pivot To (N - 1) a(row)(col) = a(row)(col) - a(pivot)(col) * s Next b(row) = b(row) - b(pivot) * s Next Next End Sub '後退代入 Private Sub backward_substitution(ByVal a, ByRef b) Dim row, col For row = (N - 1) To 0 Step -1 For col = (N - 1) To (row + 1) Step - 1 b(row) = b(row) - a(row)(col) * b(col) Next b(row) = b(row) / a(row)(row) Next End Sub 'ピボット選択 Private Sub pivoting(ByRef a, ByRef b) Dim pivot, row, col For pivot = 0 To (N - 1) '各列で 一番値が大きい行を 探す Dim max_row: max_row = pivot Dim max_val: max_val = 0 For row = pivot To (N - 1) If Abs(a(row)(pivot)) > max_val Then '一番値が大きい行 max_val = Abs(a(row)(pivot)) max_row = row End If Next '一番値が大きい行と入れ替え If max_row <> pivot Then Dim tmp For col = 0 To (N - 1) tmp = a(max_row)(col) a(max_row)(col) = a(pivot)(col) a(pivot)(col) = tmp Next tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp End If Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '2次元配列を表示 Private Sub disp_matrix(ByVal matrix) Dim row, col For Each row In matrix For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine Next End Sub
Option Explicit Private Const N = 4 Private a: a = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1)) Private b: b = Array(8,17,20,16) 'ピボット選択 pivoting a,b WScript.StdOut.WriteLine "pivoting" WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '前進消去 forward_elimination a,b WScript.StdOut.WriteLine "forward elimination" WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '後退代入 backward_substitution a,b WScript.StdOut.WriteLine "X" disp_vector b '前進消去 Private Sub forward_elimination(ByRef a, ByRef b) Dim pivot, row, col For pivot = 0 To (N - 1) For row = 0 To (N - 1) If row <> pivot Then Dim s: s = a(row)(pivot) / a(pivot)(pivot) For col = pivot To (N - 1) a(row)(col) = a(row)(col) - a(pivot)(col) * s Next b(row) = b(row) - b(pivot) * s End If Next Next End Sub '後退代入 Private Sub backward_substitution(ByVal a, ByRef b) Dim pivot For pivot = 0 To (N -1) b(pivot) = b(pivot) / a(pivot)(pivot) Next End Sub 'ピボット選択 Private Sub pivoting(ByRef a, ByRef b) Dim pivot, row, col For pivot = 0 To (N - 1) '各列で 一番値が大きい行を 探す Dim max_row: max_row = pivot Dim max_val: max_val = 0 For row = pivot To (N - 1) If Abs(a(row)(pivot)) > max_val Then '一番値が大きい行 max_val = Abs(a(row)(pivot)) max_row = row End If Next '一番値が大きい行と入れ替え If max_row <> pivot Then Dim tmp For col = 0 To (N - 1) tmp = a(max_row)(col) a(max_row)(col) = a(pivot)(col) a(pivot)(col) = tmp Next tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp End If Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '2次元配列を表示 Private Sub disp_matrix(ByVal matrix) Dim row, col For Each row In matrix For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine Next End Sub
Option Explicit Private Const N = 4 Private a: a = Array(Array(-1,-2,7,-2),Array(1,-1,-2,6),Array(9,2,1,1),Array(2,8,-2,1)) Private b: b = Array(8,17,20,16) 'ピボット選択 pivoting a,b WScript.StdOut.WriteLine "pivoting" WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '前進消去 forward_elimination a,b WScript.StdOut.WriteLine "LU" disp_matrix a 'Ly=b から y を求める (前進代入) Private y: y = Array(0.0,0.0,0.0,0.0) forward_substitution a,b,y WScript.StdOut.WriteLine "Y" disp_vector y 'Ux=y から x を求める (後退代入) Private x: x = Array(0.0,0.0,0.0,0.0) backward_substitution a,y,x WScript.StdOut.WriteLine "X" disp_vector x '前進消去 Private Sub forward_elimination(ByRef a, ByVal b) Dim pivot, row, col For pivot = 0 To (N - 2) For row = (pivot + 1) To (N - 1) Dim s: s = a(row)(pivot) / a(pivot)(pivot) For col = pivot To (N - 1) a(row)(col) = a(row)(col) - a(pivot)(col) * s 'これが 上三角行列 Next a(row)(pivot) = s 'これが 下三角行列 'b(row) = b(row) - b(pivot) * s 'この値は変更しない Next Next End Sub '前進代入 Private Sub forward_substitution(ByVal a, ByVal b, ByRef y) Dim row, col For row = 0 To (N - 1) For col = 0 To row b(row) = b(row) - a(row)(col) * y(col) Next y(row) = b(row) Next End Sub '後退代入 Private Sub backward_substitution(ByVal a, ByVal y, ByRef x) Dim row, col For row = (N - 1) To 0 Step -1 For col = (N - 1) To (row + 1) Step - 1 y(row) = y(row) - a(row)(col) * x(col) Next x(row) = y(row) / a(row)(row) Next End Sub 'ピボット選択 Private Sub pivoting(ByRef a, ByRef b) Dim pivot, row, col For pivot = 0 To (N - 1) '各列で 一番値が大きい行を 探す Dim max_row: max_row = pivot Dim max_val: max_val = 0 For row = pivot To (N - 1) If Abs(a(row)(pivot)) > max_val Then '一番値が大きい行 max_val = Abs(a(row)(pivot)) max_row = row End If Next '一番値が大きい行と入れ替え If max_row <> pivot Then Dim tmp For col = 0 To (N - 1) tmp = a(max_row)(col) a(max_row)(col) = a(pivot)(col) a(pivot)(col) = tmp Next tmp = b(max_row) b(max_row) = b(pivot) b(pivot) = tmp End If Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '2次元配列を表示 Private Sub disp_matrix(ByVal matrix) Dim row, col For Each row In matrix For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine Next End Sub
Option Explicit Private Const N = 4 Private a: a = Array(Array(5.0,2.0,3.0,4.0),Array(2.0,10.0,6.0,7.0),Array(3.0,6.0,15.0,9.0),Array(4.0,7.0,9.0,20.0)) Private b: b = Array(34.0,68.0,96.0,125.0) WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '前進消去 forward_elimination a,b WScript.StdOut.WriteLine "LL^T" disp_matrix a 'Ly=b から y を求める (前進代入) Private y: y = Array(0.0,0.0,0.0,0.0) forward_substitution a,b,y WScript.StdOut.WriteLine "Y" disp_vector y 'L^Tx=y から x を求める (後退代入) Private x: x = Array(0.0,0.0,0.0,0.0) backward_substitution a,y,x WScript.StdOut.WriteLine "X" disp_vector x '前進消去 Private Sub forward_elimination(ByRef a, ByVal b) Dim pivot, row, col For pivot = 0 To (N - 1) Dim s: s = 0 For col = 0 To (pivot - 1) s = s + a(pivot)(col) * a(pivot)(col) Next 'ここで根号の中が負の値になると計算できない! a(pivot)(pivot) = Sqr(a(pivot)(pivot) - s) For row = (pivot + 1) To (N - 1) s = 0 For col = 0 To (pivot - 1) s = s + a(row)(col) * a(pivot)(col) Next a(row)(pivot) = (a(row)(pivot) - s) / a(pivot)(pivot) a(pivot)(row) = a(row)(pivot) Next Next End Sub '前進代入 Private Sub forward_substitution(ByVal a, ByVal b, ByRef y) Dim row, col For row = 0 To (N - 1) For col = 0 To row b(row) = b(row) - a(row)(col) * y(col) Next y(row) = b(row) / a(row)(row) Next End Sub '後退代入 Private Sub backward_substitution(ByVal a, ByVal y, ByRef x) Dim row, col For row = (N - 1) To 0 Step -1 For col = (N - 1) To (row + 1) Step - 1 y(row) = y(row) - a(row)(col) * x(col) Next x(row) = y(row) / a(row)(row) Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '2次元配列を表示 Private Sub disp_matrix(ByVal matrix) Dim row, col For Each row In matrix For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine Next End Sub
Option Explicit Private Const N = 4 Private a: a = Array(Array(5.0,2.0,3.0,4.0),Array(2.0,10.0,6.0,7.0),Array(3.0,6.0,15.0,9.0),Array(4.0,7.0,9.0,20.0)) Private b: b = Array(34.0,68.0,96.0,125.0) WScript.StdOut.WriteLine "A" disp_matrix a WScript.StdOut.WriteLine "B" disp_vector b WScript.StdOut.WriteLine '前進消去 forward_elimination a,b WScript.StdOut.WriteLine "LDL^T" disp_matrix a 'Ly=b から y を求める (前進代入) Private y: y = Array(0.0,0.0,0.0,0.0) forward_substitution a,b,y WScript.StdOut.WriteLine "Y" disp_vector y 'DL^Tx=y から x を求める (後退代入) Private x: x = Array(0.0,0.0,0.0,0.0) backward_substitution a,y,x WScript.StdOut.WriteLine "X" disp_vector x '前進消去 Private Sub forward_elimination(ByRef a, ByVal b) Dim pivot, row, col, k For pivot = 0 To (N - 1) 'pivot < k の場合 Dim s: s = 0 For col = 0 To (pivot - 1) s = a(pivot)(col) For k = 0 To (col - 1) s = s - a(pivot)(k) * a(col)(k) * a(k)(k) Next a(pivot)(col) = s / a(col)(col) a(col)(pivot) = a(pivot)(col) Next 'pivot == k の場合 s = a(pivot)(pivot) For k = 0 To (pivot - 1) s = s - a(pivot)(k) * a(pivot)(k) * a(k)(k) Next a(pivot)(pivot) = s Next End Sub '前進代入 Private Sub forward_substitution(ByVal a, ByVal b, ByRef y) Dim row, col For row = 0 To (N - 1) For col = 0 To row b(row) = b(row) - a(row)(col) * y(col) Next y(row) = b(row) Next End Sub '後退代入 Private Sub backward_substitution(ByVal a, ByVal y, ByRef x) Dim row, col For row = (N - 1) To 0 Step -1 For col = (N - 1) To (row + 1) Step - 1 y(row) = y(row) - a(row)(col) * a(row)(row) * x(col) Next x(row) = y(row) / a(row)(row) Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '2次元配列を表示 Private Sub disp_matrix(ByVal matrix) Dim row, col For Each row In matrix For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine Next End Sub
Option Explicit Private Const N = 3 'ベキ乗法で最大固有値を求める Call Main Private Sub Main() Dim a: a = 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)) Dim x: x = Array(1.0 ,0.0 ,0.0 ,0.0) 'ベキ乗法 Dim lambda: lambda = power(a, x) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" WScript.StdOut.WriteLine Right(Space(14) & FormatNumber(lambda, 10, -1, 0, 0), 14) WScript.StdOut.WriteLine "eigenvector" Call disp_vector(x) End Sub 'ベキ乗法 Private Function power(ByVal a, ByRef x0) Dim i, j, k Dim lambda: lambda = 0.0 '正規化 (ベクトル x0 の長さを1にする) Call normarize(x0) Dim e0: e0 = 0.0 For i = 0 To N e0 = e0 + x0(i) Next For k = 1 To 200 '1次元配列を表示 WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab Call disp_vector(x0) '行列の積 x1 = A × x0 Dim x1: x1 = Array(0.0, 0.0, 0.0, 0.0) For i = 0 To N For j = 0 To N x1(i) = x1(i) + a(i)(j) * x0(j) Next Next '内積 Dim p0: p0 = 0.0 Dim p1: p1 = 0.0 For i = 0 To N p0 = p0 + x1(i) * x1(i) p1 = p1 + x1(i) * x0(i) Next '固有値 lambda = p0 / p1 '正規化 (ベクトル x1 の長さを1にする) Call normarize(x1) '収束判定 Dim e1: e1 = 0.0 For i = 0 To N e1 = e1 + x1(i) Next If Abs(e1 - e0) < 0.00000000001 Then Exit For x0 = x1 e0 = e1 Next power = lambda End Function '1次元配列を表示 Private Sub disp_vector(ByVal row) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '正規化 (ベクトルの長さを1にする) Private Sub normarize(ByRef x()) Dim s: s = 0.0 Dim i For i = 0 To N s = s + x(i) * x(i) Next s = Sqr(s) For i = 0 To N x(i) = x(i) / s Next End Sub
Option Explicit Private Const N = 3 '逆ベキ乗法で最小固有値を求める Call Main Private Sub Main() Dim a: a = 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)) Dim x: x = Array(1.0 ,0.0 ,0.0 ,0.0) 'LU分解 Call forward_elimination(a) '逆ベキ乗法 Dim lambda: lambda = inverse(a, x) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" WScript.StdOut.WriteLine Right(Space(14) & FormatNumber(lambda, 10, -1, 0, 0), 14) WScript.StdOut.WriteLine "eigenvector" Call disp_vector(x) End Sub '逆ベキ乗法 Private Function inverse(ByVal a, ByRef x0) Dim i, j, k Dim lambda: lambda = 0.0 '正規化 (ベクトル x0 の長さを1にする) Call normarize(x0) Dim e0: e0 = 0.0 For i = 0 To N e0 = e0 + x0(i) Next For k = 1 To 200 '1次元配列を表示 WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab Call disp_vector(x0) 'Ly = b から y を求める (前進代入) Dim y: y = Array(0.0, 0.0, 0.0, 0.0) Call forward_substitution(a,y,x0) 'Ux = y から x を求める (後退代入) Dim x1: x1 = Array(0.0, 0.0, 0.0, 0.0) Call backward_substitution(a,x1,y) '内積 Dim p0: p0 = 0.0 Dim p1: p1 = 0.0 For i = 0 To N p0 = p0 + x1(i) * x1(i) p1 = p1 + x1(i) * x0(i) Next '固有値 lambda = p1 / p0 '正規化 (ベクトル x1 の長さを1にする) Call normarize(x1) '収束判定 Dim e1: e1 = 0.0 For i = 0 To N e1 = e1 + x1(i) Next If Abs(e1 - e0) < 0.00000000001 Then Exit For x0 = x1 e0 = e1 Next inverse = lambda End Function 'LU分解 Private Sub forward_elimination(ByRef a) Dim pivot, row, col For pivot = 0 To (N - 1) For row = (pivot + 1) To N Dim s: s = a(row)(pivot) / a(pivot)(pivot) For col = pivot To N a(row)(col) = a(row)(col) - a(pivot)(col) * s 'これが 上三角行列 Next a(row)(pivot) = s 'これが 下三角行列 Next Next End Sub '前進代入 (Ly = b から y を求める) Private Sub forward_substitution(ByVal a, ByRef y, ByVal b) Dim row, col For row = 0 To N For col = 0 To row b(row) = b(row) - a(row)(col) * y(col) Next y(row) = b(row) Next End Sub '後退代入 (Ux = y から x を求める) Private Sub backward_substitution(ByVal a, ByRef x, ByVal y) Dim row, col For row = N To 0 Step -1 For col = N To (row + 1) Step - 1 y(row) = y(row) - a(row)(col) * x(col) Next x(row) = y(row) / a(row)(row) Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '正規化 (ベクトルの長さを1にする) Private Sub normarize(ByRef x()) Dim s: s = 0.0 Dim i For i = 0 To N s = s + x(i) * x(i) Next s = Sqr(s) For i = 0 To N x(i) = x(i) / s Next End Sub
Option Explicit Private Const N = 3 'LR分解で固有値を求める Call Main Private Sub Main() Dim a: a = 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)) Dim l: l = Array(_ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0)) Dim u: u = Array(_ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0)) Dim i, j, k For k = 1 To 200 'LU分解 Call decomp(a, l, u) '行列の積 Call multiply(u, l, a) '対角要素を表示 WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab Call disp_eigenvalue(a) '収束判定 Dim e: e = 0.0 For i = 1 To N For j = 0 To i - 1 e = e + Abs(a(i)(j)) Next Next If e < 0.00000000001 Then Exit For Next WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" Call disp_eigenvalue(a) End Sub '対角要素を表示 Private Sub disp_eigenvalue(ByVal a) Dim i For i = 0 To N WScript.StdOut.Write Right(Space(14) & FormatNumber(a(i)(i), 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub 'LU分解 Private Sub decomp(ByRef a, ByRef l, ByRef u) Dim i For i = 0 To N Dim j For j = 0 To N l(i)(j) = 0.0 u(i)(j) = 0.0 Next Next l(0)(0) = 1.0 For j = 0 To N u(0)(j) = a(0)(j) Next For i = 1 To N u(i)(0) = 0.0 l(0)(i) = 0.0 l(i)(0) = a(i)(0) / u(0)(0) Next For i = 1 To N l(i)(i) = 1.0 Dim t: t = a(i)(i) Dim k For k = 0 To i t = t - l(i)(k) * u(k)(i) Next u(i)(i) = t For j = i + 1 To N u(j)(i) = 0.0 l(i)(j) = 0.0 t = a(j)(i) For k = 0 To i t = t - l(j)(k) * u(k)(i) Next l(j)(i) = t / u(i)(i) t = a(i)(j) For k = 0 To i t = t - l(i)(k) * u(k)(j) Next u(i)(j) = t Next Next End Sub '行列の積 Private Sub multiply(ByRef a, ByRef b, ByRef c) Dim i For i = 0 To N Dim j For j = 0 To N Dim s: s = 0.0 Dim k For k = 0 To N s = s + a(i)(k) * b(k)(j) Next c(i)(j) = s Next Next End Sub
Option Explicit Private Const N = 3 'QR分解で固有値を求める Call Main Private Sub Main() Dim a: a = 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)) Dim q: q = Array(_ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0)) Dim r: r = Array(_ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0), _ Array(0.0, 0.0, 0.0, 0.0)) Dim i, j, k For k = 1 To 200 'QR分解 Call decomp(a, q, r) '行列の積 Call multiply(r, q, a) '対角要素を表示 WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab Call disp_eigenvalue(a) '収束判定 Dim e: e = 0.0 For i = 1 To N For j = 0 To i - 1 e = e + Abs(a(i)(j)) Next Next If e < 0.00000000001 Then Exit For Next WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" Call disp_eigenvalue(a) End Sub '対角要素を表示 Private Sub disp_eigenvalue(ByVal a) Dim i For i = 0 To N WScript.StdOut.Write Right(Space(14) & FormatNumber(a(i)(i), 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub 'QR分解 Private Sub decomp(ByRef a, ByRef q, ByRef r) Dim x: x = Array(0.0, 0.0, 0.0, 0.0) Dim k For k = 0 To N Dim i For i = 0 To N x(i) = a(i)( k) Next Dim j For j = 0 To k - 1 Dim t: t = 0 For i = 0 To N t = t + a(i)(k) * q(i)(j) Next r(j)(k) = t r(k)(j) = 0 For i = 0 To N x(i) = x(i) - t * q(i)(j) Next Next t = 0 For i = 0 To N t = t + x(i) * x(i) Next r(k)(k) = Sqr(t) For i = 0 To N q(i)(k) = x(i) / r(k)(k) Next Next End Sub '行列の積 Private Sub multiply(ByRef a, ByRef b, ByRef c) Dim i For i = 0 To N Dim j For j = 0 To N Dim s: s = 0.0 Dim k For k = 0 To N s = s + a(i)(k) * b(k)(j) Next c(i)(j) = s Next Next End Sub
Option Explicit Private Const N = 3 Private Const PI = 3.14159265358979 'ヤコビ法で固有値を求める Call Main Private Sub Main() Dim a: a = 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)) Dim v: v = 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)) 'ヤコビ法 Call jacobi(a, v) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" Call disp_eigenvalue(a) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvector" Call disp_eigenvector(v) End Sub 'ヤコビ法 Private Sub jacobi(ByRef a, ByRef v) Dim i, j, k Dim max_val Dim p, q For k = 1 To 100 '最大値を探す max_val = 0.0 For i = 0 To N For j = i + 1 To N If (max_val < Abs(a(i)(j))) Then max_val = Abs(a(i)(j)) p = i q = j End If Next Next 'θ を求める Dim t If Abs(a(p)(p) - a(q)(q)) < 0.00000000001 Then 'a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = PI / 4.0 If (a(p)(p) < 0) Then t = -t End If Else 'a_{pp} ≠ a_{qq} のとき t = Atn(2.0 * a(p)(q) / (a(p)(p) - a(q)(q))) / 2.0 End If 'θ を使って 行列 U を作成し、A = U^t × A × U Dim c: c = Cos(t) Dim s: s = Sin(t) 'U^t × A Dim t1, t2 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 Next '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 Next '行列の対角要素を表示 (固有値) WScript.StdOut.Write Right(Space(3) & k, 3) & vbTab Call disp_eigenvalue(a) '収束判定 If max_val < 0.00000000001 Then Exit For Next End Sub '行列の対角要素を表示 Private Sub disp_eigenvalue(ByVal x) Dim i For i = 0 To N WScript.StdOut.Write Right(Space(14) & FormatNumber(x(i)(i), 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '固有ベクトルを表示 Private Sub disp_eigenvector(ByVal matrix) Dim row, col For Each row In matrix normarize(row) disp_vector(row) WScript.StdOut.WriteLine Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next End Sub '正規化 (ベクトルの長さを1にする) Private Sub normarize(ByRef x()) Dim s: s = 0.0 Dim i For i = 0 To N s = s + x(i) * x(i) Next s = Sqr(s) For i = 0 To N x(i) = x(i) / s Next End Sub
Option Explicit Private Const N = 3 'ハウスホルダー変換とQR分解で固有値を求める Call Main Private Sub Main() Dim a: a = 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)) Dim d: d = Array(0.0 ,0.0 ,0.0 ,0.0) Dim e: e = Array(0.0 ,0.0 ,0.0 ,0.0) 'ハウスホルダー変換 Call tridiagonalize(a, d, e) 'QR分解 Call decomp(a, d, e) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvalue" Call disp_vector(d) WScript.StdOut.WriteLine WScript.StdOut.WriteLine "eigenvector" Call disp_matrix(a) End Sub 'ハウスホルダー変換 Private Sub tridiagonalize(ByRef a, ByRef d, ByRef e) Dim v: v = Array(0.0 ,0.0 ,0.0 ,0.0) Dim i, j, k For k = 0 To N - 2 d(k) = a(k)(k) Dim t: t = 0.0 Dim w: w = Array(0.0, 0.0, 0.0, 0.0) For i = k + 1 To N w(i) = a(i)(k) t = t + w(i) * w(i) Next Dim s: s = Sqr(t) If w(k + 1) < 0 Then s = -s End If If Abs(s) < 0.00000000001 Then e(k + 1) = 0.0 Else e(k + 1) = -s w(k + 1) = w(k + 1) + s s = 1 / Sqr(w(k + 1) * s) For i = k + 1 To N w(i) = w(i) * s Next For i = k + 1 To N s = 0.0 For j = k + 1 To N If j <= i Then s = s + a(i)(j) * w(j) Else s = s + a(j)(i) * w(j) End If Next v(i) = s Next s = 0 For i = k + 1 To N s = s + w(i) * v(i) Next s = s / 2 For i = k + 1 To N v(i) = v(i) - s * w(i) Next For i = k + 1 To N For j = k + 1 To i a(i)(j) = a(i)(j) - w(i) * v(j) - w(j) * v(i) Next Next '次の行は固有ベクトルを求めないなら不要 For i = k + 1 To N a(i)(k) = w(i) Next End If Next 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 Step -1 If k < N - 1 Then For i = k + 1 To N w(i) = a(i)(k) Next For i = k + 1 To N s = 0.0 For j = k + 1 To N s = s + a(i)(j) * w(j) Next v(i) = s Next For i = k + 1 To N For j = k + 1 To N a(i)(j) = a(i)(j) - v(i) * w(j) Next Next End If For i = 0 To N a(i)(k) = 0.0 Next a(k)(k) = 1.0 Next End Sub 'QR分解 Private Sub decomp(ByRef a, ByRef d, ByRef e) Dim i, j, k e(0) = 1.0 Dim h: h = N Do While (Abs(e(h)) < 0.00000000001) h = h - 1 Loop Do While h > 0 e(0) = 0.0 Dim l: l = h - 1 Do While (Abs(e(l)) >= 0.00000000001) l = l - 1 Loop For j = 1 To 100 Dim w: w = (d(h - 1) - d(h)) / 2.0 Dim s: s = Sqr(w * w + e(h) * e(h)) If w < 0.0 Then s = - s End If Dim x: x = d(l) - d(h) + e(h) * e(h) / (w + s) Dim y: y = e(l + 1) Dim z: z = 0.0 Dim t, u For k = l To h - 1 If Abs(x) >= Abs(y) Then t = -y / x u = 1 / Sqr(t * t + 1.0) s = t * u Else t = -x / y s = 1 / Sqr(t * t + 1.0) If t < 0 Then s = -s End If u = t * s End If 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 Next If k < N - 1 Then x = e(k + 1) y = -s * e(k + 2) z = y e(k + 2) = u * e(k + 2) End If Next WScript.StdOut.Write Right(Space(3) & j, 3) & vbTab Call disp_vector(d) '収束判定 If (Abs(e(h)) < 0.00000000001) Then Exit For Next e(0) = 1.0 Do While (Abs(e(h)) < 0.00000000001) h = h - 1 Loop Loop '次の行は固有ベクトルを求めないなら不要 For k = 0 To N - 1 l = k For i = k + 1 To N If d(i) > d(l) Then l = i End If Next 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 Next Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row()) Dim col For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine End Sub '2次元配列を表示 Private Sub disp_matrix(ByVal matrix()) Dim row, col For Each row In matrix For Each col In row WScript.StdOut.Write Right(Space(14) & FormatNumber(col, 10, -1, 0, 0), 14) & vbTab Next WScript.StdOut.WriteLine Next End Sub