Curve Fitting for y = cekx

Suppose we have ( x1, y1) , (x2, y2)… (xn,yn) and we want to fit these points to an exponential curve y = CeAx.

We use the least-square procedure E(A,C) = ∑(CeAx – yk)2

We want to minimize E(A,C) . So, we obtain the partial derivatives ∂E/∂A and also ∂E/∂C and set it to 0 .

∂E/∂A = ∑2(CeAx –yk)CxkeAxk = 0

∑(2C2xke2Axk - 2CxkykeAxk) = 0

∑Cxke2Axk - ∑xkykeAxk = 0 - I

∂E/∂C = ∑2(CeAxk – yk) eAxk) = 0

∑Ce2Axk - ∑ykeAxk = 0 - II


I and II form a system of 2 non-linear equations that can be solved iteratively using methods like Newton for systems of non-linear eqns.

We can find A and C in y = CeAx that minimizes the square error without having to solve the system of non linear equations. We can transform y = CeAx into a linear equation using a transformation

Transformations for data linearization

We can transform y = CeAx into a curve of the form y = Ax+ B
If we make the substitution z = ln(y)

ln y = ln(CeAx)
= ln C + Ax
z = Ax + b

we are going to try to fit the data into z = Ax + B.

xk yk zk = lnyk xk2 xkzk

0 1.5 0.4055 0 0
1 2.5 0.9163 1 0.9163
2 3.5 1.2528 4 2.5056
3 5 1.6094 9 4.8282
4 7.5 2.0149 16 8.0596
∑xk = 10 ∑xk2 = 30 ∑zk = 6.1989 ∑xkzk = 16.3097

A∑xk2 + B∑xk = ∑xkzk I
A∑xk + BN = ∑zk II

30A + 10B = 16.3097
10A + 5B = 6.1989

Solving for A, B

A = 0.3912
B = 0.4574

z = 0.3912x + 0.4574

y = CeAx
B = lnC
C = eB
= e.4574
= 1.5799

y = 1.5799e0.3912x

If we solve the non linear least square for y = CeAx without using a substitution , instead using an iterative method , we would get y = 1.670899e0.38705x which is closer to the solution than using transformation.


Polynomial fitting

We want to find c1,c2,… cn+1 that best fits the points (x1,y1),(x2,y2)… (xn,yn) in the polynomial
f(x) = c1+ c2x + c3x2+ … + cn+1xn

The difference between Lagrange or Newton polynomials and polynomial fitting is that Lagrange/Newton polynomials will pass though all the points. In polynomial fitting, the polynomial may not pass through all the points but the error is minimized.

Lagrange/Newton approximation
Polynomial passes through the data
If we have m data points, we will obtain a polynomial of degree m+1

Polynomial fitting
Polynomial will pass close to the data points trying to reduce error.
Degree of polynomial is given in advance and it does not depend on number of points
The big disadvantage of polynomial approximation is that a large number of points will give a polynomial of larger degree i.e. it is difficult to evaluate and it can be imprecise/far away in values between points.



For example, if we want to fit to a quadratic
f(x) = Ax2 + Bx + C that fits the (x1,y1),(x2,y2)…(xn, yn)

E(A,B,C) = ∑(Axk2 + Bxk + C – yk)2

Minimize E(A,B,C)



Piecewise cubic splines

We want to fit N+1 data points (x0,y0), (x1,y1) … (xN,yN) where a=x0<x1<x2<…<xN = b into N cubic polynomials Sk(x) with coefficients Sk,0, Sk,1 , Sk,2 , Sk,3 that satisfies the following properties

Rule

 

I) S(x) = Sk(x) = Sk,0 + Sk,1(x-xk) + Sk,2(x-xk)2 + Sk,3(x-xk)3
for x = [xk,xk+1] , k = 0,1,…,N-1
We will us different cubic polynomials between each 2 points

II) S(xk) = yk for k = 0,1, … N .
The polynomial will pass through each point

III) Sk(xk) = Sk+1(xk+1) k = 0,1,…,N-2
A polynomial ends where a new one begins ( polynomials are connected)

IV) Sk’(xk+1) = Sk+1’(xk+1) k = 0,1,…,N-2
Continuity of its first derivative points.

V) Sk’’(xk+1) = Sk+1’’(xk+1)

Note: Polynomial approx builds a polynomial that passes through all points. If we have N+1 data points, the polynomial will be of degree N.

Polynomial fitting builds a polynomial that passes close to the N+1 data points. The error is minimum. The points may not touch the polynomials. Degree of polynomial < N . Polynomial can be of low degree

