Numerical Analysis (aimee)

We learned a more efficient way to code Newton’s divided difference in class last week. We began by trying the code with the equation with five equally spaced points:

f=@(x) sin(pi*x)

n=5

x=linspace(-1,1,n)’

c=zeros(n,1)

F=f(x)

c(1)=F(1)

for i=1:n-1

F = diff(F)./(x(i+1:n)-x(1:n-i))

c(i+1)=F(1)

end

m=101;

xx = linspace(-1,1,m)’;

phi = ones(m,n);

for i=2:n

phi(:,i) = phi(:,i-1).*(xx-x(i-1));

end

P=phi*c;

plot(xx,f(xx),‘-r’)

hold on

plot(xx,P,‘-ob’)

With the error:

Using the same function with Gaussian points:

With the error:

Next we used Gaussian points with the equation

%f=@(x) sin(pi*x)

f=@(x) 1./(1+4.*x.^2)

n=7

%x=linspace(-1,1,n)’

x = -cos(pi*(0:n-1)/(n-1))

c=zeros(n,1)

F=f(x)

c(1)=F(1)

for i=1:n-1

F = diff(F)./(x(i+1:n)-x(1:n-i))

c(i+1)=F(1)

end

m=101;

xx = linspace(-1,1,m)’;

phi = ones(m,n);

for i=2:n

phi(:,i) = phi(:,i-1).*(xx-x(i-1));

end

P=phi*c;

plot(xx,f(xx),‘-r’)

hold on

plot(xx,P,‘-ob’)

With the error:

Using the same function with equidistant points:

With the error:

I thought it would be interesting to do a little research into finding out some more information on the Runge Phenomenon. The Runge phenomenon occurs with polynomial interpolation for polynomials of high degree. The higher the order of the polynomial used, the larger the bound of error becomes, showing that instead of moving towards zero, the error is increasing.

Using equally spaced points, as the degree of the polynomial increases, the error near the endpoints gets increasingly larger. It can be proven that the error of the Runge phenomenon using equally spaced points will approach infinity as the degree of the polynomial increases. This contradicts the Weierstrass approximation theorem that explains that there exists a sequence of approximating polynomials with the error converging to zero.

The following picture demonstrates this:

The red curve is the Runge function.

The blue curve is a 5th order interpolating polynomial.

The green curve is a 9th order interpolating polynomial

At the interpolating points, the error between the function and the polynomial is zero. Between the interpolating points, especially near the endpoints, as the degree of the polynomial increases, the error increases.

The following table shows the Runge function, its derivatives, and their evaluations at

Since the magnitude of higher order derivatives of the Runge function gets larger, the bound for error for higher order polynomial interpolations is larger.

The Runge phenomenon shows that using equally spaced points with interpolation of higher degree polynomials can cause problems. Using Gaussian points, the error appears to be more regular, oscillations are minimized, and the error decreases as the degree of the polynomial interpolation increases.

Using the cubic spline interpolation method can help to avoid this since it uses piecewise polynomials.

When looking for a method that can be more useful in cases where the Lagrange polynomial may not work as well, Newton’s divided differences is one possibility. Starting with the interpolating polynomial:

Imposing the conditions:

And finding:

Etc…

One advantage this method has over the Lagrange method is that the degree of p can be increased from to by adding the term:

There would only need to be a calculation done for this specific added point, rather than in the Lagrange formula where every term would change.

Since the term depends only on , the nth divided difference can be expressed as

With some algebraic manipulation, the divided difference can be represented as:

for any value of n. Newton’s divided difference can be easier to use than the Vandermonde matrix method because it uses a lower or upper triangular matrix, rather than the more complicated Vandermonde matrix.

I found a website: http://www.math.ucla.edu/~ronmiech/YAN/ndivdiff.html that will generate the upper triangular matrix for Newton’s divided difference. Using this program and 5 values on the interval [-1,1] with the function sin(pi*x), I received the following

-1.0000   -0.5000         0        0.5000    1.0000

0       -1.0000         0        1.0000         0

-2.0000    2.0000    2.0000   -2.0000         0

4.0000         0       -4.0000         0             0

-2.6667   -2.6667         0             0             0

I tried to write a MATLAB code to compute the matrix. With some knowledge of MATLAB and compiling various sources of MATLAB help, I came up with the following code:

x=[-1 -.5 0 .5 1];

y=[sin(-pi) sin(-.5.*pi) 0 sin(0.5*pi) sin(pi)];

n=length(x);

c=zeros(n,n+1);

c(:,1)=x’;

c(:,2)=y’;

for j=3:n+1;

for i=j-1:n;

i1=i-j+2;

c(i,j)=( c(i,j-1)-c(i-1,j-1) )/( c(i,1)-c(i1,1) );

end;

end;

c’

This code yielded the following matrix

-1.0000   -0.5000         0    0.5000    1.0000

-0.0000   -1.0000         0    1.0000    0.0000

0   -2.0000    2.0000    2.0000   -2.0000

0         0    4.0000         0   -4.0000

0         0         0   -2.6667   -2.6667

The two matrices are relatively the same, except that the columns are in a different order, so it seems as though the code worked well.

The interpolation polynomial found using Newton’s method is the same as the one found using Lagrange method as well as Vandermonde method. As a theorem:

Given distinct points and corresponding values there exists a unique polynomial of degree in for which

In order to code the cubic spline interpolation in MATLAB, Katie found the book Numerical Computing with MATLAB by Cleve B. Moler in the library and it gave a good description on what to do.

The book explained to create the following m files:

function d = splineslopes(h, delta)

n = length(h) +1;
a= zeros(size(h)); b=a; c=a; r=a;

a(1:n-2) = h(2:n-1);
a(n-1) = h(n-2) + h(n-1);
b(1) = h(2);
b(2:n-1) = 2*(h(2:n-1)+h(1:n-2));
b(n) = h(n-2);
c(1) = h(1) + h(2);
c(2:n-1) = h(1:n-2);

r(1)=((h(1)+2*c(1))*h(2)*delta(1)+...
    h(1)^2*delta(2))/c(1);
r(2:n-1) = 3*(h(2:n-1).*delta(1:n-2)+...
    h(1:n-2).*delta(2:n-1));
r(n) = (h(n-1)^2*delta(n-2)+...
    (2*a(n-1)+h(n-1))*h(n-2)*delta(n-1))/a(n-1);

d = tridisolve(a,b,c,r);

function v = splinetx(x,y,u)
%SPLINETX  Textbook spline function.
%  v = splinetx(x,y,u) finds the piecewise cubic interpolatory
%  spline S(x), with S(x(j)) = y(j), and returns v(k) = S(u(k)).
%
%  See SPLINE, PCHIPTX.

%  First derivatives

   h = diff(x);
   delta = diff(y)./h;
   d = splineslopes(h,delta);

%  Piecewise polynomial coefficients

   n = length(x);
   c = (3*delta - 2*d(1:n-1) - d(2:n))./h;
   b = (d(1:n-1) - 2*delta + d(2:n))./h.^2;

%  Find subinterval indices k so that x(k) <= u < x(k+1)

   k = ones(size(u));
   for j = 2:n-1
      k(x(j) <= u) = j;
   end

%  Evaluate spline

   s = u - x(k);
   v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k)));

function x = tridisolve(a,b,c,d)
%   TRIDISOLVE  Solve tridiagonal system of equations.
%     x = TRIDISOLVE(a,b,c,d) solves the system of linear equations
%     b(1)*x(1) + c(1)*x(2) = d(1),
%     a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1,
%     a(n-1)*x(n-1) + b(n)*x(n) = d(n).
%
%   The algorithm does not use pivoting, so the results might
%   be inaccurate if abs(b) is much smaller than abs(a)+abs(c).
%   More robust, but slower, alternatives with pivoting are:
%     x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)
%     x = S\d where S = spdiags([[a; 0] b [0; c]],[-1 0 1],n,n)

