User:Egm6341.s11.team3/sub6

From Wikiversity
Jump to navigation Jump to search

= 6.1 Proof of HOTRE theorem additional terms derivation Part b= [1] [2]

Given: p4 and p5 values along the error expression and residual integral term[edit | edit source]

From Proof of HOTRE notes by Prof. Vu-Quoc, functions expressions are given as:

with expression for at term is given as

Find: Expression for p6, p7, E at p6 and p7, the expression for global-to-local variable conversion and values for functions 'd'[edit | edit source]

(a) Expression for and

(b) Find

(c) Find

Solution[edit | edit source]

(a)[edit | edit source]

Having this expression for error

, we integrate by parts the residual part once again

by noting that

then by taking integral over the interval, for this more specific case, in interval , we get

Then re-writing the equations above

2 things must be noted:

  1. the rewritten term is the same as the residual term in error as well as the residual next term is of the same format as term.
  2. note that by the definition of integration by parts, that was done in deriving expression for , we expressed or

Therefore,

, where

in same fashion we expand it once more

, where we have now

in order to reduce some error in residual function and noting that this is an error that generated from integration, we want to make 1 function odd such that negative and positive contribution to error would cancel out. Since we note that are odd functions, we will force with constants for this function to cancel out at the interval

That is:

They both end at the same point, namely '0'

as well as the function will pass though origin, in order to enforce the oddity of the function (remember that odd functions are symmetric with origin and, therefore, they pass though origin)

then by noting that constants and were already found for lower order 'p' functions and

we plug in the final boundary value to find

Therefore

(b)[edit | edit source]

From the equation (3) given in Prof. Vu-Quoc's notes

, where is the size of an interval, namely,

By re-arranging members

(c)[edit | edit source]

From Prof. Vu-Quoc's notes on expression for function "d"

Miscellaneous/Matlab Code for caluclations[edit | edit source]

--Egm6341.s11.team3.JA 03:07, 4 April 2011 (UTC)

6.2 Derivation of Trapezoidal Rule Error [3][edit | edit source]

 REMARK: We did the solution on our own 

Given: The Error Expression of the Trapezoidal Rule[edit | edit source]

where:

Find: Derive the Given Expression[edit | edit source]

Derive equations 6.2.1 and 6.2.2.

Solution[edit | edit source]

From class lecture, we have the following error expression:

[4]

Where:

[5]

By recognizing that x(+1)=b and x(-1)=a, equation 6.2.4 may be used to solve for three of the parameters in equation 6.2.3 as follows:

Substituting these three expressions (Equations 6.2.5, 6.2.6, and 6.2.7) into the error expression in Equation 6.2.3 yields:

By substituting the expression given by Equation 6.2.2 into Equation 6.2.8, the final answer is achieved:

However, this answer is different than the one given in the problem statement. Unless the solution given in the problem statement is incorrect, I believe this is likely due to a misprint in the lecture notes. Had the expression in Equation 6.2.3 been given instead as:

Then, by following the same exact steps shown above, the derivation would instead have led to:

Which, after making the substitution given by Equation 6.2.2, agrees with the desired final solution. So, by assuming the aforementioned misprint in Lecture 31 Page 3 Equation 2, the error may be expressed as:

Egm6341.s11.team3.russo 19:35, 28 March 2011 (UTC)

References[edit | edit source]

  1. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg31.djvu&page=1
  2. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg31.djvu&page=2
  3. (Lecture 31 page 3)
  4. (meeting 31 page 3).
  5. (meeting 30 page 3).

6.3 Bernoulli Number[1][edit | edit source]

Given[edit | edit source]










Find[edit | edit source]


Solution[edit | edit source]

by Eq(6.3.1)









Therefore



by Eq(6.3.2)
Thus







Therefore




Egm6341.s11.team3.sungsikkim 03:22, 29 March 2011 (UTC)

6.4 To derive the Higher Order Trapezoidal Rule Error equation by canceling odd order derivative terms[edit | edit source]

Given[edit | edit source]

Refer Lecture slide http://en.wikiversity.org/w/index.php?title=File%3ANm1.s11.mtg30.djvu&page=2 for problem statement.

.

where can be written as,

.

.

where .

After successive integration by parts,

Find[edit | edit source]

Try to derive the Higher Order Trapezoidal Rule Error equation by canceling odd order derivative terms of , i.e., make terms to zero at .

Solution[edit | edit source]

We know,

At ,

We also know,

And,

At ,

.

Functions are odd functions.

Since are odd functions, we know .

And henceforth becomes

Using the below relationship

.

In Equation (4) can be transferred to as shown below:

.

Egm6341.s11.team3.rakesh 05:16, 2 April 2011 (UTC)

6.5 Determining polynomial pairs using recurrence relationship[edit | edit source]

Given[edit | edit source]

Recurrence Relationship

Where;

Find[edit | edit source]

Calculate the coefficients and the pairs of polynomials using the Recurrence formula

Solution[edit | edit source]

Using the recurrence relation we can determine the constants to evaluate the pairs of polynomials,

with i varying from 1 to 4 ,we have 4 equations and 4 unknowns which can be solved simultaneously to find the constants.

Given that,

when

when

when

substitute the values and solve for c_{9} ,

Now we calculate the pair polynomials using the expressions given below for odd and even terms , and then the constants are substituted in them ,

The expressions are generated as follows,

Now we get the final Pair Polynomials ,by just substituting the constants in the maple generated equations and the pairs are listed below :

6.6 Understand Kessler's code[edit | edit source]

Given: Kessler's code[edit | edit source]

Kessler's code is given below

Find: understand Kessler's code line by line[edit | edit source]

Using to understand Kessler's code line by line, starting with the best of Spring 2010

Solution[edit | edit source]

The Kessler's code is to compute the coefficients and the constants where , according to the following algorithm[2]:

1) :

2) :

  • function traperror
function [c,p]=traperror(n)
The input of this function  is a integer that ranges from 1 to . The function returns two arrays  and , corresponding to the constants  and the coefficients .
f=1; g=2;
cn=-1; cd=1;
Set initial values for f, g, cn and cd, where the initial value of cn and cd are the numerator and denominator of , respectively.
for k=1:n
Start iteration from  to 
   f=f.*g.*(g+1);
   % k=1  g=[4;2]
   % k=2  g=[6;4;2]
   % k=3  g=[8;6;4;2]
  Update vector f according to , where  is the denominator of the factors of  in (Eqn 6.6.2) as 
   [newcn,newcd]=fracsum(-1*cn,cd.*f);
   cn=[cn;newcn]; cd=[cd;newcd];
   % k=1 cn=[-1;1], cd=[1;6]
   % k=2 cn=[-1;1;-7], cd=[1;6;360]
   % k=3 cn=[-1;1;-7;31], cd=[1;6;360;15120
  Calculate  according to (Eqn 6.6.2)
   f=[f;1];
   g=[2+g(1);g];
   % k=1  f=[6;1]  g=[4;2]
   % k=2  f=[120;6;1]  g=[6;4;2]
   % k=3  f=[5040;120;6;1]  g=[8;6;4;2]
  Update vector f and g, where  and  are respectively the numerator and denominator of the factors of  in (Eqn 6.6.2) as 
   [newpn,newpd]=fracsum((g-1).*cn,f.*cd);
   % k=1  newpn=-1, newpd=6
   % k=2  newpn=1, newpd=360
   % k=3  newpn=-2, newpd=15120
  Calculate  according to (Eqn 6.6.3)
   pn(k,1)=newpn; pd(k,1)=newpd;
   % k=1  pn(1,1)=-1, pd(1,1)=6
   % k=2  pn(2,1)=1, pd(2,1)=360
   % k=3  pn(3,1)=-2, pd(3,1)=15120
end
c=int64([cn cd]);
p=int64([pn pd]);
% k=1 c=[-1 1; 1 6], p=[-1 6]
% k=2 c=[-1 1; 1 6; -7 360],p=[-1 6; 1 360]
% k=3 c=[-1 1; 1 6; -7 360; 31 15120], p=[-1 6; 1 360; -2 15120]
Add  and  to the output arrays c and p
  • function fracsum
function [nsum,dsum]=fracsum(n,d)
This function is to calculate the summation of  from vectors  and . The outputs nsum and dsum returns the numerator and denominator,respectively. 
div=gcd(round(n),round(d));
The function round is used to round off a number to its nearest integer. The function gcd returns an array containing the greatest common divisors of the corresponding elements of the two input arrays.
n=round(n./div);
d=round(d./div);
This is to remove the common divisors from  and , while the ratio  is preserved.
dsum=1;
for k=1:length(d)
dsum=lcm(dsum,d(k));
end
The function lcm returns the least common multiple of corresponding elements of two input arrays
nsum=dsum*sum(n./d);
div=gcd(round(nsum),round(dsum));
nsum=nsum/div;

Egm6341.s11.team3.ren 06:31, 30 March 2011 (UTC)

References[edit | edit source]

  1. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg31.djvu&page=3
  2. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg32.djvu&page=3

6.7 Compute Circumferences [1] [2] [3][edit | edit source]

 REMARK: We did the solution on our own 

Given: Bifolium and Ellipse Expressions[edit | edit source]

The expression of a bifolium:

with

And an ellipse with :

Find: The Circumferences[edit | edit source]

Compute the circumferences of the bifolium and the ellipse by integrations:

with the ellipse represented by:

Also compute the circumferences of the ellipse using:

Computer the integrations using the following methods:

a) Comp. Trap.

b) Romberg.

c) Clenshaw-Curtis.

d) Fast Clenshaw-Curtis.

Solution:[edit | edit source]

(1) Circumference for Bifolium:[edit | edit source]

From (Eq. 6.7.1), (Eq. 6.7.2), we derive the Integration as:

Using WolframAlpha, the exact solution is calculated to be [4]:

a)

The Comp. Trap. rule [5] is written as:

b)

The Romberg method is implemented following the Richardson extrapolation [6] :

c)

The Clenshaw-Curtis quadrature method could be realized with the Trefethen's clencurt code [7]. The original integration could be written as:

where , and determined by,

d)

The fast Clenshaw-Curtis method involves fast fft, and the code is obtained at [8].

The matlab code for computing the circumference of the bifolium is attached below:

Results:[edit | edit source]

The integration results for a), c) and d) are shown and compared in the following table:

a) Comp. Trap. c) Clenshaw-Curtis d) Fast C-C
2 3.14159 0.43618 2.09440 1.48337 3.14159 0.43618
4 3.14159 0.43618 2.13113 1.44664 3.49066 0.08711
8 3.49462 0.08315 3.34207 0.23570 3.67745 0.09968
16 3.55150 0.02627 3.51396 0.06381 3.61754 0.03977
32 3.57134 0.00643 3.56185 0.01592 3.58628 0.00851
64 3.57617 0.00160 3.57381 0.00396 3.57982 0.00205
128 3.57737 0.00040 3.57678 0.00099 3.57828 0.00051
256 3.57767 0.00010 3.57753 0.00024 3.57790 0.00013
512 3.57775 0.00002 3.57771 0.00006 3.57780 0.00003
1024 3.57777 0.00000 3.57776 0.00001 3.57778 0.00001
2048 3.57777 0.00000 3.57778 0.00001
Elapsed time(s) 0.023071 0.443672 0.024568


The Romberg table for b) is:

2 3.14159 0.43618 3.14159 0.43618 3.64368 0.06591 3.56646 0.01131 3.57867 0.00090 3.57775 0.00002 3.57777 0.00000
4 3.14159 0.43618 3.61230 0.03453 3.56767 0.01010 3.57862 0.00085 3.57775 0.00002 3.57777 0.00000
8 3.49462 0.08315 3.57046 0.00731 3.57845 0.00068 3.57775 0.00002 3.57777 0.00000
16 3.55150 0.02627 3.57795 0.00018 3.57776 0.00001 3.57777 0.00000
32 3.57134 0.00643 3.57778 0.00001 3.57777 0.00000
64 3.57617 0.00160 3.57777 0.00000
128 3.57737 0.00040


(2) Circumference for Ellipse: Method1[edit | edit source]

From definition of eccentricity of ellipse [9]:

From (Eq. 6.7.2), (Eq. 6.7.3), we derive the Integration as:

Using WolframAlpha, the exact solution is calculated to be [10]:

Following the same approaches as in (Eq. 6.7.7), (Eq. 6.7.8) and (Eq. 6.7.9), the integration is computed using the following code:

Results:[edit | edit source]

The integration results for a), c) and d) are shown and compared in the following table:

a) Comp. Trap. c) Clenshaw-Curtis d) Fast C-C
2 62.8319 18.9728 42.8526 1.0065 7.8152 36.0439
4 35.3235 8.5356 24.0379 19.8212 20.5879 23.2712
8 28.5676 15.2915 48.6865 4.8274 48.3473 4.4882
16 36.4533 7.4058 43.7775 0.0816 43.6728 0.1863
32 42.8790 0.9801 43.8597 0.0006 43.8598 0.0007
64 43.8038 0.0553 43.8591 0.0000 43.8591 0.0000
128 43.8582 0.0009
256 43.8591 0.0000
Elapsed time(s) 0.016347 0.022058 0.017077


The Romberg table for b) is:

2 62.8319 18.9728 26.1541 17.7050 26.3264 17.5327 40.1489 3.7102 45.5249 1.6658 44.0225 0.1634 43.8568 0.0023 43.8582 0.0009
4 35.3235 8.5356 26.3156 17.5435 39.9329 3.9262 45.5039 1.6448 44.0240 0.1649 43.8568 0.0023 43.8582 0.0009
8 28.5676 15.2915 39.0818 4.7773 45.4169 1.5578 44.0298 0.1707 43.8570 0.0021 43.8582 0.0009
16 36.4533 7.4058 45.0210 1.1619 44.0514 0.1923 43.8577 0.0014 43.8582 0.0009
32 42.8790 0.9801 44.1120 0.2529 43.8607 0.0016 43.8582 0.0009
64 43.8038 0.0553 43.8764 0.0173 43.8582 0.0009
128 43.8582 0.0009 43.8594 0.0003
256 43.8591 0.0000


(3) Circumference for Ellipse: Method2[edit | edit source]

The 2nd method of circumference for ellipse follows (Eq. 6.7.4). Similar to previous calculations, the function integrated here is defined as:

Using WolframAlpha, the exact solution is calculated to be [11]:

Which is identical to (Eq. 6.7.13), following the same approaches as in (Eq. 6.7.7), (Eq. 6.7.8) and (Eq. 6.7.9), the integration is computed using the following code:

Results:[edit | edit source]

The integration results for a), c) and d) are shown and compared in the following table:

a) Comp. Trap. c) Clenshaw-Curtis d) Fast C-C
2 43.6129 0.2462 44.5369 0.6778 44.0921 0.2330
4 43.8525 0.0066 43.8395 0.0196 43.7309 0.1282
8 43.8591 0.0000 43.8591 0.0000 43.8595 0.0004
16 43.8591 0.0000
Elapsed time(s) 0.014058 0.019758 0.015220


The Romberg table for b) is:

2 43.6129 0.2462 43.9324 0.0733 43.8565 0.0026
4 43.8525 0.0066 43.8613 0.0022
8 43.8591 0.0000


Egm6341.s11.team3.Xia 18:00, 31 March 2011 (UTC)

References[edit | edit source]

  1. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg33.djvu&page=1
  2. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg33.djvu&page=2
  3. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg33.djvu&page=3
  4. http://www.wolframalpha.com/input/?i=integration+%28%282*sinx*cosx%5E2%29%5E2%2B%282*cosx%5E3-4*sinx%5E2*cosx%29%5E2%29%5E%280.5%29+for+x+from+0+to+pi
  5. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg7.djvu&page=4
  6. http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg29.djvu&page=5
  7. http://people.maths.ox.ac.uk/trefethen/clencurt.m
  8. http://www.mathworks.com/matlabcentral/fileexchange/6911-fast-clenshaw-curtis-quadrature
  9. http://mathworld.wolfram.com/Ellipse.html
  10. http://www.wolframalpha.com/input/?i=integrate+%289%2F10%29%281-%28sqrt%2891%29%2F10%29cosx%29%5E%28-2%29sqrt%28191%2F100-%28sqrt%2891%29%2F5%29cosx%29+from+0+to+2pi
  11. http://www.wolframalpha.com/input/?i=integrate+40sqrt%281-0.91sin%5E2x%29+for+x+from+0+to+pi%2F2

6.8 Derivation of Arc Length Expression [1][edit | edit source]

 REMARK: We did the solution on our own 

Given: A picture of Triangle AOB and the expression for Arc Length PQ [2][edit | edit source]

TriangleAOB

Find: Derive the Given Expression[edit | edit source]

Use Triangle AOB and the Law of Cosines to derive Equation 6.8.1

Solution[edit | edit source]

The Law of Cosines [3] is as follows:

Where a, b and c are the sides of a triangle and is the angle between sides a and b.

Upon inspection it is clear that the Law of Cosines may be applied to Triangle OAB (pictured above):

In order to achieve the desired results, we must denote the triangular components in terms of elliptical components:

Where the desired arc length PQ is the infinitesimal arc length dl integrated over the angle POQ. Thus, we must find an expression for dl. Rewriting Equation 6.8.3 in terms of elliptical components yields:

Solving for our desired parameter dl yields:

Thus an expression governing the infinitesimal arc length dl has been developed. As stated previously, the desired arc length PQ is equivalent to this infinitesimal arc length dl integrated over the angle POQ. By denoting the angles that points Q and P make with the origin and axis as and , respectively, the arc length may be found:

Which is exactly equivalent to the expression given in the problem statement by Equation 6.8.1.

Egm6341.s11.team3.russo 22:38, 28 March 2011 (UTC)

References[edit | edit source]

  1. (Lecture 33 page 2)
  2. (Lecture 32 page 4)
  3. (Wikipedia Article on Law of Cosines)

6.9 Determine 3rd and 4th rows of the Matrix A, verify the inverse of the Matrix,Identify the basis functions , . Plot 's=[edit | edit source]

[1],[2]

Given[edit | edit source]

6.9.a) The matrix is given to be:

.

Find[edit | edit source]

.

and

.

Solution[edit | edit source]

We know that

.

and

.

thus using the above relations where in both cases we then have

.

and

.

Hence we can put forth that the 3rd and 4th rows are and respectively.

Given[edit | edit source]

6.9.b)

From the matrix

.

Find[edit | edit source]

Verify the inverse of the above given matrix

.

Solution[edit | edit source]

A=[1 0 0 0;0 1 0 0;1 1 1 1;0 1 2 3];
Ainv=inv(A);
Ainv

The above sequence of commands generates a matrix of the elements in it as shown below

Ainv =

     1     0     0     0
     0     1     0     0
    -3    -2     3    -1
     2     1    -2     1

6.9.c) Identify the basis functions , . Plot 's[edit | edit source]

Solution[edit | edit source]







































Note for Verifying that :









Hv
Plot of N1, N2, N3, N4 as a function s

Egm6341.s11.team3.rakesh 19:21, 4 April 2011 (UTC)

Contributing Members[edit | edit source]

Team Contribution Table
Problem Number Assigned To Solved By Typed By Proofread
6.1 Janulevicius, Arunas Janulevicius, Arunas Janulevicius, Arunas Russo, Jr Charles T
6.2 Russo, Jr Charles T Russo, Jr Charles T Russo, Jr Charles T rakesh
6.3 Kim-1, Sungsik Kim-1, Sungsik Kim-1, Sungsik Balu, Elango
6.4 Yerramsetti, Surya Rakesh rakesh rakesh Russo, Jr Charles T
6.5 Balu, Elango Balu, Elango Balu, Elango Xia, Xi
6.6 Ren,Zheng Xia, Xi Xia, Xi Janulevicius, Arunas
6.7 Xia, Xi Yerramsetti, Surya Rakesh Yerramsetti, Surya Rakesh rakesh
6.8 Russo, Jr Charles T Russo, Jr Charles T Russo, Jr Charles T Kim-1, Sungsik
6.9 Yerramsetti, Surya Rakesh Yerramsetti, Surya Rakesh Yerramsetti, Surya Rakesh Xia, Xi