Cubic Splines
Instead of using single polynomials for all points, we use M polynomials for M+1 points
For each polynomial Sk(x) there are 4 coefficients to be determined Sk,0 , Sk,1, Sk,2 , Sk,3 (4 degrees of freedom) . There are N polynomials bringing the number of coefficients to be to be determined to 4N. From rule II, and with N+1 data points, we have N+1 constraints. From III, IV and V, each supplies N-1 constraints. So, they give 3(N-1) constraints.
N+1 + 3(N-1) = 4N-2. This leaves 2 additional degrees of freedom i.e. need 2 more equations. We will call them end point constraints and they involve S’(x) and S’’(x) at x0 and xN

Since S(x) is cubic and S’(x) is spline, we obtain equation of this line with lagrange.

S’’k(x) = S’’k(xk)(x-xk+1)/(xk-xk+1) + S’’k(xk+1)(x-xk)/(xk+1 – xk)
So we have that Sk’’(x) = S’’(xk) when x = xk and Sk’’(x) = S’’(xk+1) when x = xk+1

We are going to use mk = S’’(xk) and mk+1 = S’’(xk+1) and hk = mk+1 - mk

When we substitute Sk’’(x) = mk(xk+1 – xk)/hk + mk+1(x-xk)/hk

Sk(x) = ∫ ∫ Sk’’(x) = ∫ ∫ (xk+1 – xk)/hk + mk+1(x-xk)/hk dx

= ∫ -mk(xk+1 – xk)2/2hk + mk+1(x-xk)2/2hk + c1

= mk(xk+1-xk)3/6hk + mk+1(x-xk)3/6hk +c1x +c2

We can express c1x + c2 = Pk(xk+1 –x) + qk(x-xk)
Since Pk and qk are constants, we have

Sk(x) = mk(xk+1-xk)3/6hk + mk+1(x-xk)3/6hk + Pk(xk+1 –x) + qk(x-xk) - I

From rule II, each polynomial starts and ends in (xk, yk) and (xk+1, yk+1). So, substituting x = xk in I , we get

Sk(xk) = yk =  mk(x_k+1 - xk)³/6hk + Pk(x_k+1 - xk)                  ; hk = x_k+1 - xk

yk = mk.hk²/6 + Pk.hk

Pk = (yk - mk.hk²/6)/hk    - II

do same for x = x_k+1

Sk(xk+1) = yk+1 = m_k+1.hk²/6 + qk,hk

qk = (y_k+1 - m_k+1.hk²/6)/hk      - III

substituting II and III into I,

Sk(x) = mk(x_k+1 - xk)³/6hk + m_k+1(x-xk)³/6hk + (yk - mk.hk/6)(x_k+1 - xk) + y_k+1 - m_k+1.hk/6       - IV

This equation is closer to the spline that we are looking for. So, we need to find m_i , i = 0,1,...,N

To obtain these values, we differentiate IV

S'k(x) = -mk(x_k+1 - x)²/2hk + m_k+1(x - xk)²/2 - (yk/hk - mkhk/6) + y_k+1/hk - m_k+1hk/6           - V

From rule IV, we know tat k = k-1 substituted to V, will give S_k-1'(xk) = Sk'(xk)     . So, we evaluate V for Sk' and S_k-1'

Sk'(xk) = -mk.hk/3 + (Y_k+1 - Yk)/hk  - m_k+1.hk/6

lets say dk = (y_k+1 - yk)/hk

Sk'(xk) = -mk.hk/3 +dk - m_k+1.hk/6           - VI

We value V for S'k-1

S_k-1'(xk) =  mk.h_k-1/3  + m_k-1.hk-1/6 + (yk - y_k-1)/h_k-1      

d_k-1 = (yk - y_k-1)/h_k-1

S_k-1'(xk) = mk.h_k-1/3 + m_k-1.h_k-1/6 + d_k-1/h_k-1       - VII

From rule IV, S_k-1'(xk) = Sk'(xk)  i.e VI = VII

-mk.hk/3 - m_k+1.hk/6 + dk = mk.h_k-1/3 + m_k-1.h_k-1/6 + d_k-1 

6(dk - d_k-1) = 2.mk(hk + h_k-1 ) + m_k+1.hk + m_k-1.h_k-1

if uk = 6(dk- d_k-1)

m_k+1.hk + 2mk(hk + h_k-1) + m_k-1.h_k-1 = uk  ; this is for k = 1,2,...,N-1
This gives  N-1 equations to solve
m0,m1,...,mN ( N+1 variables). Hence we need 2 more equations. The other 2 equations will be given with conditions at the end points