# 10. Time dependent problem in 1D

we consider the time dependent problem
\begin{align*}
  -k_{\perp}^2 \partial_t \phi &= \partial_s J
  \\
  \partial_t A + \mu \partial_t J &= \partial_s \left( n - \phi \right)
  \\
  \partial_t n &= \partial_s J
  \\
  \beta J &= k_{\perp}^2 A
%  \label{}
\end{align*}
where $s \in [-\pi, \pi]$, $\beta \sim 10^{-3}$, $\mu \sim 10^{-4}$ and $k_{\perp} \in [10^{-2},10^{1}]$.
\\
It is easy to check that the eigenvalues related to the previous system are $\{-V_a k_{\parallel}, 0, V_a k_{\parallel} \}$ with $V_a := \frac{1+k_{\perp}^2}{\beta + \mu k_{\perp}^2}$. 


## Wave equation for $A$

Multiplying the equation on $n$ by $k_{\perp}^2$ then adding it to the equation on $\phi$, we get
$$
\left( \beta + \mu k_{\perp}^2 \right) \partial_{tt} A = \left( 1 + k_{\perp}^2 \right) \partial_{ss} A
$$
which leads to the wave equation
$$
\partial_{tt} A = \frac{1+k_{\perp}^2}{\beta + \mu k_{\perp}^2} \partial_{ss} A
$$

## Time discretization

Let's define $\gamma := \frac{k_{\perp}^2}{\beta}$, and replace $J$ in the equation on $A$. We get
\begin{align*}
  \partial_t \phi &= - \frac{1}{\beta} \partial_s A
  \\
  \partial_t A &= \frac{1}{1 + \mu \gamma} \partial_s \left( n - \phi \right)
  \\
  \partial_t n &= \gamma \partial_s A
%  \label{}
\end{align*}

Using a full implicit time scheme, we have,
\begin{align*}
   \frac{\phi^{k+1} - \phi^{k}}{\Delta t} &= - \frac{1}{\beta} \partial_s A^{k+1}
  \\
  \frac{A^{k+1} - A^{k}}{\Delta t} &= \frac{1}{1 + \mu \gamma} \partial_s \left( n^{k+1} - \phi^{k+1} \right)
  \\
  \frac{n^{k+1} - n^{k}}{\Delta t} &= \gamma \partial_s A^{k+1}
%  \label{}
\end{align*}

finally, 

\begin{align*}
    \phi^{k+1} + \frac{\Delta t}{\beta} \partial_s A^{k+1}  &= \phi^k
    \\
    \frac{\Delta t}{1+\mu \gamma} \phi^{k+1} + A^{k+1} - \frac{\Delta t}{1+\mu \gamma} \partial_s n^{k+1} &= A^k
    \\
    -\Delta t \gamma \partial_s A^{k+1} + n^{k+1} &= n^k
\end{align*}



## Weak formulation

Let $v$ denote a test function, in a Finite Elements space $V \subset H^1(\Omega)$. Multiplying all the equations by $v$, then integrating over the whole domain, we get

\begin{align*}
    \langle \phi^{k+1}, v \rangle 
    + \frac{\Delta t}{\beta} \langle \partial_s A^{k+1}, v \rangle  &= \langle \phi^k, v \rangle
    \\
    \frac{\Delta t}{1+\mu \gamma} \langle \phi^{k+1}, v \rangle + \langle A^{k+1}, v \rangle - \frac{\Delta t}{1+\mu \gamma} \langle \partial_s n^{k+1}, v \rangle &= \langle A^k, v \rangle
    \\
    -\Delta t \gamma \langle \partial_s A^{k+1}, v \rangle + \langle n^{k+1}, v \rangle &= \langle n^k, v \rangle
\end{align*}

We use a symmetrized weak formulation, where we assume having periodic boundary conditions:




\begin{align*}
    \langle \phi^{k+1}, v \rangle 
    + \frac{\Delta t}{2 \beta} \langle \partial_s A^{k+1}, v \rangle 
    - \frac{\Delta t}{2 \beta} \langle A^{k+1}, \partial_s v \rangle
    &= \langle \phi^k, v \rangle
    \\
    \frac{\Delta t}{1+\mu \gamma} \langle \phi^{k+1}, v \rangle 
    + \langle A^{k+1}, v \rangle 
    - \frac{\Delta t}{2+2\mu \gamma} \langle \partial_s n^{k+1}, v \rangle 
    + \frac{\Delta t}{2+2\mu \gamma} \langle n^{k+1}, \partial_s v \rangle 
    &= \langle A^k, v \rangle
    \\
    -\frac{\Delta t \gamma}{2} \langle \partial_s A^{k+1}, v \rangle 
    +\frac{\Delta t \gamma}{2} \langle A^{k+1}, \partial_s v \rangle 
    + \langle n^{k+1}, v \rangle 
    &= \langle n^k, v \rangle
\end{align*}

In order to simplify the notation, we introduction the following bilinear form 
$$
b( v,u ) := \frac{1}{2} \left( \langle \partial_s u, v \rangle - \langle u, \partial_s v \rangle \right)
$$
then our weak formulation writes

\begin{align*}
    \langle \phi^{k+1}, v \rangle 
    + \frac{\Delta t}{\beta} b(v,A^{k+1}) 
    &= \langle \phi^k, v \rangle
    \\
    \frac{\Delta t}{1+\mu \gamma} \langle \phi^{k+1}, v \rangle 
    + \langle A^{k+1}, v \rangle 
    - \frac{\Delta t}{1+\mu \gamma} b(v,n^{k+1}) 
    &= \langle A^k, v \rangle
    \\
    -\Delta t \gamma b(v,A^{k+1}) 
    + \langle n^{k+1}, v \rangle 
    &= \langle n^k, v \rangle
\end{align*}

Finally, let's introduce the weak formulation related to the mass matrix $a_m(v,u) := \langle u, v \rangle$, 

\begin{align*}
    a_m(v, \phi^{k+1}) 
    + \frac{\Delta t}{\beta} b(v, A^{k+1})  
    &= a_m(v, \phi^k)
    \\
    \frac{\Delta t}{1+\mu \gamma} a_m(v, \phi^{k+1}) 
    + a_m(v, A^{k+1}) 
    - \frac{\Delta t}{1+\mu \gamma} b(v,  n^{k+1}) 
    &= a_m(v, A^k)
    \\
    - \Delta t \gamma b(v, A^{k+1}) 
    + a_m(v, n^{k+1}) 
    &= a_m(v, n^k)
\end{align*}


```python
import numpy as np
from numpy import linspace, zeros, pi

from sympy.core.containers import Tuple
from sympy import symbols
from sympy import Symbol
from sympy import Lambda
from sympy import Function

from gelato.glt import glt_symbol
from gelato.calculus   import (Dot, Cross, Grad, Curl, Rot, Div, dx)
from gelato.calculus   import Constant
from gelato.fem.assembly import assemble_matrix
from gelato.fem.utils    import compile_kernel
from gelato.fem.utils    import compile_symbol

from spl.fem.splines import SplineSpace
from spl.fem.vector  import VectorFemSpace

from IPython.display import Math
from sympy import latex
```


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

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

a_m = lambda v,u: u*v
b = lambda v,u: 0.5*(dx(u)*v - u*dx(v))
```

\begin{align*}
    a_m(v, \phi^{k+1}) 
    + \frac{\Delta t}{\beta} b(v, A^{k+1})  
    &= a_m(v, \phi^k)
    \\
    \frac{\Delta t}{1+\mu \gamma} a_m(v, \phi^{k+1}) 
    + a_m(v, A^{k+1}) 
    - \frac{\Delta t}{1+\mu \gamma} b(v,  n^{k+1}) 
    &= a_m(v, A^k)
    \\
    - \Delta t \gamma b(v, A^{k+1}) 
    + a_m(v, n^{k+1}) 
    &= a_m(v, n^k)
\end{align*}


```python
phi, A, n = symbols('phi A n')
v0, v1, v2 = symbols('v0 v1 v2')

dt = Constant('dt')
beta = Constant('beta')
mu = Constant('mu')
gamma = Constant('gamma')

a = Lambda((x,v0,v1,v2,phi,A,n), 
             a_m(v0, phi) + dt/beta * b(v0, A) 
           + dt/(1+mu * gamma) * a_m(v1, phi) 
           + a_m(v1, A) 
           - dt/(1+mu * gamma) * b(v1, n) 
           - dt * gamma * b(v2, A) + a_m(v2, n))
```


```python
# create a finite element space
p  = 3
ne = 64

grid = linspace(0., 1., ne+1)

W = SplineSpace(p, grid=grid)
V = VectorFemSpace(W, W, W)
```


```python
symbol = glt_symbol(a, space=V)
```


```python
Math(latex(symbol))
```




$$\left[\begin{matrix}m_{1} & \frac{1.0 a_{1}}{\beta} dt & 0\\\frac{dt m_{1}}{\gamma \mu + 1} & m_{1} & - \frac{1.0 dt a_{1}}{\gamma \mu + 1}\\0 & - 1.0 dt \gamma a_{1} & m_{1}\end{matrix}\right]$$




```python
eigen = symbol.eigenvals()
eigen = list(eigen.keys())
```


```python
Math(latex(eigen))
```




$$\left [ m_{1}, \quad \frac{1}{\beta \gamma \mu + \beta} \left(\beta \gamma \mu m_{1} + \beta m_{1} + dt \sqrt{\beta^{2} \gamma^{2} \mu a_{1}^{2} + \beta^{2} \gamma a_{1}^{2} + \beta \gamma \mu a_{1} m_{1} + \beta a_{1} m_{1}}\right), \quad \frac{1}{\beta \gamma \mu + \beta} \left(\beta \gamma \mu m_{1} + \beta m_{1} - dt \sqrt{\beta^{2} \gamma^{2} \mu a_{1}^{2} + \beta^{2} \gamma a_{1}^{2} + \beta \gamma \mu a_{1} m_{1} + \beta a_{1} m_{1}}\right)\right ]$$




```python
# compute the symbol of the mass
symbol_m = glt_symbol(Lambda((x,v,u), u*v), space=V)
symbol_a = glt_symbol(Lambda((x,v,u), dx(u)*v), space=V)
```


```python
eigen_normalized = [e/symbol_m for e in eigen]
```


```python
Math(latex(eigen_normalized))
```




$$\left [ 1, \quad \frac{1}{m_{1} \left(\beta \gamma \mu + \beta\right)} \left(\beta \gamma \mu m_{1} + \beta m_{1} - dt \sqrt{\beta^{2} \gamma^{2} \mu a_{1}^{2} + \beta^{2} \gamma a_{1}^{2} - \beta \gamma \mu a_{1} m_{1} - \beta a_{1} m_{1}}\right), \quad \frac{1}{m_{1} \left(\beta \gamma \mu + \beta\right)} \left(\beta \gamma \mu m_{1} + \beta m_{1} + dt \sqrt{\beta^{2} \gamma^{2} \mu a_{1}^{2} + \beta^{2} \gamma a_{1}^{2} - \beta \gamma \mu a_{1} m_{1} - \beta a_{1} m_{1}}\right)\right ]$$




```python
e = eigen_normalized[1]
```


```python
from sympy import simplify, cancel, collect, expand
```


```python
Math(latex(cancel(e-1)))
```




$$- \frac{dt}{\beta \gamma \mu m_{1} + \beta m_{1}} \sqrt{\beta^{2} \gamma^{2} \mu a_{1}^{2} + \beta^{2} \gamma a_{1}^{2} - \beta \gamma \mu a_{1} m_{1} - \beta a_{1} m_{1}}$$




```python
Math(latex(symbol_a))
```




$$a_{1}$$




```python
print(symbol_a.is_complex)
```

    None



```python

```


```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

```

1e6b6770f67a1d562d7c00fdb2b5fdaa9f0073af
ipynb
Jupyter Notebook
lessons/Chapter8/20_time_dependent_1d.ipynb
ratnania/IGA-Python
a9d7aa9bd14d4b3f1b12cdfbc2f9bf3c0a68fff4
MIT
2018-04-27T15:40:17.000Z
2020-08-13T08:45:35.000Z
lessons/Chapter8/20_time_dependent_1d.ipynb
GabrielJie/IGA-Python
a9d7aa9bd14d4b3f1b12cdfbc2f9bf3c0a68fff4
MIT
2021-06-08T22:59:19.000Z
2022-01-17T20:36:56.000Z
lessons/Chapter8/20_time_dependent_1d.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