Check analytical solution for FiPy Robin BC example

gistfile1.py

c, P, D, r, C1, C2 = var('c P D r C1 C2')
# Characteristic equation to solve ODE
f1(r) = r**2 - P*r- D
# Solution for roots
print solve(f1(r) == 0, r)
# (answer)
r1 = 0.5*(P - sqrt(P^2 + 4*D))
r2 = 0.5*(P + sqrt(P^2 + 4*D))
# Function and first derivative
csol(x, C1, C2) = C1*exp(r1*x) + C2*exp(r2*x)
dcsoldx(x, C1, C2) = derivative(csol, x)
# Solve for constants given BC's
C1C2 = solve(
    [P == -dcsoldx(0, C1, C2) + P*csol(0, C1, C2),
     dcsoldx(1, C1, C2) == 0], C1, C2)
print C1C2[0][0]
print
print C1C2[0][1]
A = sqrt(P^2 + 4*D)
BA = P^4 + 4*D*P^2 - D^2*(e^(-A) + e^(A) - 2)
C1tA = P/2*(P^3 + (e^(A) + 3)*P*D + A*(P^2 - D*(e^(A) - 1)))
C2tA = P/2*(P^3 + (e^(-A) + 3)*P*D - A*(P^2 - D*(e^(-A) - 1)))

csln(x, P, D) = (C1tA*exp(r1*x) + C2tA*exp(r2*x))/BA
dcslndx(x, P, D) = derivative(csln(x, P, D), x)

Pt = 3.0
Dt = 2.0
At = sqrt(Pt + 4*Dt^2)
xt = 0.4
# Not same as original version!
print "new solution test point:", csln(xt, Pt, Dt)
print "original solution test point:", (2 * Pt * exp(Pt * xt / 2) * ((Pt + At) * exp(At / 2 * (1 - xt))
             - (Pt - At) * exp(-At / 2 *(1 - xt)))/
             ((Pt + At)**2*exp(At / 2)- (Pt - At)**2 * exp(-At / 2)))

print "New solution boundary conditions (LHS, and RHS, should be equal):"
print Pt, -dcslndx(0, Pt, Dt) + Pt*csln(0, Pt, Dt)
print 0, dcslndx(1, Pt, Dt)