LPSolutionviaLUDecomp.py

#-------------------------------------------------------------------------------
# Name:        Numerical Project Assignment #1
# Purpose:     Given an integer n and a function f, this program will
#              solve the n points finite-difference approximation of the
#              Poisson-Laplace problem using LU decomposition.
#
# Author:      Addison Euhus
#
# Created:     24/03/2014
# Copyright:   (c) Addison Euhus 2014
#-------------------------------------------------------------------------------

from numpy import *
from pylab import *

#The function f in the LP Method
def f(x):
    return 0*x+2

#The true solution, to be used as a comparison function to be plotted along
#with solution, if desired
def truef(x):
    return (-x**2+x)

#The number of steps (integer n) in the LP Method
n = 10

#Variable assignment
h = 1/(n+1.0)
W = linspace(0,1,n+2)

#Generates the right hand side of the matrix equation for the LP problem
Y = h**2*f(W)[1:n+1]

#Constructs the left hand side of the matrix equation for the LP problem
#The lower triangular matrix
L = eye(n) - diag(array([k/(k+1.) for k in range(1,n)]),-1)
#The upper triangular matrix
U = diag(array([(k+1.)/k for k in range(1,n+1)])) - diag(ones(n-1),1)

#Solving for Z in LZ = Y
Z = [Y[0]]
for i in range(1,n):
    Z.append(Y[i]-L[i][i-1]*Z[i-1])

#Solving for X in UX = Z, starting from the end of X to the beginning and
#reversing the order of X upon completion
X = [float(Z[n-1])/U[n-1][n-1]]
for i in range(1,n):
    X.append(float(Z[n-1-i]-U[n-1-i][n-i]*X[i-1])/U[n-1-i][n-1-i])
X.reverse()

#Adjusts solution for endpoints
solution = zeros(n+2)
solution[1:n+1] = X

#Plots the solution
plot(W,solution,'b')

#Plots the true solution
plot(linspace(0,1,500),truef(linspace(0,1,500)),'r')

show()