(c) Juan Gomez 2019. Thanks to Universidad EAFIT for support. This material is part of the course Introduction to Finite Element Analysis

# SOLUTION: Introduction to One-Dimensional Lagrange Interpolation

## By Juan Gomez


```python
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
import sympy as sym
```

### Class activity

**Follow the steps indicated next to implement, in an independent notebook, an interpolation scheme typical of the ones used in FEM codes. The resulting scheme would correspond to the interpolation within a single element**.

## Shape functions

* Assuming that a constant interpolation domain is defined as $x \in \left[ {{-1.0},{+1.0}} \right]$ with nodal points corresponding to $x^0 = -1.0$, $x^1 = +1.0$, $x^2 = 0.0$, $x^3 = -0.25$ and $x^4 = +0.25$ use the function `LagrangPoly()` to generate the interpolation functions associated to these 5 nodal points. In a single graph plot the resulting polynomials.


```python
def LagrangPoly(x, order, i, xi=None):
    if xi == None:
        xi = sym.symbols('x:%d' % (order+1))
    index = list(range(order+1))
    index.pop(i)
    return sym.prod([(x-xi[j])/(xi[i]-xi[j]) for j in index])
```


```python
npts = 200
yy = np.zeros((npts))
xx = np.linspace(-1, 1, npts)
```


```python
x = sym.symbols('x')
pol = []
pol.append(sym.simplify(LagrangPoly(x, 4, 0, [-1, 1, 0 , -0.25 , 0.25])))
pol.append(sym.simplify(LagrangPoly(x, 4, 1, [-1, 1, 0 , -0.25 , 0.25])))
pol.append(sym.simplify(LagrangPoly(x, 4, 2, [-1, 1, 0 , -0.25 , 0.25])))
pol.append(sym.simplify(LagrangPoly(x, 4, 3, [-1, 1, 0 , -0.25 , 0.25])))
pol.append(sym.simplify(LagrangPoly(x, 4, 4, [-1, 1, 0 , -0.25 , 0.25])))
#print(pol)
```


```python
plt.figure(0)
plt.grid()
for k in range(5):
    for i in range(npts):
        yy[i] = pol[k].subs([(x, xx[i])])
    plt.plot(xx, yy)
```


    <IPython.core.display.Javascript object>






## Interpolating polynomial

* Use the interpolation functions found in the previous step, to define a matrix $\left[N(x)\right]$, which will be called in what follows the `function interpolation matrix`, in such a way that the interpolation operation can be conducted using matrix operations like:

$$\hat f(x) = \left[N(x)\right]\left\{F\right\}$$

and where $F$ is a vector storing the nodal values of the function. In your Notebook print the symbolic version of the interpolation matrix.

* In your Notebook plot a function $f(x)$ (representing temperature in a bar) and its interpolated version $\hat f(x)$ using the matrix scheme developed in the previous step. The code should plot also the first order derivative of the function.


```python
fx = lambda x: x**4 + 4.0*x**3 + 3.0*x**2 - 10.0
```


```python
npts = 200
yy = np.zeros((npts))
xx = np.linspace(-1, 1, npts)
```


```python
fd = np.array([fx(-1.0), fx(1.0), fx(0.0) , fx(-0.25) , fx(+0.25)])
print(fd)
```

    [-10.          -2.         -10.          -9.87109375  -9.74609375]



```python
print(pol[0]*fd[0] + pol[1]*fd[1] + pol[2]*fd[2] + pol[3]*fd[3] + pol[4]*fd[4])
```

    -160.0*x**4 + 170.0*x**2 - 5.0*x*(0.8*x - 0.2)*(x - 1)*(1.33333333333333*x + 0.333333333333333) - 1.0*x*(0.8*x + 0.2)*(x + 1)*(1.33333333333333*x - 0.333333333333333) + 42.1166666666667*x*(x - 1)*(x + 1)*(2.0*x - 0.5) + 41.5833333333333*x*(x - 1)*(x + 1)*(2.0*x + 0.5) - 10.0



