Option Explicit Public Const m As Integer = 3 Public Const s As Integer = 3 Public Const n As Integer = 9 Dim eps As Double, epd As Double, q As Double, f As Double Sub ƒ{ƒ^ƒ“1_Click() Dim a(s, s) As Double Dim b(s) As Double Dim c(s) As Double Dim En(n, n) As Double, Em(m, m) As Double Dim jacob(m, m) As Double Dim indx(n) As Integer Dim z(n) As Double Dim mat_temp(n, n) As Double, lumat(n, n) As Double, vec(n) As Double Dim ff(n) As Double Dim d As Double Dim i As Integer, j As Integer, k As Integer, kk As Integer Dim x1 As Double, x2 As Double, x3 As Double Dim t As Double, dt As Double Application.ScreenUpdating = False eps = 0.0099 epd = 0.0000198 q = 0.0000762 f = 1# For i = 1 To n For j = 1 To n If i = j Then En(i, j) = 1# Else En(i, j) = 0# End If Next j Next i For i = 1 To m For j = 1 To m If i = j Then Em(i, j) = 1# Else Em(i, j) = 0# End If Next j Next i a(1, 1) = (88# - 7# * Sqr(6#)) / 360# a(1, 2) = (296# - 169# * Sqr(6#)) / 1800# a(1, 3) = (-2# + 3# * Sqr(6#)) / 225# a(2, 1) = (296# + 169# * Sqr(6#)) / 1800# a(2, 2) = (88# + 7# * Sqr(6#)) / 360# a(2, 3) = (-2# - 3# * Sqr(6#)) / 225# a(3, 1) = (16# - Sqr(6#)) / 36# a(3, 2) = (16# + Sqr(6#)) / 36# a(3, 3) = 1# / 9# b(1) = (16# - Sqr(6#)) / 36# b(2) = (16# + Sqr(6#)) / 36# b(3) = 1# / 9# c(1) = (4# - Sqr(6#)) / 10# c(2) = (4# + Sqr(6#)) / 10# c(3) = 1# t = 0# dt = 0.002 x1 = 1# x2 = 1# x3 = 1# Worksheets("Sheet1").Cells(2, 1) = t Worksheets("Sheet1").Cells(2, 2) = x1 Worksheets("Sheet1").Cells(2, 3) = x2 Worksheets("Sheet1").Cells(2, 4) = x3 For k = 1 To 30000 For i = 1 To n z(i) = 0# Next i jacob(1, 1) = f11(t, x1, x2, x3) jacob(1, 2) = f12(t, x1, x2, x3) jacob(1, 3) = f13(t, x1, x2, x3) jacob(2, 1) = f21(t, x1, x2, x3) jacob(2, 2) = f22(t, x1, x2, x3) jacob(2, 3) = f23(t, x1, x2, x3) jacob(3, 1) = f31(t, x1, x2, x3) jacob(3, 2) = f32(t, x1, x2, x3) jacob(3, 3) = f33(t, x1, x2, x3) Call tensorproduct(a, jacob, mat_temp) For i = 1 To n For j = 1 To n lumat(i, j) = En(i, j) - dt * mat_temp(i, j) Next j Next i Call ludcmp(lumat, indx, d) Call tensorproduct(a, Em, mat_temp) For kk = 1 To 20 ff(1) = f1(t + c(1) * dt, x1 + z(1), x2 + z(2), x3 + z(3)) ff(2) = f2(t + c(1) * dt, x1 + z(1), x2 + z(2), x3 + z(3)) ff(3) = f3(t + c(1) * dt, x1 + z(1), x2 + z(2), x3 + z(3)) ff(4) = f1(t + c(2) * dt, x1 + z(4), x2 + z(5), x3 + z(6)) ff(5) = f2(t + c(2) * dt, x1 + z(4), x2 + z(5), x3 + z(6)) ff(6) = f3(t + c(2) * dt, x1 + z(4), x2 + z(5), x3 + z(6)) ff(7) = f1(t + c(3) * dt, x1 + z(7), x2 + z(8), x3 + z(9)) ff(8) = f2(t + c(3) * dt, x1 + z(7), x2 + z(8), x3 + z(9)) ff(9) = f3(t + c(3) * dt, x1 + z(7), x2 + z(8), x3 + z(9)) For i = 1 To n vec(i) = -z(i) For j = 1 To n vec(i) = vec(i) + dt * mat_temp(i, j) * ff(j) Next j Next i Call lubksb(lumat, indx, vec) For i = 1 To n z(i) = z(i) + vec(i) Next i Next kk x1 = x1 + z(7) x2 = x2 + z(8) x3 = x3 + z(9) t = t + dt Worksheets("Sheet1").Cells(k + 2, 1) = t Worksheets("Sheet1").Cells(k + 2, 2) = x1 Worksheets("Sheet1").Cells(k + 2, 3) = x2 Worksheets("Sheet1").Cells(k + 2, 4) = x3 Next k Application.ScreenUpdating = True End Sub Function f1(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f1 = (q * x2 - x1 * x2 + x1 * (1# - x1)) / eps End Function Function f2(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f2 = (-q * x2 - x1 * x2 + f * x3) / epd End Function Function f3(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f3 = x1 - x3 End Function Function f11(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f11 = (-x2 + 1# - 2# * x1) / eps End Function Function f12(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f12 = (q - x1) / eps End Function Function f13(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f13 = 0# End Function Function f21(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f21 = (-x2) / epd End Function Function f22(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f22 = (-q - x1) / epd End Function Function f23(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f23 = f / epd End Function Function f31(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f31 = 1# End Function Function f32(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f32 = 0# End Function Function f33(t As Double, x1 As Double, x2 As Double, x3 As Double) As Double f33 = -1# End Function Sub ludcmp(a() As Double, indx() As Integer, d As Double) Dim i As Integer, imax As Integer, j As Integer, k As Integer Dim big As Double, dum As Double, sum As Double, temp As Double Dim vv(n) As Double d = 1# For i = 1 To n big = 0# For j = 1 To n temp = Abs(a(i, j)) If temp > big Then big = temp End If Next j vv(i) = 1# / big Next i For j = 1 To n For i = 1 To j - 1 sum = a(i, j) For k = 1 To i - 1 sum = sum - a(i, k) * a(k, j) Next k a(i, j) = sum Next i big = 0# For i = j To n sum = a(i, j) For k = 1 To j - 1 sum = sum - a(i, k) * a(k, j) Next k a(i, j) = sum dum = vv(i) * Abs(sum) If dum >= big Then big = dum imax = i End If Next i If j <> imax Then For k = 1 To n dum = a(imax, k) a(imax, k) = a(j, k) a(j, k) = dum Next k d = -d vv(imax) = vv(j) End If indx(j) = imax If a(j, j) = 0# Then a(j, j) = 1E-20 End If If j <> n Then dum = 1# / a(j, j) For i = j + 1 To n a(i, j) = a(i, j) * dum Next i End If Next j End Sub Sub lubksb(a() As Double, indx() As Integer, b() As Double) Dim i As Integer, ii As Integer, ip As Integer, j As Integer Dim sum As Double For i = 1 To n ip = indx(i) sum = b(ip) b(ip) = b(i) If ii <> 0 Then For j = ii To i - 1 sum = sum - a(i, j) * b(j) Next j Else If sum <> 0# Then ii = i End If End If b(i) = sum Next i For i = n To 1 Step -1 sum = b(i) For j = i + 1 To n sum = sum - a(i, j) * b(j) Next j b(i) = sum / a(i, i) Next i End Sub Sub invmat(a() As Double) Dim y(n, n) As Double Dim indx(n) As Integer Dim col(n) As Double Dim d As Double Dim i As Integer, j As Integer Call ludcmp(a, indx, d) For j = 1 To n For i = 1 To n col(i) = 0# Next i col(j) = 1# Call lubksb(a, indx, col) For i = 1 To n y(i, j) = col(i) Next i Next j For i = 1 To n For j = 1 To n a(i, j) = y(i, j) Next j Next i End Sub Function det(a() As Double) As Double Dim j As Integer Dim indx(n) As Integer Dim d As Double Call ludcmp(a, indx, d) For j = 1 To n d = d * a(j, j) Next j det = d End Function Sub linsolve(a() As Double, b() As Double) Dim indx(n) As Integer Dim d As Double Call ludcmp(a, indx, d) Call lubksb(a, indx, b) End Sub Sub tensorproduct(a() As Double, b() As Double, c() As Double) Dim i As Integer, j As Integer Dim i1 As Integer, j1 As Integer Dim i2 As Integer, j2 As Integer For i1 = 1 To s For j1 = 1 To s For i2 = 1 To m For j2 = 1 To m i = (i1 - 1) * m + i2 j = (j1 - 1) * m + j2 c(i, j) = a(i1, j1) * b(i2, j2) Next j2 Next i2 Next j1 Next i1 End Sub