# L-38 MCS 471 Fri 18 Nov 2022 : heatbackward.jl

# Application of backward differences to a heat diffusion example,
# taken from section 8.1 of Numerical Analysis by Timothy Sauer.

using Printf
using LinearAlgebra

Base.show(io::IO, f::Float64) = @printf(io, "%.3e", f)

"""
    heatbackward(D::Float64,f0::Function,lb::Function,rb::Function,
                 a::Float64,b::Float64,T::Float64,M::Int,N::Int)

returns the M-by-N matrix of approximated values over [a,b]
with time running from 0 to T.

The diffusion coefficient is given by D.

The function f0 is a function in the space variable x
and gives the initial temperature distribution over [a,b].

The function lb is a function in the time variable
and defines the temperature at the left bound of [a,b].

The function rb is a function in the time variable 
and defines the temperature at the right bound of [a,b].

EXAMPLE :
    f(x) = (sin(2*pi*x))^2
    lb(t) = 0
    rb(t) = 0
    sol = heatbackward(1.0,f,lb,rb,0.0,1.0,1.0,10,250)
"""
function heatbackward(D::Float64,f0::Function,lb::Function,rb::Function,
                      a::Float64,b::Float64,T::Float64,M::Int,N::Int)
    result = zeros(M,N)
    offset = zeros(M)
    h = (b-a)/(M-1)
    k = T/(N-1)
    sigma = (D*k)/h^2
    Adiag = (1+2*sigma)*ones(M)
    Asubd = -sigma*ones(M-1)
    A = diagm(Adiag) + diagm(+1 => Asubd) + diagm(-1 => Asubd)
    result[:,1] = Array([f0(i*h) for i=1:M])
    for j=2:N
        offset[1] = sigma*lb((j-1)*k)
        offset[M] = sigma*rb((j-1)*k)
        rhs = result[:,j-1] + offset
        result[:,j] = A\rhs
    end
    return result
end

"""
    Runs a particular heat diffusion example.
"""
function main()
    f(x) = (sin(2*pi*x))^2
    lb(t) = 0
    rb(t) = 0
    sol = heatbackward(1.0,f,lb,rb,0.0,1.0,1.0,6,6)
    println("The solution :")
    show(stdout, "text/plain", sol); println("")
end

main()