Unlike other suggestions, this is short and doesn't use external libraries like numpy. (Not that using other libraries is bad...it's nice not need to, especially for such a simple problem.)

def line_intersection(line1, line2):    xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])    ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1])    def det(a, b):        return a[0] * b[1] - a[1] * b[0]    div = det(xdiff, ydiff)    if div == 0:       raise Exception('lines do not intersect')    d = (det(*line1), det(*line2))    x = det(d, xdiff) / div    y = det(d, ydiff) / div    return x, yprint line_intersection((A, B), (C, D))

And FYI, I would use tuples instead of lists for your points. E.g.

A = (X, Y)

EDIT: Initially there was a typo. That was fixed Sept 2014 thanks to @zidik.

This is simply the Python transliteration of the following formula, where the lines are (a1, a2) and (b1, b2) and the intersection is p. (If the denominator is zero, the lines have no unique intersection.)

Can't stand aside,

So we have linear system:

A₁ * x + B₁ * y = C₁
A₂ * x + B₂ * y = C₂

let's do it with Cramer's rule, so solution can be found in determinants:

x = D_x/D
y = D_y/D

where D is main determinant of the system:

A₁ B₁
A₂ B₂

and D_x and D_y can be found from matricies:

C₁ B₁
C₂ B₂

and

A₁ C₁
A₂ C₂

(notice, as C column consequently substitues the coef. columns of x and y)

So now the python, for clarity for us, to not mess things up let's do mapping between math and python. We will use array L for storing our coefs A, B, C of the line equations and intestead of pretty x, y we'll have [0], [1], but anyway. Thus, what I wrote above will have the following form further in the code:

for D

L1[0] L1[1]
L2[0] L2[1]

for D_x

L1[2] L1[1]
L2[2] L2[1]

for D_y

L1[0] L1[2]
L2[0] L2[2]

Now go for coding:

line - produces coefs A, B, C of line equation by two points provided,
intersection - finds intersection point (if any) of two lines provided by coefs.

from __future__ import division def line(p1, p2):    A = (p1[1] - p2[1])    B = (p2[0] - p1[0])    C = (p1[0]*p2[1] - p2[0]*p1[1])    return A, B, -Cdef intersection(L1, L2):    D  = L1[0] * L2[1] - L1[1] * L2[0]    Dx = L1[2] * L2[1] - L1[1] * L2[2]    Dy = L1[0] * L2[2] - L1[2] * L2[0]    if D != 0:        x = Dx / D        y = Dy / D        return x,y    else:        return False

Usage example:

L1 = line([0,1], [2,3])L2 = line([2,3], [0,4])R = intersection(L1, L2)if R:    print "Intersection detected:", Relse:    print "No single intersection point detected"

Here is a solution using the Shapely library. Shapely is often used for GIS work, but is built to be useful for computational geometry. I changed your inputs from lists to tuples.

Problem

# Given these endpoints#line 1A = (X, Y)B = (X, Y)#line 2C = (X, Y)D = (X, Y)# Compute this:point_of_intersection = (X, Y)

Solution

import shapelyfrom shapely.geometry import LineString, Pointline1 = LineString([A, B])line2 = LineString([C, D])int_pt = line1.intersection(line2)point_of_intersection = int_pt.x, int_pt.yprint(point_of_intersection)