diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml new file mode 100644 index 0000000..38562be --- /dev/null +++ b/.github/workflows/Documenter.yml @@ -0,0 +1,24 @@ +name: Documenter +on: + push: + branches: + - main + tags: '*' + pull_request: +jobs: + build: + runs-on: ubuntu-latest # [ubuntu-latest, windows-latest, macOS-latest] + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@latest + with: + version: '1.8.5' + - name: Check Git installation + run: git --version + - name: Install dependencies + run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - name: Build and deploy + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key + run: julia --project=docs/ docs/make.jl \ No newline at end of file diff --git a/.github/workflows/Runtests.yml b/.github/workflows/Runtests.yml new file mode 100644 index 0000000..769ca3e --- /dev/null +++ b/.github/workflows/Runtests.yml @@ -0,0 +1,21 @@ +name: Runtests +on: + push: + pull_request: +jobs: + test: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: ['1.8.5'] + julia-arch: [x64] + os: [ubuntu-latest] # [ubuntu-latest, windows-latest, macOS-latest] + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.julia-version }} + - uses: julia-actions/julia-buildpkg@latest + - name: Run tests + uses: julia-actions/julia-runtest@latest + \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..be0b405 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,11 @@ +{ + "spellright.language": [ + "en", + "sl" + ], + "spellright.documentTypes": [ + "markdown", + "latex", + "plaintext" + ] +} \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..4bbcfd1 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2023 jakob + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml index cddf415..9d8569d 100644 --- a/Project.toml +++ b/Project.toml @@ -4,4 +4,6 @@ authors = ["lovc21 "] version = "0.1.0" [deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md new file mode 100644 index 0000000..0126a17 --- /dev/null +++ b/README.md @@ -0,0 +1,28 @@ +# Interpolation using Barycentric Formula + +[![Documenter](https://github.com/lovc21/Interpolation_with_the_barycentric_formula.jl/actions/workflows/Documenter.yml/badge.svg)](https://github.com/lovc21/Interpolation_with_the_barycentric_formula.jl/actions/workflows/Documenter.yml) +[![Runtests](https://github.com/lovc21/Interpolation_with_the_barycentric_formula.jl/actions/workflows/Runtests.yml/badge.svg)](https://github.com/lovc21/Interpolation_with_the_barycentric_formula.jlin/actions/workflows/Runtests.yml) + +In this repository, you can finde the code for Homework 3 of the Numerical Methods course.The code is written in Julia, and the main implementation can be found in the file `src/Interpolation_with_the_barycentric_formula.jl`. The code is tested using the file `test/runtests.jl`, and it is documented using the file docs/make.jl, you can test the code for the specific by running the `Scripts/script.jl` file. + +To run the code, it is necessary to have Julia installed on your computer. Once downloaded, you can run the code for the following equations: e^(-x^2) on \([-1,1]\) , sin(x)/x on \([0,10]\) and |x^2-2x| on \([1,3]\).To run the code, you can simply start the script `Scripts/script.jl`.The script will compute the interpolation for the three equations and will return the reuslts to the precision of 1e-6. + +The code and the Mathematical backround is documented using Documenter.jl. To see the documentation, you can run the docs/make.jl file. The documentation is also available online at [documentation](https://lovc21.github.io/Interpolation_with_the_barycentric_formula.jl/dev/). + +--- +Here are the results of the interpolation for the three equations for the nodes 1 ,3 the number so that the error is less than 1e-6: + +1. e^(-x^2) +

+ +

+ +2. sin(x)/x +

+ +

+ +3. |x^2-2x| +

+ +

\ No newline at end of file diff --git a/Scripts/Ploting.jl b/Scripts/Ploting.jl new file mode 100644 index 0000000..d051fac --- /dev/null +++ b/Scripts/Ploting.jl @@ -0,0 +1,19 @@ + +using Plots + +x1 = -1:0.01:1 +x2 = 0.01:0.1:10 +x3 = 1:0.01:3 + +f1(x) = exp(-x^2) +f2(x) = sin(x) / x +#f2(x) = x == 0 ? 1 : sin(x)/x +f3(x) = abs(x^2 - 2x) + +# Plot the functions +plot(x1, f1.(x1), label="e^{-x^2}", legend=:topright) + +plot(x2, f2.(x2), label="sin(x)/x") + +plot(x3, f3.(x3), label="|x^2-2x|") + diff --git a/Scripts/script.jl b/Scripts/script.jl new file mode 100644 index 0000000..f77b636 --- /dev/null +++ b/Scripts/script.jl @@ -0,0 +1,39 @@ +using Interpolation_with_the_barycentric_formula + +precision = 1e-6 +n = 10 + +f1(x) = exp.(-x.^2) +#f2(x) = sin.(x) ./ x to ne dela +f2(x) = x == 0 ? 1 : sin(x)/x # to dela +f3(x) = abs.(x.^2 - 2 .* x) + +functions = [f1, f2, f3] +intervals = [(-1, 1), (0, 10), (1, 3)] +names = ["exp(-x^2)", "sin(x)/x", "|x^2 - 2x|"] + + +# Calculate and print the optimal degree for each function +for (i, f) in enumerate(functions) + degree = get_optimal_degree(f, intervals[i][1], intervals[i][2], n, precision) + println("Optimal degree for $(names[i]) on [$(intervals[i][1]), $(intervals[i][2])]: $degree") +end + +# Compute the maximum error for each function at the optimal degrees and print the results +degrees = [10, 13, 1567] +for (i, f) in enumerate(functions) + error = compute_max_error(f, intervals[i][1], intervals[i][2], degrees[i]) + println("Maximum error for $(names[i]) on [$(intervals[i][1]), $(intervals[i][2])]: $error") +end + +# Display the plots for each function for increasing n values up to their max degrees +for (i, f) in enumerate(functions) + for n in [1, 3, degrees[i]] + p = plot_interpolation(f, intervals[i][1], intervals[i][2], n, names[i]) + display(p) + sleep(2) + end +end + + + diff --git a/docs/make.jl b/docs/make.jl index 21e35f7..7c503a0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,6 +10,8 @@ makedocs( # Documenter can also automatically deploy documentation to gh-pages. # See "Hosting Documentation" and deploydocs() in the Documenter manual # for more information. -#=deploydocs( - repo = "" -)=# +deploydocs( + repo = "github.com/lovc21/Interpolation_with_the_barycentric_formula.jl.git", + push_preview = true, + devbranch = "main", +) diff --git a/docs/src/images/ezgif-4-53e4254853.gif b/docs/src/images/ezgif-4-53e4254853.gif new file mode 100644 index 0000000..33393ed Binary files /dev/null and b/docs/src/images/ezgif-4-53e4254853.gif differ diff --git a/docs/src/images/ezgif-4-cfcfa1457c.gif b/docs/src/images/ezgif-4-cfcfa1457c.gif new file mode 100644 index 0000000..cad12c1 Binary files /dev/null and b/docs/src/images/ezgif-4-cfcfa1457c.gif differ diff --git a/docs/src/images/ezgif-4-e9a28c1c2b.gif b/docs/src/images/ezgif-4-e9a28c1c2b.gif new file mode 100644 index 0000000..4d27fb3 Binary files /dev/null and b/docs/src/images/ezgif-4-e9a28c1c2b.gif differ diff --git a/docs/src/index.md b/docs/src/index.md index 3498c49..7032f38 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,3 +1,160 @@ -# Interpolation_with_the_barycentric_formula.jl +# Interpolation with the barycentric formula -Documentation for Interpolation_with_the_barycentric_formula.jl +This is the documentation for the repository of homework 3 Interpolation with the barycentric formula of the course Numerical Methods. This documentation is devided in two parts: the first one is the documentation of the mathematical backround of the task and the second part is the documentation of the code we used to solve the task. + +## Mathematical backround + +### Introduction +The point of the given task is to numerically interpolate the following function: e^(-x^2) on \([-1,1]\) , sin(x)/x on \([0,10]\) and |x^2-2x| on \([1,3]\). The interpolation is done with the barycentric formula and Chebyshev points. + +### Chebyshev points +Chebyshev points are the roots of the Chebyshev polynomials. The Chebyshev polynomials are defined as a sequence of orthogonal polynomials defined on the interval [−1,1]. + +Chebyshev points of the first kind in the interval [−1,1] are given by the formula: + + + + + +
+ +
+ +And for the Chebyshev points on the interval [a, b]: + + + + + +
+ +
+ +### Barycentric formula +The barycentric formula is a method for numerical interpolation. It is based on the Lagrange interpolation formula. The Lagrange interpolation formula is the following: + + + + + +
+ +
+ +But this formula is numerically unstable, so we use the barycentric formula instead. The barycentric formula is the following: + + + + + +
+ +
+ +This formula offers a numerically stable way to interpolate a function and it is also faster than the Lagrange interpolation formula. + +### The results of the interpolation +The following pictures show the results of the interpolation of the given functions with the barycentric formula and Chebyshev points. + +1. e^(-x^2) +

