# 1. High order pathology for the Stiffness matrix

For this chapter, we will need the [GeLaTo](https://github.com/pyccel/GeLaTo) library.

In this example, we show the dependency of the **stiffness** symbol with respect to the spline degree in 1d, we will also give some asymptotic expansion of the basic operators: **stiffness**, **mass** and **advection**.


```python
# ... needed imports
from numpy import linspace, zeros, pi
import matplotlib.pyplot as plt
%matplotlib inline

from IPython.display import Math

from sympy import Symbol
from sympy import Lambda

from gelato.glt import glt_symbol
from gelato.glt import glt_lambdify
from gelato.calculus   import Dot, Grad, dx

from spl.fem.splines import SplineSpace
```

We define the weak form, as:


```python
x = Symbol('x')

u = Symbol('u')
v = Symbol('v')

a_s = Lambda((x,v,u), Dot(Grad(u), Grad(v)))
```


```python
V = SplineSpace(2, grid=linspace(0., 1., 32+1))
```


```python
from IPython.display import Math
from gelato.printing.latex import glt_latex

symbol = glt_symbol(a_s, space=V, evaluate=True)
Math(glt_latex(symbol))
```




$$\left( \left ( x, \quad \theta_1\right ) \mapsto - \frac{64}{3} \cos{\left (\theta_1 \right )} - \frac{32}{3} \cos{\left (2 \theta_1 \right )} + 32 \right)$$




```python
from gelato.fem.utils    import compile_symbol
symbol_s = compile_symbol('symbol_1', a_s, V, backend='fortran')
```


```python
# ... make the symbol a numeric function, that can be evaluated
f = glt_lambdify(symbol)

nx = 100
x = linspace(  0, 1, nx)
t = linspace(-pi,pi, nx)
w = f(x,t)
# ...

# ... example of symbol evaluation
e = zeros(nx)
symbol_s(x, t, e)
# ...

```


```python
from numpy import allclose

allclose(e,w)
```




    True




```python

```

The following function allows to plot the eigenvalues for a given *grid* and different spline **degrees**:


```python
def plot_symbols(expr, degrees=[2, 3], ne=64):
    
    colors = ["r", "b", "g", "y", "m", "k", "c"]

    for i,p in enumerate(degrees):
        # ...  create a finite element space
        V = SplineSpace(p, grid=linspace(0., 1., ne+1))
        # ...

        # ... compute the GLT symbol
        symbol = glt_symbol(expr, space=V, evaluate=True)
        # ...

        # ... make the symbol a numeric function, that can be evaluated
        f = glt_lambdify(symbol)

        nx = 100
        x1 = linspace(  0, 1, nx)
        t1 = linspace(-pi,pi, nx)
        w = f(x1,t1)
        # ...
        
        plt.plot(t1, w.real, "-"+colors[i], label="$p=" + str(p) + "$")
        plt.legend(loc=9);
```


```python
plot_symbols(a_s, degrees=[2, 3, 4, 5])
```

As expected, we notice the convergence to 0 at $| \theta | = \pi$.

We consider now the weak formulation related to the **mass** matrix:


```python
x = Symbol('x')

u = Symbol('u')
v = Symbol('v')

# weak form for the mass matrix
a_m = Lambda((x,v,u), u * v)
```

We will use the following function to plot the division of the symbol of the **stiffness** by the symbol of the **mass**  of degree $p-1$


```python
def plot_symbols(expr_s, expr_m, degrees=[2, 3], ne=64):
    
    colors = ["r", "b", "g", "y", "m", "k", "c"]

    for i,p in enumerate(degrees):
        # ...  create a finite element space
        V = SplineSpace(p, grid=linspace(0., 1., ne+1))
        W = SplineSpace(p-1, grid=linspace(0., 1., ne+1))
        # ...

        # ... compute the GLT symbol
        symbol_s = glt_symbol(expr_s, space=V, evaluate=True)
        symbol_m = glt_symbol(expr_m, space=W, evaluate=True)
        # ...

        # ... make the symbol a numeric function, that can be evaluated
        f_s = glt_lambdify(symbol_s)
        f_m = glt_lambdify(symbol_m)

        nx = 100
        x = linspace(  0, 1, nx)
        t = linspace(-pi,pi, nx)
        w = f_s(x,t) / f_m(x,t)
        # ...
        
        plt.plot(t, w.real, "-"+colors[i], label="$p=" + str(p) + "$")
        plt.legend(loc=9);
```


```python
plot_symbols(a_s, a_m, degrees=[2, 3, 4, 5])
```

We see now that the *new symbol* is not sensitive to the *B-Spline* degree.

In fact, this symbol corresponds to the $S M^{-1}$ (and $M^{-1} S$), if we denote by $S$ and $M$ the *stiffness* and *mass* matrices respectively.

## Symbolic expression and asymptotic expansion 

In this section, we give the asymptotic expansion of the symbol associated to the **mass**, **stiffness** and **advection** matrices.

We start by defining our weak formulations:


```python
x = Symbol('x')

u = Symbol('u')
v = Symbol('v')

a_m = Lambda((x,v,u), u * v)
a_s = Lambda((x,v,u), Dot(Grad(u), Grad(v)))
a_a = Lambda((x,v,u), dx(v) * u)
```


```python
# ...
p  = 3
ne = 16
# ...

# ...  create a finite element space
V = SplineSpace(p, grid=linspace(0., 1., ne+1))
# ...

# ... compute the GLT symbol
symbol_m = glt_symbol(a_m, space=V, evaluate=True)
symbol_s = glt_symbol(a_s, space=V, evaluate=True)
symbol_a = glt_symbol(a_a, space=V, evaluate=True)
# ...
```

We can now evaluate, symbolicaly, the symbols at the sampling points, using the following function, and print them using *Latex*:


```python
from gelato.printing.latex import glt_latex
from IPython.display import Math

t1 = Symbol("t1")

# ...
def to_latex(expr):    
    m = Symbol("m")
    N = Symbol("N")
    pi = Symbol("pi")
    return glt_latex(expr.subs({t1: 2 * pi * m / N}))
# ...
```

---
**Note**:
since **expr** is a sympy **Lambda** expression, we need to use its *expr* member

---


```python
# symbol of the mass matrix
Math(to_latex(symbol_m.expr))
```




$$\frac{397}{13440} \cos{\left (\frac{2 m}{N} \pi \right )} + \frac{1}{336} \cos{\left (\frac{4 m}{N} \pi \right )} + \frac{1}{40320} \cos{\left (\frac{6 m}{N} \pi \right )} + \frac{2396825387}{79999999672}$$




```python
# symbol of the stiffness matrix
Math(to_latex(symbol_s.expr))
```




$$- 4 \cos{\left (\frac{2 m}{N} \pi \right )} - \frac{32}{5} \cos{\left (\frac{4 m}{N} \pi \right )} - \frac{4}{15} \cos{\left (\frac{6 m}{N} \pi \right )} + \frac{32}{3}$$




```python
# symbol of the advection matrix
Math(to_latex(symbol_a.expr))
```




$$- i \left(- \frac{49}{72} \sin{\left (\frac{2 m}{N} \pi \right )} - \frac{7}{45} \sin{\left (\frac{4 m}{N} \pi \right )} - \frac{1}{360} \sin{\left (\frac{6 m}{N} \pi \right )}\right)$$



We now give the asymptotic expansion of the GLT symbol arount $0$, up to the order $8$:


```python
# asymptotic expansion of the mass symbol
Math(to_latex(symbol_m.expr.series(t1,0,8)))
```




$$- \frac{4 m^{6} \pi^{6}}{189 N^{6}} + \frac{19 m^{4} \pi^{4}}{360 N^{4}} - \frac{m^{2} \pi^{2}}{12 N^{2}} + \frac{1574999993543}{25199999896680} + \mathcal{O}\left(\frac{m^{8} \pi^{8}}{N^{8}}; \left ( N, \quad m, \quad \pi\right )\rightarrow \left ( 0, \quad 0, \quad 0\right )\right)$$




```python
# asymptotic expansion of the stiffness symbol
Math(to_latex(symbol_s.expr.series(t1,0,8)))
```




$$\frac{2432 m^{6} \pi^{6}}{45 N^{6}} - \frac{256 m^{4} \pi^{4}}{3 N^{4}} + \frac{64 m^{2}}{N^{2}} \pi^{2} + \mathcal{O}\left(\frac{m^{8} \pi^{8}}{N^{8}}; \left ( N, \quad m, \quad \pi\right )\rightarrow \left ( 0, \quad 0, \quad 0\right )\right)$$




```python
# asymptotic expansion of the advection symbol
Math(to_latex(symbol_a.expr.series(t1,0,8)))
```




$$- \frac{128 i m^{7} \pi^{7}}{189 N^{7}} + \frac{76 i m^{5} \pi^{5}}{45 N^{5}} - \frac{8 i m^{3} \pi^{3}}{3 N^{3}} + \frac{2 i}{N} m \pi + \mathcal{O}\left(\frac{m^{8} \pi^{8}}{N^{8}}; \left ( N, \quad m, \quad \pi\right )\rightarrow \left ( 0, \quad 0, \quad 0\right )\right)$$




```python

```


```python

```


```python
fs = glt_lambdify(symbol_s)
fm = glt_lambdify(symbol_m)
fa = glt_lambdify(symbol_a)

nx = 100
x1 = linspace(  0, 1, nx)
t1 = linspace(-pi,pi, nx)

s1 = fs(x1,t1)
m1 = fm(x1,t1)
a1 = fa(x1,t1)
```


```python
w = a1**2 + m1 * s1
plt.plot(t1, w.real);
```


```python
s2 = s1
m2 = m1
a2 = a1

w = (m1*s2 + m2*s1)**2 + 4*(a1**2*m2**2 + a2**2*m1**2)

plt.plot(t1, w.real);




```


```python
# ...
p1  = 3 ; p2 = 3
ne1 = 64 ; ne2 = 64
# ...

# ...  create a finite element space
V1 = SplineSpace(p1, grid=linspace(0., 1., ne1+1))
V2 = SplineSpace(p2, grid=linspace(0., 1., ne2+1))
# ...

# ... compute the GLT symbol
symbol_m1 = glt_symbol(a_m, space=V1, evaluate=True)
symbol_s1 = glt_symbol(a_s, space=V1, evaluate=True)
symbol_a1 = glt_symbol(a_a, space=V1, evaluate=True)

symbol_m2 = glt_symbol(a_m, space=V2, evaluate=True)
symbol_s2 = glt_symbol(a_s, space=V2, evaluate=True)
symbol_a2 = glt_symbol(a_a, space=V2, evaluate=True)
# ...

# ...
fs1 = glt_lambdify(symbol_s1)
fm1 = glt_lambdify(symbol_m1)
fa1 = glt_lambdify(symbol_a1)

fs2 = glt_lambdify(symbol_s2)
fm2 = glt_lambdify(symbol_m2)
fa2 = glt_lambdify(symbol_a2)
# ...

x1 = x2 = linspace(  0, 1, 100)
t1 = t2 = linspace(-pi,pi, 100)

# ...
s1 = fs1(x1,t1)
m1 = fm1(x1,t1)
a1 = fa1(x1,t1)

s2 = fs2(x2,t2)
m2 = fm2(x2,t2)
a2 = fa2(x2,t2)
# ...
```


```python
w = (m1*s2 + m2*s1)**2 + 4*(a1**2*m2**2 + a2**2*m1**2)
print(min(w.real))
plt.plot(t1, w.real);
```


```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../../styles/custom.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>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }
    div.cell{
        width:600px;
        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:#054BCD;
        margin-bottom: 0.5em;
        margin-top: 0.5em;
        display: block;
    }	
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 16pt;
        color: #054BCD;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }  
</style>






```python

```


```python

```

f57d4acca9b58f98dd5947859a7ad908a3cfa987
ipynb
Jupyter Notebook
lessons/Chapter8/01_stiffness_pathology.ipynb
ratnania/IGA-Python
a9d7aa9bd14d4b3f1b12cdfbc2f9bf3c0a68fff4
MIT
2018-04-27T15:40:17.000Z
2020-08-13T08:45:35.000Z
lessons/Chapter8/01_stiffness_pathology.ipynb
GabrielJie/IGA-Python
a9d7aa9bd14d4b3f1b12cdfbc2f9bf3c0a68fff4
MIT
2021-06-08T22:59:19.000Z
2022-01-17T20:36:56.000Z
lessons/Chapter8/01_stiffness_pathology.ipynb
GabrielJie/IGA-Python
a9d7aa9bd14d4b3f1b12cdfbc2f9bf3c0a68fff4
MIT
2018-10-06T01:30:20.000Z
2021-12-31T02:42:05.000Z
Qwen/Qwen-72B
1. YES
2. YES