```python
# 그래프, 수학 기능 추가
# Add graph and math features
import pylab as py
import numpy as np


```

# 1차 적분<br>First Order Numerical Integration



[](https://www.youtube.com/watch?v=1p0NHR5w0Lc)



다시 면적이 1인 반원을 생각해 보자.<br>Again, let's think about the half circle with area of 1.



$$
\begin{align}
    \pi r^2 &= 2 \\
    r^2 &= \frac{2}{\pi} \\
    r &= \sqrt{\frac{2}{\pi}}
\end{align}
$$




```python
r = py.sqrt(2.0 / py.pi)


```


```python
def half_circle(x):
    return py.sqrt(r**2 - x**2)


```

$$
y = \sqrt{r^2 - x^2}
$$




```python
import plot_num_int as pi


```


```python
pi.plot_a_half_circle_of_area(1)
pi.axis_equal_grid_True()


```

이번에는 사다리꼴 규칙을 이용해서 구해 보기로 하자.<br>
This time, let's use the trapezoid rule to find its area.



## 사다리꼴 규칙<br>Trapezoid Rule



다음과 같은 사다리꼴을 생각해 보자.<br>Let's think about a trapezoid as follows.




```python
x_array = (0, 1)
y_array = (1, 2)

py.fill_between(x_array, y_array)
py.axis('equal')
py.axis('off')

py.text(-0.25, 0.5, '$y_i$')
py.text(1.15, 1, '$y_{i+1}$')
py.text(0.5, -0.3, '$\Delta x$');


```

사다리꼴의 면적은 다음과 같다.<br>
Area of a trapezoid is as follows.



$$
a_i=\frac{1}{2} \left( y_i + y_{i+1} \right) \Delta x
$$



## 1차 적분<br>First order numerical integration



마찬가지로 일정 간격으로 $x$ 좌표를 나누어 보자.<br>
Same as before, let's divide $x$ coordinates in a constant interval.




```python
n = 10

pi.plot_half_circle_with_stems(n, 1)

# 사다리꼴의 좌표를 나눔 Find coordinates for the trapezoids
x_array_bar = py.linspace(-r, r, n+1)
y_array_bar = half_circle(x_array_bar)

# 각 사다리꼴의 폭 Width of each trapezoid
delta_x = x_array_bar[1] - x_array_bar[0]

# 일련의 사다리꼴을 그림 Plot a series of the trapezoids
xp, yp = x_array_bar[0], y_array_bar[0]
for x, y in zip(x_array_bar[1:], y_array_bar[1:]):
    py.fill_between((xp, x), (yp, y), alpha=0.5, color=py.random((1, 3)))
    xp, yp = x, y

py.axis('equal')
py.grid(True)


```

사다리꼴의 면적을 하나씩 구해서 더해보자.<br>Let's accumulate the area of trapezoids.



$$
    Area = \sum_{k=0}^{n-1} F_k
$$



$$
    F_k = \frac{\Delta x}{2}\left[f(x_k)+f(x_{k+1})\right]
$$



$$
    Area = \sum_{k=0}^{n-1}  \frac{1}{2}\left[f(x_k)+f(x_{k+1})\right] \Delta x
$$




```python
def get_delta_x(xi, xe, n):
    return (xe - xi) / n


```


```python
def num_int_1(f, xi, xe, n, b_verbose=False):
    x_array = py.linspace(xi, xe, n+1)
    delta_x = x_array[1] - x_array[0]

    integration_result = 0.0
    x_k = x_array[0]
    y_k = f(x_k)

    for k, x_k_plus_1 in enumerate(x_array[1:]):
        
        y_k_plus_1 = f(x_k_plus_1)
        F_k = 0.5 * (y_k + y_k_plus_1) * (x_k_plus_1 - x_k)
        if b_verbose: print('i = %2d, F_k = %g' % (k, F_k))
        integration_result += F_k
        x_k, y_k = x_k_plus_1, y_k_plus_1
    
    return integration_result


```


```python
n = 10
result = num_int_1(half_circle, -r, r, n, b_verbose=True)
print('result =', result)


```

예상한 값 1에 더 비슷한 값을 얻기 위해 더 잘게 나누어 보자<br>
To obtain the result closer to the expected value of 1, let's divide with a narrower interval.




```python
n = 100
result = num_int_1(half_circle, -r, r, n)
print('result =', result)


```


```python
%timeit -n 100 result = num_int_1(half_circle, -r, r, n)


```

### $cos \theta$의 반 주기<br>Half period of $cos \theta$




```python
n = 10
result_cos = num_int_1(py.cos, 0, py.pi, n, b_verbose=True)
print('result =', result_cos)

```


```python
n = 100
result_cos = num_int_1(py.cos, 0, py.pi, n)
print('result =', result_cos)

```

### 1/4 원<br>A quarter circle




```python
n = 10
result_quarter = num_int_1(half_circle, -r, 0, n, b_verbose=True)
print('result =', result_quarter)

```


```python
n = 100
result_quarter = num_int_1(half_circle, -r, 0, n)
print('result =', result_quarter)

```

## 연습문제<br>Exercises



도전 과제 1 : 다른 조건이 같을 때 0차 적분과 사다리꼴 적분의 오차를 비교해 보시오. 필요하면 해당 파이썬 함수를 복사하시오.<br>Try this 1 : Compare the errors of the zeroth and first order integrations of the half circle example above using the same conditions. Duplicate the python function if necessary.




```python

```

도전 과제 2 : 길이 $L=3[m]$ 인 외팔보가 분포 하중 $\omega=50sin\left(\frac{1}{2L}\pi x\right)[N/m]$을 받고 있을 때 전단력과 굽힘모멘트 선도를 구하시오.<br>
Try this 2 : Plot diagrams of shear force and bending moment of a cantilever with length $L=3m$ under distributed load $\omega=50sin\left(\frac{1}{2L}\pi x\right)[N/m]$. <br>
(ref : C 4.4, Pytel, Kiusalaas & Sharma, Mechanics of Materials, 2nd Ed, SI, Cengage Learning, 2011.)




```python

```

## 함수형 프로그래밍<br>Functional programming



간격이 일정하다면 면적의 근사값을 다음과 같이 바꾸어 쓸 수 있다.<br>
If the interval $\Delta x$ is constant, we may rewrite the approximation of the area as follows.



$$
\begin{align}
    Area &= \sum_{k=0}^{n-1}  \frac{1}{2}\left[f(x_k)+f(x_{k+1})\right] \Delta x \\
   &= \Delta x \sum_{k=0}^{n-1}  \frac{1}{2}\left[f(x_k)+f(x_{k+1})\right] 
\end{align}
$$



$$
\begin{align}
   \sum_{k=0}^{n-1}  \frac{1}{2}\left[f(x_k)+f(x_{k+1})\right] &= \frac{1}{2}\left[f(x_0)+f(x_1)\right] \\
   &+ \frac{1}{2}\left[f(x_1)+f(x_2)\right] \\
   &+ \frac{1}{2}\left[f(x_2)+f(x_3)\right] \\
   & \ldots \\
   &+ \frac{1}{2}\left[f(x_{n-2})+f(x_{n-1})\right] \\
   &+ \frac{1}{2}\left[f(x_{n-1})+f(x_{n})\right] \\
   &= \frac{1}{2}f(x_0) + \sum_{k=1}^{n-1}  f(x_k) + \frac{1}{2}f(x_{n}) \\
   &= \frac{1}{2}\left[f(x_0) + f(x_{n})\right] + \sum_{k=1}^{n-1}  f(x_k)
\end{align}
$$



$$
\begin{align}
    Area &= \Delta x \sum_{k=0}^{n-1}  \frac{1}{2}\left[f(x_k)+f(x_{k+1})\right] \\
    &= \Delta x \left[\frac{1}{2}\left[f(x_0) + f(x_{n})\right] + \sum_{k=1}^{n-1}  f(x_k)\right]
\end{align}
$$



할당문 없이 `sum()` 과 `map()` 함수로 구현해 보자.<br>
Instead of assignments, let's implement using `sum()` and `map()` functions.




```python
def num_int_1_functional(f, xi, xe, n):
    # get_delta_x() 함수 호출 횟수를 줄이기 위해 함수 안의 함수를 사용
    # To reduce the number of calling function get_delta_x(), define inner functions
    def with_delta_x(f, xi, n, delta_x=get_delta_x(xi, xe, n)):

        return delta_x * (
            0.5 * (f(xi) + f(xe))
            + sum(
                map(
                    f,
                    py.arange(xi + delta_x, xe - delta_x*0.1, delta_x),
                )
            )
        )

    return with_delta_x(f, xi, n)


```


```python
n = 100
result_func = num_int_1_functional(half_circle, -r, r, n)
print('result_func =', result_func)


```


```python
assert 1e-7 > abs(result - result_func), f"result = {result}, result_func = {result_func}"


```


```python
%timeit -n 100 result_func = num_int_1_functional(half_circle, -r, r, n)


```

## NumPy 벡터화<br>Vectorization of NumPy




```python
import pylab as py


```


```python
def num_int_1_vector_with_delta_x(f, xi, xe, n, delta_x):
    return delta_x * (
        f(py.arange(xi+delta_x, xe-delta_x*0.5, get_delta_x(xi, xe, n))).sum()
        + 0.5 * f(py.array((xi, xe))).sum()
    )


def num_int_1_vector(f, xi, xe, n):
    
    return num_int_1_vector_with_delta_x(f, xi, xe, n, get_delta_x(xi, xe, n))


```


```python
n = 100
result_vect = num_int_1_vector(half_circle, -r, r, n)
print('result_vect =', result_vect)


```


```python
assert 1e-7 > abs(result - result_vect), f"result = {result}, result_vect = {result_vect}"


```


```python
%timeit -n 100 result_func = num_int_1_vector(half_circle, -r, r, n)


```

## 시험<br>Test



아래는 함수가 맞게 작동하는지 확인함<br>
Following cells verify whether the functions work correctly.




```python
import pylab as py
r = py.sqrt(1.0 / py.pi)
n = 10
delta_x = r/n


def half_circle(x):
    return py.sqrt(r**2 - x ** 2)


assert 0.25 > num_int_1(half_circle, -r, 0, n)
assert 0.25 > num_int_1(half_circle, 0, r, n)
assert 0.25 > num_int_1_functional(half_circle, -r, 0, n)
assert 0.25 > num_int_1_functional(half_circle, 0, r, n)
assert 0.25 > num_int_1_vector(half_circle, -r, 0, n)
assert 0.25 > num_int_1_vector(half_circle, 0, r, n)


```


```python
assert 0.1 > (abs(num_int_1(half_circle, -r, 0, n) - 0.25) * 4)
assert 0.1 > (abs(num_int_1(half_circle, 0, r, n) - 0.25) * 4)
assert 0.1 > (abs(num_int_1_functional(half_circle, -r, 0, n) - 0.25) * 4)
assert 0.1 > (abs(num_int_1_functional(half_circle, 0, r, n) - 0.25) * 4)
assert 0.1 > (abs(num_int_1_vector(half_circle, -r, 0, n) - 0.25) * 4)
assert 0.1 > (abs(num_int_1_vector(half_circle, 0, r, n) - 0.25) * 4)


```

## Final Bell<br>마지막 종




```python
# stackoverfow.com/a/24634221
import os
os.system("printf '\a'");


```


```python

```

0e8f506ecc5d87f7a16efa560b159b7a7661a5b0
ipynb
Jupyter Notebook
10_first_order.ipynb
kangwonlee/19ECA-30-num-int
20f45cb6cde958ebb901341a0253abfb95efe92d
BSD-3-Clause
10_first_order.ipynb
kangwonlee/19ECA-30-num-int
20f45cb6cde958ebb901341a0253abfb95efe92d
BSD-3-Clause
10_first_order.ipynb
kangwonlee/19ECA-30-num-int
20f45cb6cde958ebb901341a0253abfb95efe92d
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES