User:Egm6341.s10.Team4.yunseok/HW5

Problem 5: Integration of logistic equation with $h=2^kh$

Given

Refer Lecture slide 42-2 for problem statement

The logistic equation

$\displaystyle \frac{d\overline{x}}{dt} = rx\left(1- \overline{x}\right)$

.

where,

$\displaystyle \overline{x} = \frac{x}{x_{max}}$

Find

Integrate logistic equation already found $\hat{h}$ such that Hermit-simpson algorithm yields error($10^{-6}$)

1) Run Hermit-Simpson with $\displaystyle h= 2^k \hat{h}$

2) Develop and Run Forward Euler with $\displaystyle h= 2^k \hat{h}$

3) Develop and Run Backward Euler with $\displaystyle h= 2^k \hat{h}$

Solution

1) Run Hermit-Simpson algo $\displaystyle h= 2^k \hat{h}$

See problem 3

2) Develop and Run Forward Euler with $\displaystyle h= 2^k \hat{h}$

$\displaystyle \begin{align} &\frac{d\overline{x}}{dt} = r\overline{x}\left(1- \overline{x}\right)\\ &\frac{\overline{x}_{i+1}-\overline{x}_i}{h}= r\overline{x_i}\left(1- \overline{x_i}\right)\\ &\overline{x}_{i+1}= hr\overline{x}_i\left(1- \overline{x}_i\right)+\overline{x}_i\\ &\overline{x}_{i+1}= (hr+1)\overline{x}_i-hr\overline{x}^2_i\\ \end{align}$ .

MATLAB CODE:

clc;
clear all;
k=1;
h=0.01*2^k;      %Step-size
r=1.2;
xmax=10;
t=0:h:10; %Range of time,t 
xbar_0 = 2/xmax;     %Initial Condition
tol = 1e-6;  %Absolute Tolerance
 
for i=1:length(t)-1;
    xbar(1)=xbar_0;
    xbar(i+1)= h*r*xbar(i)*(1-xbar(i))+xbar(i);
end
 x=xmax*xbar;
plot(t,double(x));                  %Plot the function
xlabel('Time(t)');
ylabel('Population(x)');
title('Population vs Time');
grid on;

3) Develop and Run Backward Euler with $\displaystyle h= 2^k \hat{h}$

$\displaystyle \begin{align} &\frac{d\overline{x}}{dt} = r\overline{x}\left(1- \overline{x}\right)\\ &\frac{\overline{x}_{i+1}-\overline{x}_i}{h}= r\overline{x}_{i+1}\left(1- \overline{x}_{i+1}\right)\\ &-\overline{x}_{i}= hr\overline{x}_{i+1}\left(1- \overline{x}_{i+1}\right)-\overline{x}_{i+1}\\ &=>hr\overline{x}^2_{i+1}+(1-hr)\overline{x}_{i+1}-\overline{x}_{i}= 0\\ \end{align}$ .

$\displaystyle \begin{align} &\overline{x}_{i+1}=\frac{hr-1 + \sqrt{(1-hr)^2+4hrx_i}}{2hr}\\ &\overline{x}_{i+1}=\frac{hr-1 - \sqrt{(1-hr)^2+4hrx_i}}{2hr}\\ \end{align}$ .

MATLAB CODE:

 
 
clc;
clear all;
k=1;  
h=0.01*2^k;      %Step-size
r=1.2;
xmax=10;
t=0:h:10; %Range of time,t 
xbar_0 = 2/xmax;     %Initial Condition
tol = 1e-6;  %Absolute Tolerance
 
for i=1:length(t)-1;
    xbar(1)=xbar_0;
    xbar(i+1)= ((h*r-1)+sqrt((1-h*r)^2+4*h*r*xbar(i)))/(2*h*r);
end
 x=xmax*xbar;
plot(t,double(x));                  %Plot the function
xlabel('Time(t)');
ylabel('Population(x)');
title('Population vs Time');
grid on;

Conclusion

As the step size($2^kh$) increase, the Forwaord Euler method shows the unstable properties.

on the other hands, The Hermit simpson and Backworkd Euler shows the relatively stable properties.

Problem 11: Circumference of Ellipse

Given

Lecture notes p.42-2

Arc length

$\displaystyle dl = [dx+dy]^{1/2}\,dt$

