Dự kiến ​​Gauss-Seidel cho LCP

Question

Tôi đang tìm triển khai C# của thuật toán Gauss-Seidel được chiếu để giải quyết linear complementarity problem. Cho đến nay tôi đã tìm thấy một trong những viết bằng C + + trong thư viện Bullet, nhưng tiếc là nó được tối ưu hóa cao (vì vậy nó sẽ rất khó để dịch nó thành C#).

Trong câu hỏi similar được đề xuất để xem xét numerical libraries for .NET. Tất cả chúng chỉ chứa các thuật toán để giải quyết systems of linear equations.

Chỉnh sửa: Ngay cả khi tôi đã tìm thấy, nó có vẻ không hoàn chỉnh, do đó câu hỏi vẫn mở.

Answer 1

Bạn đã thực hiện Gauss Seidel mà không cần chiếu. Đối với Gauss Seidel dự kiến, bạn cần phải dự án (hoặc kẹp) giải pháp trong giới hạn dưới và trên:

public static double[] Solve (double[,] matrix, double[] right, 
           double relaxation, int iterations, 
           double[] lo, double[] hi) 
{ 
    // Validation omitted 
    var x = right; 
    double delta; 

    // Gauss-Seidel with Successive OverRelaxation Solver 
    for (int k = 0; k < iterations; ++k) { 
     for (int i = 0; i < right.Length; ++i) { 
      delta = 0.0f; 

      for (int j = 0; j < i; ++j) 
       delta += matrix [i, j] * x [j]; 
      for (int j = i + 1; j < right.Length; ++j) 
       delta += matrix [i, j] * x [j]; 

      delta = (right [i] - delta)/matrix [i, i]; 
      x [i] += relaxation * (delta - x [i]); 
    // Project the solution within the lower and higher limits 
      if (x[i]<lo[i]) 
       x[i]=lo[i]; 
      if (x[i]>hi[i]) 
       x[i]=hi[i]; 
     } 
    } 
    return x; 
} 

Đây là một sửa đổi nhỏ. Dưới đây là một ý chính cho thấy cách trích xuất ma trận A và vector b từ Thư viện Vật lý Bullet và giải quyết nó bằng cách sử dụng dự kiến ​​Gauss Seidel: https://gist.github.com/erwincoumans/6666160

Answer 2

Sau một tuần tìm kiếm tôi cuối cùng đã tìm thấy this ấn phẩm (bằng tiếng Nga, dựa trên công việc của Kenny Erleben). Một thuật toán Gauss-Seidel dự kiến ​​được mô tả ở đó và sau đó mở rộng với SOR và điều kiện chấm dứt. Tất cả những điều đó với các ví dụ trong C++, mà tôi đã sử dụng cho việc triển khai C# cơ bản này:

public static double[] Solve (double[,] matrix, double[] right, 
           double relaxation, int iterations) 
{ 
    // Validation omitted 
    var x = right; 
    double delta; 

    // Gauss-Seidel with Successive OverRelaxation Solver 
    for (int k = 0; k < iterations; ++k) { 
     for (int i = 0; i < right.Length; ++i) { 
      delta = 0.0f; 

      for (int j = 0; j < i; ++j) 
       delta += matrix [i, j] * x [j]; 
      for (int j = i + 1; j < right.Length; ++j) 
       delta += matrix [i, j] * x [j]; 

      delta = (right [i] - delta)/matrix [i, i]; 
      x [i] += relaxation * (delta - x [i]); 
     } 
    } 

    return x; 
}