User:Egm6341.s10.Team4.riherd
Mark Riherd
pmriherd@gmail.com
7.6 - Hermite-Simpson Integration of Coupled ODE's
[edit | edit source]Given
[edit | edit source]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 , we find that our time stepping method takes the form
Additionally, using our interpolation .
However, no exact solution exists in order to determine . Given , the Newton-Rhapson method can be employed in order to find . By converting to , the Newton-Rhapson method takes the form
where
Equations of Motion:
Let
The equations of motion that we are concerned with are
The physical parameters for which we will solve using are
where
The initial conditions are
The flight is controlled using the angle of attack and the thrust, which are shown below
Find
[edit | edit source]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
[edit | edit source]In performing this analysis, we see that our Jacobian, can be described as
Using the interpolation for z_{i+1/2} as a function of z_{i} and z_{i+1} described above, this becomes
In evaluating this equation, our Jacobian can be described as
Thus, as long as can be calculated, the Jacobian can be found.
The obvious next step is to calculate . The full expansion of can be seen to be
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
[edit | edit source]Given
[edit | edit source]
Find
[edit | edit source]We are asked to show that
Solution
[edit | edit source]Let and .
Putting this into our equation,
The function is periodic over a period of and in this period the function is even around the point . Thus, by symmetry, this integral can be broken up into 4 parts of equal magnitude such that
6.6 - Inversion of Local Domain Transformation
[edit | edit source]Given
[edit | edit source]In using our Hermit interpolation we use the transformation from to using the transformation
Find
[edit | edit source]We are asked to invert this transformation from the form of to
Solution
[edit | edit source]Inverting this transformation to find s in terms of t is a simple algebra problem.
5.6 - Evaluation of d2i
[edit | edit source]Given
[edit | edit source]The hyperbolic cotangent is defined as .
The function , where Bi are the Bernoulli numbers.
Find
[edit | edit source]We are asked to verify the values of di for i=2,4, and 6, where di is the constant . This is done previously given the values d2=-1/12, d4=1/720 and d6=-1/30240.
We are also asked to determine the values of d8 and d10.
Additionally, we are asked to compare these values of di to another relationship for di used in our error of the trapezoidal rule, where .
Solution
[edit | edit source]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
Substituting all of this into the function , we see that
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
This series expansion of xcoth(x) in hand, we are able to determine the values of 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 d2i using a definition given by the polynomial functions defined in our trapezoidal error analysis, where .
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
With these coefficients in hand, we are able to calculate the values for d2i using the formula as given above relating p2i(1) and d2i.
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 d2i are consistent with each other.
Problem 4.9
[edit | edit source]Given
[edit | edit source]We have been given various integration methods thus far throughout the class, we are now asked to compare their efficiency and accuracy.
Find
[edit | edit source]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 over the interval 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
[edit | edit source]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 | |
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,
Breaking this into and we can see that , 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
[edit | edit source]Given
[edit | edit source]Find
[edit | edit source]We are asked to integrate the function over the interval 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
[edit | edit source]Operating in the same manner as above, only for a different 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 | |
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
[edit | edit source]Given
[edit | edit source]Rhomberg interpolation for the composite trapezoidal rule is defined as
Find
[edit | edit source]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 .
Solution
[edit | edit source]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.
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 |
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 | 0.000162753774536 | 0.000040211006324 | 0.000010024753971 | 0.000002504460147 | ||
32 | 0.000040689689013 | 0.000010142669019 | 0.000002533836294 | 0.000000633345056 | ||
64 | 0.000010172500094 | 0.000002541266795 | 0.000000635201004 | 0.000000158793027 | ||
128 | 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
[edit | edit source]Given
[edit | edit source]We know that E(x) is defined as .
Find
[edit | edit source]We are asked to show that the derivative of our error is .
Solution
[edit | edit source]Assuming that is n+1 times differentiable, then .
Let us assume because , takes the form .
We also know that . Taking the derivative of , we see that , a constant. Taking the derivative of , we see that .
Substituting this back into our equation for , we see that .
Problem 2.12
[edit | edit source]Given
[edit | edit source]where .
Find
[edit | edit source]Show that .
Solution
[edit | edit source]We are asked to show that for , where .
Using The Fundamental Theorem of Calculus (Second Part), we know that for some differentiable function with with derivative , .
Allowing and , then .
If the derivative of is now taken, the .
Substituting in our definitions we see that
.
Homework #1
[edit | edit source]Problem 1.8
[edit | edit source]We are asked to evaluate the integral for which using an approximation of using Taylor Series, and also using the trapezoidal rule and Simpson's Rule.
Taylor Series
The function can be approximated using Taylor series to find that Inserting this into our function, , and simplifying the expression, we find that While this form of our function, , presents an exact solution, we need to truncate this in order to solve it numerically. Therefore, let be our approximation of and let This function, is able to be integrated analytically, when the results of that are presented below in Table 1.8.1.
While can be integrated numerically, it still has some degree of error because it is only an approximation of . The error in this approximation can be analyzed by integrating the remainder of , where for some .
Integrating the remainder of our function, , we come to see that our error is . This can be evaluated to find that . Thus, our maximum error using the Taylor series is . 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, .
For two points on a curve, where x=a,b. The trapezoidal rule can be expressed as
For n+1 number of points, the domain of can be broken up into n segments of equal length, h, and the trapezoidal rule works our such that , where .
For the case that , 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, .
For three points on a curve, for x=a, (a+b)/2,b. .
For n+1 number of points, the domain of can be broken into n segments of equal length, h, and Simpson's Rule works out such that , where .
For the case that , 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.
n | |||||||
---|---|---|---|---|---|---|---|
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 | 0.000162753774536 | 0.000000026561305 | ||
32 | 1.317902151454404 | 1.317942841143417 | 1.317902153114908 | 0.000040689689013 | 0.000000001660504 | ||
64 | 1.317902151454404 | 1.317912323954498 | 1.317902151558192 | 0.000010172500094 | 0.000000000103787 | ||
128 | 1.317902151454404 | 1.317904694584294 | 1.317902151460890 | 0.000002543129890 | 0.000000000006486 |
Problem 1.11
[edit | edit source]From interpolation theory, we know that .
Additionally, we define .
For :
Evaluating these functions at the interpolation points and , we find that:
and
and
and
Reinserting these values of back into , we can see that
Thus, we can see that .