数值分析C程序
// ----------------------------------------------------------------------------
// 多项式求值的秦九韶算法
// 功能： 求 n-1 次多项式在 x 处的函数值
// 参数： a － 存放多项式系数的数组
// n － 多项式系数数组的大小
// x － 函数的自变量的指定值
// 返回值：n-1 次多项式在 x 处的函数值 y = Pn-1(x)
// ----------------------------------------------------------------------------
double Polynomial(double a[], int n, double x)
{
int i;
double t = a[n-1];
for (i=n-2; i>=0; i--) t = t * x + a[i];
return t;
}
///////////////////////////////////////////////////////////////////////////////
// ----------------------------------------------------------------------------
// 带有列主元的高斯消元法
// 功能： 求解线性方程组 Ax = b
// 参数： A － 指向n*n系数矩阵的指针
// b － 常数向量的指针
// n － 方程组的维数
// 返回值：0 － 如果成功。线性方程组的解保存在 b 中
// 1 － 求解失败
// ----------------------------------------------------------------------------
int LEquations_Gauss(double A[], double b[], int n)
{
int i, j, k;
double s, t;
// 高斯列主元消元,使 Ax=b -> Ux=b'
for(k=0; k
// 选取列主元素
j = k; t = fabs(A[k*n+k]);
for(i=k+1; i
if((s = fabs(A[i*n+k])) > t)
{
t = s;
j = i;
}
}
if(t < 1.0e-30) return 1; // 方程为病态,或无穷多解
if(j != k) // 交换方程的第 j 行和 k 行
{
for(i=k; i
t = A[j*n+i];
A[j*n+i] = A[k*n+i];
A[k*n+i] = t;
}
t = b[j]; b[j] = b[k]; b[k] = t;
}
t = 1.0 / A[k*n+k]; // 先将对角线元素变为 1.0
for(i=k+1; i
for(i=k+1; i
t = A[i*n+k];
for(j=k+1; j
}
}
// 回代,求解 Ux=b'
for(i=n-2; i>=0; i--)
for(j=i+1; j
return 0;
}
// ----------------------------------------------------------------------------
// 带有列主元的约当消元法
// 功能： 求解线性方程组 Ax = b
// 参数： A － 指向n*n系数矩阵的指针
// b － 常数向量的指针
// n － 方程组的维数
// 返回值：0 － 如果成功。线性方程组的解保存在 b 中
// 1 － 求解失败
// ----------------------------------------------------------------------------
int LEquations_Jordan(double A[], double b[], int n)
{
int i, j, k;
double s, t;
// 消元主循环,使 A -> I, 从而 b -> b'为所求的解
for(k=0; k
// 选取列主元素
j = k; t = fabs(A[k*n+k]);
for(i=k+1; i
if((s = fabs(A[i*n+k])) > t)
{
t = s;
j = i;
}
}
if(t < 1.0e-30) return 1; // 方程为病态,或无穷多解
if(j != k) // 交换方程的第 j 行和 k 行
{
for(i=k; i
t = A[j*n+i];
A[j*n+i] = A[k*n+i];
A[k*n+i] = t;
}
t = b[j]; b[j] = b[k]; b[k] = t;
}
t = 1.0 / A[k*n+k]; // 先将对角线元素变为 1.0
for(i=k+1; i
for(i=0; i
if(i != k && A[i*n+k] != 0.0)
{
t = A[i*n+k];
for(j=k+1; j
}
}
}
return 0;
}
// ----------------------------------------------------------------------------
// 求矩阵 a 的逆矩阵 b
// 功能： 求矩阵 a 的逆矩阵 b = a-1
// 参数： a － 指向n*n矩阵的指针
// b － 指向n*n矩阵的指针,用于返回a的逆矩阵
// n － 维数
// 返回值：0 － 如果成功。解保存在 b 中
// 1 － 求解失败
// ----------------------------------------------------------------------------
int Matrix_Invert(double a[], double b[], int n)
{
int i, j, k, m;
double d, t;
for(i=0; i
else b[i*n+j] = 0.0;
for(k=0; k
j = k; t = fabs(a[k*n+k]);
for(i=k+1; i
if((d = fabs(a[i*n+k])) > t)
{
t = d;
j = i;
}
}
if(t < 1.0e-30) return 1;
if(j != k)
{
for(i=k; i
t = a[j*n+i];
a[j*n+i] = a[k*n+i];
a[k*n+i] = t;
}
for(i=0; i
t = b[j*n+i];
b[j*n+i] = b[k*n+i];
b[k*n+i] = t;
}
}
// 为了使a的第k行的第k元素为1,a和b的第k行都除以a(k,k)
t = 1.0 / a[k*n+k];
for(i=k; i
for(i=0; i
if(i != k && a[i*n+k] != 0.0)
{
t = a[i*n+k];
for(j=k+1; j
}
}
return 0;
}
// ----------------------------------------------------------------------------
// Doolittle分解法(直接三角分解)
// 功能： 求解线性方程组 Ax = b
// 参数： A － 指向n*n系数矩阵的指针
// b － 常数向量的指针
// n － 方程组的维数
// 返回值：0 － 如果成功。线性方程组的解保存在 b 中
// 1 － 不能直接三角分解
// ----------------------------------------------------------------------------
int LEquations_LU(double A[], double b[], int n)
{
int i, j, k;
double s, t;
// 进行直接三角分解,L和U分别保存在A的下三角和上三角
for(k=0; k
for(i=k+1; i
第k+1行
{
t = 0.0; // 计算下三角矩阵L的第k列的第i行元素
for(j=0; j
t = 0.0; // 计算上三角矩阵U的第k+1行的第i+1列元素
for(j=0; j<=k; j++) t += A[(k+1)*n+j] * A[j*n+i];
A[(k+1)*n+i] -= t;
}
}
// 求解 Ly = b 并将所求的结果(向量y)保存在b中
for(i=1; i
t = 0.0; k = i * n;
for(j=0; j
}
// 求解 Ux = y 并将所求的结果(向量x)保存在b中
for(i=n-1; i>=0; i--)
{
t = 0.0; k = i * n + i + 1;
for(j=i+1; j
}
return 0;
}
// ----------------------------------------------------------------------------
// 带列主元素的三角分解法
// 功能： 求解线性方程组 Ax = b
// 参数： A － 指向n*n系数矩阵的指针,求解成功时保存PA的三角分解LU
// b － 常数向量的指针,求解成功时保存线性方程组的解
// n － 方程组的维数
// 返回值：0 － 如果成功。线性方程组的解保存在 x 中
// 1 － 不能直接三角分解
// ----------------------------------------------------------------------------
int LEquations_GLU(double A[], double b[], int n)
{
int i, j, k;
double s, t;
// 高斯列主元消元,使 PA -> LU, b -> Pb
for(k=0; k
// 选取列主元素
j = k; t = fabs(A[k*n+k]);
for(i=k+1; i
if((s = fabs(A[i*n+k])) > t)
{
t = s;
j = i;
}
}
if(t < 1.0e-30) return 1; // 病态
if(j != k) // 交换方程的第 j 行和 k 行
{
for(i=0; i
t = A[j*n+i];
A[j*n+i] = A[k*n+i];
A[k*n+i] = t;
}
t = b[j]; b[j] = b[k]; b[k] = t;
}
s = 1.0 / A[k*n+k];
for(i=k+1; i
A[i*n+k] *= s; t = A[i*n+k];
for(j=k+1; j
}
}
// 求解 Ly = b 并将所求的结果(向量y)保存在b中
for(i=1; i
t = 0.0; k = i * n;
for(j=0; j
}
// 求解 Ux = y 并将所求的结果(向量x)保存在b中
for(i=n-1; i>=0; i--)
{
t = 0.0; k = i * n + i + 1;
for(j=i+1; j
}
return 0;
}
// ----------------------------------------------------------------------------
// 平方根法(LLT法)
// 功能： 求解对称正定线性方程组 Ax = b
// 参数： A － 指向n*n系数矩阵的指针
// b － 常数向量的指针
// n － 方程组的维数
// 返回值：无。
// 线性方程组的解保存在 b 中
// ----------------------------------------------------------------------------
int LEquations_LT(double *A, doubl
e *b, int n)
{
int i, j, k;
double t;
for(i=0; i
// 先求对角线元素
t = 0.0;
for(k=0; k
if (A[i*n+i] < 1.0e-30) return 1;
// 再求对角线以下的 i 列各元素
for(j=i+1; j
t = 0.0;
for(k=0; k
A[j*n+i] = A[i*n+j] = (A[i*n+j] - t) / A[i*n+i];
}
}
for(i=0; i
t = 0.0;
for(j=0; j
}
for(i=n-1; i>=0; i--) // 求解方程组LTx=y
{
t = 0.0;
for(j=i+1; j
}
return 0;
}
// ----------------------------------------------------------------------------
// 改进的平方根法(LDLT法)
// 功能： 求解对称正定线性方程组 Ax = b
// 参数： A － 指向n*n系数矩阵的指针
// b － 常数向量的指针
// n － 方程组的维数
// 返回值：无。
// 线性方程组的解保存在 b 中
// ----------------------------------------------------------------------------
int LEquations_LDT(double *A, double *b, int n)
{
int i, j, k;
double s, t;
for(i=0; i
// 先求对角线元素
t = 0.0;
for(k=0; k
s = fabs(A[i*i+i]);
if (s < 1.0e-30) return 1;
s = 1.0 / A[i*n+i];
// 再求对角线以下的 i 列各元素
for(j=i+1; j
t = 0.0;
for(k=0; k
A[j*n+i] = A[i*n+j] = (A[i*n+j] - t) * s;
}
}
for(i=1; i
t = 0.0;
for(j=0; j
}
for(i=n-1; i>=0; i--) // 求解方程组DLTx=y
{
t = 0.0;
for(j=i+1; j
}
}
// ----------------------------------------------------------------------------
// 追赶法
// 功能： 求解三对角线性方程组 Ax = b
// 参数： A － 指向3*n系数矩阵的指针(存放顺序为ai,bi,ci行)
// b － 常数向量的指针
// n － 方程组的维数
// 返回值：无。
// 线性方程组的解保存在 b 中
// ----------------------------------------------------------------------------
int LEquations_Pursue(double *A, double *b, int n)
{
int i, j, k;
double t;
// 计算 LU 分解,和求解 Ly = b
for(i=1, j=n+1, k=n+n; i
if (A[j-1] == 0.0) return 1;
A[i] /= A[j-1];
b[i] -= b[i-1] * A[i];
A[j] -= A[i] * A[k];
}
// 解U*x =y(此时的上三角矩阵的对角线元素已化为全1)
k = n + n - 2;
if (A[k+1] == 0.0) return 1;
b[n-1] /= A[k+1];
for(i=n-2; i>=0; i--, k--)
{
if (A[n+i] == 0.0) return 1;
b[i] =
(b[i] - A[k] * b[i+1]) / A[n+i];
}
}
// 该函数用于计算两个向量差的无穷范数的值
double deVec_Inf(double *x1, double *x2, int n)
{
int i;
double t, p=0.0;
for(i=0; i
t = fabs(x1[i] - x2[i]);
if(t > p) p = t;
}
return p;
}
// ----------------------------------------------------------------------------
// 调整方程 Ax=b 试图使系数矩阵对角占优
// 功能： 用全主元素法调整方程 Ax=b 为对角线最大，并判断调整后的方程是否对角占优
// 参数： A － 指向n阶系数矩阵的指针
// b － 常数向量的指针
// n － 方程组的维数
// 返回值：0 － 如果成功。A已对角占优且b也随之改变
// 1 － 失败，A和b均未还原
// -1 － 没有足够的自由空间
// ----------------------------------------------------------------------------
int SuperDiagonal(double *A, double *b, int n)
{
int i, j, k, u, v;
int *p = (int *)malloc(n*sizeof(int));
double s, t;
if (p == NULL) return -1;
for (i=0; i
t = 0.0; u = v = 0;
for (i=0; i
if (p[i]) continue;
for (j=0; j
if (p[j]) continue;
if (t < (s=fabs(A[i*n+j])))
{
t = s;
u = i; v = j;
}
}
}
// 测试该行能否主元素占优
s = 0.0;
for (i=0; i
if (i != v) s += fabs(A[u*n+i]);
}
if (t <= s) { free(p); return 1; } // 不可能对角占优
// 交换主元素到对角线
if (u != v)
{
for (i=0; i
s = A[u*n+i];
A[u*n+i] = A[v*n+i];
A[v*n+i] = s;
}
s = b[u]; b[u] = b[v]; b[v] = s;
}
// 标记 Auv 从主元素已变成对角元素 Avv
p[v] = 1;
}
free(p);
return 0;
}
// ----------------------------------------------------------------------------
// 雅可比迭代法
// 功能： 求解线性方程组 Ax = b
// 将方程写成迭代公式xn+1 = D-1(D-a)xn + D-1b.
// 参数： A － 指向n阶系数矩阵的指针
// b － 常数向量的指针
// x － 用于保存方程近似解的n阶向量的指针
// d － 求解的精度控制
// n － 方程组的维数
// 返回值：0 － 如果成功。所求线性方程组的解保存在 x 中
// 1 － 求解失败
// -1 － 没有足够的自由空间
// ----------------------------------------------------------------------------
int LEquations_JacobiAlt(double A[], double b[], double x[], double d, int n)
{
int i, j, m=100/*最大迭代次数100*/;
double temp, *p;
// 用全主元素法调整方程 Ax = b，并测试是否对角占优
if ((i = SuperDiagonal(A, b, n))) return i;
// 为x(k-1)申请内存
if((p = (double *)malloc(n * sizeof(double))) == NULL) return -1;
for(i=
0; i
p[i] = 0.0; // 给出x0的初值
if (fabs(A[i*n+i]) < 1.0e-30) { free(p); return 1; }
temp = 1.0 / A[i*n+i]; // 将D-1预先乘上矩阵a,以便后续的计算简捷。
for(j=0; j
x[i] = b[i]; // 给出x1的初值
}
// 迭代过程
while(deVec_Inf(x, p, n) > d)
{
if(--m < 0){ free(p); return 1; } // 计算了100次迭代还没收敛(求解失败)
for(i=0; i
for(i=0; i
temp = 0.0;
for(j=0; j
x[i] = b[i] - temp; // 新的x分量的值保存在向量p中
}
}
free(p);
return 0; // 求解成功,所求的解保存在向量x中！
}
// ----------------------------------------------------------------------------
// 高斯-赛德尔迭代法
// 功能： 求解线性方程组 Ax = b
// 将a = L+D+U,xn+1 = -(D-L)-1Uxn + (D+L)-1b在程序中处理成：
// x = D-1(b - (a-D)x),即迭代中的新值依次替代旧值x.
// 参数： A － 指向n*n系数矩阵的指针
// b － 常数向量的指针
// x － 用于保存方程近似解的n阶向量的指针
// d － 求解的精度控制
// n － 方程组的维数
// 返回值：0 － 如果成功。所求线性方程组的解保存在 x 中
// 1 － 求解失败
// ----------------------------------------------------------------------------
int LEquations_GaussSeidelAlt(double A[], double b[], double x[], double d, int n)
{
int i, j, k, m=100/*最大迭代次数100*/;
double s, t, temp;
// 用全主元素法调整方程 Ax = b，并测试是否对角占优
if ((i = SuperDiagonal(A, b, n))) return i;
// 用 x0=(0,0,...,0)T 求x1,并计算 x1-x0 的范数
for (i=0; i
for(i=0; i
if (fabs(A[i*n+i]) < 1.0e-30) return 1;
t = 1.0 / A[i*n+i];
A[i*n+i] = temp = 0.0;
for(j=0; j
A[i*n+j] *= t;
temp += A[i*n+j] * x[j];
}
b[i] *= t;
x[i] = b[i] - temp;
t = fabs(x[i]);
if(t > s) s = t; // 计算x1的无穷向量范数
}
// 迭代过程
while(s > d)
{
if(--m < 0) return 1; // 计算了100次迭代还没收敛,认为是求解失败
s = 0.0;
for(i=0; i
temp = 0.0;
for(j=0; j
if(i != j) temp += a[i*n+j] * x[j];
}
t = b[i] - temp;
temp = fabs(t - x[i]);
if(temp > s) s = temp; // 计算(xn+1-xn)的无穷范数
x[i] = t; // 新的x分量的值仍保存在向量x中
}
}
return 0; // 求解成功,所求的解保存在向量x中！
}
// ---------
-------------------------------------------------------------------
// 超松弛法(SOR)
// 功能： 求解线性方程组 Ax = b
// 将a = L+D+U,xn+1 = (D+ωL)-1((1+ω)D-ωU)xn + ω(D+ωL)-1b在程序中处理成：
// x = (I-ωD-1a)x + ωD-1b,即迭代中的新值依次替代旧值x.
// 参数： A － 指向n*n系数矩阵的指针
// b － 常数向量的指针
// x － 用于保存方程近似解的n阶向量的指针
// w － 松弛因子,0<ω<2
// d － 求解的精度控制
// n － 方程组的维数
// 返回值：0 － 如果成功。所求线性方程组的解保存在 x 中
// 1 － 求解精度没有达到要求
// ----------------------------------------------------------------------------
int LEquations_SORAlt(double A[], double b[], double x[], double w, double d, int n)
{
int i, j, k, m=100/*最大迭代次数100*/;
double t, temp, s;
// 此处应处理矩阵a,使其主对角线的元素的绝对值在相应的列中最大,
// 以确保D-1的存在和求解精度。
// 限制松弛因子在0<ω<2的范围内,使迭代过程收敛。
if(w < 1.0e-10) w = 1.0e-10;
else if(w > 2.0 - 1.0e-10) w = 2.0 - 1.0e-10;
for (i=0; i
s = 0.0; // 给出x0的范数和x1的范数的初值
for(i=0; i
if (fabs(Ai*n+i]) < 1.0e-30) return 1;
t = 1.0 / A[i*n+i];
// 若取迭代初值x0=0,则x1 = (D+L)-1b的计算可以从(D+L)x1=b中解出。
A[i*n+i] = temp = 0.0;
for(j=0; j
A[i*n+j] *= t;
temp += A[i*n+j] * x[j];
}
b[i] *= t;
x[i] = (b[i] - temp) * w + (w - 1.0) * x[i];
t = fabs(x[i]);
if(t > s) s = t; // 计算x1的无穷向量范数
}
while(s > d)
{
if(--m < 0) return 1; // 计算了100次迭代还没收敛(求解精度没有达到要求)
s = 0.0;
for(i=0; i
temp = 0.0;
for(j=0; j
temp = fabs(t);
if(temp > s) s = temp; // 计算(xn+1-xn)的无穷范数
x[i] = t; // 新的x分量的值仍保存在向量x中
}
}
return 0; // 求解成功,所求的解保存在向量x中！
}
// ----------------------------------------------------------------------------
// Lagrange插值
// 功能： 用n个插值点 (xi,yi) i=(0,n-1) 确定n-1次多项式Pn-1(x)
// 参数： x － 插值点的x坐标
// y － 插值点的y坐标
// n － 插值点的个数
// a － 返回插值多项式Pn-1(x)系数的数组
// 返回值：0 － 如果成功
// 1 －(错误)插值点有重点
// -1 －没有足够的自由空间
// ----------------------------------------------------------------------------
int Interpolate_Lagrange(doubl
e *x, double *y, int n, double *a)
{
int i, j;
double s, t, *p;
if ((p = (double *)malloc((n+1)*sizeof(double))) == NULL) return -1;
// p=ωn(x), a=(0,0,...,0)T
p[0] = -x[0]; p[1] = 1.0; a[0] = 0.0;
for(i=1; i
p[i+1] = p[i]; a[i] = 0.0;
for (j=i; j>=0; j--) p[j] = p[j-1] - p[j] * x[i];
p[0] *= -x[i];
}
// 计算 a += ∑[yi * li(x)]
for (i=0; i
t = 1.0;
for(j=0; j
if (fabs(t) < 1.0e-30) { free(p); return 1; }
t = y[i] / t;
s = p[n];
for (j=n-1; j>=0; j--)
{
a[j] += t * s; // 用于计算 a += li(x)*yi
s = p[j] + s * x[i]; // 用于计算 ωn(x) / (x-xi)
}
}
free(p);
return 0;
}
// ----------------------------------------------------------------------------
// Newton插值
// 功能： 用n个插值点 (xi,yi) i=(0,n-1) 确定n-1次多项式Pn-1(x)
// 参数：
// x[] － 一维数组,长度为n,存放给定的n个结点的值x(i)
// y[] － 一维数组,长度为n,存放给定的n个结点的函数值y(i),
// y(i) = f(x(i)), i=0,1,...,n-1
// n － 插值点的个数
// a － 返回插值多项式Pn-1(x)系数的数组
// 返回值：0 － 如果成功
// 1 －(错误)插值点有重点
// -1 －没有足够的自由空间
// ----------------------------------------------------------------------------
int Interpolate_Newton(double *x, double *y, int n, double *a)
{
int i, j;
double t, *p, *q;
if((p = (double *)malloc((n+n)*sizeof(double))) == NULL) return -1;
q = p + n;
p[0] = 1.0; a[0] = y[0];
for(i=0; i
for(i=1; i
for(j=n-1; j>=i; j--) // 计算各阶差商
{
t = x[j] - x[j-i];
if (fabs(t) < 1.0e-30) { free(p); return 1; }
q[j] = (q[j] - q[j-1]) / t;
}
// 计算 q = ωi(x) = (x-xi-1)乘(x-x0)...(x-xi-2) 的积,
// 同时在 a 中累加 f[x0,...,xi] * ωi(x)
p[i] = p[i-1]; a[i] = q[i];
for(j=i-1; j>0; j--)
{
p[j] = p[j-1] - x[i-1] * p[j];
a[j] += p[j] * q[i];
}
p[0] *= -x[i-1];
a[0] += p[0] * q[i];
}
free(p);
return 0;
}
// ----------------------------------------------------------------------------
// Hermite(埃尔米特)插值
// 功能： 用n个插值点 {(x(i),y(i),y'(i)), i=(0,n-1)} 确定2n-1次多项式P2n-1(x)
// 参数：
// x[] － 一维数组,长度为n,存放给定的n个结点的值x(i)
// y[] － 一维数组,长度为n,存放给定的n个结点的函数值y(i),
// y(i) = f(x(i)), i=0,1,...,n-1
// dy[]－ 一维数组,长度为n,存放给定的n个结点的函数导数值y'(i),
// y'(i) = f'(x(i)), i=0,1,...,n-1
// n － 结点的个数
// a[] － 一维数组,长度为2n,存放插值多项式的系数
// 返回值：double 型,指
定的插值点t的函数近似值f(t)
//////////////////////////////////////////////////////////////////////
double Interpolate_Hermite(double x[], double y[], double dy[], int n, double a[])
{
int i, j;
double r, s, t, *p, *q, u, v;
if ((p = (double *)malloc((4*n+2)*sizeof(double))) == NULL) return -1;
q = p + 2*n;
// 计算q = ωn(x)^2, 多项式系数初始值 a = (0,0,...,0)T
a[0] = a[1] = 0.0;
q[0] = x[0] * x[0]; q[1] = -(x[i] + x[i]); q[2] = 1.0;
for(i=1; i
k = i + i; a[i+i] = a[i+i+1] = 0.0;
s = 2.0 * x[i]; t = x[i] * x[i];
q[k+2] = q[k];
q[k+1] = q[k-1] - s * q[k];
for (j=k; j>1; j--) q[j] = q[j-2] - s * q[j-1] + t * q[j];
q[1] = t * q[1] - s * q[0];
q[0] *= t;
}
// 计算 ∑[(αi(x)+βi(x))*li(x)^2]
for (i=0; i
// 先计算l'i(xi)
r = 1.0;
for(j=0; j
if (fabs(r) < 1.0e-30) { free(p); return 1; }
s = dy[i] - 2.0 * y[i] * r;
t = y[i] - x[i] * s;
// 计算 p = ωn(x)^2 / (x - xi)^2
k = n + n - 2;
u = q[k+2]; v = q[k+1];
for (j=k; j>=0; j--)
{
u = v + 2.0 * x[i];
v = q[j] - x[i] * x[i];
p[j] = u; // 保存 ωn(x)^2 / (x-xi)^2
}
u = 1.0 / (r * r);
// 计算 a += p * (s * x + t) / r^2
a[0] += t * p[0] * u;
for (j=1; j
a[k] += s * p[k] * u;
}
free(p);
return 0;
}
// ----------------------------------------------------------------------------
// 分段线性插值
// 功能： 用n个插值点 {(xi,yi), i=(0,n-1)} 确定n-1个形如aix+bi的一次多项式
// 参数： x － 插值点的x坐标
// y － 插值点的y坐标
// xp－ 自变量值(求在此处的函数值)
// n － 插值点的个数
// ----------------------------------------------------------------------------
void Spline1_GetValue(double *x, double *y, double xp, int n)
{
int i;
double t;
// 确定求值函数区间
if (xp < x[0]) i = 0;
else if (xp > x[n-1]) i = n - 2;
{
for(i=0; i
}
// 计算第i区间上在xp处的函数值
t = (y[i+1] - y[i]) / (x[i+1] - x[i]);
return (y[i] - t * (xp - x[i]));
}
// ----------------------------------------------------------------------------
// 三次样条插值(第一类边界条件)
// 功能： 用n个插值点 (xi,yi) i=(0,n-1) 确定n-1个3次插值多项式Si(x),参数用二次导数给出
// 参数： x － 插值点的x坐标数组
// y － 插值点的y坐标数组
// m － 存放结果(n维向量)
// ya－ 使用左端点的一阶导数值作为一个边界条件
// yb－ 使用右端点的一阶导数值作为另一个边界条件
// n － 插值点的个数
// 返回值：0 － 如果成功
// -1 －(错误)插值点有重点
// ---------------------------------------------------------
-------------------
int Interpolate_Spline3_1(double *x, double *y, double *m, double ya, double yb, int n)
{
int i, k=n+n;
double t, h, *a;
if((a = (double *)malloc(3*n*sizeof(double))) == NULL) return -1;
// 先处理首和尾方程的系数
a[0] = a[k+n-1] = 0.0; a[n] = a[k-1] = 2.0; a[n-1] = a[k] = 1.0;
h = x[1] - x[0]; m[0] = 6.0 * ((y[1] - y[0]) / h - ya) / h;
h = x[n-1] - x[n-2]; m[n-1] = 6.0 * (yb - (y[n-1] - y[n-2]) / h) / h;
// 计算三对角线性方程组的系数矩阵和右端向量
for(i=1; i
h = x[i] - x[i-1]; t = x[i+1] - x[i-1];
a[n+i] = 2.0; a[i] = h / t; a[k+i] = 1.0 - a[i];
m[i] = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / h) / t;
}
// 用追赶法解三对角线性方程组,并把求得的解通过m返回
LEquations_Pursue(a, m, n);
free(a);
return 0;
}
// ----------------------------------------------------------------------------
// 三次样条插值(第二类边界条件)
// 功能： 用n个插值点 (xi,yi) i=(0,n-1) 确定n-1个3次插值多项式Si(x),参数用二次导数给出
// 参数： x － 插值点的x坐标数组
// y － 插值点的y坐标数组
// m － 存放结果(n维向量)
// ya－ 使用左端点的二阶导数值作为一个边界条件
// yb－ 使用右端点的二阶导数值作为另一个边界条件
// n － 插值点的个数
// 返回值：0 － 如果成功
// -1 －(错误)插值点有重点
// ----------------------------------------------------------------------------
int Interpolate_Spline3_2(double *x, double *y, double *m, double ya, double yb, int n)
{
int i, k=n+n;
double t, h, *a;
if((a = (double *)malloc(3*n*sizeof(double))) == NULL) return -1;
// 先处理首和尾方程的系数
a[0] = a[n-1] = a[k] = a[k+j-1] = 0.0;
a[n] = a[k-1] = 1.0;
m[0] = ya; m[n-1] = yb;
// 计算三对角线性方程组的系数矩阵和右端向量
for(i=1; i
h = x[i] - x[i-1]; t = x[i+1] - x[i-1];
a[n+i] = 2.0; a[i] = h / t; a[k+i] = 1.0 - a[i];
m[i] = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / h) / t;
}
// 用追赶法解三对角线性方程组,并把求得的解通过m返回
LEquations_Pursue(a, m, n);
free(a);
return 0;
}
// ----------------------------------------------------------------------------
// 三次样条插值(第三类边界条件)
// 功能： 用n个插值点 (xi,yi) i=(0,n-1) 确定n-1个3次插值多项式Si(x),参数用二次导数给出
// 参数： x － 插值点的x坐标数组
// y － 插值点的y坐标数组
// m － 存放结果(n维向量)
// n － 插值点的个数
// 返回值：0 － 如果成功
// -1 －(错误)插值点有重点
// ----------------------------------------------------------------------------
int Interpolate_Spline3_3(double
*x, double *y, double *m, int n)
{
int i, k=n*n;
double t, h, *a;
if((a = (double *)malloc(k*sizeof(double))) == NULL) return -1;
for (i=0; i
// 先处理首和尾方程的系数
a[0] = 1.0; a[n-1] = -1.0; a[k-1] = 2.0;
h = x[1] - x[0]; t = x[n-1] - x[n-2];
a[k-n+1] = h / (h + t); a[k-2] = t / (h + t);
m[0] = 0.0;
m[n-1] = 6.0 * ((y[1] - y[0]) / h - (y[n-1] - y[n-2]) / t) / (h + t);
// 计算三对角线性方程组的系数矩阵和右端向量
for(i=1; i
h = x[i] - x[i-1]; t = x[i+1] - x[i-1];
a[i*n+i] = 2.0;
a[i*n+i-1] = h / t; a[i*n+i+1] = 1.0 - a[i*n+i-1];
m[i] = 6.0 * ((y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / h) / t;
}
// 用Gauss列主元素法解线性方程组,并把求得的解通过m返回
LEquations_Gauss(a, m, n);
free(a);
return 0;
}
// ----------------------------------------------------------------------------
// 功能： 三次样条插值多项式的求值函数。
// 参数： x － 插值点的x坐标数组
// y － 插值点的y坐标数组
// m － 三次样条插值的(插值点)二阶导数数组
// xp－ 求值点的x坐标
// n － 数组的长度(或指出是n-1次插值多项式)
// 返回值：0 － 成功
// -1 － 插值失败
// ----------------------------------------------------------------------------
double Spline3_GetValue(double *x, double *y, double *m, double xp, iint n)
{
int i, k;
double h, t, p, s;
// 确定求值函数区间
if (xp < x[0]) i = 0;
else if (xp > x[n-1]) i = n - 2;
{
for(i=0; i
}
k = i + 1;
// 使用si(xp)计算三次样条的函数值
h = x[k] - x[i];
t = xp - x[i]; s = m[k] * t * t * t;
p = xp - x[k]; s -= m[i] * p * p * p;
s += (6.0 * y[k] - m[k] * h * h) * t;
s -= (6.0 * y[i] - m[i] * h * h) * p;
s /= 6.0 * h;
return s;
}
// ----------------------------------------------------------------------------
// 功能： 对给定观测点数组(xi,yi) i=(0,m-1)用最小二乘法求一个不超过n次的多项式。
// 参数： x － 观测点的x坐标
// y － 观测点的y坐标
// m － 观测点的个数
// n － 回归多项式的次数+1
// a － 返回最小二乘多项式Pn-1(x)系数的数组
// 返回值：0 － 如果成功
// 1 － 用于回归的数据太少,所求的结果不会有太大意义,因而本程序对此不做处理
// -1 － 系统没有足够的可供计算的自由空间
// ----------------------------------------------------------------------------
int minQ(double *x, double *y, int m, int n, double *a)
{
int i, j, k = n*n;
double t, *A, *b=a, *p;
if(m >= n) return 1;
if((A = p = (double *)malloc((k+m)*sizeof(double))) == NULL) return -1;
p += k;
for(i=0; i
方程组的系数矩阵A
for(i=0; i
a[0] = (double)m;
for(j=0; j
for(i=1; i
for(j=0; j
p[j] *= x[j];
a[i] += p[j];
b[j] += p[j] * y[j];
}
}
for(i=1; i
for(j=0; j
p[j] *= x[j];
a[i*n+n-1] += p[j];
}
}
// 从a的部分值出之中产生整个系数矩阵a
for(i=1; i
for(j=0; j
k = n + i - 1;
for(j=i; j
}
LDT(A, b, n); // 用LDLT法求解法方程组
free(A);
return 0;
}
// ----------------------------------------------------------------------------
// 功能： 对给定观测点数组(xi,yi) i=(0,n-1)用最小二乘法求f(x)=abx中的参数a和b
// 参数： x － 观测点的x坐标
// y － 观测点的y坐标
// a － 获取回归方程的系数a
// b － 获取回归方程的参数b
// n － 观测点的个数
// ----------------------------------------------------------------------------
void eMinQ(double *x, double *y, double *a, double *b, int n)
{
int i;
double xt, yt, xx, xy, t;
xt = yt = xx = xy = 0.0; // 初始化为0用于求和
for(i=0; i
xt += x[i];
t = log(y[i]);
yt += t;
xx += x[i] * x[i];
xy += x[i] * t;
}
*b = ((double)n * xy - xt * yt) / ((double)n * xx - xt * xt);
*a = exp((yt - *b * xt) / (double)n);
}
// ----------------------------------------------------------------------------
// 功能： 对给定观测点数组(xi,yi) i=(0,m-1)用线性最小二乘法的一般形式求解。
// 参数： x － 观测点的x坐标
// y － 观测点的y坐标
// f － 函数(向量)f=(f0,f1,...,fm)的指针数组
// w － 权系数向量(数组的指针)
// m － 观测点的个数
// n － 回归多项式的次数+1
// a － 返回最小二乘多项式Pn-1(x)系数的数组
// 返回值：0 － 如果成功
// 1 － 用于回归的数据太少,所求的结果不会有太大意义,因而本程序对此不做处理
// -1 － 系统没有足够的可供计算的自由空间
// ----------------------------------------------------------------------------
int aMinQ(double *x, double *y, double (*f[])(double), double *w, int m, int n, double *a)
{
int i, j, k = n*n;
double t, *A, *b=a;
if((A = (double *)malloc(k*sizeof(double))) == NULL) return -1;
for(i=0; i
for(j=i; j
t = 0.0;
for(k=0; k
[j*(n+1)+i] = A[i*(n+1)+j] = t;
}
b[i] = 0.0
for(k=0; k
LDT(A, b, n); // 用LDLT法求解法方程组
free(A);
return 0;
}
// 线性最小二乘法一般形式的通用算法
// ----------------------------------------------------------------------------
// 功能： 对给定观测点数组(xi,yi) i=(0,m-1)用最小二乘法求一个不超过n次的多项式。
// 参数： x － 观测点的x坐标
// y － 观测点的y坐标
// w － 目标函数的全因子
// m － 观测点的个数
// n － 最小二乘多项式的次数+1
// a － 用于返回最小二乘多项式的系数数组
// 返回值：0 － 如果成功,结果保存在数组 a 中
// -1 － 系统没有足够的可供计算的自由空间
// ----------------------------------------------------------------------------
int minQs(double *x, double *y, double *w, int m, int n, double *a)
{
int i, k;
double *r, *u, *v, *p;
double ak, algha, bet, t, t0, t1, temp;
if(n < n ) return 1;
if((r = u = v = (double *)malloc((n+n+n) * sizeof(double))) == NULL) return -1;
u += n; v += n + n; // u和v分别存放φk-1和φk-2多项式的系数
// 初始化k=0,1时的各项参数
for(i=0; i
t0 = t1 = ak = u[0] = 0.0;
for(i=0; i
t0 += w[i]; // 计算a0和α1的分母部分的值
u[0] += w[i] * x[i]; // 计算α1的分子部分的值
r[0] += w[i] * y[i]; // 计算a0的分子部分的值
}
u[0] = -u[0] / t0; u[1] = v[0] = 1.0; // u=φ1(x), v=φ0(x)
r[0] /= t0; // r=a0*φ0(x)
for(i=0; i
temp = (x[i] + u[0]);
ak += w[i] * y[i] * temp;
t1 += w[i] * temp * temp;
}
ak /= t1; r[0] += u[0] * ak; // (a2) = ak = (y,φ) / (φ1,φ1)
r[1] = ak; // r = a0*φ0(x) + a1*φ1(x)
// 主循环,u=(φk-1)和v(=φk-2)及合并后的多项式r(=P)的系数
for(k=2; k
alpha = 0.0;
for(i=0; i
temp = Polynomial(u, k, x[i]);
alpha += w[i] * x[i] * temp * temp;
}
alpha /= t1; // 计算αk的值
bet = t1 / t0; // 计算βk的值
// 求φk多项式的系数,并保存到v中
v[k] = u[k-1];
v[k-1] = u[k-2] - alpha * u[k-1];
for(i=k-2; i>0; i--)
v[i] = u[i-1] - alpha * u[i] - bet * v[i];
v[0] = -alpha * u[0] - bet * v[0];
p = u; u = v; v = p; // 保持u和v仍分别存放φk-1和φk-2多项式的系数
t = ak = 0.0;
for(i=0; i
temp = Polynomial(u, k+1, x[i]);
t += w[i] * temp * temp;
ak += w[i] * y[i] * temp;
}
ak /= t; t0 = t1; t1 = t;
r[k] = ak;
for(i=0; i
for(i=0; i
free(r);
return 0;
}
// -------------------------------------------
---------------------------------
// 用Legendre多项式求最佳平方逼近多项式
// 功能： 求f(x)在[a,b]上的最佳平方逼近n-1次多项式P(x)
// 参数： f － 函数f(x)(的指针)
// a － 区间的左端点
// b － 区间的右端点
// c － 返回多项式系数数组
// n － 返回多项式系数数组的长度
// d － 求解精度控制参数
// 返回值：0 － 如果成功,多项式P(x)保存在数组 c 中
// -1 － 系统没有足够的可供计算的自由空间
// ----------------------------------------------------------------------------
int Approach_Legendre(double (*f)(double), double a, double b, double d, double *c, int n)
{
int i, k;
double *u, *v, *p, s, t, r, ak, r1, r2;
if ((u = (double *)malloc((n+n)*sizeof(double))) == NULL) return -1;
v = u + n;
r = 0.5;
s = 2.0 / (b - a); t = (a + b) / (a - b);
v[0] = 1.0; // P0(x)
u[0] = t; u[1] = s; // P1(x)
c[0] = r * Integral_Romberg(f, a, b, d); // φ0(x)
r += 1.0;
ak = r * Integral_Legendre(f, u, 2, a, b, d);
c[0] += t * ak; c[1] = s * ak; // φ(x) += φ1(x)
for (k=2; k
r2 = 1.0 / (double)k; r1 = 1.0 - r2; r2 = 2.0 - r2;
// 计算 Pk(x)
v[k] = u[k-1] * r2 * s;
v[k-1] = r2 * (u[k-2] * s + u[k-1] * t);
for (i=k-2; i>0; i++)
v[i] = r2 * (u[i-1] * s + u[i] * t) + r1 * v[i];
v[0] = r2 * u[0] * t + r1 * v[0];
p = v; v = u; u = p; // 交换, 使 u = φk(x); v = φk-1(x)
// 计算ak
r += 1.0;
ak = r * Integral_Legendre(f, u, k+1, a, b, d);
// φ(x) += φk(x)
for (i=0; i
c[k] = ak * u[k];
}
if (v < u) u = v;
free(u);
return 0;
}
// ----------------------------------------------------------------------------
// 功能： Legendre(勒让德)多项式求值。
// 参数： x － 自变量x
// n － 多项式的次数
// 返回值：Pn(x)的值
// ----------------------------------------------------------------------------
double Legendre(double x, int n)
{
int i;
double p, p0=1.0, p1=x;
if(n <= 0) return 1.0;
else if(n == 1) return x;
for(i=2; i
p = ((double)(i+i-1) * x * p1 - (double)(i-1) * p0) / (double)i;
p0 = p1; p1 = p;
}
return p;
}
// ----------------------------------------------------------------------------
// 功能： Chebyshev(切比雪夫)多项式求值。
// 参数： x － 自变量x
// n － 多项式的次数
// 返回值：Tn(x)的值
// ----------------------------------------------------------------------------
double Chebyshev(double x, int n)
{
int i;
double t, t0=1.0, t1=x;
if(n <= 0) return 1.0;
else if(n == 1) return x;
for(i=2; i
t = 2.0 * x * t1 - t0;
t0 = t1; t1 = t;
}
return t;
}
// -------------------------------------------
---------------------------------
// 功能： Laguerre(拉盖尔)多项式求值。
// 参数： x － 自变量x
// n － 多项式的次数
// 返回值：Ln(x)的值
// ----------------------------------------------------------------------------
double Laguerre(double x, int n)
{
int i;
double L, L0 = 1.0, L1 = 1.0 - x;
if(n <= 0) return 1.0;
else if(n == 1) return x;
for(i=2; i
L = ((double)(i+i-1)-x) * L1 - x * x * L0;
L0 = L1; L1 = L;
}
return p;
}
// ----------------------------------------------------------------------------
// 功能： Hermite(艾尔米特)多项式求值。
// 参数： x － 自变量x
// n － 多项式的次数
// 返回值：Hn(x)的值
// ----------------------------------------------------------------------------
double Laguerre(double x, int n)
{
int i;
double H, H0 = 1.0, H1 = x + x;
if(n <= 0) return 1.0;
else if(n == 1) return x;
for(i=2; i
H = 2.0 * (x * H1 - (double)(i-1) * H0;
H0 = H1; H1 = H;
}
return p;
}
// ----------------------------------------------------------------------------
// 功能： Chebyshev最佳一致逼近n次多项式求法。
// 参数： f － [a,b]上的逼近目标函数的指针
// x － (返回值)插值点的x坐标,实参必须保证能存放n个double数据
// y － (返回值)Newton插值多项式的系数,实参必须保证能存放n个double数据
// xa－ f定义的(用于逼近的)区间左端点
// xb－ f定义的(用于逼近的)区间右端点
// n － 逼近多项式的次数
// a － 返回插值多项式Pn-1(x)系数的数组
// 返回值：0 － 如果成功。结果通过x和y返回插值点列，a返回插值多项式系数
// 1 －(错误)插值点有重点
// -1 －没有足够的自由空间
// ----------------------------------------------------------------------------
int ChebyshevNewton(double (*f)(double), double *x, double *y,
double xa, double xb, int n, double a[])
{
int i, j;
double t=(xa+xb)* 0.5, p=(xb-xa)* 0.5;
for(i=0; i
x[i] = t + p * cos((double)(i+i+1) * 3.14159265358979323846 / (double)(n+n+2));
y[i] = (*f)(x[i]);
}
return (Interpolate_Newton(x, y, n, a));
}
// ----------------------------------------------------------------------------
// 功能： 求多项式函数的导函数。
// 参数： a － 多项式系数数组(也存放求导结果)
// n － 数组的长度(或指出是n-1次多项式)
// 返回值：无
// ----------------------------------------------------------------------------
void diff(double *a, int n)
{
int i;
for(i=1; i
a[i-1] = a[i] * (double)i;
}
}
(2) 给定一组离散点(xi,yi)(i=0,1,...,n-1),通过插值多项
式求导函数
extern int Newton(double *x, double *y, int n);
extern int tNewton(double *x, double *y, int n);
// ----------------------------------------------------------------------------
// 功能： 通过Newton插值多项式求导函数。
// 参数： x － 插值点的x坐标数组
// y － 插值点的y坐标数组,(也存放求导结果多项式的系数数组)
// n － 数组的长度(或指出是n-1次插值多项式)
// 返回值：0 － 成功
// -1 － 插值失败
// ----------------------------------------------------------------------------
int nDiff(double *x, double *y, int n)
{
if(!Newton(x,y,n)) return -1;
if(!tNewton(x,y,n)) return -1;
diff(y, n);
return 0;
}
extern int spline(double *x, double *y, double *m, double ya, double yb, int n);
// ----------------------------------------------------------------------------
// 功能： 三次样条插值多项式导函数的求值函数。
// 参数： x － 插值点的x坐标数组
// y － 插值点的y坐标数组
// m － 三次样条插值的(插值点)二阶导数数组
// xp－ 求值点的x坐标
// n － 数组的长度(或指出是n-1次插值多项式)
// 返回值：0 － 成功
// -1 － 插值失败
// ----------------------------------------------------------------------------
double Val_dSpline(double *x, double *y, double *m, double xp, iint n)
{
int i, k;
double h, t, s;
// 确定求值函数区间
for(i=1; i
k = i - 1;
h = x[i] - x[k];
t = xp - x[k]; s = m[i] * t * t;
t = xp - x[i]; s -= m[k] * t * t;
t = y[i] - y[k]; s += t + t;
s /= h + h;
s -= h * (m[i] - m[k]) / 6.0;
return s;
}
// ----------------------------------------------------------------------------
// 复化梯形求积法
// 功能： 求f(x)在[a,b]上的定积分
// 参数： f － 积分函数(的指针)
// a － 积分区间的左端点
// b － 积分区间的右端点
// d－ 求解精度控制参数
// 返回值：f(x)在[a,b]上的积分
// ----------------------------------------------------------------------------
double Integral_Trapezoid(double (*f)(double), double a, double b, double d)
{
int i, n = 2;
double s, ss, t, x, h = b - a;
t = ((*f)(a) + (*f)(b)) * 0.5; s = h * t;
h *= 0.5;
t += (*f)((a+b)*0.5); ss = h * t;
while(fabs(ss-s) > d)
{
s = ss;
x = a + h * 0.5;
for(i=0; i
t += (*f)(x);
x += h;
}
h *= 0.5; n <<= 1;
ss = t * h;
}
return ss;
}
// ----------------------------------------------------------------------------
// 复化辛普森法求积法
// 功能： 求f(x)在[a,b]上的定积分
// 参数： f － 积分函数(的指针)
// a － 积分区间的左端点
// b － 积分区间的右端点
//
d－ 求解精度控制参数
// 返回值：f(x)在[a,b]上的积分
// ----------------------------------------------------------------------------
double Integral_Simpson(double (*f)(double), double a, double b, double d)
{
int i, n=2;
double s, ss, t, tt, p, h, x;
h = b - a;
t = (*f)(a) + (*f)(b);
p = (*f)( (a+b)* 0.5 );
tt = p + p;
ss = (t + tt + tt) * h / 6.0;
s = ss + d + d;
while(fabs(ss-s) > d)
{
p = 0; t += tt;
h *= 0.5; x = a + h * 0.5;
for(i=0; i
tt = p + p; n <<= 1;
s = ss;
ss = (t + tt + tt) * h / 6.0;
}
return ss;
}
// ----------------------------------------------------------------------------
// 龙贝格求积法
// 功能： 求f(x)在[a,b]上的定积分
// 参数： f － 积分函数(的指针)
// a － 积分区间的左端点
// b － 积分区间的右端点
// d － 求解精度控制参数
// 返回值：f(x)在[a,b]上的积分
// ----------------------------------------------------------------------------
double Integral_Romberg(double (*f)(double), double a, double b, double d)
{
int i, n=4;
double buff[8], *t, *tt, *r, temp, p, h, x;
t = buff; tt = &buff[4];
p = (*f)(a) + (*f)(b);
h = (b - a) * 0.5;
tt[0] = p * h; ; 计算 T1
p += 2 * (*f)(a+h); h *= 0.5;
t[0] = p * h; ; 计算 T2
t[1] = (4.0 * t[0] - tt[0]) / 3.0; ; 计算 S1
p += 2 * ((*f)(a+h) + (*f)(b-h));
h /= 2.0;
tt[0] = p * h; ; 计算 T4
tt[1] = ( 4.0 * tt[0] - t[0]) / 3.0; ; 计算 S2
tt[2] = (16.0 * tt[1] - t[1]) / 15.0; ; 计算 C1
tt[3] = tt[2]; t[3] = t[1]; ; 为计算和收敛判断准备初值
while(fabs(tt[3]-t[3]) > d && n)
{
temp = 0; x = a + h;
for(i=0; i
n <<= 1; h /= 2.0;
t[0] = p * h; ; 计算 T8k, k=n/8
t[1] = ( 4.0 * t[0] - tt[0]) / 3.0; ; 计算 S4k
t[2] = (16.0 * t[1] - tt[1]) / 15.0; ; 计算 C2k
t[3] = (64.0 * t[2] - tt[2]) / 63.0; ; 计算 Rk
r = t; t = tt; tt = r; ; 交换列
}
return rr[0];
}
// ----------------------------------------------------------------------------
// 含有多项式的龙贝格求积法
// 功能： 求P(x)*f(x)在[a,b]上的定积分(其中P(x)是n-1次多项式)
// 参数： f － 积分函数(的指针)
// c － 多项式系数数组
// n － 多项式系数数组的长度
// a － 积分区间的左端点
// b － 积分区间的右端点
// d － 求解精度控制参数
// 返回值：f(x)在[a,b]上的积分
// ----------------------------------------------------------------------------
double Integral_Legendre(double (*f)(double), double *c, int n, double a, double b, double d)
{
int i, m=4;
double buff[8], *t, *tt, *r, temp, p, h, x;
t = buff; tt = &buff[4];
p = (*f)(a)
* Polynomial(c, n, b) + (*f)(b) * Polynomial(c, n, b);
h = (b - a) * 0.5;
tt[0] = p * h; ; 计算 T1
p += 2.0 * (*f)(a+h) * Polynomial(c, n, a+h);
h *= 0.5;
t[0] = p * h; ; 计算 T2
t[1] = (4.0 * t[0] - tt[0]) / 3.0; ; 计算 S1
p += 2 * ((*f)(a+h) * Polynomial(c, n, a+h) + (*f)(b-h) * Polynomial(c, n, b-h));
h /= 2.0;
tt[0] = p * h; ; 计算 T4
tt[1] = ( 4.0 * tt[0] - t[0]) / 3.0; ; 计算 S2
tt[2] = (16.0 * tt[1] - t[1]) / 15.0; ; 计算 C1
tt[3] = tt[2]; t[3] = t[1]; ; 为计算和收敛判断准备初值
while(fabs(tt[3]-t[3]) > d && m)
{
temp = 0; x = a + h;
for(i=0; i
temp += (*f)(x) * Polynomial(c, n, x);
x += h + h;
}
p += temp + temp;
m <<= 1; h /= 2.0;
t[0] = p * h; ; 计算 T8k (k=n/8)
t[1] = ( 4.0 * t[0] - tt[0]) / 3.0; ; 计算 S4k
t[2] = (16.0 * t[1] - tt[1]) / 15.0; ; 计算 C2k
t[3] = (64.0 * t[2] - tt[2]) / 63.0; ; 计算 Rk
r = t; t = tt; tt = r; ; 交换列
}
return rr[0];
}
// ----------------------------------------------------------------------------
// 二分法
// 功能： 求解非线性方程f(x)=0在[a,b]区间上的根,并返回所求的结果。
// 参数： double f － 指向方程f(x)=0中的函数f的指针
// double a － 求解区间的一个端点
// double b － 求解区间的另一个端点
// double d － 求解的精度(误差限)
// 返回值：double,成功时返回方程的根,否则返回最大的双精度数值。
// ----------------------------------------------------------------------------
double NLEquation_Split(double (*f)(double), double a, double b, double d)
{
double r, s, t, p;
s = (*f)(a); t = (*f)(b);
// 让a对应函数值为负值的端点坐标,b对应函数值为正数的端点坐标,
// 以便使二分法主循环内部的计算逻辑更简单些。
if(s < 0.0)
{
if(t < 0.0) return DBL_MAX; // 求解失败
// DBL_MAX 为最大的双精度数：1.7976931348623158e+308
// 该常量(宏)定义在头文件 double.h 中
} else if(s > 0.0) {
if(t > 0.0) return DBL_MAX; // 求解失败
r = s; s = t; t = r;
r = a; a = b; b = r;
} else return a;
// 二分法主循环
while(fabs(b-a) > d)
{
p = (a + b) * 0.5; r = (*f)(p);
if(p < 0.0) { a = p; s = r; }
else { b = p; t = r; }
}
// 返回最终所确定区间的中值
return ((a + b) * 0.5);
}
// ----------------------------------------------------------------------------
// 简单迭代法
// 功能： 求解非线性方程f(x)=0在x附近的一个根,并返回所求的结果。
// 参数： double f － 指向方程f(x)=0中的函数f的指针
// double x － 初始迭代点
// double d － 求解的精度(误差限)
// int N － 最大迭代次数限制
// 返回值：double,成功时返回方程的根,否则返回最大的双精度数值。
// ----------------------------------------------------------------------------
double NLEquation_Alternate(double (*f)(double), double x, double d, int N)
{
double t;
while (N-- > 0)
{
t = (*f)(x);
if (fabs(t-x) < d) return t; // 返回方程的根
x = t;
}
return DBL_MAX; // 达到了最大迭代次数
}
// ----------------------------------------------------------------------------
// 牛顿迭代法
// 功能： 求解非线性方程f(x)=0在x附近的一个根,并返回所求的结果。
// 参数： double f － 指向方程f(x)=0中的函数f的指针
// double f1－ 指向方程f(x)=0中的函数f的导函数f'的指针
// double x － 初始迭代点
// double d － 求解的精度(误差限)
// int N － 最大迭代次数限制
// 返回值：double,成功时返回方程的根,否则返回最大的双精度数值。
// ----------------------------------------------------------------------------
double NLEquation_Newton(double (*f)(double), double (*f1)(double), double x, double d, int N)
{
double t, s;
while (N-- > 0)
{
s = (*f1)(x);
if (fabs(s) < 1.0e-30) break;
t = x - (*f)(x) / s;
if (fabs(t-x) < d) return t; // 返回方程的根
x = t;
}
return DBL_MAX; // 达到了最大迭代次数
}
// ----------------------------------------------------------------------------
// 简化的牛顿迭代法
// 功能： 求解非线性方程f(x)=0在x附近的一个根,并返回所求的结果。
// 参数： double f － 指向方程f(x)=0中的函数f的指针
// double f1－ 1 / f'(x0)
// double x － 初始迭代点
// double d － 求解的精度(误差限)
// int N － 最大迭代次数限制
// 返回值：double,成功时返回方程的根,否则返回最大的双精度数值。
// ----------------------------------------------------------------------------
double NLEquation_NewtonSimple(double (*f)(double), double f1, double x, double d, int N)
{
double t;
while (N-- > 0)
{
t = x - (*f)(x) * f1;
if (fabs(t-x) < d) return t; // 返回方程的根
x = t;
}
return DBL_MAX; // 达到了最大迭代次数
}
// ----------------------------------------------------------------------------
// 牛顿下山法
// 功能： 求解非线性方程f(x)=0在x附近的一个根,并返回所求的结果。
// 参数： double f － 指向方程f(x)=0中的函数f的指针
// double f1－ 指向方程f(x)=0中的函数f的导函数f'的指针
// double x － 初始迭代点
// double d － 求解的精度(误差限)
// int N － 最大迭代次数限制
// 返回值：double,成功时返回方程的根,否则返回最大的双精度数值。
// ----------------------------------------------------------------------------
double NLEquation_NewtonDown(dou
ble (*f)(double), double (*f1)(double), double x, double d, int N)
{
double t, s, r, fx, alpha;
fx = (*f)(x); s = (*f1)(x);
if (fabs(s) < 1.0e-30) return DBL_MAX;
alpha = -fx / s;
while (N-- > 0)
{
t = x + alpha;
r = (*f)(t);
if (fabs(r) > fabs(fx)) alpha *= 0.5;
else
{
if (fabs(t-x) < d)
{
if (fabs(r) < d) return t; // 返回方程的根
break;
}
x = t; fx = r;
s = (*f1)(x);
if (fabs(s) < 1.0e-30) break;
alpha = -fx / s;
}
}
return DBL_MAX; // 达到了最大迭代次数
}
// ----------------------------------------------------------------------------
// 埃特金加速方法(也叫Steffensen方法)
// 功能： 求解非线性方程f(x)=0在x附近的一个根。
// 参数： double f － 指向方程f(x)=0中的函数f的指针
// double x － 求解区间的一个端点
// double d － 求解的精度(误差限)
// int N － 最大迭代次数限制
// 返回值：double,成功时返回方程的根,否则返回最大的双精度数值。
// ----------------------------------------------------------------------------
double NLEquation_Aitken(double (*f)(double), double x, double d, int N)
{
double s, t, r;
while (N-- > 0)
{
s = (*f)(x); t = (*f)(s);
r = t - s - s + x;
if (fabs(r) < 1.0e-30) break;
t = x - (s-x)*(s-x) / r; // 埃特金加速
if (fabs(t-x) < d) return t; // 返回方程的根
x = t;
}
return DBL_MAX; // 达到了最大迭代次数
}
// ----------------------------------------------------------------------------
// 弦截法
// 功能： 求解非线性方程f(x)=0在[a,b]区间上的根,并返回所求的结果。
// 参数： double f － 指向方程f(x)=0中的函数f的指针
// double a － 求解区间的一个端点
// double b － 求解区间的另一个端点
// double d － 求解的精度(误差限)
// int N － 最大迭代次数限制
// 返回值：double,成功时返回方程的根,否则返回最大的双精度数值。
// ----------------------------------------------------------------------------
double NLEquation_Bowcut(double (*f)(double), double a, double b, double d, int N)
{
double t, s, fa, fb;
fb = (*f)(a);
while (N-- > 0)
{
fa = fb; fb = (*f)(b);
s = fb - fa;
if (fabs(s) < 1.0e-30) break;
t = b - (b - a) * fb / s;
a = b; b = t;
if (fabs(b-a) < d) return t; // 返回方程的根
}
return DBL_MAX; // 达到了最大迭代次数
}
// ----------------------------------------------------------------------------
// Muller法(抛物线法)
// 功能： 求解非线性方程f(x)=0在x附近的一个根,并返回所求的结果。
// 参数： double f － 指向方程f(x)=0中的函数f的指针
// double x0－ 初始第一个迭代点
// double x1－ 初始第二个迭代点
//