Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- def read_inputs(filename):
- with open(filename) as f:
- x, y = [], []
- for line in f.readlines():
- coord = [float(n) for n in line.strip().split()]
- x.append(coord[0])
- y.append(coord[1])
- return x, y
- def read_from_file(file_name):
- f=open("table.txt","r")
- lines=f.readlines()
- x=[]
- y=[]
- for q in lines:
- temp = q.strip().split()
- x.append(float(temp[0]))
- y.append(float(temp[1]))
- x = np.array(x)
- y = np.array(y)
- return x,y
- x, y = read_from_file("table.txt")
- def diff(x, y):
- '''x : array of data points
- y : array of f(x) '''
- n = len(x)
- a = []
- for i in range(n):
- a.append(y[i])
- for j in range(1, n):
- for i in range(n-1, j-1, -1):
- a[i] = float(a[i]-a[i-1])/float(x[i]-x[i-j])
- return np.array(a)
- the_values = diff(x, y)
- def write_to_file(the_array, file_name):
- my_list = the_array.tolist()
- with open(file_name, 'w') as f:
- for item in my_list:
- f.write("%s\n" % item)
- write_to_file(diff(x,y), "divided_differences.txt")
- def Horner(abscissaes, newton_diff, x):
- if len(abscissaes) == 1:
- # Exit condition.
- return newton_diff[0]
- term, newton_diff = newton_diff[0], newton_diff[1:]
- abscissae, abscissaes = abscissaes[0], abscissaes[1:]
- next = Horner(abscissaes, newton_diff, x)
- return term + (x - abscissae) * next
- def Horner_not_recusrive(abscissaes, newton_diff, x):
- # Reverse iteration, starting from the end makes it easier.
- value = newton_diff[-1]
- for i, abscissae in reversed(list(enumerate(abscissaes[:-1]))):
- term = newton_diff[i]
- value = term + (x - abscissae) * value
- return value
- abscissaes, ordinates = read_inputs("table.txt")
- diff = newton_diff(abscissaes, ordinates)
- x = 20
- print(Horner_not_recusrive(np.array(abscissaes), diff, x))
- print(Horner(np.array(abscissaes), diff, x))
- print(-0.5 + (x - 1) * (0.9 + (x - 2) * (-0.216667 + (x - 4) * (0.0285714 + (x - 8) * (-0.00199562)))))
- def write_to_files( shortfilename, longfilename, npoints, listx, diff, a, b ):
- shortfile = open(shortfilename, "w")
- longfile = open(longfilename, "w")
- shortfile.writelines(["%.20f\t" % item for item in diff])
- shortfile.close()
- # we create a mesh of npoints abscissae equidistant from each other by distance step
- x = a
- step = (b-a)/float(npoints-1)
- for i in range(npoints):
- y = Horner(listx,diff,x)
- longfile.write("{:.20f} \t {:.20f} \n".format(x,y))
- x = x+step
- longfile.close()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement