# 2.1

def get_soln(x0, t0, tmax, deltat, A, spacing=1, x_vals=None):
    if x_vals is None:
        # Pass a dict as x_vals to have that dict update in place
        # This is helpful when computations that might not finish
        # You don't want to do x_vals={} because python will keep the reference to that dict
        x_vals = {}
    x_vals[t0] = x0
    x = x0
    t = t0
    counter = 0
    while t < tmax:
        dxdt = A*x
        new_t = t + deltat
        new_x = x + dxdt * deltat
        t = new_t
        x = new_x
        if (counter%spacing) == 0: # This is to cut down on RAM use
            x_vals[t] = x
        counter += 1
        # Alternatively, if you want to do matrix-multiplication by hand, here's something that should work for 2x2
        # dx1dt = A[0][0]*x[0] + A[0][1]*x[1]
        # dx2dt = A[1][0]*x[0] + A[1][1]*x[1]
        # new_t = t + deltat
        # new_x1 = x[0] + dx1dt * deltat
        # new_x2 = x[1] + dx2dt * deltat
        # t = new_t
        # x = (new_x1, new_x2)
        # if (counter%spacing) == 0:
        #     x_vals[t] = x
        # counter += 1
    x_vals[t] = x
    return x_vals


# 2.2

# I overlooked something when assigning this problem, sorry!
# The matrix A := [[2900994527*I - 6408044856, 10152071064*I - 4456031192], [3222022554*I - 11196078372, 22879134373*I + 6408044856]] / 6445
# diagonalizes as P*D*P^{-1}, with
# D = [1000003*I         0]
#     [        0 3000017*I]
# P = [7*I - 10      8*I]
#     [ 7*I + 4  9*I + 9]
# Setting w(t) := P^{-1}x(t), so that the differential equation x' = Ax becomes w' = Dw,
# we have w_1(t) = w_1(0) * exp(lambda1*t), and similarly for w_2(t).
# This gives w_1(t+deltat) = w_1(t) * (1 + lambda1*deltat + lambda1^2 deltat^2 / 2 + ...)
# and the linear approximation w_1(t) * (1 + lambda1*deltat).
# When lambda1 is purely imaginary and deltat << 1/|lambda1|,
# one can compute using the Taylor series for sqrt(1 + x) around x = 0 that
# |1 + lambda1*deltat| = 1 + |lambda1|^2*deltat^2 + O(lambda1^4*deltat^4).
# This means that, after n steps, Euler's method will give
# |w(t)| ~= |w(0)| * (1 + |lambda1|^2*deltat^2)^n.
# It takes tmax/deltat steps to get from 0 to tmax, so n can be up to that big.
#    (1 + |lambda1|^2*deltat^2)^(tmax/deltat)
#  = (1 + |lambda1|^2*deltat^2)^( (1/(|lambda1|^2*deltat^2)) * tmax*|lambda1|^2*deltat )
# ~= e^(tmax*|lambda1|^2*deltat)
# so in fact one must take deltat to be substantially smaller than 1/|lambda1|^2,
# not the 1/|lambda1| I was expecting.
# The thing I overlooked is at every step the norm of w(t) is biased to increase,
# unlike a more "random" system where sometimes it might increase and sometimes it might decrease.
# In the latter case, one has "square root cancellation", like how if you flip a coin X times,
# the deviation from exactly 50/50 heads/tails will be about size sqrt(X), not size X.
# Here the biggest eigenvalue is about 3*10^6, so tmax*|lambda|^2 ~= 10^14,
# and taking deltat much smaller than 10^(-14) is not realistic.
# There are probably ways to fix this, but that's too big of a detour I think.
# If your code works for smaller systems, e.g. the ones from problem 1, great!
# If you considered, especially after doing problem 2.3, that deltat would need to be small, excellent


# 2.3

# As mentioned above, A = PDP^{-1}, where
# D = [1000003*I         0]
#     [        0 3000017*I]
# P = [7*I - 10      8*I]
#     [ 7*I + 4  9*I + 9]
# The general solution to the differential equation is
# x(t) = c1*exp(lambda1*t)*v1 + c2*exp(lambda2*t)*v2
# where lambda1 = 1000003*I, lambda2 = 3000017*I, v1 = (7*I - 10, 7*I + 4)^T, v2 = (8*I, 9*I + 9)^T (the columns of P)
# and c1, c2 are arbitrary complex constants.
# To make it so that x(0) = (2,1)^T, we need c1*v1 + c2*v2 = (2,1)^T, i.e. (c1,c2)^T = P^{-1}*(2,1)^T = (46/6445*I - 1168/6445, -383/12890*I + 2159/12890)^T
# So x(t) = (46/6445*I - 1168/6445)*exp(1000003*I*t)*(7*I - 10, 7*I + 4)^T + (-383/12890*I + 2159/12890)*exp(3000017*I*t)*(8*I, 9*I + 9)^T
# and substituting t=10 gives x(10) ~= (0.546582360963995 + 0.736469857104057*I, -1.06390905292510 - 0.367953811974633*I).