{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simple example of importing a function on a non uniform grid\n",
"\n",
"ApproxFun generally needs, or needs to be able to evaluate, function values on the Chebyshev grid. There are other methods, but here is one to use when you have data that is not on a Chebyshev grid. It uses [The AAA algorithm for rational approximation](http://arxiv.org/abs/1612.00337) (from the registered package **BaryRational**) to first form a rational approximation to the data. Then it uses that approximation to evaluate on a Chebyshev grid and then transforms those values to Chebyshev expansion coefficients. Finally, it uses those coefficients to create an ApproxFun Fun.\n",
"\n",
"This example first constructs the Fun from a random grid and plots it against the true value. Next the first and second derivatives are now calculated by using ApproxFun operators and plotted against their true values. Finally, a more complex function of all three Funs is constructed $G(x) = (f(x) + 4*f'(x))^2 + sin(f''(x))$ and then integrated and plotted.\n",
"\n",
"Note that we create a new project called ApproxAAA and minimally populate it so as to not depend on you setting up the environment correctly or having this notebook mess with your default environment. When no longer wanted, remember to delete the ApproxAAA directory.\n",
"\n",
"#### Credits\n",
"Initital example (with Splines) by [Juan Ignacio Polanco](https://jipolanco.gitlab.io/) from the [Julia Discourse](https://discourse.julialang.org/t/best-way-to-take-derivatives-of-unevenly-spaced-data-with-interpolations-discrete-derivatives/54097)\n",
"\n",
"We also changed the domain from [0, 1] to [-1, 1] (to match the ApproxFun defaults) and sampled 25 random points instead of 20 (but that should make it a more challenging problem)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m environment at `~/ghub/ApproxFunExamples/Extras/ApproxAAA/Project.toml`\n",
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `~/ghub/ApproxFunExamples/Extras/ApproxAAA/Project.toml`\n",
"\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `~/ghub/ApproxFunExamples/Extras/ApproxAAA/Manifest.toml`\n",
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `~/ghub/ApproxFunExamples/Extras/ApproxAAA/Project.toml`\n",
"\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `~/ghub/ApproxFunExamples/Extras/ApproxAAA/Manifest.toml`\n",
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `~/ghub/ApproxFunExamples/Extras/ApproxAAA/Project.toml`\n",
"\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `~/ghub/ApproxFunExamples/Extras/ApproxAAA/Manifest.toml`\n"
]
}
],
"source": [
"using Pkg\n",
"Pkg.activate(\"ApproxAAA\")\n",
"Pkg.add(\"ApproxFun\")\n",
"Pkg.add(\"BaryRational\")\n",
"Pkg.add(\"Plots\")\n",
"\n",
"using ApproxFun\n",
"using BaryRational\n",
"using Plots\n",
"using Random"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Actual function and its derivatives\n",
"f(x) = sin(2π * x)\n",
"f′(x) = 2π * cos(2π * x)\n",
"f′′(x) = -4π^2 * sin(2π * x)\n",
"\n",
"# Generate uneven grid of ``N`` points in ``[-1, 1]``\n",
"N = 25\n",
"rng = MersenneTwister(42)\n",
"x = sort!(2rand(rng, N) .- 1)\n",
"x[1] = -1; x[end] = 1; # make sure endpoints are included\n",
"y = f.(x) # actual function at the random grid points\n",
"\n",
"# Now create the rational approximation\n",
"modl = aaa(x, y)\n",
"\n",
"M = 129\n",
"pts = chebyshevpoints(M)\n",
"S = Chebyshev()\n",
"\n",
"# The Fun and its first two derivatives\n",
"nfun = Fun(S, ApproxFun.transform(S, modl.(pts)))\n",
"D′ = nfun'\n",
"D′′ = nfun''\n",
"\n",
"xplot = -1:0.01:1\n",
"plt = plot(f, 0, 1; label = \"f(x) exact\", xlabel = \"x\", title = \"Function interpolation\")\n",
"plot!(plt, xplot, nfun.(xplot); label = \"f(x) Interpolation\")\n",
"scatter!(plt, x, y; label = \"Input data\")"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plt = plot(f′, -1, 1; label = \"f'(x) exact\", xlabel = \"x\", legend = :top, title = \"1st derivative approximation\")\n",
"plot!(plt, D′, -1, 1; label = \"f'(x) approx\")\n",
"scatter!(plt, x, D′.(x); label = \"Approximation at grid points\")"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plt = plot(f′′, -1, 1; label = \"f''(x) exact\", xlabel = \"x\", legend = :top, title = \"2nd derivative approximation\")\n",
"plot!(plt, xplot, D′′.(xplot); label = \"f''(x) approx\")\n",
"scatter!(plt, x, D′′.(x); label = \"Approximation at grid points\")"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# A 'complicated' function of f(x), f'(x), and f''(x) from the random grid approximation\n",
"G = (nfun + 4*D′)^2 + sin(D′′)\n",
"\n",
"# and its integral\n",
"Gint = cumsum(G)\n",
"\n",
"# Now compare to the 'exact' integral\n",
"fz = sin(2pi * Fun())\n",
"GE = (fz + 4*fz')^2 + sin(fz'')\n",
"GEint = cumsum(GE)\n",
"\n",
"plt = plot(xplot, Gint.(xplot); label = \"Integral approximation\", legend = :top, title=\"Integral Approximation of (f(x) + 4f'(x))^2 + sin(f''(x))\")\n",
"plot!(plt, xplot, GEint.(xplot); label = \"Exact Integral\")\n",
"scatter!(plt, x, Gint.(x); label = \"Approximation at grid points\")"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.4179153218719304e-6"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# We know that \"converged to the eye\" is not good enough so check the error norm too\n",
"norm(GEint - Gint)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So not too bad after all."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.6.0-beta1",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}