using LinearAlgebra, PyPlot

A = [ 0.1 -0.1
      0.5 -1   ]

t = range(0, 20, length=1000)
# find solution by the brute-force eᴬᵗ [1,10]:
x = [exp(A*t)*[1,10] for t in t]
plot(t, [x[1] for x in x], "r-")
plot(t, [x[2] for x in x], "b-")
xlabel(L"time $t$")
ylabel(L"solution $x(t)$")
legend([L"1st $x$ component", L"2nd $x$ component"])
title("solution of 2 coupled ODEs")
grid()

t = range(0, 100, length=1000)
# find solution by the brute-force eᴬᵗ [1,10]:
x = [exp(A*t)*[1,10] for t in t]
plot(t, [x[1] for x in x], "r-")
plot(t, [x[2] for x in x], "b-")
xlabel(L"time $t$")
ylabel(L"solution $x(t)$")
legend([L"1st $x$ component", L"2nd $x$ component"])
title("solution of 2 coupled ODEs")
grid()

λ, X = eigen(A)
λ

X

c = X \ [1, 10]

A = [ 0     1 
     -0.01  0 ]

t = range(0, 20π*5, length=1000)
# find solution by the brute-force eᴬᵗ [0,1]:
plot(t, [(exp(A*t)*[0,1])[1] for t in t], "r-")
xlabel(L"time $t$")
ylabel(L"solution $x(t)$")
title("motion of a frictionless mass on a spring")
grid()

λ, X = eigen(A)
λ

X

c = X \ [0, 1]

2*real(c[1]*X[:,1])

ξ = 2 * c[1] * X[1,1]

# polar form of α:
r = abs(ξ)
ϕ = angle(ξ)
r, ϕ/π

t = range(0, 20π*5, length=1000)
# find solution by the brute-force eᴬᵗ [0,1]:
plot(t, [(exp(A*t)*[0,1])[1] for t in t], "r-")
plot(t, [r*cos(-0.1*t + ϕ) for t in t], "k--")
xlabel(L"time $t$")
ylabel(L"solution $x(t)$")
title("motion of a frictionless mass on a spring")
grid()
legend(["numerical", "analytical"])

2π/0.1

B = [ 0      1 
     -0.01  -0.02 ]

t = range(0, 20π*5, length=1000)
# find solution by the brute-force eᴬᵗ [0,1]:
plot(t, [(exp(B*t)*[0,1])[1] for t in t], "r-")
xlabel(L"time $t$")
ylabel(L"solution $x(t)$")
title("motion of a mass on a spring + drag")
grid()

eigvals(B)

sqrt(0.01 - (0.02/2)^2)

t = range(0, 20π*5, length=1000)
# find solution by the brute-force eᴬᵗ [0,1]:
plot(t, [(exp(B*t)*[0,1])[1] for t in t], "r-")
plot(t, 10*exp.(-0.01 * t), "k--")
xlabel(L"time $t$")
ylabel(L"solution $x(t)$")
title("motion of a mass on a spring + drag")
legend(["solution \$x(t)\$", "exponential decay \$e^{-\\alpha t}\$"])
grid()

eigvals([ 0      1 
       -0.01  -0.3 ])

βs = [0.02, 0.05, 0.1, 0.2, 0.3, 0.4]
t = range(0, 20π*5, length=1000)
for β in βs
    Bd = [ 0      1 
          -0.01  -β ]
    plot(t, [(exp(Bd*t)*[0,1])[1] for t in t], "-")
end
xlabel(L"time $t$")
ylabel(L"solution $x(t)$")
title("mass on a spring for different drags")
legend(["\$\\beta=$β\$" for β in βs])
grid()

using Interact
t = range(0, 20π*5, length=1000)
fig = figure()
@manipulate for β in 0.0:0.01:0.3
    Bd = [ 0      1 
          -0.01  -β ]
    withfig(fig) do
        c = β < 0.1999 ? "blue" : β > 0.2001 ? "red" : "black"
        plot(t, [(exp(Bd*t)*[0,1])[1] for t in t], "-", color=c)
        xlabel(L"time $t$")
        ylabel(L"solution $x(t)$")
        T = c == "blue" ? "underdamped" : c == "red" ? "overdamped" : "critically damped"
        title("$T mass on a spring for drag \$\\beta=$β\$")
        ylim(-4,8)
        grid()
    end
end

C = [ 0     0    1 0
      0     0    0 1
     -0.02  0.01 0 0
      0.01 -0.02 0 0 ]

t = range(0, 20π*5, length=1000)
# find solution by the brute-force eᴬᵗ [0,1]:
x = [(exp(C*t)*[0,0,1,0]) for t in t]
plot(t, [x[1] for x in x], "r-")
plot(t, [x[2]+20 for x in x], "b-")
xlabel(L"time $t$")
ylabel(L"solutions $x(t)$")
title("motion of 2 masses on springs")
legend([L"first mass $x$", L"second mass $x+20$"])
grid()

eigvals(C)

λ, X = eigen(C)
t = range(0, 20π*5, length=1000)
# initial condition from the real part of the first eigenvector:
x = [(exp(C*t)*X[:,1]) for t in t]
plot(t, [real(x[1]) for x in x], "r-")
plot(t, [real(x[2])+2 for x in x], "b-")
xlabel(L"time $t$")
ylabel(L"solutions $x(t)$")
title("first normal mode \$\\omega_1 = 0.1\$")
legend([L"first mass $x$", L"second mass $x+2$"])
grid()

# initial condition from the real part of the eigenvector:
x = [(exp(C*t)*X[:,3]) for t in t]
plot(t, [real(x[1]) for x in x], "r-")
plot(t, [real(x[2])+2 for x in x], "b-")
xlabel(L"time $t$")
ylabel(L"solutions $x(t)$")
title("second normal mode \$\\omega_2 = $(imag(λ[1])) = 0.1\\sqrt{2}\$")
legend([L"first mass $x$", L"second mass $x+2$"])
grid()