$\displaystyle x = a cost$

$\displaystyle y = b sint$

Find

Show that,
$\displaystyle C = \int dl = a \int_{0}^{2\pi} [1-e^2cos^2t]^{1/2}\,dt$
where,
Eccentricity, $\displaystyle e = \Bigg[1-\frac{b^2}{a^2}\Bigg]^{1/2}$

Solution

$\displaystyle \frac {dx}{dt} = -a sin t$

$\displaystyle \frac {dy}{dt} = b cos t$

and

$\displaystyle \begin{align} dl &= [dx^2 + dy^2 ]^{1/2}\\ &= \Bigg[\bigg(\frac{dx}{dt}\bigg)^2 + \bigg(\frac{dy}{dt}\bigg)^2 \Bigg]^{1/2} dt\\ &=\Bigg[(-a sin t)^2 + (b cos t)^2\Bigg]^{1/2} dt\\ &= \Bigg[(a^2 sin^2 t) + (b^2 cos^2 t)\Bigg]^{1/2} dt\\ &= a \Bigg[sin^2 t + \left(\frac{b^2}{a^2} cos^2 t\right)\Bigg]^{1/2} dt\\ &= a \Bigg[sin^2 t + cos^2 t + \frac{b^2}{a^2} cos^2 t- cos^2 t \Bigg]^{1/2} dt\\ &= a \Bigg[1 -\left(1-\frac{b^2}{a^2}\right) cos^2 t \Bigg]^{1/2} dt\\ &= a \Bigg[1 - e^2 cos^2 t \Bigg]^{1/2} dt\\ \end{align}$

.

So,

$\displaystyle C = \int_{0}^{2\pi} dl = a \int_{0}^{2\pi} [1-e^2cos^2t]^{1/2}\,dt$

Problem 11: Verification and comaprision of Problem 10 of HW5

Given

The given condition of No 10 problem of HW5

Refer Lecture slide 30-4 for problem statement Polar form relative to focus . Eq (3) from Page 30-3 $\displaystyle \ \ r(\theta) = \frac{a(1-e^2)}{1-e \, cos(\theta)}$ In our problem , a=1 where, eccentricity $\displaystyle e = \sin(\frac{\pi}{12})$ First Method to compute the Arc Length for an elliptical curve is using Eq (4) from Page 30-3 Since for the Ellipse $\displaystyle \theta$ range is from 0 to 2 $\displaystyle\pi$ Circumference $\displaystyle I( \theta) = \int_0^{2\pi} dl$ Where Eq (1) $\displaystyle dl= \sqrt{ {r}^{2}+ {\left( \frac{dr}{d \theta}\right)}^{2}} d\theta$. The second method uses the elliptic integral of the second kind which is defined as: $E(e^2) = \int_0^{\pi/2}\sqrt {1-e^2 \sin^2\theta}\ d\theta\$     Eq(6)   in page   30-4 The circumference of an ellipse is $\displaystyle C = 4 a E(e^2)$     Eq(5)   in page   30-4, Giving the integral: $\displaystyle I= 4 E(e^2) = 4\int_0^{\pi/2}\sqrt {1-e^2 \sin^2\theta}\ d\theta$ Find Arc length of ellipse using the two Integrals and the integration methods below. compute $\displaystyle I_n$ to the error $\displaystyle 10^{-10}$ , compute time using tic/toc matlab command and error estimate 1) Composit Trapozoidal rule 2) Romberg Table 3) Clencurt

.

Find

verify the result of pb10 of HW5(Team3) and compare with the own code.

Solution

first method $\displaystyle I( \theta) = \int_0^{2\pi} \sqrt{ {r}^{2}+ {\left( \frac{dr}{d \theta}\right)}^{2}} d\theta$

When run the Team 3's clencult code, they does'nt porvide answer(while loop continuously run.) the reason is that they use the I_exact value as Quad's solution(6.17660193898874),

but this value is not correct. when they calculate the Quad solution, they didn't designate the Toll error level(e.g. $\displaystyle 10^{-10}$). However their clencult code itself is correct,so when I change the I_exact value correctly, they provide correct answer.

<Romberg Table>

Team3

Own work

Romberg Table is correct.

Second method $\displaystyle I= 4 E(e^2) = 4\int_0^{\pi/2}\sqrt {1-e^2 \sin^2\theta}\ d\theta$

