using LinearAlgebra 
using Plots 
using IterativeSolvers

function func(y,t,ϵ)
    return -(1/ϵ)*(y-sin(t)); 
end

function actual_solution(t,ϵ)
    return (sin(t) - ϵ*cos(t))/(1+ϵ^2) + (1 + ϵ/(1+ϵ^2))*exp(-t/ϵ);  
end

function compute_solution(t,ϵ)
    n = length(t);
    y = zeros(n);
    for i in 1:n
       y[i] = actual_solution(t[i],ϵ); 
    end
    return y
end

function explicit_euler(t,y0,ϵ)
    n = length(t);
    y = zeros(n);
    y[1] = y0;
    dt = t[2]-t[1];
    for i in 2:n
       y[i] = y[i-1] + dt*func(y[i-1],t[i-1],ϵ);  
    end
    return y;
end

function sor_implicit_euler(t,y0,ϵ) 
    n = length(t);
    y = zeros(n);
    y[1] = y0;
    dt = t[2]-t[1];
    res = 100;
    τ = 1e-6;
    iter = 1e6; 
    κ2 = (dt/ϵ); 
    κ1 = 1/(1 + κ2);
    κ3 = 1/κ1;
    it = 0;
    while (it<iter && res>τ)
        for i in 2:n 
           y[i] = κ1*(y[i-1] + κ2*sin(t[i])); 
        end
        for i in 2:n 
           res+= κ2*sin(t[i]) + y[i-1] - κ3*y[i];
        end
        res = sqrt(abs(res/n^2));
        it += 1;
    end
    return y;
end 

function sor_trapezoidal_method(t,y0,ϵ)
    n = length(t);
    y = zeros(n);
    y[1] = y0;
    dt = t[2]-t[1];
    res = 100;
    τ = 1e-6;
    iter = 1e6; 
    κ1 = (dt/(2*ϵ)); 
    κ2 = 1 - κ1;
    κ3 = 1 + κ1;
    κ4 = 1/κ3;
    it = 0;
    while (it<iter && res>τ)
        for i in 2:n 
            y[i] = κ4*(κ1*(sin(t[i]) + sin(t[i-1])) + κ2*y[i-1]);
        end 
        for i in 2:n 
            res += κ1*(sin(t[i]) + sin(t[i-1])) + κ2*y[i-1] - κ3*y[i];
        end
        res = sqrt(abs(res/n^2));
        it += 1;
    end
    return y;
end

function norm(y1,y2) 
    n = length(y1);
    norm = 0;
    for i = 1:n 
        norm += abs(y1[i]-y2[i]);
    end
    norm = sqrt(norm/n*n);
    println("The norm is ",norm)
end

x0 = 0;
xl = pi; 
N = 99;
y0 = 1;
x = LinRange(x0,xl,N)

ϵ = 1;
y_a = compute_solution(x,ϵ);
y = sor_trapezoidal_method(x,y0,ϵ);
plot(x,y)
p1 = plot!(x,y_a,title="Symplectic ϵ=1")
norm(y,y_a);

ϵ = 0.1;
y_a = compute_solution(x,ϵ);
y = sor_trapezoidal_method(x,y0,ϵ);
plot(x,y)
p2 = plot!(x,y_a,title="Symplectic  ϵ=0.1")
norm(y,y_a);

ϵ = 0.01;
y_a = compute_solution(x,ϵ);
y = sor_trapezoidal_method(x,y0,ϵ);
plot(x,y)
p3 = plot!(x,y_a,title="Symplectic  ϵ=0.01")
norm(y,y_a);

ϵ = 0.001;
y_a = compute_solution(x,ϵ);
y = sor_trapezoidal_method(x,y0,ϵ);
plot(x,y)
p4 = plot!(x,y_a,title="Symplectic  ϵ=0.001")
norm(y,y_a);

plot(p1,p2,p3,p4,layout=(2,2),legend=false)

x0 = 0;
xl = pi; 
N = 99;
y0 = 1;
x = LinRange(x0,xl,N)

ϵ = 1;
y_a = compute_solution(x,ϵ);
y = explicit_euler(x,y0,ϵ);
plot(x,y)
p1 = plot!(x,y_a,title="Explcit Euler ϵ=1")
norm(y,y_a);

ϵ = .1;
y_a = compute_solution(x,ϵ);
y = explicit_euler(x,y0,ϵ);
plot(x,y)
p2 = plot!(x,y_a,title="Explcit Euler ϵ=0.1")
norm(y,y_a);

ϵ = .01;
y_a = compute_solution(x,ϵ);
y = explicit_euler(x,y0,ϵ);
plot(x,y)
p3 = plot!(x,y_a,title="Explcit Euler ϵ=0.01")
norm(y,y_a);

ϵ = .001;
y_a = compute_solution(x,ϵ);
y = explicit_euler(x,y0,ϵ);
plot(x,y)
p4 = plot!(x,y_a,title="Explcit Euler ϵ=0.001")
norm(y,y_a);

plot(p1,p2,p3,p4,layout=(2,2),legend=false)

x0 = 0;
xl = pi; 
N = 99;
y0 = 1;
x = LinRange(x0,xl,N);

ϵ = 1;
y_a = compute_solution(x,ϵ);
y = sor_implicit_euler(x,y0,ϵ);
plot(x,y)
p1 = plot!(x,y_a,title="Implicit Euler ϵ=1")
norm(y,y_a);

ϵ = 0.1;
y_a = compute_solution(x,ϵ);
y = sor_implicit_euler(x,y0,ϵ);
plot(x,y)
p2 = plot!(x,y_a,title="Implicit Euler ϵ=0.1")
norm(y,y_a);

ϵ = 0.01;
y_a = compute_solution(x,ϵ);
y = sor_implicit_euler(x,y0,ϵ);
plot(x,y)
p3 = plot!(x,y_a,title="Implicit Euler ϵ=0.01")
norm(y,y_a);

ϵ = 0.001;
y_a = compute_solution(x,ϵ);
y = sor_implicit_euler(x,y0,ϵ);
plot(x,y)
p4 = plot!(x,y_a,title="Implicit Euler ϵ=0.001")
norm(y,y_a);

plot(p1,p2,p3,p4,layout=(2,2),legend=false)


