horners method in 2-D

polyDer2D.m

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\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)}$$

%% 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

polyDer2D.py

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

polyFit2D.m

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\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)}$$

%% 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;

polyVal2D.html

\[\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>

polyVal2D.m

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\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)}$$

%% 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

polyVal2D.py

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