How to implement Finite Difference Method ODE Boundary Value Problem in Python?

Question

I implemented the Finite Differences Method for an ODE with Boundary Value Problem. Here is the approximations I used for the FDM:

And here is the balk problem:

with u(0) = u(L) = 0 (attached on both edges)!

I could choose another approximation formula for u'', which has the order 8:

Which approximation formula should I choose? Formula (1) or forumla (4)? Which one is better for this problem? I really don't know why wouldn't formula 4 be more accurate here and can't explain it.

And does anyone know how should I change my actual implementation (using formula (1)), so that it works with formula (4)?

Thank you in advance for any explanation!

Here is my actual implementation, which works as it should (I think at least! If anyone thinks it doesn't please let me know):

import numpy

def fd( x, v, w, t, a, b ):
    """Implements the shooting method to solve linear second order BVPs

    Compute finite difference solution to the BVP

        u'' = x(t) + v(t) u + w(t) u'
        u(t[0]) = a, u(t[n-1]) = b

    t should be passed in as an n element array.   x, v, and w should be
    either n element arrays corresponding to x(t), v(t) and w(t) or
    scalars, in which case an n element array with the given value is
    generated for each of them.

    USAGE:
        u'' = x(t) + v(t) u + w(t) u'
        u = fd(x, v, w, t, a, b)

    INPUT:
        x,v,w - arrays containing x(t), v(t), and w(t) values.  May be
                specified as Python lists, NumPy arrays, or scalars.  In
                each case they are converted to NumPy arrays.
        t     - array of n time values to determine u at
        a     - solution value at the left boundary: a = u(t[0])
        b     - solution value at the right boundary: b = u(t[n-1])

    OUTPUT:
        u     - array of solution function values corresponding to the
                values in the supplied array t.
    """

    # Get the dimension of t and make sure that t is an n-element vector

    if type( t ) != numpy.ndarray:
        if type( t ) == list:
            t = numpy.array( t )
        else:
            t = numpy.array( [ float( t ) ] )

    n = len( t )

    # Make sure that x, v, and w are either scalars or n-element vectors.
    # If they are scalars then we create vectors with the scalar value in
    # each position.

    if type( x ) == int or type( x ) == float:
        x = numpy.array( [ float( x ) ] * n )

    if type( v ) == int or type( v ) == float:
        v = numpy.array( [ float( v ) ] * n )

    if type( w ) == int or type( w ) == float:
        w = numpy.array( [ float( w ) ] * n )

    # Compute the stepsize.  It is assumed that all elements in t are
    # equally spaced.

    h = t[1] - t[0];

    print('h: ', h)

    # Construct tridiagonal system; boundary conditions appear as first and
    # last equations in system.

    A = -( 1.0 + ((w[1:n] * h) / 2.0 ))
    A[-1] = 0.0
    #print('before A: ', A)

    C = -( 1.0 - ((w[0:n-1] * h) / 2.0 ))
    C[0] = 0.0
    #print('before C: ', C)

    print('v: ', v)
    B = 2.0 + (h * h * v)
    B[0] = B[n-1] = 1.0
    #print('before B: ', B)

    print('x: ', x)
    D = - (h * h * x)
    D[0] = a
    D[n-1] = b
    #print('before D: ', D)

    # Solve tridiagonal system

    for i in range( 1, n ):
        umult = A[i-1] / B[i-1]
        D[i] = D[i] - (umult * D[i-1])
        B[i] = B[i] - (umult * C[i-1])
        #print('umult {}: {}'.format(i, umult))
        #print('D {}: {}'.format(i, D))
        #print('B {}: {}'.format(i, B))

    u = numpy.zeros( n )
    u[n-1] = D[n-1] / B[n-1]

    for i in range( n - 2, -1, -1 ):
        u[i] = ( D[i] - (C[i] * u[i+1])) / B[i]
    print('u: ', u)
    return u

#-----------------------------------------------------------------------------#
# Solves u'' = (T/K)*u - (w/2K)* (x^2 - x*L), u(0)=0, u(L) = u(9) = 0 using the
# finite difference method.
#-----------------------------------------------------------------------------#

# Set up constants

T = 500
K = 5 * pow(10,9)
W = 100
L = 9

# Set up interval.  We will solve the problem for both n=64 and n=128.

a = 0.0
b = L
n1 = 4
#n2 = 20
t1 = linspace( a, b, n1 )
#t2 = linspace( a, b, n2 )


# Compute finite difference solutions
#fd( x, v, w, t, a, b )

ufd1 = fd( (W/(2*K)) * (pow(t1,2) - (L*t1)), (T/K), 0, t1, 0, 0)
#ufd2 = fd( (W/(2*K)) * (pow(t2,2) - (L*t2)), (T/K), 0, t2, 0, 0)


# Prepare for display; set interactive mode to true so each plot
# is shown as it is generated

# Plot solutions

plot( t1, ufd1, 'ro')
title( 'Finite Difference Method' )
xlabel( '$t$' )
ylabel( '$x$' )
legend( ( '%3d points' % n1, '%3d points' % n2 ), loc='lower right' )
draw()

Answer 1

Higher order finite difference schemes come with greater accuracy, but also with stricter assumptions on the regularity of your solution $u$. So what is better is hard to say in advance, since (unless you perform a detailed study) it is hard to tell what kind of regularity your solution suffices.

In your code you would then have to adapt the part where you set the matrix $A$ since it is not tridiagonal any longer, but has 3 sub-diagonals. As your code is currently written, it is only working for linear problems of the kind $u'' = x(t) + v(t) u + w(t) u'$ - the problem formulation you give (3), however, is clearly non-linear in $u$ and its derivatives.