{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Romberg Integration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The application of the composite Trapezoidal rule to $\\displaystyle \\int_0^1 \\sqrt{1-x^2} dx$ to approximate $\\pi$ requires too many evaluations of $\\sqrt{1-x^2}$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\n", "(process:12248): GLib-GIO-WARNING **: 14:19:57.456: Unexpectedly, UWP app `Clipchamp.Clipchamp_2.5.1.0_neutral__yxz26nhyzhsrt' (AUMId `Clipchamp.Clipchamp_yxz26nhyzhsrt!App') supports 41 extensions but has no verbs\n", "\n", "(process:12248): GLib-GIO-WARNING **: 14:19:57.518: Unexpectedly, UWP app `Evernote.Evernote_10.47.7.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n" ] } ], "source": [ "using Plots" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "f (generic function with 1 method)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x) = sqrt(1-x^2)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r = 0.0:0.01:1.0\n", "plot(f,r,aspect_ratio=1, label=\"sqrt(1-x^2)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We approximate the area of the quarter circle by the application of the composite trapezoidal rule, defined in an adaptive fashion." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "adaptrap" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", " function adaptrap(f::Function,a::Float64,b::Float64,n::Int64)\n", "\n", "Returns a vector of n approximations for the definite integral of\n", "the function f over the interval [a,b], the i-th entry in the vector\n", "uses 2^i function evaluations.\n", "\n", "Example:\n", "\n", " t = adaptrap(cos,0,pi/2,10)\n", "\"\"\"\n", "function adaptrap(f::Function,a::Float64,b::Float64,n::Int64)\n", " t = zeros(n)\n", " h = (b-a) # size of subinterval\n", " m = 1 # number of subintervals\n", " t[1] = (f(a) + f(b))*h/2\n", " for i = 2:n\n", " h = h/2\n", " for j=0:m-1\n", " t[i] = t[i] + f(a+h+j*2*h)\n", " end;\n", " t[i] = t[i-1]/2 + h*t[i]\n", " m = 2*m\n", " end\n", " return t\n", "end" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "20-element Vector{Float64}:\n", " 0.5\n", " 0.6830127018922193\n", " 0.7489272670256102\n", " 0.7724547860892934\n", " 0.7808132594569352\n", " 0.7837756057192828\n", " 0.784824228194921\n", " 0.7851951980991537\n", " 0.7853263957393075\n", " 0.7853727881799138\n", " 0.7853891916347545\n", " 0.785394991352862\n", " 0.7853970419019385\n", " 0.7853977668874246\n", " 0.7853980232097235\n", " 0.7853981138335553\n", " 0.7853981458739578\n", " 0.7853981572019572\n", " 0.7853981612070116\n", " 0.7853981626230124" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = adaptrap(f,0.0,1.0,20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us compute the errors, subtracting all elements in t by pi/4." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "20-element Vector{String}:\n", " \"-2.85e-01\"\n", " \"-1.02e-01\"\n", " \"-3.65e-02\"\n", " \"-1.29e-02\"\n", " \"-4.58e-03\"\n", " \"-1.62e-03\"\n", " \"-5.74e-04\"\n", " \"-2.03e-04\"\n", " \"-7.18e-05\"\n", " \"-2.54e-05\"\n", " \"-8.97e-06\"\n", " \"-3.17e-06\"\n", " \"-1.12e-06\"\n", " \"-3.97e-07\"\n", " \"-1.40e-07\"\n", " \"-4.96e-08\"\n", " \"-1.75e-08\"\n", " \"-6.20e-09\"\n", " \"-2.19e-09\"\n", " \"-7.74e-10\"" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Printf\n", "e = [@sprintf(\"%.2e\",T - pi/4) for T in t]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At one million function evaluations, we have about 10 decimal places correct." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Romberg extrapolation does not apply to this problem because the function $\\sqrt(1-x^2)$ is not sufficiently many times continuous differentiable everywhere in $[0,1]$, because already the first derivative has $\\sqrt{1-x^2}$ in the denominator and thus is not defined at $x = 1$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore, we must reformulate the problem." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "s = 0.01:0.01:sqrt(2.0)/2.0\n", "L(x) = sqrt(2.0)/2.0\n", "D(x) = x\n", "plot(f,r,aspect_ratio=1, label=\"sqrt(1-x^2)\")\n", "plot!(L,s,color=\"red\", label=\"y=sqrt(2)/2\")\n", "plot!(D,s,label=\"y=x\")\n", "plot!([sqrt(2.0)/2.0, sqrt(2.0)/2], [0.01, sqrt(2.0)/2.0],color=\"red\", label=\"x=sqrt(2)/2\")\n", "savefig(\"romberg4pi.png\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The problem is then reduced to $\\displaystyle \\int_0^{\\sqrt{2}/2} \\sqrt{1-x^2} - \\frac{\\sqrt{2}}{2} dx$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "romberg" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", " function romberg(t::Array{Float64,1})\n", "\n", "Applies extrapolation to the approximations in t,\n", "returned by the composite Trapezoidal rule.\n", "\n", "Example:\n", "\n", " t = adaptrap(cos,0,pi/2,6)\n", " et = romberg(t);\n", "\"\"\"\n", "function romberg(t::Array{Float64,1})\n", " n = length(t)\n", " et = zeros(n,n)\n", " et[:,1] = t\n", " for i = 2:n\n", " for j = 2:i\n", " r = 4^(j-1)\n", " et[i,j] = (r*et[i,j-1] - et[i-1,j-1])/(r - 1);\n", " end\n", " end\n", " return et\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before the application of ``romberg`` we must prepare the array of the results of the composite Trapezoidal rule." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The composite trapezium rule :\n", " 2 1.0355339059327372e-01 3.91e-02\n", " 4 1.3249560917971068e-01 1.02e-02\n", " 8 1.4011017603181600e-01 2.59e-03\n", " 16 1.4204903931769053e-01 6.50e-04\n", " 32 1.4253638456870030e-01 1.63e-04\n", " 64 1.4265839556359677e-01 4.07e-05\n" ] } ], "source": [ "f(x) = sqrt(1-x^2) - sqrt(2.0)/2.0\n", "exact = (pi - 2.0)/8.0\n", "n = 6\n", "t = adaptrap(f,0.0,sqrt(2.0)/2.0,n)\n", "println(\"The composite trapezium rule :\")\n", "N = 2\n", "for i=1:n\n", " strnbr = @sprintf(\"%7d\", N)\n", " strerr = @sprintf(\"%.2e\", abs(exact - t[i]))\n", " strapp = @sprintf(\"%.16e\", t[i])\n", " println(\"$strnbr $strapp $strerr\")\n", " N = 2*N\n", "end" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Romberg integration :\n", "1.0355339059327372e-01 3.91e-02\n", "1.4214301537518967e-01 5.56e-04\n", "1.4268205495633965e-01 1.70e-05\n", "1.4269871825008892e-01 3.63e-07\n", "1.4269907778110696e-01 3.92e-09\n", "1.4269908168053008e-01 1.82e-11\n" ] } ], "source": [ "et = romberg(t)\n", "println(\"Romberg integration :\")\n", "for i=1:n\n", " strerr = @sprintf(\"%.2e\", abs(exact - et[i,i]))\n", " strapp = @sprintf(\"%.16e\", et[i,i])\n", " println(\"$strapp $strerr\")\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that with 64 function evaluations we have about 11 correct decimal places. Let us compute the error table." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Errors in the table \n", " 3.91e-02\n", " 1.02e-02 5.56e-04\n", " 2.59e-03 5.07e-05 1.70e-05\n", " 6.50e-04 3.75e-06 6.24e-07 3.63e-07\n", " 1.63e-04 2.49e-07 1.50e-08 5.32e-09 3.92e-09\n", " 4.07e-05 1.58e-08 2.76e-10 4.27e-11 2.20e-11 1.82e-11\n" ] } ], "source": [ "println(\"Errors in the table \")\n", "for i=1:n\n", " for j=1:i\n", " strerr = @sprintf(\"%.2e\", abs(exact - et[i,j]))\n", " print(\" $strerr\")\n", " end\n", " println(\"\")\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The errors in the table shows the progression of the extrapolation." ] } ], "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 }