Polynomial fit passing through specified points
Jiro's pick this week is polyfix by Are Mjaavatten.
Have you ever wanted to fit a polynomial to your data and have the line go through some specified points? What about specifying the slope at a certain point? Let's take a look at some options, including Are's entry.
Contents
Polynomial fitting
The function polyfit lets you fit a polynomial to your data. We'll first create some data.
t = linspace(0,2,31)';
y = sin(pi*t)+ 0.1*randn(size(t));
plot(t,y,'o')
We'll fit a 3rd order polynomial to the data.
$$ y \approx x_1t^3 + x_2t^2 + x_3t + x_4$$
x_polyfit = polyfit(t,y,3) y1 = polyval(x_polyfit,t); plot(t,y,'o',t,y1) legend('data','polyfit')
x_polyfit = 2.8843 -8.6941 6.0419 -0.2050
Constrain to go through certain points
What if you want this polynomial to go through certain points. Perhaps, you want the curve to cross (0, 0) and (2, 0). This is where Are's entry comes into play. But first, let me talk about a different method. I found this question on MATLAB Answers. There are several ways to deal with this, and one of them is to use a function like lsqlin from Optimization Toolbox. lsqlin solves the following least-squares curve fitting problem.
$$\min_{x}\frac{1}{2}||C\cdot x - d||_2^2 \quad \rm {such \; that} \quad \Bigg \{ \begin{array}{l} A \cdot x \leq b, \\ Aeq \cdot x = beq, \\ lb \leq x \leq ub \end{array} $$
For a 3rd order polynomial, $C$ and $d$ are defined as the following.
$$C = \left[ t^3 \quad t^2 \quad t \quad 1 \right],\quad d = y$$
$x$ represents the coefficients for the polynomial terms.
$$x : \left[ \begin{array}{c} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array} \right]$$
C = [t.^3 t.^2 t ones(size(t))]; d = y;
There are no inequality constraints, so we set $A$ and $b$ to empty.
A = []; % No inequality constraint b = []; % No inequality constraint
For the equality constraints, we need to set the $y$ values at $t$ = 0 and $t$ = 2 to zero. This means
$$\begin{array}{rccccc} x_1(0)^3 + x_2(0)^2 + x_3(0) + x_4 = 0 & \Rightarrow & [0\quad0\quad0\quad1] & \cdot & x = & 0 \\ x_1(2)^3 + x_2(2)^2 + x_3(2) + x_4 = 0 & \Rightarrow & \underbrace{[2^3\quad2^2\quad2\quad1]}_{Aeq} & \cdot & x = & \underbrace{0}_{beq} \end{array}$$
Putting these into $Aeq$ and $beq$, we get the following.
Aeq = [0 0 0 1 ; % t = 0 2^3 2^2 2 1]; % t = 2 beq = [0 ; % y = 0 (when t = 0) 0]; % y = 0 (when t = 2)
We then call lsqlin to solve for the polynomial coefficients.
x_lsqlin1 = lsqlin(C,d,A,b,Aeq,beq)
Warning: The trust-region-reflective algorithm can handle bound constraints only; using active-set algorithm instead. Optimization terminated. x_lsqlin1 = 2.5524 -7.6807 5.1518 0
plot(t,y,'o',t,C*x_lsqlin1) legend('data','lsqlin') grid on
Constrain to have a certain derivative
Let's add more constraint for the curve to go through (0.5,1) and (1.5,-1). Also, we want the derivatives at those points to be zero. That's additional 4 constraints. With so many constraints, we won't be able to satisfy them with a 3rd order polynomial. Let's try a 7th order.
$$ y \approx x_1t^7 + x_2t^6 + x_3t^5 + x_4t^4 + x_5t^3 + x_6t^2 + x_7t + x_8$$
C = [t.^7 t.^6 t.^5 t.^4 t.^3 t.^2 t ones(size(t))];
The derivative for the 7th order polynomial is defined as
$$\frac{dy}{dt} \approx 7x_1t^6 + 6x_2t^5 + 5x_3t^4 + 4x_4t^3 + 3x_5t^2 + 2x_6t + x_7$$
For example, the constraint "the derivative at $t$ = 0.5 is zero" can be represented as
$$ \begin{array}{c} { 7x_1(0.5)^6 + 6x_2(0.5)^5 + 5x_3(0.5)^4 + 4x_4(0.5)^3 + 3x_5(0.5)^2 + 2x_6(0.5) + x_7 = 0 } \\ \\ { \Downarrow } \\ \\ { \underbrace{\left[ 7(0.5)^6\quad6(0.5)^5\quad5(0.5)^4\quad4(0.5)^3\quad3(0.5)^2\quad2(0.5)\quad1\quad0\right]}_{Aeq}\cdot x = \underbrace{0}_{beq} } \end{array}$$
Aeq = [0 0 0 0 0 0 0 1 ; % t = 0 2^7 2^6 2^5 2^4 2^3 2^2 2 1 ; % t = 2 0.5^7 0.5^6 0.5^5 0.5^4 0.5^3 0.5^2 0.5 1 ; % t = 0.5 1.5^7 1.5^6 1.5^5 1.5^4 1.5^3 1.5^2 1.5 1 ; % t = 1.5 7*0.5^6 6*0.5^5 5*0.5^4 4*0.5^3 3*0.5^2 2*0.5 1 0 ; % t = 0.5 (for dy/dt) 7*1.5^6 6*1.5^5 5*1.5^4 4*1.5^3 3*1.5^2 2*1.5 1 0]; % t = 1.5 (for dy/dt) beq = [ 0 ; % y = 0 (when t = 0) 0 ; % y = 0 (when t = 2) 1 ; % y = 1 (when t = 0.5) -1 ; % y = -1 (when t = 1.5) 0 ; % dy/dt = 0 (when t = 0.5) 0]; % dy/dt = 0 (when t = 1.5)
Let's fit and see if we've accomplished our goal.
x_lsqlin2 = lsqlin(C,d,A,b,Aeq,beq) plot(t,y,'o',t,C*x_lsqlin2) legend('data','lsqlin (with derivative constraint)') grid on set(gca,'XTick',0:0.5:2) % Adjust tick marks to show points of interest
Warning: The trust-region-reflective algorithm can handle bound constraints only; using active-set algorithm instead. Optimization terminated. x_lsqlin2 = -0.2925 2.1030 -7.8142 17.6810 -19.6794 5.7234 2.2753 0.0000
polyfix
In the MATLAB Answers post I mentioned above, Are actually posted a response mentioning polyfix. This entry achieves the goal of performing a polynomial fit with constraints to pass through specific points with specific derivatives. Let's solve the same problem using polyfix, 7th order polynomial to fit through (0,0), (2,0), (0.5,1), (1.5,-1) and derivative of zero at $t$ = 0.5 and $t$ = 1.5.
x_polyfix = polyfix(t,y,7,[0 2 0.5 1.5],[0 0 1 -1],[0.5 1.5],[0 0]) plot(t,y,'o',t,polyval(x_polyfix,t)) legend('data','polyfix') grid on set(gca,'XTick',0:0.5:2) % Adjust tick marks to show points of interest
x_polyfix = Columns 1 through 7 -0.2925 2.1030 -7.8142 17.6810 -19.6794 5.7234 2.2753 Column 8 0
If we compare this result with that from lsqlin, we see that it's essentially identical.
norm(x_lsqlin2-x_polyfix')
ans = 2.1194e-10
Optimization Toolbox
Are's polyfit is great for performing polynomial fits with constraints around passing points. If you need other sophisticated constraints, you would want to check out Optimization Toolbox.
Comments
Give this a try and let us know what you think here or leave a comment for Are.
- 범주:
- Picks
댓글
댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.