Team 3 only calculated the quad method.

<Romberg Table> Own work

In conclusion, Team 3's work is almost correct except for

1. quad function -- the error order$\displaystyle (10^{-8})$ was bigger than criterion $\displaystyle 10^{-10}$ because of not designating error tolerance level in the function

2. seconde method was not calculated except for quad function.

Also, when compare the result bwtween first and second mehtod, the results were same. and in the Romberg table, the second method was faster than first method.(see the upper table)

Matlab Code: method 1-Quad

clc;
clear all;
format long;
e=sin(pi/12); % Lecture Note 25-2
func=@(theta)((1-e.^2).^2./(1-e.*cos(theta)).^2+(1-e.^2)^2./(1-e.*cos(theta)).^4.*e.^2.*sin(theta).^2).^(1/2); % Defined dl from above
C=quad(func,0,2*pi,10^-11) % Integrate dl

Matlab Code: method 1-Composit Trap.Rule

% Trapezoidal Integration HW 5-10
% Modified from HW 4-11
clc
clear
format long
 
%limit
lowlimit=0.;
highlimit=2*pi;
 
%exact I (calculated by quad comment with Tol = 10^-11)
e=sin(pi/12);
F = @(x) sqrt(1-e^2*sin(x).*sin(x));
I = 4*quad(F,0,pi/2,10^-11);
 
tic
n=1;
E_n=1;
 
while E_n >= 10^(-10)
%make x(i)
x(1)=lowlimit;
h=(highlimit-lowlimit)/n;
for i=1:n;
    x(i+1)=x(i)+h;
end
 
%make I_n
sum=0.;
for i=2:n;
    sum=sum+ff5_10(x(i));
end
I_n=h*(sum+0.5*(ff5_10(x(1))+ff5_10(x(n+1))));
E_n=abs(I-I_n);
n=n+1;
end
n
toc

subfunction ff5_10

function y = ff5_10(x)
e=sin(pi/12);
r=(1-e^2)/(1-e*cos(x));
dr=-r*e*sin(x)/(1-e*cos(x));
%dr=-((1-e^2)*e)*sin(x)/((1-e*cos(x))^2);
y=sqrt(r^2+dr^2);
end

Matlab Code: method 1 - clencult

% ClenCurt spectral integration (modified from Trefethen p30.m)
% Modified from HW 4-11 By Roni Plachta for HW 5-10
% Computation: various values of N needed to get Error
% of integration < 1E-10 for (1-e^2*sin(x)^2)^.5 x=0 to x=pi/2
 
clear all
close all
format long
 
Nmax = 23; 
 
%limits
a=0;
b=2*pi;
 
%exact I (calculated by quad comment with Tol = 10^-11)
e=sin(pi/12);
F = @(x) sqrt(1-e^2*sin(x).*sin(x));
I = 4*quad(F,0,pi/2,1E-10);
 
tic
  for N = 1:Nmax; 
   [x,w] = clencurt(N);
 
   clear f;
   f(N+1,1)=0;
 
 % ReScale   
  x = 0.5 * ( ( x + 1.0 ) * a - ( x - 1.0 ) * b );
  w = 0.5 * ( b - a ) * w;
 
  for i = 1:N+1
        f(i,1)= ff5_10(x(i));
  end
 
  E(N)=0;
  S(N)=0;
  for i = 1:N+1
        S(N)= S(N)+ (w(i)*f(i));
  end
 
        E(N) = abs(I-S(N));
  end
 toc 
 
 xlswrite('5_10_7b.xls', E);
 S(N)

Matlab Code: method 1 - Romberg Table

clear all
close all
format long
clc
 
% integration Range
lowlimit=0.;
highlimit=2*pi;%pi/2;
%exact I (calculated by quad comment with Tol = 10^-11)
e=sin(pi/12);
F = @(x) sqrt(1-e^2*sin(x).*sin(x));
I = quad(F,0,pi/2,10^-11);
 
tic
E_n(1)=1;
k=1;
I_n(1)= ((highlimit-lowlimit)/2*(ff(lowlimit)+ff(highlimit)));
 