+ +

+ +2. sin(x)/x +

+ +

+ +3. |x^2-2x| +

+ +

+ +## Code documentation +First we calculate the Chebyshev nodes using the second formula for the Chebyshev points. We use the following function to calculate the Chebyshev nodes: + +```julia +function chebyshev_nodes(n, a, b) + + # Check the nodes calculation + nodes = [cos(π * (i / n)) for i in 0:n] + + # Check the mapping to the interval + mapped_nodes = [(a + b) / 2 + (b - a) / 2 * x for x in nodes] + + return mapped_nodes +end +``` +Next we calculate the weights for the barycentric formula. We use the following function to calculate the weights: + +```julia +function chebyshev_weights(n) + weights = [Float64((-1)^i) for i in 0:n] + weights[1] *= 0.5 + weights[end] *= 0.5 + return weights +end +``` +After that we can caculate the barycentric interpolation using the barycentric formula. We use the following function to calculate the barycentric interpolation: +```julia +function barycentric_interpolation(nodes,values, weights, x) + if x in nodes + return values[nodes .== x][1] + else + # Calculate the numerator + numerator = sum(values[i] * weights[i] ./ (x .- nodes[i]) for i in 1:length(nodes)) + # Calculate the denominator + denominator = sum(weights[i] ./ (x .- nodes[i]) for i in 1:length(nodes)) + end + + return numerator / denominator +end +``` + +After we have the barycentric interpolation we can run the code so that we are in the exaptible level of precision. We use the following function to run the code: +```julia +function get_optimal_degree(func, a, b,n,precision) + while true + nodes = chebyshev_nodes(n, a, b) + values = [func(node) for node in nodes] + weights = chebyshev_weights(n) + + test_points = range(a, stop=b, length=1000) + max_error = 0.0 + + for x in test_points + interpolated_value = barycentric_interpolation(nodes, values, weights, x) + error = abs(func(x) - interpolated_value) + max_error = max(max_error, error) + end + + if max_error < precision + return n + end + + n += 1 + + end +end +``` + +Finally we can plot the results of the interpolation. We use the following function to plot the results: +```julia +function plot_interpolation(f, a, b, n, name) + nodes = chebyshev_nodes(n, a, b) + values = [f(node) for node in nodes] + weights = chebyshev_weights(n) + + x_vals = range(a, stop=b, length=1000) + y_true = [f(x) for x in x_vals] + y_interp = [barycentric_interpolation(nodes, values, weights, x) for x in x_vals] + + p = plot(x_vals, y_true, label="Original $name", lw=2) + plot!(p, x_vals, y_interp, label="Interpolated $name n=$n", linestyle=:dash, lw=2) + title!("Interpolation of $name with degree $n") + xlabel!("x") + ylabel!("f(x)") + return p +end +``` \ No newline at end of file diff --git a/images/ezgif-4-53e4254853.gif b/images/ezgif-4-53e4254853.gif new file mode 100644 index 0000000..33393ed Binary files /dev/null and b/images/ezgif-4-53e4254853.gif differ diff --git a/images/ezgif-4-cfcfa1457c.gif b/images/ezgif-4-cfcfa1457c.gif new file mode 100644 index 0000000..cad12c1 Binary files /dev/null and b/images/ezgif-4-cfcfa1457c.gif differ diff --git a/images/ezgif-4-e9a28c1c2b.gif b/images/ezgif-4-e9a28c1c2b.gif new file mode 100644 index 0000000..4d27fb3 Binary files /dev/null and b/images/ezgif-4-e9a28c1c2b.gif differ diff --git a/src/Interpolation_with_the_barycentric_formula.jl b/src/Interpolation_with_the_barycentric_formula.jl index d7a0232..c8392b7 100644 --- a/src/Interpolation_with_the_barycentric_formula.jl +++ b/src/Interpolation_with_the_barycentric_formula.jl @@ -1,5 +1,208 @@ module Interpolation_with_the_barycentric_formula -greet() = print("Hello World!") +using LinearAlgebra +using Plots + +export barycentric_interpolation, chebyshev_nodes, chebyshev_weights, get_optimal_degree, compute_max_error, plot_interpolation + +""" +This function calculates the Chebyshev nodes for the given interval. + Parameters: + n: the number of nodes + a: the left endpoint of the interval + b: the right endpoint of the interval + Returns: + mapped_nodes: the Chebyshev nodes +""" +function chebyshev_nodes(n, a, b) + + # Check the nodes calculation + nodes = [cos(π * (i / n)) for i in 0:n] + + # Check the mapping to the interval + mapped_nodes = [(a + b) / 2 + (b - a) / 2 * x for x in nodes] + + return mapped_nodes +end + +"""chebyshev_nodes(10, -1, 1) +11-element Vector{Float64}: + 1.0 + 0.9510565162951535 + 0.8090169943749475 + 0.5877852522924731 + 0.30901699437494745 + 6.123233995736766e-17 + -0.30901699437494734 + -0.587785252292473 + -0.8090169943749473 + -0.9510565162951535 + -1.0 +""" + +""" +This function calculates the barycentric weights for the Chebyshev nodes. + Parameters: + n: the number of nodes + Returns: + weights: the barycentric weights +""" +function chebyshev_weights(n) + weights = [Float64((-1)^i) for i in 0:n] + weights[1] *= 0.5 + weights[end] *= 0.5 + return weights +end + +""" +chebyshev_weights(10) +11-element Vector{Float64}: + 0.5 + -1.0 + 1.0 + -1.0 + 1.0 + -1.0 + 1.0 + -1.0 + 1.0 + -1.0 + 0.5 +""" + +""" +This function calculates the barycentric interpolation for the given function. + Parameters: + nodes: the nodes + values: the values + weights: the weights + x: the point + Returns: + numerator / denominator: the barycentric interpolation +""" +function barycentric_interpolation(nodes,values, weights, x) + if x in nodes + return values[nodes .== x][1] + else + # Calculate the numerator + numerator = sum(values[i] * weights[i] ./ (x .- nodes[i]) for i in 1:length(nodes)) + # Calculate the denominator + denominator = sum(weights[i] ./ (x .- nodes[i]) for i in 1:length(nodes)) + end + + return numerator / denominator +end + +""" +n = 10 +a, b = -1, 1 +nodes = chebyshev_nodes(n, a, b) +values = [exp(-node^2) for node in nodes] +weights = chebyshev_weights(n) +x = 1e-6 +interpolated_value = barycentric_interpolation(nodes, values, weights, x) +0.999999999999 +""" + +""" +This function calculates the optimal degree for the given function. + Parameters: + func: the function + a: the left endpoint of the interval + b: the right endpoint of the interval + n: the number of nodes + precision: the precision + Returns: + n: the optimal degree +""" + +function get_optimal_degree(func, a, b,n,precision) + while true + nodes = chebyshev_nodes(n, a, b) + values = [func(node) for node in nodes] + weights = chebyshev_weights(n) + + test_points = range(a, stop=b, length=1000) + max_error = 0.0 + + for x in test_points + interpolated_value = barycentric_interpolation(nodes, values, weights, x) + error = abs(func(x) - interpolated_value) + max_error = max(max_error, error) + end + + if max_error < precision + return n + end + + n += 1 + + end +end + +""" +n = 1 +precision = 1e-6 + +f1(x) = exp(-x^2) +f2(x) = x == 0 ? 1 : sin(x)/x +f3(x) = abs(x^2 - 2 * x) + +println("e^(-x^2) [-1, 1]: ", get_optimal_degree(f1, -1, 1, n, precision)) +println("sin(x)/x on [0, 10]: ", get_optimal_degree(f2, 0, 10, n, precision)) +println("|x^2 - 2x| on [1, 3]: ", get_optimal_degree(f3, 1, 3, n, precision)) +e^(-x^2) [-1, 1]: 10 +sin(x)/x on [0, 10]: 13 +|x^2 - 2x| on [1, 3]: 1567 +""" + +""" +This function computes the maximum error for the given function. + Parameters: + f: the function + a: the left endpoint of the interval + b: the right endpoint of the interval + n: the number of nodes + Returns: + max_error: the maximum error +""" +function compute_max_error(f, a, b, n) + nodes = chebyshev_nodes(n, a, b) + values = [f(node) for node in nodes] + weights = chebyshev_weights(n) + + x_vals = range(a, stop=b, length=1000) + errors = [abs(f(x) - barycentric_interpolation(nodes, values, weights, x)) for x in x_vals] + + return maximum(errors) +end + +""" +This function plots the original function and the interpolated function. + Parameters: + f: the original function + a: the left endpoint of the interval + b: the right endpoint of the interval + n: the number of nodes + name: the name of the function + Returns: + p: the plot +""" +function plot_interpolation(f, a, b, n, name) + nodes = chebyshev_nodes(n, a, b) + values = [f(node) for node in nodes] + weights = chebyshev_weights(n) + + x_vals = range(a, stop=b, length=1000) + y_true = [f(x) for x in x_vals] + y_interp = [barycentric_interpolation(nodes, values, weights, x) for x in x_vals] + + p = plot(x_vals, y_true, label="Original $name", lw=2) + plot!(p, x_vals, y_interp, label="Interpolated $name n=$n", linestyle=:dash, lw=2) + title!("Interpolation of $name with degree $n") + xlabel!("x") + ylabel!("f(x)") + return p +end end # module Interpolation_with_the_barycentric_formula diff --git a/test/runtests.jl b/test/runtests.jl index 5b84ae8..d01ad4a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,48 @@ using Test +using Interpolation_with_the_barycentric_formula +# Test 1: Chebyshev Nodes +@testset "Test 1: Chebyshev Nodes" begin + result = chebyshev_nodes(10, -1, 1) + expected = [1.0, 0.9510565162951535, 0.8090169943749475, 0.5877852522924731, 0.30901699437494745, 6.123233995736766e-17, -0.30901699437494734, -0.587785252292473, -0.8090169943749473, -0.9510565162951535, -1.0] + @test result ≈ expected +end -@test true +# Test 2: Chebyshev Weights +@testset "Test 2: Chebyshev Weights" begin + result = chebyshev_weights(10) + expected = [0.5, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 0.5] + @test result ≈ expected +end + +# Test 3: Barycentric Interpolation +@testset "Test 3: Barycentric Interpolation" begin + nodes = chebyshev_nodes(10, -1, 1) + values = [exp(-node^2) for node in nodes] + weights = chebyshev_weights(10) + result = barycentric_interpolation(nodes, values, weights, 1e-6) + @test result ≈ 0.999999999999 +end + +# Test 4: Get Optimal Degree +@testset "Test 4: Get Optimal Degree" begin + f1 = x -> exp(-x^2) + result1 = get_optimal_degree(f1, -1, 1, 1, 1e-6) + @test result1 == 10 + + f2 = x -> x == 0 ? 1 : sin(x)/x + result2 = get_optimal_degree(f2, 0, 10, 1, 1e-6) + @test result2 == 13 + + f3 = x -> abs(x^2 - 2 * x) + result3 = get_optimal_degree(f3, 1, 3, 1, 1e-6) + @test result3 == 1567 +end + +# Test 5: Compute Max Error +@testset "Test 5: Compute Max Error" begin + f1 = x -> exp(-x^2) + result1 = compute_max_error(f1, -1, 1, 10) + @test 0 <= result1 <= 1e-6 +end