Step by Step solution about How to resolve the algorithm Runge-Kutta method step by step in the Julia programming language

RK4 Method for Solving Ordinary Differential Equations in Julia

Introduction: This code solves a system of ordinary differential equations (ODEs) using the Runge-Kutta 4th order (RK4) method in the Julia programming language.

ODE Function: The first line defines the ODE function f(x, y) that takes in two variables x and y and returns their product multiplied by the square root of y.

Theoretical Solution Function: The theoric(t) function returns the analytical solution for the ODE at time t.

RK4 Algorithm Lambda Functions: The rk4(f) function implements the RK4 algorithm as a series of nested lambda functions:

1. The outermost lambda function (1st lambda) takes in f, t, y, and δt.
2. The next lambda function (2nd lambda) calculates δy1 based on f(t, y).
3. The subsequent lambda functions (3rd, 4th, and 5th lambda) calculate δy2, δy3, and δy4 respectively, using the midpoint method and the current values of y and t.
4. The deepest lambda function (5th lambda) uses these δy values to calculate y_{n+1}.

RK4 Algorithm Main Function: The δy function assigns the rk4(f) lambda function to a variable, which can be used to calculate the change in y for a given t, y, and δt.

Solving the ODE: The code initializes the initial conditions t₀, y₀, and tmax, and sets the time step δt. It iteratively applies the RK4 algorithm to calculate the solution at each time step, and prints the error at each integer time value.

Additional RK4 Function for General ODEs: The rk4(f, x₀, y₀, x₁, n) function solves a general ODE system of the form dy/dt = f(x, y) using the RK4 method. It returns vectors of x and y values for the solution.

Example Usage: The code also provides an example of using the rk4 function to solve the ODE f(x, y) = x * √(y) with initial conditions x₀ = 0, y₀ = 1, x₁ = 10, and n = 100. It prints the solution and the error at 10 time points.

f(x, y) = x * sqrt(y)
theoric(t) = (t ^ 2 + 4.0) ^ 2 / 16.0

rk4(f) = (t, y, δt) ->  # 1st (result) lambda
         ((δy1) ->      # 2nd lambda
         ((δy2) ->      # 3rd lambda
         ((δy3) ->      # 4th lambda
         ((δy4) -> ( δy1 + 2δy2 + 2δy3 + δy4 ) / 6 # 5th and deepest lambda: calc y_{n+1}
         )(δt * f(t + δt, y + δy3))         # calc δy₄
         )(δt * f(t + δt / 2, y + δy2 / 2)) # calc δy₃
         )(δt * f(t + δt / 2, y + δy1 / 2)) # calc δy₂
         )(δt * f(t, y))                    # calc δy₁

δy = rk4(f)
t₀, δt, tmax = 0.0, 0.1, 10.0
y₀ = 1.0

t, y = t₀, y₀
while t ≤ tmax
    if t ≈ round(t) @printf("y(%4.1f) = %10.6f\terror: %12.6e\n", t, y, abs(y - theoric(t))) end
    y += δy(t, y, δt)
    t += δt
end


function rk4(f::Function, x₀::Float64, y₀::Float64, x₁::Float64, n)
    vx = Vector{Float64}(undef, n + 1)
    vy = Vector{Float64}(undef, n + 1)
    vx[1] = x = x₀
    vy[1] = y = y₀
    h = (x₁ - x₀) / n
    for i in 1:n
        k₁ = h * f(x, y)
        k₂ = h * f(x + 0.5h, y + 0.5k₁)
        k₃ = h * f(x + 0.5h, y + 0.5k₂)
        k₄ = h * f(x + h, y + k₃)
        vx[i + 1] = x = x₀ + i * h
        vy[i + 1] = y = y + (k₁ + 2k₂ + 2k₃ + k₄) / 6
    end
    return vx, vy
end

vx, vy = rk4(f, 0.0, 1.0, 10.0, 100)
for (x, y) in Iterators.take(zip(vx, vy), 10)
    @printf("%4.1f %10.5f %+12.4e\n", x, y, y - theoric(x))
end