📄 nr.cs
字号:
anorm = Math.Max(anorm, (Math.Abs(w[i]) + Math.Abs(rv1[i])));
}
for (i = n - 1; i >= 0; i--)
{
if (i < n - 1)
{
if (g != 0.0)
{
for (j = l; j < n; j++)
v[j, i] = (a[i, j] / a[i, l]) / g;
for (j = l; j < n; j++)
{
for (s = 0.0, k = l; k < n; k++) s += a[i, k] * v[k, j];
for (k = l; k < n; k++) v[k, j] += s * v[k, i];
}
}
for (j = l; j < n; j++) v[i, j] = v[j, i] = 0.0;
}
v[i, i] = 1.0;
g = rv1[i];
l = i;
}
for (i = Math.Min(m, n) - 1; i >= 0; i--)
{
l = i + 1;
g = w[i];
for (j = l; j < n; j++) a[i, j] = 0.0;
if (g != 0.0)
{
g = 1.0 / g;
for (j = l; j < n; j++)
{
for (s = 0.0, k = l; k < m; k++) s += a[k, i] * a[k, j];
f = (s / a[i, i]) * g;
for (k = i; k < m; k++) a[k, j] += f * a[k, i];
}
for (j = i; j < m; j++) a[j, i] *= g;
}
else for (j = i; j < m; j++) a[j, i] = 0.0;
++a[i, i];
}
for (k = n - 1; k >= 0; k--)
{
for (its = 0; its < 30; its++)
{
flag = true;
for (l = k; l >= 0; l--)
{
nm = l - 1;
if (Math.Abs(rv1[l]) + anorm == anorm)
{
flag = false;
break;
}
if (Math.Abs(w[nm]) + anorm == anorm) break;
}
if (flag)
{
c = 0.0;
s = 1.0;
for (i = l - 1; i < k + 1; i++)
{
f = s * rv1[i];
rv1[i] = c * rv1[i];
if (Math.Abs(f) + anorm == anorm) break;
g = w[i];
h = SQRT(f, g);
w[i] = h;
h = 1.0 / h;
c = g * h;
s = -f * h;
for (j = 0; j < m; j++)
{
y = a[j, nm];
z = a[j, i];
a[j, nm] = y * c + z * s;
a[j, i] = z * c - y * s;
}
}
}
z = w[k];
if (l == k)
{
if (z < 0.0)
{
w[k] = -z;
for (j = 0; j < n; j++) v[j, k] = -v[j, k];
}
break;
}
if (its == 29) throw new NRException("no convergence in 30 svdcmp iterations");
x = w[l];
nm = k - 1;
y = w[nm];
g = rv1[nm];
h = rv1[k];
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
g = SQRT(f, 1.0);
f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;
c = s = 1.0;
for (j = l; j <= nm; j++)
{
i = j + 1;
g = rv1[i];
y = w[i];
h = s * g;
g = c * g;
z = SQRT(f, h);
rv1[j] = z;
c = f / z;
s = h / z;
f = x * c + g * s;
g = g * c - x * s;
h = y * s;
y *= c;
for (jj = 0; jj < n; jj++)
{
x = v[jj, j];
z = v[jj, i];
v[jj, j] = x * c + z * s;
v[jj, i] = z * c - x * s;
}
z = SQRT(f, h);
w[j] = z;
if (z != 0)
{
z = 1.0 / z;
c = f * z;
s = h * z;
}
f = c * g + s * y;
x = c * y - s * g;
for (jj = 0; jj < m; jj++)
{
y = a[jj, j];
z = a[jj, i];
a[jj, j] = y * c + z * s;
a[jj, i] = z * c - y * s;
}
}
rv1[l] = 0.0;
rv1[k] = f;
w[k] = x;
}
}
W = new NRVec(w);
V = new NRMatrix(v);
}
public static void SVbksb(NRMatrix u, NRVec w, NRMatrix v, NRVec b, out NRVec x)
//u,w,v为A的奇异值分解,b为列向量，x为Ax=b的解向量。
{
int jj, j, i;
double s;
double wmax = 0.0;
//找出最大奇异值
for (int k = 0; k < w.Rows; k++)
if (w[k] > wmax) wmax = w[k];
// 设置奇异值域值
double wmin = wmax * (1.0e-16);
// 零化"small" 奇异值
for (int k = 0; k < w.Rows; k++)
if (w[k] < wmin) w[k] = 0.0;
int m = u.Rows;
int n = u.Cols;
NRVec tmp = new NRVec(n);
for (j = 0; j < n; j++)
{
s = 0.0;
if (w[j] != 0.0)
{
for (i = 0; i < m; i++) s += u[i, j] * b[i];
s /= w[j];
}
tmp[j] = s;
}
x = new NRVec(b.Rows);
for (j = 0; j < n; j++)
{
s = 0.0;
for (jj = 0; jj < n; jj++) s += v[j, jj] * tmp[jj];
x[j] = s;
}
}
static double update(NRMatrix u, NRMatrix Flag, double h, int i, int j) //FastSweeping的更新程序
{
double uxmin = 0.0, uymin = 0.0, result = 0.0;
if (i > 0 && i < u.Rows - 1 && j > 0 && j < u.Cols - 1)
{
uxmin = Math.Min(u[i - 1, j], u[i + 1, j]);
uymin = Math.Min(u[i, j - 1], u[i, j + 1]);
}
if (i == 0 && j > 0 && j < u.Cols - 1)
{
uxmin = u[i + 1, j];
uymin = Math.Min(u[i, j - 1], u[i, j + 1]);
}
if (i == 0 && j == 0)
{
uxmin = u[i + 1, j];
uymin = u[i, j + 1];
}
if (i == 0 && j == u.Cols - 1)
{
uxmin = u[i + 1, j];
uymin = u[i, j - 1];
}
if (i == u.Rows - 1 && j < u.Cols - 1 && j > 0)
{
uxmin = u[i - 1, j];
uymin = Math.Min(u[i, j - 1], u[i, j + 1]);
}
if (i == u.Rows - 1 && j == 0)
{
uxmin = u[i - 1, j];
uymin = u[i, j + 1];
}
if (i == u.Rows - 1 && j == u.Cols - 1)
{
uxmin = u[i - 1, j];
uymin = u[i, j - 1];
}
if (i > 0 && i < u.Rows - 1 && j == 0)
{
uxmin = Math.Min(u[i - 1, j], u[i + 1, j]);
uymin = u[i, j + 1];
}
if (i > 0 && i < u.Rows - 1 && j == u.Cols - 1)
{
uxmin = Math.Min(u[i - 1, j], u[i + 1, j]);
uymin = u[i, j - 1];
}
if (Flag[i, j] != 1)
if (Math.Abs(uxmin - uymin) >= h)
result = Math.Min(uxmin, uymin) + h;
else
result = 0.5 * (uxmin + uymin + Math.Sqrt(2 * h * h - (uxmin - uymin) * (uxmin - uymin)));
return result;
}
public static int FastSweeping2D(NRMatrix u, Point[] BP, double h) //距离函数生成器，
//u为包围所有边界点BP的方形计算区域，h为网格步长
//返回值为扫描迭代次数step
{
NRMatrix Flag = new NRMatrix(u.Rows, u.Cols);
//标记边界点
for (int k = 0; k < BP.Length; k++)
Flag[(int)BP[k].X, (int)BP[k].Y] = 1;
//初始化边界点的距离值为零，其余点的距离值为无穷大。
for (int i = 0; i < u.Rows; i++)
for (int j = 0; j < u.Cols; j++)
if (Flag[i, j] == 1) u[i, j] = 0.0;
else
u[i, j] = 1.0E+20;
double maxerror = 1.0, temp = 0.0;
int step = 0;
while (maxerror > 1.0E-08)
{
//扫描开始
maxerror = 0.0;
for (int i = 0; i < u.Rows; i++)
for (int j = 0; j < u.Cols; j++)
//update
{
temp = update(u, Flag, h, i, j);
if (Math.Abs(u[i, j] - temp) > maxerror) maxerror = Math.Abs(u[i, j] - temp);
u[i, j] = temp;
}
for (int i = u.Rows - 1; i != 0; i--)
for (int j = 0; j < u.Cols; j++)
//update
{
temp = update(u, Flag, h, i, j);
if (Math.Abs(u[i, j] - temp) > maxerror) maxerror = Math.Abs(u[i, j] - temp);
u[i, j] = temp;
}
for (int i = u.Rows - 1; i != 0; i--)
for (int j = u.Cols - 1; j != 0; j--)
//update
{
temp = update(u, Flag, h, i, j);
if (Math.Abs(u[i, j] - temp) > maxerror) maxerror = Math.Abs(u[i, j] - temp);
u[i, j] = temp;
}
for (int i = 0; i < u.Rows; i++)
for (int j = u.Cols - 1; j != 0; j--)
//update
{
temp =update(u, Flag, h, i, j);
if (Math.Abs(u[i, j] - temp) > maxerror) maxerror = Math.Abs(u[i, j] - temp);
u[i, j] = temp;
}
step++;
if (step > 10000) break;//一万次扫描，不能满足精度要求，扫描结束。
}
return step;
}
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -