Merge pull request #24 from JuliaTurkuDataScience/2nd-order-example

2nd order example

examples/SimpleHarmonicMotion.jl

@@ -4,25 +4,25 @@ using Plots
 ## inputs
 tSpan = [0, 10]     # [intial time, final time]
 x0 = [1, 1]             # intial value
-β = 1            # order of the derivative
+β = 2            # order of the derivative
 par = [16.0, 4.0] # [spring constant for a mass on a spring, inertial mass]
-h = 0.1
+h = 0.01
 
 ## Equation
 function F(t, n, β, x, par)
 
       K = par[1]
       m = par[2]
 
-      K ./ m .* x[n]
+      - K ./ m .* x[n]
 
 end
 
 ## Numerical solution
 t, Yapp = FDEsolver(F, tSpan, x0, β, nothing, par, h = h)
 
 #plot
-plot(t, Yapp, linewidth = 5, title = "Solution of a 1D with order > 1",
+plot(t, Yapp, linewidth = 5, title = "Simple harmonic motion (order=2)",
      xaxis = "Time (t)", yaxis = "x(t)", label = "Approximation")
 a = x0[1] .* map(cos, sqrt(par[1] / par[2]) .* t) .+ x0[2] ./ sqrt(par[1] / par[2]) .* map(sin, sqrt(par[1] / par[2]) .* t)
 plot!(t, a, lw = 3, ls = :dash, label = "Exact solution")

src/SupFuns.jl

@@ -100,23 +100,3 @@ end
 
 ## Gamma function for vectors ##
 Γ(b) = map(gamma, b)
-
-a_n0(n, β) = ((n - 1) .^ (β .+ 1) .- n .^ β .* (n .- β .- 1)) ./ Γ(β .+ 2)
-
-## function α ##
-function α(n, β)
-
-    if n == 0
-
-        return 1 ./ Γ(β .+ 2)
-
-    else
-
-        return ((n .- 1) .^ (β .+ 1) .- 2 .* n .^ (β .+ 1) .+ (n .+ 1) .^ (β .+ 1)) ./ Γ(β .+ 2)
-
-    end
-
-end
-
-## Gamma function for vectors ##
-Γ(b) = map(gamma, b)

src/main.jl

@@ -9,7 +9,7 @@ function FDEsolver(F, tSpan, y0, β, ::Nothing, par...; h = 0.01, nc = 3, StopIt
 
     # Calculate T with taylor expansion
     m::Int64 = ceil(β[1])
-
+    
     for n in 1:N
 
         T0 = taylor_expansion(tSpan[1], t[n], y0, β, m)

test/SimpleHarmonicMotion_test.jl

@@ -0,0 +1,30 @@
+using FdeSolver
+using Test
+using Statistics
+
+@testset "FdeSolver.jl" begin
+
+tSpan = [0, 10]     # [intial time, final time]
+x0 = [1, 1]             # intial value
+β = 2            # order of the derivative
+par = [16.0, 4.0] # [spring constant for a mass on a spring, inertial mass]
+h = 0.01
+
+function F(t, n, β, x, par)
+
+      K = par[1]
+      m = par[2]
+
+      - K ./ m .* x[n]
+
+end
+
+t, Yapp = FDEsolver(F, tSpan, x0, β, nothing, par, h = h)
+
+a = x0[1] .* map(cos, sqrt(par[1] / par[2]) .* t) .+ x0[2] ./ sqrt(par[1] / par[2]) .* map(sin, sqrt(par[1] / par[2]) .* t)
+
+@test @isdefined(t)
+@test @isdefined(Yapp)
+@test abs((mean(Yapp) - mean(a))) < 0.05
+
+end

test/runtests.jl

@@ -6,4 +6,6 @@ using SafeTestsets
 
 @safetestset "solving 3 parametric FDEs with inner function" begin include("3FDE_par_test.jl") end
 
+@safetestset "solving  Simple Harmonic Motion (order=2) " begin include("SimpleHarmonicMotion_test.jl") end
+
 @safetestset "sample test" begin include("sample_test.jl") end