From 331c9c7f019ab5489c038d1746ae4d6334035653 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Tue, 26 Mar 2024 13:17:56 +0000 Subject: [PATCH 1/2] fix: CubicSpline's natural boundary conditions --- src/interpolation_caches.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index bf293d7d..b5173feb 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -268,10 +268,10 @@ function QuadraticSpline(u::uType, t; extrapolate = false) where {uType <: Abstr end """ - QuadraticSpline(u, t; extrapolate = false) + CubicSpline(u, t; extrapolate = false) It is a spline interpolation using piecewise cubic polynomials between each pair of data points. Its first and second derivative is also continuous. -Extrapolation extends the last cubic polynomial on each side. +Second derivative on both ends are zero, which are also called "natural" boundary conditions. Extrapolation extends the last cubic polynomial on each side. ## Arguments @@ -303,9 +303,9 @@ function CubicSpline(u::uType, u, t = munge_data(u, t) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) - dl = h[2:(n + 1)] + dl = vcat(h[2:n], zero(eltype(h))) d_tmp = 2 .* (h[1:(n + 1)] .+ h[2:(n + 2)]) - du = h[2:(n + 1)] + du = vcat(zero(eltype(h)), h[3:(n + 1)]) tA = Tridiagonal(dl, d_tmp, du) # zero for element type of d, which we don't know yet @@ -324,9 +324,9 @@ function CubicSpline(u::uType, t; extrapolate = false) where {uType <: AbstractV u, t = munge_data(u, t) n = length(t) - 1 h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0) - dl = h[2:(n + 1)] + dl = vcat(h[2:n], zero(eltype(h))) d_tmp = 2 .* (h[1:(n + 1)] .+ h[2:(n + 2)]) - du = h[2:(n + 1)] + du = vcat(zero(eltype(h)), h[3:(n + 1)]) tA = Tridiagonal(dl, d_tmp, du) d_ = map( i -> i == 1 || i == n + 1 ? zeros(eltype(t), size(u[1])) : From 60ad06a860d5d6b7beb36b8425d0ed2b40b3cb7d Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Tue, 26 Mar 2024 13:18:16 +0000 Subject: [PATCH 2/2] test: update CubicSpline tests --- test/interpolation_tests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index ebfe1947..5be8c566 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -498,8 +498,8 @@ end A = CubicSpline(u, t; extrapolate = true) # Solution - P₁ = x -> 1 + 1.5x + x^2 + 0.5x^3 # for x ∈ [-1.0, 0.0] - P₂ = x -> 1 + 1.5x + x^2 - 0.5x^3 # for x ∈ [0.0, 1.0] + P₁ = x -> 1 + 1.5x + 0.75 * x^2 + 0.25 * x^3 # for x ∈ [-1.0, 0.0] + P₂ = x -> 1 + 1.5x + 0.75 * x^2 - 0.25 * x^3 # for x ∈ [0.0, 1.0] for (_t, _u) in zip(t, u) @test A(_t) == _u @@ -534,8 +534,8 @@ end u = [0.0, 1.0, 3.0] t = [-1.0, 0.0, 1.0] A = CubicSpline(u, t; extrapolate = true) - @test A(-2.0) ≈ -2.0 - @test A(2.0) ≈ 4.0 + @test A(-2.0) ≈ -1.0 + @test A(2.0) ≈ 5.0 A = CubicSpline(u, t) @test_throws DataInterpolations.ExtrapolationError A(-2.0) @test_throws DataInterpolations.ExtrapolationError A(2.0)