Блог

Я написал код, который решает систему линейных уравнений сервал раз подряд во время моделирования методом Монте-Карло.

Для каждого запуска входные данные слегка изменяются, и решение должно вычисляться каждый раз заново. Целью этого является получение функции распределения вероятностей результатов (решений линейной системы).

Итак, мой вопрос таков:
Есть ли способ решить систему линейных уравнений только один раз и сохранить общее решение, чтобы для каждого запуска по методу Монте-Карло решения можно было вычислять напрямую?

Это позволило бы очень сэкономить время, поскольку для правильного моделирования мне нужно не менее 20 тыс. запусков, и даже для небольших систем с тремя неизвестными это занимает много времени. Мой код решает это линейное уравнение каждый раз по-новому, поскольку в его исходной версии количество переменных и, следовательно, количество входных величин должно быть замкнутым, и, следовательно, общие решения неизвестны.

Вот мой алгоритм исключения по Гауссу.

 Function gaussian_elimination(w As Variant, mm As Variant, R As Variant, rb As Variant, n_iso As Integer) As Variant()

'initializing running indexes
    Dim i As Integer
    Dim j As Integer
    Dim h As Integer
    Dim n As Integer

    n = n_iso

'runing variables for Gauss elimination
    Dim ip As Integer
    Dim q As Integer
    Dim p As Integer
    Dim z As Double
    Dim temp1(1, 1) As Variant
    Dim temp2(1, 1) As Variant
    Dim sum As Variant


'initializing b vector
    Dim b() As Variant
    ReDim b(1 To n - 1, 0 To 1)

'initializing k vector
    Dim k() As Variant
    ReDim k(1 To n - 1, 0 To 1)

'initializing A matrix
    Dim a() As Variant
    ReDim a(1 To n - 1, 1 To n - 1)

'initializing X matrix
    Dim x() As Variant
    ReDim x(1 To n - 1, 1 To n)


' calculating b vector

    For i = 1 To (n - 1) Step 1
        b(i, 0) = mm(1, 0) / (w(i   n - 1, 0) * (rb(i, 0) - R(i * n, 0))) - mm(1, 0) / (w(i, 0) * (R(i, 0) - rb(i, 0)))
    Next i

'calculating A matrix
    For i = 1 To (n - 1) Step 1
        For j = 1 To (n - 1) Step 1
            a(i, j) = mm(j   1, 0) * ((R(j   i * (n - 1), 0) / (w(i   n - 1, 0) * (rb(i, 0) - R(i * n, 0)))) - (R(j, 0) / (w(i, 0) * (R(i, 0) - rb(i, 0)))))
        Next j
    Next i



    'using on board solving routine

    Dim A_Inv As Variant
    Dim k_vec As Variant
    Dim b_dummy As Variant





'filling X matrix
    For i = 1 To n - 1 Step 1
        For j = 1 To n Step 1
                If j = (n) Then
                    x(i, j) = b(i, 0)
                Else: x(i, j) = a(i, j)
                End If
        Next j
    Next i

   'Gaussian elimination
    For i = 1 To (n - 2) Step 1
        For j = i   1 To (n - 1) Step 1
            If (Abs(x(j, i)) > Abs(x(i, i))) Then
                For h = 1 To n
                    temp1(1, 1) = x(i, h)
                    temp2(1, 1) = x(j, h)
                    x(i, h) = temp2(1, 1)
                    x(j, h) = temp1(1, 1)
                Next h
            End If
        Next

        For p = i   1 To n - 1
            z = x(p, i) / x(i, i)
            For q = i   1 To n
                x(p, q) = x(p, q) - z * x(i, q)
            Next q
            x(p, i) = 0
        Next p
    Next i






'calculatiing k factors backwards
    If Abs(x(UBound(x, 1), UBound(x, 2) - 1)) <= 0 Then
        MsgBox "Equation system can not be solved! Solving for k factors faild", vbExclamation, "Warning!"
        Exit Function
    End If

    k((UBound(x, 1)), 0) = x((UBound(x, 1)), UBound(x, 2)) / x((UBound(x, 1)), (UBound(x, 2) - 1))

    For i = ((UBound(x, 1) - 1)) To (LBound(x, 1)) Step -1
        sum = x(i, UBound(x, 2))
        For j = i   1 To (UBound(x, 2) - 1) Step 1
            sum = sum - x(i, j) * k(j, 0)
        Next j
        k(i, 0) = sum / x(i, i)
    Next i

    For i = 1 To n - 1
        k(i, 0) = (-1) * k(i, 0)
    Next i

    gaussian_elimination = k



End Function