The Ritz method is a direct variational method to find an approximate solution for boundary value problems. The method is named after Walter Ritz.
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, 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$$
In [1]:
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()
In [2]:
x = symbols('x')
In [3]:
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
In [4]:
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
In [5]:
u2, c = approx_sol(2)
In [6]:
J = funct_J(u2, x)
Thus, the approximated solution is
In [7]:
sol = solve([diff(J, c[0]), diff(J, c[1])], c)
u_ap = u2.subs(sol)
u_ap
Out[7]:
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}$).
In [8]:
u_anal = sin(x)/sin(1) - x
In [9]:
%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
In [10]:
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
In [11]:
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
In [12]:
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()
In [12]:
In [13]:
from IPython.core.display import HTML
def css_styling():
styles = open('../styles/custom_barba.css', 'r').read()
return HTML(styles)
css_styling()
Out[13]: