Python. Recursive matrix product with result as a set of vectors

Question

As an application of Eulers method, I'm trying to implement a code which would compute the recursive matrix product Yn = Yn-1 + A(Yn-1), where Y is a vector and A is a matrix such that the product is defined. This is the current code I have

def f(A, y):
    return A.dot(y)


def euler(f, t0, y0, T, dt):
    t = np.arange(t0, T + dt, dt)
    y = [0,0,0,0]*len(t)
    y[0] = y0
    for i in range(1, len(t)):
        y[i] = y[i - 1] + f(A, y[i - 1])*dt
        return t, y


# Define problem specific values


A = np.array([[0, 0, 1, 0],
              [0, 0, 0, 1],
              [-2, -3, 0, 0],
              [-3, -2, 0, 0]])

y1_0 = 1
y2_0 = 2
y3_0 = 0
y4_0 = 0
y0 = [y1_0, y2_0, y3_0, y4_0]

t,y = euler(f,0,y0,2,1)

print(t,y)

For example, the result for points in the range t0 = 0, T = 2 should be the vectors Y1 and Y2. Instead I have

[0 1 2] [[1, 2, 0, 0], array([ 1, 2, -8, -7]), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Something is wrong here. While Y1 = [1, 2, -8, -7 ] does show up, there is all of this unnecessary stuff. And Y2 is not printed at all. I suspect this is due to how I define the variable y. For every point in the range of t, I need a vector of 4 zeros - which is then filled up by the function euler, I think. How should correct this?

Answer 1

The computer always does what you tell it to do. In your case y is constructed by repeating 4 zeros len(t) times, giving a list of 12 zeros. The first list entry is replaced by the list y0. The second list entry is replaced by the result of the numpy operations which is a numpy.array. Then the return statement at the level of the loop instructions breaks the loop and returns the t and y arrays. y still contains 10 zeros from its construction that were not replaced.

So construct

    y = np.zeros([len(t), len(y0)])

and repair the indentation level.