{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook introduces the Galerkin method with an example of a linear boundary value problem."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using Printf\n",
"Base.show(io::IO, f::Float64) = @printf(io, \"%.7e\", f)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider $y'' = 4 y$, for $x$ in $[0, 1]$, $y(0) = 1$ and $y(1) = 3$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This problem has an exact solution, evaluated in the function below."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"exactsol"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" exactsol(x::Float64)\n",
"\n",
"returns the value for the exact solution at x.\n",
"\"\"\"\n",
"function exactsol(x::Float64)\n",
" d = exp(2.0) - exp(-2.0)\n",
" a = (3.0 - exp(-2.0))*exp(2*x)\n",
" b = (exp(2.0) - 3.0)*exp(-2*x)\n",
" return (a + b)/d\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 1. Collocation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The idea of collocation is to impose a candidate solution, a solution of a particular form.\n",
"The first type of solution is a polynomial, expressed in the standard basis 1, $x$, $x^2$, etc.\n",
"The coefficients of this polynomial are obtained via the solution of a linear system, defined by the function below."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"setup1"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" A, b = setup1(dim::Int)\n",
"\n",
"returns the linear system to solve a BVP\n",
"with collocation over the interval [a,b]\n",
"with dim points, with a step size h of (b-a)/(dim-1).\n",
"\"\"\"\n",
"function setup1(dim::Int, a::Float64, b::Float64)\n",
" h = (b - a)/(dim-1)\n",
" A = zeros(dim,dim); A[1,1] = 1\n",
" A[dim,:] = ones(dim)\n",
" xi = h\n",
" for i=2:dim-1\n",
" A[i,1] = -4.0\n",
" A[i,2] = -4.0*xi\n",
" A[i,3] = -4*xi^2 + 2\n",
" if dim > 3\n",
" A[i,4] = -4*xi^3 + 6*xi\n",
" for j=5:dim\n",
" A[i,j] = -4*xi^(j-1) + (j-1)*(j-2)*xi^(j-3)\n",
" end\n",
" end\n",
" xi = xi + h\n",
" end\n",
" rhs = zeros(dim); rhs[1] = 1; rhs[dim] = 3\n",
" return A, rhs\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us run a particular example, with ``dim`` points."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dim = 3 # the dimension of the problem"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"the matrix A :\n",
"3×3 Matrix{Float64}:\n",
" 1.0000000e+00 0.0000000e+00 0.0000000e+00\n",
" -4.0000000e+00 -2.0000000e+00 1.0000000e+00\n",
" 1.0000000e+00 1.0000000e+00 1.0000000e+00\n",
"the right hand side b :\n",
"3-element Vector{Float64}:\n",
" 1.0000000e+00\n",
" 0.0000000e+00\n",
" 3.0000000e+00\n"
]
}
],
"source": [
"A, b = setup1(dim, 0.0, 1.0)\n",
"println(\"the matrix A :\")\n",
"show(stdout, \"text/plain\", A); println(\"\")\n",
"println(\"the right hand side b :\")\n",
"show(stdout, \"text/plain\", b); println(\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution is obtained by solving a linear system ``A*x = b``, solved with the backslash operator."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 1.0000000e+00\n",
" -6.6666667e-01\n",
" 2.6666667e+00"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cff = A\\b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution ``cff`` represents the coefficients of a polynomial.\n",
"The first element ``cff[1]`` is the constant term in the polynomial.\n",
"That the value of ``cff[1]`` equals ``1.0`` is because of the first boundary condition: ``y(0) = 1.0``."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To compare with the exact solution, we evaluate the polynomial with coefficients in ``cff``, for some step size ``h``."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"approximation : exact value : error\n",
"1.0000000e+00 : 1.0000000e+00 : 0.00e+00\n",
"1.3333333e+00 : 1.2961085e+00 : 3.72e-02\n",
"3.0000000e+00 : 3.0000000e+00 : 4.44e-16\n"
]
}
],
"source": [
"h = 1.0/(dim-1)\n",
"exactval = Array([exactsol((i-1)*h) for i=1:dim])\n",
"println(\"approximation : exact value : error\")\n",
"for i=1:dim\n",
" valsol = evalpoly((i-1)*h, cff)\n",
" show(stdout, \"text/plain\", valsol); print(\" : \")\n",
" show(stdout, \"text/plain\", exactval[i]); print(\" : \")\n",
" err = abs(valsol - exactval[i])\n",
" strerr = @sprintf(\"%.2e\", err)\n",
" println(strerr)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While the errors are acceptable in low dimension, working with high degree polynomials in the standard basis does not lead to accurate results for larger problems."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:108884): GLib-GIO-WARNING **: 19:24:41.403: Unexpectedly, UWP app `Evernote.Evernote_10.48.4.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n",
"\n",
"(process:108884): GLib-GIO-WARNING **: 19:24:41.403: Unexpectedly, UWP app `Clipchamp.Clipchamp_2.5.4.0_neutral__yxz26nhyzhsrt' (AUMId `Clipchamp.Clipchamp_yxz26nhyzhsrt!App') supports 41 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xr = 0.0:0.01:1.0\n",
"yr = [exactsol(t) for t in xr]\n",
"px = [evalpoly(t, cff) for t in xr]\n",
"plot(xr,yr,label=\"y(x)\",legend=:topleft)\n",
"plot!(xr,px,label=\"p(x)\", legend=:topleft)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"figcollocation.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 2. The Galerkin Method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve the same example, but now with an orthogonal basis of spline functions.\n",
"This will result in a tridiagonal system, defined by the function below."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"setup2"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" A, b = setup2(dim::Int)\n",
"\n",
"returns the linear system to solve a BVP\n",
"with the Galerkin method over the interval [a,b]\n",
"with dim points, with a step size h of (b-a)/(dim+1).\n",
"\"\"\"\n",
"function setup2(dim::Int, a::Float64, b::Float64)\n",
" h = (b - a)/(dim+1)\n",
" alpha = 2.0/h + 8.0*h/3.0\n",
" beta = -1.0/h + 4.0*h/6.0\n",
" d = alpha*ones(dim)\n",
" e = beta*ones(dim-1)\n",
" A = diagm(d) + diagm(+1 => e) + diagm(-1 => e)\n",
" rhs = zeros(dim);\n",
" rhs[1] = -beta; rhs[dim] = -3*beta\n",
" return A, rhs\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us solve this for an example, of dimension ``dim``."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dim = 3"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The matrix A :\n",
"3×3 Matrix{Float64}:\n",
" 8.6666667e+00 -3.8333333e+00 0.0000000e+00\n",
" -3.8333333e+00 8.6666667e+00 -3.8333333e+00\n",
" 0.0000000e+00 -3.8333333e+00 8.6666667e+00\n",
"The vector b :\n",
"3-element Vector{Float64}:\n",
" 3.8333333e+00\n",
" 0.0000000e+00\n",
" 1.1500000e+01\n"
]
}
],
"source": [
"A, b = setup2(dim, 0.0, 1.0)\n",
"println(\"The matrix A :\")\n",
"show(stdout, \"text/plain\", A); println(\"\")\n",
"println(\"The vector b :\")\n",
"show(stdout, \"text/plain\", b); println(\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Observe that the matrix ``A`` is tridiagonal."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 1.0109122e+00\n",
" 1.2855407e+00\n",
" 1.8955276e+00"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol = A\\b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution ``sol`` contains the values of the solution at a grid with step size ``h`` defined next."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.5000000e-01"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h = 1.0/(dim+1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us compare the approximations in ``sol`` to the exact solution, defined by the function ``exactsol``."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"approximation : exact value : error\n",
"1.0109122e+00 : 1.0181162e+00 : 7.20e-03\n",
"1.2855407e+00 : 1.2961085e+00 : 1.06e-02\n",
"1.8955276e+00 : 1.9049351e+00 : 9.41e-03\n"
]
}
],
"source": [
"exactval = Array([exactsol(i*h) for i=1:dim])\n",
"println(\"approximation : exact value : error\")\n",
"for i=1:dim\n",
" show(stdout, \"text/plain\", sol[i]); print(\" : \")\n",
" show(stdout, \"text/plain\", exactval[i]); print(\" : \")\n",
" err = abs(sol[i] - exactval[i])\n",
" strerr = @sprintf(\"%.2e\", err)\n",
" println(strerr)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Even with a large step size we obtain two to three decimal places of accuracy."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To plot the approximation, we must define the spline function."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"splinefun"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" splinefun(x::Float64, i::Int, n::Int, a::Float64, b::Float64)\n",
"\n",
"returns the value of the i-th basis function at x,\n",
"on the interval [a,b] divided in n equal subintervals.\n",
"\"\"\"\n",
"function splinefun(x::Float64, i::Int, n::Int, a::Float64, b::Float64)\n",
" h = (b-a)/n\n",
" xi = a + i*h\n",
" if(i == 0)\n",
" if((x >= a) && (x < xi + h))\n",
" return (xi + h - x)/h\n",
" else\n",
" return 0.0\n",
" end\n",
" end\n",
" if(i == n)\n",
" if((x > xi - h) && (x <= xi))\n",
" return (x - xi + h)/h\n",
" else\n",
" return 0.0\n",
" end\n",
" end\n",
" if((x > xi - h) && (x <= xi))\n",
" return (x - xi + h)/h\n",
" end\n",
" if((x >= xi) && (x < xi + h))\n",
" return (xi + h - x)/h\n",
" end\n",
" return 0.0\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To check, let us make the tent plot."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xr = 0.0:0.01:1.0\n",
"y0 = [splinefun(x,0,4,0.0,1.0) for x in xr]\n",
"y1 = [splinefun(x,1,4,0.0,1.0) for x in xr]\n",
"y2 = [splinefun(x,2,4,0.0,1.0) for x in xr]\n",
"y3 = [splinefun(x,3,4,0.0,1.0) for x in xr]\n",
"y4 = [splinefun(x,4,4,0.0,1.0) for x in xr]\n",
"plot(xr,y0,label=\"B0\")\n",
"plot!(xr,y1,label=\"B1\")\n",
"plot!(xr,y2,label=\"B2\")\n",
"plot!(xr,y3,label=\"B3\")\n",
"plot!(xr,y4,label=\"B4\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The piecewise linear B-spline ``S`` is defined using the spline functions."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Sfun"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" Sfun(x::Float64)\n",
"\"\"\"\n",
"function Sfun(x::Float64)\n",
" result = 1.0*splinefun(x,0,4,0.0,1.0)\n",
" for i=1:3\n",
" result = result + sol[i]*splinefun(x,i,4,0.0,1.0)\n",
" end\n",
" result = result + 3.0*splinefun(x,4,4,0.0,1.0)\n",
" return result\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we are ready to plot the solution ``S(x)``, the piecewise linear B-spline."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xr = 0.0:0.01:1.0\n",
"yr = [exactsol(t) for t in xr]\n",
"Sx = [Sfun(t) for t in xr]\n",
"plot(xr,yr,label=\"y(x)\",legend=:topleft)\n",
"plot!(xr,Sx,label=\"S(x)\", legend=:topleft)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"figGalerkin.png\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.0",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}