while E_n(k) >= 10^(-10)
    k=k+1;
    n=2^(k-1);
    x(1)=lowlimit;
    h=(highlimit-lowlimit)/n;
    for i=1:n/2;
        x(2*i)=x(1)+(2*i-1)*h;
    end
 
    %make I_n
    sum=0.;
    for j=1:n/2;
        sum=sum+ff(x(2*j));
    end
    I_n(k)=0.5*I_n(k-1) + h*sum;
%    E_n(k)=abs(I-I_n(k));
 
%% Richardson Interpolation %%
T1=zeros(length(n),1);
T2=zeros(length(n),1);
 
for i=1:k-1
     % T1(n) Terms
     T1(i+1)=(4*I_n(i+1)-I_n(i))/(4-1);
 end
 
 for i=1:k-2
     % T2(n) Terms
     T2(i+2)=(4^2*T1(i+2)-T1(i+1))/(4^2-1);
 end
 
  for i=1:k-3
     % T3(n) Terms
     T3(i+3)=(4^3*T2(i+3)-T2(i+2))/(4^3-1);
  end
 
  for i=1:k-4
     % T3(n) Terms
     T4(i+4)=(4^4*T3(i+4)-T3(i+3))/(4^4-1);
  end
 
    E_n(k)=abs(4*I-T2(k-1));
end
toc
k
n
I=4*I;
I_n=I_n;
T1=T1;
T2=T2;
T3=T3;
T4=T4;
Rombergtable=[I_n' T1' T2' T3' T4']

subfuntion ff:

function y = ff(x)
e=sin(pi/12);
%y=sqrt(1-e^2*(sin(x))^2);
y=((1-e.^2).^2./(1-e.*cos(x)).^2+(1-e.^2)^2./(1-e.*cos(x)).^4.*e.^2.*sin(x).^2).^(1/2);
end

<method 2>

Matlab Code: method 2 - Quad

clc
clear
e=sin(pi/12);
F = @(x) sqrt(1-e^2*sin(x).*sin(x));
tic
I = 4*quad(F,0,pi/2,10^-11)
toc

Matlab Code: method 2 - Compsit Trap. rule

clc
clear
 
%limit
lowlimit=0.;
highlimit=2*pi; %pi/2;
 
%exact I (calculated by quad comment with Tol = 10^-11)
e=sin(pi/12);
F = @(x) sqrt(1-e^2*sin(x).*sin(x));
I = 4*quad(F,0,pi/2,10^-11);
 
tic
n=1;
E_n=1;
while E_n >= 10^(-10)
%make x(i)
x(1)=lowlimit;
h=(highlimit-lowlimit)/n;
for i=1:n;
    x(i+1)=x(i)+h;
end
 
%make I_n
sum=0.;
for i=2:n;
    sum=sum+ff(x(i));
end
I_n=h*(sum+0.5*(ff(x(1))+ff(x(n+1))))
%I_n=4*h*(sum+0.5*(ff(x(1))+ff(x(n+1))))
E_n=abs(I-I_n);
n=n+1;
end
n
toc

subfunction ff

function y = ff(x)
e=sin(pi/12);
y=sqrt(1-e^2*(sin(x))^2);
%y=((1-e.^2).^2./(1-e.*cos(x)).^2+(1-e.^2)^2./(1-e.*cos(x)).^4.*e.^2.*sin(x).^2).^(1/2);
end

Matlab Code: method 2 - Clencurt

% Computation: various values of N needed to get Error
% of integration < 1E-10
 
clear
Nmax = 9; 
 
  %limit
a=0;
b=pi/2;
 
%exact I (calculated by quad comment with Tol = 10^-11)
e=sin(pi/12);
F = @(x) sqrt(1-e^2*sin(x).*sin(x));
I = 4*quad(F,0,pi/2,1E-10);
 
tic
  for N = 1:Nmax; i = 1:N;
    [x,w] = clencurt(N);
 
 % ReScale   
  x = 0.5 * ( ( x + 1.0 ) * a - ( x - 1.0 ) * b );
  w = 0.5 * ( b - a ) * w;
 
    f = 4*sqrt(1-e^2*(sin(x).*sin(x)));
    E(N) = abs(I-abs(w*f));
  end
 toc

Matlab Code: method 2 - RombergTable

clear all
close all
format long
clc
 
% integration Range
lowlimit=0.;
highlimit=pi/2;
%exact I (calculated by quad comment with Tol = 10^-11)
e=sin(pi/12);
F = @(x) sqrt(1-e^2*sin(x).*sin(x));
I = quad(F,0,pi/2,10^-11);
 
tic
E_n(1)=1;
k=1;
I_n(1)= ((highlimit-lowlimit)/2*(ff(lowlimit)+ff(highlimit)));
 
while E_n(k) >= 10^(-10)
    k=k+1;
    n=2^(k-1);
    x(1)=lowlimit;
    h=(highlimit-lowlimit)/n;
    for i=1:n/2;
        x(2*i)=x(1)+(2*i-1)*h;
    end
 
    %make I_n
    sum=0.;
    for j=1:n/2;
        sum=sum+ff(x(2*j));
    end
    I_n(k)=0.5*I_n(k-1) + h*sum;
%    E_n(k)=abs(I-I_n(k));
 
%% Richardson Interpolation %%
T1=zeros(length(n),1);
T2=zeros(length(n),1);
 
for i=1:k-1
     % T1(n) Terms
     T1(i+1)=(4*I_n(i+1)-I_n(i))/(4-1);
 end
 
 for i=1:k-2
     % T2(n) Terms
     T2(i+2)=(4^2*T1(i+2)-T1(i+1))/(4^2-1);
 end
 
  for i=1:k-3
     % T3(n) Terms
     T3(i+3)=(4^3*T2(i+3)-T2(i+2))/(4^3-1);
  end
 
  for i=1:k-4
     % T3(n) Terms
     T4(i+4)=(4^4*T3(i+4)-T3(i+3))/(4^4-1);
  end
 
    E_n(k)=abs(I-T2(k-1));
end
toc
k
n
I=4*I;
I_n=4*I_n;
T1=4*T1;
T2=4*T2;
T3=4*T3;
T4=4*T4;
Rombergtable=[I_n' T1' T2' T3' T4']

Problem 5: Identification of basis function $\displaystyle \overline{N_i}$

Given

Refer Lecture slide 35-2,4 for problem statement

The relationship between $\displaystyle Z(s)$ and $\displaystyle d_i$(degree of freedom)

$\displaystyle Z(s)=\sum_{i=0}^{3}c_i s^i = \sum_{i=1}^{4}N_i(s) d_i = \sum_{i=1}^{4} \overline{N_i}(s) \, \overline{d_i}$

$\displaystyle \longrightarrow (1)$

.

where,

$\displaystyle N_i(s)$ = basis function

$\displaystyle d_i(s)$ = degree of freedom

$\displaystyle d_1=Z_i$ $\displaystyle \,\,\,\,\, ,d_2=\dot{Z_i}$ $\displaystyle \,\,\,\,\, ,d_3=Z_{i+1}$ $\displaystyle \,\,\,\,\, ,d_4=\dot{Z}_{i+1}$

Find

Identify the the basis function $\displaystyle \left\{ N_i(s) \right\}\,\,\, ,i=1,2,3,4$, and plot $\displaystyle \overline{N_i}(s)$'s

Solution

From the equation (1) of 35-4(lecture slide)

$\displaystyle \begin{Bmatrix} c_0 \\ c _1 \\ c _2 \\ c_3 \end{Bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -3 & -2 & 3 & -1 \\ 2 & 1 & -2 & 1 \end{bmatrix} \begin{Bmatrix} Z_i \\ Z_i^' \\ Z_{i+1} \\ Z_{i+1}^' \end{Bmatrix} =\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -3 & -2 & 3 & -1 \\ 2 & 1 & -2 & 1 \end{bmatrix} \begin{Bmatrix} \overline{d}_1 \\ \overline{d}_2 \\ \overline{d}_3 \\ \overline{d}_4 \\ \end{Bmatrix}$

.

where,

$\displaystyle \overline{d}_2= h \, d_2$ $\displaystyle , \,\,\,\,\, \overline{d}_4= h \, d_4$ $\displaystyle , \,\,\,\,\, \overline{d}_1=d_1$ $\displaystyle , \,\,\,\,\, \overline{d}_3=d_3$

Expand the matrix equation.

$\displaystyle \begin{align} &c_0= \overline{d}_1 \\ &c_1= \overline{d}_2 \\ &c_2= -3\overline{d}_1 -2\overline{d}_2+3\overline{d}_3-\overline{d}_4 \\ &c_3= 2\overline{d}_1 +1\overline{d}_2-2\overline{d}_3+\overline{d}_4 \\ \end{align}$

.

plug in $\displaystyle c_i$ to equation (1)

$\displaystyle \begin{align} Z(s)&=\sum_{i=0}^{3}c_i s^i = c_0 + c_1 s^1 + c_2 s^2 + c_3 s^3 \\ &= \overline{d}_1 + \overline{d}_2 s^1 + (-3\overline{d}_1 -2\overline{d}_2+3\overline{d}_3-\overline{d}_4) s^2 + (2\overline{d}_1 +1\overline{d}_2-2\overline{d}_3+\overline{d}_4) s^3 \end{align}$

.

arrange in terms of $\displaystyle \overline{d}_i$

$\displaystyle \begin{align} & = (2s^3-3s^2+1)\overline{d}_1+(s^3-2s^2+s)\overline{d}_2+(-2s^3+3s^2)\overline{d}_3+(s^3-s^2)\overline{d}_4 \\ &=\overline{N}_1(s) \, \overline{d}_1+\overline{N}_2(s) \, \overline{d}_2+\overline{N}_3(s) \, \overline{d}_3+\overline{N}_4(s) \, \overline{d}_4 \end{align}$

.

So,

$\displaystyle \begin{align} & \overline{N}_1= 2s^3-3s^2+1 \\ & \overline{N}_2= s^3-2s^2+s \Rightarrow N_2= h(s^3-2s^2+s)\\ & \overline{N}_3= -2s^3+3s^2 \\ & \overline{N}_4= s^3-s^2 \Rightarrow N_4= h(s^3-s^2) \end{align}$

.

Plot of $\displaystyle N_i$

Problem 5: Derivation of Compsit Trap. Rule Error

Given

Refer Lecture slide 28-3 for problem statement

The Error of Composit Trap. Rule

$\displaystyle E^1_n= \sum_{k=0}^{n-1} \left[ \int_{x_{k}}^{x_{k+1}}f(x)dx-\frac{h}{2} \left\{ f(x_k)+f(x_{k+1}) \right\} \right]$

.

Find

Derive the Compsit Trap.Rule Error equation.

$\displaystyle E=\sum_{r=1}^{\ell} 2 \, \overline{d}_{2r} \, h^{2r-1}\Bigg[f^{(2r-1)}(b)-f^{(2r-1)}(a)\Bigg]-\frac{h^{2\ell}}{2^{2\ell}} \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}} \, p_{2\ell}(t_{k}(x)) \, f^{(2\ell)}(x)dx$

.

Solution

The Error of Composit Trap. Rule

$\displaystyle E^1_n= \sum_{k=0}^{n-1} \left[ \int_{x_{k}}^{x_{k+1}}f(x)dx-\frac{h}{2} \left\{ f(x_k)+f(x_{k+1}) \right\} \right]$

.

Transfer int.$\displaystyle [x_k, x_{k+1}]$ to $\displaystyle [1,+1]$

where,

$\displaystyle x(t):=t\frac{h}{2}+\frac{x_k+x_{k+1}}{2}, \, \, \, \, \, \, \, t \in [-1,1]\, \, \,\, \, \, , h:=(b-a)/n$

.

$\displaystyle g_k(t) = f(x(t)) \,\,\,\,\,\,such \,\,\,that\,\,\,\,\,x \in [x_k, x_{k+1}]$

.

$\displaystyle E^1_n= \frac{h}{2}\sum_{k=0}^{n-1} \left[ \int_{-1}^{+1}g_k(t)dt- \left\{ g_k(-1)+g_k(+1) \right\} \right]$

.

through the successive integration by parts, the below equation can be obtained21-2,326-3

$\displaystyle \begin{align} E^1_n & = \frac{h}{2} \left[ \sum_{k=0}^{n-1}\Big[p_2(t)g_k^{(1)}(t)+p_4(t)g_k^{(3)}(t_k)+p_6(t)g_k^{(5)}(t)+......+p_{2\ell}(t)g_k^{(2\ell-1)}(t)\Big]_{-1}^{+1}- \sum_{k=0}^{n-1} \left[ \int_{-1}^{1}p_{2\ell}(t)g_k^{(2\ell)}(t) dt \right]\right]\\ &= \frac{h}{2} \left[ \sum_{r=1}^{\ell}\sum_{k=0}^{n-1}\Big[p_{2r}(t)g_k^{(2r-1)}(t)\Big]_{-1}^{+1} - \sum_{k=0}^{n-1} \left[ \int_{-1}^{1}p_{2\ell}(t)g_k^{(2\ell)}(t) dt \right]\right] \end{align}$

.

Extend the first summation term.

$\displaystyle \begin{align} E^1_n &= \frac{h}{2} \left[ \sum_{r=1}^{\ell}\sum_{k=0}^{n-1}\Big[p_{2r}(+1)g_k^{(2r-1)}(+1)-p_{2r}(-1)g_k^{(2r-1)}(-1)\Big]- \sum_{k=0}^{n-1} \left[ \int_{-1}^{1}p_{2l}(t)g_k^{(2\ell)}(t) dt \right]\right]\\ & because \,\,\,\,\,\,p_{2r}(+1)=p_{2r}(-1),\\ &= \frac{h}{2} \left[ \sum_{r=1}^{\ell}\sum_{k=0}^{n-1} p_{2r}(+1) \Big[g_k^{(2r-1)}(+1)-g_k^{(2r-1)}(-1)\Big]- \sum_{k=0}^{n-1} \left[ \int_{-1}^{1}p_{2l}(t)g_k^{(2\ell)}(t) dt\right]\right] \end{align}$

.

Using below relationship,

$\displaystyle g_k^{(i)}(t)=\left(\frac{h}{2}\right)^if^{(i)}(x(t)),,\,\,\,\,\ x \in [x_k,x_{k+1}]$

.

$\displaystyle g_k(t)$can be changed to $\displaystyle f_(x_k)$ and below was obtained.

$\displaystyle \begin{align} &g_k^{(2r-1)}(+1)=\left(\frac{h}{2}\right)^{2r-1}f^{(2r-1)}(x_{k+1}),\\ &g_k^{(2r-1)}(-1)=\left(\frac{h}{2}\right)^{2r-1}f^{(2r-1)}(x_{k}),\\ &g_k^{(2\ell)}(t)=\left(\frac{h}{2}\right)^{2\ell}f^{(2\ell)}(x_{t}) \end{align}$

.

Plug in thess relationships to the equation $\displaystyle E_n^1$

$\displaystyle \begin{align} E^1_n & = \left(\frac{h}{2}\right) \left[ \sum_{r=1}^{\ell} p_{2r}(+1) \left(\frac{h}{2}\right)^{2r-1} \sum_{k=0}^{n-1} \Big[f^{(2r-1)}(x_{k+1})-f^{(2r-1)}(x_k)\Big]- \sum_{k=0}^{n-1} \left[ \int_{x_k}^{x_{k+1}}p_{2\ell}(t_k(x))\left(\frac{h}{2}\right)^{2\ell}f^{(2\ell)}(x) dx \right] \right]\\ &= \sum_{r=1}^{\ell} h^{2r} \frac{p_{2r}(+1)}{2^{2r}} \Big[f^{(2r-1)}(x_{n})-f^{(2r-1)}(x_0)\Big]- \left(\frac{h}{2}\right)^{2\ell+1} \sum_{k=0}^{n-1} \left[ \int_{x_k}^{x_{k+1}}p_{2\ell}(t_k(x))f^{(2\ell)}(x) dx \right] \end{align}$

.

by the definition of $\displaystyle \overline{d}_{2r}$

$\displaystyle \overline{d}_{2r} = \frac{p_{2r}(+1)}{2^{2r}}$

$\displaystyle x_n=b, x_0=a$

$\displaystyle E^1_n = \sum_{r=1}^{\ell} h^{2r} \overline{d}_{2r} \Big[f^{(2r-1)}(b)-f^{(2r-1)}(a)\Big]- \underbrace{\left(\frac{h}{2}\right)^{2\ell+1}}_{(?)} \sum_{k=0}^{n-1} \left[ \int_{x_k}^{x_{k+1}}p_{2\ell}(t_k(x))f^{(2\ell)}(x) dx \right]$

This solution is little bit different from the given solution of the lecture.

the coefficient of second summation term that marked $\displaystyle (?)$

we got $\displaystyle \left(\frac{h}{2}\right)^{2\ell+1}$ instead of $\displaystyle \left(\frac{h}{2}\right)^{2\ell}$

Problem 11: Derivation of the equation of arclength using the cosine law

Given

Refer Lecture slide [1] for problem statement

Figure 1

$\displaystyle r(\theta) = \frac{1-e^2}{1-e \, cos(\theta)}$

where, eccentricity $\displaystyle e = \sin(\frac{\pi}{12})$

Find

Derive the equation of arclength using triangle OAB and cosine law.

$d\ell=d\theta \left[ r^2+ \left( \frac{dr}{d\theta}\right)^2\right]^{\frac{1}{2}}$

Solution

[reference] the law of cosine $c^2 = a^2 + b^2 - 2ab\cos(\gamma)\,$

.

Apply the cosine law to triabgle OAB in the Figure1 (See Figure 1)

$\displaystyle dl^2=\overline{AB}^2 = \overline{OA}^2 + \overline{OB}^2 - 2\overline{OA}\,\,\, \overline{OB}\cos(d\theta)\,$

$\displaystyle dl = \sqrt{\overline{OA}^2 + \overline{OB}^2 - 2\overline{OA}\,\,\, \overline{OB}\cos(d\theta)}\,$

$\displaystyle dl = \sqrt{r(\theta)^2 + r(\theta+d\theta)^2 - 2r(\theta)r(\theta+d\theta)\cos(d\theta)}\,$

$\displaystyle dl = \sqrt{\underbrace{\Big[r(\theta) - r(\theta+d\theta)\Big]^2}_{(a)} + 2r(\theta)r(\theta+d\theta)\underbrace{(1-\cos(d\theta))}_{(b)}}\,$

Part$\displaystyle (a)$ is approximately

$\displaystyle (r(\theta) - r(\theta+d\theta))^2 \simeq (dr)^2 \,$

For part$\displaystyle (b)$

Using the other cosine summation law

$\displaystyle cos(\alpha + \beta)=cos(\alpha)cos(\beta)-sin(\alpha)sin(\beta) \,$

$\displaystyle cos(d\theta)=cos\left(\frac{d\theta}{2}+\frac{d\theta}{2}\right) =cos\left(\frac{d\theta}{2}\right)cos\left(\frac{d\theta}{2}\right)-sin\left(\frac{d\theta}{2}\right)sin\left(\frac{d\theta}{2}\right)= cos^2\left(\frac{d\theta}{2}\right)-sin^2\left(\frac{d\theta}{2}\right) \,$

because, $\displaystyle cos^2\left(\frac{d\theta}{2}\right)= 1-sin^2\left(\frac{d\theta}{2}\right)\,$

$\displaystyle = \left(1-sin^2\left(\frac{d\theta}{2}\right)\right)-sin^2\left(\frac{d\theta}{2}\right) = 1-2sin^2\left(\frac{d\theta}{2}\right)\ \,$

In small $\displaystyle d\theta$ approximately,

$\displaystyle sin(\frac{d\theta}{2}) \simeq \frac{d\theta}{2} \,$

using this fact,

$\displaystyle 1-cos(d\theta)= 1-\left(1-2sin^2\left(\frac{d\theta}{2}\right)\right) = 2\left(\frac{d\theta}{2}\right)^2 =\frac{(d\theta)^2}{2}$

So, Plug in part(a) and Part(b)

$\displaystyle dl = \sqrt{(r(\theta) - r(\theta+d\theta))^2 + 2r(\theta)r(\theta+d\theta)(1-\cos(d\theta))}\,$

$\displaystyle = \sqrt{(dr)^2 + 2r(\theta)r(\theta+d\theta)(\frac{(d\theta)^2}{2})} = \sqrt{(dr)^2 + r(\theta)r(\theta+d\theta)(d\theta)^2} \,$

$\displaystyle = d\theta \sqrt{\left(\frac{dr}{d\theta}\right)^2 + r(\theta)r(\theta+d\theta)} \,$

when $\displaystyle d\theta$ is very small, approximately, $\displaystyle r(\theta)r(\theta+d\theta) \simeq r(\theta)^2 \,$

$\displaystyle \simeq d\theta \sqrt{\left(\frac{dr}{d\theta}\right)^2 + r(\theta)^2} \,$