Module VB0101 Public Sub Main() Console.WriteLine(3 + 5) Console.WriteLine(3 - 5) Console.WriteLine(3 * 5) Console.WriteLine(3 ^ 5) Console.WriteLine(5 / 3) Console.WriteLine(5 \ 3) Console.WriteLine(5 Mod 3) Console.Write(3 * 5 & vbNewLine) Console.WriteLine(String.Format("{0,3:D}", 3 * 5)) Console.WriteLine(String.Format("{0,23:F20}", 5 / 3)) End Sub End Module
Module VB0102 Public Sub Main() Dim i As Integer = 3 * 5 Console.WriteLine("3 * 5 = " & i) Console.WriteLine(String.Format("3 * 5 = {0}", i)) End Sub End Module
Module VB0103 Public Sub Main() For i As Integer = 1 To 9 Console.Write(i & ", ") Next Console.WriteLine() End Sub End Module
Module VB0104 Public Sub Main() For i As Integer = 1 To 9 If i Mod 3 = 0 Then Console.Write(i & ", ") End If Next Console.WriteLine() End Sub End Module
Module VB0105 Public Sub Main() Dim sum As Integer = 0 For i As Integer = 1 To 99 If i Mod 3 = 0 Then sum += i End If Next Console.WriteLine(sum) End Sub End Module
Module VB0301 Public Sub Main() ' 3 の倍数の合計 Console.WriteLine( sn(3, 999) ) End Sub ' 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数 Private Function sn(ByVal a As Integer, ByVal lim As Integer) As Integer Dim n As Integer = lim \ a ' 項数:n = 上限:lim / 公差:a Dim l As Integer = n * a ' 末項:l = 項数:n * 公差:a Return (a + l) * n \ 2 ' 総和:sn = (初項:a + 末項:l) * 項数:n / 2 End Function End Module
Module VB0302 Sub Main() ' 10000 までの 自然数の和 ' 項目数 n = 10000 Dim n As Integer = 10000 Console.WriteLine( n * (n + 1) \ 2 ) End Sub End Module
Module VB0303 Sub Main() ' 10000 までの 偶数の和 ' 項目数 n = 5000 Dim n As Integer = 10000 \ 2 Console.WriteLine( n * (n + 1) ) End Sub End Module
Module VB0304 Sub Main() ' 10000 までの 奇数の和 ' 項目数 n = 5000 Dim n As Integer = 10000 \ 2 Console.WriteLine( n ^ 2 ) End Sub End Module
Module VB0305 Sub Main() ' 1000 までの 自然数の2乗の和 Dim n As Integer = 1000 Console.WriteLine( n * (n + 1) * (2 * n + 1) \ 6 ) End Sub End Module
Module VB0306 Sub Main() ' 100 までの 自然数の3乗の和 Dim n As Integer = 100 Console.WriteLine( (n ^ 2) * ((n + 1) ^ 2) \ 4 ) End Sub End Module
Module VB0307 Sub Main() ' 初項 2, 公比 3, 項数 10 の等比数列の和 Dim n As Integer = 10 Dim a As Integer = 2 Dim r As Integer = 3 Console.WriteLine( (a * ((r ^ n) - 1)) \ (r - 1) ) End Sub End Module
Option Explicit Module VB0401 Sub Main() Dim a As Integer = 5 '初項 5 Dim d As Integer = 3 '公差 3 Dim n As Integer = 10 '項数 10 Dim p As Long = 1 '積 For i As Integer = 1 To n Dim m As Integer = a + (d * (i - 1)) p *= m Next Console.WriteLine(p) End Sub End Module
Option Explicit Module VB0402 Sub Main() '初項 5, 公差 3, 項数 10 の数列の積を表示する Console.WriteLine(product(5, 3, 10)) End Sub Private Function product(ByVal m As Integer, ByVal d As Integer, ByVal n As Integer) As Long If n = 0 Then Return 1 Else Return m * product(m + d, d, n - 1) End If End Function End Module
Module VB0403 '階乗を求める関数 Private Function Fact(ByVal n As Integer) As Integer If n <= 1 Then Return 1 Else Return n * Fact(n - 1) End If End Function Sub Main() '10の階乗 Console.WriteLine(Fact(10)) Console.WriteLine(10 * 9 * 8 * 7 * 6 * 5 * 4 * 3 * 2 * 1) End Sub End Module
Module VB0404 '下降階乗冪 Private Function FallingFact(ByVal x As Integer, ByVal n As Integer) As Integer If n <= 1 Then Return x Else Return x * FallingFact(x - 1, n - 1) End If End Function Sub Main() '10 から 6 までの 総乗 Console.WriteLine(FallingFact(10, 5)) Console.WriteLine(10 * 9 * 8 * 7 * 6) End Sub End Module
Module VB0405 '上昇階乗冪 Private Function RisingFact(ByVal x As Integer, ByVal n As Integer) As Integer If n <= 1 Then Return x Else Return x * RisingFact(x + 1, n - 1) End If End Function Sub Main() '10 から 14 までの 総乗 Console.WriteLine(RisingFact(10, 5)) Console.WriteLine(10 * 11 * 12 * 13 * 14) End Sub End Module
Module VB0406 '階乗 Private Function Fact(ByVal n As Integer) As Integer If n <= 1 Then Return 1 Else Return n * Fact(n - 1) End If End Function '下降階乗冪 Private Function FallingFact(ByVal x As Integer, ByVal n As Integer) As Integer If n <= 1 Then Return x Else Return x * FallingFact(x - 1, n - 1) End If End Function Sub Main() '順列 (異なる 10 個のものから 5 個取ってできる順列の総数) Dim n As Integer = 10 Dim r As Integer = 5 Console.WriteLine(Fact(n) / Fact(n - r)) Console.WriteLine(FallingFact(n, r)) End Sub End Module
Module VB0407 Sub Main() '重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) Dim n As Integer = 10 Dim r As Integer = 5 Console.WriteLine(n ^ r) End Sub End Module
Module VB040101 '組合せ Private Function Comb(ByVal n As Integer, ByVal r As Integer) As Integer If (r = 0) OrElse (r = n) Then Return 1 ElseIf r = 1 Then Return n Else Return Comb(n - 1, r - 1) + Comb(n - 1, r) End If End Function Sub Main() '組合せ (異なる 10 個のものから 5 個取ってできる組合せの総数) Dim n As Integer = 10 Dim r As Integer = 5 Console.WriteLine(Comb(n, r)) End Sub End Module
Module VB0409 '組合せ Private Function Comb(ByVal n As Integer, ByVal r As Integer) As Integer If (r = 0) OrElse (r = n) Then Return 1 ElseIf r = 1 Then Return n Else Return Comb(n - 1, r - 1) + Comb(n - 1, r) End If End Function Sub Main() '重複組合せ (異なる 10 個のものから重複を許して 5 個とる組合せの総数) Dim n As Integer = 10 Dim r As Integer = 5 Console.WriteLine(Comb(n + r - 1, r)) End Sub End Module
Module VB0501 Public Sub Main() For degree As Integer = 0 To 360 Step 15 If (degree Mod 30 = 0 OrElse degree Mod 45 = 0) Then Dim radian As Double = degree * Math.PI / 180.0 '自作の正弦関数 Dim d1 As Double = mySin(radian, 1, False, radian, 1.0, radian) '標準の正弦関数 Dim d2 As Double = Math.Sin(radian) '標準関数との差異 Console.WriteLine(String.Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", degree, d1, d2, d1 - d2)) End If Next End Sub '自作の正弦関数 Private Function mySin(ByVal x As Double, ByVal n As Integer, ByVal nega As Boolean, ByVal numerator As Double, ByVal denominator As Double, ByVal y As Double) As Double Dim m As Integer = 2 * n denominator = denominator * (m + 1) * m numerator = numerator * x * x Dim a As Double = numerator / denominator '十分な精度になったら処理を抜ける If (a <= 0.00000000001) Then Return y Else Return y + mySin(x, n + 1, Not nega, numerator, denominator, If(nega, a, -a)) End If End Function End Module
Module VB0502 Public Sub Main() For degree As Integer = 0 To 360 Step 15 If (degree Mod 30 = 0 OrElse degree Mod 45 = 0) Then Dim radian As Double = degree * Math.PI / 180.0 '自作の余弦関数 Dim d1 As Double = myCos(radian, 1, False, 1.0, 1.0, 1.0) '標準の余弦関数 Dim d2 As Double = Math.Cos(radian) '標準関数との差異 Console.WriteLine(String.Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", degree, d1, d2, d1 - d2)) End If Next End Sub '自作の余弦関数 Private Function myCos(ByVal x As Double, ByVal n As Integer, ByVal nega As Boolean, ByVal numerator As Double, ByVal denominator As Double, ByVal y As Double) As Double Dim m As Integer = 2 * n denominator = denominator * m * (m - 1) numerator = numerator * x * x Dim a As Double = numerator / denominator '十分な精度になったら処理を抜ける If (a <= 0.00000000001) Then Return y Else Return y + myCos(x, n + 1, Not nega, numerator, denominator, If(nega, a, -a)) End If End Function End Module
Module VB0503 Public Sub Main() For i As Integer = 0 To 180 Step 15 If (i Mod 180 <> 0) Then Dim degree As Integer = i - 90 Dim radian As Double = degree * Math.PI / 180.0 Dim x2 As Double = radian * radian '自作の正接関数 Dim d1 As Double = myTan(radian, x2, 15, 0.0) '15:必要な精度が得られる十分大きな奇数 '標準の正接関数 Dim d2 As Double = Math.Tan(radian) '標準関数との差異 Console.WriteLine(String.Format("{0,3:D} : {1,13:F10} - {2,13:F10} = {3,13:F10}", degree, d1, d2, d1 - d2)) End If Next End Sub '自作の正接関数 Private Function myTan(ByVal x As Double, ByVal x2 As Double, ByVal n As Integer, ByVal t As Double) As Double t = x2 / (n - t) n -= 2 If (n <= 1) Then Return x / (1 - t) Else Return myTan(x, x2, n, t) End If End Function End Module
Module VB0504 Public Sub Main() For i As Integer = 0 To 20 Dim x As Double = (i - 10) / 4.0 '標準の指数関数 Dim d1 As Double = Math.Exp(x) '自作の指数関数 Dim d2 As Double = myExp(x, 1, 1.0, 1.0, 1.0) '標準関数との差異 Console.WriteLine(String.Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", x, d1, d2, d1 - d2)) Next End Sub '自作の指数関数 Private Function myExp(ByVal x As Double, ByVal n As Integer, ByVal numerator As Double, ByVal denominator As Double, ByVal y As Double) As Double denominator = denominator * n numerator = numerator * x Dim a As Double = numerator / denominator '十分な精度になったら処理を抜ける If (Math.Abs(a) <= 0.00000000001) Then Return y Else Return y + myExp(x, n + 1, numerator, denominator, a) End If End Function End Module
Module VB0505 Public Sub Main() For i As Integer = 0 To 20 Dim x As Double = (i - 10) / 4.0 '標準の指数関数 Dim d1 As Double = Math.Exp(x) '自作の指数関数 Dim x2 As Double = x * x Dim d2 As Double = myExp(x, x2, 30, 0.0) '30:必要な精度が得られるよう, 6から始めて4ずつ増加させる '標準関数との差異 Console.WriteLine(String.Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", x, d1, d2, d1 - d2)) Next End Sub '自作の指数関数 Private Function myExp(ByVal x As Double, ByVal x2 As Double, ByVal n As Integer, ByVal t As Double) As Double t = x2 / (n + t) n -= 4 If (n < 6) Then Return 1 + ((2 * x) / (2 - x + t)) Else Return myExp(x, x2, n, t) End If End Function End Module
Module VB0506 Public Sub Main() For i As Integer = 1 To 20 Dim x As Double = i / 5.0 '標準の対数関数 Dim d1 As Double = Math.Log(x) '自作の対数関数 Dim x2 As Double = (x - 1) / (x + 1) Dim d2 As Double = 2 * myLog(x2, x2, 1.0, x2) '標準関数との差異 Console.WriteLine(String.Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", x, d1, d2, d1 - d2)) Next End Sub '自作の対数関数 Private Function myLog(ByVal x2 As Double, ByVal numerator As Double, ByVal denominator As Double, ByVal y As Double) As Double denominator = denominator + 2 numerator = numerator * x2 * x2 Dim a As Double = numerator / denominator '十分な精度になったら処理を抜ける If (Math.Abs(a) <= 0.00000000001) Then Return y Else Return y + myLog(x2, numerator, denominator, a) End If End Function End Module
Module VB0507 Public Sub Main() For i As Integer = 1 To 20 Dim x As Double = i / 5.0 '標準の対数関数 Dim d1 As Double = Math.Log(x) '自作の対数関数 Dim d2 As Double = myLog(x - 1, 27, 0.0) '27:必要な精度が得られる十分大きな奇数 '標準関数との差異 Console.WriteLine(String.Format("{0,5:F2} : {1,13:F10} - {2,13:F10} = {3,13:F10}", x, d1, d2, d1 - d2)) Next End Sub '自作の対数関数 Private Function myLog(ByVal x As Double, ByVal n As Integer, ByVal t As Double) As Double Dim n2 As Integer = n Dim x2 As Double = 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 Return x / (1 + t) Else Return myLog(x, n - 1, t) End If End Function End Module
Module VB0508 Public Sub Main() For i As Integer = 0 To 20 Dim x As Integer = i - 10 '自作の双曲線正弦関数 Dim d1 As Double = mySinh(x, 1, x, 1.0, x) '標準の双曲線正弦関数 Dim d2 As Double = Math.Sinh(x) '標準関数との差異 Console.WriteLine(String.Format("{0,3:D} : {1,17:F10} - {2,17:F10} = {3,13:F10}", x, d1, d2, d1 - d2)) Next End Sub '自作の双曲線正弦関数 Private Function mySinh(ByVal x As Double, ByVal n As Integer, ByVal numerator As Double, ByVal denominator As Double, ByVal y As Double) As Double Dim m As Integer = 2 * n denominator = denominator * (m + 1) * m numerator = numerator * x * x Dim a As Double = numerator / denominator '十分な精度になったら処理を抜ける If (Math.Abs(a) <= 0.00000000001) Then Return y Else Return y + mySinh(x, n + 1, numerator, denominator, a) End If End Function End Module
Module VB0509 Public Sub Main() For i As Integer = 0 To 20 Dim x As Integer = i - 10 '自作の双曲線余弦関数 Dim d1 As Double = myCosh(x, 1, 1.0, 1.0, 1.0) '標準の双曲線余弦関数 Dim d2 As Double = Math.Cosh(x) '標準関数との差異 Console.WriteLine(String.Format("{0,3:D} : {1,17:F10} - {2,17:F10} = {3,13:F10}", x, d1, d2, d1 - d2)) Next End Sub '自作の双曲線余弦関数 Private Function myCosh(ByVal x As Double, ByVal n As Integer, ByVal numerator As Double, ByVal denominator As Double, ByVal y As Double) As Double Dim m As Integer = 2 * n denominator = denominator * m * (m - 1) numerator = numerator * x * x Dim a As Double = numerator / denominator '十分な精度になったら処理を抜ける If (Math.Abs(a) <= 0.00000000001) Then Return y Else Return y + myCosh(x, n + 1, numerator, denominator, a) End If End Function End Module
Module VB0601 Public Sub Main() Const a As Double = 0 Const b As Double = 1 '台形則で積分 Dim n As Integer = 2 For j As Integer = 1 To 10 Dim h As Double = (b - a) / n Dim s As Double = 0 Dim x As Double = a For i As Integer = 1 To n - 1 x += h s += f(x) Next s = h * ((f(a) + f(b)) / 2 + s) n *= 2 '結果を π と比較 Console.WriteLine(String.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI)) Next End Sub Private Function f(ByVal x As Double) As Double Return 4 / (1 + x * x) End Function End Module
Module VB0602 Public Sub Main() Const a As Double = 0 Const b As Double = 1 '中点則で積分 Dim n As Integer = 2 For j As Integer = 1 To 10 Dim h As Double = (b - a) / n Dim s As Double = 0 Dim x As Double = a + (h / 2) For i As Integer = 1 To n s += f(x) x += h Next s = h * s n *= 2 '結果を π と比較 Console.WriteLine(String.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI)) Next End Sub Private Function f(ByVal x As Double) As Double Return 4 / (1 + x * x) End Function End Module
Module VB0603 Public Sub Main() Const a As Double = 0 Const b As Double = 1 'Simpson則で積分 Dim n As Integer = 2 For j As Integer = 1 To 5 Dim h As Double = (b - a) / n Dim s2 As Double = 0 Dim s4 As Double = 0 Dim x As Double = a + h For i As Integer = 1 To n \ 2 s4 += f(x) x += h s2 += f(x) x += h Next s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 Dim s As Double = (s2 + s4) * h / 3 n *= 2 '結果を π と比較 Console.WriteLine(String.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI)) Next End Sub Private Function f(ByVal x As Double) As Double Return 4 / (1 + x * x) End Function End Module
Module VB0604 Public Sub Main() Const a As Double = 0 Const b As Double = 1 Dim t(6, 6) As Double '台形則で積分 Dim n As Integer = 2 For i As Integer = 1 To 6 Dim h As Double = (b - a) / n Dim s As Double = 0 Dim x As Double = a For j As Integer = 1 To n - 1 x += h s += f(x) Next '結果を保存 t(i,1) = h * ((f(a) + f(b)) / 2 + s) n *= 2 Next 'Richardsonの補外法 n = 4 For j As Integer = 2 To 6 For i As Integer = 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 '結果を π と比較 Console.WriteLine(String.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, t(i,j), t(i,j) - Math.PI)) End If Next n *= 4 Next End Sub Private Function f(ByVal x As Double) As Double Return 4 / (1 + x * x) End Function End Module
Module VB0701 'データ点の数 - 1 Private Const N As Integer = 6 Public Sub Main() Dim x(N) As Double Dim y(N) As Double '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット For i As Integer = 0 To N Dim d As Double = i * 1.5 - 4.5 x(i) = d y(i) = f(d) Next '0.5刻みで 与えられていない値を補間 For i As Integer = 0 To 18 Dim d As Double = i * 0.5 - 4.5 Dim d1 As Double = f(d) Dim d2 As Double = lagrange(d, x, y) '元の関数と比較 Console.WriteLine(String.Format("{0,5:F2}{4}{1,8:F5}{4}{2,8:F5}{4}{3,8:F5}", d, d1, d2, d1 - d2, vbTab)) Next End Sub '元の関数 Private Function f(ByVal x As Double) As Double Return x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function 'Lagrange (ラグランジュ) 補間 Private Function lagrange(ByVal d As Double, ByVal x() As Double, ByVal y() As Double) As Double Dim sum As Double = 0 For i As Integer = 0 To N Dim prod As Double = y(i) For j As Integer = 0 To N If j <> i Then prod *= (d - x(j)) / (x(i) - x(j)) End If Next sum += prod Next Return sum End Function End Module
Module VB0702 'データ点の数 - 1 Private Const N As Integer = 6 Public Sub Main() Dim x(N) As Double Dim y(N) As Double '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット For i As Integer = 0 To N Dim d As Double = i * 1.5 - 4.5 x(i) = d y(i) = f(d) Next '0.5刻みで 与えられていない値を補間 For i As Integer = 0 To 18 Dim d As Double = i * 0.5 - 4.5 Dim d1 As Double = f(d) Dim d2 As Double = neville(d, x, y) '元の関数と比較 Console.WriteLine(String.Format("{0,5:F2}{4}{1,8:F5}{4}{2,8:F5}{4}{3,8:F5}", d, d1, d2, d1 - d2, vbTab)) Next End Sub '元の関数 Private Function f(ByVal x As Double) As Double Return x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function 'Neville (ネヴィル) 補間 Private Function neville(ByVal d As Double, ByVal x() As Double, ByVal y() As Double) As Double Dim w(N, N) As Double For i As Integer = 0 To N w(0,i) = y(i) Next For j As Integer = 1 To N For i As Integer = 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 Return w(N,0) End Function End Module
Module VB0703 'データ点の数 - 1 Private Const N As Integer = 6 Public Sub Main() Dim x(N) As Double Dim y(N) As Double '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット For i As Integer = 0 To N Dim d1 As Double = i * 1.5 - 4.5 x(i) = d1 y(i) = f(d1) Next '差分商の表を作る Dim d(N, N) As Double For j As Integer = 0 To N d(0,j) = y(j) Next For i As Integer = 1 To N For j As Integer = 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(N) As Double For j As Integer = 0 To N a(j) = d(j,0) Next '0.5刻みで 与えられていない値を補間 For i As Integer = 0 To 18 Dim d1 As Double = i * 0.5 - 4.5 Dim d2 As Double = f(d1) Dim d3 As Double = newton(d1, x, a) '元の関数と比較 Console.WriteLine(String.Format("{0,5:F2}{4}{1,8:F5}{4}{2,8:F5}{4}{3,8:F5}", d1, d2, d3, d2 - d3, vbTab)) Next End Sub '元の関数 Private Function f(ByVal x As Double) As Double Return x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function 'Newton (ニュートン) 補間 Private Function newton(ByVal d As Double, ByVal x() As Double, ByVal a() As Double) As Double Dim sum As Double = a(0) For i As Integer = 1 To N Dim prod As Double = a(i) For j As Integer = 0 To (i - 1) If j <> i Then prod *= (d - x(j)) End If Next sum += prod Next Return sum End Function End Module
Module VB0704 'データ点の数 - 1 Private Const N As Integer = 6 Private Const Nx2 As Integer = 13 Public Sub Main() Dim x(N) As Double Dim y(N) As Double Dim yd(N) As Double '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット For i As Integer = 0 To N Dim d1 As Double = i * 1.5 - 4.5 x(i) = d1 y(i) = f(d1) yd(i) = fd(d1) Next '差分商の表を作る Dim z(Nx2) As Double Dim d(Nx2, Nx2) As Double For i As Integer = 0 To Nx2 Dim j As Integer = i \ 2 z(i) = x(j) d(0,i) = y(j) Next For i As Integer = 1 To Nx2 For j As Integer = 0 To (Nx2 - i) If i = 1 AndAlso 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(Nx2) As Double For j As Integer = 0 To Nx2 a(j) = d(j,0) Next '0.5刻みで 与えられていない値を補間 For i As Integer = 0 To 18 Dim d1 As Double = i * 0.5 - 4.5 Dim d2 As Double = f(d1) Dim d3 As Double = hermite(d1, z, a) '元の関数と比較 Console.WriteLine(String.Format("{0,5:F2}{4}{1,8:F5}{4}{2,8:F5}{4}{3,8:F5}", d1, d2, d3, d2 - d3, vbTab)) Next End Sub '元の関数 Private Function f(ByVal x As Double) As Double Return x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function '導関数 Private Function fd(ByVal x As Double) As Double Return 1 - (x ^ 2) / 2 + (x ^ 4) / (4 * 3 * 2) End Function 'Hermite (エルミート) 補間 Private Function hermite(ByVal d As Double, ByVal z() As Double, ByVal a() As Double) As Double Dim sum As Double = a(0) For i As Integer = 1 To Nx2 Dim prod As Double = a(i) For j As Integer = 0 To (i - 1) prod *= (d - z(j)) Next sum += prod Next Return sum End Function End Module
Module VB0705 'データ点の数 - 1 Private Const N As Integer = 6 Public Sub Main() Dim x(N) As Double Dim y(N) As Double '1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット For i As Integer = 0 To N Dim d1 As Double = i * 1.5 - 4.5 x(i) = d1 y(i) = f(d1) Next '3項方程式の係数の表を作る Dim a(N) As Double Dim b(N) As Double Dim c(N) As Double Dim d(N) As Double For i As Integer = 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(N) As Double Dim s(N) As Double g(1) = b(1) s(1) = d(1) For i As Integer = 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(N) As Double z(0) = 0 z(N) = 0 z(N-1) = s(N-1) / g(N-1) For i As Integer = N - 2 To 1 Step -1 z(i) = (s(i) - c(i) * z(i+1)) / g(i) Next '0.5刻みで 与えられていない値を補間 For i As Integer = 0 To 18 Dim d1 As Double = i * 0.5 - 4.5 Dim d2 As Double = f(d1) Dim d3 As Double = spline(d1, x, y, z) '元の関数と比較 Console.WriteLine(String.Format("{0,5:F2}{4}{1,8:F5}{4}{2,8:F5}{4}{3,8:F5}", d1, d2, d3, d2 - d3, vbTab)) Next End Sub '元の関数 Private Function f(ByVal x As Double) As Double Return x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2) End Function 'Spline (スプライン) 補間 Private Function spline(ByVal d As Double, ByVal x() As Double, ByVal y() As Double, ByVal z() As Double) As Double '補間関数値がどの区間にあるか Dim k As Integer = -1 For i As Integer = 1 To N If d <= x(i) Then k = i - 1 Exit For End If Next If (k < 0) Then k = N Dim d1 As Double = x(k+1) - d Dim d2 As Double = d - x(k) Dim d3 As Double = x(k+1) - x(k) Return (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 End Module
Option Explicit Module VB0801 '重力加速度 Private Const g As Double = -9.8 '空気抵抗係数 Private Const k As Double = -0.01 '時間間隔(秒) Private Const h As Double = 0.01 Public Sub Main() '角度 Const degree As Double = 45 Dim radian As Double = degree * Math.PI / 180.0 '初速 250 km/h -> 秒速に変換 Dim v As Double = 250 * 1000 \ 3600 '水平方向の速度 Dim vx(1) As Double vx(0) = v * Math.Cos(radian) '鉛直方向の速度 Dim vy(1) As Double vy(0) = v * Math.Sin(radian) '経過秒数 Dim t As Double = 0.0 '位置 Dim x As Double = 0.0 Dim y As Double = 0.0 'Euler法 Dim i As Integer = 1 Do While (y >= 0.0) '経過秒数 t = i * h '位置 x += h * vx(0) y += h * vy(0) Console.WriteLine(String.Format("{0,4:F2}{5}{1,8:F5}{5}{2,9:F5}{5}{3,9:F5}{5}{4,8:F5}", t, vx(0), vy(0), x, y, vbTab)) '速度 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 += 1 Loop End Sub '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx As Double, ByVal vy As Double) As Double Return k * Math.Sqrt(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx As Double, ByVal vy As Double) As Double Return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) End Function End Module
Option Explicit Module VB0802 '重力加速度 Private Const g As Double = -9.8 '空気抵抗係数 Private Const k As Double = -0.01 '時間間隔(秒) Private Const h As Double = 0.01 Public Sub Main() '角度 Const degree As Double = 45 Dim radian As Double = degree * Math.PI / 180.0 '初速 250 km/h -> 秒速に変換 Dim v As Double = 250 * 1000 \ 3600 '水平方向の速度 Dim vx(2) As Double vx(0) = v * Math.Cos(radian) '鉛直方向の速度 Dim vy(2) As Double vy(0) = v * Math.Sin(radian) '経過秒数 Dim t As Double = 0.0 '位置 Dim x(2) As Double x(0) = 0.0 Dim y(2) As Double y(0) = 0.0 'Heun法 Dim i As Integer = 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) Console.WriteLine(String.Format("{0,4:F2}{5}{1,8:F5}{5}{2,9:F5}{5}{3,9:F5}{5}{4,8:F5}", t, vx(0), vy(0), x(0), y(0), vbTab)) i += 1 Loop End Sub '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx As Double, ByVal vy As Double) As Double Return k * Math.Sqrt(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx As Double, ByVal vy As Double) As Double Return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) End Function End Module
Option Explicit Module VB0803 '重力加速度 Private Const g As Double = -9.8 '空気抵抗係数 Private Const k As Double = -0.01 '時間間隔(秒) Private Const h As Double = 0.01 Public Sub Main() '角度 Const degree As Double = 45 Dim radian As Double = degree * Math.PI / 180.0 '初速 250 km/h -> 秒速に変換 Dim v As Double = 250 * 1000 \ 3600 '水平方向の速度 Dim vx(1) As Double vx(0) = v * Math.Cos(radian) '鉛直方向の速度 Dim vy(1) As Double vy(0) = v * Math.Sin(radian) '経過秒数 Dim t As Double = 0.0 '位置 Dim x(1) As Double x(0) = 0.0 Dim y(1) As Double y(0) = 0.0 '中点法 Dim i As Integer = 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 As Double = vx(0) + vx(1) / 2.0 Dim wy As Double = 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 Console.WriteLine(String.Format("{0,4:F2}{5}{1,8:F5}{5}{2,9:F5}{5}{3,9:F5}{5}{4,9:F5}", t, vx(0), vy(0), x(0), y(0), vbTab)) i += 1 Loop End Sub '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx As Double, ByVal vy As Double) As Double Return k * Math.Sqrt(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx As Double, ByVal vy As Double) As Double Return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) End Function End Module
Option Explicit Module VB0804 '重力加速度 Private Const g As Double = -9.8 '空気抵抗係数 Private Const k As Double = -0.01 '時間間隔(秒) Private Const h As Double = 0.01 Public Sub Main() '角度 Const degree As Double = 45 Dim radian As Double = degree * Math.PI / 180.0 '初速 250 km/h -> 秒速に変換 Dim v As Double = 250 * 1000 \ 3600 '水平方向の速度 Dim vx(5) As Double vx(0) = v * Math.Cos(radian) '鉛直方向の速度 Dim vy(5) As Double vy(0) = v * Math.Sin(radian) '経過秒数 Dim t As Double = 0.0 '位置 Dim x(5) As Double x(0) = 0.0 Dim y(5) As Double y(0) = 0.0 'Runge-Kutta法 Dim i As Integer = 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 As Double = vx(0) + vx(1) / 2.0 Dim wy As Double = 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(1) + x(2) * 2 + x(3) * 2 + x(4)) / 6 y(0) += ( y(1) + y(2) * 2 + y(3) * 2 + y(4)) / 6 vx(0) += (vx(1) + vx(2) * 2 + vx(3) * 2 + vx(4)) / 6 vy(0) += (vy(1) + vy(2) * 2 + vy(3) * 2 + vy(4)) / 6 Console.WriteLine(String.Format("{0,4:F2}{5}{1,8:F5}{5}{2,9:F5}{5}{3,9:F5}{5}{4,9:F5}", t, vx(0), vy(0), x(0), y(0), vbTab)) i += 1 Loop End Sub '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx As Double, ByVal vy As Double) As Double Return k * Math.Sqrt(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx As Double, ByVal vy As Double) As Double Return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) End Function End Module
Option Explicit Module VB0805 '重力加速度 Private Const g As Double = -9.8 '空気抵抗係数 Private Const k As Double = -0.01 '時間間隔(秒) Private Const h As Double = 0.01 Public Sub Main() '角度 Const degree As Double = 45 Dim radian As Double = degree * Math.PI / 180.0 '初速 250 km/h -> 秒速に変換 Dim v As Double = 250 * 1000 \ 3600 '水平方向の速度 Dim vx(5) As Double vx(0) = v * Math.Cos(radian) '鉛直方向の速度 Dim vy(5) As Double vy(0) = v * Math.Sin(radian) '経過秒数 Dim t As Double = 0.0 '位置 Dim x(5) As Double x(0) = 0.0 Dim y(5) As Double y(0) = 0.0 'Runge-Kutta-Gill法 Dim i As Integer = 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 As Double = vx(0) + vx(1) / 2.0 Dim wy As Double = 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) * ((Math.Sqrt(2.0) - 1) / 2) + vx(2) * (1 - 1 / Math.Sqrt(2.0)) wy = vy(0) + vy(1) * ((Math.Sqrt(2.0) - 1) / 2) + vy(2) * (1 - 1 / Math.Sqrt(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) / Math.Sqrt(2.0) + vx(3) * (1 + 1 / Math.Sqrt(2.0)) wy = vy(0) - vy(2) / Math.Sqrt(2.0) + vy(3) * (1 + 1 / Math.Sqrt(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(1) + x(2) * (2 - Math.Sqrt(2.0)) + x(3) * (2 + Math.Sqrt(2.0)) + x(4)) / 6 y(0) += ( y(1) + y(2) * (2 - Math.Sqrt(2.0)) + y(3) * (2 + Math.Sqrt(2.0)) + y(4)) / 6 vx(0) += (vx(1) + vx(2) * (2 - Math.Sqrt(2.0)) + vx(3) * (2 + Math.Sqrt(2.0)) + vx(4)) / 6 vy(0) += (vy(1) + vy(2) * (2 - Math.Sqrt(2.0)) + vy(3) * (2 + Math.Sqrt(2.0)) + vy(4)) / 6 Console.WriteLine(String.Format("{0,4:F2}{5}{1,8:F5}{5}{2,9:F5}{5}{3,9:F5}{5}{4,9:F5}", t, vx(0), vy(0), x(0), y(0), vbTab)) i += 1 Loop End Sub '空気抵抗による水平方向の減速分 Private Function fx(ByVal vx As Double, ByVal vy As Double) As Double Return k * Math.Sqrt(vx * vx + vy * vy) * vx End Function '重力と空気抵抗による鉛直方向の減速分 Private Function fy(ByVal vx As Double, ByVal vy As Double) As Double Return g + (k * Math.Sqrt(vx * vx + vy * vy) * vy) End Function End Module
Option Explicit Module VB0901 Public Sub Main() Dim a As Double = 1.0 Dim b As Double = 2.0 Console.WriteLine(String.Format("{0,12:F10}", bisection(a, b))) End Sub Private Function bisection(ByVal a As Double, ByVal b As Double) Dim c As Double Do While(True) '区間 (a, b) の中点 c = (a + b) / 2 c = (a + b) / 2 Console.WriteLine(String.Format("{0,12:F10}{2}{1,13:F10}", c, c - Math.Sqrt(2), vbTab)) Dim fc As Double = f(c) If Math.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 Return c End Function Private Function f(ByVal x As Double) As Double Return x * x - 2.0 End Function End Module
Option Explicit Module VB0902 Public Sub Main() Dim a As Double = 1.0 Dim b As Double = 2.0 Console.WriteLine(String.Format("{0,12:F10}", falseposition(a, b))) End Sub Private Function falseposition(ByVal a As Double, ByVal b As Double) Dim c As Double Do While(True) '点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c = (a * f(b) - b * f(a)) / (f(b) - f(a)) Console.WriteLine(String.Format("{0,12:F10}{2}{1,13:F10}", c, c - Math.Sqrt(2), vbTab)) Dim fc As Double = f(c) If Math.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 Return c End Function Private Function f(ByVal x As Double) As Double Return x * x - 2.0 End Function End Module
Option Explicit Module VB0903 Public Sub Main() Dim x As Double = 1.0 Console.WriteLine(String.Format("{0,12:F10}", iterative(x))) End Sub Private Function iterative(ByVal x0 As Double) Dim x1 As Double Do While(True) x1 = g(x0) Console.WriteLine(String.Format("{0,12:F10}{2}{1,13:F10}", x1, x1 - Math.Sqrt(2), vbTab)) If Math.Abs(x1 - x0) < 0.0000000001 Then Exit Do x0 = x1 Loop Return x1 End Function Private Function g(ByVal x As Double) As Double Return (x / 2) + (1 / x) End Function End Module
Option Explicit Module VB0904 Public Sub Main() Dim x As Double = 2.0 Console.WriteLine(String.Format("{0,12:F10}", newton(x))) End Sub Private Function newton(ByVal x0 As Double) Dim x1 As Double Do While(True) x1 = x0 - (f0(x0) / f1(x0)) Console.WriteLine(String.Format("{0,12:F10}{2}{1,13:F10}", x1, x1 - Math.Sqrt(2), vbTab)) If Math.Abs(x1 - x0) < 0.0000000001 Then Exit Do x0 = x1 Loop Return x1 End Function Private Function f0(ByVal x As Double) As Double Return x * x - 2 End Function Private Function f1(ByVal x As Double) As Double Return 2 * x End Function End Module
Option Explicit Module VB0905 Public Sub Main() Dim x As Double = 2.0 Console.WriteLine(String.Format("{0,12:F10}", bailey(x))) End Sub Private Function bailey(ByVal x0 As Double) Dim x1 As Double Do While(True) x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0))))) Console.WriteLine(String.Format("{0,12:F10}{2}{1,13:F10}", x1, x1 - Math.Sqrt(2), vbTab)) If Math.Abs(x1 - x0) < 0.0000000001 Then Exit Do x0 = x1 Loop Return x1 End Function Private Function f0(ByVal x As Double) As Double Return x * x - 2 End Function Private Function f1(ByVal x As Double) As Double Return 2 * x End Function Private Function f2(ByVal x As Double) As Double Return 2 End Function End Module
Option Explicit Module VB0906 Public Sub Main() Dim x0 As Double = 1.0 Dim x1 As Double = 2.0 Console.WriteLine(String.Format("{0,12:F10}", secant(x0, x1))) End Sub Private Function secant(ByVal x0 As Double, ByVal x1 As Double) Dim x2 As Double Do While(True) x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) Console.WriteLine(String.Format("{0,12:F10}{2}{1,13:F10}", x2, x2 - Math.Sqrt(2), vbTab)) If Math.Abs(x2 - x1) < 0.0000000001 Then Exit Do x0 = x1 x1 = x2 Loop Return x2 End Function Private Function f(ByVal x As Double) As Double Return x * x - 2 End Function End Module
Option Explicit Module VB1001 Private Const N As Integer = 3 Public Sub Main() Dim a(,) As Double = {{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}} Dim b() As Double = {20,16,8,17} Dim c() As Double = {0,0,0,0} 'ヤコビの反復法 jacobi(a,b,c) Console.WriteLine("X") disp_vector(c) End Sub '1次元配列の表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub 'ヤコビの反復法 Private Sub jacobi(ByVal a(,) As Double, ByVal b() As Double, ByVal x0() As Double) Do While(True) Dim x1() As Double = {0,0,0,0} Dim finish As Boolean = True For i As Integer = 0 To N x1(i) = 0 For j As Integer = 0 To N If j <> i Then x1(i) += a(i,j) * x0(j) End If Next x1(i) = (b(i) - x1(i)) / a(i,i) If (Math.Abs(x1(i) - x0(i)) > 0.0000000001) Then finish = False Next For i As Integer = 0 To N x0(i) = x1(i) Next If finish Then Return disp_vector(x0) Loop End Sub End Module
Option Explicit Module VB1002 Private Const N As Integer = 3 Public Sub Main() Dim a(,) As Double = {{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}} Dim b() As Double = {20,16,8,17} Dim c() As Double = {0,0,0,0} 'ガウス・ザイデル法 gauss(a,b,c) Console.WriteLine("X") disp_vector(c) End Sub '1次元配列の表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub 'ガウス・ザイデル法 Private Sub gauss(ByVal a(,) As Double, ByVal b() As Double, ByVal x0() As Double) Do While(True) Dim finish As Boolean = True For i As Integer = 0 To N Dim x1 As Double = 0 For j As Integer = 0 To N If j <> i Then x1 += a(i,j) * x0(j) End If Next x1 = (b(i) - x1) / a(i,i) If (Math.Abs(x1 - x0(i)) > 0.0000000001) Then finish = False x0(i) = x1 Next If finish Then Return disp_vector(x0) Loop End Sub End Module
Option Explicit Module VB1003 Private Const N As Integer = 3 Public Sub Main() Dim a(,) As Double = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}} Dim b() As Double = {8,17,20,16} 'ピボット選択 pivoting(a,b) Console.WriteLine("pivoting") Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '前進消去 forward_elimination(a,b) Console.WriteLine("forward elimination") Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '後退代入 backward_substitution(a,b) Console.WriteLine("X") disp_vector(b) End Sub '1次元配列の表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '2次元配列の表示 Private Sub disp_matrix(ByVal matrix(,) As Double) Dim i As Integer = 0 For Each col As Double In matrix Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) i += 1 If i > N Then i = 0 Console.WriteLine() End If Next End Sub '前進消去 Private Sub forward_elimination(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To (N - 1) For row As Integer = (pivot + 1) To N Dim s As Double = a(row,pivot) / a(pivot,pivot) For col As Integer = pivot To N a(row,col) -= a(pivot,col) * s Next b(row) -= b(pivot) * s Next Next End Sub '後退代入 Private Sub backward_substitution(ByVal a(,) As Double, ByRef b() As Double) For row As Integer = N To 0 Step -1 For col As Integer = N To (row + 1) Step - 1 b(row) -= a(row,col) * b(col) Next b(row) /= a(row,row) Next End Sub 'ピボット選択 Private Sub pivoting(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N '各列で 一番値が大きい行を 探す Dim max_row As Integer = pivot Dim max_val As Double = 0 For row As Integer = pivot To N If Math.Abs(a(row,pivot)) > max_val Then '一番値が大きい行 max_val = Math.Abs(a(row,pivot)) max_row = row End If Next '一番値が大きい行と入れ替え If max_row <> pivot Then Dim tmp As Double For col As Integer = 0 To N 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 End Module
Option Explicit Module VB1004 Private Const N As Integer = 3 Public Sub Main() Dim a(,) As Double = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}} Dim b() As Double = {8,17,20,16} 'ピボット選択 pivoting(a,b) Console.WriteLine("pivoting") Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '前進消去 forward_elimination(a,b) Console.WriteLine("forward elimination") Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '後退代入 backward_substitution(a,b) Console.WriteLine("X") disp_vector(b) End Sub '1次元配列の表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '2次元配列の表示 Private Sub disp_matrix(ByVal matrix(,) As Double) Dim i As Integer = 0 For Each col As Double In matrix Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) i += 1 If i > N Then i = 0 Console.WriteLine() End If Next End Sub '前進消去 Private Sub forward_elimination(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N For row As Integer = 0 To N If row = pivot Then Continue For Dim s As Double = a(row,pivot) / a(pivot,pivot) For col As Integer = pivot To N a(row,col) -= a(pivot,col) * s Next b(row) -= b(pivot) * s Next Next End Sub '後退代入 Private Sub backward_substitution(ByVal a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N b(pivot) /= a(pivot,pivot) Next End Sub 'ピボット選択 Private Sub pivoting(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N '各列で 一番値が大きい行を 探す Dim max_row As Integer = pivot Dim max_val As Double = 0 For row As Integer = pivot To N If Math.Abs(a(row,pivot)) > max_val Then '一番値が大きい行 max_val = Math.Abs(a(row,pivot)) max_row = row End If Next '一番値が大きい行と入れ替え If max_row <> pivot Then Dim tmp As Double For col As Integer = 0 To N 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 End Module
Option Explicit Module VB1005 Private Const N As Integer = 3 Public Sub Main() Dim a(,) As Double = {{-1,-2,7,-2},{1,-1,-2,6},{9,2,1,1},{2,8,-2,1}} Dim b() As Double = {8,17,20,16} 'ピボット選択 pivoting(a,b) Console.WriteLine("pivoting") Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '前進消去 forward_elimination(a,b) Console.WriteLine("LU") disp_matrix(a) 'Ly=b から y を求める (前進代入) Dim y() As Double = {0,0,0,0} forward_substitution(a,b,y) Console.WriteLine("Y") disp_vector(y) 'Ux=y から x を求める (後退代入) Dim x() As Double = {0,0,0,0} backward_substitution(a,y,x) Console.WriteLine("X") disp_vector(x) End Sub '1次元配列の表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '2次元配列の表示 Private Sub disp_matrix(ByVal matrix(,) As Double) Dim i As Integer = 0 For Each col As Double In matrix Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) i += 1 If i > N Then i = 0 Console.WriteLine() End If Next End Sub '前進消去 Private Sub forward_elimination(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To (N - 1) For row As Integer = (pivot + 1) To N Dim s As Double = a(row,pivot) / a(pivot,pivot) For col As Integer = pivot To N a(row,col) -= a(pivot,col) * s 'これが 上三角行列 Next a(row,pivot) = s 'これが 下三角行列 'b(row) -= b(pivot) * s 'この値は変更しない Next Next End Sub '前進代入 Private Sub forward_substitution(ByVal a(,) As Double, ByVal b() As Double, ByRef y() As Double) For row As Integer = 0 To N For col As Integer = 0 To row b(row) -= a(row,col) * y(col) Next y(row) = b(row) Next End Sub '後退代入 Private Sub backward_substitution(ByVal a(,) As Double, ByVal y() As Double, ByRef x() As Double) For row As Integer = N To 0 Step -1 For col As Integer = N To (row + 1) Step - 1 y(row) -= a(row,col) * x(col) Next x(row) = y(row) / a(row,row) Next End Sub 'ピボット選択 Private Sub pivoting(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N '各列で 一番値が大きい行を 探す Dim max_row As Integer = pivot Dim max_val As Double = 0 For row As Integer = pivot To N If Math.Abs(a(row,pivot)) > max_val Then '一番値が大きい行 max_val = Math.Abs(a(row,pivot)) max_row = row End If Next '一番値が大きい行と入れ替え If max_row <> pivot Then Dim tmp As Double For col As Integer = 0 To N 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 End Module
Option Explicit Module VB1006 Private Const N As Integer = 3 Public Sub Main() Dim a(,) As Double = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}} Dim b() As Double = {34,68,96,125} Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '前進消去 forward_elimination(a,b) Console.WriteLine("LL^T") disp_matrix(a) 'Ly=b から y を求める (前進代入) Dim y() As Double = {0,0,0,0} forward_substitution(a,b,y) Console.WriteLine("Y") disp_vector(y) 'L^Tx=y から x を求める (後退代入) Dim x() As Double = {0,0,0,0} backward_substitution(a,y,x) Console.WriteLine("X") disp_vector(x) End Sub '1次元配列の表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '2次元配列の表示 Private Sub disp_matrix(ByVal matrix(,) As Double) Dim i As Integer = 0 For Each col As Double In matrix Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) i += 1 If i > N Then i = 0 Console.WriteLine() End If Next End Sub '前進消去 Private Sub forward_elimination(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N Dim s As Double = 0 For col As Integer = 0 To (pivot - 1) s += a(pivot,col) * a(pivot,col) Next 'ここで根号の中が負の値になると計算できない! a(pivot,pivot) = Math.Sqrt(a(pivot,pivot) - s) For row As Integer = (pivot + 1) To N s = 0 For col As Integer = 0 To (pivot - 1) 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(,) As Double, ByVal b() As Double, ByRef y() As Double) For row As Integer = 0 To N For col As Integer = 0 To 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(,) As Double, ByVal y() As Double, ByRef x() As Double) For row As Integer = N To 0 Step -1 For col As Integer = N To (row + 1) Step - 1 y(row) -= a(row,col) * x(col) Next x(row) = y(row) / a(row,row) Next End Sub End Module
Option Explicit Module VB1007 Private Const N As Integer = 3 Public Sub Main() Dim a(,) As Double = {{5,2,3,4},{2,10,6,7},{3,6,15,9},{4,7,9,20}} Dim b() As Double = {34,68,96,125} Console.WriteLine("A") disp_matrix(a) Console.WriteLine("B") disp_vector(b) Console.WriteLine '前進消去 forward_elimination(a,b) Console.WriteLine("LDL^T") disp_matrix(a) 'Ly=b から y を求める (前進代入) Dim y() As Double = {0,0,0,0} forward_substitution(a,b,y) Console.WriteLine("Y") disp_vector(y) 'DL^Tx=y から x を求める (後退代入) Dim x() As Double = {0,0,0,0} backward_substitution(a,y,x) Console.WriteLine("X") disp_vector(x) End Sub '1次元配列の表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '2次元配列の表示 Private Sub disp_matrix(ByVal matrix(,) As Double) Dim i As Integer = 0 For Each col As Double In matrix Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) i += 1 If i > N Then i = 0 Console.WriteLine() End If Next End Sub '前進消去 Private Sub forward_elimination(ByRef a(,) As Double, ByRef b() As Double) For pivot As Integer = 0 To N 'pivot < k の場合 Dim s As Double = 0 For col As Integer = 0 To (pivot - 1) s = a(pivot,col) For k As Integer = 0 To (col - 1) 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 As Integer = 0 To (pivot - 1) s -= a(pivot,k) * a(pivot,k) * a(k,k) Next a(pivot,pivot) = s Next End Sub '前進代入 Private Sub forward_substitution(ByVal a(,) As Double, ByVal b() As Double, ByRef y() As Double) For row As Integer = 0 To N For col As Integer = 0 To row b(row) -= a(row,col) * y(col) Next y(row) = b(row) Next End Sub '後退代入 Private Sub backward_substitution(ByVal a(,) As Double, ByVal y() As Double, ByRef x() As Double) For row As Integer = N To 0 Step -1 For col As Integer = N To (row + 1) Step - 1 y(row) -= a(row,col) * a(row,row) * x(col) Next x(row) = y(row) / a(row,row) Next End Sub End Module
Option Explicit Module VB1101 Private Const N As Integer = 3 'ベキ乗法で最大固有値を求める Public Sub Main() Dim a(,) As Double = {{5.0, 4.0, 1.0, 1.0}, _ {4.0, 5.0, 1.0, 1.0}, _ {1.0, 1.0, 4.0, 2.0}, _ {1.0, 1.0, 2.0, 4.0}} Dim x() As Double = {1.0 ,0.0 ,0.0 ,0.0} 'ベキ乗法 Dim lambda As Double = power(a, x) Console.WriteLine() Console.WriteLine("eigenvalue") Console.WriteLine(String.Format("{0,14:F10}", lambda)) Console.WriteLine("eigenvector") disp_vector(x) End Sub 'ベキ乗法 Private Function power(ByVal a(,) As Double, ByRef x0() As Double) As Double Dim lambda As Double = 0.0 '正規化 (ベクトル x0 の長さを1にする) normarize(x0) Dim e0 As Double = 0.0 For i As Integer = 0 To N e0 += x0(i) Next For k As Integer = 1 To 100 '1次元配列を表示 Console.Write(String.Format("{0,3:D}{1}", k, vbTab)) disp_vector(x0) '行列の積 x1 = A × x0 Dim x1() As Double = {0.0, 0.0, 0.0, 0.0} For i As Integer = 0 To N For j As Integer = 0 To N x1(i) += a(i, j) * x0(j) Next Next '内積 Dim p0 As Double = 0.0 Dim p1 As Double = 0.0 For i As Integer = 0 To N p0 += x1(i) * x1(i) p1 += x1(i) * x0(i) Next '固有値 lambda = p0 / p1 '正規化 (ベクトル x1 の長さを1にする) normarize(x1) '収束判定 Dim e1 As Double = 0.0 For i As Integer = 0 To N e1 += x1(i) Next If Math.Abs(e1 - e0) < 0.00000000001 Then Exit For For i As Integer = 0 To N x0(i) = x1(i) Next e0 = e1 Next Return lambda End Function '1次元配列を表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '正規化 (ベクトルの長さを1にする) Private Sub normarize(ByRef x() As Double) Dim s As Double = 0.0 For i As Integer = 0 To N s += x(i) * x(i) Next s = Math.Sqrt(s) For i As Integer = 0 To N x(i) /= s Next End Sub End Module
Option Explicit Module VB1102 Private Const N As Integer = 3 '逆ベキ乗法で最小固有値を求める Public Sub Main() Dim a(,) As Double = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} Dim x() As Double = {1.0 ,0.0 ,0.0 ,0.0} 'LU分解 forward_elimination(a) '逆ベキ乗法 Dim lambda As Double = inverse(a, x) Console.WriteLine() Console.WriteLine("eigenvalue") Console.WriteLine(String.Format("{0,14:F10}", lambda)) Console.WriteLine("eigenvector") disp_vector(x) End Sub '逆ベキ乗法 Private Function inverse(ByVal a(,) As Double, ByRef x0() As Double) As Double Dim lambda As Double = 0.0 '正規化 (ベクトル x0 の長さを1にする) normarize(x0) Dim e0 As Double = 0.0 For i As Integer = 0 To N e0 += x0(i) Next For k As Integer = 1 To 100 '1次元配列を表示 Console.Write(String.Format("{0,3:D}{1}", k, vbTab)) disp_vector(x0) 'Ly = b から y を求める (前進代入) Dim b() As Double = New Double(N){} Dim y() As Double = New Double(N){} For i As Integer = 0 To N b(i) = x0(i) Next forward_substitution(a,y,b) 'Ux = y から x を求める (後退代入) Dim x1() As Double = New Double(N){} backward_substitution(a,x1,y) '内積 Dim p0 As Double = 0.0 Dim p1 As Double = 0.0 For i As Integer = 0 To N p0 += x1(i) * x1(i) p1 += x1(i) * x0(i) Next '固有値 lambda = p1 / p0 '正規化 (ベクトル x1 の長さを1にする) normarize(x1) '収束判定 Dim e1 As Double = 0.0 For i As Integer = 0 To N e1 += x1(i) Next If Math.Abs(e1 - e0) < 0.00000000001 Then Exit For For i As Integer = 0 To N x0(i) = x1(i) Next e0 = e1 Next Return lambda End Function 'LU分解 Private Sub forward_elimination(ByRef a(,) As Double) For pivot As Integer = 0 To (N - 1) For row As Integer = (pivot + 1) To N Dim s As Double = a(row, pivot) / a(pivot, pivot) For col As Integer = pivot To N 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(,) As Double, ByRef y() As Double, ByVal b() As Double) For row As Integer = 0 To N For col As Integer = 0 To 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(,) As Double, ByRef x() As Double, ByVal y() As Double) For row As Integer = N To 0 Step -1 For col As Integer = N To (row + 1) Step - 1 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() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '正規化 (ベクトルの長さを1にする) Private Sub normarize(ByRef x() As Double) Dim s As Double = 0.0 For i As Integer = 0 To N s += x(i) * x(i) Next s = Math.Sqrt(s) For i As Integer = 0 To N x(i) /= s Next End Sub End Module
Option Explicit Module VB1103 Private Const N As Integer = 3 'LR分解で固有値を求める Public Sub Main() Dim a(,) As Double = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} Dim l(,) As Double = New Double(N,N){} Dim u(,) As Double = New Double(N,N){} For k As Integer= 1 To 200 'LU分解 decomp(a, l, u) '行列の積 multiply(u, l, a) '対角要素を表示 Console.Write(String.Format("{0,3:D}{1}", k, vbTab)) disp_eigenvalue(a) '収束判定 Dim e As Double = 0.0 For i As Integer = 1 To N For j As Integer = 0 To i - 1 e += Math.Abs(a(i, j)) Next Next If e < 0.00000000001 Then Exit For Next Console.WriteLine() Console.WriteLine("eigenvalue") disp_eigenvalue(a) End Sub 'LU分解 Private Sub decomp(ByRef a(,) As Double, ByRef l(,) As Double, ByRef u(,) As Double) For i As Integer = 0 To N For j As Integer = 0 To N l(i, j) = 0.0 u(i, j) = 0.0 Next Next l(0, 0) = 1.0 For j As Integer = 0 To N u(0, j) = a(0, j) Next For i As Integer = 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 As Integer = 1 To N l(i, i) = 1.0 Dim t As Double = a(i, i) For k As Integer = 0 To i t -= l(i, k) * u(k, i) Next u(i, i) = t For j As Integer = i + 1 To N u(j, i) = 0.0 l(i, j) = 0.0 t = a(j, i) For k As Integer = 0 To i t -= l(j, k) * u(k, i) Next l(j, i) = t / u(i, i) t = a(i, j) For k As Integer = 0 To i t -= l(i, k) * u(k, j) Next u(i, j) = t Next Next End Sub '行列の積 Private Sub multiply(ByRef a(,) As Double, ByRef b(,) As Double, ByRef c(,) As Double) For i As Integer = 0 To N For j As Integer = 0 To N Dim s As Double = 0.0 For k As Integer = 0 To N s += a(i, k) * b(k, j) Next c(i, j) = s Next Next End Sub '対角要素を表示 Private Sub disp_eigenvalue(ByVal a(,) As Double) For i As Integer = 0 To N Console.Write(String.Format("{0,14:F10}{1}", a(i, i), vbTab)) Next Console.WriteLine() End Sub End Module
Option Explicit Module VB1104 Private Const N As Integer = 3 'QR分解で固有値を求める Public Sub Main() Dim a(,) As Double = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} Dim q(,) As Double = New Double(N,N){} Dim r(,) As Double = New Double(N,N){} For k As Integer= 1 To 100 'QR分解 decomp(a, q, r) '行列の積 multiply(r, q, a) '対角要素を表示 Console.Write(String.Format("{0,3:D}{1}", k, vbTab)) disp_eigenvalue(a) '収束判定 Dim e As Double = 0.0 For i As Integer = 1 To N For j As Integer = 0 To i - 1 e += Math.Abs(a(i, j)) Next Next If e < 0.00000000001 Then Exit For Next Console.WriteLine() Console.WriteLine("eigenvalue") disp_eigenvalue(a) End Sub 'QR分解 Private Sub decomp(ByRef a(,) As Double, ByRef q(,) As Double, ByRef r(,) As Double) Dim x() As Double = New Double(N){} For k As Integer = 0 To N For i As Integer = 0 To N x(i) = a(i, k) Next For j As Integer = 0 To k - 1 Dim t As Double = 0.0 For i As Integer = 0 To N t += a(i, k) * q(i, j) Next r(j, k) = t r(k, j) = 0.0 For i As Integer = 0 To N x(i) -= t * q(i, j) Next Next Dim s As Double = 0.0 For i As Integer = 0 To N s += x(i) * x(i) Next r(k, k) = Math.Sqrt(s) For i As Integer = 0 To N q(i, k) = x(i) / r(k, k) Next Next End Sub '行列の積 Private Sub multiply(ByRef a(,) As Double, ByRef b(,) As Double, ByRef c(,) As Double) For i As Integer = 0 To N For j As Integer = 0 To N Dim s As Double = 0.0 For k As Integer = 0 To N s += a(i, k) * b(k, j) Next c(i, j) = s Next Next End Sub '対角要素を表示 Private Sub disp_eigenvalue(ByVal a(,) As Double) For i As Integer = 0 To N Console.Write(String.Format("{0,14:F10}{1}", a(i, i), vbTab)) Next Console.WriteLine() End Sub End Module
Option Explicit Module VB1105 Private Const N As Integer = 3 'ヤコビ法で固有値を求める Public Sub Main() Dim a(,) As Double = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} Dim v(,) As Double = {{1.0 ,0.0 ,0.0 ,0.0}, {0.0 ,1.0 ,0.0 ,0.0}, {0.0 ,0.0 ,1.0 ,0.0}, {0.0 ,0.0 ,0.0 ,1.0}} 'ヤコビ法 jacobi(a, v) Console.WriteLine() Console.WriteLine("eigenvalue") disp_eigenvalue(a) Console.WriteLine() Console.WriteLine("eigenvector") disp_eigenvector(v) End Sub 'ヤコビ法 Private Sub jacobi(ByRef a(,) As Double, ByRef v(,) As Double) For k As Integer = 1 To 100 '最大値を探す Dim p, q As Integer Dim max_val As Double = 0.0 For i As Integer = 0 To N For j As Integer = i + 1 To N If (max_val < Math.Abs(a(i, j))) Then max_val = Math.Abs(a(i, j)) p = i q = j End If Next Next 'θ を求める Dim t As Double = 0.0 If Math.Abs(a(p, p) - a(q, q)) < 0.00000000001 Then 'a_{pp} = a_{qq} のとき、回転角tをπ/4にする t = Math.PI / 4.0 If (a(p, p) < 0) Then t = -t End If Else 'a_{pp} ≠ a_{qq} のとき t = Math.Atan(2.0 * a(p, q) / (a(p, p) - a(q, q))) / 2.0 End If 'θ を使って 行列 U を作成し、A = U^t × A × U Dim c As Double = Math.Cos(t) Dim s As Double = Math.Sin(t) 'U^t × A Dim t1 As Double = 0.0 Dim t2 As Double = 0.0 For i As Integer = 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 As Integer = 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 '対角要素を表示 Console.Write(String.Format("{0,3:D}{1}", k, vbTab)) disp_eigenvalue(a) '収束判定 If max_val < 0.00000000001 Then Exit For Next End Sub '対角要素を表示 Private Sub disp_eigenvalue(ByVal a(,) As Double) For i As Integer = 0 To N Console.Write(String.Format("{0,14:F10}{1}", a(i, i), vbTab)) Next Console.WriteLine() End Sub '固有ベクトルを表示 Private Sub disp_eigenvector(ByVal matrix(,) As Double) For i As Integer = 0 To N Dim row() As Double = New Double(N){} For j As Integer = 0 To N row(j) = matrix(i, j) Next normarize(row) disp_vector(row) Next End Sub '1次元配列を表示 Private Sub disp_vector(ByVal row() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '正規化 (ベクトルの長さを1にする) Private Sub normarize(ByRef x() As Double) Dim s As Double = 0.0 For i As Integer = 0 To N s += x(i) * x(i) Next s = Math.Sqrt(s) For i As Integer = 0 To N x(i) /= s Next End Sub End Module
Option Explicit Module VB1106 Private Const N As Integer = 3 'ハウスホルダー変換とQR分解で固有値を求める Public Sub Main() Dim a(,) As Double = {{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} Dim d() As Double = New Double(N){} Dim e() As Double = New Double(N){} 'ハウスホルダー変換 tridiagonalize(a, d, e) 'QR分解 decomp(a, d, e) Console.WriteLine() Console.WriteLine("eigenvalue") disp_vector(d) Console.WriteLine() Console.WriteLine("eigenvector") disp_matrix(a) End Sub 'ハウスホルダー変換 Private Sub tridiagonalize(ByRef a(,) As Double, ByRef d() As Double, ByRef e() As Double) Dim w() As Double = New Double(N){} Dim v() As Double = New Double(N){} For k As Integer = 0 To N - 2 d(k) = a(k, k) Dim t As Double = 0.0 For i As Integer = k + 1 To N w(i) = a(i, k) t += w(i) * w(i) Next Dim s As Double = Math.Sqrt(t) If w(k + 1) < 0 Then s = -s End If If Math.Abs(s) < 0.00000000001 Then e(k + 1) = 0.0 Else e(k + 1) = -s w(k + 1) += s s = 1 / Math.Sqrt(w(k + 1) * s) For i As Integer = k + 1 To N w(i) *= s Next For i As Integer = k + 1 To N s = 0.0 For j As Integer = k + 1 To N If j <= i Then s += a(i, j) * w(j) Else s += a(j, i) * w(j) End If Next v(i) = s Next s = 0.0 For i As Integer = k + 1 To N s += w(i) * v(i) Next s /= 2.0 For i As Integer = k + 1 To N v(i) -= s * w(i) Next For i As Integer = k + 1 To N For j As Integer = k + 1 To i a(i, j) -= w(i) * v(j) + w(j) * v(i) Next Next '次の行は固有ベクトルを求めないなら不要 For i As Integer = 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 As Integer = N To 0 Step -1 w = {0.0 ,0.0 ,0.0 ,0.0} If k < N - 1 Then For i As Integer = k + 1 To N w(i) = a(i, k) Next For i As Integer = k + 1 To N Dim s As Double = 0.0 For j As Integer = k + 1 To N s += a(i, j) * w(j) Next v(i) = s Next For i As Integer = k + 1 To N For j As Integer = k + 1 To N a(i, j) -= v(i) * w(j) Next Next End If For i As Integer = 0 To N a(i, k) = 0.0 Next a(k, k) = 1.0 Next End Sub 'QR分解 Private Sub decomp(ByRef a(,) As Double, ByRef d() As Double, ByRef e() As Double) e(0) = 1.0 Dim h As Integer = N Do While (Math.Abs(e(h)) < 0.00000000001) h -= 1 Loop Do While h > 0 e(0) = 0.0 Dim l As Integer = h - 1 Do While (Math.Abs(e(l)) >= 0.00000000001) l -= 1 Loop For j As Integer = 1 To 100 Dim w As Double = (d(h - 1) - d(h)) / 2.0 Dim s As Double = Math.Sqrt(w * w + e(h) * e(h)) If w < 0.0 Then s = -s End If Dim x As Double = d(l) - d(h) + e(h) * e(h) / (w + s) Dim y As Double = e(l + 1) Dim z As Double = 0.0 Dim t As Double = 0.0 Dim u As Double = 0.0 For k As Integer = l To h - 1 If Math.Abs(x) >= Math.Abs(y) Then 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 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 As Integer = 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 Console.Write(String.Format("{0,3:D}{1}", j, vbTab)) Call disp_vector(d) '収束判定 If (Math.Abs(e(h)) < 0.00000000001) Then Exit For Next e(0) = 1.0 Do While (Math.Abs(e(h)) < 0.00000000001) h -= 1 Loop Loop '次の行は固有ベクトルを求めないなら不要 For k As Integer = 0 To N - 1 Dim l As Integer = k For i As Integer = k + 1 To N If d(i) > d(l) Then l = i End If Next Dim t As Double = d(k) d(k) = d(l) d(l) = t For i As Integer = 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() As Double) For Each col As Double In row Console.Write(String.Format("{0,14:F10}{1}", col, vbTab)) Next Console.WriteLine() End Sub '2次元配列を表示 Private Sub disp_matrix(ByVal matrix(,) As Double) For row As Integer = 0 To N For col As Integer = 0 To N Console.Write(String.Format("{0,14:F10}{1}", matrix(row, col), vbTab)) Next Console.WriteLine() Next End Sub End Module