The format of p, the polynomial coefficients, is similar to 1-D polynomials:
f(x,y)=p1xnym+p2x(n−1)ym+…+pn+1ym+…pn+2xny(m−1)+pn+3x(n−1)y(m−1)+…+p2(n+1)y(m−1)+……pm(n+1)+1∗xn+pm(n+1)+2∗x(n−1)+…+p(n+1)(m+1)
Gists:
BTW the double nested loop here is about 25% faster than 2 smaller loops with a lot of vectorization. Not sure why, but I occasionally find that loops are actually more efficient than big matrix calculations.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| function [fx,fy] = polyDer2D(p,x,y,n,m) | |
| %POLYDER2D Evaluate derivatives of 2-D polynomial using Horner's method. | |
| % F = POLYDER2D(P,X,Y,N,M) evaluates the derivatives of 2-D polynomial P at | |
| % the points specified by X and Y, which must have the same dimensions. The | |
| % outputs FX and FY will have the same dimensions as X and Y. N and M specify | |
| % the order of X and Y respectively. Polynomial coefficients are in the | |
| % following order. | |
| % | |
| % F(X,Y) = P_1 * X^N * Y^M + P_2 * X^{N-1} * Y^M + ... + P_{N+1} * Y^M + ... | |
| % P_{N+2} * X^N * Y^{M-1} + P_{N+3} * X^{N-1} * Y^{M-1} + ... + P_{2*(N+1)} * Y^{M-1} + ... | |
| % ... | |
| % P_{M*(N+1)+1} * X^N + P_{M*(N+1)+2} * X^{N-1} + ... + P_{(N+1)*(M+1)} | |
| % | |
| % See also: POLYFITN by John D'Errico on MathWorks MATLAB Central FEX | |
| % http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn | |
| % | |
| % Mark Mikofski, http://poquitopicante.blogspot.com | |
| % Version 1-0, 2013-04-17 | |
| %% LaTex | |
| % | |
| % f(x,y)=p1xnym+p2x(n−1)ym+…+pn+1ym+… | |
| % | |
| % pn+2xny(m−1)+pn+3x(n−1)y(m−1)+…+p2(n+1)y(m−1)+… | |
| % | |
| % … | |
| % | |
| % pm(n+1)+1∗xn+pm(n+1)+2∗x(n−1)+…+p(n+1)(m+1) | |
| %% check input args | |
| validateattributes(p,{'numeric'},{'2d','nonempty','real','finite'}, ... | |
| 'polyDer2D','p',1) | |
| validateattributes(x,{'numeric'},{'nonempty','real','finite'}, ... | |
| 'polyDer2D','x',2) | |
| validateattributes(y,{'numeric'},{'nonempty','real','finite'}, ... | |
| 'polyDer2D','y',3) | |
| assert(all(size(x)==size(y)),'polyDer2D:sizeMismatch', ... | |
| 'X and Y must be the same size.') | |
| % use size of p to set n & m | |
| pdims = size(p); | |
| if nargin<4 && all(pdims>1) | |
| n = pdims(1)-1; | |
| m = pdims(2)-1; | |
| end | |
| validateattributes(n,{'numeric'},{'scalar','integer','positive','<',10}, ... | |
| 'polyDer2D','n',4) | |
| validateattributes(m,{'numeric'},{'scalar','integer','positive','<',10}, ... | |
| 'polyDer2D','m',5) | |
| if all(pdims>1) | |
| assert(pdims(1)==n+1,'polyDer2D:xOrderMismatch', ... | |
| 'The number of x coefficients doesn''t match the order n.') | |
| assert(pdims(2)==m+1,'polyDer2D:yOrderMismatch', ... | |
| 'The number of y coefficients doesn''t match the order m.') | |
| end | |
| p = p(:); | |
| Np = numel(p); | |
| assert(Np==(n+1)*(m+1),'OrderMismatch', ... | |
| ['The number of x & y coefficients doesn''t match the orders ', ... | |
| 'n & m.']) | |
| %% fx = df/dx | |
| fx = n*p(1); | |
| for ni = 1:n-1 | |
| fx = fx.*x+(n-ni)*p(1+ni); | |
| end | |
| for mi = 1:m | |
| mj = (n+1)*mi+1; | |
| gx = n*p(mj); | |
| for ni = 1:n-1 | |
| gx = gx.*x+(n-ni)*p(mj+ni); | |
| end | |
| fx = fx.*y+gx; | |
| end | |
| %% fy = df/dy | |
| fy = p(1); | |
| for ni = 1:n | |
| fy = fy.*x+p(1+ni); | |
| end | |
| fy = m*fy; | |
| for mi = 1:m-1 | |
| mj = (n+1)*mi+1; | |
| gy = p(mj); | |
| for ni = 1:n | |
| gy = gy.*x+p(mj+ni); | |
| end | |
| fy = fy.*y+(m-mi)*gy; | |
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| def polyDer2D(p, x, y, n, m): | |
| """ | |
| polyDer2D(p, x, y, n, m) | |
| Evaluate derivatives of a 2-D polynomial using Horner's method. | |
| Evaluates the derivatives of 2-D polynomial `p` at the points | |
| specified by `x` and `y`, which must be the same dimensions. The | |
| outputs `(fx, fy)` will have the same dimensions as `x` and `y`. | |
| The order of `x` and `y` are specified by `n` and `m`, respectively. | |
| Parameters | |
| ---------- | |
| p : array_like | |
| Polynomial coefficients in order specified by polyVal2D.html. | |
| x : array_like | |
| Values of 1st independent variable. | |
| y : array_like | |
| Values of 2nd independent variable. | |
| n : int | |
| Order of the 1st independent variable, `x`. | |
| m : int | |
| Order of the 2nd independent variable, `y`. | |
| Returns | |
| ------- | |
| fx : ndarray | |
| Derivative with respect to x. | |
| fy : ndarray | |
| Derivative with respect to y. | |
| See Also | |
| -------- | |
| numpy.polynomial.polynomial.polyval2d : Evaluate a 2-D polynomial at points (x, y). | |
| Example | |
| -------- | |
| >>> print polyVal2D([1,2,3,4,5,6],2,3,2,1) | |
| >>> (39, 11) | |
| >>> 1*2*2*3 + 2*3 + 4*2*2 + 5 | |
| 39 | |
| >>> 1*(2**2) + 2*2 + 3 | |
| 11 | |
| >>> print polyVal2D([1,2,3,4,5,6,7,8,9],2,3,2,2) | |
| >>> (153, 98) | |
| >>> 1*2*2*(3**2) + 2*(3**2) + 4*2*2*3 + 5*3 + 7*2*2 + 8 | |
| 153 | |
| >>> 1*2*(2**2)*3 + 2*2*2*3 + 3*2*3 + 4*(2**2) + 5*2 + 6 | |
| 98 | |
| """ | |
| # TODO: check input args | |
| p = np.array(p) | |
| x = np.array(x) | |
| y = np.array(y) | |
| n = np.array(n) | |
| m = np.array(m) | |
| # fx = df/dx | |
| fx = n * p[0] | |
| for ni in np.arange(n - 1): | |
| fx = fx * x + (n - ni - 1) * p[1 + ni] | |
| for mi in np.arange(m): | |
| mj = (n + 1) * (mi + 1) | |
| gx = n * p[mj] | |
| for ni in np.arange(n - 1): | |
| gx = gx * x + (n - ni - 1) * p[mj + 1 + ni] | |
| fx = fx * y + gx | |
| # fy = df/dy | |
| fy = p[0] | |
| for ni in np.arange(n): | |
| fy = fy * x + p[1 + ni] | |
| fy = m * fy | |
| for mi in np.arange(m - 1): | |
| mj = (n + 1) * (mi + 1) | |
| gy = p[mj] | |
| for ni in np.arange(n): | |
| gy = gy * x + p[mj + 1 + ni] | |
| fy = fy * y + (m - mi - 1) * gy | |
| return fx, fy |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| function p = polyFit2D(f,x,y,n,m) | |
| %POLYFIT2D Fit data with a 2-D polynomial. | |
| % P = POLYFIT2D(F,X,Y,N,M) determines the 2-D polynomial P that best fit the | |
| % data F(X,Y) in the least squares sense. Polynomial coefficients are in the | |
| % following order. | |
| % | |
| % F(X,Y) = P_1 * X^N * Y^M + P_2 * X^{N-1} * Y^M + ... + P_{N+1} * Y^M + ... | |
| % P_{N+2} * X^N * Y^{M-1} + P_{N+3} * X^{N-1} * Y^{M-1} + ... + P_{2*(N+1)} * Y^{M-1} + ... | |
| % ... | |
| % P_{M*(N+1)+1} * X^N + P_{M*(N+1)+2} * X^{N-1} + ... + P_{(N+1)*(M+1)} | |
| % | |
| % See also: POLYFITN by John D'Errico on MathWorks MATLAB Central FEX | |
| % http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn | |
| % | |
| % Mark Mikofski, http://poquitopicante.blogspot.com | |
| % Version 1-0, 2013-03-13 | |
| %% LaTex | |
| % | |
| % f(x,y)=p1xnym+p2x(n−1)ym+…+pn+1ym+… | |
| % | |
| % pn+2xny(m−1)+pn+3x(n−1)y(m−1)+…+p2(n+1)y(m−1)+… | |
| % | |
| % … | |
| % | |
| % pm(n+1)+1∗xn+pm(n+1)+2∗x(n−1)+…+p(n+1)(m+1) | |
| %% check input args | |
| validateattributes(f,{'numeric'},{'nonempty','real','finite'}, ... | |
| 'polyFit2D','f',1) | |
| validateattributes(x,{'numeric'},{'nonempty','real','finite'}, ... | |
| 'polyFit2D','x',2) | |
| validateattributes(y,{'numeric'},{'nonempty','real','finite'}, ... | |
| 'polyFit2D','y',3) | |
| assert(all(size(x)==size(y)),'polyFit2D:sizeMismatch', ... | |
| 'X and Y must be the same size.') | |
| assert(all(size(x)==size(f)),'polyFit2D:sizeMismatch', ... | |
| 'X and F must be the same size.') | |
| validateattributes(n,{'numeric'},{'scalar','integer','positive','<',10}, ... | |
| 'polyFit2D','n',4) | |
| validateattributes(m,{'numeric'},{'scalar','integer','positive','<',10}, ... | |
| 'polyFit2D','m',5) | |
| Npoints = numel(f); | |
| assert(all((n+1)*(m+1)<Npoints),'polyFit2D:degreeTooBig', ... | |
| 'Degree must be smaller than number of data poiints.') | |
| %% fit data | |
| f = f(:); | |
| x = x(:)*ones(1,n+1); | |
| y = y(:)*ones(1,(n+1)*(m+1)); | |
| z = ones(Npoints,(n+1)*(m+1)); | |
| for ni = 1:n | |
| z(:,1:ni) = z(:,1:ni).*x(:,1:ni); | |
| end | |
| for mi = 1:m | |
| mj = (n+1)*mi; | |
| z(:,1:mj) = z(:,1:mj).*y(:,1:mj); | |
| for ni = 1:n | |
| z(:,mj+1:mj+ni) = z(:,mj+1:mj+ni).*x(:,1:ni); | |
| end | |
| end | |
| p = z\f; |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| \[\begin{aligned} | |
| f\left(x,y\right)=&p_1 x^n y^m+p_2 x^{\left(n-1\right)} y^m+\ldots+p_{n+1} y^m+\ldots \\ | |
| &p_{n+2}x^ny^{\left(m-1\right)}+p_{n+3}x^{\left(n-1\right)}y^{\left(m-1\right)}+\ldots+p_{2\left(n+1\right)}y^{\left(m-1\right)}+\ldots \\ | |
| &\ldots \\ | |
| &p_{m\left(n+1\right)+1}*x^n+p_{m\left(n+1\right)+2}*x^{\left(n-1\right)}+\ldots+p_{\left(n+1\right)\left(m+1\right)} | |
| \end{aligned} \] | |
| <script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"></script> |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| function f = polyVal2D(p,x,y,n,m) | |
| %POLYVAL2D Evaluate a 2-D polynomial using Horner's method. | |
| % F = POLYVAL2D(P,X,Y) evaluates the 2-D polynomial P at the points | |
| % specified by X and Y, which must have the same dimensions. The output F | |
| % will have the same dimensions as X and Y. N and M specify the order of X | |
| % and Y respectively. Polynomial coefficients are in the following order. | |
| % | |
| % F(X,Y) = P_1 * X^N * Y^M + P_2 * X^{N-1} * Y^M + ... + P_{N+1} * Y^M + ... | |
| % P_{N+2} * X^N * Y^{M-1} + P_{N+3} * X^{N-1} * Y^{M-1} + ... + P_{2*(N+1)} * Y^{M-1} + ... | |
| % ... | |
| % P_{M*(N+1)+1} * X^N + P_{M*(N+1)+2} * X^{N-1} + ... + P_{(N+1)*(M+1)} | |
| % | |
| % | |
| % See also: POLYFITN by John D'Errico on MathWorks MATLAB Central FEX | |
| % http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn | |
| % | |
| % Mark Mikofski, http://poquitopicante.blogspot.com | |
| % Version 1-0, 2013-03-13 | |
| %% LaTex | |
| % | |
| % f(x,y)=p1xnym+p2x(n−1)ym+…+pn+1ym+… | |
| % | |
| % pn+2xny(m−1)+pn+3x(n−1)y(m−1)+…+p2(n+1)y(m−1)+… | |
| % | |
| % … | |
| % | |
| % pm(n+1)+1∗xn+pm(n+1)+2∗x(n−1)+…+p(n+1)(m+1) | |
| %% check input args | |
| validateattributes(p,{'numeric'},{'2d','nonempty','real','finite'}, ... | |
| 'polyVal2D','p',1) | |
| validateattributes(x,{'numeric'},{'nonempty','real','finite'}, ... | |
| 'polyVal2D','x',2) | |
| validateattributes(y,{'numeric'},{'nonempty','real','finite'}, ... | |
| 'polyVal2D','y',3) | |
| assert(all(size(x)==size(y)),'polyVal2D:sizeMismatch', ... | |
| 'X and Y must be the same size.') | |
| % use size of p to set n & m | |
| pdims = size(p); | |
| if nargin<4 && all(pdims>1) | |
| n = pdims(1)-1; | |
| m = pdims(2)-1; | |
| end | |
| validateattributes(n,{'numeric'},{'scalar','integer','positive','<',10}, ... | |
| 'polyVal2D','n',4) | |
| validateattributes(m,{'numeric'},{'scalar','integer','positive','<',10}, ... | |
| 'polyVal2D','m',5) | |
| if all(pdims>1) | |
| assert(pdims(1)==n+1,'polyVal2D:xOrderMismatch', ... | |
| 'The number of x coefficients doesn''t match the order n.') | |
| assert(pdims(2)==m+1,'polyVal2D:yOrderMismatch', ... | |
| 'The number of y coefficients doesn''t match the order m.') | |
| end | |
| %% evaluate polynomial P | |
| f = p(1); | |
| for ni = 1:n | |
| f = f.*x+p(1+ni); | |
| end | |
| for mi = 1:m | |
| mj = (n+1)*mi+1; | |
| g = p(mj); | |
| for ni = 1:n | |
| g = g.*x+p(mj+ni); | |
| end | |
| f = f.*y+g; | |
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import numpy as np | |
| def polyVal2D(p, x, y, n, m): | |
| """ | |
| polyVal2D(p, x, y, n, m) | |
| Evaluate a 2-D polynomial using Horner's method. | |
| Evaluates the 2-D polynomial `p` at the points specified by `x` | |
| and `y`, which must be the same dimensions. The output `f` will | |
| have the same dimensions as `x` and `y`. The order of `x` and `y` | |
| are specified by `n` and `m`, respectively. | |
| Parameters | |
| ---------- | |
| p : array_like | |
| Polynomial coefficients in order specified by polyVal2D.html. | |
| x : array_like | |
| Values of 1st independent variable. | |
| y : array_like | |
| Values of 2nd independent variable. | |
| n : int | |
| Order of the 1st independent variable, `x`. | |
| m : int | |
| Order of the 2nd independent variable, `y`. | |
| Returns | |
| ------- | |
| f : ndarray | |
| Values of the evaluated 2-D polynomial. | |
| See Also | |
| -------- | |
| numpy.polynomial.polynomial.polyval2d : Evaluate a 2-D polynomial at points (x, y). | |
| Example | |
| -------- | |
| >>> (5**2 * 6**2 + 2 * 5 * 6**2 + 3 * 6**2 + 4 * 5**2 * 6 + | |
| ... 5 * 5 * 6 + 6 * 6 + 7 * 5**2 + 8 * 5 + 9) | |
| 2378 | |
| >>> polyVal2D([1,2,3,4,5,6,7,8,9],5,6,2,2) | |
| 2378 | |
| >>> 5 * 6 + 2 * 6 + 3 * 5 + 4 | |
| 61 | |
| >>> polyVal2D([1,2,3,4],5,6,1,1) | |
| 61 | |
| """ | |
| # TODO: check input args | |
| p = np.array(p) | |
| x = np.array(x) | |
| y = np.array(y) | |
| n = np.array(n) | |
| m = np.array(m) | |
| # evaluate f(x,y) | |
| f = p[0] | |
| for ni in np.arange(n): | |
| f = f * x + p[1 + ni] | |
| for mi in np.arange(m): | |
| mj = (n + 1) * (mi + 1) | |
| g = p[mj] | |
| for ni in np.arange(n): | |
| g = g * x + p[mj + 1 + ni] | |
| f = f * y + g | |
| return f |