```python
plt.figure(2)
plt.grid()
for i in range(npts):
    yy[i] = fd[0]*pol[0].subs([(x, xx[i])]) + fd[1]*pol[1].subs([(x, xx[i])]) \
        + fd[2]*pol[2].subs([(x, xx[i])]) + fd[3]*pol[3].subs([(x, xx[i])]) \
        + fd[4]*pol[4].subs([(x, xx[i])])

zz = fx(xx)
plt.plot([-1, 1, 0 , -0.25 , 0.25], fd, 'ko')
plt.plot(xx, yy, 'r--')
plt.plot(xx, zz)
```


    <IPython.core.display.Javascript object>









    [<matplotlib.lines.Line2D at 0xb1561e5f8>]



## Derivative


```python
fdx = lambda x: 4*x**3 + 12.0*x**2 + 6.0*x

fc = np.array([fdx(-1.0), fdx(1.0), fdx(0.0) , fdx(-0.25) , fdx(0.25)])
```

The interpolatted derivatives are then obtained according to:


\begin{equation}
\hat f'(x) = \frac{dL^0(x)}{dx} f^0 + \frac{dL^1(x)}{dx}(x) f^1 + \frac{dL^2(x)}{dx} f^2 + + \frac{dL^3(x)}{dx} f^3 + + \frac{dL^4(x)}{dx} f^4
\end{equation}


```python
dpol = []
dpol.append(sym.diff(pol[0], x))
dpol.append(sym.diff(pol[1], x))
dpol.append(sym.diff(pol[2], x))
dpol.append(sym.diff(pol[3], x))
dpol.append(sym.diff(pol[4], x))
#print(dpol)
#print(dpol[0]*fd[0] + dpol[1]*fd[1] + dpol[2]*fd[2] + dpol[3]*fd[3] + dpol[4]*fd[4])
#
plt.figure(3)
yy = fdx(xx)
plt.plot(xx, yy, 'r--')
plt.plot([-1, 1, 0 , -0.25 , 0.25], fc, 'ko')
#
for i in range(npts):
    yy[i] = fd[0]*dpol[0].subs([(x, xx[i])]) + fd[1] * \
        dpol[1].subs([(x, xx[i])]) + fd[2]*dpol[2].subs([(x, xx[i])]) \
    + fd[3]*dpol[3].subs([(x, xx[i])]) + fd[4]*dpol[4].subs([(x, xx[i])])
plt.plot(xx, yy)
```


    <IPython.core.display.Javascript object>









    [<matplotlib.lines.Line2D at 0xb1590ae48>]



## Matrix scheme like in the FEM

* To conduct the interpolation in the Notebook create a Python function or subroutine that takes as input parameters the coordinate $x$ where the function is to be approximated and the vector of nodal values of the temperature function and returns upon execution the interpolated value of the function.


```python
def shape1d(x , un):
    N = sym.Matrix([[x*(0.8*x - 0.2)*(x - 1)*(1.33333333333333*x + 0.333333333333333)/2,
         x*(0.8*x + 0.2)*(x + 1)*(1.33333333333333*x - 0.333333333333333)/2,
         16.0*x**4 - 17.0*x**2 + 1.0 , -4.26666666666667*x*(x - 1)*(x + 1)*(2.0*x - 0.5) ,
         -4.26666666666667*x*(x - 1)*(x + 1)*(2.0*x + 0.5) ]])
    u_int = N*un
#
    return u_int
```


```python
uu = sym.Matrix([-10.0, -2.0 , -10.0 , -9.87 , -9.74])

xx = 0.50

u_int = shape1d(xx , uu)
print(u_int)
```

    Matrix([[-8.67200000000003]])



```python
from IPython.core.display import HTML
def css_styling():
    styles = open('./nb_style.css', 'r').read()
    return HTML(styles)
css_styling()
```





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

<style>

/*
Template for Notebooks for Modelación computacional.

Based on Lorena Barba template available at:

    https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css
*/

/* Fonts */
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

/* Text */
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}	
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}	
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}
</style>









```python

```

ee84a77552c4aec484803ffd4049bdda6266a2b7
ipynb
Jupyter Notebook
solutions/SOL_01_lagrange_1d_principles.ipynb
AppliedMechanics-EAFIT/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2019-11-26T13:28:30.000Z
2022-02-16T17:57:11.000Z
solutions/SOL_01_lagrange_1d_principles.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
solutions/SOL_01_lagrange_1d_principles.ipynb
jgomezc1/Introductory-Finite-Elements
a4b44d8bf29bcd40185e51ee036f38102f9c6a72
MIT
2020-02-17T07:24:59.000Z
2022-03-02T07:54:28.000Z
Qwen/Qwen-72B
1. YES
2. YES