x = d;
n = length(x);

for j = 1:n-1
   mu = a(j)/b(j);
   b(j+1) = b(j+1) - mu*c(j);
   x(j+1) = x(j+1) - mu*x(j);
end

x(n) = x(n)/b(n);
for j = n-1:-1:1
   x(j) = (x(j)-c(j)*x(j+1))/b(j);
end

Then typing in the MATLAB command window: 

>> n=5;
>> x=linspace(-1,1,n)';
y= sin(pi*x);
m=101;
xx = linspace(-1,1,m)';
>> P= splinetx(x,y,xx);
>> plot(xx, P)
>> y1 = sin(pi*xx);
>> plot(xx, y1)
>> plot(xx,y1,xx,P)
>> y2= abs(y1 - P);
>> plot(xx, y2)

We used the code and changed the functions and amount of points to generate the following graphs today in class:

sin(pi*x) with 5 points

error for sin(pi*x) with 5 points

sin(pi*x) w 10 points

error for sin(pi*x) w 10 points

error for sin(pi*x) w 100 points

error for sin(pi*x) w 1000 points

We tried to use the Runge function with Gaussian points and it was able to graph the polynomial interpolation but we couldn’t get it to graph the actual function.

Cubic Spline

DEFINITION: Suppose that are n+1 points, where . The function S(x) is called a cubic spline if there exists cubic polynomials with coefficients that satisfy the properties:

(i) for and

(ii) The spline passes through each data point: for

(iii) The spline forms a continuous function over [a,b]: for

(iv) The spline forms a smooth function: for

(v) The second derivative is continous:

Natural Spline: There exists a unique cubic spline with the free boundary conditions and

Clamped Spline: There exists a unique cubic spline with the first derivative boundary conditions and

Theorem: Minimum Property of Clamped cubic splines: Assume and is the unique clamped spline interpolant for which passes through and satisfies the clamped end conditions and . Then

Given the function f(x) = sin(pi*x) interpolate a polynomial for the interval [-1, 1]

syms t L P

n=5

x=linspace(-1,1,n)

Y=sin(pi*x)

for j=1:5

L(j)=1;

for i=1:5

if (i~=j)

L(j) = L(j)*(t-x(i))/(x(j)-x(i));

end

end

end

P=L*Y’

expand(P)

ans =

-27043212804868892242677135040851/10141204801825835211973625643008*t^3+108172851219475573938466140184915/40564819207303340847894502572032*t

To graph the polynomial versus the function:

x=[-1:.01:1]

y=(108172851219475573938466140184915*x)/40564819207303340847894502572032 – (27043212804868892242677135040851*x.^3)/10141204801825835211973625643008

y1=sin(pi*x)

plot(x,y,x,y1)

From the graph, it can be seen that the polynomial interpolation (which used 5 points) is the same as the actual function in exactly five places. At all other points it is very similar.

To graph the error:

y3=abs((108172851219475573938466140184915*x)/40564819207303340847894502572032 – (27043212804868892242677135040851*x.^3)/10141204801825835211973625643008-(sin(pi*x)))

plot(x,y3)

Given the function f(x) = exp(x) interpolate a polynomial for the interval [-1, 1]

syms t L P

n=5

x=linspace(-1,1,n)

Y=exp(x)

for j=1:5

L(j)=1;

for i=1:5

if (i~=j)

L(j) = L(j)*(t-x(i))/(x(j)-x(i));

end

end

end

P=L*Y’

expand(P)

ans =

48904249067013/1125899906842624*t^4+149756602621819/844424930131968*t^3+2250200748318013/4503599627370496*t^2+13481801331389657/13510798882111488*t+1

To graph the polynomial versus the function:

x=[-1:.01:1]

y=48904249067013/1125899906842624*t^4+149756602621819/844424930131968*t^3+2250200748318013/4503599627370496*t^2+13481801331389657/13510798882111488*t+1

y2=exp(x)

plot(x,y,x,y2)

Zooming in closer at a point that is not an interpolation point:

If zooming in close to an interpolation point, the graphs are the same.

It appears from far away that the polynomial interpolation and the function are the same, but after zooming it, it is clear they are very close but not the same.

When zooming in around zero, or one of the points that was used in interpolation, it appears that the graphs are the same, similar to what happened with the sin(pi*x) error from above.

When graphing the error:

y2=abs((48904249067013*x.^4)/1125899906842624 + (149756602621819*x.^3)/844424930131968 + (2250200748318013*x.^2)/4503599627370496 + (13481801331389657*x)/13510798882111488 + 1-(exp(x)))

plot(x,y2)

This proves that the error is zero at the points that were interpolated.

In an attempt to find examples that can cause the Vandermonde method or Lagrange method to fail, Katie and I made up various examples. The first example was chosen at random.

Example 1:

Using 6 points: (0, -5) (0.2, 0) (0.2, 0.1) (0.5, 1) (0.52, -0.1) (1,8)

Works in both Vandermonde and Lagrange

We then decided to try to make the interval between points smaller.

Example 2:

Using 4 points: (0,0) (0.001, 0.005) (0.002, 0.004) (0.005, -0.001)

Works in Vandermonde and Lagrange

In the previous example, we saw that there was no constant term. We thought this to be because the point (0,0) was chosen. Next, we decided to choose points that were very far apart.

Example 3:

Using 5 points: (-100, 89) (75, 200) (0,1) (1000,1000) (-20, 2)

Works in Vandermonde and Lagrange

After not being able to find any sort of error when choosing points that were close together then points that were far apart, we decided to try to choose points that had x values that were very close together and y values that were very far apart.

Example 4:

Using 3 points: (0.0001, 1000) (0.0002, -2) (0.0003, 200)

Worked in Vandermonde and Lagrange

Next, instead of working with points, we decided to use a specific function.

Example 5:

Y=cos(x) with 20 points chosen between -1 and 1

Works in Vandermonde and Lagrange

Since the cos function seemed to work efficiently in both methods, we thought to use the cos function of x to a negative power and add more points

Example 6:

n=50

x=linspace(-1,1,n)

Y=cos(x.^-1)

A=vander(x);

B=Y’;

Ax=B;

x=A\B

Result: Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.

In the Lagrange method, it seemed to compute ok, but took a very long time.

Overall, our results were inconclusive. We thought we may have found something that caused either one or both methods to fail (by putting x^-1) or increasing the amount of points but we were unable to find any concrete results.

When given various points in the x-y plane, the task to attempt to find a function that satisfies these points can sometimes be very difficult. One approach that can be taken uses the idea of the Vandermonde matrix and solving a linear system.

Suppose the polynomial is of the form . A linear system of the form Ax=B can be formed using the Vandermonde matrix of the points as A, the unknown coefficients as x, and the known function values as B.

A linear system of this form can be solved by hand using linear algebra, or much easier by using MATLAB. Taking the following example using data given in class:

A = vander(-1:0.5:1)

A =

1.0000   -1.0000    1.0000   -1.0000    1.0000

0.0625   -0.1250    0.2500   -0.5000    1.0000

0         0             0             0            1.0000

0.0625    0.1250    0.2500    0.5000    1.0000

1.0000    1.0000    1.0000    1.0000    1.0000

B=[0;-.5;-1;.75;-0]

B =

0

-0.5000

-1.0000

0.7500

0

x=A\B

x =

-4.6667

-1.6667

5.6667

1.6667

-1.0000

This data shows the polynomial which satisfies these points is of the form:

There is also another approach that can be used to find a polynomial which satisfies a set of points. This is known as Lagrange interpolation. Lagrange polynomials are interpolating polynomials which equal zero in at all points given, except one. For example: will equal one at and zero at

The polynomial interpolation can be defined as …

After computing the solution by hand, we received help with constructing the MATLAB code for the Lagrange interpolation and saw there were two different ways it could be done. The first used the code

syms t L P

x=[-1, -1/2, 0, 1/2, 1]

L(1) = ((t-x(2))*(t-x(3))*(t-x(4))*(t-x(5)))/((x(1)-x(2))*(x(1)-x(3))*(x(1)-x(4))*(x(1)-x(5)))

L(2) = ((t-x(1))*(t-x(3))*(t-x(4))*(t-x(5)))/((x(2)-x(1))*(x(2)-x(3))*(x(2)-x(4))*(x(2)-x(5)))

L(3) = ((t-x(2))*(t-x(1))*(t-x(4))*(t-x(5)))/((x(3)-x(2))*(x(3)-x(1))*(x(3)-x(4))*(x(3)-x(5)))

L(4) = ((t-x(2))*(t-x(3))*(t-x(1))*(t-x(5)))/((x(4)-x(2))*(x(4)-x(3))*(x(4)-x(1))*(x(4)-x(5)))

L(5) = ((t-x(2))*(t-x(3))*(t-x(4))*(t-x(1)))/((x(5)-x(2))*(x(5)-x(3))*(x(5)-x(4))*(x(5)-x(1)))

Y=[0; 3/4; -1; -1/2; 0]

P=L*Y

expand(P)

We received the answer:

P =

-2*(t+1)*t*(t-1/2)*(t-1)-4*(t+1/2)*(t+1)*(t-1/2)*(t-1)+4/3*(t+1/2)*t*(t+1)*(t-1)

ans =

-14/3*t^4+5/3*t^3+17/3*t^2-5/3*t-1

The second code, although more complex at first, will be better for future use.

syms t L P

x=[-1, -1/2, 0, 1/2, 1];

Y=[0; 3/4; -1; -1/2; 0];

for j=1:5

L(j)=1;

for i=1:5

if (i~=j)

L(j) = L(j)*(t-x(i))/(x(j)-x(i));

end

end

end

P=L*Y;

expand(P)

ans =

-14/3*t^4+5/3*t^3+17/3*t^2-5/3*t-1

It can be seen that both of these codes yield the same answer, which is also the same as the answer found using the Vandermonde matrix approach. The next step is to look more into these two approaches and find when one method will be more efficient than the other and whether or not both methods will work for all cases.

Introduction:

After understanding that Newton’s method generally has second order convergence, I wondered if it was possible to have 3^rd-order, 4^th-order, or even nth order of convergence.

A third order algorithm will triple the number of correct digits in each iteration. To find a third order iteration, one must attempt to determine an expression such that will lead to the desired 3^rd-order iterative formula. Reqriting as

The choice of letting will eliminate from the formula, which is required.

Iteration formula:

This yields the third-order formula which has the property

In all cases of third order iterations,

Since , the typical third-order behavior is

Example:

Using this third-order formula to estimate radical 2 with the following MATLAB code:

x(1)=1.0;
for i=1 : 6
    x(i+1) = sqrt(2) + ((x(i)-sqrt(2))^3)/(3*((x(i)))^2 +2)
end

x =
  Columns 1 through 3
   1.000000000000000   1.400000000000000   1.414213197969543
 Columns 4 through 6
   1.414213562373095   1.414213562373095   1.414213562373095
 Column 7
   1.414213562373095


Graphing the error,

Halley’s method:

The method of convergence with higher order than Newton’s method that is often looked at is Halley’s method. Halley is most famous for analyzing the orbit of the comet in 1682 and making the correct prediction of its return 76 years later.

This method uses the quadratic term of the taylor series expansion of f and requires that the second derivative is continuous.

is defined as the route of the quadratic expansion:

Each iteration would be of the form

Since there are two signs in the denominator, the sign is chosen to try to minimize the difference between and

The order of convergence of Halley’s method is three. Again, Halley’s function involves the evaluation of the second derivative and each iteration requires three function evaluations.

Thoughts:

