User:Egm6341.s10.Team4.riherd

Mark Riherd

pmriherd@gmail.com

7.6 - Hermite-Simpson Integration of Coupled ODE's

Given

This is a flight control problem. The dynamics of flight are well documented and will not be described in detail here.

General Method:

By using collocation to force our Hermite-Simpson interpolation to exactly solve the given ODE's at the mid point $\displaystyle z_{i+1/2}$, we find that our time stepping method takes the form

$\displaystyle \mathbf {z_{i+1}} -\mathbf {z_{i}} ={\frac {\tfrac {h}{2}}{3}}(\mathbf {f_{i}} +4\mathbf {f_{i+1/2}+\mathbf {f_{i_{1}}} } )$

Additionally, using our interpolation $\displaystyle \mathbf {f_{i+1/2}} ={\tfrac {1}{2}}(\mathbf {z_{i}} +\mathbf {z_{i=1}} +{\tfrac {h}{8}}(\mathbf {f_{i}} -\mathbf {f_{i+1}} ))$.

However, no exact solution exists in order to determine $\displaystyle \mathbf {z_{i+1}}$. Given $\displaystyle \mathbf {z_{i}}$, the Newton-Rhapson method can be employed in order to find $\displaystyle \mathbf {z_{i+1}}$. By converting $\displaystyle \mathbf {z_{i+1}} -\mathbf {z_{i}} ={\frac {h}{6}}(\mathbf {f_{i}} +4\mathbf {f_{i+1/2}+\mathbf {f_{i_{1}}} } )$ to $\displaystyle \mathbf {F(z_{i},z_{i+1})} =\mathbf {0}$, the Newton-Rhapson method takes the form

$\displaystyle \mathbf {z_{i+1}^{k+1}} =\mathbf {z_{i+1}^{k}} -{\begin{bmatrix}\mathbf {\frac {dF(z_{i},z_{i+1}^{k})}{dz_{i+1}}} \end{bmatrix}}^{-1}\mathbf {F(z_{i},z_{i+1}^{k})}$

where

$\displaystyle \mathbf {F(z_{i},z_{i+1})} =\mathbf {z_{i+1}} -\mathbf {z_{i}} -{\frac {h}{6}}(\mathbf {f_{i}} +4\mathbf {f_{i+1/2}+\mathbf {f_{i_{1}}} } )=0$

Equations of Motion:

Let

$\displaystyle {\begin{bmatrix}z\end{bmatrix}}={\begin{bmatrix}x\\y\\v\\\gamma \end{bmatrix}}$

The equations of motion that we are concerned with are

$\displaystyle {\begin{bmatrix}{\dot {z}}\end{bmatrix}}={\begin{bmatrix}{\dot {x}}\\{\dot {y}}\\{\dot {v}}\\{\dot {\gamma }}\end{bmatrix}}={\begin{bmatrix}vcos\gamma \\vsin\gamma \\{\frac {1}{m}}((T-D)cos\alpha -Lsin\alpha -gsin\gamma )\\{\frac {T-D}{m}}sin\alpha +{\frac {L}{mv}}cos\alpha -{\frac {g}{v}}cos\gamma \end{bmatrix}}$

The physical parameters for which we will solve using are

$\displaystyle D={\tfrac {1}{2}}C_{d}\rho v^{2}S_{ref}$

$\displaystyle L={\tfrac {1}{2}}C_{l}\rho v^{2}S_{ref}$

where

$\displaystyle C_{d}=a_{1}\alpha ^{2}+a_{2}\alpha +a_{3}$

$\displaystyle C_{l}=b_{1}\alpha +b_{2}$

$\displaystyle \rho =c_{1}y^{2}+c_{2}y+c_{3}$

$\displaystyle a_{1}=-1.9431{\tfrac {1}{rad^{2}}}$

$\displaystyle a_{2}=-0.1499{\tfrac {1}{rad}}$

$\displaystyle a_{3}=0.2359$

$\displaystyle b_{1}=21.9{\tfrac {1}{rad}}$

$\displaystyle b_{2}=0.0$

$\displaystyle c_{1}=3.312E-9{\tfrac {kg}{m}}$

$\displaystyle c_{2}=-1.142E-4{\tfrac {kg}{m^{2}}}$

$\displaystyle c_{3}=1.224{\tfrac {kg}{m^{3}}}$

The initial conditions are

$\displaystyle x(t=0)=0.0m$

$\displaystyle y(t=0)=30.0m$

$\displaystyle v(t=0)=272.0{\tfrac {m}{s}}$

$\displaystyle \gamma (t=0)=0.0rad$

The flight is controlled using the angle of attack and the thrust, which are shown below

Find

We are asked to determine x,y,v, and \gamma using the Hermite-Simpson/Newton-Rhapson Method and then verify this with MATLAB's ode solver ode45.

Solution

In performing this analysis, we see that our Jacobian, $\displaystyle |\mathbf {\tfrac {dF}{dz_{i+1}}} |$ can be described as

$\displaystyle \mathbf {\frac {dF(z_{i},z_{i+1})}{dz_{i+1}}} =\mathbf {\frac {dz_{i+1}}{dz_{i+1}}} -\mathbf {\frac {dz_{i}}{dz_{i+1}}} +{\frac {h}{6}}(\mathbf {\frac {df_{i}}{dz_{i+1}}} +4\mathbf {\frac {df_{i+1/2}}{dz_{i+1}}} +\mathbf {\frac {df_{i+1}}{dz_{i+1}}} )$

Using the interpolation for z_{i+1/2} as a function of z_{i} and z_{i+1} described above, this becomes

$\displaystyle \mathbf {\frac {dF(z_{i},z_{i+1})}{dz_{i+1}}} =\mathbf {\frac {dz_{i+1}}{dz_{i+1}}} -\mathbf {\frac {dz_{i}}{dz_{i+1}}} +{\frac {h}{6}}(\mathbf {\frac {df_{i}}{dz_{i+1}}} +4\mathbf {{\frac {df_{i+1/2}}{dz_{i+1/2}}}{\frac {dz_{i+1/2/2}}{dz_{i+1}}}} +\mathbf {\frac {df_{i+1}}{dz_{i+1}}} )$

In evaluating this equation, our Jacobian can be described as

$\displaystyle \mathbf {\frac {dF(z_{i},z_{i+1})}{dz_{i+1}}} =\mathbf {I} -{\frac {h}{6}}(3\mathbf {I} -{\frac {h}{2}}\mathbf {\frac {df_{i+1}}{dz_{i+1}}} )\mathbf {\frac {df_{i+1}}{dz_{i+1}}}$

Thus, as long as $\displaystyle {\tfrac {df}{dz}}$ can be calculated, the Jacobian can be found.

The obvious next step is to calculate $\displaystyle {\tfrac {df}{dz}}$. The full expansion of $\displaystyle {\tfrac {df}{dz}}$ can be seen to be

$\displaystyle {\frac {\mathbf {df} }{\mathbf {dz} }}={\begin{bmatrix}0&0&cos(\gamma )&-vsin(\gamma )\\0&0&sin(\gamma )&vcos(\gamma )\\0&{\frac {1}{m}}(-{\frac {dD}{dy}}cos\alpha -{\frac {dL}{dy}}cos\alpha )&{\frac {1}{m}}(-{\frac {dD}{dv}}cos\alpha -{\frac {dv}{dy}}cos\alpha &{\frac {1}{m}}(-gcos\gamma )\\0&{\frac {1}{mv}}({\frac {dD}{dy}}sin\alpha +{\frac {dL}{dy}}cos\alpha )&-{\frac {T}{mv^{2}}}sin\alpha -{\frac {1}{m}}{\frac {d}{dv}}({\frac {D}{v}})cos\alpha +{\frac {1}{m}}{\frac {d}{dv}}({\frac {L}{v}})cos\alpha +{\frac {g}{v^{2}}}cos\gamma &{\frac {g}{v}}sin\gamma \\\end{bmatrix}}$

In performing the evaluations we can see that the two forms of numerical integration produce a very similar output.

Source Codes

Main Program

% Numerical Methods
% Homework 7 - Problem 6
% Mark Riherd

% This problem will solve the optimal control problem described in class
% using a Forward Euler finite difference scheme

close all
clear all
clc

nn=2^12+1;
t_span=40;

% Constants
m=1005;
g=9.81;
S_ref=0.3376;
a1=-1.9431;
a2=-0.1499;
a3=0.2359;
b1=21.9;
b2=0;
c1=3.312E-9;
c2=-1.142E-4;
c3=1.224;

% Initial conditions
x(1)=0.0;
y(1)=30.0;
v(1)=272;
gamma(1)=0.0;

% Parameters put into a vector
p=[m;g;S_ref;a1;a2;a3;b1;b2;c1;c2;c3];

% Time span
tt=0:t_span/(nn-1):t_span;

dt=tt(2)-tt(1);
h=dt;

alpha=zeros(1,nn);
thrust=zeros(1,nn);

% Control vector
for i=1:nn
   if (tt(i)<21)
      alpha(i)=0.03;
      thrust(i)=6000;
   elseif (tt(i)<27)
       alpha(i)=0.09+(0.13-0.09)*(tt(i)-21)/(27-21);
       thrust(i)=6000;
   elseif (tt(i)<33)
       alpha(i)=-0.13+(-0.2+0.13)*(tt(i)-27)/(33-27);
       thrust(i)=1000;
   elseif (tt(i)<40.01)
       alpha(i)=-0.2+(-0.13+0.2)*(tt(i)-33)/(40-33);
       thrust(i)=6000;
   end
end

% Hermite-Simpson/Newton-Rhapson method

% First point
z(:,1)=[x(1);y(1);v(1);gamma(1)];

for i=1:(nn-1)

    % For first iteration of Newton-Rhapson method
     zzp1=z(:,i);

    % End points
    f_i=z_dot(z(:,i),p,alpha(i),thrust(i));
    f_i1=z_dot(zzp1,p,alpha(i+1),thrust(i+1));

    % Mid-point
    g=1/2*(zzp1+z(:,i))-h/8*(f_i-f_i1);
    f_ih=z_dot(g,p,(alpha(i)+alpha(i+1))/2.0,(thrust(i)+thrust(i+1))/2.0);

    for j=1:20
        ff=zzp1-z(:,i)-h/6*(f_i+4*f_ih+f_i1);

        dfdz_j=dfdz(zzp1,p,alpha(i+1),thrust(i+1));
        dffdz=eye(4)-h/6*dfdz_j*(3*eye(4)+h/2*dfdz_j);
        dffdz_inv=inv(dffdz);

        % New values for z
        zz=zzp1;
        zzp1=zz-dffdz*ff;

        % Calculate f's for next iteration
        % End points
        f_i=z_dot(z(:,i),p,alpha(i),thrust(i));
        f_i1=z_dot(zzp1,p,alpha(i+1),thrust(i+1));

        % Mid-point
        g=1/2*(z(:,i)+zzp1)+h/8*(f_i-f_i1);
        f_ih=z_dot(g,p,(alpha(i)+alpha(i+1))/2.0,(thrust(i)+thrust(i+1))/2.0);

    end

    z(:,(i+1))=zzp1;

end

x=z(1,:);
y=z(2,:);
v=z(3,:);
gamma=z(4,:);

% ODE45 portion
tt_span=[0 40];
z0=[0; 30; 272; 0];

[tout,yout] = ode45(@flight45,tt_span,z0);

figure
subplot(4,1,1), plot(tt,v,tout,yout(:,3))
title('Velocity')
xlabel('t (s)')
ylabel('v(t) (m/s)')
legend('Hermite-Simpson/Newton-Rhapson Method','ode45')
subplot(4,1,2), plot(tt,y,tout,yout(:,2))
title('Altitude')
xlabel('t (s)')
ylabel('y(t) (m)')
subplot(4,1,3), plot(tt,gamma,tout,yout(:,4))
title('Flight angle \gamma(t)')
xlabel('t (s)')
ylabel('\gamma(t) (rad)')
subplot(4,1,4), plot(x,y,yout(:,1),yout(:,2))
title('Flight Trajectory')
xlabel('x (m)')
ylabel('y(x) m')

figure
subplot(2,1,1), plot(tt,alpha)
title('Angle of Attack')
xlabel('t (s)')
ylabel('\alpha(t) (rad)')
axis([0 40 -0.2 0.2])
subplot(2,1,2), plot(tt,thrust)
title('Thrust')
xlabel('t (s)')
ylabel('T(t) (rad)')
axis([0 40 0 7000])

Calculate the Jacobian

function [J]=dfdz(z,p,alpha,T)

% Extracting parameters
m=p(1);
g=p(2);
S_ref=p(3);
a1=p(4);
a2=p(5);
a3=p(6);
b1=p(7);
b2=p(8);
c1=p(9);
c2=p(10);
c3=p(11);

% Extracting state variable
x=z(1);
y=z(2);
v=z(3);
gamma=z(4);

% Calculating intermediate terms and derivatives of L and D
Cd=a1*alpha^2+a2*alpha+a3;
Cl=b1*alpha+b2;
rho=c1*y^2+c2*y+c3;

dDdy=0.5*Cd*v^2*S_ref*(2*c1*y+c2);
dLdy=0.5*Cl*v^2*S_ref*(2*c1*y+c2);

dDdv=rho*Cd*v*S_ref;
dLdv=rho*Cl*v*S_ref;

dDovdv=.5*rho*Cd*S_ref;
dLovdv=.5*rho*Cl*S_ref;

% Calculate the Jacobian
% x_dot terms
J(1,1)=0.0;
J(1,2)=0.0;
J(1,3)=cos(gamma);
J(1,4)=-v*sin(gamma);

% y_dot terms
J(2,1)=0.0;
J(2,2)=0.0;
J(2,3)=sin(gamma);
J(2,4)=v*cos(gamma);

% v_dot terms
J(3,1)=0.0;
J(3,2)=1/m*(-dDdy*cos(alpha)-dLdy*sin(alpha));
J(3,3)=1/m*(-dDdv*cos(alpha)-dLdv*sin(alpha));
J(3,4)=(-g*cos(gamma));

% gamma_dot terms
J(4,1)=0.0;
J(4,2)=1/(m*v)*(-dDdy*sin(alpha)+dLdy*cos(alpha));
J(4,3)=-T/(m*v^2)*sin(alpha)-1/m*dDovdv*sin(alpha)+1/m*dLovdv*cos(alpha)+...
    g/v^2*cos(gamma);
J(4,4)=g/v*sin(gamma);

end

Calculate f

function [f] = z_dot(z,p,alpha,T)

% Extracting parameters
m=p(1);
g=p(2);
S_ref=p(3);
a1=p(4);
a2=p(5);
a3=p(6);
b1=p(7);
b2=p(8);
c1=p(9);
c2=p(10);
c3=p(11);

% Extracting state variable
y=z(2);
v=z(3);
gamma=z(4);

% Calculating intermediate terms and derivatives of L and D
Cd=a1*alpha^2+a2*alpha+a3;
Cl=b1*alpha+b2;
rho=c1*y^2+c2*y+c3;

D=.5*Cd*v^2*rho*S_ref;
L=.5*Cl*v^2*rho*S_ref;

f(1)=v*cos(gamma);
f(2)=v*sin(gamma);
f(3)=1/m*((T-D)*cos(alpha)-L*sin(alpha))-g*sin(gamma);
f(4)=(T-D)/(m*v)*sin(alpha)+L/(m*v)*cos(alpha)-g/v*cos(gamma);

f=f';

end

Function for ode 45 to use

function dy = flight45(t,y)

% Extracting parameters
% Constants
m=1005;
g=9.81;
S_ref=0.3376;
a1=-1.9431;
a2=-0.1499;
a3=0.2359;
b1=21.9;
b2=0;
c1=3.312E-9;
c2=-1.142E-4;
c3=1.224;

% alpha and thrust
if (t<21)
  alpha=0.03;
  thrust=6000;
elseif (t<27)
   alpha=0.09+(0.13-0.09)*(t-21)/(27-21);
   thrust=6000;
elseif (t<33)
   alpha=-0.13+(-0.2+0.13)*(t-27)/(33-27);
   thrust=1000;
elseif (t<40.01)
   alpha=-0.2+(-0.13+0.2)*(t-33)/(40-33);
   thrust=6000;
end

% Calculating intermediate terms and derivatives of L and D
Cd=a1*alpha^2+a2*alpha+a3;
Cl=b1*alpha+b2;
rho=c1*y(2)^2+c2*y(2)+c3;

D=.5*Cd*y(3)^2*rho*S_ref;
L=.5*Cl*y(3)^2*rho*S_ref;

dy = zeros(4,1);    % a column vector
dy(1) = y(3)*cos(y(4));
dy(2) = y(3)*sin(y(4));
dy(3) = 1/m*((thrust-D)*cos(alpha)-L*sin(alpha))-g*sin(y(4));
dy(4) = (thrust-D)/(m*y(3))*sin(alpha)+L/(m*y(3))*cos(alpha)-g/y(3)*cos(y(4));

7.12 - Manipulation of Elliptical Integration

Given

$\displaystyle C=a\int _{t=0}^{2\pi }{\sqrt {1+e^{2}cos^{2}t}}dt$

Find

We are asked to show that $\displaystyle C=a\int _{t=0}^{2\pi }{\sqrt {1+e^{2}cos^{2}t}}dt=4a\int _{t=0}^{\tfrac {\pi }{2}}{\sqrt {1+e^{2}sin^{2}\alpha }}d\alpha$

Solution

Let $\displaystyle t=cos^{-1}(sin\alpha )$ and $\displaystyle dt=cos\alpha {\frac {1}{\sqrt {1-sin^{2}\alpha }}}d\alpha =cos\alpha {\frac {1}{cos\alpha }}d\alpha =d\alpha$.

Putting this into our equation,

$\displaystyle C=a\int _{t=0}^{2\pi }{\sqrt {1+e^{2}sin^{2}\alpha }}d\alpha$

The function $\displaystyle {\sqrt {1+e^{2}sin^{2}\alpha }}$ is periodic over a period of $\displaystyle \pi$ and in this period the function is even around the point $\displaystyle \alpha ={\tfrac {\pi }{2}}$. Thus, by symmetry, this integral can be broken up into 4 parts of equal magnitude such that

$\displaystyle C=4a\int _{t=0}^{\tfrac {\pi }{2}}{\sqrt {1+e^{2}sin^{2}\alpha }}d\alpha$

6.6 - Inversion of Local Domain Transformation

Given

In using our Hermit interpolation we use the transformation from $\displaystyle t\epsilon [t_{i},t_{i+1}]$ to $\displaystyle s\epsilon [0,1]$ using the transformation $\displaystyle t(s)=(1-s)t_{i}+st_{i+1}$

Find

We are asked to invert this transformation from the form of $\displaystyle t=t(s)$ to $\displaystyle s=s(t)$

Solution

Inverting this transformation to find s in terms of t is a simple algebra problem.

$\displaystyle t=(1-s)t_{i}+st_{i+1}$

$\displaystyle t=t_{i}+s(t_{i+1}-t_{i})$

$\displaystyle s={\frac {t-t_{i}}{t_{i+1}-t_{i}}}$

5.6 - Evaluation of d_2i

Given

The hyperbolic cotangent is defined as $\displaystyle \ coth(x)={\frac {cosh(x)}{sin(x)}}={\frac {e^{x}+e^{-x}}{e^{x}-e^{-x}}}.$.

The function $\displaystyle \ f(x)=xcosh(x)=\sum _{i=0}^{\infty }{\frac {2^{2i}B_{2i}}{(2i)!}}x^{2i}$, where B_i are the Bernoulli numbers.

Find

We are asked to verify the values of d_i for i=2,4, and 6, where d_i is the constant $\displaystyle \ d_{2i}={\frac {-B_{2i}}{(2i)!}}$. This is done previously given the values d₂=-1/12, d₄=1/720 and d₆=-1/30240.

We are also asked to determine the values of d₈ and d₁₀.

Additionally, we are asked to compare these values of d_i to another relationship for d_i used in our error of the trapezoidal rule, where $\displaystyle \ d_{2i}={\frac {p_{2i}(1)}{2^{2i}}}$.

Solution

Expanding cosh(x) and sinh(x) in terms of the exponential function and then expanding those exponential functions as series around the point x=0, we see that

$\displaystyle \ f(x)=xcosh(x)=\sum _{i=1}^{\infty }-2^{2i}d_{2i}x^{2i}$

$\displaystyle \ cosh(x)=1+{\frac {1}{2}}x^{2}-{\frac {1}{24}}x^{4}+{\frac {1}{720}}x^{6}+{\frac {1}{40320}}x^{8}+{\frac {1}{3628800}}x^{10}+...$

$\displaystyle \ sinh(x)=x+{\frac {1}{6}}x^{3}+{\frac {1}{120}}x^{5}+{\frac {1}{5040}}x^{7}+{\frac {1}{362880}}x^{9}+{\frac {1}{39916800}}x^{11}+...$

Substituting all of this into the function $\displaystyle \ f(x)=xcosh(x)$, we see that

$\displaystyle \ f(x)={\frac {x+{\frac {1}{2}}x^{3}-{\frac {1}{24}}x^{5}+{\frac {1}{720}}x^{7}+{\frac {1}{40320}}x^{9}+{\frac {1}{3628800}}x^{11}...}{x+{\frac {1}{6}}x^{3}+{\frac {1}{120}}x^{5}+{\frac {1}{5040}}x^{7}+{\frac {1}{362880}}x^{9}+{\frac {1}{39916800}}x^{11}+...}}$

Performing the appropriate polynomial division (which will not be shown here due to length and no clear way to show long division in Latex), we see that

$\displaystyle \ f(x)=xcosh(x)=1+{\frac {1}{3}}x^{2}-{\frac {1}{45}}x^{4}+{\frac {2}{975}}x^{6}-{\frac {1}{4725}}x^{8}+{\frac {2}{93555}}x^{10}-...$

This series expansion of xcoth(x) in hand, we are able to determine the values of $\displaystyle \ 2^{2i}d_{2i}$ to be -

2i	d2i (calculated from f(x)=xcosh(x))	d2i (calculated using known Bernoulli numbers)	d2i (given)
0	1	1	1
2	-1/12	-1/12	-1/2
4	1/720	1/720	1/720
6	-1/30240	-1/30240	-1/30240
8	1/209600	1/209600	n/a
10	-1/95800320	-1/95800320	n/a

Additionally, we are asked to calculate d_2i using a definition given by the polynomial functions defined in our trapezoidal error analysis, where $\displaystyle \ d_{2i}={\frac {p_{2i}(1)}{2^{2i}}}$.

From analysis done in class and coefficients as given by Kessler, P. (Trapezoidal Rule Error, http://www.mechanicaldust.com/UCB/math128a/extrap.pdf, 2008), we know that

$\displaystyle \ p_{2}(t)=-{\frac {t^{2}}{2}}+{\frac {1}{6}},p_{2}(1)=-{\frac {1}{3}}$

$\displaystyle \ p_{4}(t)=-{\frac {t^{4}}{24}}+{\frac {t^{2}}{12}}-{\frac {7}{360}},p_{4}(1)={\frac {1}{45}}$

$\displaystyle \ p_{6}(t)=-{\frac {t^{6}}{720}}+{\frac {t^{4}}{144}}-{\frac {7t^{2}}{720}}+{\frac {31}{15120}},p_{6}(1)=-{\frac {2}{945}}$

$\displaystyle \ p_{8}(t)=-{\frac {t^{8}}{40320}}+{\frac {t^{6}}{4320}}-{\frac {7t^{4}}{8640}}+{\frac {31t^{2}}{30240}}-{\frac {127}{604800}},p_{8}(1)={\frac {1}{4725}}$

$\displaystyle \ p_{10}=-{\frac {t^{1}0}{3628800}}+{\frac {t^{8}}{241920}}-{\frac {7t^{6}}{259200}}+{\frac {31t^{4}}{362880}}-{\frac {127t^{2}}{1209600}}+{\frac {73}{3421440}},p_{10}(t)={\frac {-2}{93555}}$

With these coefficients in hand, we are able to calculate the values for d_2i using the formula as given above relating p_2i(1) and d_2i.

2i	d2i (calculated from f(x)=xcosh(x))	d2i (calculated using polynomial coefficients)
0	1	1
2	-1/12	-1/12
4	1/720	1/720
6	-1/30240	-1/30240
8	1/209600	1/209600
10	-1/95800320	-1/95800320

Comparing these values, we see that the formulas given to find d_2i are consistent with each other.

Problem 4.9

Given

We have been given various integration methods thus far throughout the class, we are now asked to compare their efficiency and accuracy.

Find

We are asked to compare the second column of the Rhomberg Table (T1(n)) to the Composite Simpson's Rule, and derive any relationship that exists.

We are asked to integrate the function $\displaystyle \ f(x)={\frac {exp(x)-1}{x}}$ over the interval $\displaystyle \ x\epsilon [0,1]$ using the following methods and to compare the computational cost as a function of run-time -

• Composite Trapezoidal (optimized)
• Composite Simpson's Rule
• Rhomberg Interpolation
• chebfun sum command - a function that uses Chebyshev polynomials to integrate
• the matlab functions trapz, quad, and clencurt.

Solution

For the first three methods, pre-existing codes written by the authors were used in order to evaluate the computational time required to perform the integration. For the Rhomberg Extrapolation, the computation time was determine by the time required to compute T0(n) plus the time to perform the extrapolation.

	Trapezoidal Rule (optimized)	Simpson's Rule	Rhomberg Extrapolation	Matlab (trapz)	Matlab (quad)	chebfun
Integration time	0.077443875638587	0.0000884298553269	0.002205027039298	0.110061293768379	0.029048107601096	0.401365075398550
Error	0.000000000620887	0.000000000103787	0.000000000103788 (T1)	0.000000000620881	0.0000000005025	$\displaystyle O(10^{-16})$
n	8192	256	512	8192	n/a	n/a

When the Simpson's Rule and T1 (Richardsen Extrapolation) are compared to each other, the results are identical out to all of the significant digits calculated, thus we believe that the two are equivalents processes. Furthermore,

$\displaystyle T_{1}(2n)={\frac {4T_{0}(2n)-T_{0}(n)}{4-1}}$

Breaking this into $\displaystyle T_{0}(n)={\frac {h}{2}}({\frac {1}{2}}f_{0}+f_{1}+f_{2}+...{\frac {1}{2}}f_{n}$ and $\displaystyle T_{0}(n/2)=h({\frac {1}{2}}f_{0}+f_{2}+f_{4}+...{\frac {1}{2}}f_{n}$ we can see that $\displaystyle T_{1}(n)={\frac {h}{6}}(f_{0}+4f_{1}+2f_{2}+4f_{3}+f_{n})$, which is equal to the Simpson's Rule for 2n points.

For the Trapezoidal Rule and Rhomberg Extrapolation

% Numerical Methods HW4
% Mark Riherd
%
% This is a program to estimate the integration of (exp(x)-1)/x from x=0 to
% x=1 using the composite trapezoidal rule and Richardsen Interpolation.

clear all
close all
format long

% From previous Taylor Series work
I_exact=1.317902151454404;

%% Composite Trapezoidal Series %%

tic

n=[1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192];

% For T(2)
h=1./n(1);
I_n(1)=h*(1/2*(1+(exp(1)-1)));

% For n=4,8,16,32,...
for i=2:length(n)
    h=1./n(i);

    I_n(i)=.5*I_n(i-1);

    for j=1:n(i)
        xx(j)=j/n(i);
    end

    % Interior points
    for j=1:n(i)/2
        % jj-th point pinpoints points out of the finer domain
        jj=(2*j-1);
        I_n(i)=h*((exp(xx(jj))-1)/xx(jj))+I_n(i);
    end
time_trap_opt(i)=toc

end

%% Rhomberg Interpolation %%
tic

T1=zeros(length(n),1);
T2=zeros(length(n),1);
T3=zeros(length(n),1);
T4=zeros(length(n),1);
T5=zeros(length(n),1);

 for i=1:(length(n)-1)
     T1(i+1)=(4*I_n(i+1)-I_n(i))/(4-1);
     T2(i+1)=(4^2*T1(i+1)-T1(i))/(4^2-1);
     T3(i+1)=(4^3*T2(i+1)-T2(i))/(4^3-1);
     T4(i+1)=(4^4*T3(i+1)-T3(i))/(4^4-1);
     T5(i+1)=(4^5*T4(i+1)-T4(i))/(4^5-1);
 time_rhom(i)=toc

 end

 I_n=I_n'
 T1=T1
 T2=T2
 T3=T3

 E0=abs(I_exact-I_n)
 E1=abs(I_exact-T1)
 E2=abs(I_exact-T2)
 E3=abs(I_exact-T3)
 E4=abs(I_exact-T4)
 E5=abs(I_exact-T5)

For the Simpson's Rule

% Numerical Methods HW1
% Mark Riherd
%
% This is a program to estimate the integration of (exp(x)-1)/x from x=0 to
% x=1 using the composite Simpson's rule.

clear all
close all

tic
I_exact=1.317902151454404;

n=[2,4,8,16,32,64,128,256];

for i=1:length(n)
    h=1./n(i);

    for j=1:(n(i)+1)
       xx(j)=(j-1)/n(i);
   end

   I_n(i)=h/3*(1+(exp(1)-1));

   for j=2:(n(i))
    if (mod(j,2)==0)
          I_n(i)=h/3*4*((exp(xx(j))-1)/xx(j))+I_n(i);
    else
          I_n(i)=h/3*2*((exp(xx(j))-1)/xx(j))+I_n(i);
    end
   end

   E_n(i)=I_n(i)-I_exact;

time_simp(i)=toc

end

I_n=I_n'
E_n=E_n'

For Matlab functions and chebshev

% Numerical Methods HW4
% Mark Riherd
%
% This is a program to estimate the integration of (exp(x)-1)/x from x=0 to
% x=1 using the Matlab tools.

clear all
close all

% trapz
I_exact=1.317902151454404;

n=[2,4,8,16,32,64,128,256,512,1024,2048,4096,8192];

for i=1:length(n)
    tic

    xx(1)=0;
    f(1)=1;

    for j=2:(n(i)+1)
       xx(j)=(j-1)/n(i);
       f(j)=(exp(xx(j))-1)/xx(j);
    end

    I_n_t(i)=trapz(xx,f);
time_trapz(i,1)=toc
end

I_n_t=I_n_t'
E_n_t=abs(I_exact-I_n_t)

% quad
tic
F = @(x)(exp(x)-1)./x;
Q = quad(F,0,1);
time_quad=toc

% chebyshev
tic
f=chebfun('(exp(x)-1)./x',[0,1]);
Q = sum(f);
E_n_C=I_exact-Q
time_cheb=toc

Problem 4.10

Find

We are asked to integrate the function $\displaystyle \ f(x)={\frac {1}{1+x^{2}}}$ over the interval $\displaystyle \ x\epsilon [-5,5]$ using the following methods and to compare the computational cost as a function of run-time -

• Composite Trapezoidal (optimized)
• Composite Simpson's Rule
• Rhomberg Interpolation
• chebfun sum command - a function that uses Chebyshev polynomials to integrate
• the matlab functions trapz, quad, and clencurt.

Solution

Operating in the same manner as above, only for a different $\displaystyle f(x)$ and bounds of integration.

	Trapezoidal Rule (optimized)	Simpson's Rule	Rhomberg Extrapolation	Matlab (trapz)	Matlab (quad)	chebfun
Integration time	1.026093406564399	0.000155742133262	0.002772121982166	1.494371285726102	0.027459449851913	0.376844048450762
Error	0.000000000229628	0.000000000163010	0.000000000163012 (T1)	0.000000000229616	0.0000000004441	$\displaystyle O(10^{-16})$
n	32768	256	512	32768	n/a	n/a

For the Trapezoidal and Rhomberg Extrapolation

% Numerical Methods HW4
% Mark Riherd
%
% This is a program to estimate the integration of 1/(1+x^2) from x=-5 to
% x=5 using the composite trapezoidal rule and Richardsen Interpolation.

clear all
close all
format long

% From analytical
I_exact=atan(5)-atan(-5);

%% Composite Trapezoidal Series %%
tic

n=[1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768];%,65536,131072,262144,524288,1048576,2097152,4194304,8388608,16777216];
x_span=10;
x_start=-5;

% For T(2)
h=x_span*1./n(1);
I_n(1)=h/2*(1/(1+(-5)^2)+1/(1+5^2));

% For n=4,8,16,32,...
for i=2:length(n)

    h=x_span./n(i);

    I_n(i)=.5*I_n(i-1);

    for j=1:n(i)-1
        xx(j)=(j*x_span)/n(i)+x_start;
    end

    % Interior points
    for j=1:n(i)/2
        % jj-th point pinpoints points out of the finer domain
        jj=(2*j-1);
        I_n(i)=h*(1/(1+xx(jj)^2))+I_n(i);
    end

    time_trap_opt(i)=toc

end

%% Richardson Interpolation %%
tic

T1=zeros(length(n),1);
T2=zeros(length(n),1);
T3=zeros(length(n),1);
T4=zeros(length(n),1);
T5=zeros(length(n),1);
T6=zeros(length(n),1);

 for i=1:(length(n)-1)
     T1(i+1)=(4*I_n(i+1)-I_n(i))/(4-1);
     T2(i+1)=(4^2*T1(i+1)-T1(i))/(4^2-1);
     T3(i+1)=(4^3*T2(i+1)-T2(i))/(4^3-1);
     T4(i+1)=(4^4*T3(i+1)-T3(i))/(4^4-1);
     T5(i+1)=(4^5*T4(i+1)-T4(i))/(4^5-1);
     T6(i+1)=(4^6*T5(i+1)-T5(i))/(4^6-1);

     time_rhom(i)=toc

 end

 T2(2)=0;
 T3(2)=0;
 T3(3)=0;

 I_n=I_n'
 T1=T1
 T2=T2
 T3=T3
 T4=T4
 T5=T5
 T6=T6

 E0=abs(I_exact-I_n)
 E1=abs(I_exact-T1)
 E2=abs(I_exact-T2)
 E3=abs(I_exact-T3)
 E4=abs(I_exact-T4)
 E5=abs(I_exact-T5)
 E6=abs(I_exact-T6)

For the Composite Simpson's Rule

% Numerical Methods HW4
% Mark Riherd
%
% This is a program to estimate the integration of 1/(1+x^2) from x=-5 to
% x=5 using the composite Simpson's Rule.

clear all
close all
format long

I_exact=atan(5)-atan(-5);

n=[2,4,8,16,32,64,128,256];%,32,64,128,256,512,1024,2048,4096,8192,16384,32768];
x_span=10;
x_start=-5;

for i=1:length(n)
    tic
    h=x_span./n(i);

    for j=1:(n(i)+1)
       xx(j)=x_span*(j-1)/n(i)+x_start;
    end

   I_n(i)=h/3*(1/(1+5^2)+1/(1+(-5)^2));

   for j=2:(n(i))
    if (mod(j,2)==0)
          I_n(i)=h/3*4*(1/(1+xx(j)^2))+I_n(i);
    else
          I_n(i)=h/3*2*(1/(1+xx(j)^2))+I_n(i);
    end
   end

   time_simp(i)=toc;

end

time_simp=time_simp'
I_n=I_n'
E_n=abs(I_exact-I_n)

For the Matlab functions

% Numerical Methods HW4
% Mark Riherd
%
% This is a program to estimate the integration of 1/(1+x^2) from x=-5 to
% x=5 using the Matlab tools.

clear all
close all

% trapz
I_exact=atan(5)-atan(-5);

n=[2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768];
x_span=10;
x_start=-5;

for i=1:length(n)
    tic

    xx(1)=-5;
    f(1)=1/(1+5^2);

    for j=2:(n(i)+1)
       xx(j)=x_span*(j-1)/n(i)+x_start;
       f(j)=1/(1+xx(j)^2);
    end

    I_n_t(i)=trapz(xx,f);
time_trapz(i,1)=toc
end

I_n_t=I_n_t'
E_n_t=abs(I_exact-I_n_t)

% quad
tic
F = @(x)(1./(1+x.^2));
Q = quad(F,-5,5,1.e-8);
E_n_Q=I_exact-Q

time_quad=toc

% chebyshev
tic
f=chebfun('(1./(1+x.^2))',[-5,5]);
Q = sum(f);
E_n_C=I_exact-Q
time_cheb=toc

Problem 3.last

Given

Rhomberg interpolation for the composite trapezoidal rule is defined as $\displaystyle \ T_{k}(n)={\frac {2^{2k}T_{k-1}(2n)-T_{k-1}(n)}{2^{2k}-1}}$

Find

We are asked to develop an efficient code that calculates and determines the efficiency of Rhomberg interpolation and compares it to the results to the Taylor Series integration method, using the function $\displaystyle \ f(x)={\frac {e^{x}-1}{x}}$.

Solution

The method for applying the trapezoidal rule has been well defined in previous work, thus it will not be shown again here.

The Rhomberg interpolation is described above, thus we will move right to out results.

In performing the Rhomberg interpolation, we found that with each additional interpolation of the composite trapezoidal series, an additional level of accuracy was obtained, with a linear return on the level of interpolation. However, even with this increased convergence to the exact solution, Rhomberg interpolation still does not converge nearly as quickly as the Taylor Series approximation of this integral does.

Thus, while the Rhomberg interpolation technique is helpful, the Taylor series method still converges much more quickly as a method of numerical integration.

Integration approximations

n	$I_{n,Taylor}$	$T_{0}(n)$	$T_{1}(n)$	$T_{2}(n)$	$T_{3}(n)$
2	1.250000000000000	1.328291727814890	1.317908916683140	1.318594340699574	1.317845576773824
4	1.315972222222222	1.320504619466078	1.318422984695466	1.318032767755261	1.317934391841460
8	1.317901815239985	1.318553086868629	1.318057156314024	1.317940540336073	1.317911725984127
16	1.317902151454404	1.318064905228940	1.317942362460728	1.317912176208376	1.317904655914551
32	1.317902151454404	1.317942841143417	1.317912294123424	1.317904685290699	1.317902784799460
64	1.317902151454404	1.317912323954498	1.317904692721200	1.317902786655409	1.317902310247431
128	1.317902151454404	1.317904694584294	1.317902787120757	1.317902310363742	1.317902191181285

Error Analysis

n	$E_{n,Taylor}$	$E_{0}(n)$	$E_{1}(n)$	$E_{2}(n)$	$E_{3}(n)$
2	0.151015657136614	0.010389576360486	0.000006765228736	0.000692189245170	0.000056574680580
4	0.004530469714098	0.002602468011674	0.000520833241061	0.000130616300857	0.000032240387056
8	0.000000832317334	0.000650935414225	0.000155004859619	0.000038388881668	0.000009574529722
16	$<O(10^{-15})$	0.000162753774536	0.000040211006324	0.000010024753971	0.000002504460147
32	$<O(10^{-15})$	0.000040689689013	0.000010142669019	0.000002533836294	0.000000633345056
64	$<O(10^{-15})$	0.000010172500094	0.000002541266795	0.000000635201004	0.000000158793027
128	$<O(10^{-15})$	0.000002543129890	0.000000635666352	0.000000158909337	0.000000039726880

Figure - Error of Various Integration Methods as a Function of the Number of Sub-Intervals Integrated Over

Code to run calculations -

% Numerical Methods HW3
% Mark Riherd
%
% This is a program to estimate the integration of (exp(x)-1)/x from x=0 to
% x=1 using the composite trapezoidal rule and Richardsen Interpolation.

clear all
close all

% From previous Taylor Series
I_exact=1.317902151454404;

I_Taylor=[1.250000000000000; 1.315972222222222; 1.317901815239985;...
          1.317902151454404; 1.317902151454404; 1.317902151454404];
E_Taylor=[0.151015657136614; 0.004530469714098; 0.000000832317334;...
          0.000000000000001; 0.000000000000000; 0.000000000000000];

%% Composite Taylor Series %%

% Values to work at
n=[2,4,8,16,32,64,128,256,512,1024,2048];

% Domain to work over - in order for code to work most efficiently
nn=n(length(n));
for j=1:(nn)
       xx(j)=j/nn;
end

for i=1:length(n)
    h=1./n(i);

    % First and last points
    I_n(i)=h/2*(1+(exp(1)-1));

    % Interior points
    for j=1:(n(i)-1)
        % jj-th point pinpoints points out of the finer domain
        jj=(j)*nn/n(i);
        I_n(i)=h*((exp(xx(jj))-1)/xx(jj))+I_n(i);
    end

    E_n(i)=I_n(i)-I_exact;
end

%% Richardson Interpolation %%

for i=1:(length(n)-1)
    % T1(n) Terms
    T1(i)=(2^(2*i)*I_n(i+1)-I_n(i))/(2^(2*i)-1);

   for j=1:i-1
       % T2(n) Terms
       T2(j)=(2^(2*j)*T1(j+1)-T1(j))/(2^(2*j)-1);

       for k=1:j-1
           % T3(n) Terms
           T3(k)=(2^(2*k)*T2(k+1)-T2(k))/(2^(2*k)-1);
       end
   end
end

I_n=I_n'
T1=T1'
T2=T2'
T3=T3'

E_n=E_n'
E1=I_exact-T1
E2=I_exact-T2
E3=I_exact-T3

% Plotting Error
loglog(n,E_n,n,[abs(E1);0],n,[abs(E2);0;0],n,[abs(E3);0;0;0],...
    n,[0;E_Taylor; 0; 0; 0; 0],n,10^-6*ones(length(n),1),':')
xlabel('n')
ylabel('Error(n)')
title('Error of Composite Trapezoial Rules and Rhomberg Interpolations')
legend('T0(n) - Error','T1(n) - Error','T2(n) - Error','T3(n) - Error','Taylor Series- Error')
axis([4 1024 10^-7.5 10^-1.5])

Problem 2.6

Given

We know that E(x) is defined as $\displaystyle \ E(x)=f(x)-f_{n}(x)$.

Find

We are asked to show that the $(n+1)^{th}$ derivative of our error $\displaystyle \ E(x)=f(x)-f_{n}(x)$ is $\displaystyle \ E^{(n+1)}(x)=f^{(n+1)}(x)$.

Solution

Assuming that $\displaystyle \ f_{n}(x)$ is n+1 times differentiable, then $\displaystyle \ E_{n}^{(n+1)}(x)=f^{(n+1)}(x)-f_{n}^{(n+1)}(x)$.

Let us assume because $\displaystyle \ f_{n}(x)\epsilon P_{n}$, $\displaystyle \ f_{n}(x)$ takes the form $\displaystyle \ f_{n}(x)=a_{n}x^{n}+a_{n-1}x^{n-1}+...+a_{1}x+a_{0}$.

We also know that $\displaystyle \ f_{n}(x)\epsilon P_{n}$. Taking the $n^{th}$ derivative of $\displaystyle \ f_{n}(x)$, we see that $\displaystyle \ f_{n}^{(n)}(x)=(n!)a_{n}$, a constant. Taking the $(n+1)^{th}$ derivative of $\displaystyle \ f_{n}(x)$, we see that $\displaystyle \ f_{n}^{(n+1)}(x)=0$.

Substituting this back into our equation for $\displaystyle \ E_{n}^{(n+1)}(x)$, we see that $\displaystyle \ E^{(n+1)}(x)=f^{(n+1)}-f_{n}^{(n+1)}(x)=f^{(n+1)}(x)-(0)=f^{(n+1)}(x)$.

Problem 2.12

Given

$\displaystyle \ \alpha (t)=\int _{-t}^{t}f(x(t'))dt'$ where $\displaystyle \ f(x(t))=F(t)$.

Find

Show that $\displaystyle \ \alpha (t)=F(-t)+F(t)$.

Solution

We are asked to show that for $\displaystyle \ \alpha (t)=\int _{-t}^{t}f(x(t'))dt'=\int _{-t}^{t}F(t')dt',\alpha ^{(1)}(t)=F(-t)+F(t)$, where $\displaystyle \ f(x(t))=F(t)$.

Using The Fundamental Theorem of Calculus (Second Part), we know that for some differentiable function $\displaystyle \ G(t)$ with with derivative $\displaystyle \ g(t)=G'(t)$, $\displaystyle \ \int _{a}^{b}g(t)dt=G(b)-G(a)$.

Allowing $\displaystyle \ f(x(t))=F(t)$ and $\displaystyle \ G'(t)=F(t)$, then $\displaystyle \ \alpha (t)=\int _{-t}^{t}F(t')dt'=G(t)-G(-t)$.

If the derivative of $\displaystyle \alpha$ is now taken, the $\displaystyle \ \alpha ^{(1)}(t)={\frac {d}{dt}}(G(t)-G(-t))=G'(t)+G'(-t)$.

Substituting in our definitions $\displaystyle \ F(t)=G'(t)$ we see that

$\displaystyle \ \alpha ^{(1)}(t)=F(-t)+F(t)$.

Homework #1

Problem 1.8

We are asked to evaluate the integral $\int _{0}^{1}f(x)dx$ for which $f(x)={\frac {e^{x}-1}{x}}$ using an approximation of $f(x)$ using Taylor Series, and also using the trapezoidal rule and Simpson's Rule.

Taylor Series

The function $e^{x}$ can be approximated using Taylor series to find that $e^{x}=\displaystyle \sum _{i=0}^{\infty }{\frac {x^{i}}{i!}}.$ Inserting this into our function, $f(x)$, and simplifying the expression, we find that $f(x)=\displaystyle \sum _{i=1}^{\infty }{\frac {x^{(i-1)}}{i!}}.$ While this form of our function, $f(x)$, presents an exact solution, we need to truncate this in order to solve it numerically. Therefore, let $f_{n}(x)$ be our approximation of $f(x)$ and let $f_{n}(x)=\displaystyle \sum _{i=1}^{n}{\frac {x^{(i-1)}}{i!}}.$ This function, $f_{n}(x)$ is able to be integrated analytically, when $\displaystyle \ \int _{0}^{1}{\frac {e^{x}-1}{x}}dx=\sum _{i=1}^{\infty }{\frac {x^{i}}{(i!)i}}$ the results of that are presented below in Table 1.8.1.

While $f_{n}(x)$ can be integrated numerically, it still has some degree of error because it is only an approximation of $f(x)$. The error in this approximation can be analyzed by integrating the remainder of $f_{n}(x)$, where $R\displaystyle \ _{n}(x)=f(x)-f_{n}(x)={\frac {x^{n}}{(n+1)!}}exp(\xi ),$ for some $\displaystyle \ \xi \epsilon [0,x]$.

Integrating the remainder of our function, $f(x)$, we come to see that our error is $\displaystyle \ E_{n}=I-I_{n}=\int _{0}^{1}[f(x)-f_{n}(x)]dx=\int _{0}^{1}{\frac {x^{n}}{(n+1)}}exp(\xi )dx$. This can be evaluated to find that $\displaystyle \ E_{n}=exp(\xi ){\frac {1}{(n+1)!(n+1)}}$. Thus, our maximum error using the Taylor series is $\displaystyle \ E_{n,max}=\max |exp(\xi )|{\frac {1}{(n+1)!(n+1)}}$. The results of this are also shown in Figure 1.8.1 below.

Trapezoidal Rule

The trapezoidal rule is a method of numerical integration that uses linear interpolation between two points in order to calculate the area under a given curve, $f(x)$.

For two points on a curve, where x=a,b. The trapezoidal rule can be expressed as $\displaystyle \ I_{1}=\int _{a}^{b}f_{n}(x)dx={\frac {b-a}{2}}(f(a)+f(b)).$

For n+1 number of points, the domain of $\displaystyle \ x\epsilon [a,b]$ can be broken up into n segments of equal length, h, and the trapezoidal rule works our such that $\displaystyle \ I_{n}={\frac {h}{2}}(f_{0}+2f_{1}+2f_{2}+...+2f_{n-1}+f_{n})$, where $\displaystyle \ f_{i}=f(x_{i})$.

For the case that $f(x)={\frac {e^{x}-1}{x}}$, results of this method are shown below in Table 1.8.1.

The error of this case, based on the result of using the Taylor series approximation of f(x) for a large value of n, is also available in Table 1.8.1 below.

Simpsons Rule

Whereas the trapezoidal rule uses a linear method of interpolation, Simpson's Rule uses a 2nd order polynomial to determine the area under a given curve, $f(x)$.

For three points on a curve, for x=a, (a+b)/2,b. $\displaystyle \ I_{2}=\int _{a}^{b}f_{n}(x)dx={\frac {b-a}{6}}(f(a)+4f({\frac {a+b}{2}})+f(b))$.

For n+1 number of points, the domain of $\displaystyle \ x\epsilon [a,b]$ can be broken into n segments of equal length, h, and Simpson's Rule works out such that $\displaystyle \ I_{n}={\frac {h}{6}}(f_{0}+4f_{1}+2f_{2}+4f_{3}+...2f_{n-2}+4f_{n-1}+f_{n})$, where $\displaystyle \ f_{i}=f(x_{i})$.

For the case that $f(x)={\frac {e^{x}-1}{x}}$, results of this method are shown below in Table 1.8.1.

The error of this case, based on the result of using the Taylor series approximation of f(x) for a large value of n, is also available in Table 1.8.1 below.

Table 1.8.1 - Integration approximations and error analysis

n	$I_{n,Taylor}$	$I_{n,trap}$	$I_{n,Simpson}$	$E_{n,Taylor}$	$E_{n,trap}$	$E_{n,Simpson}$
2	1.250000000000000	1.328291727814890	1.318008665676679	0.151015657136614	0.010389576360486	0.000106514222274
4	1.315972222222222	1.320504619466078	1.317908916683141	0.004530469714098	0.002602468011674	0.000006765228736
8	1.317901815239985	1.318553086868629	1.317902576002813	0.000000832317334	0.000650935414225	0.000000424548408
16	1.317902151454404	1.318064905228940	1.317902178015710	$<O(10^{-15})$	0.000162753774536	0.000000026561305
32	1.317902151454404	1.317942841143417	1.317902153114908	$<O(10^{-15})$	0.000040689689013	0.000000001660504
64	1.317902151454404	1.317912323954498	1.317902151558192	$<O(10^{-15})$	0.000010172500094	0.000000000103787
128	1.317902151454404	1.317904694584294	1.317902151460890	$<O(10^{-15})$	0.000002543129890	0.000000000006486

Problem 1.11

From interpolation theory, we know that $\displaystyle l_{i}=\prod _{j=0,j\neq i}^{n}{\frac {x-x_{j}}{x_{i}-x_{j}}}$.

Additionally, we define $\displaystyle P_{2}(x)=\sum _{i=0}^{2}l_{i}(x)f(x_{i})=f(x)$.

For $n=2$:

$\displaystyle l_{0}=\prod _{j=0,j\neq i}^{n=2}{\frac {x-x_{j}}{x_{i}-x_{j}}}={\frac {x-x_{1}}{x_{0}-x_{1}}}{\frac {x-x_{2}}{x_{0}-x_{2}}}$

$\displaystyle l_{1}=\prod _{j=0,j\neq i}^{n=2}{\frac {x-x_{j}}{x_{i}-x_{j}}}={\frac {x-x_{0}}{x_{1}-x_{0}}}{\frac {x-x_{2}}{x_{1}-x_{2}}}$

$\displaystyle l_{2}=\prod _{j=0,j\neq i}^{n=2}{\frac {x-x_{j}}{x_{i}-x_{j}}}={\frac {x-x_{0}}{x_{2}-x_{0}}}{\frac {x-x_{1}}{x_{2}-x_{1}}}$

Evaluating these functions $\displaystyle l_{0}(x),l_{1}(x)andl_{2}(x)$ at the interpolation points $\displaystyle x_{0},x_{1}$ and $\displaystyle x_{2}$, we find that:

$\displaystyle l_{0}(x_{0})=1,l_{0}(x_{1})=0$ and $\displaystyle l_{0}(x_{2})=0$

$\displaystyle l_{1}(x_{0})=0,l_{1}(x_{1})=1$ and $\displaystyle l_{1}(x_{2})=0$

$\displaystyle l_{2}(x_{0})=0,l_{2}(x_{1})=0$ and $\displaystyle l_{2}(x_{2})=1$

Reinserting these values of $\displaystyle l_{i}(x_{j})$ back into $\displaystyle P_{2}(x_{j})$, we can see that

$\displaystyle P_{2}(x_{0})=(1)f(x_{0})+(0)f(x_{1})+(0)f(x_{2})=f(x_{0})$

$\displaystyle P_{2}(x_{1})=(0)f(x_{0})+(1)f(x_{1})+(0)f(x_{2})=f(x_{1})$

$\displaystyle P_{2}(x_{2})=(0)f(x_{0})+(0)f(x_{1})+(1)f(x_{2})=f(x_{2})$

Thus, we can see that $\displaystyle P_{2}(x_{j})_{j=0,1,2}=\sum _{i=0}^{n=2}l_{i}(x_{j})f(x_{i})=f(x_{j})$.