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]; }