Curve fitting for y = ce^(Ax)
- Given: (x[1],y[1]), (x[2],y[2]), ... (x[n],y[n])
- Fitting these points to an exponential curve
- E(A,C) = sum[k=1;n](ce^(Ax[k]) - y[k])^2
- Want to minimize error, so we take the derivative and solve
c sum[k=1;n](x[k]e^(2Ax[k])) - sum[k=1;n](y[k]x[k]e^(Ax[k]))
= 0
c sum[k=1;n](e^(2Ax[k])) - sum[k=1;n](y[k]e^(Ax[k])) = 0
- These equations can be solved iteratively using the Newton method
Transformations for data linearization
- Attempting to make y = ce^(Ax) look like y = Ax + B
- ln y = Ax + ln c by log of both sides, which we write as z = Ax + B
Example
-------
(x[k], y[k]) (0, 1.5) (1,
2.5) (2, 3.5) (3,
5.0) (4, 7.5) | 10
z[k] 0.4055 0.9163 1.2528 1.6094 2.0149 | 6.1989
x[k]^2 0 1 4 9 16 | 30
x[k]z[k] 0 0.9113 2.5036 4.8282 8.0596 | 16.3097
Substituting these values into A and B, we get
A(30) + B(10) = 16.3097
A(10) + B(5) = 6.1989
Solving this system results in A = 0.3912, B = .4574
y = 1.5799e^(0.3912x), which was 1.6108995e^(0.3835705x) via Newton method
------------------------------------------------------------------------------
Polynomial fitting
- Given: (x[1],y[1]), (x[2],y[2]), ... (x[n],y[n])
- Want to find c[1], c[2], ... c[m+1] that best fit the given points
- Polynomial is f(x) = c[1] + c[2]x + c[3]x^2 + ... c[m+1]x^m
Polynomial approx. vs. fitting
- Lagrange/Newton take m points, give m-1 degree polynomial
- Fitting is given degree in advance, which is independent of the points
- Approximation results in a curve which passes through each point
- Fitting results in a curve which will pass close/reduce error
Other transformations used to fit curves
- y = A/x + B becomes y = Az + B with z = 1/x
- y = Cx^A becomes w = Az + B, z = ln(x), w = ln(y)
- y = Ce^(Ax) becomes z = Ax + B, z = ln(y), B = ln(C)
Polynomial fitting
- We want to fit (x[1],y[1]), (x[2],y[2]), etc. into y = Ax^2 + Bx + C
- Again, minimizing least square error
A*sum[k=1;n]x[k]^4+B*sum[k=1;n]x[k]^3+C*sum[k=1;n]x[k]^2=sum[k=1;n]y[k]x[k]^2
A*sum[k=1;n]x[k]^3+B*sum[k=1;n]x[k]^2+C*sum[k=1;n]x[k] =sum[k=1;n]y[k]x[k]
A*sum[k=1;n]x[k]^2+B*sum[k=1;n]x[k] +C =sum[k=1;n]y[k]
Questions about homework
- HW1: mkarahan@cs.purdue.edu
- HW2: gpark@cs.purdue.edu
Piecewise Cubic Splines
- Want to fit n+1 data points into n cubic polynomials S[k](x) that have coefficients
of S[k,0], S[k,1], S[k,2], S[k,3] from k=1 to n
- These methods, when carried to 3D, are used for CGI
- Rules
* S[k](x) = S[k,0] + S[k,1]x + S[k,2]x^2 + S[k,3]x^3
for [x[k],x[k+1]]
* S[k](x[k]) = y[k]
* S[k](x[k+1]) = S[k+1](x[k+1])
* S'[k](x[k+1]) = S'[k+1](x[k+1])
* S"[k](x[k+1]) = S"[k+1](x[k+1])
- Polynomial approximation gets "wiggly" as n gets large
- Curve fitting may not be as "wiggly" but doesn't pass through all
points
- Cubic splines don't get "wiggly" and pass through all points
- For each polynomial S[k](x), there are 4 coefficients (degrees of freedom)
- Total, there are 4n coefficients
- There are n+1 constraints from rule #2, n-1 each from #3, 4, 5
- Total, there are 4n-2 constraints
- Therefore, there are 2 additions equations needed for end-point constraints
(using S'(x) and S"(x) at x[0] and x[n]), to be discussed later
- Since the spline is cubic, we can obtain equations using Lagrange approx.
- m[k] = S"(x[k]), m[k+1] = S"(x[k+1]) and h[k] = x[k+1] - x[k]
- S"[k](x) = m[k]((x[k+1] - x[k])/h[k]) + m[k+1]((x - x[k])/h[k])
7-15-04
Splines (cont.)
- Better than Bezier curves due to continuity in the second derivative
- The second derivative is a line (since the spline is cubic)
S"[k](x)=S"(x[k])((x-x[k+1])/(x[k]-x[k+1]))+S"(x[k+1])((x-x[k])/(x[k+1]-x[k]))
S[k](x) =∫[x[k];x[k+1]]∫[x[k];x[k+1]] S"[k](x)
S[k](x) =∫[x[k];x[k+1]]∫[x[k];x[k+1]]((m[k](x[k+1]-x)/(h[k]))+(m[k+1](x-x[k])/(h[k])))dx
S[k](x) =(m[k]/6h[k])(x[k+1]-x)^3+(m[k+1]/6h[k])(x-x[k])^3+c[1]x+c[2]
S[k](x) =(m[k]/6h[k])(x[k+1]-x)^3+(m[k+1]/6h[k])(x-x[k])^3+p[k](x[k+1]-x)+q[k](x-x[k])
S[k](x[k]) = y[k] = (m[k]/6)(h[k])^2 + P[k](h[k])
p[k] = (y[k] - (m[k]/6)(h[k])^2)(1/h[k])
p[k] = y[k]/h[k] - (m[k]/6)h[k]
Do the same for x = x[k+1]
q[k] = y[k+1]/h[k] - (m[k+1]/6)h[k]
Substituting back, we get:
S[k](x) = (m[k]/6h[k])(x[k+1]-x)^3+(m[k+1]/6h[k])(x-x[k])^3+(y[k]/h[k] - (m[k]/6)h[k])(x[k+1]-x)+(y[k+1]/h[k]
- (m[k+1]/6)h[k])(x-x[k])
This equation is closer to the spline we're looking for, though we need m[i]
S'[k](x) = -(m[k]/2h[k])(x[k+1]-x)^2+(m[k]/2h[k])(x-x[k])^2-((y[k]/h[k])-(m[k]h[k])/6)+((y[k+1]/h[k])-(m[k+1]h[k])/6)
d[k] = (y[k+1]-y[k])/h[k]
S'[k](x[k]) = (-1/3)m[k]h[k] - (m[k+1]h[k])/6) + d[k]
S'[k-1](x) = (-1/2)(m[k-1]/h[k-1])(x[k]-x)^2+(1/2)(m[k]/h[k-1])(x-x[k-1])^2-((y[k-1]/h[k-1])-(m[k-1]h[k-1])/6)+((y[k]/h[k-1])-(m[k]h[k-1])/6)
S'[k-1](x[k]) = (1/3)(m[k]h[k-1])+(1/6)(m[k-1]h[k-1])+d[k-1]
6(d[k]-d[k-1]) = m[k-1]h[k-1] + 2m[k](h[k-1]+h[k]) + m[k+1]h[k]
U[k] = m[k-1]h[k-1] + 2m[k](h[k-1]+h[k]) + m[k+1]h[k]
This gives us n-1 equations to solve for m[0], m[1], ... m[n], needing the
endpoint equations in order to finally solve for all m[i]
7-16-04
Splines (cont. again)
- The equations give us a strictly diagonal dominant matrix
- s[k,0] = y[k]
- s[k,1] = d[k] - (h[k](2m[k]+m[n+1])/6
- s[k,2] = m[k]/2
- s[k,3] = (m[k+1]-m[k])/(6m[k])
- Endpoint constraints
* Clamped spline = first derivatives at end points are
given
* Natural spline = second derivatives at end points are
zero
* Extrapolated spline = S[0] and S[1] are the same spline,
S[n-2] and S[n-1]
are the same spline (used in order to connect multiple splines)
* Parabolic terminated spline = S[0] and S[n-1] are quadratic
polynomials
* End-point curvature adjusted spline = First/second
derivatives given
Example
-------
Find the natural spline that approximates y = sin(x) in the interval [0,π/2]
with 4 points. These points will be 0, π/6, π/3, and π/2.
k x[k] y[k]
------------------------------
0 0 0
1 π/6 = 0.3236 0.5
2 π/3 = 1.0472 0.8660
3 π/2 = 1.5708 1
h[0] = h[1] = h[2] = π/6
= 0.5236
d[0] = (y[1]-y[0])/h[0] = 0.9549
d[1] = (y[2]-y[1])/h[1] = 0.6990
d[2] = (y[2]-y[1])/h[2] = 0.2559
u[1] = 6(d[1]-d[0]) = -1.5354
u[2] = 6(d[2]-d[1]) = -2.6586
m[0] = 0, m[3] = 0 (from end-point constraints)
2(h[0]+h[1])m[1]+h[1]m[2]=u[1]
h[n-2]m[n-2]+2(h[n-2]+h[n-1])m[n-1]=u[n-1]
2.0944 m[1] + 0.5236 m[2] = -1.5354
0.5326 m[1] + 2.0944 m[2] = -2.6386
m[1] = -0.4435
m[2] = -1.1585
S[0,0] = 0
S[0,1] = 0.9936
S[0,2] = 0
S[0,3] = -0.1412
S[0](x) = 0.9936 x - 0.1412 x^3, 0 ≤ x ≤ π/6
S[1](x) = 0.5 + 0.8775(x-0.5236) + 0.2218(x-0.5236)^2 - 0.3276(x-0.5236)^3
S[2](x) = 0.8660 + 0.4581(x-1.0472) + 0.5743(x-1.0472)^2 + 0.3688(x-1.0472)^3