# Weighted-residual method

Let us consider the equation

$$A u = f\quad \text{in } \Omega$$

For an approximation $u_N$ of $u$, the residual, $R_N$, is defined by

$$R_N \equiv Au_N - f$$

When the residual is made orthogonal to the subspace spanned by a base $\{\psi_k\}$, we have a weighted-residual method, i.e.,

$$\langle R_N, \psi_k\rangle =0 \quad k=1, 2, \cdots, N$$


In [1]:

from __future__ import division, print_function
import numpy as np
from sympy import *
from sympy.plotting import plot3d
from scipy.linalg import eigh
from scipy.special import jn_zeros as Jn_zeros, jn as Jn
import matplotlib.pyplot as plt




In [2]:

init_session()
%matplotlib inline
plt.style.use("seaborn-notebook")




IPython console for SymPy 1.2 (Python 3.6.6-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.2/



## Ritz method: Axisymmetric modes in a circular membrane



In [3]:

def u_fun(r, m):
""" Trial function. """
c = symbols('c0:%i' % m)
w = (1 - r**2) *sum(c[k]*r**(2*k) for k in range (0, m))
return w, c




In [4]:

r = symbols('r')
m = 7
u, coef = u_fun(r, m)




In [5]:

T_inte = u**2
U_inte = diff(u, r)**2




In [6]:

display(U_inte)
display(T_inte)




$$\left(- 2 r \left(c_{0} + c_{1} r^{2} + c_{2} r^{4} + c_{3} r^{6} + c_{4} r^{8} + c_{5} r^{10} + c_{6} r^{12}\right) + \left(- r^{2} + 1\right) \left(2 c_{1} r + 4 c_{2} r^{3} + 6 c_{3} r^{5} + 8 c_{4} r^{7} + 10 c_{5} r^{9} + 12 c_{6} r^{11}\right)\right)^{2}$$

$$\left(- r^{2} + 1\right)^{2} \left(c_{0} + c_{1} r^{2} + c_{2} r^{4} + c_{3} r^{6} + c_{4} r^{8} + c_{5} r^{10} + c_{6} r^{12}\right)^{2}$$




In [7]:

U = integrate(expand(r*U_inte), (r, 0, 1))
T = integrate(expand(r*T_inte), (r, 0, 1))




In [8]:

K = Matrix(m, m, lambda ii, jj: diff(U, coef[ii], coef[jj]))
K




Out[8]:

$$\left[\begin{matrix}2 & \frac{2}{3} & \frac{1}{3} & \frac{1}{5} & \frac{2}{15} & \frac{2}{21} & \frac{1}{14}\\\frac{2}{3} & \frac{2}{3} & \frac{7}{15} & \frac{1}{3} & \frac{26}{105} & \frac{4}{21} & \frac{19}{126}\\\frac{1}{3} & \frac{7}{15} & \frac{2}{5} & \frac{34}{105} & \frac{11}{42} & \frac{3}{14} & \frac{8}{45}\\\frac{1}{5} & \frac{1}{3} & \frac{34}{105} & \frac{2}{7} & \frac{31}{126} & \frac{19}{90} & \frac{2}{11}\\\frac{2}{15} & \frac{26}{105} & \frac{11}{42} & \frac{31}{126} & \frac{2}{9} & \frac{98}{495} & \frac{29}{165}\\\frac{2}{21} & \frac{4}{21} & \frac{3}{14} & \frac{19}{90} & \frac{98}{495} & \frac{2}{11} & \frac{71}{429}\\\frac{1}{14} & \frac{19}{126} & \frac{8}{45} & \frac{2}{11} & \frac{29}{165} & \frac{71}{429} & \frac{2}{13}\end{matrix}\right]$$




In [9]:

M = Matrix(m, m, lambda ii, jj: diff(T, coef[ii], coef[jj]))
M




Out[9]:

$$\left[\begin{matrix}\frac{1}{3} & \frac{1}{12} & \frac{1}{30} & \frac{1}{60} & \frac{1}{105} & \frac{1}{168} & \frac{1}{252}\\\frac{1}{12} & \frac{1}{30} & \frac{1}{60} & \frac{1}{105} & \frac{1}{168} & \frac{1}{252} & \frac{1}{360}\\\frac{1}{30} & \frac{1}{60} & \frac{1}{105} & \frac{1}{168} & \frac{1}{252} & \frac{1}{360} & \frac{1}{495}\\\frac{1}{60} & \frac{1}{105} & \frac{1}{168} & \frac{1}{252} & \frac{1}{360} & \frac{1}{495} & \frac{1}{660}\\\frac{1}{105} & \frac{1}{168} & \frac{1}{252} & \frac{1}{360} & \frac{1}{495} & \frac{1}{660} & \frac{1}{858}\\\frac{1}{168} & \frac{1}{252} & \frac{1}{360} & \frac{1}{495} & \frac{1}{660} & \frac{1}{858} & \frac{1}{1092}\\\frac{1}{252} & \frac{1}{360} & \frac{1}{495} & \frac{1}{660} & \frac{1}{858} & \frac{1}{1092} & \frac{1}{1365}\end{matrix}\right]$$




In [10]:

Kn = np.array(K).astype(np.float64)
Mn = np.array(M).astype(np.float64)




In [11]:

vals, vecs = eigh(Kn, Mn, eigvals=(0, m-1))
np.sqrt(vals)




Out[11]:

array([ 2.40482556,  5.52007811,  8.65373016, 11.79598495, 15.24615171,
21.46269268, 41.26282741])




In [23]:

lam = Jn_zeros(0, m)




In [13]:

r_vec = np.linspace(0, 1, 60)
plt.figure(figsize=(14, 5))
ax1 = plt.subplot(1, 2, 1)
ax2 = plt.subplot(1, 2, 2)
for num in range(5):
u_num = lambdify((r), u.subs({coef[kk]: vecs[kk, num] for kk in range(m)}),
"numpy")
ax1.plot(r_vec, u_num(r_vec)/u_num(0))
ax1.set_title("Approximated solution")
ax2.plot(r_vec, Jn(0, lam[num]*r_vec), label=r"$m=%i$"%num)
ax2.set_title("Exact solution")
plt.legend(loc="best", framealpha=0);






## Bubnov-Galerkin method

The Bubnov-Galerkint methood is a generalization of the Ritz method. As in the Ritz method, it seeks an approximate solution as a linear combination of base functions

$$u_N = \sum_{i=1}^{N} c_i \phi_i\, .$$

In this case, the coefficients $c_i$ are determined from the condition that the residual $R_N$ is orthogonal to the basis functions $\phi_1, \phi_2, \cdots, \phi_N$:

$$\langle R_N, \phi_k\rangle =0\quad k=1, 2, \cdots, N \, .$$

If the operator is positive-definite then $A=T^*T$ and the problem can be rewritten identically to the Ritz method.



In [14]:

def u_fun(r, m):
""" Trial function. """
c = symbols('c0:%i' % m)
w = (1 - r**2) *sum(c[k]*r**(2*k) for k in range (0, m))
return w, c




In [15]:

r = symbols('r')
m = 7
u, coef = u_fun(r, m)




In [16]:

u




Out[16]:

$$\left(- r^{2} + 1\right) \left(c_{0} + c_{1} r^{2} + c_{2} r^{4} + c_{3} r^{6} + c_{4} r^{8} + c_{5} r^{10} + c_{6} r^{12}\right)$$




In [17]:

U = -integrate(diff(r*diff(u, r), r)*u, (r, 0, 1))
T = integrate(r*u**2, (r, 0, 1))




In [18]:

K = Matrix(m, m, lambda ii, jj: diff(U, coef[ii], coef[jj]))
K




Out[18]:

$$\left[\begin{matrix}2 & \frac{2}{3} & \frac{1}{3} & \frac{1}{5} & \frac{2}{15} & \frac{2}{21} & \frac{1}{14}\\\frac{2}{3} & \frac{2}{3} & \frac{7}{15} & \frac{1}{3} & \frac{26}{105} & \frac{4}{21} & \frac{19}{126}\\\frac{1}{3} & \frac{7}{15} & \frac{2}{5} & \frac{34}{105} & \frac{11}{42} & \frac{3}{14} & \frac{8}{45}\\\frac{1}{5} & \frac{1}{3} & \frac{34}{105} & \frac{2}{7} & \frac{31}{126} & \frac{19}{90} & \frac{2}{11}\\\frac{2}{15} & \frac{26}{105} & \frac{11}{42} & \frac{31}{126} & \frac{2}{9} & \frac{98}{495} & \frac{29}{165}\\\frac{2}{21} & \frac{4}{21} & \frac{3}{14} & \frac{19}{90} & \frac{98}{495} & \frac{2}{11} & \frac{71}{429}\\\frac{1}{14} & \frac{19}{126} & \frac{8}{45} & \frac{2}{11} & \frac{29}{165} & \frac{71}{429} & \frac{2}{13}\end{matrix}\right]$$




In [19]:

M = Matrix(m, m, lambda ii, jj: diff(T, coef[ii], coef[jj]))
M




Out[19]:

$$\left[\begin{matrix}\frac{1}{3} & \frac{1}{12} & \frac{1}{30} & \frac{1}{60} & \frac{1}{105} & \frac{1}{168} & \frac{1}{252}\\\frac{1}{12} & \frac{1}{30} & \frac{1}{60} & \frac{1}{105} & \frac{1}{168} & \frac{1}{252} & \frac{1}{360}\\\frac{1}{30} & \frac{1}{60} & \frac{1}{105} & \frac{1}{168} & \frac{1}{252} & \frac{1}{360} & \frac{1}{495}\\\frac{1}{60} & \frac{1}{105} & \frac{1}{168} & \frac{1}{252} & \frac{1}{360} & \frac{1}{495} & \frac{1}{660}\\\frac{1}{105} & \frac{1}{168} & \frac{1}{252} & \frac{1}{360} & \frac{1}{495} & \frac{1}{660} & \frac{1}{858}\\\frac{1}{168} & \frac{1}{252} & \frac{1}{360} & \frac{1}{495} & \frac{1}{660} & \frac{1}{858} & \frac{1}{1092}\\\frac{1}{252} & \frac{1}{360} & \frac{1}{495} & \frac{1}{660} & \frac{1}{858} & \frac{1}{1092} & \frac{1}{1365}\end{matrix}\right]$$




In [20]:

Kn = np.array(K).astype(np.float64)
Mn = np.array(M).astype(np.float64)




In [21]:

vals, vecs = eigh(Kn, Mn, eigvals=(0, m-1))
np.sqrt(vals)




Out[21]:

array([ 2.40482556,  5.52007811,  8.65373016, 11.79598495, 15.24615171,
21.46269268, 41.26282741])



The following cell change the style of the notebook.



In [22]:

from IPython.core.display import HTML
def css_styling():
return HTML(styles)
css_styling()




Out[22]:

/* 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 )
}

MathJax.Hub.Config({
TeX: {
extensions: ["AMSmath.js"]
},
tex2jax: {
inlineMath: [ ['$','$'], ["\$","\$"] ],
displayMath: [ ['$$','$$'], ["\$","\$"] ]
},
displayAlign: 'center', // Change this to 'center' to center equations.
"HTML-CSS": {
styles: {'.MathJax_Display': {"margin": 4}}
}
});




In [ ]: