fixed typo

0ecdfbc6 · Isabel Slingerland · 625d3bf8 · 0ecdfbc6
Commit 0ecdfbc6 authored 6 months ago by Isabel Slingerland
--- a/book/numerical_methods/4-numerical-integration.ipynb
+++ b/book/numerical_methods/4-numerical-integration.ipynb
@@ -338,7 +338,7 @@
    "$$\n",
    "\n",
    "\n",
-    "```{admonition} ** Derivation**\n",
+    "```{admonition} **Derivation**\n",
    ":class: tip, dropdown\n",
    "\n",
    "Start by using a Taylor series expansion:\n",

 %% Cell type:markdown id: tags:

 # Numerical Integration

 %% Cell type:markdown id: tags:

 NOTE
 ```{note}
 **Important things to retain from this block:**
 * Knowing how to apply different discrete integration techniques
 * Having an idea on how the truncation errors change from technique to technique
 ```

 %% Cell type:markdown id: tags:

 ## Motivation
 Numerical integration is used often in diverse fields of engineering. For example to develop advanced numerical model techniques or to analyze field or laboratory data where a desired physical quantity is an integral of other measured quantitites. So, the analytic integral is not always known and has to be computed (approximated) numerically. For example, the variance density spectrum (in ocean waves) is a discrete signal obtained from the water level variations and its **integral** is a quantity that characterizes the wave field.

 %% Cell type:markdown id: tags:

 ## Discrete Integration

 %% Cell type:markdown id: tags:




 In a 1 dimension situation, the integral

 $$
 I=\int_a^b f(x)dx
 $$

 is defined between points and represents the area under the curve $f$ limited by $[a,b]$. A visual way to approximate that area would be to use a squared page and count the number of squares that fall under that curve. In the example figure, you can count 11 squares.

 **The numerical integration is an approximation** to that area done in a strict way. There are several paths to follow  using polynomial representation. These widely used formulas are known as Newton-Cotes formulas.


 * First, we split the interval of integration (a,b) on several partial intervals.

 * Then, we replace $f(x)$ by the interpolation polynomial in each interval and calculate the integral over it.

 * Finally, sum all the areas to find (approximate) the integral.


 ```{figure} figs/integral_illustration.png
 :name: integral_illustration

 two integral illustration showing the area under the curve for interval [a,b]
 ```

 %% Cell type:markdown id: tags:

 ### Left Riemann Integration

 ```{figure} figs/left_riemann_new.png
 :name: left_riemann_new.png

 Illustrating the left Riemann Sum over interval [a,b]
 ```

 Using the lowest polynomial order leads to the Riemann integral also known as the rectangle rule. Basically, the function is divided in segments with a known $\Delta x$ and the area of each segment is represented by a rectangle. Its height, related to the discretized points of the function, can be defined by the left or by the right.

 Therefore, we get the left Riemann integration formula given by

 $$
 \int_{a}^{b} f(x) \, dx \approx  \sum_{i=0}^{n-1} f(x_i)\Delta x, \hspace{5mm} \Delta x = \frac{b - a}{n}
 $$

 where $n$ is the number of segments.

 ### Right Riemann Integration



 ```{figure} figs/right_riemann_new.png
 :name: right_riemann_new.png

 Illustrating the right Riemann Sum over interval [a,b]
 ```

 In an analogous way, the formula for the right Riemann integration is:

 $$
 \int_{a}^{b} f(x) \, dx \approx  \sum_{i=1}^{n} f(x_i)\Delta x \approx \sum_{i=0}^{n-1} f(x_{i+1})\Delta x
 $$

 Note the difference respect to the left Riemann integration, specially the starting and end points.

 ### Midpoint Rule


 ```{figure} figs/midpoint_new.png
 :name: midpoint_new.png

 Illustrating the right Midpoint Rule over interval [a,b]
 ```
 The midpoint rule defines the height of the rectangle at the midpoint between the two end points of each interval:

 $$
 \int_{a}^{b} f(x) \, dx \approx \sum_{i=0}^{n-1} f(\frac{x_{i}+x_{i+1}}{2}) \Delta x
 $$



 ### Trapezoidal Rule


 ```{figure} figs/trapezoidal_new.png
 :name: trapezoidal_new.png

 Illustrating the Trapezoidal Rule over interval [a,b]
 ```

 The trapezoidal rule uses a polynomial of the first order. As a result the area to be computed is a trapezoid (the top of the rectangle has a slope now) where the bases are the values assumed by the function in the end points of each of the partial intervals, and its height is just the difference between two consecutive points. Mathematically it reduces to the average of two rectangles:

 $$
 \int_{a}^{b} f(x) \, dx \approx \sum_{i=0}^{n-1} \frac{f(x_{i}) + f(x_{i+1})}{2} \Delta x
 $$

 ### Simpson’s Rule



 ```{figure} figs/Simpson_s.png
 :name: Simpson_s.png

 Illustrating the Simpson Rule over interval [a,b]
 ```
 Simpson’s rule uses a second order polynomial. It computes the integral of the function $f(x)$ by fitting parabolas needing three consecutive function points and adding the areas under them (see $p_1(x)$ and $p_2(x)$ in the figure). In formula form:

 $$
 \int_{a}^{b} f(x) \, dx \approx \sum_{i=1}^{n/2} \frac{f(x_{2i-2}) + 4f(x_{2i-1}) + f(x_{2i})}{6} 2\Delta x
 $$

 Note that the interval length is multiplied with two, since three points are needed.




 %% Cell type:markdown id: tags:

 :::{card} ***Exercise***

 Estimate the integral $\int_{0.5}^1 x^4$ with a total number of steps n=2:
 - Left Rieman Sum
 - Midpoint rule
 - Trapezoidal

 and then compute the real error $\epsilon= | f(x)- \hat{f(x)|}$, where f(x) is the exact integral and $\hat{f(x)}$ the approximated integral

 `````{admonition} Tip
 :class: tip
 Tip: feel free to use the code cell below as a calculator.
 `````


 ```{admonition} **Solution**
 :class: tip, dropdown


 - find $\Delta x$,

 $$
 \frac{b-a}{n}=\frac{1-0.5}{2}=0.25
 $$

 **Left Rieman Sum**:

 $$
 \sum_{i=0}^{n-1} f(x_i)\Delta x = f(x_0) \Delta x + f(x_1) \Delta x
 $$


 $$
 0.5^4 \cdot 0.25 + 0.75^4 \cdot 0.25 =  0.0947265625
 $$

 **Midpoint Rule**:

 $$
 \sum_{i=0}^{n-1} f(\frac{x_{i}+x_{i+1}}{2}) \Delta x = f(\frac{x_0+x_1}{2}) \Delta x + f(\frac{x_1+x_2}{2}) \Delta x
 $$



 $$
 0.625^4 0.25 + 0.875^4 0.25= 0.1846923828125
 $$

 **Trapezoidal Rule**:

 $$
 \sum_{i=0}^{n-1} \frac{f(x_{i}) + f(x_{i+1})}{2} \Delta x =  \frac{f(x_{0}) + f(x_{1})}{2} \Delta x + \frac{f(x_{1}) + f(x_{2})}{2} \Delta x
 $$

 $$
 \frac{f(0.5) + f(0.75)}{2} \Delta x + \frac{f(0.75) + f(1)}{2} \Delta x = \frac{(0.5^4 +0.75^4)}{2} \cdot 0.25 + \frac{(0.75^4 +1^4)}{2} \cdot 0.25= 0.2119140625
 $$

 **Exact**

 $$
 \int_{0.5}^1 x^4$ = (\frac{1}{5}x^5)\bigg|_{0.5}^{1} = \frac{1}{5}1^5- \frac{1}{5}0.5^5= 0.19375
 $$

 **Error**


 $$
 \epsilon= | f(x)- \hat{f(x)}|
 $$

 error Left Rieman = | 0.19375 - 0.0947265625| = 0.0990234375

 error Midpoint Rule = | 0.19375 - 0.1846923828125| = 0.009057617187500006

 error Trapezoidal Rule = | 0.19375 - 0.2119140625| = 0.018164062499999994

 ```
 :::

 %% Cell type:code id: tags:

 ``` python
 #feel free to use this as your calculator for the exercise above

 def f(x):
    return x**4

 f_left_rieman=
 f_Midpoint=
 f_trapezoidal=
 f_exact=

 print(f"f_left_rieman = {f_left_rieman}")
 print(f"f_Midpoint = {f_Midpoint}")
 print(f"f_trapezoidal = {f_trapezoidal}")
 print(f"f_exact = {f_exact}")
 print(f"Error left Riemann = {abs(f_exact-f_left_rieman)}")
 print(f"Error Midpoint = {abs(f_exact-f_Midpoint)}")
 print(f"Error trapezoidal = {abs(f_exact-f_trapezoidal)}")
 ```

 %% Cell type:markdown id: tags:

 ### Truncation Errors

 The truncation error for each method is the difference between the exact integral and the approximated integral. In the case of the left Riemann method:

 $$
 \text{Left Riemann Error }=  \int_a^b f(x) dx - \sum_{i=0}^{n-1} f(x_i)\Delta x
 $$


 Now, we can use the Taylor Series Expansion of $f(x)$ instead of $f(x)$ itself and solve that integral. If you are curious about the mathematical steps, further below you can find the detailed derivation. For the moment, let's see the absolute errors for the different techniques.

 $$
 \text{Left Riemann Error }=  \left| \bar f'(b-a)\Delta x /2  \right| \text { therefore } \mathcal O(\Delta x)
 $$

 $$
 \text{Right Riemann Error }=  \left| \bar f'(b-a)\Delta x /2  \right| \text { therefore } \mathcal O(\Delta x)
 $$

 $$
 \text{Midpoint Error }=  \left| \bar f''(b-a)\Delta x^2 /2  \right| \text { therefore } \mathcal O(\Delta x^2)
 $$

 $$
 \text{Trapezoidal Error }=  \left| \bar f''(b-a)\Delta x^2 /2  \right| \text { therefore } \mathcal O(\Delta x^2)
 $$


 $$
 \text{Simpsons Error }=  \left| \bar f''''(b-a)\Delta x^4 /2  \right| \text { therefore } \mathcal O(\Delta x^4)
 $$

 #### Analysis

 As we can see, for the left and right Riemann integration techniques we find errors $\sim \dfrac{1}{n}$, while for the midpoint and trapezoidal rules, we find them $\sim\dfrac{1}{n^2}$, and with Simpson's rule $\sim\dfrac{1}{n^4}$.

 %% Cell type:markdown id: tags:

 #### Truncation Error Derivation for Left Riemann integration

 ```{note}

 You don't need to know the derivation for the exam.
 ```






 %% Cell type:markdown id: tags:

 :::{card} Derivation


 $$
 \text{Left Riemann Error }= | \int_a^b f(x) dx - \sum_{i=0}^{n-1} f(x_i)\Delta x | \hspace{5mm} (1)
 $$


-```{admonition} ** Derivation**
+```{admonition} **Derivation**
 :class: tip, dropdown

 Start by using a Taylor series expansion:


 $$
 f(x)= f(x_i)+ (x-x_i)f'(\xi(x)) + \mathcal{O}(\Delta x)^2 \hspace{5mm}(2)
 $$


 Substitute this Taylor expansion (2) in the term $\int_a^b f(x)dx$ in equation (1):


 $$
 \int_a^b \bigg[ f(x_i)+ (x-x_i)f'(x_i)+ \mathcal{O}(\Delta x)^2 \bigg] dx \hspace{5mm} (3)
 $$

 Rewrite this as the sum of the integrals over the sub-intervals $[x_i,x_{i+1}]$:

 $$
 = \sum_{i=0}^{n-1} \int_{x_i}^{x_{x+1}} \bigg[  f(x_i)+ (x-x_i)f'(x_i)+ \frac{1}{2}(x-x_i)^2f''(x)+ \mathcal{O}(\Delta x)^2 \bigg]dx \hspace{5mm} (4)
 $$

 Using integration by substitution, $u= x-x_i$, therefore $\frac{du}{dx}=1$ meaning du=dx

 $$
 =\sum_{i=0}^{n-1} \int_{x_i}^{x_{x+1}} \bigg[  f(x_i)+ (u)f'(x_i) + \mathcal{O}(\Delta x)^2 \bigg]du \hspace{5mm} (5)
 $$


 $$
 =\sum_{i=0}^{n-1}  \bigg[  (u)f(x_i)+ \frac{(u)^2}{2}f'(x_i)+ \mathcal{O}(\Delta x)^2  \bigg] \bigg|^{x_{i+1}}_{x_{i}} \hspace{5mm} (6)
 $$

 $$
 =\sum_{i=0}^{n-1}  \bigg[  (x_{i+1}-x_i)f(x_i)+ \frac{(x_{i+1}-x_i)^2}{2}f'(x_i))+ \mathcal{O}(\Delta x)^2  \bigg] \hspace{5mm} (7)
 $$

 $$
 =\sum_{i=0}^{n-1}  \bigg[  \Delta xf(x_i)+ \frac{\Delta x^2}{2}f'(x_i))+ \mathcal{O}\Delta x^2  \bigg] \hspace{5mm} (8)
 $$

 Plug equation (8) in equation (1):

 $$
 |\sum_{i=0}^{n-1}  \bigg[  (\Delta x)f(x_i)+ \frac{\Delta x^2}{2}f'(x_i))+ \mathcal{O}(\Delta x)^2\bigg]  - \sum_{i=0}^{n-1} f(x_i)\Delta x | \hspace{5mm} (9)
 $$

 $$
 =|\sum_{i=0}^{n-1}  \bigg[  \Delta xf(x_i)+ \frac{\Delta x^2}{2}f'(x_i))+ \mathcal{O}(\Delta x)^2  - f(x_i)\Delta x\bigg]| \hspace{5mm} (10)
 $$

 The term $\Delta xf(x_i)$ falls away:

 $$
 =|\sum_{i=0}^{n-1}  \bigg[  \frac{\Delta x^2}{2}f'(x_i))+ \mathcal{O}(\Delta x)^2 \bigg]| \hspace{5mm} (11)
 $$


 $$
 = n \cdot (b-a)\frac{\Delta x^2}{2}f'(x_i) \hspace{5mm} (12)
 $$

 $\Delta x = \frac{b-a}{n}$ this gives:

 ```

 $$
 \text{Left Riemann Error }= | \int_a^b f(x) dx - \sum_{i=0}^{n-1} f(x_i)\Delta x | = | \bar{f}'(b-a) \frac{\Delta x}{2} |
 $$

 The resulting left Riemann error is proportional to the average of the derivatives, the length of the domain (b-a) and the step size (\Delta x). This result tells us that for steeper functions the error will be higher than for milder ones. For example, a function of a straight horizontal line with average derivatives equal to 0 would have an error equal to 0! Also, the error will be larger if the domain is larger, which is kind of obvious. As those first two terms (average derivative and domain length) are inherent to the problem requirements, we say that the error is of the order \Delta x.


 :::

 %% Cell type:markdown id: tags:

 :::{card} Exercise


 Estimate the integral $\int_{1}^{3} \frac{5}{x^4} dx$ with a step-size $\Delta x= 0.5$

 - Right Riemann Sum
 - Simpson Rule

 `````{admonition} Tip
 :class: tip
 Tip: feel free to use the code cell below as a calculator.
 `````

 ```{admonition} **Solution**
 :class: tip, dropdown

 Find $n$

 $$
 n = \frac{b-a}{\Delta x}= \frac{3-1}{0.5}= 4
 $$

 ***Right Riemann Sum***

 $$
 \int_{a}^{b} f(x) \, dx \approx  \sum_{i=1}^{n} f(x_i)\Delta x \approx \sum_{i=0}^{n-1} f(x_{i+1})\Delta x
 $$

 $$
 =  \frac{5}{1.5^4} \cdot 0.5 +\frac{5}{2^4} \cdot 0.5 + \frac{5}{2.5^4} \cdot 0.5 + \frac{5}{3^4} \cdot 0.5
 $$

 $$
 = 0.7449413580246914
 $$

 ***Simpsons Rule***

 $$
 \int_{a}^{b} f(x) \, dx \approx \sum_{i=1}^{n/2} \frac{f(x_{2i-2}) + 4f(x_{2i-1}) + f(x_{2i})}{6} 2\Delta x
 $$

 $$
 \frac{ \frac{5}{1^4}+ 4\frac{5}{1.5^4} + \frac{5}{2^4}}{6} 2\cdot 0.5 +\frac{ \frac{5}{2^4}+ 4\frac{5}{2.5^4} + \frac{5}{3^4}}{6} 2\cdot 0.5 = 1.69155761316872
 $$


 ***Exact***

 $$
 \int_{1}^{3} \frac{5}{x^4} \, dx =  -\frac{5}{3x^3}\bigg|_1^3 = 1.60493827160
 $$

 ***error***

 $$
 \text{Right Riemann error} = | 1.6049382716 - 0.744941358| = 0.8599969
 $$

 $$
 \text{Simpsons rule error} = | 1.6049382716 - 0.744941358| =  0.0866193
 $$

 ```
 :::

 %% Cell type:code id: tags:

 ``` python
 #feel free to use this as your calculator for the exercise above
 def f(x):
    return 5/(x**4)

 f_right_riemann =
 f_simps =
 f_exact =

 print(f"f_right_riemann = {f_right_riemann}")
 print(f"f_simps = {f_simps}")
 print(f"f_exact = {f_exact}")
 print(f"Error right Riemann = {abs(f_exact-f_right_riemann)}")
 print(f"Error Simpsons = {abs(f_exact-f_simps)}")
 ```

 %% Cell type:markdown id: tags: