University of Florida/Egm6341/s10.Team2/HW6

Problem 1: Arc Length Computation and Comparisons

Statement

P.32-1 1. Compute the arc length to O(1e-10), using error analysis of composite trapezoidal rule, and verify the results by Romberg quadrature rule, Clenshaw-Curtis rule and chebfun/sum command.

2. Compute it again by means of the elliptical quadrature formula.

Solution

   clear;clc
   % function Definition
   syms e x
   r=(1-e^2)/(1-e*cos(x));
   drdx=diff(r,'x'); %Define dr/dx
   dl=(r^2+drdx^2)^(1/2); %Define dl
   d2ldx2=diff(dl,'x',2); %Define 2nd derivative of dl
   format long
   e=sin(pi/12);
   x=linspace(0,pi/3,2^14+1);
   d2l = -(1/4)./((1-e^2).^2./(1-e.*cos(x)).^2+(1-e^2).^2./(1-e.*cos(x)).^4.*e^2.*sin(x).^2).^(3/2).*(-2*(1-e^2).^2./(1-e.*cos(x)).^3.*e.*sin(x)-4*(1-e^2).^2./(1-e.*cos(x)).^5.*e^3.*sin(x).^3+2*(1-e^2).^2./(1-e.*cos(x)).^4.*e^2.*sin(x).*cos(x)).^2+(1/2)./((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).*(4*(1-e^2).^2./(1-e.*cos(x)).^4.*e^2.*sin(x).^2-2*(1-e^2).^2./(1-e.*cos(x)).^3.*e.*cos(x)+20*(1-e^2).^2./(1-e.*cos(x)).^6.*e^4.*sin(x).^4-20*(1-e^2).^2./(1-e.*cos(x)).^5.*e^3.*sin(x).^2.*cos(x)+2*(1-e^2).^2./(1-e.*cos(x)).^4.*e^2.*cos(x).^2);
   %plot(x,d2l)
   M2=max(abs(d2l));
   n_theory=ceil((pi^3*M2*1e10./(27*12)).^(1/2))
   n_computed=length(x)-1

   %Composite Trapezoidal
   [I,step]=ctrap('((1-sin(pi/12)^2)^2/(1-sin(12/pi)*cos(x))^2+(1-sin(pi/12)^2)^2/(1-sin(pi/12)*cos(x))^4*sin(pi/12)^2*sin(x)^2)^(1/2)',0,pi/3,1e-4)

   %Romberg Quadrature
   e=sin(pi/12);
   z = @(x) ((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);
   I_romb=rombquad(z,0,pi/3,10)

   %Clenshaw-Curtis
   e=sin(pi/12);
   [xx,w]=fclencurt(2^10+1,0,pi/3);
   yy=((1-e.^2).^2./(1-e.*cos(xx)).^2+(1-e.^2)^2./(1-e.*cos(xx)).^4.*e.^2.*sin(xx).^2).^(1/2);
   I_cc=w'*yy

   %Chebfun/Sum
   yy = chebfun('((1-sin(pi/12).^2).^2./(1-sin(pi/12).*cos(x)).^2+(1-sin(pi/12).^2)^2./(1-sin(pi/12).*cos(x)).^4.*sin(pi/12).^2.*sin(x).^2).^(1/2)',[0 pi/3]);
   II=sum(yy)

   %Elliptical Integral of second kind, using the adaptive simpson's rule
   clear;clc
   format long
   e=sin(pi/12);
   b=@(x) sqrt(1-e.^2.*sin(x).^2);
   I=quad(b,0,pi/3,1e-10)

   %function ctrap
   function [I,step] = ctrap(f,a,b,eps)
   n=1;
   h=(b-a)/2;
   I1=0;
   I2=(subs(sym(f),findsym(sym(f)),a) + subs(sym(f),findsym(sym(f)),b))*h;
   while abs(I2-I1)>eps
   n=2*n;
   h=(b-a)/n;
   I1=I2;
   I2=0;
      for i=0:n-1
          x=a+h*i;
          x1=x+h;
          I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+...
                       subs(sym(f),findsym(sym(f)),x1));
      end
   end
   I=I2;
   step=n;

   %function rombquad
   function I_romb = rombquad(f,a,b,m)
   h = 2.^((1:m)-1)*(b-a)/2^(m-1);             % These are the intervals used.
   k1 = 2.^((m-2):-1:-1)*2+1;                  % Index into the intervals.
   f1 = feval(f,a:h(1):b);                     % Function evaluations.
   R = zeros(1,m);                             % Pre-allocation.
   % Define the starting vector:
   for ii = 1:m
   R(ii) = 0.5*h(ii)*(f1(1)+2*...
   sum(f1(k1(end-ii+1):k1(end-ii+1)-1:end-1)) + f1(end));
   end
   % Interpolations:
   for jj = 2:m
   jpower = (4^(jj-1)-1);
   for kk = 1:(m-jj+1)
   R(kk) = R(kk)+(R(kk)-R(kk+1))/jpower;
   end
   end
   I_romb = R(1);

   %function fclencurt
   function [x,w]=fclencurt(N1,a,b)
   N=N1-1; bma=b-a;
   c=zeros(N1,2);
   c(1:2:N1,1)=(2./[1 1-(2:2:N).^2 ])'; c(2,2)=1;
   f=real(ifft([c(1:N1,:);c(N:-1:2,:)]));
   w=bma*([f(1,1); 2*f(2:N,1); f(N1,1)])/2;
   x=0.5*((b+a)+N*bma*f(1:N1,2));

Results Discussion

1) For error analysis, the n required for O(1e-10) is 16756.

2) For actual performace of trapezoidal rule, n is 16384.

3) For the first method, I = 1.263690485847868, I_romb = 1.263690485847868, I_cc = 1.263690485847865, II = 1.263690485847865.

4) For the second method, I = 1.036826643842352. This is because theta is defined differently in both cases.

Author

Egm6341.s10.team2.lee 20:58, 7 April 2010 (UTC)

Problem 2: Derivation of the Linear Momentum Term in the Y' Direction

Statement

P.33-2

At an instant t+dt, consider v+dv to show that

$dP_{y'}=mVd\gamma \quad$

after neglecting the higher order terms.

Solution

The Above figure represents the original axes and the new axes. It also shows the vector representations of the velocities at time=t and time=t+dt.

Our goal is to identify the following:

$dPy'=mV_{y'}$

The velocity at time t+dt is represented by V+dV. This means an infinitesimally small amount of time has elapsed. With this small change in time the angle $d\gamma$ is considered to be small.

Using trigonometry to express the component of V+dV that is in the direction of Y', we write the following

$V_{y'}=\sin(d\gamma )[V+dV]$

Using Taylor Series expansion to approximate the sine of the small angle (as described in Small-Angle Approximation), we can identify that for small values of the angle:

$\sin(d\gamma )=d\gamma$

Using this approximation the following can be written:

$V_{y'}=d\gamma [V+dV]=d\gamma V+d\gamma dV$

By ignoring the higher order term of the previous result we are able to show:

$dPy'=mV_{y'}=m(Vd\gamma )$

Author

Solved by--Guillermo Varela (GV)

---

Problem 3: Computation of the co-efficient matrix to express $C_{i}$ interms of the degrees of freedom

Statement

P.35-3

From the discussion in the above mentioned slides it we have the following equation:

${\begin{bmatrix}A_{1,1}&A_{1,2}&A_{1,3}&A_{1,4}\\A_{2,1}&A_{2,2}&A_{2,3}&A_{2,4}\\A_{3,1}&A_{3,2}&A_{3,3}&A_{3,4}\\A_{4,1}&A_{4,2}&A_{4,3}&A_{4,4}\end{bmatrix}}{\begin{Bmatrix}c_{0}\\c_{1}\\c_{2}\\c_{3}\end{Bmatrix}}={\begin{Bmatrix}Z_{i}\\Z_{i}'\\Z_{i+1}\\Z_{i+1}'\end{Bmatrix}}$

The first two rows of matrix A, is given. Calculate the rows 3 and 4.So that the co-efficients can be expressed as below:

${\begin{Bmatrix}c_{0}\\c_{1}\\c_{2}\\c_{3}\end{Bmatrix}}={\begin{bmatrix}A_{1,1}&A_{2,1}&A_{3,1}&A_{4,1}\\A_{1,2}&A_{2,2}&A_{3,2}&A_{4,2}\\A_{1,3}&A_{2,3}&A_{3,3}&A_{4,3}\\A_{1,4}&A_{2,4}&A_{3,4}&A_{4,4}\end{bmatrix}}{\begin{Bmatrix}Z_{i}\\Z_{i}'\\Z_{i+1}\\Z_{i+1}'\end{Bmatrix}}$

Solution

From the discussion in P.35-2 and P.35-3 we have the following equations:

$\displaystyle Z(s)=\sum _{i=0}^{3}c_{i}s^{i}=\sum _{i=1}^{3}N_{i}(s)d_{i}$

(2 p35-2)

$\displaystyle d_{1}=Z_{i}$ $\displaystyle d_{3}=Z_{i+1}$ $\displaystyle d_{2}={\dot {Z_{i}}}$ $\displaystyle d_{4}={{\dot {Z}}_{i+1}}$

(3 p35-2)

$\displaystyle {\dot {Z}}={\frac {dZ}{dt}}=\underbrace {\frac {dZ}{ds}} _{Z'}.{\frac {ds}{dt}}$

(4 p35-2)

Using ${\frac {ds}{dt}}={\frac {1}{h}}$ and the above equations we have shown the first two rows of the matrix A to be ${\begin{bmatrix}1&0&0&0\\0&1&0&0\end{bmatrix}}$

Using the equations above we have

$\displaystyle d_{3}=Z_{i+1}=Z(s=1)=c_{0}+c_{1}+c_{2}+c_{3}$

(1)

$\displaystyle d_{4}={\dot {Z}}_{i+1}={\frac {1}{h}}Z'(s=1)=c_{1}+2c_{2}+3c_{3}$

(2)

Putting the results in matrix form we have eqn 7 P.35-3.

${\begin{bmatrix}1&0&0&0\\0&1&0&0\\1&1&1&1\\0&1&2&3\end{bmatrix}}{\begin{Bmatrix}c_{0}\\c_{1}\\c_{2}\\c_{3}\end{Bmatrix}}={\begin{Bmatrix}Z_{i}\\Z'_{i}\\Z_{i+1}\\Z'_{i+1}\end{Bmatrix}}$

Author

--Egm6341.s10.team2.niki 21:45, 5 April 2010 (UTC)

--Proofread by Egm6341.s10.team2.lee 00:38, 6 April 2010 (UTC)

Problem 4: Inverse of Matrix

Statement

P. 35-4

Solution

The Inverse of the Matrix $A={\begin{bmatrix}1&0&0&0\\0&1&0&0\\1&1&1&1\\0&1&2&3\end{bmatrix}}\;\;$ can be found out by doing row operations as follows:

The augmented matrix of A and I is represented by

$X={\begin{bmatrix}1&0&0&0&|\;1&0&0&0\\0&1&0&0&|\;0&1&0&0\\1&1&1&1&|\;0&0&1&0\\0&1&2&3&|\;0&0&0&1\end{bmatrix}}$

By carrying out row transformation: $R_{3}\rightarrow 3R_{3}-R_{4}$ on the augmented matrix 'X', it becomes:

$X={\begin{bmatrix}1&0&0&0&|\;1&0&0&0\\0&1&0&0&|\;0&1&0&0\\3&2&1&0&|\;0&0&3&-1\\0&1&2&3&|\;0&0&0&1\end{bmatrix}}$

By carrying out row transformation: $R_{3}\rightarrow R_{3}-3R_{1}-2R_{2}$ on the augmented matrix 'X', it becomes:

$X={\begin{bmatrix}1&0&0&0&|\;+1&0&0&0\\0&1&0&0&|\;+0&1&0&0\\0&0&1&0&|\;-3&-2&3&-1\\0&1&2&3&|\;+0&0&0&1\end{bmatrix}}$

By carrying out row transformation: $R_{4}\rightarrow R_{4}-R_{2}-2R_{3}$ on the augmented matrix 'X', it becomes:

$X={\begin{bmatrix}1&0&0&0&|\;+1&0&0&0\\0&1&0&0&|\;+0&1&0&0\\0&0&1&0&|\;-3&-2&3&-1\\0&0&0&3&|\;+6&3&-6&3\end{bmatrix}}$

Finally, by carrying out row transformation: $R_{4}\rightarrow R_{4}/3$ on the augmented matrix 'X', it becomes:

$X={\begin{bmatrix}1&0&0&0&|\;+1&0&0&0\\0&1&0&0&|\;+0&1&0&0\\0&0&1&0&|\;-3&-2&3&-1\\0&0&0&1&|\;+2&1&-2&1\end{bmatrix}}$

Here the Left part of the augmented matrix is reduced to a unit matrix $\;I_{(4X4)}$ by carrying out the row operations and as a result we obtained a transformed matrix on the right part of the augmented matrix which is $\;A^{-1}$. Therefore,

$A^{-1}={\begin{bmatrix}+1&0&0&0\\+0&1&0&0\\-3&-2&3&-1\\+2&1&-2&1\end{bmatrix}}$

The same result can be verified by using the Matlab code below and 'B' is the value of $\;A^{-1}$ in the code.

Matlab code to verify the inverse of matrix- 'A':

clear;
clc;
A=[1,0,0,0;0,1,0,0;1,1,1,1;0,1,2,3];
B=inv(A);
syms Zi Zip Zi1 Zip1;
Z=[Zi;Zip;Zi+1;Zip+1];
C=B*Z;

Author

Solved by--Srikanth Madala (SM)

---

Problem 5:Identify basis functions

Statement

P. 35-4

To identify the set of basis functions: $\;N_{1}(s),N_{2}(s),N_{3}(s)\;and\;N_{4}(s)$ and plot them as a function of 's' varying between $\;s=0\;and\;s=1$

Solution

Using the $\;A^{-1}$ calculated in the above problem, we can find out the values of the coefficients- $\;C_{0},C_{1},C_{2}\;and\;C_{3}$ in terms of $Z_{i},Z'_{i},Z_{i+1}\;and\;Z'_{i+1}$ by using the following equation:

${\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}}$

Therefore,
${\begin{Bmatrix}C_{0}=Z_{i}\\C_{1}=Z'_{i}\\C_{2}=-3Z_{i}-2Z'_{i}+3Z_{i+1}-1Z'_{i+1}\\C_{3}=2Z_{i}+1Z'_{i}-2Z_{i+1}+1Z'_{i+1}\end{Bmatrix}}$

But, we know that:

$Z\left(s\right)=\sum _{i=1}^{4}N_{i}\left(s\right)d_{i}=\sum _{i=0}^{3}C_{i}s^{i}$

$N_{1}\left(s\right)d_{1}+N_{2}\left(s\right)d_{2}+N_{3}\left(s\right)d_{3}+N_{4}\left(s\right)d_{4}=C_{0}s^{0}+C_{1}s^{1}+C_{2}s^{2}+C_{3}s^{3}$

By definition: $d_{1}:=Z_{i},\;d_{2}:=Z'_{i},\;d_{3}:=Z_{i+1}\;and\;d_{4}:=Z'_{i+1}$ and also we already know the values of the coefficients- $\;C_{0},C_{1},C_{2}\;and\;C_{3}$ in terms of $Z_{i},Z'_{i},Z_{i+1}\;and\;Z'_{i+1}.$ Therefore, plugging these values in the above equation, gives us:

$N_{1}\left(s\right)Z_{i}+N_{2}\left(s\right)Z'_{i}+N_{3}\left(s\right)Z_{i+1}+N_{4}\left(s\right)Z'_{i+1}=Z_{i}+Z'_{i}s+\left[-3Z_{i}-2Z'_{i}+3Z_{i+1}-Z'_{i+1}\right]s^{2}+\left[2Z_{i}+Z'_{i}-2Z_{i+1}+Z'_{i+1}\right]s^{3}$

$N_{1}\left(s\right)Z_{i}+N_{2}\left(s\right)Z'_{i}+N_{3}\left(s\right)Z_{i+1}+N_{4}\left(s\right)Z'_{i+1}=Z_{i}\left(1-3s^{2}+2s^{3}\right)+Z'_{i}\left(s-2s^{2}+s^{3}\right)+Z_{i+1}\left(3s^{2}-2s^{3}\right)+Z'_{i+1}\left(-s^{2}+s^{3}\right)$

By comparing LHS and RHS, we can deduce the following:

$N_{1}\left(s\right)=\left(1-3s^{2}+2s^{3}\right)$ $N_{2}\left(s\right)=\left(s-2s^{2}+s^{3}\right)$ $N_{3}\left(s\right)=\left(3s^{2}-2s^{3}\right)$ $N_{4}\left(s\right)=\left(-s^{2}+s^{3}\right)$

Matlab code to generate the following plots:

clear;
clc;
syms s;
N1=1-(3*(s^2))+(2*(s^3));
N2=s-(2*(s^2))+(s^3);
N3=(3*(s^2))-(2*(s^3));
N4=-(s^2)+(s^3);
subplot(2,2,1); ezplot(N1,[0,1])
subplot(2,2,2); ezplot(N2,[0,1])
subplot(2,2,3); ezplot(N3,[0,1])
subplot(2,2,4); ezplot(N4,[0,1])

Author

Solved by--Srikanth Madala (SM)
Proofread by--Egm6341.s10.team2.lee 02:13, 8 April 2010 (UTC)

Problem 6: Transformation of Variables

Statement

P.36-1 Express s in terms of t

Solution

Since

$t(s)=t_{i}+{\frac {t_{i+1}-t_{i}}{1-0}}S=t_{i}+hS$

Solving for S,

$s(t)={\frac {t-t_{i}}{h}}$

Verify:

$s(t_{i})=0\quad and\quad s(t_{i+1})=1$

Author

Egm6341.s10.team2.lee 00:10, 5 April 2010 (UTC)

Problem 7:Enforcing compliance with the DE at the point (s=1/2)

Statement

P. 36-2 By enforcing compliance with the DE at the point $Z_{i+1}=Z(s=1/2)$ show that

$\displaystyle Z_{i+1/2}=Z(s=1/2)={\frac {1}{2}}[Z_{i}+Z_{i+1}]+{\frac {h}{8}}[f_{i}-f_{i+1}]$

Solution

We have from eqn 1 P. 35-4 and Prob 4 of HW6 [1] the following

$\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}}$

(1 p 35-4)

In eqn 2 P. 35-2 at the point $s={\frac {1}{2}}$ we get

$\displaystyle Z(s={\frac {1}{2}})=Z_{i}+Z'_{i}({\frac {1}{2}})+[-3Z_{i}-2Z'_{i}+3Z_{i+1}-Z_{i+1}]({\frac {1}{4}})+[2Z_{i}+Z'_{i}-2Z_{i+1}+Z'_{i+1}]({\frac {1}{8}})$

(1)

$\displaystyle Z(s={\frac {1}{2}})=Z_{i}[1-{\frac {3}{4}}+{\frac {1}{4}}]+Z'_{i}[{\frac {1}{2}}-{\frac {1}{2}}+{\frac {1}{8}}]+Z_{i+1}[{\frac {3}{4}}-{\frac {1}{4}}]+Z'_{i+1}[{\frac {-1}{4}}+{\frac {1}{8}}]$

(2)

$\displaystyle Z(s={\frac {1}{2}})=Z_{i}[{\frac {1}{2}}]+Z'_{i}[{\frac {1}{8}}]+Z_{i+1}[{\frac {1}{2}}]+Z'_{i+1}[-{\frac {1}{8}}]$

(3)

$\displaystyle Z(s={\frac {1}{2}})={\frac {1}{2}}[Z_{i}+Z_{i+1}]+{\frac {1}{8}}[Z'_{i}-Z'_{i+1}]$

(4)

But we know from eqn 1P. 36-1 that $Z'=h{\dot {Z}}$. By imposing collocation we have $Z_{i}=f_{i}$ and $Z_{i+1}=f_{i+1}$. Using these results in the above equation we get

$\displaystyle Z_{i+1/2}=Z(s=1/2)={\frac {1}{2}}[Z_{i}+Z_{i+1}]+{\frac {h}{8}}[f_{i}-f_{i+1}]$

Author

--Egm6341.s10.team2.niki 23:03, 5 April 2010 (UTC)

--Proofread by Egm6341.s10.team2.lee 00:37, 6 April 2010 (UTC)

Problem 8:Enforcing compliance with the DE at (s=1/2):Part B

Statement

P.36-2 By enforcing compliance with the DE at the point $Z_{i+1}=Z(s=1/2)$ show that

$\displaystyle Z'_{i+1/2}=Z'(s=1/2)=-{\frac {3}{2}}[Z_{i}-Z_{i+1}]-{\frac {1}{2}}[Z'_{i}+Z'_{i+1}]$

Solution

We have from eqn 1 P.35-4 and Prob 4 of HW6 [2] the following

$\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}}$ (1 p 35-4)

We know from eqn 3 P.35-3

$\displaystyle Z'(s)=\sum _{i=1}^{3}ic_{i}s^{i-1}$

Expanding and substituting for the constants we get:

$\displaystyle Z'(s=1/2)=c_{1}+2c_{2}({\frac {1}{2}})+3c_{3}({\frac {1}{4}})$

(1)

$\displaystyle Z'(s=1/2)=Z_{i}+2[-3Z_{i}-2Z'_{i}+3Z_{i+1}-Z'_{i+1}]({\frac {1}{2}})+3[2Z_{i}+Z'_{i}-2Z_{i+1}+Z'_{1+1}]({\frac {1}{4}})$

(2)

Rearranging terms we get

$\displaystyle Z'(s=1/2)=Z_{i}[-3+{\frac {3}{2}}]+Z'_{i}[1-2+{\frac {3}{4}}]+Z_{i+1}[3-{\frac {3}{2}}]+Z'_{i+1}[-1+{\frac {3}{4}}]$

(3)

$\displaystyle Z'_{i+1/2}=Z'(s=1/2)=-{\frac {3}{2}}[Z_{i}-Z_{i+1}]-{\frac {1}{2}}[Z'_{i}+Z'_{i+1}]$

which is the required result.

Author

--Egm6341.s10.team2.niki 23:29, 5 April 2010 (UTC)

--Proofread by Egm6341.s10.team2.lee 00:37, 6 April 2010 (UTC)

Problem 9: Condition for collocation at $t_{i+1/2}$

Statement

P. 36-3 To prove that the condition for collocation at $\;t_{i+1/2}$ is:

$Z_{i+1}=Z_{i}+{\frac {\frac {h}{2}}{3}}\left[f_{i}+4f_{i+{\frac {1}{2}}}+f_{i+1}\right]$

Solution

From P.36-2 we know that:

$Z'_{i+{\frac {1}{2}}}={\frac {-3}{2}}\left(Z_{i}-Z_{i+1}\right)-{\frac {1}{4}}\left(Z'_{i}+Z'_{i+1}\right)$

using ${\dot {Z}}={\frac {1}{h}}Z'\;$, we get:

${\dot {Z}}_{i+{\frac {1}{2}}}={\frac {-3}{2h}}\left(Z_{i}-Z_{i+1}\right)-{\frac {1}{4}}\left({\dot {Z}}_{i}+{\dot {Z}}_{i+1}\right)$

${\dot {Z}}_{i+{\frac {1}{2}}}={\frac {-3}{2h}}\left(Z_{i}-Z_{i+1}\right)-{\frac {1}{4}}\left(f_{i}+f_{i+1}\right)$

Subtracting $\;f_{i+{\frac {1}{2}}}$ from both sides:

${\dot {Z}}_{i+{\frac {1}{2}}}-f_{i+{\frac {1}{2}}}={\frac {-3}{2h}}\left(Z_{i}-Z_{i+1}\right)-{\frac {1}{4}}\left(f_{i}+f_{i+1}\right)-f_{i+{\frac {1}{2}}}$

We know that in order for the collocation to be satisfied at point $t_{i+1/2}$, $\;{\dot {Z}}_{i+{\frac {1}{2}}}-f_{i+{\frac {1}{2}}}=0$

$\therefore {\frac {-3}{2h}}\left(Z_{i}-Z_{i+1}\right)-{\frac {1}{4}}\left(f_{i}+f_{i+1}+4f_{i+{\frac {1}{2}}}\right)=0$

$\left(Z_{i+1}-Z_{i}\right)={\frac {2h}{3*4}}\left[f_{i}+f_{i+1}+4f_{i+{\frac {1}{2}}}\right]$

$\therefore \;Z_{i+1}=Z_{i}+{\frac {\frac {h}{2}}{3}}\left[f_{i}+4f_{i+{\frac {1}{2}}}+f_{i+1}\right]$

Author

Solved by-Srikanth Madala (SM)

---

Problem 10: Commentary on the matlab code available on Kessler's paper for the evaluation of polynomials and their coefficients in the trapezoidal rule error

Statement

P. 37-1

Run Kessler's code to reproduce his table, and complete the HW of P.30-2, i.e., follow Kessler's code line by line to produce the pairs $\;[P_{2}(x),P_{3}(x)],\;[P_{4}(x),P_{5}(x)],\;[P_{6}(x),P_{7}(x)]$ and write comments on his code.

Solution

Matlab code

function [c,p]= traperror(n)
%traperror is a user-defined function which takes in 'n' as the input argument
%and returns: 1.The values of the coefficients c_1, c_3, ...in the polynomial p.
%2. The values of the polynomials p_2, p_4, ... , p_{2n} at t=+1

f=1; g=2;
cn=-1; cd=1;
%defining variables 'f,g,cn and cd'
%cn is the numerator of the coefficients c_1, c_3, ...
%cd is the denominator of the coefficients c_1, c_3, ...
%cn is initially set to 1 and cd to -1 because c_1 = -1,  i.e,(-1/1).

for k=1:n
%beginning of 'for' repeatition loop with 'n' as the input value of 'traperror'
%function signifying the desired number of polynomials p_2, p_4, ... , p_{2n} at t=+1.

f=f.*g.*(g+1);
%This is an iterative step for 'f' to generates the odd-numbered factorials like 1! =1 , 3! =6, 5!=120 and so on.
%In each iteration, the value of 'f' is reassigned by using element-wise multiplication of matrices f,g and (g+1)

[newcn,newcd]=fracsum(-1*cn,cd.*f);
%calling the user-defined function 'fracsum' in the process of computing
%the numerator and denominator of the coefficients c_1, c_3, ...that are used in the polynomials p_2, p_4, ... , p_{2n} in each iteration
%of the for loop. In each for loop iteration, we call the function
%'fracsum' and pass the arguments n=-1*cn and d=cd.*f to the function fracsum(n,d).
%The function fracsum returns the values [nsum,dsum] which are stored in
%[newcn,newcd]. It should be noted that the values of cn,cd and f all
%incremented with an additional row element in each successive iteration.

cn=[cn;newcn]; cd=[cd;newcd];
% the 'cn' (numerator of the coefficient)  and 'cd' (denominator of the coefficient) values are reassigned in
%the execution of the For loop. In each iteration of the for loop, the new values of 'cn' and 'cd' matrices
%are incremented by one additional row element 'newcn' and 'newcd' respectively

f=[f;1];
%the 'f' matrix is incremented by an additional row element=1, in each iteration of the for loop

g=[2+g(1);g];
% the 'g' matrix is incremented by an additional row element=even number, in each iteration of the for loop

[newpn,newpd]=fracsum((g-1).*cn,f.*cd);
%calling the user-defined function 'fracsum' in the process of evaluating the numerator and denominator of the polynomials p_2, p_4, ... , p_{2n}
%computed at t=+1. calling the function 'fracsum' and passing the arguments n=(g-1).*cn and d=f.*cd to the function fracsum(n,d).
%The function fracsum returns the values [nsum,dsum] which are stored in [newpn,newpd]

pn(k,1)=newpn;
%In each 'k'th ('k' varying from 1 to n) iteration of the for loop, numerator element of each polynomial
%function-P_2k(t=1) is stored in the 'k'th row of a column matrix 'pn'

pd(k,1)=newpd;
%In each 'k'th ('k' varying from 1 to n) iteration of the for loop, denominator element of each polynomial
%function-P_2k(t=1) is stored in the 'k'th row of a column matrix 'pd'

end
% end statement of the for loop to check and see if k=n; the condition k=n ends the loop
% or else the loop is repeated until the statement is true.

c=int64([cn cd]);
% converts the augmented matrix 'c' containing the numerator and the denominator of the coefficients- c_1, c_3, ... to signed 64-bit integers

p=int64([pn pd]);
% converts the augmented matrix 'p' containing the numerator and the denominator of the polynomial functions- p_2, p_4, ... , p_{2n}
% evaluated at t=+1 to signed 64-bit integers

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [nsum,dsum]= fracsum(n,d)
%user-defined function which takes in 'n,d' as the input arguments and returns: 'nsum' and 'dsum'.
%This function is called twice in the execution of the function traperror(n). First time, the arguments 'n=-1*cn' and 'd=cd.*f'
%are passed to this function and the second time, the arguments 'n=(g-1).*cn'and 'd=f.*cd'
%are passed to this function. During the first time call of each iteration, this function
%successively returns the the numerator and denominator of the coefficients
%c_1, c_3, ...that are used in the polynomials p_2, p_4, ... , p_{2n}. During the second time call of each iteration, this function
%successively returns the the numerator and denominator of the the polynomials p_2, p_4, ... , p_{2n}computed at t=+1
%[c1/(2k-1)!] + [c3/(2k-3)!] + [c5/(2k-5)!] +....+ [c2k-5/5!] + [c2k-3/3!]+ [c2k-1/1!] = 0 is the recurssive formula used in the computation of
%c3,c5,...c2k-5, c2k-3 and c2k-1 starting with c1=-1. 'f' in the argument
%passed is essentially the odd number factorial in the above formula.
%Therefore 'n' and 'd' that are passed to the function are numerator and
%denominator of each term that is calculated as the iteration progresses.

div=gcd(round(n),round(d));
%'round' is a matlab function to rounding off the fraction to the nearest integer and 'gcd' is a matlab function that
%returns the greatest common divisor of round(n) and round(d). This value of gcd is stored in a variable 'div'.
%In the first time call of the function in each iteration, 'div' matrix
%elements are assigned the elemental gcd's of the coefficient numerator-  (cn)*-1 and
%coefficient denominator*odd number factorial- cd.*f

n=round(n./div);
% 'n' matrix is reassigned by carrying out element wise division with 'div' matrix and then rounding off

d=round(d./div);
% 'd' matrix is reassigned by carrying out element wise division with 'div' matrix and then rounding off

dsum=1;
% a new variable 'dsum' is initially assigned a value of 1

for k=1:length(d)
% a new 'for' loop is started; the 'length' function measures the row size of
% the column matrix 'd'

dsum=lcm(dsum,d(k));
%'lcm' is a matlab function that returns the lowest common multiple of the values in the parentheses and in this case dsum and d(k)
% The lcm of the denominator is taken successively to calculate the new denominator sum and the value is stored in the function
% return variable-'dsum'

end
% end of for loop to check and see if k=length(d)?

nsum=dsum*sum(n./d);
%'nsum' variable is defined by carrying out the arithemetic operation between 'dsum' and sum(n./d)
% calculation of the sum of successive coefficients- c3,c5,...c2k-5, c2k-3
% and c2k-1 using sum function which is further multiplied by 'dsum' to yield an pseudo 'nsum'

div=gcd(round(nsum),round(dsum));
%'div' is reassigned a different value which is the GCD of rounded values of 'nsum' and 'dsum'

nsum=nsum/div;
%'nsum' is redefined
% the actual new numerator sum is found and stored in the function return variable-'nsum'

Author

Solved by-Srikanth Madala (SM)

---

Problem 11: Verify the results of Problem 10, HW #5 by using Team 4's code to debug

Statement

P. 37-1

An Ellipse was given:

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

where:

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

$a=1$

Calculate the Circumference of the elipse as follows:

Method 1:

$C=\int _{0}^{2\pi }\left[r^{2}+({\frac {dr}{d\theta }})^{2}\right]^{\frac {1}{2}}$

Method 2:

$C=4aE(e^{2})$

$E(e^{2})=\int _{0}^{\frac {\pi }{2}}\left[1-e^{2}\sin(\theta )\right]^{\frac {1}{2}}d\theta$

Use Team 4's Solution to debug own code and compare results.

Solution

Composite Trapezoidal Rule

After looking at Team 4's codes I was able to gain knowledge into checking for the error of the numerical integration. Like our team, they also used the 'Quad' command as the theoretical answer to compare their results to. In addition they used an extra argument to the 'Quad' command which provides for more accurate results by providing the order of the error.

In addition the way the derivative of the radius relative to the angle is calculated differently. We used Matlab's Symbolic toolbox to find the integration results while team 4 used the theoretical results then plugged in values for their respective variables.

Regarding Method 2, the same applies.

Below are our results

Method 1

Composite Trapezoidal Rule
Team	n	Value	Error	Team Difference
Team 2	18	6.176517900343786	6.568967592102126e-012	0.000084086656214
Team 4	18	6.176601987	2.89E-11

Method 2

Composite Trapezoidal Rule
Team	n	Value	Error	Team Difference
Team 2	4	6.176601987658692	5.062616992290714e-014	0.000000000658692
Team 4	4	6.176601987	6.81E-12

Matlab Code Used

Modified Matlab code for composite Trapezoidal calculation

function [I,E]=circompos(a,b,n)
format long
h=(b-a)/n;
pi2=2*pi();
i=0;
Itot=0;
It=0;
It2=0;
 %replace elarc with function calculator
while i<=n
    if i==0
        Itot1=(1/2)*elarc(a);
    else if i<n
            u=h*i;
            It(i)= elarc(u);
        else
            It2=(1/2)*elarc(b);
        end
    end
Itot=Itot1+sum(It)+It2;
i=i+1;
end

I=Itot*(h);

%Error Calculation, replace elarc with function calculator
Ir=quad(@elarc,0,pi2,10^-11) %added the extra argument for the error
E=Ir-I;

Clencurt Command

Comparing the two codes used for using the 'clencurt' command were different. Our approach took a simple way to just use the weights to calculate our results with no error bounding, as it was done by Team 4.

Method 1

Composite Trapezoidal Rule
Team	n	Value	Error	Team Difference
Team 2	23	6.176517900343786	7.505640553517878e-011	0.000084086656214
Team 4	23	6.176601987	7.55E-11

Method 2

Composite Trapezoidal Rule
Team	n	Value	Error	Team Difference
Team 2	9	6.176517900343786	9.930172963734663e-005
Team 4	9	6.176601987	1.59E-11

Modified Matlab Code For Clencurt

Matlab code

function [I,E] = runclencurt(n)
[xx,ww] = clencurt(n);                               %Calls the clencurt function
    xx=pi+xx.*pi;
    ww=ww.*pi;
    F=elarc(xx);  %Evaluation of the function
    I_cl=ww*F;
    I=6.176517900343786;
    E=abs(I_cl-I);

Romberg Table

The solution approaches to creating the Romberg Table were different from both teams. Team 4's approach uses the simple trapezoidal rule formula as compared to Our team's approach of using the matlab 'trapz' command. In addition the method used to calculate row by row of the Romberg Table were different.

When running Team 4's code I'm not able to get all of the Romberg Table to be completed. At the moment Team 2 is continiously working on figuring out the cause for the difference of the solutions as shown below:

Method 1

Composite Trapezoidal Rule
Team	n	Value	Error	Team Difference
Team 2	128	6.176517900343941	0.486692881551676	0.000084086656059
Team 4	128	6.176601987	1.10E-11

Method 2

Composite Trapezoidal Rule
Team	n	Value	Error	Team Difference
Team 2	32	6.176601987658692	-0.295956574802858	0.000084086656214
Team 4	32	6.176601987	5.30E-14

Author

Solved by--Guillermo Varela (GV)

---

Signatures

Solved problems 3,7,8 and proofread --Niki Nachappa (NN) 23:31, 5 April 2010 (UTC)

Solved problems 4,5,9,10 --Srikanth Madala (SM)

Solved problems 2, 11 and Proofread 4, 8 --Guillermo Varela (GV) 20:50, 7 April 2010 (UTC)

Solved problems 1,6 and Proofread 3,7,8 --Egm6341.s10.team2.lee 02:12, 8 April 2010 (UTC)

Proofread problems 4,9,5,10 --Patrick O'Donoughue(PO)

Solution Assignments and Reviewers

Problem Assignments
Problem #	Solution	Reviewed
Problem 1	JP	NN,GV
Problem 2	GV	PO,JP
Problem 3	NN	SM,PO
Problem 4	PO	GV,SM
Problem 5	SM	JP,GV
Problem 6	JP	PO,GV
Problem 7	NN	SM,JP
Problem 8	NN	GV,PO
Problem 9	PO	JP,SM
Problem 10	SM	NN,GV
Problem 11	GV	PO,JP