User:EML4500.f08.JAMAMA/NM/Xia

= Numerical Methods =

Extra Solution
From Taylor series :

$$ e^x = 1 + x + \frac{x^2}{2} + \frac{x^3}{6} + \frac{x^4}{24} + \frac{x^5}{120} + ... $$

Then the original fuction could be expressed as:

$$ f(x) = 1 + \frac{x}{2} + \frac{x^2}{6} + \frac{x^3}{24} + \frac{x^4}{120} + ... $$

$$ \Rightarrow \lim_{x\to 0} f(x) = 1 + \frac{0}{2} + \frac{0^2}{6} + \frac{0^3}{24} + \frac{0^4}{120} + ... =1 $$

The original function is plotted from 0 to 1. For comparison, serie-form of f(x) of order 2nd and 3rd are also shown in the following plot.



Extra Solution
From Taylor series :


 * {| style="width:100%" border="0"

$$ e^x = 1 + \sum_{i=1}^n \frac{x^i}{i!} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }

Then the original fuction could be expressed as:


 * {| style="width:100%" border="0"

$$ f(x) = \frac{e^{x}-1}{x} = \sum_{i=1}^n \frac{x^{i-1}}{i!} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle
 * }


 * {| style="width:100%" border="0"

$$ \Rightarrow P_n(x) = \sum_{i=0}^{n} \frac{x^{i}}{(i+1)!} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 2.12)
 * }

To find $$ R_n(x) $$, we need to solve for $$ f^{(n+1)}(x) $$ which is difficult. So, here we compute $$ R_n(x) $$ from the previous function $$ e^{x} $$. From (Eq. 2.7), we could write $$ e^{x} $$ as:


 * {| style="width:100%" border="0"

$$ e^{x} = 1+ \sum_{i=1}^{n+1} \frac{x^i}{i!} + \frac{x^{n+2}}{(n+2)!} e^\xi, \quad \xi \in [0, x] $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle $$
 * }


 * {| style="width:100%" border="0"

$$ \Rightarrow f(x) = \frac{e^{x}-1}{x} = \sum_{i=0}^{n} \frac{x^i}{(i+1)!} + \frac{x^{n+1}}{(n+2)!} e^\xi, \quad \xi \in [0, x] $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 2.13)

$$
 * }

Comparing (Eq. 2.12) and (Eq. 2.13) yields,


 * {| style="width:100%" border="0"

$$ R_n(x) = \frac{x^{n+1}}{(n+2)!} e^\xi, \quad \xi \in [0, x] $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle

$$
 * }

Solution
Noting from equation 3.10


 * {| style="width:100%" border="0"

$$ f_n(x) = {P}_{n}(x) = \sum_{i=0}^{n} \frac{x^{i}}{(i+1)!} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 4.1)
 * }

So the integration could be derived as


 * {| style="width:100%" border="0"

$$ I_n = \int_{-1}^{1} f_n(x)dx = \int_{-1}^{1} \sum_{i=0}^{n} \frac{x^{i}}{(i+1)!}dx = \sum_{i=0}^{n} \int_{-1}^{1}\frac{x^{i}}{(i+1)!}dx $$ $$ substitute $$ \quad n $$ with $$ \quad 2k $$ yields,
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle (Eq. 4.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I_n = \sum_{i=0}^{k} \frac{2}{(2i+1)(2i+1)!} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.3)
 * }

As $$ \quad n = 2,4,8,16,... $$, use matlab to solve for the integration and the code is attached below:

The results are shown in the following table:

Solution
From Simpson's rule :

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I_n(x) = \frac {h}{3} [f_0 + 4f_1 + 2\sum_{i=2}^{n-2} f_i + 4f_{n-1} + f_n] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Matlab code is written to solve the problem and attached below:

The solution and iteration history is shown in following table:

Solution
For Gauss-Legendre quad, the integration could be written as :

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I_n(x) = \sum_{i=1}^{n} w_if(x_i) $$ for $$ i=1,2,...,n $$ and $$ x_i $$ are the roots for (Eq) $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

The weight functions are given by :

<span id="(1)">
 * {| style="width:100%" border="0"

$$ w_i = \frac {-2}{(n+1)P_{n}^{(1)}(x_i)P_{n+1}(x_i)} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ P_{n}(x) = \sum_{i=0}^{[n/2]} (-1)^i\frac{(2n-2i)!x^{n-2i}}{2^ni!(n-i)!(n-2i)!} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Combining the above equations, we can obtain the roots and weight function values for n=2 and n=4 as presented in lecture notes. The following matlab code uses those values directly and compute the integrations for n=2 and n=4.

The results are shown in the following table:

Obviously from the above table, the error for n=4 already smaller than the order of 0.000001, so there is no need of further computaions.

Problem
Given from lecture notes:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ p_2(x)=\sum_{i=0}^{n=2} l_{i,2}(x)f(x_i) $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.1)
 * }
 * {| style="width:100%" border="0"

$$ l_{i,n}=\prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.2)
 * }

Derive: Simple simpson's rule.

Solution
From definition:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I_2(f) = \int_{a}^{b} f_2(x)dx = \sum_{i=0}^{2} \int_{a}^{b} l_{i,2}(x)dx f(x_i) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.3)
 * }

With

<span id="(1)">
 * {| style="width:100%" border="0"

$$ x_0=a;x_1=\frac {a+b}{2};x_2=b. $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

from (Eq. 11.2), we can obtain:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{0,2}(x) = \frac {(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} = \frac {2}{(a-b)^2} [x^2-\frac{a+3b}{2} x +\frac{(a+b)b}{2}] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.4)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{1,2}(x) = \frac {(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} = -\frac {4}{(a-b)^2} [x^2-(a+b) x +ab] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.5)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{2,2}(x) = \frac {(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} = \frac {2}{(b-a)^2} [x^2-\frac{3a+b}{2} x +\frac{(a+b)a}{2}] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.6)
 * }

Integrate (Eq. 11.4)(Eq. 11.5)(Eq. 11.6), we can solve

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \int_{a}^{b}l_{0,2}(x) = \frac {2}{(b-a)^2} [-\frac{1}{12} (a-b)^3] = \frac{1}{6}(b-a) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.7)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \int_{a}^{b}l_{1,2}(x) = \frac {4}{(a-b)^2} [-\frac{1}{6} (a-b)^3] = \frac{4}{6}(b-a) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.8)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \int_{a}^{b}l_{2,2}(x) = \frac {2}{(b-a)^2} [-\frac{1}{12} (a-b)^3] = \frac{1}{6}(b-a) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.9)
 * }

With

<span id="(1)">
 * {| style="width:100%" border="0"

$$ h=\frac {b-a}{2} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Plug (Eq. 11.7)(Eq. 11.8)(Eq. 11.9) into (Eq. 11.3), yields

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I_2(f) = \int_{a}^{b} l_{0,2}(x)dx f(x_0) + \int_{a}^{b} l_{1,2}(x)dx f(x_1) + \int_{a}^{b} l_{2,2}(x)dx f(x_2) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 11.10)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow I_2(f) = \frac {h}{3} (f(x_0) + 4f(x_1) + f(x_2)) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Which is right the simple simpson's rule.

Problem
Given:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f(x)=\frac {e^x-1}{x} on [-1,1], x_0=-1,x_n=1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.1)
 * }

Find: 1) Construct $$ f_n(x) $$ by lagrange interpolation functions for $$ n=1,2,4,8,16 $$ and then plot $$ f(x) $$ and $$ f_n(x) $$. 2) Compute weight functions for $$ n=1,2,4,8,16 $$ and then form tables. 3) Compute $$ I_n(x) $$ for $$ n=1,2,4,8 $$ and compare to $$ I $$ that is solved using WolframAlpha with more digits. 4) For $$ n=5 $$, plot $$ l_0,l_1,l_2 $$ and then predict the shapes of $$ l_3,l_4,l_5 $$ based on those.

Solution
1) According to the problem, we construct the fuctions from following equations:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_n(x)=p_n(x)=\sum_{i=0}^{n} l_{i,n}(x)f(x_i) $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.2)
 * }
 * {| style="width:100%" border="0"

$$ l_{i,n}=\prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.3)
 * }

For n=1,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ x_0=-1,x_1=1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{0,1}(x) = \frac {x-x_1}{x_0-x_1} = -\frac{1}{2}(x-1) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.4)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{1,1}(x) = \frac {x-x_0}{x_1-x_0} = \frac{1}{2}(x+1) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.5)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f(x_0) = \frac {e^{x_0}-1}{x_0} = 1-e^{-1} \displaystyle $$ $$ f(x_1) = \frac {e^{x_1}-1}{x_1} = e-1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow f_1(x) = l_{0,1}(x)f(x_0) + l_{1,1}(x)f(x_1) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.6)
 * }

For n=2,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ x_0=-1,x_1=0,x_2=1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{0,2}(x) = \frac {(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} = \frac{1}{2}x(x-1) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.7)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{1,2}(x) = \frac {(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} = -(x+1)(x-1) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.8)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{1,2}(x) = \frac {(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} = \frac{1}{2}(x+1)x $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.9)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f(x_0) = \frac {e^{x_0}-1}{x_0} = 1-e^{-1} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f(x_1) = \frac {e^{x_1}-1}{x_1} = 1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f(x_2) = \frac {e^{x_2}-1}{x_2} = e-1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Rightarrow f_2(x) = l_{0,2}(x)f(x_0) + l_{1,2}(x)f(x_1) +  l_{2,2}(x)f(x_2)$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.10)
 * }

For n=4,8,16, we use matlab for plotting with (Eq. 12.2) and (Eq. 12.3). The code is attached below:

The plots are shown as following:



2) For n=1, integrate previous (Eq. 12.4) and (Eq. 12.5) from [-1,1], we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ w_{0,1} = \int_{-1}^{1} l_{0,1}(x)dx = 1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ w_{1,1} = \int_{-1}^{1} l_{1,1}(x)dx = 1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

For n=2, integrate previous (Eq. 12.7), (Eq. 12.8) and (Eq. 12.9) from [-1,1], we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ w_{0,2} = \int_{-1}^{1} l_{0,2}(x)dx = \frac{1}{3} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ w_{1,2} = \int_{-1}^{1} l_{1,2}(x)dx = \frac{4}{3} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ w_{2,2} = \int_{-1}^{1} l_{2,2}(x)dx = \frac{1}{3} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

For n=4,8,16, by apply LIET, we obtain the linear system equations:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \sum_{i=0}^{n}w_{i,n}(x_i)^j = \frac{b^{j+1}-a^{j+1}}{j+1} $$ when $$ j=0,1,2,...,n $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.11)
 * }

Solving (Eq. 12.11) for n=4,8,16 by matlab and the code is attached below:

The results are shown in the following tables:

3) Integrate using Newton-Cotes Method, The integration could be expressed as:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I_n(f) = \int_{a}^{b} f_n(x)dx = \sum_{i=0}^{n} w_{i,n}f(x_i) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.12)
 * }

With previously calculated value of weight functions, The results of the integrations are shown in the following table and compared with that obtained from WolframAlpha.

4) For n=5, the lagrange functions are computed from (Eq. 12.3) using matlab, the code is attached below:

The first three fuctions are plotted as follows:



Now, we will analyze the shapes of the lagrange functions. First it is obvious to identify:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ x_i = -x_{n-i} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

From (Eq. 12.3), we could derive

<span id="(1)">
 * {| style="width:100%" border="0"

$$ l_{n-i,n}(x) = \prod_{j=0,j \neq n-i}^{n} \frac{x-x_j}{x_{n-i}-x_j} = \prod_{j=0,j \neq n-i}^{n} \frac{x+x_{n-j}}{x_{n-i}+x_{n-j}} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \prod_{j=0,n-j \neq i}^{n} \frac{-x-x_{n-j}}{x_{i}-x_{n-j}} = \prod_{k=0,k \neq i}^{n} \frac{-x-x_{k}}{x_{i}-x_{k}} = l_{i,n}(-x) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 12.13)
 * }

