/************************* Program 2.2 ****************************/ /* */ /************************************************************************/ /* Please Note: */ /* */ /* (1) This computer program is written by Tao Pang in conjunction with */ /* his book, "An Introduction to Computational Physics," published */ /* by Cambridge University Press in 1997. */ /* */ /* (2) No warranties, express or implied, are made for this program. */ /* */ /************************************************************************/ #define NMAX 21 void pfit (n,m,x,f,a,u) /* Subroutine generating orthonormal polynomials U(M,N) up to (M-1)th order and coefficients A(M), for the least squares approximation of the function F(N) at X(N). Other variables used: G(K) for g_k, H(K) for h_k, S(K) for . Copyright (c) Tao Pang 1997. */ int n,m; double x[],f[]; double a[],u[NMAX][NMAX]; { int i,j; double stmp; double s[NMAX],g[NMAX],h[NMAX]; if (n > NMAX) { printf ("Too many points\n"); exit (1); } if (m > n) { printf ("Order too high\n"); exit (2); } /* Set up the zeroth order polynomial u_0 */ s[0] = 0; g[0] = 0; a[0] = 0; for (i = 0; i < n; ++i) { u[0][i] = 1; stmp = u[0][i]*u[0][i]; s[0] = s[0]+stmp; g[0] = g[0]+x[i]*stmp; a[0] = a[0]+u[0][i]*f[i]; } g[0] = g[0]/s[0]; a[0] = a[0]/s[0]; h[0] = 0; /* Set up the first order polynomial u_1 */ s[1] = 0; g[1] = 0; h[1] = 0; a[1] = 0; for (i = 0; i < n; ++i) { u[1][i] = x[i]*u[0][i]-g[0]*u[0][i]; s[1] = s[1]+u[1][i]*u[1][i]; g[1] = g[1]+x[i]*u[1][i]*u[1][i]; h[1] = h[1]+x[i]*u[1][i]*u[0][i]; a[1] = a[1]+u[1][i]*f[i]; } g[1] = g[1]/s[1]; a[1] = a[1]/s[1]; h[1] = h[1]/s[0]; /* Higher order polynomials u_k from the recursive relation */ if (m >= 3) { for (i = 1; i < m-1; ++i) { s[i+1] = 0; g[i+1] = 0; h[i+1] = 0; a[i+1] = 0; for (j = 0; j < n; ++j) { u[i+1][j] = x[j]*u[i][j]-g[i]*u[i][j]-h[i]*u[i-1][j]; s[i+1] = s[i+1]+u[i+1][j]*u[i+1][j]; g[i+1] = g[i+1]+x[j]*u[i+1][j]*u[i+1][j]; h[i+1] = h[i+1]+x[j]*u[i+1][j]*u[i][j]; a[i+1] = a[i+1]+u[i+1][j]*f[j]; } g[i+1] = g[i+1]/s[i+1]; a[i+1] = a[i+1]/s[i+1]; h[i+1] = h[i+1]/s[i]; } } }