Undetermined Coefficients

In previous sections, we approximated the definite integral and derivatives of any given function $y=f(x)$ based on its Lagrange interpolation $L_n(x)$:

$$\frac{d}{dx}f(x_j) \approx \frac{d}{dx}L_n(x_j) =\frac{d}{dx}\left[\sum_{i=0}^n f(x_i)l_i(x_j)\right]$$

$$= \sum_{i=0}^n f(x_i)\frac{d}{dx}l_i(x_j) =\sum_{i=0}^n d_{ij} f(x_i),\;\;\;\;\;(j=0,\cdots,n)$$ (79)

$$\int_a^b f(x)\,dx \approx \int_a^b L_n(x)\,dx =\int_a^b \left[\sum_{i=0}^n f(x_i)l_i(x)\right]\,dx$$

$$= \sum_{i=0}^n f(x_i)\int_a^b l_i(x)\,dx =\sum_{i=0}^n c_if(x_i)$$ (80)

whre

$$d_{ij}=d_i(x)\big\vert _{x=x_j}=\frac{d}{dx}l_i(x)\big\vert _{x=x_j} =\frac{d}{dx}l_i(x_j)$$ (81)

is the coefficient $d_i(x)=l'_i(x)$ in Eq. (12) evaluated at $x=x_j,\;(i,j=0,\cdots,n)$, and $c_i=\int_a^b l_i(x)\,dx$ is given in Eq. (25).

These coefficients can also be obtained by an alternative method of undetermined coefficients. As both $c_i$ and $d_{ij}$ are independent of the specific function $f(x)$ in question, we can find them based on a set of particular functions $f(x)=x^k,\;(k=0,\cdots,n)$:

$$I[f]=\int_a^b f(x)\,dx = \int_a^b x^k\,dx=\frac{b^{k+1}-a^{k+1}}{k+1} =\sum_{i=0}^n c_if(x_i)=\sum_{i=0}^n c_ix_i^k \;\;\;\;\;\;(k=0,\cdots,n)$$ (82)

$$D[f]=\frac{d}{dx}f(x_j) = \frac{d}{dx}x^k_j=k x^{k-1}_j =\sum_{i=0}^n d_{ij} f(x_i)=\sum_{i=0}^n d_{ij}x_i^k \;\;\;\;\;\;(k,j=0,\cdots,n)$$ (83)

We consider $x_i^k$ as the component in the kth row and ith column of a $(n+1)\times (n+1)$ matrix ${\bf A}$ so that the two equations can be expressed in matrix forms as:

$$\left[\begin{array}{cccc}1&1&\cdots&1\\ x_0&x_1&\cdots&x_n\\ x_0... ...b^3-a^3}{3}\\ \vdots\\ \frac{b^{n+1}-a^{n+1}}{n+1}\end{array}\right] \;\;\;\;\; \;\;\;\;\;\; {\bf Ac}={\bf b}$$ (84)

$$\left[\begin{array}{cccc}1&1&\cdots&1\\ x_0&x_1&\cdots&x_n\\ x_0... ...dots&\ddots&\vdots\\ nx_0^{n-1}&nx_1^{n-1}&\cdots&nx_n^{n-1}\end{array}\right]$$

$$\;\;\;\;\; \;\;\;\;\;\; {\bf AD}={\bf A}[{\bf d}_0,\cdots,{\bf d}_n]={\bf B} =[{\bf b}_0,\cdots,{\bf b}_n]$$ (85)

In the second equation, the jth column ${\bf d}_j=[d_{0j},\cdots,d_{nj}]^T$ of ${\bf D}=[{\bf d}_0,\cdots,{\bf d}_n]$ is composed of the $n+1$ coefficients for approximating $f'(x_j)$. Solving these equation systems we get the coefficients needed:

$${\bf c} = {\bf A}^{-1}{\bf b}$$ (86)

$${\bf D} = {\bf A}^{-1}{\bf B}, \;\;\;\; \;\;\;\; {\bf d}_j={\bf A}^{-1}{\bf b}_j\;\;\;\;\;\;(j=0,\cdots,n)$$ (87)

Once ${\bf c}$ and ${\bf D}$ are obtained, the integral and derivatives of any function $f(x)$ represented by its $n+1$ samples in ${\bf f}=[f(x_0),\cdots,f(x_n)]^T$ can be approximated as

