Showing posts with label polynomial segment. Show all posts
Showing posts with label polynomial segment. Show all posts

Wednesday, January 16, 2013

Length of a Polynomial Segment - Part 1; and the HP 50g Program MPFIT


The next two blog entries are going to deal with finding the length of a polynomial segment.

Polynomial Segments

What is a polynomial segment? I could not find a formal definition, so I'll describe it is this: similar to a line segment where a pair of points are connected with a line, a polynomial segment is a set of points that are connected by polynomial curve. The degree of the polynomial is determined by the number of points, less one. For instance, if we have 3 points, the polynomial that connects them the order 2 (quadratic). A cubic polynomial (order 3) is used to connect a set of 4 points.

The Goal

We have a set of points {(x0, y0), (x1, y1), (x2, y2),... (xn, yn)}. Connect the points using a polynomial. What is the length of that segment?

General Method

1. Find the polynomial to connect the set of points. Order of the elements is important. This is determined by using the matrix calculation (X^T X)^-1 (X^T y) to determine the coefficients. The resulting matrix contains the coefficients of the polynomial; with each term being a coefficient of x of increasing power, from 0 to n-1.

2. Use the arc length formula ∫ ( √ (1 + (f'(x))^2 ) dx, a, b) where a is the minimum of the x values and b is the maximum of the x values.

We recommend using a calculator or math software that can handle matrices and integrals for this procedure.

But what is X and y?

y represents a vector of y-values of the set of points.

X is a Vandermonde Matrix. A Vandermonde matrix is a n × n matrix is formed with each row containing a geometric progression of each term. The first column, each term is raised to the 0th power, the second column each term to the 1st, the third column each term is raised to the second power, and so on until the nth row where each term is raised to the n-1 power.

For instance, a Vandermonde matrix consisting of the set [a, b, c, d] would look like this:

Example 1

Find a polynomial segment that fits the points (0,0), (4,8), and (8,6). Find the length of that segment. (See the diagram above).

The x values are [0, 4, 8] with xmin = 0 and xmax = 8. The y values are [0, 8, 6]. Let X be the Vandermonde matrix with the set [0, 4, 8] and the vector y be represented by the corresponding y-values.

Now we need to calculate (X^T X)^-1 (X^T y).

Then:

The resulting matrix contains the coefficients of the required polynomial. Since the vector has 3 entries, the order of the polynomial is 3-1=2.

Going top-down:
0 is the constant (coefficient of x^0),
3.25 is the coefficient of x,
and -0.3125 is the coefficient of x^2.

Our polynomial is f(x) = 3.25x - 0.03125x^2.

To find the length, use the arc length formula using end points xmin = 0 and xmax = 8. We find that the length of the polynomial segment connecting (0,0), (4,8), and (8,6) is about 14.23920.

Example 2

What is the length of the polynomial segment connected with the points (-2,5), (0,-6), (3,2), and (4,1). We are working with 4 points, the required polynomial has the order 4-1=3. (Cubic polynomial)

The Vandermonde matrix X with the vector y are as follows:

Calculating (X^T X)^-1 (X^T y) yields (approximately):

[ [ -6 ]
[ 2.98333 ]
[ 2.725 ]
[ -0.75833 ] ]

To form the cubic polynomial:
-6 is the constant,
2.98333 is the coefficient of x,
2.725 is the coefficient of x^2, and
-0.75833 is the coefficient of x^3.

With xmin = -2 and xmax = 4, we find the length of this segment is approximately 32.60952. (See below).

HP50g Program MPFIT

The HP 50g program MPFIT fits a polynomial to the set of points {(x0, y0), (x1, y1), (x2, y2), ... , (x_(n-1), y_(n-1))}.

Input:
2: vector of x coordinates
1: vector of y coordinates

Output:
2: Coefficients in matrix form [ [a0] [a1] [a2] ... [a_(n-1)] ]
1: Polynomial of X in the form a0 + a1 * X + a2 * X^2 + ... + a_(n-1) * X^(n-1)

Program MPFIT (size 183.5 bytes)

<< DUP SIZE OBJ→ DROP → MX MY N
<< 'X' PURGE
MY OBJ→ 1 + →ARRY
MX VANDERMONDE
LSQ → MR
<< MR DUP
{1, 1} GET
2 N FOR K
MR K 1 2 →LIST GET
'X' K 1 - ^ * +
NEXT >> >> >>

Example:
Input:
2: [1, 3, 6]
1: [6, 3, 5]

Output:
2: [[8.8] [-3.23333333333] [0.43333333333]]
1: 8.8+-3.23333333333*X+0.43333333333*X^2

Which means the polynomial fit for the given points is (approximately):

8.8 - 3.23333x + 0.43333x^2

(Screen shots of this example are shown below)

To get the length of the segment with the HP 50g, first note the xmin and xmax values, which are 1 and 6, respectively. This keystroke (in RPN Mode) should do the trick:

[ ' ] [ X ] [ENTER]
[RS] [COS] ( δ )
[LS] [ √ ] (x²)
1 [ + ] [ √ ]
1 [ENTER]
6 [ENTER]
3 [LS] [EVAL] (PRG) [F1] (STACK) [NXT] [F1] (ROLL)
[ ' ] [ X ] [ENTER]
[RS] [TAN] ( ∫ )

Length = 7.75662832829



My next blog entry will contain programs for the TI-84+ and Casio Prizm (fx series) calculators.

Have a great day everyone,

Eddie


This blog is property of Edward Shore. 2013

Length of a Polynomial Segment - Part 2 (TI-84+/Casio); From Coffee Klatch in San Dimas

Today I am coming to you from Coffee Klatch in San Dimas, CA. This is one of my most favorite places to hang out, and have great coffee and food.

Today's blog is the second part of the polynomial segment blog. Here is a link to Part 1.


To Recall:

A polynomial segment is a set of points that is connected by a polynomial.

To find a perfectly fit polynomial to fit n points:

1. Let X be the Vandermonde matrix of order (n-1) × (n-1) for the following:

2. Let y be the vector of y values.
3. The coefficient matrix is calculated as (X^T X)^-1 (X^T y). The result is a vector of coefficients of x^k where k = 0 to n-1.

On most graphing (and some scientific calculators), you can accomplish this task in two ways:

1. Use matrix operations. This is good for any size polynomial.
2. Take advantage of the curve fitting functions LinReg (to fit 2 points), QuadReg (to fit 3 points), CubicReg (to fit 4 points), and QuartReg (to fit 5 points). *Regression models and names vary for each calculator.

Once the polynomial, f(x) is found, find the arc length is found by:

The limits of the integral are the minimum and maximum x values, respectively.




Here are two programs that accomplish the task of finding a polynomial fit for 3, 4, or 5 points for the TI-84+ (and its family), and the Casio Prizm (and its fx-9xxx family).


TI-84+

POLYARC
1/12/2013
Fit a polynomial to a set of three, four, or five points.
Find the arc length of the polynomial from the end points.

Example 1: (0,0), (4,8), (8,6)
Results: y = -0.3125x^2 + 3.25
Length = 14.23920

Example 2: (0,3), (2,-1), (4,2), (6,0)
Results: y = -0.25x^3 + 2.375x^2 - 5.75x + 3
Length = 12.46834

Example 3: (0,0), (1,1), (2,0), (3,-1), (4,0.5)
Results: y ≈ -0.02083333x^4 + 0.20833333x^3 - 1.77083333x^2 + 2.54166666x
Length = 6.36142413

Program POLYARC (327 bytes)
Input "NO OF PTS(3-5)", N
If N<3 or N>5
1/0
N → dim(L1)
N → dim(L2)
For(K,1,N)
Disp K
Input "X=", X
Input "Y=", Y
X → L1(K)
Y → L2(K)
End
If N=3
Then
QuadReg L1, L2
fnInt(√(1+(2aX+b)²),X,min(L1),max(L1)) → L
Pause {a, b, c}
End
If N=4
Then
CubicReg L1, L2
Pause {a, b, c, d}
fnInt(√(1+(3aX² + 2bX + c)²), X, min(L1), max(L1)) → L
End
If N=5
Then
QuartReg L1, L2
Pause {a, b, c, d, e}
fnInt(√(1 + (4aX³ + 3bX² + 2cX + d)² , X, min(L1), max(L2)) → L
End
Disp "LENGTH="
Pause L

Where to find:
a, b, c, d, e: VARS, 5 (Statistics), right, right, [2, 3, 4, 5, and 6 respectively]
You need the lower case a, b, c, d, and e!

Note: No "bending backwards" is allowed


Casio Prizm

POLYARC
1/12/2013
Fit a polynomial to a set of three or four points.
Find the arc length of the polynomial from the end points.

Example 1: (0,0), (4,8), (8,6)
Results: y = -0.3125x^2 + 3.25
Length = 14.23920

Example 2: (0,3), (2,-1), (4,2), (6,0)
Results: y = -0.25x^3 + 2.375x^2 - 5.75x + 3
Length = 12.46834

Example 3: (0,0), (1,1), (2,0), (3,-1), (4,0.5)
Results: y ≈ -0.02083333x^4 + 0.20833333x^3 - 1.77083333x^2 + 2.54166666x
Length = 6.36142413

Program POLYARC (352 bytes):
"NO. OF POINTS(3,4,5)"? → N
N<3 Or N>5 ⇒ 1 ÷ 0
N → Dim List 1
N → Dim List 2
For 1 → K To N
K ◢
"X"? → List 1[K]
"Y"? → List 2[K]
Next
If N = 3
Then QuadReg List 1, List 2
∫( √(1 + (2aX + b)² ), Min(List 1), Max(List 1)) → L
{a, b, c} ◢
EndIf
If N=4
Then CubicReg List 1, List 2
∫( √(1 + (3aX ² + 2bX + c)² ), Min(List 1), Max(List 1)) → L
{a, b, c, d} ◢
EndIf
If N=5
Then QuartReg List 1, List 2
∫ ( √(1 + (4aX^3 + 3bX ² + 2cX + d)² ), Min(List 1), Max(List 2)) → L
{a, b, c, d, e} ◢
EndIf
"LENGTH="
L

How to Access:
a, b, c, d, and e: VARS, STAT (F3), GRAPH (F3), F1 (F2, F3, F4, and F5 respectively)
QuadReg: back to default menu, MENU, STAT, CALC, X^2
CubicReg: back to default menu, MENU, STAT, CALC, X^3
QuartReg: back to default menu, MENU, STAT, CALC, X^4

Note: No "bending backwards" allowed



Enjoy the day and I'll talk to you soon!

Eddie


This blog is property of Edward Shore. 2013

Python – Earth’s Radius and Gravity in US Units

Python – Earth’s Radius and Gravity in US Units Introduction The following script, gravus2.py, estimates the Earth’s gravity i...