5 、 数值计算

5.1 定积分计算(Romberg)

/* Romberg 求定积分
输入：积分区间[a,b]，被积函数 f(x,y,z)
输出：积分结果
f(x,y,z)示例：
double f0( double x, double l, double t )
{
return sqrt(1.0+l*l*t*t*cos(t*x)*cos(t*x));
}
*/
double Integral(double a, double b, double (*f)(double x, double y, double z), double eps,
double l, double t)
double Romberg (double a, double b, double (*f)(double x, double y, double z), double eps,
double l, double t)
{
#define MAX_N 1000
int i, j, temp2, min;
double h, R[2][MAX_N], temp4;
for (i=0; i<MAX_N; i++) {
R[0][i] = 0.0;
R[1][i] = 0.0;
}
h = b-a;
min = (int)(log(h*10.0)/log(2.0)); //h should be at most 0.1
R[0][0] = ((*f)(a, l, t)+(*f)(b, l, t))*h*0.50;
i = 1;
temp2 = 1;
while (i<MAX_N){
i++;
R[1][0] = 0.0;
for (j=1; j<=temp2; j++)
R[1][0] += (*f)(a+h*((double)j-0.50), l, t);
R[1][0] = (R[0][0] + h*R[1][0])*0.50;
temp4 = 4.0;
for (j=1; j<i; j++) {
R[1][j] = R[1][j-1] + (R[1][j-1]-R[0][j-1])/(temp4-1.0);
temp4 *= 4.0;
}
if ((fabs(R[1][i-1]-R[0][i-2])<eps)&&(i>min))
return R[1][i-1];
h *= 0.50;
temp2 *= 2;
for (j=0; j<i; j++)
R[0][j] = R[1][j];
}
return R[1][MAX_N-1];
}
double Integral(double a, double b, double (*f)(double x, double y, double z), double eps,
double l, double t)
{
#define pi 3.1415926535897932
int n;
double R, p, res;
n = (int)(floor)(b * t * 0.50 / pi);
p = 2.0 * pi / t;
res = b - (double)n * p;
if (n)
R = Romberg (a, p, f0, eps/(double)n, l, t);
R = R * (double)n + Romberg( 0.0, res, f0, eps, l, t );
return R/100.0;
}

5.2 多项式求根（牛顿法）

/* 牛顿法解多项式的根
输入：多项式系数 c[]，多项式度数 n，求在[a,b]间的根
输出：根
要求保证[a,b]间有根
*/
double fabs( double x )
{
return (x<0)? -x : x;
}
double f(int m, double c[], double x)
{
int i;
double p = c[m];
for (i=m; i>0; i--)
p = p*x + c[i-1];
return p;
}
int newton(double x0, double *r,
double c[], double cp[], int n,
double a, double b, double eps)
{
int MAX_ITERATION = 1000;
int i = 1;
double x1, x2, fp, eps2 = eps/10.0;
x1 = x0;
while (i < MAX_ITERATION) {
x2 = f(n, c, x1);
fp = f(n-1, cp, x1);
if ((fabs(fp)<0.000000001) && (fabs(x2)>1.0))
return 0;
x2 = x1 - x2/fp;
if (fabs(x1-x2)<eps2) {
if (x2<a || x2>b)
73
return 0;
*r = x2;
return 1;
}
x1 = x2;
i++;
}
return 0;
}
double Polynomial_Root(double c[], int n, double a, double b, double eps)
{
double *cp;
int i;
double root;
cp = (double *)calloc(n, sizeof(double));
for (i=n-1; i>=0; i--) {
cp[i] = (i+1)*c[i+1];
}
if (a>b) {
root = a; a = b; b = root;
}
if ((!newton(a, &root, c, cp, n, a, b, eps)) &&
(!newton(b, &root, c, cp, n, a, b, eps)))
newton((a+b)*0.5, &root, c, cp, n, a, b, eps);
free(cp);
if (fabs(root)<eps)
return fabs(root);
else
return root;
}

5.3 周期性方程（追赶法）

/* 追赶法解周期性方程
周期性方程定义：| a1 b1 c1 ... | = x1
| a2 b2 c2 ... | = x2
| ... | * X = ...
| cn-1 ... an-1 bn-1 | = xn-1
| bn cn an | = xn
输入：a[],b[],c[],x[]
输出：求解结果 X 在 x[]中
*/
void run()
{
c[0] /= b[0]; a[0] /= b[0]; x[0] /= b[0];
for (int i = 1; i < N - 1; i ++) {
double temp = b[i] - a[i] * c[i - 1];
c[i] /= temp;
x[i] = (x[i] - a[i] * x[i - 1]) / temp;
a[i] = -a[i] * a[i - 1] / temp;
}
a[N - 2] = -a[N - 2] - c[N - 2];
for (int i = N - 3; i >= 0; i --) {
a[i] = -a[i] - c[i] * a[i + 1];
x[i] -= c[i] * x[i + 1];
}
x[N - 1] -= (c[N - 1] * x[0] + a[N - 1] * x[N - 2]);
x[N - 1] /= (c[N - 1] * a[0] + a[N - 1] * a[N - 2] + b[N - 1]);
for (int i = N - 2; i >= 0; i --)
x[i] += a[i] * x[N - 1];
}