$$\int_a^b f(x)\,dx = \sum_{i=0}^nc_if(x_i)={\bf c}^T{\bf f}$$ (88)

$$\frac{d}{dx}f(x_j) = \sum_{i=0}^nd_{ij}f(x_i)={\bf d}_j^T{\bf f} \;\;\;\;\;(j=0,\cdots,n)$$ (89)

This method of undetermined coefficients is implemented by the Matlab code below:

    A=zeros(n+1);     % coefficient matrix on the left
    b=zeros(n+1,1);   % vector on the right for integral
    B=zeros(n+1);     % matrix on the right for derivatives
    A(1,:)=1;         % first row of A
    B(1,:)=0;         % first row of B
    b(1,:)=b-a;       % first element of b
    for i=1:n         % for the remaining n rows
        A(i+1,:)=x.^i;
        B(i+1,:)=i*x.^(i-1);
        b(i+1)=(b^(i+1)-a^(i+1))/(i+1);
    end
    c=inv(A)*b;       % find coefficients for approximating integral
    D=inv(A)*B;       % find coefficients for approximating derivatives
    f(x)*c;           % approximated integral
    f(x)*D(:,i+1);    % approximated derivatives

Example: Find $f'(x_i)$ and $\int_0^1 f(x)\,dx$ of $f(x)=e^x$, with groun truth:

$$\int_0^1 e^x\,dx=e^x\big\vert _0^1=e-1=1.718281828,\;\;\;\;\;\; \frac{d}{dx}e^x=e^x=2.718281828^x$$ (90)

When $n=5$, we have $a=0,\,b=1,\,n=5,\;h=(b-a)/n=1/5$, and the approximation of the integral generated by the code above is $1.71828231$ with relative error $2.82e-07$. The approximated derivatives at $x_0=a,\cdots,x_n=x_5=1$ are:

