Skip to content

Instantly share code, notes, and snippets.

@mikofski
Last active Mar 13, 2020
Embed
What would you like to do?
horners method in 2-D
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
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
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;
\[\begin{aligned}
f\left(x,y\right)=&amp;p_1 x^n y^m+p_2 x^{\left(n-1\right)} y^m+\ldots+p_{n+1} y^m+\ldots \\
&amp;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 \\
&amp;\ldots \\
&amp;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>
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
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment