构建Ax=b时的方法

GAMES201 高级物理引擎课程专区
#1

大家好！目前正在做Homework 1，碰到的老问题是构建系统的Ax=b时，A的写法往往很繁复。比如欧拉视角的流场有nx * ny个cell，那么A就是一个(nx * ny, nx * ny)的五对角矩阵，加上各种边界条件以后代码就变得非常僵硬。请问这方面有什么值得参考的优雅的实现吗？谢谢

2 个赞
#2

Good question. The space complicity would be O(N^6) if using this attempt in 3D, which is apperantly impratical in production. But don’t worry, there’re some useful tricks to resolve this problem.

In fact, A is often highly-sparsed. Take tridiagonal matrix as example:

A = \begin{bmatrix} 2 &-1& 0& 0 &0&0&-1 \\ -1& 2& -1& 0 &0&0&0 \\ 0 &-1& 2& -1 &0 &0&0 \\ 0 &0& -1& 2 &-1 &0&0 \\ 0& 0& 0 &-1& 2 &-1&0 \\ 0& 0& 0 &0 &-1 &2&-1 \\ -1 &0& 0& 0 &0 &-1&2 \\ \end{bmatrix}

See? Despite the matrix have N^2 elements, only 3N of them are non-zero.

So instead of storing the huge&sparse A with size N\times N,
we may use a small&dense A_{diag} with size N\times3 to store it’s non-zero elements:

A_{diag} = \begin{bmatrix} -1 & 2 & -1 \\ -1 & 2 & -1 \\ -1 & 2 & -1 \\ -1 & 2 & -1 \\ -1 & 2 & -1 \\ -1 & 2 & -1 \\ -1 & 2 & -1 \\ \end{bmatrix}

And use a function A(i, j) to replace A[i, j]:

def A(i, j):
  ret = 0.0
  if abs(i - j) <= 1:  # is around diagonal?
    ret = A_diag[i, j - i + 1]

To sum up:
Instead of storing the whole matrix in memory, only store the minimal required elements to reproduce the whole matrix.

This idea is called on-fly, which is one of the @yuanming -styled term used in our course.

3 个赞