🍠超定和欠定线性方程组C代码数值算法

关键词

C代码 | 数值 | 算法 | 超定线性方程组 | 欠定线性方程组 | 对偶单纯形法 | 三角分解 | 线性规划 | 范数 | 平面曲线 | 有界变量 | 残差 | 切比雪夫 | 近似 | 不等式 | 数学 | 离散点 | 容差 | 分段线性 | 欧拉常数

超定和欠定线性方程组 C代码数值算法 | Notion亞圖跨際 on Notion

🏈指点迷津 | Brief

要点

1. 超定线性方程组问题：函数1使用对偶单纯形法和基矩阵的三角分解求解 $L_1$范数。函数2使用对偶单纯形法来制定问题的对偶线性规划，跳过某些中间单纯形迭代。
2. 线性单侧 $L_1$ 近似问题：函数从超定线性方程组的上方或下方计算单侧 $L_1$解，使用改进的单纯形法来制定问题的线性规划。
3. 有界变量的 $L_1$ 近似问题：在解向量的元素限制在 -1 和 1 之间的条件下，函数使用改进的单纯形法，并跳过某些中间单纯形迭代、计算超定线性方程组的 $L_1$ 解。
4. 平面曲线的 $L_1$ 多边形近似：代码计算离散点集 {x,y} 的直线多边形。 该近似值使得任何线的 L_1误差范数不超过预先指定的容差。 多边形中的线数事先是未知的，可以是开放的或封闭的。
5. 平面曲线的分段 $L_1$ 近似：代码1 计算离散点集 {x,y} 的线性分段 $L_1$ 近似。 该近似值使得任何段的 $L_1$ 残差（误差）范数都不会超过给定的容差。 代码2 计算由平面曲线 $y = f(x)$ 离散化所得的给定数据点集 {x,y} 的“接近平衡”分段线性 $L_1​$近似值。
6. 线性切比雪夫近似：该程序使用改进的单纯形法来解决问题的线性规划问题，求解切比雪夫 $\mathrm{L}_{\infty}$ 范数。
7. 单侧切比雪夫近似：C代码使用改进的单纯形法来解决问题的线性规划问题，从超定线性方程组的上方或下方计算单侧切比雪夫。

8. 有界变量的切比雪夫近似：在解向量的元素介于 -1 和 1 之间的条件下，C代码计算超定线性方程组的切比雪夫解。

9. 受限切比雪夫近似：C代码使用改进的单纯形法来解决问题的线性规划问题，计算超定线性方程组的受限切比雪夫解，出于数值稳定性的目的，该代码对基础矩阵使用三角分解。

10. 严格切比雪夫近似：C代码使用改进的单纯形法来解决问题的线性规划问题，计算超定线性方程组的“严格”切比雪夫解。

11. 分段切比雪夫近似：C代码计算离散点集 {x,y} 的线性分段切比雪夫$（\mathrm{L}_{\infty} ）$近似。 该近似值使得任何段的切比雪夫残差（误差）范数都不会超过给定的容差。
12. 线性不等式的解：包含从系统下方计算单侧“切比雪夫”解和从系统下方计算单侧$\mathrm{L}_1$解。

C代码数值计算示例

余弦和正弦积分

余弦和正弦积分定义为:

$\begin{aligned} \operatorname{Ci}(x) & =\gamma+\ln x+\int_0^x \frac{\cos t-1}{t} d t \\ \operatorname{Si}(x) & =\int_0^x \frac{\sin t}{t} d t \end{aligned}$

这里 $\gamma \approx 0.5772 \ldots$ 是欧拉常数。我们只需要一种方法来计算 $x>0$ 的函数，因为

$\operatorname{Si}(-x)=-\operatorname{Si}(x), \quad \operatorname{Ci}(-x)=\operatorname{Ci}(x)-i \pi$

我们可以再次通过幂级数和复连分数的明智组合来评估这些函数。 此序列有

$\begin{aligned} \operatorname{Si}(x) & =x-\frac{x^3}{3 \cdot 3 !}+\frac{x^5}{5 \cdot 5 !}-\cdots \\ \operatorname{Ci}(x) & =\gamma+\ln x+\left(-\frac{x^2}{2 \cdot 2 !}+\frac{x^4}{4 \cdot 4 !}-\cdots\right) \end{aligned}$

指数积分 $E_1(ix)$ 的连分数为

$\begin{aligned} E_1(i x) & =-\operatorname{Ci}(x)+i[\operatorname{Si}(x)-\pi / 2] \\ & =e^{-i x}\left(\frac{1}{i x+} \frac{1}{1+} \frac{1}{i x+} \frac{2}{1+} \frac{2}{i x+} \cdots\right) \\ & =e^{-i x}\left(\frac{1}{1+i x-} \frac{1^2}{3+i x-} \frac{2^2}{5+i x-} \cdots\right) \end{aligned}$

连分数的“偶”形式在最后一行给出，对于大约相同的计算量，收敛速度是原来的两倍。在这种情况下，从交替级数到连分数的一个很好的交叉点是 $x=2$。至于菲涅尔积分，对于较大的 $x$，精度可能会受到正弦和余弦例程的精度的限制。

 #include <math.h>
 ​
 #include "complex.h"
 ​
 #define EPS 6.0e-8
 #define EULER 0.57721566
 #define MAXIT 100
 #define PIBY2 1.5707963 π / 2.
 #define FPMIN 1.0e-30
 #define TMIN 2.0
 #define TRUE 1
 #define ONE Complex(1.0, 0.0)
 void cisi(float x, float * ci, float * si) {
   void nrerror(char error_text[]);
   int i, k, odd;
   float a, err, fact, sign, sum, sumc, sums, t, term;
   fcomplex h, b, c, d, del;
   t = fabs(x);
   if (t == 0.0) {
     * si = 0.0;
     * ci = -1.0 / FPMIN;
     return;
   }
   if (t > TMIN) {
     b = Complex(1.0, t);
     c = Complex(1.0 / FPMIN, 0.0);
     d = h = Cdiv(ONE, b);
     for (i = 2; i <= MAXIT; i++) {
       a = -(i - 1) * (i - 1);
       b = Cadd(b, Complex(2.0, 0.0));
       d = Cdiv(ONE, Cadd(RCmul(a, d), b));
       c = Cadd(b, Cdiv(Complex(a, 0.0), c));
       del = Cmul(c, d);
       h = Cmul(h, del);
       if (fabs(del.r - 1.0) + fabs(del.i) < EPS) break;
     }
     if (i > MAXIT) nrerror("cf failed in cisi");
     h = Cmul(Complex(cos(t), -sin(t)), h);
     * ci = -h.r;
     * si = PIBY2 + h.i;
   } else {
     if (t < sqrt(FPMIN)) {
       sumc = 0.0;
       sums = t;
     } else {
       sum = sums = sumc = 0.0;
       sign = fact = 1.0;
       odd = TRUE;
       for (k = 1; k <= MAXIT; k++) {
         fact *= t / k;
         term = fact / k;
         sum += sign * term;
         err = term / fabs(sum);
         if (odd) {
           sign = -sign;
           sums = sum;
           sum = sumc;
         } else {
           sumc = sum;
           sum = sums;
         }
         if (err < EPS) break;
         odd = !odd;
       }
       if (k > MAXIT) nrerror("maxits exceeded in cisi");
     }
     * si = sums;
     * ci = sumc + log(t) + EULER;
   }
   if (x < 0.0) * si = -( * si);
 }

道森积分

道森积分 $F(x)$ 定义为

$F(x)=e^{-x^2} \int_0^x e^{t^2} d t$

该函数还可以通过以下方式与复杂的误差函数相关：

$F(z)=\frac{i \sqrt{\pi}}{2} e^{-z^2}[1-\operatorname{erfc}(-i z)]$

$F(z)$​ 的一个显着近似是

$F(z)=\lim _{h \rightarrow 0} \frac{1}{\sqrt{\pi}} \sum_{n \text { odd }} \frac{e^{-(z-n h)^2}}{n}$

上面等式的不同寻常之处在于，随着 $h$ 变小，其精度呈指数级增长，因此相当适中的 $h$ 值（以及相应的级数相当快的收敛）给出了非常准确的近似值。首先方便的是移动求和索引以将其大约集中在指数项的最大值上。定义 $n_0$ 为最接近 x / h 的偶数，且 $x_0 \equiv n_0 h、x^{\prime} \equiv x-x_0$和 $n^{\prime} \equiv n- n_0$，于是

$F(x) \approx \frac{1}{\sqrt{\pi}} \sum_{\substack{n^{\prime}=-N \\ n^{\prime} \text { odd }}}^N \frac{e^{-\left(x^{\prime}-n^{\prime} h\right)^2}}{n^{\prime}+n_0}$

其中，当 h 足够小且 N 足够大时，近似等式是准确的。如果我们注意到这一点，这个公式的计算可以大大加快

$e^{-\left(x^{\prime}-n^{\prime} h\right)^2}=e^{-x^{\prime 2}} e^{-\left(n^{\prime} h\right)^2}\left(e^{2 x^{\prime} h}\right)^{n^{\prime}}$

第一个因子计算一次，第二个因子是要存储的常量数组，第三个因子可以递归计算，因此只需要计算两个指数。通过将总和分解为 $n^{\prime}$ 的正值和负值，还可以利用系数 $e^{-\left(n^{\prime} h\right)^2}​​$ 的对称性。

在下面的例程中，选择 $h=0.4$ 和 $N=11$。由于求和的对称性以及对 n 奇数值的限制，for 循环的限制为 1 到 6 。此浮点版本结果的精度约为 $2 \times 10^{-7}$。为了保持 x=0 附近的相对精度（其中 $F(x)$ 消失），程序分支到评估 $F(x)$ 的幂级数 ，其中 $|x|<0.2$ ​。

 #include <math.h>
 ​
 #include "nrutil.h"
 ​
 #define NMAX 6
 #define H 0.4
 #define A1(2.0 / 3.0)
 #define A2 0.4
 #define A3(2.0 / 7.0)
 float dawson(float x)
 {
   int i, n0;
   float d1, d2, e1, e2, sum, x2, xp, xx, ans;
   static float c[NMAX + 1];
   static int init = 0;
   
   if (init == 0) {
     init = 1;
     for (i = 1; i <= NMAX; i++) c[i] = exp(-SQR((2.0 * i - 1.0) * H));
   }
   if (fabs(x) < 0.2) {
 ​
     x2 = x * x;
     ans = x * (1.0 - A1 * x2 * (1.0 - A2 * x2 * (1.0 - A3 * x2)));
   } else {
     
     xx = fabs(x);
     n0 = 2 * (int)(0.5 * xx / H + 0.5);
     xp = xx - n0 * H;
     e1 = exp(2.0 * xp * H);
     e2 = e1 * e1;
     d1 = n0 + 1;
     d2 = d1 - 2.0;
     sum = 0.0;
     for (i = 1; i <= NMAX; i++, d1 += 2.0, d2 -= 2.0, e1 *= e2)
       sum += c[i] * (e1 / d1 + 1.0 / (d2 * e1));
     ans = 0.5641895835 * SIGN(exp(-xp * xp), x) * sum;
   }
   return ans;
 }