# polyfit2.py3 general polygons of 2 variables # # p = polyfit2(nx, ny, pwr, x, y, z) p computed for least square fit # # x y x y x y x p # i=0, i<=nx z[i,j] = p[0,0] + p[1,0]*x[i] + p[2,0]*x[i]^2 + p[3,0]*x[i]^3 + # j=0, j<=ny p[0,1]*y[j] + p[0,2]*y[j]^2 + p[0,3]*y[j]^3 + # p[1,1]*x[i]*y[j] + p[2,1]*x[i]^2*y[j] + # p[1,2]*x[i]*y[j]^2 # # The problem is stated as follows : # Given measured data for values of z based on values of x and y e.g. # # z i j z_actual x y z value for testing here # - - - ------- --- --- 10 z values needed for nx=3 ny=3 # 0 0 0 1.0 0.0 0.0 p for testing 1, 2, 3, ... 9, 10 # 1 1 0 10.0 1.0 0.0 # 2 2 0 49.0 2.0 0.0 # 3 3 0 169.0 3.2 0.0 # 4 0 1 23.077 0.0 1.1 # 5 0 2 102.787 0.0 2.1 # 6 0 3 282.697 0.0 3.1 # 7 1 1 191.587 1.0 1.1 # 8 2 1 152.477 2.0 1.1 # 9 1 2 191.587 1.0 2.1 # # Find p[] such that # z_approximate = p[0,0] + p[1,0]*x[i] + p[2,0]*x[i]^2 + p[3,0]*x[i]^3 + # p[0,1]*y[j] + p[0,2]*y[j]^2 + p[0,3]*y[j]^3 + # p[1,1]*x[i]*y[j] + p[2,1]*x[i]^2*y[j] + # p[1,2]*x[i]*y[j]^2 # and such that the sum of (z_actual - z_approximate) squared is minimized. # To find the coefficients p solve a system of 10 linear equations # # | SUM(1*1) SUM(1*x) SUM(1*x^2) SUM(1*x^3) # SUM(1*y) SUM(1*y^2) SUM(1*y^3) # SUM(1*x*y) SUM(1*x^2*y) SUM(1*x*y^2) | # | p[0,0] | | SUM(1*z) | # | SUM(x*1) SUM(x*x) SUM(x*x^2) SUM(x*x^3) # SUM(x*y) SUM(x*y^2) SUM(x*y^3) # SUM(x*x*y) SUM(x*x^2*y) SUM(x*x*y^2) # | | p[1,0] | | SUM(x*z) | # | SUM(x^2*1) SUM(x^2*x) SUM(x^2*x^2) SUM(x^2*x^3) # SUM(x^2*y # SUM(x^2*x*y) # | | p[2,0] | | SUM(x^2*z) | # | SUM(x^3*1) SUM(x^3*x) SUM(x^3*x^2) SUM(x^3*x^3) # SUM(x^3*y # SUM(x*3*x*y # | | p[3,0] | | SUM(x^3*z) | # | SUM(y*1) SUM(y*x) SUM(y*x^2) SUM(y*x^3) # Sum(y*y) # SUM(y*x*y) # | | p[0,1] | | SUM(y*z) | # | SUM(y^2*1) SUM(y^2*x) SUM(y^2*y^2) SUM(y^2*y^3) # SUM(y*2*y) # SUM(y^2*x*y) # | | p[0,2] | | SUM(y^2*z) | # | SUM(y^3*1) SUM(y^3*x) SUM(y^3*y^2) SUM(y^3*y^3) # SUM(y^3*y) # SUM(y^3*x*y) # | | p[0,3] | | Sum(y^3*z) | # | SUM(x*y*1) SMU(x*y*x) SUM(y^3*y^2) SUM(y^3*y^3) # SUM(x*y*y) # SUM(x*y*x*y) # | | p[1,1] | | Sum(x*y*z) | # | SUM(x^2*y*1) SUM(x^2*y*x) SUM(x^2*y*x^2) SUM(x^2*y*x^3) # SUM(x^2*y*y) SUM(y^3*y*y^2) SUM(y^3*y^3) # SUM(x^2*y*x*y) SUM(x^2*y*x^2*y) SUM(x^2*y*x*y^2) | # | p[2,1] | | Sum(x^2*y*z) | # | SUM(x*y^2*1) SUM(x*y^2*x) SUM(x*y^2*x^2) SUM(x*y^2*x^3) # SUM(x*y^2*y) SUM(x*y^2*y^2) Sum(x*y^2*y^3) # SUM(x*y^2*x*y) SUM(x*y^2*x*y^2) SUM(x*y^2*x*y^2) | # | p[1,2] | | Sum(x*y^2*z) | # solve |A|*|p|=|Y| using simeq or solve # import numpy from numpy import array from numpy import power from numpy.linalg import solve def polyval1(x, y, p): # equation written z = p[0,0] + p[1,0]*x + p[2,0]*x*x + p[3,0]*x*x*x + p[0,1]*y + p[0,2]*y*y + p[0,3]*y*y*y + p[1,1]*x*y + p[2,1]*x*x*y + p[1,2]*x*y*y return z # end polyval1 def poly(nx, ny, pwr, x, y, p): # return z i<=nx j<=ny i+j<=pwr z = p[0,0] # x=0.0 y=0.0 constant term for i in range(nx+1): for j in range(ny+1): if (i+j > pwr) or (i+j < 1): continue # end if xp = power(x,i) yp = power(y,j) z = z + p[i,j]*xp*yp # print("i=",i," xp=",xp," j=",j," yp=",yp," z=",z) # end j # end i return z # end poly def f0(x,y): # functions to build A matrix return 1.0 # end f0 def f1(x,y): return x # end f1 def f2(x,y): return x*x # end f2 def f3(x,y): return x*x*x # end f3 def f4(x,y): return y # end f4 def f5(x,y): return y*y # end f5 def f6(x,y): return y*y*y # end f6 def f7(x,y): return x*y # end f7 def f8(x,y): return x*x*y # end f8 def f9(x,y): return x*y*y # end f9 print("polyfit2.py3 running") nx = 3 print("nx=",nx) x = array([0.0, 1.0, 2.0, 3.2]) print("x=",x) ny = 3 print("ny=",ny) y = array([0.0, 1.1, 2.1, 3.1]) print("y=",y) pwr = 3 print("pwr=",pwr) p = array([[(0.0) for j in range(ny+1)] for i in range(nx+1)]) p[0,0] = 1.0 # easy to check p[1,0] = 2.0 p[2,0] = 3.0 p[3,0] = 4.0 p[0,1] = 5.0 p[0,2] = 6.0 p[0,3] = 7.0 p[1,1] = 8.0 p[2,1] = 9.0 p[1,2] = 10.0 print("p=",p) print(" ") z = array([(0.0) for k in range(16)]) z1 = array([(0.0) for k in range(16)]) z[0] = poly(nx, ny, pwr, x[0], y[0], p) print("z[0] at x[0],y[0]=",z[0]," ",x[0]," ",y[0]) z1[0] = polyval1(x[0], y[0], p) print("z1[0] at x[0],y[0]=",z[0]," ",x[0]," ",y[0]) z[1] = poly(nx, ny, pwr, x[1], y[0], p) print("z[1] at x[1],y[0]=",z[1]," ",x[1]," ",y[0]) z1[1] = polyval1(x[1], y[0], p) print("z1[1] at x[1],y[0]=",z1[1]," ",x[1]," ",y[0]) z[2] = poly(nx, ny, pwr, x[2], y[0], p) print("z[2] at x[2],y[0]=",z[2]," ",x[2]," ",y[0]) z1[2] = polyval1(x[2], y[0], p) print("z1[2] at x[2],y[0]=",z1[2]," ",x[2]," ",y[0]) z[3] = poly(nx, ny, pwr, x[3], y[0], p) print("z[3] at x[3],y[0]=",z[3]," ",x[3]," ",y[0]) z1[3] = polyval1(x[3], y[0], p) print("z1[3] at x[3],y[0]=",z1[3]," ",x[3]," ",y[0]) z[4] = poly(nx, ny, pwr, x[0], y[1], p) print("z[4] at x[0],y[1]=",z[4]," ",x[0]," ",y[1]) z1[4] = polyval1(x[0], y[1], p) print("z1[4] at x[0],y[1]=",z1[4]," ",x[0]," ",y[1]) z[5] = poly(nx, ny, pwr, x[0], y[2], p) print("z[5] at x[0],y[2]=",z[5]," ",x[0]," ",y[2]) z1[5] = polyval1(x[0], y[2], p) print("z1[5] at x[0],y[2]=",z1[5]," ",x[0]," ",y[2]) z[6] = poly(nx, ny, pwr, x[0], y[3], p) print("z[6] at x[0],y[3]=",z[6]," ",x[0]," ",y[3]) z1[6] = polyval1(x[0], y[3], p) print("z1[6] at x[0],y[3]=",z1[6]," ",x[0]," ",y[3]) z[7] = poly(nx, ny, pwr, x[1], y[2], p) print("z[7] at x[1],y[1]=",z[7]," ",x[1]," ",y[1]) z1[7] = polyval1(x[1], y[1], p) print("z1[7] at x[1],y[1]=",z1[7]," ",x[1]," ",y[1]) z[8] = poly(nx, ny, pwr, x[2], y[1], p) print("z[8] at x[2],y[1]=",z[8]," ",x[2]," ",y[1]) z1[8] = polyval1(x[2], y[1], p) print("z1[8] at x[2],y[1]=",z1[8]," ",x[2]," ",y[1]) z[9] = poly(nx, ny, pwr, x[1], y[2], p) print("z[9] at x[1],y[2]=",z[9]," ",x[1]," ",y[2]) z1[9] = polyval1(x[1], y[2], p) print("z1[9] at x[1],y[2]=",z1[9]," ",x[1]," ",y[2]) print(" ") w = array([[(0.0) for j in range(ny+1)] for i in range(nx+1)]) ij = array([[0,1,2,3,0,0,0,1,2,1],[0,0,0,0,1,2,3,1,1,2]]) # i=ij[0,k] j=ij[1,k] 10 x[i],y[j] to SUM for A and w for k in range(10): i = ij[0,k] j = ij[1,k] w[i,j] = poly(nx, ny, pwr, x[i], y[j], p) print("w[",i,",",j,"]=",w[i,j]) # end k print(" ") A = array([[(0.0) for i in range(10)] for j in range(10)]) Y = array([(0.0) for k in range(10)]) P = array([(0.0) for k in range(10)]) # print("buid A matrix and Y vector to compute p, A * P = Y") for k in range(10): i = ij[0,k] j = ij[1,k] A[0,0] = A[0,0] + f0(x[i],y[j])*f0(x[i],y[j]) A[0,1] = A[0,1] + f0(x[i],y[j])*f1(x[i],y[j]) A[0,2] = A[0,2] + f0(x[i],y[j])*f2(x[i],y[j]) A[0,3] = A[0,3] + f0(x[i],y[j])*f3(x[i],y[j]) A[0,4] = A[0,4] + f0(x[i],y[j])*f4(x[i],y[j]) A[0,5] = A[0,5] + f0(x[i],y[j])*f5(x[i],y[j]) A[0,6] = A[0,6] + f0(x[i],y[j])*f6(x[i],y[j]) A[0,7] = A[0,7] + f0(x[i],y[j])*f7(x[i],y[j]) A[0,8] = A[0,8] + f0(x[i],y[j])*f8(x[i],y[j]) A[0,9] = A[0,9] + f0(x[i],y[j])*f9(x[i],y[j]) A[1,0] = A[1,0] + f1(x[i],y[j])*f0(x[i],y[j]) A[1,1] = A[1,1] + f1(x[i],y[j])*f1(x[i],y[j]) A[1,2] = A[1,2] + f1(x[i],y[j])*f2(x[i],y[j]) A[1,3] = A[1,3] + f1(x[i],y[j])*f3(x[i],y[j]) A[1,4] = A[1,4] + f1(x[i],y[j])*f4(x[i],y[j]) A[1,5] = A[1,5] + f1(x[i],y[j])*f5(x[i],y[j]) A[1,6] = A[1,6] + f1(x[i],y[j])*f6(x[i],y[j]) A[1,7] = A[1,7] + f1(x[i],y[j])*f7(x[i],y[j]) A[1,8] = A[1,8] + f1(x[i],y[j])*f8(x[i],y[j]) A[1,9] = A[1,9] + f1(x[i],y[j])*f9(x[i],y[j]) A[2,0] = A[2,0] + f2(x[i],y[j])*f0(x[i],y[j]) A[2,1] = A[2,1] + f2(x[i],y[j])*f1(x[i],y[j]) A[2,2] = A[2,2] + f2(x[i],y[j])*f2(x[i],y[j]) A[2,3] = A[2,3] + f2(x[i],y[j])*f3(x[i],y[j]) A[2,4] = A[2,4] + f2(x[i],y[j])*f4(x[i],y[j]) A[2,5] = A[2,5] + f2(x[i],y[j])*f5(x[i],y[j]) A[2,6] = A[2,6] + f2(x[i],y[j])*f6(x[i],y[j]) A[2,7] = A[2,7] + f2(x[i],y[j])*f7(x[i],y[j]) A[2,8] = A[2,8] + f2(x[i],y[j])*f8(x[i],y[j]) A[2,9] = A[2,9] + f2(x[i],y[j])*f9(x[i],y[j]) A[3,0] = A[3,0] + f3(x[i],y[j])*f0(x[i],y[j]) A[3,1] = A[3,1] + f3(x[i],y[j])*f1(x[i],y[j]) A[3,2] = A[3,2] + f3(x[i],y[j])*f2(x[i],y[j]) A[3,3] = A[3,3] + f3(x[i],y[j])*f3(x[i],y[j]) A[3,4] = A[3,4] + f3(x[i],y[j])*f4(x[i],y[j]) A[3,5] = A[3,5] + f3(x[i],y[j])*f5(x[i],y[j]) A[3,6] = A[3,6] + f3(x[i],y[j])*f6(x[i],y[j]) A[3,7] = A[3,7] + f3(x[i],y[j])*f7(x[i],y[j]) A[3,8] = A[3,8] + f3(x[i],y[j])*f8(x[i],y[j]) A[3,9] = A[3,9] + f3(x[i],y[j])*f9(x[i],y[j]) A[4,0] = A[4,0] + f4(x[i],y[j])*f0(x[i],y[j]) A[4,1] = A[4,1] + f4(x[i],y[j])*f1(x[i],y[j]) A[4,2] = A[4,2] + f4(x[i],y[j])*f2(x[i],y[j]) A[4,3] = A[4,3] + f4(x[i],y[j])*f3(x[i],y[j]) A[4,4] = A[4,4] + f4(x[i],y[j])*f4(x[i],y[j]) A[4,5] = A[4,5] + f4(x[i],y[j])*f5(x[i],y[j]) A[4,6] = A[4,6] + f4(x[i],y[j])*f6(x[i],y[j]) A[4,7] = A[4,7] + f4(x[i],y[j])*f7(x[i],y[j]) A[4,8] = A[4,8] + f4(x[i],y[j])*f8(x[i],y[j]) A[4,9] = A[4,9] + f4(x[i],y[j])*f9(x[i],y[j]) A[5,0] = A[5,0] + f5(x[i],y[j])*f0(x[i],y[j]) A[5,1] = A[5,1] + f5(x[i],y[j])*f1(x[i],y[j]) A[5,2] = A[5,2] + f5(x[i],y[j])*f2(x[i],y[j]) A[5,3] = A[5,3] + f5(x[i],y[j])*f3(x[i],y[j]) A[5,4] = A[5,4] + f5(x[i],y[j])*f4(x[i],y[j]) A[5,5] = A[5,5] + f5(x[i],y[j])*f5(x[i],y[j]) A[5,6] = A[5,6] + f5(x[i],y[j])*f6(x[i],y[j]) A[5,7] = A[5,7] + f5(x[i],y[j])*f7(x[i],y[j]) A[5,8] = A[5,8] + f5(x[i],y[j])*f8(x[i],y[j]) A[5,9] = A[5,9] + f5(x[i],y[j])*f9(x[i],y[j]) A[6,0] = A[6,0] + f6(x[i],y[j])*f0(x[i],y[j]) A[6,1] = A[6,1] + f6(x[i],y[j])*f1(x[i],y[j]) A[6,2] = A[6,2] + f6(x[i],y[j])*f2(x[i],y[j]) A[6,3] = A[6,3] + f6(x[i],y[j])*f3(x[i],y[j]) A[6,4] = A[6,4] + f6(x[i],y[j])*f4(x[i],y[j]) A[6,5] = A[6,5] + f6(x[i],y[j])*f5(x[i],y[j]) A[6,6] = A[6,6] + f6(x[i],y[j])*f6(x[i],y[j]) A[6,7] = A[6,7] + f6(x[i],y[j])*f7(x[i],y[j]) A[6,8] = A[6,8] + f6(x[i],y[j])*f8(x[i],y[j]) A[6,9] = A[6,9] + f6(x[i],y[j])*f9(x[i],y[j]) A[7,0] = A[7,0] + f7(x[i],y[j])*f0(x[i],y[j]) A[7,1] = A[7,1] + f7(x[i],y[j])*f1(x[i],y[j]) A[7,2] = A[7,2] + f7(x[i],y[j])*f2(x[i],y[j]) A[7,3] = A[7,3] + f7(x[i],y[j])*f3(x[i],y[j]) A[7,4] = A[7,4] + f7(x[i],y[j])*f4(x[i],y[j]) A[7,5] = A[7,5] + f7(x[i],y[j])*f5(x[i],y[j]) A[7,6] = A[7,6] + f7(x[i],y[j])*f6(x[i],y[j]) A[7,7] = A[7,7] + f7(x[i],y[j])*f7(x[i],y[j]) A[7,8] = A[7,8] + f7(x[i],y[j])*f8(x[i],y[j]) A[7,9] = A[7,9] + f7(x[i],y[j])*f9(x[i],y[j]) A[8,0] = A[8,0] + f8(x[i],y[j])*f0(x[i],y[j]) A[8,1] = A[8,1] + f8(x[i],y[j])*f1(x[i],y[j]) A[8,2] = A[8,2] + f8(x[i],y[j])*f2(x[i],y[j]) A[8,3] = A[8,3] + f8(x[i],y[j])*f3(x[i],y[j]) A[8,4] = A[8,4] + f8(x[i],y[j])*f4(x[i],y[j]) A[8,5] = A[8,5] + f8(x[i],y[j])*f5(x[i],y[j]) A[8,6] = A[8,6] + f8(x[i],y[j])*f6(x[i],y[j]) A[8,7] = A[8,7] + f8(x[i],y[j])*f7(x[i],y[j]) A[8,8] = A[8,8] + f8(x[i],y[j])*f8(x[i],y[j]) A[8,9] = A[8,9] + f8(x[i],y[j])*f9(x[i],y[j]) A[9,0] = A[9,0] + f9(x[i],y[j])*f0(x[i],y[j]) A[9,1] = A[9,1] + f9(x[i],y[j])*f1(x[i],y[j]) A[9,2] = A[9,2] + f9(x[i],y[j])*f2(x[i],y[j]) A[9,3] = A[9,3] + f9(x[i],y[j])*f3(x[i],y[j]) A[9,4] = A[9,4] + f9(x[i],y[j])*f4(x[i],y[j]) A[9,5] = A[9,5] + f9(x[i],y[j])*f5(x[i],y[j]) A[9,6] = A[9,6] + f9(x[i],y[j])*f6(x[i],y[j]) A[9,7] = A[9,7] + f9(x[i],y[j])*f7(x[i],y[j]) A[9,8] = A[9,8] + f9(x[i],y[j])*f8(x[i],y[j]) A[9,9] = A[9,9] + f9(x[i],y[j])*f9(x[i],y[j]) print("k=",k," i=",i," j=",j," x[i]=",x[i]," y[j]=",y[j]) # end k print(" ") print("A=") print(A) print(" ") for k in range(10): i = ij[0,k] j = ij[1,k] Y[0] = Y[0] + polyval1(x[i],y[j],p) Y[1] = Y[1] + polyval1(x[i],y[j],p)*x[i] Y[2] = Y[2] + polyval1(x[i],y[j],p)*x[i]*x[i] Y[3] = Y[3] + polyval1(x[i],y[j],p)*x[i]*x[i]*x[i] Y[4] = Y[4] + polyval1(x[i],y[j],p)*y[j] Y[5] = Y[5] + polyval1(x[i],y[j],p)*y[j]*y[j] Y[6] = Y[6] + polyval1(x[i],y[j],p)*y[j]*y[j]*y[j] Y[7] = Y[7] + polyval1(x[i],y[j],p)*x[i]*y[j] Y[8] = Y[8] + polyval1(x[i],y[j],p)*x[i]*x[i]*y[j] Y[9] = Y[9] + polyval1(x[i],y[j],p)*x[i]*y[j]*y[j] # end k print("Y=") print(Y) print(" ") P = solve(A,Y) print("P=") print(P) print(" ") print("polyfit2.py3 finished")