# Calculus of variations: Ritz method

The [Ritz method](http://en.wikipedia.org/wiki/Ritz_method) is a direct variational method to find an approximate solution for boundary value problems. The method is named after [Walter Ritz](http://en.wikipedia.org/wiki/Walther_Ritz).

## Example 1: Ordinary differential equation

Let us consider the differential equation
$$\frac{d^2u}{dx^2} + u + x = 0 \enspace , $$
with boundary conditions $u(0)=u(1)=0$.

The variational form that is equivalent to the differential equation is
$$J[u] = \int\limits_{0}^{1} F(u, u')dx = \int\limits_{0}^{1} \left[\frac{1}{2}\left( \frac{du}{dx}\right)^2 - \frac{1}{2}u^2 - ux \right] dx \enspace , $$
we can verify that computing the [Euler-Lagrange equations](http://en.wikipedia.org/wiki/Euler%E2%80%93Lagrange_equation), i.e.,
$$F_u - \frac{d}{dx}F_{u'} = 0 \enspace ,$$
where the subindices refer to derivatives with respect to $u$ and $u'$. In this case
$$\begin{align*}
F_u =& -u - x\\
F_{u'} =& u' \enspace ,
\end{align*}$$
hence
$$ -u'' - u - x = 0 \enspace ,$$
or
$$\frac{d^2 u}{dx^2} + u + x = 0 \enspace .$$

The approximating functions need to satisfy the boundary conditions, then, we propose solutions of the form
$$u_n(x) = x (1-x) \sum_{i=0}^n c_i x^i$$



```
from __future__ import division
from sympy import *
from sympy import init_printing
from sympy.solvers import solve
import matplotlib.pyplot as plt, mpld3
import numpy as np
init_printing()
```


```
x = symbols('x')
```


```
def approx_sol(num):
    """Compute the trial function using polynomials
    
    Parameters
    ----------
    num : int
        Number of terms in the expansion.
        
    Returns
    -------
    u_n : Sympy expression
        Trial function for the given basis.
    c : (num) Sympy symbols list
        List of coefficients.
        
    """    
    c = symbols('c0:%d'%num)
    u_n = x*(1 - x)*sum([c[k]*x**k for k in range(num)])
    return u_n, c
```

Let us compute the functional for this function


```
def funct_J(u, x):
    """Functional for a given trial function
    
    Parameters
    ----------
    u : Sympy expression
        Approximant solution.
    x : Sympy symbol
        Independent variable.
    
    Returns
    -------
    J : Sympy expression
        Functional for the given approximant function.
    
    """
    J = integrate(1/S(2)*(diff(u,x))**2 - 1/S(2)*u**2 - u*x, (x, 0, 1))
    return J


```

Now, let  us put all the ingredients together


```
u2, c = approx_sol(2)
```


```
J = funct_J(u2, x)
```

Thus, the approximated solution is


```
sol = solve([diff(J, c[0]), diff(J, c[1])], c)
u_ap = u2.subs(sol)
u_ap
```

The exact solution is
$$u(x) = \frac{\sin(x)}{\sin(1)} - x \enspace ,$$
and we can compare the two solutions that are really close (the relative error is $1.3778\times 10^{-5}$).


```
u_anal = sin(x)/sin(1) - x
```


```
%matplotlib inline
mpld3.enable_notebook()

vec_u = lambdify((x), u_ap, "numpy")
plt.figure(figsize=(8,4))
t = np.linspace(0, 1, 100)
y = vec_u(t)
y_anal = np.sin(t)/np.sin(1) - t
plt.plot(t, y_anal, lw=2)
plt.plot(t, y, 'r--', lw=1, markevery=10)
plt.grid(True, lw=1, color='black', alpha=0.2)
plt.xlabel("Position - x", size=15)
plt.ylabel("Solution - u", size=15)
plt.show()
```

We can compute the solution for different levels of approximation


```
plt.figure(figsize=(8,4))
t = np.linspace(0, 1, 100)
y_anal = np.sin(t)/np.sin(1) - t
plt.plot(t, y_anal, lw=2, label="Analytic")
plt.grid(True, lw=1, color='black', alpha=0.2)

t = np.linspace(0, 1, 10)
markers = ['o', 's', 'h', '*']
for cont,n in enumerate(range(1,8,2)):
    u, c = approx_sol(n)
    J = funct_J(u, x)
    sol = solve([diff(J, c[k]) for k in range(n)], c)
    u_ap = u.subs(sol)    
    vec_u = lambdify((x), u_ap, "numpy")
    y = vec_u(t)
    plt.plot(t, y, marker=markers[cont], ms=8, lw=0, label="n=%i"%n)

plt.legend(loc='best')
plt.show()
```

We can compute the error for different number of terms


```
N = range(1, 23, 2)
errors = np.zeros((len(N),))
for cont, n in enumerate(N):
    u, c = approx_sol(n)
    J = funct_J(u, x)
    sol = solve([diff(J, c[k]) for k in range(n)], c)
    u_ap = u.subs(sol) 
    rel_error = abs(integrate(u_ap - u_anal, (x,0,1))/integrate(u_anal, (x,0,1)))
    errors[cont] = rel_error
    
```


```
plt.figure(figsize=(8,4))
plt.plot(N, np.log(errors)/np.log(10), 'o')
plt.xlabel("Number of terms", size=15)
plt.ylabel("Logarithm of the relative error", size=15)
plt.grid(True, lw=1, color='black', alpha=0.2)
plt.show()
```


```

```


```
from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom_barba.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>
/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
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>






648a64083400447607191fb0c170585c8d679fdb
ipynb
Jupyter Notebook
notebooks/ritz_method.ipynb
jomorlier/FEM-Notes
3b81053aee79dc59965c3622bc0d0eb6cfc7e8ae
MIT
2020-04-15T01:53:14.000Z
2020-04-15T01:53:14.000Z
notebooks/ritz_method.ipynb
jomorlier/FEM-Notes
3b81053aee79dc59965c3622bc0d0eb6cfc7e8ae
MIT
notebooks/ritz_method.ipynb
jomorlier/FEM-Notes
3b81053aee79dc59965c3622bc0d0eb6cfc7e8ae
MIT
2020-05-25T17:19:53.000Z
2020-05-25T17:19:53.000Z
Qwen/Qwen-72B
1. YES
2. YES