So from (Eq. 12.13), we can predict that $$ l_3 $$ and $$ l_2 $$, $$ l_4 $$ and $$ l_1 $$, $$ l_5 $$ and $$ l_0 $$ are axis-symmetric function pairs (With respect to y axis). We further plot $$ l_3 $$, $$ l_4 $$, $$ l_5 $$ in the following figure and the prediction is verified.



= 3.7 Integration Derivation for Error of Trap. rule =

Prove
Prove the following equation

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{M_2}{2!} \int_{a}^{b} \left | (x-a)(x-b) \right | dx = \frac {(b-a)^3}{12} M_2 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 3.7.1)
 * }

Solution
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{M_2}{2!} \int_{a}^{b} \left | (x-a)(x-b) \right | dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_2}{2!} \int_{a}^{b} (x-a)(b-x) dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_2}{2} \int_{a}^{b} [-x^2+(a+b)x-ab] dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_2}{2} [-\frac{1}{3}x^3+\frac{1}{2}(a+b)x^2-abx]\mid _{a}^{b}  $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_2}{2} [-\frac{1}{3}(b^3-a^3)+\frac{1}{2}(a+b)(b^2-a^2)-ab(b-a)] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_2}{12} (b^3-3b^2a+3a^2b-a^3) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_2}{12} (b-a)^3 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Prove complete. Egm6341.s11.team3.Xia 19:00, 11 February 2011 (UTC)

= 3.8 Integration Derivation for Error of Simple Simpson's rule =

Prove
Prove the following equation

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{M_3}{3!} \int_{a}^{b} \left | (x-a)(x-\frac{a+b}{2})(x-b) \right | dx = \frac {(b-a)^4}{192} M_3 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 3.8.1)
 * }

Solution
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{M_3}{3!} \int_{a}^{b} \left | (x-a)(x-\frac{a+b}{2})(x-b) \right | dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_3}{3!} [\int_{a}^{\frac{a+b}{2}} (x-a)(x-\frac{a+b}{2})(x-b) dx + \int_{\frac{a+b}{2}}^{b} (x-a)(x-\frac{a+b}{2})(b-x) dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_3}{3!} (\int_{a}^{\frac{a+b}{2}} - \int_{\frac{a+b}{2}}^{b}) [x^3-\frac{3}{2}(a+b)x^2+\frac{a^2+4ab+b^2}{2}x- \frac{a^2b+b^2a}{2}] dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_3}{3!} [\frac{1}{4}x^4-\frac{1}{2}(a+b)x^3+\frac{a^2+4ab+b^2}{4}x^2- \frac{a^2b+b^2a}{2}x] (\mid_{a}^{\frac{a+b}{2}} + \mid_{b}^{\frac{a+b}{2}}) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{2M_3}{6} [\frac{1}{64}(a+b)^4-\frac{1}{16}(a+b)(a+b)^3+\frac{a^2+4ab+b^2}{16}(a+b)^2- \frac{a^2b+b^2a}{4}(a+b)] $$ $$ - \frac{M_3}{6} [\frac{1}{4}a^4-\frac{1}{2}(a+b)a^3+\frac{a^2+4ab+b^2}{4}a^2- \frac{a^3b+b^2a^2}{2}] $$ $$ - \frac{M_3}{6} [\frac{1}{4}b^4-\frac{1}{2}(a+b)b^3+\frac{a^2+4ab+b^2}{4}b^2- \frac{a^2b^2+b^3a}{2}] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_3}{192} (a^4-4a^3b-10a^2b^2-4ab^3+b^4) - \frac{M_3}{6} (-\frac{a^2b^2}{2}) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_3}{192} (a^4-4a^3b+6a^2b^2-4ab^3+b^4) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ = \frac{M_3}{192} (b-a)^4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Prove complete. Egm6341.s11.team3.Xia 20:00, 11 February 2011 (UTC)

= 4.6 Fix Runge Phenomenon Use Non-uniform Distribution of Nodes =

Given:
Given a function

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f(x) = \frac{1}{1+4x^2}\quad $$ in interval $$ x \in [-5,5] \quad $$, $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 3.3.1)
 * }

Find:
1) Use the roots of legendre polynomials to compute $$\quad f_n^L(x) $$, plot $$\quad f_n^L(x) $$ and $$\quad f(x) $$, increase n till $$\quad I_n $$ is accurate to $$\quad O(10^{-6}) $$.

2) Use the roots of chebyshew polynomials to compute $$\quad f_n^L(x) $$, plot $$\quad f_n^L(x) $$ and $$\quad f(x) $$, Use Newton-Cotes to find $$\quad I_n $$ to $$\quad O(10^{-6}) $$ accuracy.

3) Use Gauss-legendre quadrature to find $$\quad I_n $$ to $$\quad O(10^{-6}) $$ accuracy.

4) Compare results of 1), 2) and 3).

Solution:
For comparison, the integration result is computed from WolframAlpha:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \int_{-5}^{5} \frac{1}{1+4x^2} dx = 1.471127674303734 \quad $$, $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.1)
 * }

1) To apply legendre polynomials, we convert the variable domain from [-5 5] to [-1 1] by the following transformation:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f^*(X) = \frac{1}{1+100X^2}\quad $$ and $$ X=x/5 \in [-1,1] \quad $$, $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.1)
 * }

And the integration realtions are given by:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I(x) = \int_{-5}^{5} f(x)dx = \int_{-1}^{1} f^*(X)d(5X) = 5I^*(X) $$, $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.2)
 * }

The legendre polynomial are given by: <span id="(1)">
 * {| style="width:100%" border="0"

$$ P_{n}(X) = \sum_{i=0}^{[n/2]} (-1)^i\frac{(2n-2i)!X^{n-2i}}{2^ni!(n-i)!(n-2i)!} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.3)
 * }

From Newton-Cotes method:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f_n^{*L}(X)=\sum_{i=0}^{n+1} l_{i,n+1}(X)f^*(X_i) $$ $$ <span id="(1)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.4)
 * }
 * {| style="width:100%" border="0"

$$ l_{i,n+1}=\prod_{j=0,j \neq i}^{n+1} \frac{X-X_j}{X_i-X_j} $$ $$ where $$ {X_i, i=1,2,...,n} \quad $$ are the roots of (Eq. 4.6.1), $$ X_0=-1, X_{n+1}=1\quad $$. $$ I_n(X) \quad $$ is given by:
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.5)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I_{n+1}^*(X)=\sum_{i=0}^{n+1} w_{i,n+1} f^*(X_i) $$ $$ the $$ w_{i,n+1} \quad $$ is estimated by solving the linear system matrix:
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.6)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \sum_{i=0}^{n+1} w_{i,n+1}(X_i)^j = \frac {b^{j+1}-a^{j+1}}{j+1} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.7)
 * }

With matlab solving the roots for (Eq. 4.6.1), we can obtain $$\quad f_n^L(x) $$ and $$\quad I_n $$. The code is attached below:

The resultant plots are shown below:



The Integrations and Errors are listed in the table below:

2) To apply chebyshew polynomials, we convert the variable domain from [-5 5] to [-1 1] by the same transformation as (Eq. 4.6.1) and (Eq. 4.6.2):

The roots of chebyshew polynomials are given by: <span id="(1)">
 * {| style="width:100%" border="0"

$$ X_{i,n} = cos(i\frac{\pi}{n})$$ and $$ i=0,1,2,...,n \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.8)
 * }

Then similar to the previous section, we apply Newton-Cotes method and lagrange interpolation method as in (Eq. 4.6.4) and (Eq. 4.6.6). Using matlab to obtain $$\quad f_n^L(x) $$ and $$\quad I_n $$. The code is attached below:

The resultant plots are shown below:



The Integrations and Errors are listed in the table below:

3) To apply Gauss-legendre quadrature, we convert the variable domain from [-5 5] to [-1 1] by the same transformation as (Eq. 4.6.1) and (Eq. 4.6.2).

The legendre polynomial are given by (Eq. 4.6.3), from which the roots are obtained. The weight functions are computed by:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ w_i= \frac {-2}{(n+1)P_n^{(1)}(X_i)P_{n+1}(X_i)} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.9)
 * }

The integration $$\quad I_n^* $$ is therefore given by:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ I_n^* = \sum_{i=1}^{n} w_i f^*(X_i) $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 4.6.10)
 * }

Using matlab to solve for the integrations. The code is attached below:

The Integrations and Errors are listed in the table below:

4) Explanation for not achieving $$ O(10^{-6}) \quad $$: Basically, in this assignment, n is computed to 32. For higher n(such as 64), the matlab "roots" function would fail and generated complex values for legendre poly roots. Moreover, the matrix for computing the weights also became "ill" so that it is considered "single" by the machine. Therefore, n=32 is the best accuracy we got for this problem.

Comparing the results of the error, it is obvious that Newton-Cotes method (first 2 cases) has a better accuracy than Gauss-Legendre Quadrature method when n is small, while the Gauss-Legendre Quadrature method converges faster than the Newton-Cotes method. Comparing the different interpolation approaches for the Newton-Cotes method, they have similar behaviors of accuray. When n=32, all the three methods are of the order of $$ O(10^{-3}) \quad $$.

Egm6341.s11.team3.Xia 5:00, 1 March 2011 (UTC)

Results:
The integration results for a), c) and d) are shown and compared in the following table:

The Romberg table for b) is:

(2) Circumference for Ellipse: Method1
From definition of eccentricity of ellipse :

<span id="(1)">
 * {| style="width:100%" border="0"

$$ e = \sqrt {1-\frac{b^2}{a^2}} = \frac{\sqrt{91}}{10} \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.11)
 * }

From (Eq. 6.7.2), (Eq. 6.7.3), we derive the Integration as:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ C = I = \int_{0}^{2 \pi} a(1-e^2)(1-ecos \theta)^{-2}\sqrt{(1 +e^2 - 2e cos \theta)} d \theta \quad $$ $$ = \int_{0}^{2 \pi} \frac{9}{10} (1- \frac{\sqrt{91}}{10} cos \theta)^{-2}(\frac{191}{100} - \frac{\sqrt{91}}{5} cos \theta)^{\frac{1}{2}} d \theta \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.12)
 * }

Using WolframAlpha, the exact solution is calculated to be :

<span id="(1)">
 * {| style="width:100%" border="0"

$$ C_{exact} = I_{exact} = 43.8591 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.13)
 * }

Following the same approaches as in (Eq. 6.7.7), (Eq. 6.7.8) and (Eq. 6.7.9), the integration is computed using the following code:

Results:
The integration results for a), c) and d) are shown and compared in the following table:

The Romberg table for b) is:

(3) Circumference for Ellipse: Method2
The 2nd method of circumference for ellipse follows (Eq. 6.7.4). Similar to previous calculations, the function integrated here is defined as:

<span id="(1)">
 * {| style="width:100%" border="0"

$$ f(\theta) = 4a[1-e^2 sin^2 \theta]^{\frac{1}{2}} = 40\sqrt{1- 0.91 sin^2 \theta} \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.14)
 * }

Using WolframAlpha, the exact solution is calculated to be :

<span id="(1)">
 * {| style="width:100%" border="0"

$$ C_{exact} = I_{exact} = 43.8591 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 6.7.13)
 * }

Which is identical to (Eq. 6.7.13), following the same approaches as in (Eq. 6.7.7), (Eq. 6.7.8) and (Eq. 6.7.9), the integration is computed using the following code:

Results:
The integration results for a), c) and d) are shown and compared in the following table:

The Romberg table for b) is:

Egm6341.s11.team3.Xia 18:00, 31 March 2011 (UTC)

=References=