26.sagews

f = open( DATA+'thermistor-data.dat')

points = []
for line in f:
    l = line.split(" ")
    t = float( l[0] )
    r = float( l[1] )
    points.append( [r, t] )

f.close()

data_plot = list_plot( points )
data_plot.show()

var( "A,B,C" )
eq1 = 1/273.35 == A + B*ln(14.20) + C*ln(14.20)^3
eq2 = 1/295.15 == A + B*ln(5.60) + C*ln(5.60)^3
eq3 = 1/371.15 == A + B*ln(0.467) + C*ln(0.467)^3

sols = solve( [eq1,eq2,eq3], A,B,C, solution_dict=True)
[ [s[A].n(), s[B].n(), s[C].n()] for s in sols ]

[[0.00290584174603719, 0.000277263522107580, 9.00768857828608e-7]]

def temp_model( r ):
    value = 1/(0.00290584174603719 + 0.000277263522107580*ln(r) + 9.00768857828608e-7*ln(r)^3)
    return n(value)

temp_model( .5 )

368.547107879979

model_plot = plot( temp_model, .45, 15, rgbcolor="red" )
model_plot.show()

# Least-squares fit Alg. 4.5.1 from Applied Numerical Methods for Engineers by Schilling and Harris.

# compute the inner product of phi[i] and phi[i] over the dependent data points.
def self_inner_product( f ):
    global points
    sum = 0.0
    for pair in points:
        xi = pair[1]
        n = f( xi )
        sum = sum + n*n
    return sum


var('x')
m = 4  # The degree of the polynomial + 1
phi = [ None for n in range(m) ]  # This array will hold the orthoganal polynomials.
a = [ None for n in range(m) ]
b = [ None for n in range(m) ]
c = [ None for n in range(m) ]
d = [ None for n in range(m) ]

# Step 1
# with k = 0
phi[0] = lambda x: 1.0

d[0] = len(points)
rsum = 0.0
for pair in points:
    rsum = rsum + pair[1]

b[0] = 1.0/d[0] * rsum
phi[1] = lambda x: x - b[0]

# Step 2
# For k=1 to m-1 compute:

for k in range(1, m):

    print "k=%d" % k

    # Compute d_k
    rsum = 0.0
    print "a"
    for pair in points:
        xi = pair[1]
        rsum = rsum + float(phi[k]( xi )^2.0)
    print "b"
    d[k] = rsum                        # <---- The problem occurs here
    print "d_k=%f" % d[k]; sleep(1)

    # Compute b_k
    rsum = 0.0
    print "c"
    for pair in points:
        xi = pair[1]
        rsum = rsum + xi * phi[k]( xi )^2.0

    print "d"
    b[k] = 1.0/d[k] * rsum
    print "b_k=%f" % b[k]; sleep(1)

    # Compute c_k
    rsum = 0.0
    for pair in points:
        xi = pair[1]
        rsum = rsum + xi * phi[k](xi) * phi[k-1](xi)

    c[k] = 1.0/d[k-1] * rsum
    print "c_k=%f" % c[k]; sleep(1)

    # Compute phi_k
    print "\n"; sleep(1)
    phi[k+1] = lambda x: (x - b[k]) * phi[k](x) - c[k] * phi[k-1](x)

k=1 a b d_k=28758.440000 c d b_k=322.072183 c_k=1150.337600 k=2 a

Traceback (most recent call last): d[k] = rsum File "/home/matthew/.sage/sage_notebook/worksheets/admin/14/code/2.py", line 78, in <lambda> phi[k+Integer(1)] = lambda x: (x - b[k]) * phi[k](x) - c[k] * phi[k-Integer(1)](x) TypeError: unsupported operand type(s) for -: 'float' and 'NoneType'

phi[1](x); phi[2](x)

x - 315.190000000000 (x - 322.0721833729507)*(x - 315.190000000000) - 1150.337599999999