It seems as though the third-order iteration formula appears to converge faster than the second-order formula and is in turn more efficient. Fewer iterations are needed to obtain the same accuracy. However, the third-order formula is much more complex than previous iteration formulas, almost entirely cancelling out its positives. In actuality, as the order of convergence increases, the complexity of the computation of derivatives, as well as increasingly large formulas, make the benefits of these higher order iteration methods not worth the trouble.

Introduction:

Newton’s method is generally of second order convergence. The number of correct digits will double with each iteration. Newton’s method is of the form:

There are some points that are necessary in order to ensure convergence of Newton’s method:

• The derivative must be calculated directly, not approximated by using two points and finding the slope.
• The initial point chosen must be sufficiently close to the actual root, if it is not fairly close, the method could diverge.
• The root must be a simple root, or of multiplicity one. If it is not, Newton’s method will be no faster than the bisection method would be.

Kantorovich?:

I wanted to look more into trying to find out if there was some sort of theorem or information that determines how close the initial guess must be to the actual root for Newton’s method to converge. It is clear that if any of the iterations of Newton’s method require a tangent line that is horizontal, the algorithm will not converge. I also found that the Kantorovich theorem can be used to describe this, however when looking more into it, I couldn’t really understand it.

Rate of Convergence of Newton’s Method:

To find the rate of convergence of Newton’s method:

• Assume that converges to
• Set for 
• Let two positive constants and exist
• 
• The sequence converges to with order of convergence .
• The number is called the asymptotic error constant
• If $latexR=1$, there is linear convergence
• If , there will be quadratic convergence

Simple Root versus Multiple Root:

Assume that Newton-Raphson iteration produces a sequence that converges to the root of the function

Simple Root: Convergence is quadratic and

for k sufficiently large.

Multiple root of order m: Convergence is linear and

for k sufficiently large.

Proof of Order of Convergence:

Expand the Taylor polynomial of degree n=1 around the point x=p_k

Since is a zero of , let

.

With some rearranging,

.

Assuming that for all near

Since , you can divide by it

Since

Rearranging to

Using the fact that the Newton iteration formula is and substituting:

.

Using and when is sufficiently large

Taking the absolute values:

One of the last things to do concerning fixed point methods was to take a look at the error produced with each iteration using the following code.

x(1)=1.0;

for i=1 : 7

x(i+1) = -sqrt(1/8)*(x(i)^2 – 2) + x(i)

y=sqrt(2)-x(i+1)

end

Using the equation with radical two as the fixed point, starting with an initial point of x₀=1, and allowing the constant to equal the optimal constant, which is –sqrt(1/8), the following graph demonstrates how quickly the error decreases. After just six iterations, the root is found, with the error = 0.

Regarding my previous post on types of convergence of the fixed point method, I would like to clarify some things. Let F’ be the derivative of the function and r be the actual root.

• If |F’(r)|<1, the iterations will converge, therefore the error will decrease
• If |F’(r)|>1, the iterations will diverge, therefore the error will increase

If there is convergence,

• The error will decrease monotonically if 0 ≤ F’(r) < 1
• The error will decrease in an oscillatory manner if -1 < F’(r) < 0

I made up the following example which demonstrates this.

I designed the fixed point method, beginning at point x=0 to find the roots of following equation:

x^² + x – 2 = 0

x=

x_n=

This shows that with the guess of an initial point sufficiently close to the root, it is expected that the iterations will converge to the root one in an oscillatory manner and will not converge to the root negative two.

The following MATLAB code confirms this:

x(1) = 0;

for i=1:30

x(i+1)=2*(x(i)+1)^-1

end

Columns 1 through 12

0 2.0000 0.6667 1.2000 0.9091 1.0476 0.9767 1.0118 0.9942 1.0029 0.9985 1.0007

Columns 13 through 24

0.9996 1.0002 0.9999 1.0000

This demonstrates the oscillatory convergence to the root one.

No matter which initial point is chosen, the code will not converge to the other root which is negative two. I attempted to chose negative three and negative one, but still both converge to the root one. This is appealing that the MATLAB code actually works and supports the data I previously found regarding types of convergence.