Trying to compute the error from comparing two arrays

Question

Some context: I am working with the Black-Scholes model.. I have an explicit (Black-Scholes) formula which is the exact solution to my problem. I have written code which implements a finite-difference scheme (explicit Forward-Euler, i.e: forward march-in-time iterative scheme) to approximate the solution to a partial differential equation (PDE). I have discretized the problem.

I want to compare these two and verify that the error between the exact solution and my approximation to the solution is indeed small. I have stored the values of the exact solution and the approximation in two (NxM) arrays, for each time-step and node in space.

Here is some pseudocode for what I wish to do:

E = 0
for i = 1,...,N:
    E_new = | BS(x_i) - u_i |
    if E_new >= E
        E = E_new
    else:
        E = E
    end

Here is the code I have written to calculate the error:

# Compute the Error
E = np.zeros((N+1, M+1)) # so that E has NxM entries--one for each time-step and node in space
for m in range(1, M+1):
    for j in range(0, N+1):
        E_new = np.absolute(BS[j,m] - u[j,m])
        if (np.any(E_new >= E)): 
            E = E_new
        else:
            E = E

print 'Error = ', E

Here, I am taking N = 100, M = 10,000. If necessary, I can post the segments of my code where I calculate the arrays BS and u. Both are defined as arrays of the same size, and when I plot them, I am receiving an output which agrees with what the solution looks like. So I believe I am on the right track and that the rest of my code is working as it is supposed to.

However, when I run this particular segment in an attempt to calculate the error, I get that the error is roughly 29. I know this is wrong, the error should be small.

My attempt to fix the problem:

I think the problem possibly lies in the line containing if (np.any(E_new >= E)):. I want to check if any entry of the array E_new is greater than or equal to E. If I just run it with if (E_new >= E):, I get this ValueError: "The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()."

I also tried running this: if (np.any([E_new, E])):, but again no luck. I received this error message, "ValueError: setting an array element with a sequence". I have also tried using np.all in place of np.any in each of my attempts and this also does not seem to fix the issue.

Any and all help is very greatly appreciated. Thank you in advance! I am really stuck here and I do not know how to move forward. I don't need a direct answer necessarily, but some insight as to what is causing the problem and a push in the right direction would be immensely helpful. Thanks again!!

Answer 1

You can just do something to the effect of "$\mathtt{np.max(np.abs( BS - u ) )}$". This computes the $L_\infty$ norm of the error, which seems to be what you're after. It'll save you the heartache/headache of writing loops, and it also happens to execute about 100x faster in python, since you're offloading the work to functions written in C. This assumes that "BS" is your exact solution, "u" is your numerical solution, and that they're both stored in numpy arrays.

Answer 2

I would highly suggest that you use the norm functions in numpy via numpy.linalg.norm(BS - u, ord=2) which would give you the $\ell_2$ norm of your error or if you set ord=numpy.infty you can do the $\ell_\infty$.