From 9e5a1b8cb12264098aae550621ab585700211f15 Mon Sep 17 00:00:00 2001
From: Evan <60061381+E-W-Jones@users.noreply.github.com>
Date: Sat, 4 Nov 2023 04:29:36 +0000
Subject: [PATCH] feat: Add Runge-Kutta Integration (#216)

Add a 4th order Runge-Kutta integration function, and some corresponding tests.

Co-authored-by: Soc Virnyl S. Estela <socvirnyl.estela+github@uncomfyhalomacro.pl>
---
 src/math/Math.jl                    |  2 +
 src/math/runge_kutta_integration.jl | 83 +++++++++++++++++++++++++++++
 test/math.jl                        |  9 ++++
 3 files changed, 94 insertions(+)
 create mode 100644 src/math/runge_kutta_integration.jl

diff --git a/src/math/Math.jl b/src/math/Math.jl
index ef9fd2c9..f8c1eb0a 100644
--- a/src/math/Math.jl
+++ b/src/math/Math.jl
@@ -63,6 +63,7 @@ export permutation
 export prime_check
 export prime_factors
 export riemann_integration
+export runge_kutta_integration
 export simpsons_integration
 export sum_ap
 export sum_gp
@@ -115,6 +116,7 @@ include("permutation.jl")
 include("prime_check.jl")
 include("prime_factors.jl")
 include("riemann_integration.jl")
+include("runge_kutta_integration.jl")
 include("sieve_of_eratosthenes.jl")
 include("simpsons_integration.jl")
 include("sum_of_arithmetic_series.jl")
diff --git a/src/math/runge_kutta_integration.jl b/src/math/runge_kutta_integration.jl
new file mode 100644
index 00000000..c8e2fd36
--- /dev/null
+++ b/src/math/runge_kutta_integration.jl
@@ -0,0 +1,83 @@
+"""
+    runge_kutta_integration(f::Function, x0::Real, y0::Real, h::Real, x_stop::Real)
+
+Numerically solve initial value problems of the form ``y' = f(x, y)`` to find ``y(x)`` using the Runge-Kutta 4th order integration scheme.
+
+Starting with the differential equation ``\\frac{\\mathrm{d}y}{\\mathrm{d}x} = y' = f(x, y)`` and the initial condition ``y(x_0) = y_0``, each step calculates 4 approximations of the gradient
+```math
+\\begin{align*}
+    k_1 &= f(x_n, y_n),\\\\
+    k_2 &= f(x_n + \\frac{h}{2}, y_n + k_1\\frac{h}{2}),\\\\
+    k_3 &= f(x_n + \\frac{h}{2}, y_n + k_2\\frac{h}{2}),\\\\
+    k_4 &= f(x_n + h, y_n + k_3h),
+\\end{align*}
+```
+and uses the weighted average of them,
+
+```math
+\\bar{k} = \\frac{k_1 + 2k_2 + 2k_3 + k_4}{6},
+```
+
+to find the approximate value of ``y(x_{n+1})`` and update ``x`` and ``y`` accordingly
+
+```math
+\\begin{align*}
+    x &\\rightarrow x_{n+1} = x_n + h\\\\
+    y &\\rightarrow y_{n+1} = y_n + h\\bar{k}.
+\\end{align*}
+```
+
+# Arguments:
+- `f`: The function ``y' = f(x, y)`` to solve for ``y(x)``.
+- `x0`: The starting value of x.
+- `y0`: The starting value of y.
+- `h`: The step size, should be >0.
+- `x_stop`: The final value of x to solve to, i.e. integrate over the range `[x0, x_stop]`
+
+# Examples
+```julia
+julia> # If you have a constant slope of y' = 1, the analytic solution is y = x
+julia> runge_kutta_integration((x, y)->1, 0, 0, 1, 3)
+([0.0, 1.0, 2.0, 3.0], [0.0, 1.0, 2.0, 3.0])
+
+julia> # Consider solving y' = exp(x), which has the analytic solution y = exp(x).
+julia> runge_kutta_integration((x, y)->exp(x), 0, 1., 0.1, 0.1)
+([0.0, 0.1], [1.0, 1.105170921726329])
+
+julia> exp.([0.0, 0.1])
+2-element Vector{Float64}:
+ 1.0
+ 1.1051709180756477
+
+```
+
+# References
+ - [https://en.wikipedia.org/wiki/Initial_value_problem](https://en.wikipedia.org/wiki/Initial_value_problem)
+ - [https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods](https://en.wikipedia.org/wiki/Initial_value_problem)
+
+# Contributors:
+- [E-W-Jones](https://github.com/E-W-Jones)
+"""
+function runge_kutta_integration(f::Function, x0::Real, y0::Real, h::Real, x_stop::Real)
+    h > 0 || throw(DomainError(h, "The step size `h` should be >0."))
+
+    x = Float64(x0)
+    y = Float64(y0)
+    output_x = [x]
+    output_y = [y]
+    
+    while x < x_stop
+        k1 = f(x, y)
+        k2 = f(x + h/2, y + k1*h/2)
+        k3 = f(x + h/2, y + k2*h/2)
+        k4 = f(x + h  , y + k3*h  )
+
+        y += h * (k1 + 2k2 + 2k3 + k4) / 6
+        x += h
+
+        push!(output_x, x)
+        push!(output_y, y)
+    end
+
+    return output_x, output_y
+end
\ No newline at end of file
diff --git a/test/math.jl b/test/math.jl
index 03036754..931d5784 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -377,6 +377,15 @@ using TheAlgorithms.Math
         )
     end
 
+    @testset "Math: Runge_Kutta Integration" begin
+        @test runge_kutta_integration((x, y)->1, 0, 0, 1, 3) == ([0.0, 1.0, 2.0, 3.0], [0.0, 1.0, 2.0, 3.0])
+        @test begin
+            x, y = runge_kutta_integration((x, y)->cos(x), 0, 0, 1e-4, π/2)
+            isapprox(x[end], π/2; rtol=1e-4) && isapprox(y[end], 1; rtol=1e-4)
+        end
+        @test_throws DomainError runge_kutta_integration((x, y)->(), 0, 0, 0, 0)
+    end
+
     @testset "Math: Sum of Arithmetic progression" begin
         @test sum_ap(1, 1, 10) == 55.0
         @test sum_ap(1, 10, 100) == 49600.0