$$\begin{tabular}{c\vert\vert c\vert c\vert r}\hline j & \mbox{appr... ... & 8.31e-06\\ \hline 5 & 2.718187 & 2.718282 & -3.50e-05\\ \hline \end{tabular}$$

When $n=7$, $h=(b-a)/n=1/7$, we get the approximated integral $1.71828183$ with error $3.79e-10$, and the approximated derivatives at $x_0=a,\cdots,x_n=x_7=1$ are:

$$\begin{tabular}{c\vert\vert c\vert c\vert r}\hline j & \mbox{appr... ... & 1.59e-08\\ \hline 7 & 2.718282 & 2.718282 & -9.79e-08\\ \hline \end{tabular}$$

The method of undetermined coefficients can also be used to approximate higher-order derivatives of the given function at any point $x$, based on a set of $m+n+1$ equally spaced points around $x$, with $m$ on its left and $n$ on its right:

$$\frac{d^j}{dx^j}f(x)=f^{(j)}(x)\approx \sum_{k=-m}^n c_k f(x+kh), \;\;\;\;\;\;\; (j=1,\cdots,m+n)$$ (91)

To determine the unknown coefficients $c_i$, we replace each term $f(x+kh)$ by the first $m+n+1$ terms of its Taylor expansion:

$$f^{(j)}(x)\approx \sum_{k=-m}^n c_k f(x+kh)$$

$$\approx \sum_{k=-m}^n c_k \left[f(x)+f'(x)(kh) +\frac{f''(x)}{2}(kh)^2+\cdots+\frac{f^{(m+n)}(x)}{(m+n)!}(kh)^{m+n}\right]$$

$$(j=1,\cdots,m+n)$$ (92)

Collecting all terms for each $f^{(i)}(x),\;(i=0,\cdots,m+n)$ on the right-hand side and equating them to the corresponding terms on the left-hand side (only one non-zero term $f^{(j)}(x)$), we obtain a system of $m+n+1$ linear equations:

$$\sum_{k=-m}^n c_k\frac{(kh)^i}{i!}=\delta_{ij} \;\;\;\; \;\;\;\; \sum_{k=-m}^n k^i c_k =\left\{\begin{array}{ll}i!/h^i&i=j \\ 0&i\ne j\end{array}\right. \;\;\;\;\;\;\;\;\;\;(i=0,\cdots,m+n)$$ (93)

Here $k^i$ can be considered as the component in the ith row and kth column of a matrix ${\bf A}$, and the $m+n+1$ equations above can be written in matrix form:

$$\left[\begin{array}{lllll} (-m)^0&(-m+1)^0&\cdots&(n-1)^0&n^0\\ ... ...ight] =\left[\begin{array}{c}0\\ \vdots\\ j!/h^j\\ \vdots\\ 0\end{array}\right]$$ (94)

or

$${\bf A}{\bf c}_j={\bf b}_j$$ (95)

Here ${\bf b}_j$ on the right-hand side is a vector containing all zero elements except the jth one $b_j=j!/h^j$.

Moreover, we can get ${\bf c}_j$ for the jth derivative $f^{(j)}(x)$ $(j=1,\cdots,m+n$) by solving $m+n$ equations at the same time:

$${\bf A}[{\bf c}_1,\cdots,{\bf c}_{m+n}] =[{\bf b}_1,\cdots,{\bf b}_{m+n}]$$ (96)

to get

$$[{\bf c}_1,\cdots,{\bf c}_{m+n}]={\bf A}^{-1}[{\bf b}_1,\cdots,{\bf b}_{m+n}]$$ (97)

Note that given $m+n+1$ points, the highest derivative that can be approximated is $f^{(m+n)}(x)$. More points are needed to approximate $f^{(j)}(x)$ if $j>m+n$.

Consider the following examples:

• $m=n=1$

  Based on the $m+n+1=3$ points $f(x-h),\;f(x),\;f(x+h)$, we can find the coefficients needed for approximating $f'(x)$ and $f''(x)$ by solving the following equation systems:

  $$\left[\begin{array}{rcc}1 & 1 & 1\\ -1 & 0 & 1\\ 1 & 0 & 1\end{ar... ..._1,\,{\bf c}_2] =\left[\begin{array}{cc}0&0\\ 1/h&0\\ 0&2/h^2\end{array}\right]$$ (98)

  Solving the system we get ${\bf c}_1=[-1/2h,\;0\;1/2h]^T$ and ${\bf c}_2=[1/h^2,\;-2/h^2\;1/h^2]^T$, and
  

  $$f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}$$

  $$f''(x) \approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}$$

  
  We can also use three points all on either side of $x$. For example, if $m=0,\;n=2$, we have:

  $$f'(x) \approx c_0f(x)+c_1f(x+h)+c_2f(x+2h) \\$$

  Solving the equation systems

  $$\left[\begin{array}{rcc}1&1&1\\ 0&1&2\\ 0&1&4\end{array}\right] [... ..._1,\,{\bf c}_2] =\left[\begin{array}{cc}0&0\\ 1/h&0\\ 0&2/h^2\end{array}\right]$$ (99)

  we get ${\bf c}_1=[-3/2h,\;2/h,\;-1/2h]^T$ and ${\bf c}_2=[1/h^2,\;-2/h^2,\;1/h^2]^T$ for approximating $f'(x)$ and $f''(x)$, respectively:
  

  $$f'(x) \approx \frac{-3f(x)+4f(x+h)-f(x+2h)}{2h}$$

  $$f''(x) \approx \frac{f(x)-2f(x+h)+f(x+2h)}{h^2}$$

  
  

• $m=n=2$

  Based on $m+n+1=5$ equally spaced points $f(x+ih),\;(i=-3,\cdots,3)$, we can approximate any of $f'(x),\,f''(x),\,f'''(x),\,f^{(4)}(x)$ by

  $$f^{(j)}(x) \approx c_{-2}f(x-2h)+c_{-1}f(x-h)+c_0f(x)+c_1f(x+h)+c_2f(x+2h)$$ (100)

  The coefficients can be found by solving the following equation systems

  $$\left[\begin{array}{rrrrr}1&1&1&1&1\\ -2&-1&0&1&2\\ 4&1&0&1&4\\ ... ...&0&0\\ 1/h&0&0&0\\ 0&2/h^2&0&0\\ 0&0&6/h^3&0\\ 0&0&0&24/h^4\end{array}\right]$$ (101)

  to get

  $$[{\bf c}_1\;{\bf c}_1\;{\bf c}_1\;{\bf c}_1] =\left[\begin{array}... ...3h^2 & -1/h^3 & -4/h^4 \\ -1/12h &-1/12h^2 & 1/2h^3 & 1/h^4 \end{array}\right]$$ (102)

  and
  

  $$f'(x) \approx \frac{1}{12h}[f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)]$$

  $$f''(x) \approx \frac{1}{12h^2}[-f(x-2h)+16f(x-h)-30f(x)+16f(x+h)-f(x+2h)]$$

  $$f'''(x) \approx \frac{1}{2h^3}[-f(x-2h)+2f(x-h)-2f(x+h)+f(x+2h)]$$

  $$f^{(4)}(x) \approx \frac{1}{h^4}[f(x-2h)-4f(x-h)+6f(x)-4f(x+h)+f(x+2h)]$$ (103)

  
  

• $m=n=3$

  Based on $m+n+1=7$ equally space points $f(x+ih),\;(i=-3,\cdots,3)$, we can approximate and of $f'(x),\,f''(x),\cdots,f^{(6)}(x)$ by

  $$f^{(j)}(x)=\sum_{k=-3}^3 c_k f(x+kh)$$ (104)

  Solving this equation

  $$\left[\begin{array}{rrrrrrr} 1&1&1&1&1&1&1\\ -3&-2&-1&0&1&2&3\\ 9... ...\\ 0&0&0&24/h^4&0&0\\ 0&0&0&0&120/h^5&0\\ 0&0&0&0&0&720/h^6 \end{array}\right]$$

  we get the coefficients in the approximation:

  $$[{\bf c}_1\,{\bf c}_2\,{\bf c}_3\,{\bf c}_4\,{\bf c}_5\,{\bf c}_6... .../h^6 \\ 1/60h & 1/90h^2 &-1/8h^3 & -1/6h^4 & 1/2h^5 & 1/h^6 \end{array}\right]$$

  For example, the first two derivative are:
  

  $$f'(x) \approx \frac{1}{60h}[-f(x-3h)+9f(x-2h)-45f(x-h)$$

  $$+ 45f(x+h)-9f(x+2h)+f(x+4h)]$$

  $$f''(x) \approx \frac{1}{180h}[2f(x-3h)-27f(x-2h)+270f(x-h) -490f(x)$$

  $$+ 270f(x+h)-27f(x+2h)+2f(x+3h)]$$

  
  Alternatively, We can also choose to use seven unequally spaced points, such as $f(x-4h),\;f(x-2h),\;f(x-h),\;f(x), \;f(x+h),\;f(x+2h),\;f(x+4h)$ and get:
  

  $$f'(x) \approx \frac{1}{360h}[-f(x-4h)+40f(x-2h)-256f(x-h)$$

  $$+ 256f(x+h)-40f(x+2h)+f(x+4h)]$$

  $$f''(x) \approx \frac{1}{720h}[f(x-4h)-80f(x-2h)+1024f(x-h) -1890f(x)$$

  $$+ 1024f(x+h)-80f(x+2h)+f(x+4h)]$$

  
  

The Matlab code that implements the method of undetermined coefficients is shown below.

    p=[-m:n];        % all m+n+1 points 
    k=m+n+1;         % total number of points
    A=zeros(k);      % matrix on left-hand side
    b=zeros(k,k-1);  % vector on right-hand side
    A(1,:)=1;        % first row of A
    for i=:k-1         
        A(i+1,:)=p.^i;   
        b(i+1,i)=factorial(i)/h^i;   
    end
    c=inv(A)*b;      % solving linear system to get coefficients
    dx=f(x+p*h)*c;   % f(x): vector of all m+n+1 function values
                     % dx: vector of derivatives of order 1 to m+n

Example: Use the method of undetermined coefficients to find the derivatives of the Gaussian function $f(x)=e^x$ at $x=1$ with 3, 5, 7, and 9 points and step sizes 1, 1/2, 1/4, 1/8. As can be seen from the approximated $f'(x)$ its relative error listed in the table, higher accuracy is achieved by smaller step sizes and larger number of points used in the approximation:

$$\begin{tabular}{c\vert\vert c\vert c\vert c\vert c}\hline h & \mb... ...05&2.7182597/8.2e-08&2.7182819/2.7e-10&2.7182818/9.5e-13\\ \hline \end{tabular}$$