= 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
p
4
(
t
)
,
p
5
(
t
)
{\displaystyle p_{4}(t),p_{5}(t)\ }
expressions are given as:
p
4
(
t
)
=
c
1
t
4
4
!
+
c
3
t
2
2
!
+
c
5
{\displaystyle p_{4}(t)=c_{1}{\frac {t^{4}}{4!}}+c_{3}{\frac {t^{2}}{2!}}+c_{5}\ }
{\displaystyle \displaystyle }
p
5
(
t
)
=
c
1
t
5
5
!
+
c
3
t
5
5
!
+
c
5
t
{\displaystyle p_{5}(t)=c_{1}{\frac {t^{5}}{5!}}+c_{3}{\frac {t^{5}}{5!}}+c_{5}t\ }
{\displaystyle \displaystyle }
c
1
=
−
1
,
c
3
=
1
6
,
c
5
=
−
7
360
{\displaystyle c_{1}=-1,c_{3}={\frac {1}{6}},c_{5}=-{\frac {7}{360}}\ }
{\displaystyle \displaystyle }
with expression for
E
{\displaystyle E\ }
at
p
5
(
t
)
t
h
{\displaystyle p_{5}(t)^{th}\ }
term is given as
E
=
[
p
2
g
(
1
)
+
p
4
g
(
3
)
]
−
1
+
1
−
∫
−
1
+
1
p
5
g
(
5
)
d
t
{\displaystyle E=\left[p_{2}g^{(1)}+p_{4}g^{(3)}\right]_{-1}^{+1}-\ \int _{-1}^{+1}p_{5}g^{(5)}dt}
{\displaystyle \displaystyle }
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
p
6
,
p
7
{\displaystyle p_{6},p_{7}\ }
and
E
{\displaystyle E\ }
(b) Find
t
k
(
x
)
{\displaystyle t_{k}(x)\ }
(c) Find
d
1
=
d
¯
2.1
=
d
¯
2
=
p
2
(
1
)
2
2
{\displaystyle d_{1}={\bar {d}}_{2.1}={\bar {d}}_{2}={\frac {p_{2}(1)}{2^{2}}}\ }
{\displaystyle \displaystyle }
d
2
=
d
¯
4
{\displaystyle d_{2}={\bar {d}}_{4}\ }
{\displaystyle \displaystyle }
d
3
=
d
¯
6
{\displaystyle d_{3}={\bar {d}}_{6}\ }
{\displaystyle \displaystyle }
Having this expression for error
E
=
[
p
2
(
t
)
g
(
1
)
(
t
)
+
p
4
(
t
)
g
(
3
)
(
t
)
]
−
1
+
1
−
∫
−
1
+
1
p
5
(
t
)
g
(
5
)
(
t
)
d
t
{\displaystyle E=\left[p_{2}(t)g^{(1)}(t)+p_{4}(t)g^{(3)}(t)\right]_{-1}^{+1}-\ \int _{-1}^{+1}p_{5}(t)g^{(5)}(t)dt}
, we integrate by parts the residual part once again
by noting that
d
d
t
(
p
6
(
t
)
g
(
5
)
(
t
)
)
=
p
5
(
t
)
g
(
5
)
(
t
)
+
p
6
(
t
)
g
(
6
)
(
t
)
{\displaystyle {\frac {d}{dt}}\left(p_{6}(t)g^{(5)}(t)\right)=p_{5}(t)g^{(5)}(t)+p_{6}(t)g^{(6)}(t)\ }
{\displaystyle \displaystyle }
then by taking integral over the interval, for this more specific case, in interval
t
∈
[
−
1
,
+
1
]
{\displaystyle t\in [-1,+1]\ }
, we get
∫
−
1
+
1
d
d
t
(
p
6
(
t
)
g
(
5
)
(
t
)
)
=
∫
−
1
+
1
p
5
(
t
)
g
(
5
)
(
t
)
+
∫
−
1
+
1
p
6
(
t
)
g
(
6
)
(
t
)
{\displaystyle \int _{-1}^{+1}{\frac {d}{dt}}\left(p_{6}(t)g^{(5)}(t)\right)=\int _{-1}^{+1}p_{5}(t)g^{(5)}(t)+\int _{-1}^{+1}p_{6}(t)g^{(6)}(t)\ }
{\displaystyle \displaystyle }
[
p
6
(
t
)
g
(
5
)
(
t
)
d
t
]
−
1
+
1
=
∫
−
1
+
1
p
5
(
t
)
g
(
5
)
(
t
)
+
∫
−
1
+
1
p
6
(
t
)
g
(
6
)
(
t
)
{\displaystyle \left[p_{6}(t)g^{(5)}(t)dt\right]_{-1}^{+1}=\int _{-1}^{+1}p_{5}(t)g^{(5)}(t)+\int _{-1}^{+1}p_{6}(t)g^{(6)}(t)\ }
{\displaystyle \displaystyle }
Then re-writing the equations above
∫
−
1
+
1
p
5
(
t
)
g
(
5
)
(
t
)
=
[
p
6
(
t
)
g
(
5
)
(
t
)
d
t
]
−
1
+
1
−
∫
−
1
+
1
p
6
(
t
)
g
(
6
)
(
t
)
{\displaystyle \int _{-1}^{+1}p_{5}(t)g^{(5)}(t)=\left[p_{6}(t)g^{(5)}(t)dt\right]_{-1}^{+1}-\int _{-1}^{+1}p_{6}(t)g^{(6)}(t)\ }
{\displaystyle \displaystyle }
2 things must be noted:
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
∫
−
1
+
1
p
5
(
t
)
g
(
5
)
(
t
)
{\displaystyle \int _{-1}^{+1}p_{5}(t)g^{(5)}(t)\ }
term.
note that by the definition of integration by parts, that was done in deriving expression for
∫
−
1
+
1
p
5
(
t
)
g
(
5
)
(
t
)
{\displaystyle \int _{-1}^{+1}p_{5}(t)g^{(5)}(t)\ }
, we expressed
d
d
t
p
6
(
t
)
=
p
5
(
t
)
{\displaystyle {\frac {d}{dt}}p_{6}(t)=p_{5}(t)\ }
or
p
6
(
t
)
=
∫
p
5
(
t
)
{\displaystyle p_{6}(t)=\int p_{5}(t)\ }
Therefore,
E
=
[
p
2
(
t
)
g
(
1
)
(
t
)
+
p
4
(
t
)
g
(
3
)
(
t
)
]
−
1
+
1
+
[
p
6
(
t
)
g
(
5
)
(
t
)
d
t
]
−
1
+
1
−
∫
−
1
+
1
p
6
(
t
)
g
(
6
)
(
t
)
{\displaystyle E=\left[p_{2}(t)g^{(1)}(t)+p_{4}(t)g^{(3)}(t)\right]_{-1}^{+1}+\left[p_{6}(t)g^{(5)}(t)dt\right]_{-1}^{+1}-\int _{-1}^{+1}p_{6}(t)g^{(6)}(t)\ }
{\displaystyle \displaystyle }
, where
p
6
(
t
)
=
c
1
t
6
6
!
+
c
3
t
4
4
!
+
c
5
t
2
2
!
+
c
7
{\displaystyle p_{6}(t)=c_{1}{\frac {t^{6}}{6!}}+c_{3}{\frac {t^{4}}{4!}}+c_{5}{\frac {t^{2}}{2!}}+c_{7}\ }
in same fashion we expand it once more
E
=
[
p
2
(
t
)
g
(
1
)
(
t
)
+
p
4
(
t
)
g
(
3
)
(
t
)
]
−
1
+
1
+
[
p
6
(
t
)
g
(
5
)
(
t
)
d
t
+
p
7
(
t
)
g
(
6
)
(
t
)
d
t
]
−
1
+
1
−
∫
−
1
+
1
p
7
(
t
)
g
(
7
)
(
t
)
{\displaystyle E=\left[p_{2}(t)g^{(1)}(t)+p_{4}(t)g^{(3)}(t)\right]_{-1}^{+1}+\left[p_{6}(t)g^{(5)}(t)dt+p_{7}(t)g^{(6)}(t)dt\right]_{-1}^{+1}-\int _{-1}^{+1}p_{7}(t)g^{(7)}(t)\ }
{\displaystyle \displaystyle }
, where we have now
p
6
(
t
)
=
c
1
t
6
6
!
+
c
3
t
4
4
!
+
c
5
t
2
2
!
+
c
7
{\displaystyle p_{6}(t)=c_{1}{\frac {t^{6}}{6!}}+c_{3}{\frac {t^{4}}{4!}}+c_{5}{\frac {t^{2}}{2!}}+c_{7}\ }
{\displaystyle \displaystyle }
p
7
(
t
)
=
c
1
t
7
7
!
+
c
3
t
5
5
!
+
c
5
t
3
3
!
+
c
7
t
+
c
8
{\displaystyle p_{7}(t)=c_{1}{\frac {t^{7}}{7!}}+c_{3}{\frac {t^{5}}{5!}}+c_{5}{\frac {t^{3}}{3!}}+c_{7}t+c_{8}\ }
{\displaystyle \displaystyle }
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
p
1
,
p
3
,
p
5
,
.
.
.
.
{\displaystyle p_{1},p_{3},p_{5},....\ }
are odd functions, we will force with constants
c
1
−
c
7
{\displaystyle c_{1}-c_{7}\ }
for this function to cancel out at the interval
t
∈
[
−
1
,
+
1
]
{\displaystyle t\in [-1,+1]\ }
That is:
They both end at the same point, namely '0'
p
7
(
−
1
)
=
p
7
(
+
1
)
=
0
{\displaystyle p_{7}(-1)=p_{7}(+1)=0\ }
{\displaystyle \displaystyle }
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)
p
7
(
0
)
=
0
⇒
c
8
=
0
{\displaystyle p_{7}(0)=0\Rightarrow c_{8}=0\ }
{\displaystyle \displaystyle }
then by noting that constants
c
1
,
c
3
,
{\displaystyle c_{1},c_{3},\ }
and
c
5
{\displaystyle c_{5}\ }
were already found for lower order 'p' functions and
c
8
=
0
{\displaystyle c_{8}=0\ }
we plug in the final boundary value to find
c
7
{\displaystyle c_{7}\ }
p
7
(
1
)
=
(
−
1
)
1
7
7
!
+
1
6
1
5
5
!
+
(
−
7
360
)
1
3
3
!
+
c
7
(
1
)
+
0
=
0
{\displaystyle p_{7}(1)=(-1){\frac {1^{7}}{7!}}+{\frac {1}{6}}{\frac {1^{5}}{5!}}+\left(-{\frac {7}{360}}\right){\frac {1^{3}}{3!}}+c_{7}(1)+0=0\ }
{\displaystyle \displaystyle }
⇒
c
7
=
31
15
,
120
{\displaystyle \Rightarrow c_{7}={\frac {31}{15,120}}\ }
{\displaystyle \displaystyle }
Therefore
p
6
(
t
)
=
c
1
t
6
6
!
+
c
3
t
4
4
!
+
c
5
t
2
2
!
+
c
7
{\displaystyle p_{6}(t)=c_{1}{\frac {t^{6}}{6!}}+c_{3}{\frac {t^{4}}{4!}}+c_{5}{\frac {t^{2}}{2!}}+c_{7}\ }
p
7
(
t
)
=
c
1
t
7
7
!
+
c
3
t
5
5
!
+
c
5
t
3
3
!
+
c
7
t
{\displaystyle p_{7}(t)=c_{1}{\frac {t^{7}}{7!}}+c_{3}{\frac {t^{5}}{5!}}+c_{5}{\frac {t^{3}}{3!}}+c_{7}t\ }
E
=
[
p
2
(
t
)
g
(
1
)
(
t
)
+
p
4
(
t
)
g
(
3
)
(
t
)
+
p
6
(
t
)
g
(
5
)
(
t
)
d
t
]
−
1
+
1
−
∫
−
1
+
1
p
7
(
t
)
g
(
7
)
(
t
)
{\displaystyle E=\left[p_{2}(t)g^{(1)}(t)+p_{4}(t)g^{(3)}(t)+p_{6}(t)g^{(5)}(t)dt\right]_{-1}^{+1}-\int _{-1}^{+1}p_{7}(t)g^{(7)}(t)\ }
c
1
=
−
1
,
c
3
=
1
6
,
c
5
=
−
7
360
,
c
7
=
31
15
,
120
{\displaystyle c_{1}=-1,c_{3}={\frac {1}{6}},c_{5}=-{\frac {7}{360}},c_{7}={\frac {31}{15,120}}\ }
From the equation (3) given in Prof. Vu-Quoc's notes
x
(
t
k
)
=
x
k
+
x
k
+
1
2
+
t
k
h
2
,
t
k
∈
[
−
1
,
+
1
]
{\displaystyle x(t_{k})={\frac {x_{k}+x_{k+1}}{2}}+t_{k}{\frac {h}{2}}\ ,t_{k}\in [-1,+1]}
{\displaystyle \displaystyle }
, where
h
{\displaystyle h\ }
is the size of an interval, namely,
h
=
x
k
+
x
k
+
1
{\displaystyle h=x_{k}+x_{k+1}\ }
By re-arranging members
t
k
(
x
)
=
2
x
−
(
x
k
+
x
k
+
1
)
h
{\displaystyle t_{k}(x)={\frac {2x-(x_{k}+x_{k+1})}{h}}\ }
From Prof. Vu-Quoc's notes on expression for function "d"
d
¯
2
r
=
p
2
r
(
1
)
2
2
r
=
d
r
=
−
B
2
r
(
2
r
)
!
{\displaystyle {\bar {d}}_{2r}={\frac {p_{2r}(1)}{2^{2r}}}=d_{r}=-{\frac {B_{2r}}{(2r)!}}\ }
{\displaystyle \displaystyle }
p
2
(
t
)
=
c
1
t
2
2
!
+
c
3
{\displaystyle p_{2}(t)=c_{1}{\frac {t^{2}}{2!}}+c_{3}\ }
{\displaystyle \displaystyle }
p
4
(
t
)
=
c
1
t
4
4
!
+
c
3
t
2
2
!
+
c
5
{\displaystyle p_{4}(t)=c_{1}{\frac {t^{4}}{4!}}+c_{3}{\frac {t^{2}}{2!}}+c_{5}\ }
{\displaystyle \displaystyle }
p
6
(
t
)
=
c
1
t
6
6
!
+
c
3
t
4
4
!
+
c
5
t
2
2
!
+
c
7
{\displaystyle p_{6}(t)=c_{1}{\frac {t^{6}}{6!}}+c_{3}{\frac {t^{4}}{4!}}+c_{5}{\frac {t^{2}}{2!}}+c_{7}\ }
{\displaystyle \displaystyle }
d
1
=
p
2
(
1
)
2
2
=
−
1
1
2
2
!
+
1
6
2
2
=
−
1
12
{\displaystyle d_{1}={\frac {p_{2}(1)}{2^{2}}}={\frac {-1{\frac {1^{2}}{2!}}+{\frac {1}{6}}}{2^{2}}}=-{\frac {1}{12}}\ }
d
2
=
p
4
(
1
)
2
4
=
−
1
1
4
4
!
+
1
6
1
2
2
!
−
7
360
2
4
=
1
720
{\displaystyle d_{2}={\frac {p_{4}(1)}{2^{4}}}={\frac {-1{\frac {1^{4}}{4!}}+{\frac {1}{6}}{\frac {1^{2}}{2!}}-{\frac {7}{360}}}{2^{4}}}={\frac {1}{720}}}
d
3
=
p
6
(
1
)
2
6
=
−
1
1
6
6
!
+
1
6
1
4
4
!
−
7
360
1
2
2
!
+
31
15
,
120
2
6
=
−
1
30
,
240
{\displaystyle d_{3}={\frac {p_{6}(1)}{2^{6}}}={\frac {-1{\frac {1^{6}}{6!}}+{\frac {1}{6}}{\frac {1^{4}}{4!}}-{\frac {7}{360}}{\frac {1^{2}}{2!}}+{\frac {31}{15,120}}}{2^{6}}}=-{\frac {1}{30,240}}}
Miscellaneous/Matlab Code for caluclations [ edit | edit source ]
Matlab code that was composed for numerical output values
function as = here
clc
t = sym ( '1' )
c1 = - 1
c3 = 1 / 6
c5 = - 7 / 360
c7 = - ( c1 * t ^7 / factorial ( 7 ) + c3 * t ^5 / factorial ( 5 ) + c5 * t ^3 / factorial ( 3 ))
d1 = ( c1 * t ^2 / factorial ( 2 ) + c3 ) / ( 2 ^2 )
d2 = ( c1 * t ^4 / factorial ( 4 ) + c3 * t ^2 / factorial ( 2 ) + c5 ) / ( 2 ^4 )
d3 = ( c1 * t ^6 / factorial ( 6 ) + c3 * t ^4 / factorial ( 4 ) + c5 * t ^2 / factorial ( 2 ) + c7 ) / ( 2 ^6 )
Matlab output
t =
1
c1 =
-1
c3 =
0.1667
c5 =
-0.0194
c7 =
31/15120
d1 =
-1/12
d2 =
1/720
d3 =
-1/30240
--Egm6341.s11.team3.JA 03:07, 4 April 2011 (UTC)
REMARK: We did the solution on our own
Given: The Error Expression of the Trapezoidal Rule [ edit | edit source ]
E
n
T
=
I
−
T
0
(
n
)
=
∑
r
=
1
l
h
2
r
d
¯
2
r
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∑
k
=
0
n
−
1
∫
x
k
x
k
+
1
P
2
l
(
t
k
(
x
)
)
f
(
2
l
)
(
x
)
d
x
{\displaystyle E_{n}^{T}=I-T_{0}(n)=\sum _{r=1}^{l}h^{2r}{\bar {d}}_{2r}\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\sum _{k=0}^{n-1}\int _{x_{k}}^{x_{k+1}}P_{2l}(t_{k}(x))f^{(2l)}(x)dx}
(
E
q
.6
.2
.1
)
{\displaystyle \displaystyle (Eq.6.2.1)}
where:
d
¯
2
r
=
P
2
l
(
1
)
2
2
r
{\displaystyle {\bar {d}}_{2r}={\frac {P_{2l}(1)}{2^{2r}}}}
(
E
q
.6
.2
.2
)
{\displaystyle \displaystyle (Eq.6.2.2)}
Find: Derive the Given Expression [ edit | edit source ]
Derive equations 6.2.1 and 6.2.2.
From class lecture, we have the following error expression:
E
=
∑
r
=
1
l
P
2
r
(
+
1
)
[
g
(
2
r
−
1
)
(
+
1
)
−
g
(
2
r
−
1
)
(
−
1
)
]
−
∫
−
1
+
1
P
2
l
(
t
)
g
(
2
l
)
(
t
)
d
t
{\displaystyle E=\sum _{r=1}^{l}P_{2r}(+1)\left[g^{(2r-1)}(+1)-g^{(2r-1)}(-1)\right]-\int _{-1}^{+1}P_{2l}(t)g^{(2l)}(t)dt}
[4]
(
E
q
.6
.2
.3
a
)
{\displaystyle \displaystyle (Eq.6.2.3a)}
Where:
g
k
(
i
)
(
t
)
=
(
h
2
)
i
f
(
i
)
(
x
(
t
)
)
,
x
∈
[
x
k
,
x
k
+
1
]
{\displaystyle g_{k}^{(i)}(t)=\left({\frac {h}{2}}\right)^{i}f^{(i)}(x(t)),\quad x\in [x_{k},x_{k+1}]}
[5]
(
E
q
.6
.2
.4
)
{\displaystyle \displaystyle (Eq.6.2.4)}
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:
g
(
2
r
−
1
)
(
+
1
)
=
(
h
2
)
(
2
r
−
1
)
f
(
2
r
−
1
)
(
x
(
+
1
)
)
=
(
h
2
)
(
2
r
−
1
)
f
(
2
r
−
1
)
(
b
)
{\displaystyle g^{(2r-1)}(+1)=\left({\frac {h}{2}}\right)^{(2r-1)}f^{(2r-1)}(x(+1))=\left({\frac {h}{2}}\right)^{(2r-1)}f^{(2r-1)}(b)}
(
E
q
.6
.2
.5
)
{\displaystyle \displaystyle (Eq.6.2.5)}
g
(
2
r
−
1
)
(
−
1
)
=
(
h
2
)
(
2
r
−
1
)
f
(
2
r
−
1
)
(
x
(
−
1
)
)
=
(
h
2
)
(
2
r
−
1
)
f
(
2
r
−
1
)
(
a
)
{\displaystyle g^{(2r-1)}(-1)=\left({\frac {h}{2}}\right)^{(2r-1)}f^{(2r-1)}(x(-1))=\left({\frac {h}{2}}\right)^{(2r-1)}f^{(2r-1)}(a)}
(
E
q
.6
.2
.6
)
{\displaystyle \displaystyle (Eq.6.2.6)}
g
(
2
l
)
(
t
)
=
(
h
2
)
2
l
f
(
2
l
)
(
x
(
t
)
)
{\displaystyle g^{(2l)}(t)=\left({\frac {h}{2}}\right)^{2l}f^{(2l)}(x(t))}
(
E
q
.6
.2
.7
)
{\displaystyle \displaystyle (Eq.6.2.7)}
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:
E
=
∑
r
=
1
l
P
2
r
(
+
1
)
[
(
h
2
)
(
2
r
−
1
)
f
(
2
r
−
1
)
(
b
)
−
(
h
2
)
(
2
r
−
1
)
f
(
2
r
−
1
)
(
a
)
]
−
∫
−
1
+
1
P
2
l
(
t
)
(
h
2
)
2
l
f
(
2
l
)
(
x
(
t
)
)
d
t
{\displaystyle E=\sum _{r=1}^{l}P_{2r}(+1)\left[\left({\frac {h}{2}}\right)^{(2r-1)}f^{(2r-1)}(b)-\left({\frac {h}{2}}\right)^{(2r-1)}f^{(2r-1)}(a)\right]-\int _{-1}^{+1}P_{2l}(t)\left({\frac {h}{2}}\right)^{2l}f^{(2l)}(x(t))dt}
=
∑
r
=
1
l
(
h
2
)
(
2
r
−
1
)
P
2
r
(
+
1
)
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∫
−
1
+
1
P
2
l
(
t
)
f
(
2
l
)
(
x
(
t
)
)
d
t
{\displaystyle \quad =\sum _{r=1}^{l}\left({\frac {h}{2}}\right)^{(2r-1)}P_{2r}(+1)\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\int _{-1}^{+1}P_{2l}(t)f^{(2l)}(x(t))dt}
=
∑
r
=
1
l
(
h
(
2
r
)
2
(
2
r
)
)
(
h
(
−
1
)
2
(
−
1
)
)
P
2
r
(
+
1
)
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∫
−
1
+
1
P
2
l
(
t
)
f
(
2
l
)
(
x
(
t
)
)
d
t
{\displaystyle \quad =\sum _{r=1}^{l}\left({\frac {h^{(2r)}}{2^{(2r)}}}\right)\left({\frac {h^{(-1)}}{2^{(-1)}}}\right)P_{2r}(+1)\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\int _{-1}^{+1}P_{2l}(t)f^{(2l)}(x(t))dt}
=
∑
r
=
1
l
2
h
(
2
r
−
1
)
(
P
2
r
(
+
1
)
2
2
r
)
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∑
k
=
0
n
−
1
∫
x
k
x
k
−
1
P
2
l
(
t
k
(
x
)
)
f
(
2
l
)
(
x
)
d
x
{\displaystyle \quad =\sum _{r=1}^{l}2h^{(2r-1)}\left({\frac {P_{2r}(+1)}{2^{2r}}}\right)\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\sum _{k=0}^{n-1}\int _{x_{k}}^{x_{k-1}}P_{2l}(t_{k}(x))f^{(2l)}(x)dx}
(
E
q
.6
.2
.8
a
)
{\displaystyle \displaystyle (Eq.6.2.8a)}
By substituting the expression given by Equation 6.2.2 into Equation 6.2.8, the final answer is achieved:
∴
E
=
∑
r
=
1
l
2
h
(
2
r
−
1
)
d
¯
2
r
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∑
k
=
0
n
−
1
∫
x
k
x
k
−
1
P
2
l
(
t
k
(
x
)
)
f
(
2
l
)
(
x
)
d
x
{\displaystyle \therefore E=\sum _{r=1}^{l}2h^{(2r-1)}{\bar {d}}_{2r}\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\sum _{k=0}^{n-1}\int _{x_{k}}^{x_{k-1}}P_{2l}(t_{k}(x))f^{(2l)}(x)dx}
(
E
q
.6
.2
.9
a
)
{\displaystyle \displaystyle (Eq.6.2.9a)}
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:
E
=
(
h
2
)
∑
r
=
1
l
P
2
r
(
+
1
)
[
g
(
2
r
−
1
)
(
+
1
)
−
g
(
2
r
−
1
)
(
−
1
)
]
−
∫
−
1
+
1
P
2
l
(
t
)
g
(
2
l
)
(
t
)
d
t
{\displaystyle E=\left({\frac {h}{2}}\right)\sum _{r=1}^{l}P_{2r}(+1)\left[g^{(2r-1)}(+1)-g^{(2r-1)}(-1)\right]-\int _{-1}^{+1}P_{2l}(t)g^{(2l)}(t)dt}
(
E
q
.6
.2
.3
b
)
{\displaystyle \displaystyle (Eq.6.2.3b)}
Then, by following the same exact steps shown above, the derivation would instead have led to:
E
=
(
h
2
)
∑
r
=
1
l
P
2
r
(
+
1
)
[
(
h
2
)
(
2
r
−
1
)
f
(
2
r
−
1
)
(
b
)
−
(
h
2
)
(
2
r
−
1
)
f
(
2
r
−
1
)
(
a
)
]
−
∫
−
1
+
1
P
2
l
(
t
)
(
h
2
)
2
l
f
(
2
l
)
(
x
(
t
)
)
d
t
{\displaystyle E=\left({\frac {h}{2}}\right)\sum _{r=1}^{l}P_{2r}(+1)\left[\left({\frac {h}{2}}\right)^{(2r-1)}f^{(2r-1)}(b)-\left({\frac {h}{2}}\right)^{(2r-1)}f^{(2r-1)}(a)\right]-\int _{-1}^{+1}P_{2l}(t)\left({\frac {h}{2}}\right)^{2l}f^{(2l)}(x(t))dt}
=
(
h
2
)
∑
r
=
1
l
(
h
2
)
(
2
r
−
1
)
P
2
r
(
+
1
)
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∫
−
1
+
1
P
2
l
(
t
)
f
(
2
l
)
(
x
(
t
)
)
d
t
{\displaystyle \quad =\left({\frac {h}{2}}\right)\sum _{r=1}^{l}\left({\frac {h}{2}}\right)^{(2r-1)}P_{2r}(+1)\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\int _{-1}^{+1}P_{2l}(t)f^{(2l)}(x(t))dt}
=
(
h
2
)
∑
r
=
1
l
(
h
(
2
r
)
2
(
2
r
)
)
(
h
(
−
1
)
2
(
−
1
)
)
P
2
r
(
+
1
)
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∫
−
1
+
1
P
2
l
(
t
)
f
(
2
l
)
(
x
(
t
)
)
d
t
{\displaystyle \quad =\left({\frac {h}{2}}\right)\sum _{r=1}^{l}\left({\frac {h^{(2r)}}{2^{(2r)}}}\right)\left({\frac {h^{(-1)}}{2^{(-1)}}}\right)P_{2r}(+1)\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\int _{-1}^{+1}P_{2l}(t)f^{(2l)}(x(t))dt}
=
∑
r
=
1
l
h
(
2
r
)
(
P
2
r
(
+
1
)
2
2
r
)
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∑
k
=
0
n
−
1
∫
x
k
x
k
−
1
P
2
l
(
t
k
(
x
)
)
f
(
2
l
)
(
x
)
d
x
{\displaystyle \quad =\sum _{r=1}^{l}h^{(2r)}\left({\frac {P_{2r}(+1)}{2^{2r}}}\right)\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\sum _{k=0}^{n-1}\int _{x_{k}}^{x_{k-1}}P_{2l}(t_{k}(x))f^{(2l)}(x)dx}
(
E
q
.6
.2
.8
b
)
{\displaystyle \displaystyle (Eq.6.2.8b)}
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:
∴
E
=
∑
r
=
1
l
h
(
2
r
)
d
¯
2
r
[
f
(
2
r
−
1
)
(
b
)
−
f
(
2
r
−
1
)
(
a
)
]
−
(
h
2
)
2
l
∑
k
=
0
n
−
1
∫
x
k
x
k
−
1
P
2
l
(
t
k
(
x
)
)
f
(
2
l
)
(
x
)
d
x
{\displaystyle \therefore E=\sum _{r=1}^{l}h^{(2r)}{\bar {d}}_{2r}\left[f^{(2r-1)}(b)-f^{(2r-1)}(a)\right]-\left({\frac {h}{2}}\right)^{2l}\sum _{k=0}^{n-1}\int _{x_{k}}^{x_{k-1}}P_{2l}(t_{k}(x))f^{(2l)}(x)dx}
(
E
q
.6
.2
.9
b
)
{\displaystyle \displaystyle (Eq.6.2.9b)}
Egm6341.s11.team3.russo 19:35, 28 March 2011 (UTC)
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg31.djvu&page=1
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg31.djvu&page=2
↑ (Lecture 31 page 3)
↑ (meeting 31 page 3) .
↑ (meeting 30 page 3) .
x
c
o
t
h
x
=
∑
r
=
0
∞
B
2
r
(
2
r
)
!
x
2
r
{\displaystyle xcothx=\sum _{r=0}^{\infty }{\frac {{B}_{2r}}{(2r)!}}x^{2}r}
w
h
e
r
e
d
¯
2
r
=
−
B
2
r
(
2
r
)
!
{\displaystyle where\quad {\overline {d}}_{2r}=-{\frac {{B}_{2r}}{(2r)!}}}
(
E
q
.6
.3
.1
)
{\displaystyle \displaystyle (Eq.6.3.1)}
d
¯
2
r
=
P
2
r
(
1
)
2
2
r
{\displaystyle {\overline {d}}_{2r}={\frac {{P}_{2r}(1)}{{2}^{2r}}}}
(
E
q
.6
.3
.2
)
{\displaystyle \displaystyle (Eq.6.3.2)}
B
0
=
1
{\displaystyle {B}_{0}=1\quad }
B
1
=
±
1
2
{\displaystyle {B}_{1}=\pm {\frac {1}{2}}}
B
2
=
1
6
{\displaystyle {B}_{2}={\frac {1}{6}}}
B
4
=
−
1
30
{\displaystyle {B}_{4}=-{\frac {1}{30}}}
B
6
=
1
42
{\displaystyle {B}_{6}={\frac {1}{42}}}
B
8
=
−
1
30
{\displaystyle {B}_{8}=-{\frac {1}{30}}}
B
10
=
5
66
{\displaystyle {B}_{10}={\frac {5}{66}}}
1.
V
e
r
i
f
y
d
¯
2
,
d
¯
4
,
d
¯
6
{\displaystyle 1.Verify\quad {\overline {d}}_{2},{\overline {d}}_{4},{\overline {d}}_{6}}
2.
d
¯
8
,
d
¯
10
b
y
E
q
(
6.3.2
)
{\displaystyle 2.\quad {\overline {d}}_{8}\quad ,{\overline {d}}_{10}\quad by\quad Eq(6.3.2)}
by Eq(6.3.1)
d
¯
2
=
−
B
2
2
!
=
−
0.0833
{\displaystyle {\overline {d}}_{2}=-{\frac {{B}_{2}}{2!}}=-0.0833}
d
¯
4
=
−
B
4
4
!
=
0.00138
{\displaystyle {\overline {d}}_{4}=-{\frac {{B}_{4}}{4!}}=0.00138}
d
¯
6
=
−
B
6
6
!
=
−
0.000033068
{\displaystyle {\overline {d}}_{6}=-{\frac {{B}_{6}}{6!}}=-0.000033068}
d
¯
8
=
−
B
8
8
!
=
0.000000826
{\displaystyle {\overline {d}}_{8}=-{\frac {{B}_{8}}{8!}}=0.000000826}
d
¯
10
=
−
B
10
10
!
=
−
0.00000002
{\displaystyle {\overline {d}}_{10}=-{\frac {{B}_{10}}{10!}}=-0.00000002}
p
7
(
t
)
=
−
t
7
5040
+
t
5
720
−
7
t
3
2160
+
31
t
15120
{\displaystyle {p}_{7}(t)=-{\frac {t^{7}}{5040}}+{\frac {t^{5}}{720}}-{\frac {7t^{3}}{2160}}+{\frac {31t}{15120}}}
p
8
(
t
)
=
∫
p
7
(
t
)
{\displaystyle {p}_{8}(t)=\int {p}_{7}(t)}
=
−
t
8
40320
+
t
6
4320
−
7
t
4
8640
+
31
t
2
30240
+
α
{\displaystyle =-{\frac {t^{8}}{40320}}+{\frac {t^{6}}{4320}}-{\frac {7t^{4}}{8640}}+{\frac {31t^{2}}{30240}}+{\boldsymbol {\alpha }}}
p
9
(
t
)
=
∫
p
8
(
t
)
{\displaystyle {p}_{9}(t)=\int {p}_{8}(t)}
−
t
9
5040
+
t
7
17640
−
7
t
5
43200
+
31
t
3
90720
+
α
t
+
β
{\displaystyle -{\frac {t^{9}}{5040}}+{\frac {t^{7}}{17640}}-{\frac {7t^{5}}{43200}}+{\frac {31t^{3}}{90720}}+{\boldsymbol {\alpha }}t+{\boldsymbol {\beta }}}
S
i
n
c
e
p
9
(
t
)
=
O
d
d
,
p
9
(
1
)
=
0
,
p
9
(
−
1
)
=
0
,
p
9
(
0
)
=
0
{\displaystyle Since\quad {p}_{9}(t)=\quad Odd,\quad {p}_{9}(1)=0,\quad {p}_{9}(-1)=0,\quad {p}_{9}(0)=0}
Therefore
β
=
0
{\displaystyle {\boldsymbol {\beta }}=0}
α
=
−
0.000209986
{\displaystyle {\boldsymbol {\alpha }}=-0.000209986}
p
8
(
1
)
=
0.00021164
{\displaystyle {p}_{8}(1)=0.00021164\quad }
by Eq(6.3.2)
Thus
d
¯
8
=
0.00021164
2
8
=
0.000000826
{\displaystyle {\overline {d}}_{8}={\frac {0.00021164}{2^{8}}}=0.000000826}
{\displaystyle \displaystyle }
p
9
(
t
)
=
−
t
9
362880
+
t
7
17640
−
7
t
5
43200
+
31
t
3
90720
−
0.000209986
t
{\displaystyle {p}_{9}(t)=-{\frac {t^{9}}{362880}}+{\frac {t^{7}}{17640}}-{\frac {7t^{5}}{43200}}+{\frac {31t^{3}}{90720}}-0.000209986t}
p
10
(
t
)
=
∫
p
9
(
t
)
{\displaystyle {p}_{10}(t)=\int {p}_{9}(t)}
=
−
t
10
3628800
+
t
8
141120
−
7
t
4
259200
+
31
t
4
362880
+
0.000104993
t
2
+
α
{\displaystyle =-{\frac {t^{10}}{3628800}}+{\frac {t^{8}}{141120}}-{\frac {7t^{4}}{259200}}+{\frac {31t^{4}}{362880}}+0.000104993t^{2}+{\boldsymbol {\alpha }}}
p
11
(
t
)
=
∫
p
10
(
t
)
{\displaystyle {p}_{11}(t)=\int {p}_{10}(t)}
−
t
11
39916800
+
t
9
1270080
−
7
t
5
1814400
+
31
t
5
1814400
−
0.000034997
t
3
+
α
t
+
β
{\displaystyle -{\frac {{t}^{11}}{39916800}}+{\frac {t^{9}}{1270080}}-{\frac {7t^{5}}{1814400}}+{\frac {31t^{5}}{1814400}}-0.000034997{t}^{3}+{\boldsymbol {\alpha }}t+{\boldsymbol {\beta }}}
S
i
n
c
e
p
11
(
t
)
=
O
d
d
,
p
11
(
1
)
=
0
,
p
11
(
−
1
)
=
0
,
p
11
(
0
)
=
0
{\displaystyle Since\quad {p}_{11}(t)=\quad Odd,\quad {p}_{11}(1)=0,\quad {p}_{11}(-1)=0,\quad {p}_{11}(0)=0}
Therefore
β
=
0
{\displaystyle {\boldsymbol {\beta }}=0}
α
=
0.000021007
{\displaystyle {\boldsymbol {\alpha }}=0.000021007}
p
10
(
1
)
=
−
0.000018753
{\displaystyle {p}_{10}(1)=-0.000018753\quad }
d
¯
10
=
−
0.000018753
2
1
0
=
−
0.000000018
≃
−
0.00000002
{\displaystyle {\overline {d}}_{10}={\frac {-0.000018753}{2^{1}0}}=-0.000000018\simeq -0.00000002}
{\displaystyle \displaystyle }
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 ]
Refer Lecture slide http://en.wikiversity.org/w/index.php?title=File%3ANm1.s11.mtg30.djvu&page=2 for problem statement.
E
n
1
=
h
2
∑
k
=
0
n
−
1
[
∫
−
1
1
g
k
(
t
)
d
t
−
(
g
k
(
−
1
)
+
g
k
(
+
1
)
)
⏟
E
]
{\displaystyle E_{n}^{1}={\frac {h}{2}}\sum _{k=0}^{n-1}\left[\underbrace {\int _{-1}^{1}g_{k}(t)dt-(g_{k}(-1)+g_{k}(+1))} _{E}\right]}
.
where
E
{\displaystyle \displaystyle E}
can be written as,
E
=
∫
−
1
1
(
−
t
)
g
k
(
1
)
(
t
)
d
t
{\displaystyle E=\int _{-1}^{1}(-t)g_{k}^{(1)}(t)dt}
.
E
=
∫
−
1
1
(
−
t
)
⏟
p
1
(
t
)
g
k
(
1
)
(
t
)
d
t
=
[
p
2
(
t
)
g
k
(
1
)
(
t
)
]
−
1
+
1
−
∫
−
1
+
1
p
2
(
t
)
g
k
(
2
)
(
t
)
d
t
{\displaystyle \displaystyle {\begin{aligned}E&=\int _{-1}^{1}\underbrace {(-t)} _{p_{1}(t)}g_{k}^{(1)}(t)dt\\&=\left[p_{2}(t)g_{k}^{(1)}(t)\right]_{-1}^{+1}-\int _{-1}^{+1}p_{2}(t)g_{k}^{(2)}(t)dt\end{aligned}}}
.
where
p
2
(
t
)
=
∫
p
1
(
t
)
d
t
{\displaystyle \displaystyle p_{2}(t)=\int p_{1}(t)dt}
.
After successive integration by parts,
E
=
[
p
2
(
t
)
g
k
(
1
)
(
t
)
−
p
3
(
t
)
g
k
(
2
)
(
t
)
+
p
4
(
t
)
g
k
(
3
)
(
t
)
−
p
5
(
t
)
g
k
(
4
)
(
t
)
+
p
6
(
t
)
g
k
(
5
)
(
t
)
+
⋯
+
p
ℓ
(
t
)
g
k
(
ℓ
−
1
)
(
t
)
]
−
1
+
1
−
∫
−
1
1
p
ℓ
(
t
)
g
k
(
ℓ
)
(
t
)
d
t
{\displaystyle E=\left[p_{2}(t)g_{k}^{(1)}(t)-p_{3}(t)g_{k}^{(2)}(t)+p_{4}(t)g_{k}^{(3)}(t)-p_{5}(t)g_{k}^{(4)}(t)+p_{6}(t)g_{k}^{(5)}(t)+\cdots +p_{\ell }(t)g_{k}^{(\ell -1)}(t)\right]_{-1}^{+1}-\int _{-1}^{1}p_{\ell }(t)g_{k}^{(\ell )}(t)dt}
{\displaystyle \displaystyle }
Try to derive the Higher Order Trapezoidal Rule Error equation by canceling odd order derivative terms of
g
(
t
)
{\displaystyle \displaystyle g(t)}
, i.e., make terms
p
2
(
t
)
,
p
4
(
t
)
,
p
6
(
t
)
,
⋯
,
{\displaystyle \displaystyle p_{2}(t),p_{4}(t),p_{6}(t),\cdots ,}
to zero at
t
=
±
1
{\displaystyle \displaystyle t=\pm 1}
.
We know,
p
2
(
t
)
=
−
t
2
2
+
c
3
{\displaystyle p_{2}(t)=-{\frac {t^{2}}{2}}+c_{3}}
At
t
=
±
1
{\displaystyle \displaystyle t=\pm 1}
,
p
2
(
1
)
=
−
1
2
2
+
c
3
=
0
;
p
2
(
−
1
)
=
−
(
−
1
)
2
2
+
c
3
=
0
;
{\displaystyle p_{2}(1)=-{\frac {1^{2}}{2}}+c_{3}=0;p_{2}(-1)=-{\frac {(-1)^{2}}{2}}+c_{3}=0;}
c
3
=
1
2
{\displaystyle c_{3}={\frac {1}{2}}}
We also know,
p
3
(
t
)
=
−
t
3
3
!
+
t
2
+
c
4
{\displaystyle p_{3}(t)=-{\frac {t^{3}}{3!}}+{\frac {t}{2}}+c_{4}}
And,
p
4
(
t
)
=
−
t
4
4
!
+
t
2
4
+
c
4
t
+
c
5
{\displaystyle p_{4}(t)=-{\frac {t^{4}}{4!}}+{\frac {t^{2}}{4}}+c_{4}t+c_{5}}
At
t
=
±
1
{\displaystyle \displaystyle t=\pm 1}
,
p
4
(
1
)
=
−
(
1
)
4
4
!
+
(
1
)
2
4
+
c
4
(
1
)
+
c
5
;
p
4
(
−
1
)
=
−
(
−
1
)
4
4
!
+
(
−
1
)
2
4
+
c
4
(
−
1
)
+
c
5
;
{\displaystyle p_{4}(1)=-{\frac {(1)^{4}}{4!}}+{\frac {(1)^{2}}{4}}+c_{4}(1)+c_{5};p_{4}(-1)=-{\frac {(-1)^{4}}{4!}}+{\frac {(-1)^{2}}{4}}+c_{4}(-1)+c_{5};}
c
4
=
0
;
c
5
=
−
5
24
{\displaystyle \displaystyle c_{4}=0;c_{5}={\frac {-5}{24}}}
.
Functions
p
3
(
t
)
,
p
5
(
t
)
,
p
7
(
t
)
,
⋯
,
{\displaystyle \displaystyle p_{3}(t),p_{5}(t),p_{7}(t),\cdots ,}
are odd functions.
E
=
[
−
p
3
(
t
)
g
k
(
2
)
(
t
)
−
p
5
(
t
)
g
k
(
4
)
(
t
)
−
p
7
(
t
)
g
k
(
6
)
(
t
)
⋯
⋯
−
p
(
2
ℓ
−
1
)
(
t
)
g
k
(
2
ℓ
−
2
)
(
t
)
]
−
1
+
1
+
∫
−
1
1
p
(
2
ℓ
−
1
)
(
t
)
g
k
(
2
ℓ
−
1
)
(
t
)
d
t
=
∑
r
=
2
ℓ
[
−
p
(
2
r
−
1
)
g
k
(
2
r
−
2
)
]
−
1
+
1
+
∫
−
1
1
p
(
2
ℓ
−
1
)
(
t
)
g
k
(
2
ℓ
−
1
)
(
t
)
d
t
=
∑
r
=
2
ℓ
[
−
p
(
2
r
−
1
)
(
1
)
g
k
(
2
r
−
2
)
(
1
)
+
p
(
2
r
−
1
)
(
−
1
)
g
k
(
2
r
−
2
)
(
−
1
)
]
+
∫
−
1
1
p
(
2
ℓ
−
1
)
(
t
)
g
k
(
2
ℓ
−
1
)
(
t
)
d
t
{\displaystyle \displaystyle {\begin{aligned}E&=\left[-p_{3}(t)g_{k}^{(2)}(t)-p_{5}(t)g_{k}^{(4)}(t)-p_{7}(t)g_{k}^{(6)}(t)\cdots \cdots -p_{(2\ell -1)}(t)g_{k}^{(2\ell -2)}(t)\right]_{-1}^{+1}+\int _{-1}^{1}p_{(2\ell -1)}(t)g_{k}^{(2\ell -1)}(t)dt\\&=\sum _{r=2}^{\ell }\left[-p_{(2r-1)}g_{k}^{(2r-2)}\right]_{-1}^{+1}+\int _{-1}^{1}p_{(2\ell -1)}(t)g_{k}^{(2\ell -1)}(t)dt\\&=\sum _{r=2}^{\ell }\left[-p_{(2r-1)}(1)g_{k}^{(2r-2)}(1)+p_{(2r-1)}(-1)g_{k}^{(2r-2)}(-1)\right]+\int _{-1}^{1}p_{(2\ell -1)}(t)g_{k}^{(2\ell -1)}(t)dt\\\end{aligned}}}
Since
p
(
2
r
−
1
)
(
t
)
{\displaystyle \displaystyle p_{(2r-1)}(t)}
are odd functions, we know
p
(
2
r
−
1
)
(
1
)
+
p
(
2
r
−
1
)
(
−
1
)
=
0
{\displaystyle \displaystyle p_{(2r-1)}(1)+p_{(2r-1)}(-1)=0}
.
And henceforth
E
{\displaystyle \displaystyle E}
becomes
E
=
∑
r
=
2
ℓ
p
(
2
r
−
1
)
(
−
1
)
[
g
k
(
2
r
−
2
)
(
1
)
+
g
k
(
2
r
−
2
)
(
−
1
)
]
+
∫
−
1
1
p
(
2
ℓ
−
1
)
(
t
)
g
k
(
2
ℓ
−
1
)
(
t
)
d
t
{\displaystyle E=\sum _{r=2}^{\ell }p_{(2r-1)}(-1)\left[g_{k}^{(2r-2)}(1)+g_{k}^{(2r-2)}(-1)\right]+\int _{-1}^{1}p_{(2\ell -1)}(t)g_{k}^{(2\ell -1)}(t)dt}
Using the below relationship
g
k
(
i
)
(
t
)
=
(
h
2
)
i
f
(
i
)
(
x
(
t
)
)
;
x
∈
[
x
k
,
x
k
+
1
]
{\displaystyle \displaystyle g_{k}^{(i)}(t)=\left({\frac {h}{2}}\right)^{i}f^{(i)}(x(t));\ \ \ \ \ x\in [x_{k},x_{k+1}]}
.
In Equation (4)
g
k
(
t
)
{\displaystyle \displaystyle g_{k}(t)}
can be transferred to
f
(
x
k
)
{\displaystyle \displaystyle f(x_{k})}
as shown below:
g
k
(
2
r
−
2
)
(
+
1
)
=
(
h
2
)
(
2
r
−
2
)
f
(
2
r
−
2
)
(
x
k
+
1
)
,
g
k
(
2
r
−
2
)
(
−
1
)
=
(
h
2
)
(
2
r
−
2
)
f
(
2
r
−
2
)
(
x
k
)
,
{\displaystyle \displaystyle {\begin{aligned}&g_{k}^{(2r-2)}(+1)=\left({\frac {h}{2}}\right)^{(2r-2)}f^{(2r-2)}(x_{k+1}),\\&g_{k}^{(2r-2)}(-1)=\left({\frac {h}{2}}\right)^{(2r-2)}f^{(2r-2)}(x_{k}),\\\end{aligned}}}
.
E
=
∑
r
=
2
ℓ
[
(
h
2
)
(
2
r
−
2
)
p
(
2
r
−
1
)
(
−
1
)
[
f
(
2
r
−
2
)
(
x
k
+
1
)
+
f
(
2
r
−
2
)
(
x
k
)
]
]
+
(
h
2
)
(
2
ℓ
−
1
)
∫
x
k
x
k
+
1
p
(
2
ℓ
−
1
)
(
t
k
(
x
)
)
f
(
2
ℓ
−
1
)
(
x
)
d
x
{\displaystyle \displaystyle E=\sum _{r=2}^{\ell }\left[\left({\frac {h}{2}}\right)^{(2r-2)}p_{(2r-1)}(-1){\Big [}f^{(2r-2)}(x_{k+1})+f^{(2r-2)}(x_{k}){\Big ]}\right]+\left({\frac {h}{2}}\right)^{(2\ell -1)}\int _{x_{k}}^{x_{k+1}}p_{(2\ell -1)}(t_{k}(x))f^{(2\ell -1)}(x)dx}
E
n
1
=
∑
r
=
2
ℓ
[
(
h
2
)
(
2
r
−
1
)
p
(
2
r
−
1
)
(
−
1
)
∑
k
=
0
n
[
f
(
2
r
−
2
)
(
x
k
+
1
)
+
f
(
2
r
−
2
)
(
x
k
)
]
]
+
(
h
2
)
(
2
ℓ
)
∑
k
=
0
n
[
∫
−
1
1
p
(
2
ℓ
−
1
)
(
t
k
(
x
)
)
f
(
2
ℓ
−
1
)
(
x
)
d
x
]
=
∑
r
=
2
ℓ
[
(
h
2
)
(
2
r
−
1
)
p
(
2
r
−
1
)
(
−
1
)
[
f
(
2
r
−
2
)
(
x
n
)
+
f
(
2
r
−
2
)
(
x
n
−
1
)
+
⋯
⋯
+
f
(
2
r
−
2
)
(
x
0
)
β
]
]
+
(
h
2
)
(
2
ℓ
)
∑
k
=
0
n
[
∫
−
1
1
p
(
2
ℓ
−
1
)
(
t
k
(
x
)
)
f
(
2
ℓ
−
1
)
(
x
)
d
x
]
{\displaystyle \displaystyle {\begin{aligned}E_{n}^{1}&=\sum _{r=2}^{\ell }\left[\left({\frac {h}{2}}\right)^{(2r-1)}p_{(2r-1)}(-1)\sum _{k=0}^{n}{\Big [}f^{(2r-2)}(x_{k+1})+f^{(2r-2)}(x_{k}){\Big ]}\right]+\left({\frac {h}{2}}\right)^{(2\ell )}\sum _{k=0}^{n}\left[\int _{-1}^{1}p_{(2\ell -1)}(t_{k}(x))f^{(2\ell -1)}(x)dx\right]\\&=\sum _{r=2}^{\ell }\left[\left({\frac {h}{2}}\right)^{(2r-1)}p_{(2r-1)}(-1){\Big [}{f^{(2r-2)}(x_{n})+f^{(2r-2)}(x_{n-1})+\cdots \cdots +f^{(2r-2)}(x_{0})}_{\beta }{\Big ]}\right]+\left({\frac {h}{2}}\right)^{(2\ell )}\sum _{k=0}^{n}\left[\int _{-1}^{1}p_{(2\ell -1)}(t_{k}(x))f^{(2\ell -1)}(x)dx\right]\end{aligned}}}
(
E
q
.6
.2
.9
a
)
{\displaystyle \displaystyle (Eq.6.2.9a)}
Egm6341.s11.team3.rakesh 05:16, 2 April 2011 (UTC)
Recurrence Relationship
{\displaystyle \displaystyle }
c
2
i
+
1
=
−
c
1
(
2
i
+
1
)
!
−
c
3
(
2
i
−
1
)
!
−
c
5
(
2
i
−
3
)
!
−
.
.
.
−
c
2
i
−
1
3
!
{\displaystyle \displaystyle c_{2i+1}=-{\dfrac {c_{1}}{(2i+1)!}}-{\dfrac {c_{3}}{(2i-1)!}}-{\dfrac {c_{5}}{(2i-3)!}}-...-{\dfrac {c_{2i-1}}{3!}}}
Where;
{\displaystyle \displaystyle }
p
1
(
t
)
=
−
t
⇒
c
1
=
−
1
{\displaystyle \displaystyle p_{1}(t)=-t\Rightarrow c_{1}=-1}
{\displaystyle \ }
Calculate the coefficients
C
3
,
C
5
,
C
7
a
n
d
C
9
{\displaystyle C_{3},C_{5},C_{7}\;and\;C_{9}}
and the pairs of polynomials
(
P
2
,
P
3
)
,
(
P
4
,
P
5
)
,
(
P
6
,
P
7
)
;
(
P
8
,
P
9
)
{\displaystyle \left(P_{2},P_{3}\right),\;\left(P_{4},P_{5}\right),\;\left(P_{6},P_{7}\right);\left(P_{8},P_{9}\right)}
using the Recurrence formula
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.
Maple code
for i to 4 do sum ( c [ 2 * j + 1 ] * t ^( 2 * ( i - j )) / factorial ( 2 * ( i - j )), j = 0 .. i ) end do
Output
1
- c[1] + c[3]
6
1 1
--- c[1] + - c[3] + c[5]
120 6
1 1 1
---- c[1] + --- c[3] + - c[5] + c[7]
5040 120 6
1 1 1 1
------ c[1] + ---- c[3] + --- c[5] + - c[7] + c[9]
362880 5040 120 6
Given that,
i
=
1
{\displaystyle i=1\,}
c
1
3
!
+
c
3
=
0
{\displaystyle {\frac {c_{1}}{3!}}+c_{3}=0\,}
p
1
(
t
)
=
−
t
,
i
.
e
.
,
c
1
=
−
1
{\displaystyle p_{1}(t)=-t,i.e.,c_{1}=-1\,}
⇒
c
3
=
1
6
{\displaystyle \Rightarrow c_{3}={\frac {1}{6}}\,}
when
i
=
2
{\displaystyle i=2\,}
c
1
5
!
+
c
3
3
!
+
c
5
=
0
{\displaystyle {\frac {c_{1}}{5!}}+{\frac {c_{3}}{3!}}+c_{5}=0\,}
c
3
=
1
6
{\displaystyle c_{3}={\frac {1}{6}}\,}
⇒
c
5
=
−
7
360
{\displaystyle \Rightarrow c_{5}=-{\frac {7}{360}}\,}
when
i
=
3
{\displaystyle i=3\,}
c
1
7
!
+
c
3
5
!
+
c
5
3
!
+
c
7
=
0
{\displaystyle {\frac {c_{1}}{7!}}+{\frac {c_{3}}{5!}}+{\frac {c_{5}}{3!}}+c_{7}=0\,}
c
5
=
−
7
360
{\displaystyle c_{5}=-{\frac {7}{360}}\,}
⇒
c
7
=
31
15120
{\displaystyle \Rightarrow c_{7}={\frac {31}{15120}}\,}
when
i
=
4
{\displaystyle i=4\,}
c
1
9
!
+
c
3
7
!
+
c
5
5
!
+
c
7
3
!
+
c
9
=
0
{\displaystyle {\frac {c_{1}}{9!}}+{\frac {c_{3}}{7!}}+{\frac {c_{5}}{5!}}+{\frac {c_{7}}{3!}}+c_{9}=0\,}
substitute the values and solve for c_{9} ,
⇒
c
9
=
−
127
604800
{\displaystyle \Rightarrow c_{9}=-{\frac {127}{604800}}\,}
Now we calculate the pair polynomials using the expressions given below for odd and even terms , and then the constants are substituted in them ,
p
2
i
(
t
)
=
∑
j
=
0
i
c
2
j
+
1
t
2
(
i
−
j
)
[
2
(
i
−
j
)
]
!
{\displaystyle p_{2i}(t)=\sum _{j=0}^{i}c_{2j+1}{\frac {t^{2(i-j)}}{[2(i-j)]!}}\,}
p
2
i
+
1
(
t
)
=
∑
j
=
0
i
c
2
j
+
1
t
2
(
i
−
j
)
+
1
[
2
(
i
−
j
)
+
1
]
!
+
c
2
i
+
2
{\displaystyle p_{2i+1}(t)=\sum _{j=0}^{i}c_{2j+1}{\frac {t^{2(i-j)+1}}{[2(i-j)+1]!}}+c_{2i+2}\,}
The expressions are generated as follows,
Maple code for EVEN polynomials
for i to 4 do P [ 2 * i ] = sum ( c [ 2 * j + 1 ] * t ^( 2 * ( i - j )) / factorial ( 2 * ( i - j )), j = 0 .. i ) end do
Output
1 2
P[2] = - c[1] t + c[3]
2
1 4 1 2
P[4] = -- c[1] t + - c[3] t + c[5]
24 2
1 6 1 4 1 2
P[6] = --- c[1] t + -- c[3] t + - c[5] t + c[7]
720 24 2
1 8 1 6 1 4 1 2
P[8] = ----- c[1] t + --- c[3] t + -- c[5] t + - c[7] t + c[9]
40320 720 24 2
Maple code for ODD polynomials
for i to 4 do P [ 2 * i + 1 ] = sum ( c [ 2 * j + 1 ] * t ^( 2 * ( i - j ) + 1 ) / factorial ( 2 * ( i - j ) + 1 ), j = 0 .. i ) end do
Output
1 3
P[3] = - c[1] t + c[3] t
6
1 5 1 3
P[5] = --- c[1] t + - c[3] t + c[5] t
120 6
1 7 1 5 1 3
P[7] = ---- c[1] t + --- c[3] t + - c[5] t + c[7] t
5040 120 6
1 9 1 7 1 5 1 3
P[9] = ------ c[1] t + ---- c[3] t + --- c[5] t + - c[7] t + c[9] t
362880 5040 120 6
Now we get the final Pair Polynomials ,by just substituting the constants in the maple generated equations and the pairs are listed below :
p
2
=
−
t
2
2
+
1
6
,
p
3
=
−
t
3
6
+
1
t
6
{\displaystyle \displaystyle p_{2}=-{\frac {t^{2}}{2}}+{\frac {1}{6}},p_{3}=-{\frac {t^{3}}{6}}+{\frac {1t}{6}}\,}
{\displaystyle \displaystyle }
p
4
=
−
t
4
24
+
t
2
12
−
7
360
,
p
5
=
−
t
5
120
+
t
3
36
−
7
t
360
{\displaystyle \displaystyle p_{4}=-{\frac {t^{4}}{24}}+{\frac {t^{2}}{12}}-{\frac {7}{360}},p_{5}=-{\frac {t^{5}}{120}}+{\frac {t^{3}}{36}}-{\frac {7t}{360}}\,}
{\displaystyle \displaystyle }
p
6
=
−
t
6
720
+
t
4
144
−
7
t
2
720
+
31
15120
,
p
7
=
−
t
7
5040
+
t
5
720
−
7
t
3
2160
+
31
t
15120
{\displaystyle \displaystyle p_{6}=-{\frac {t^{6}}{720}}+{\frac {t^{4}}{144}}-{\frac {7t^{2}}{720}}+{\frac {31}{15120}},p_{7}=-{\frac {t^{7}}{5040}}+{\frac {t^{5}}{720}}-{\frac {7t^{3}}{2160}}+{\frac {31t}{15120}}\,}
{\displaystyle \displaystyle }
p
8
=
−
t
8
40320
+
t
6
4320
−
7
t
4
8640
+
31
t
2
30240
−
127
604800
,
p
9
=
−
t
9
362880
+
t
7
30240
−
7
t
5
43200
+
31
t
3
90720
−
127
t
604800
{\displaystyle \displaystyle p_{8}=-{\frac {t^{8}}{40320}}+{\frac {t^{6}}{4320}}-{\frac {7t^{4}}{8640}}+{\frac {31t^{2}}{30240}}-{\frac {127}{604800}},p_{9}=-{\frac {t^{9}}{362880}}+{\frac {t^{7}}{30240}}-{\frac {7t^{5}}{43200}}+{\frac {31t^{3}}{90720}}-{\frac {127t}{604800}}\,}
{\displaystyle \displaystyle }
Kessler's code is given below
Matlab code
function [c,p]= traperror ( n)
%compute the coefficients p_2, p_4, ... , p_{2n} associated with the
%trapezoidal rule error expansion. Also compute c_1, c_3, ...
f = 1 ; g = 2 ;
cn = - 1 ; cd = 1 ;
for k = 1 : n
f = f .* g .* ( g + 1 );
[ newcn , newcd ]= fracsum ( - 1 * cn , cd .* f );
cn =[ cn ; newcn ]; cd =[ cd ; newcd ];
f =[ f ; 1 ];
g =[ 2 + g ( 1 ); g ];
[ newpn , newpd ]= fracsum (( g - 1 ) .* cn , f .* cd );
pn ( k , 1 )= newpn ; pd ( k , 1 )= newpd ;
end
c = int64 ([ cn cd ]);
p = int64 ([ pn pd ]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [nsum,dsum]= fracsum ( n,d)
div = gcd ( round ( n ), round ( d ));
n = round ( n ./ div );
d = round ( d ./ div );
dsum = 1 ;
for k = 1 : length ( d )
dsum = lcm ( dsum , d ( k ));
end
nsum = dsum * sum ( n ./ d );
div = gcd ( round ( nsum ), round ( dsum ));
nsum = nsum / div ;
Find: understand Kessler's code line by line [ edit | edit source ]
Using
(
p
2
i
,
p
2
i
+
1
)
,
i
=
1
,
2
,
3
{\displaystyle (p_{2i},\ p_{2i+1}),\ i=1,2,3}
to understand Kessler's code line by line, starting with the best of Spring 2010
The Kessler's code is to compute the coefficients
p
2
n
{\displaystyle p_{2n}}
and the constants
c
n
{\displaystyle c_{n}}
where
n
=
1
,
2
,
3
,
.
.
.
,
n
{\displaystyle n=1,2,3,...,n}
, according to the following algorithm[2] :
1)
i
=
0
{\displaystyle i=0}
:
p
1
(
t
)
=
c
1
t
,
c
1
=
−
1
{\displaystyle p_{1}(t)=c_{1}t,\ c_{1}=-1}
(
E
q
.6
.6
.1
)
{\displaystyle (Eq.6.6.1)\ }
2)
i
=
1
,
2
,
3
,
.
.
.
{\displaystyle i=1,2,3,...}
:
c
2
i
+
1
=
−
∑
j
=
0
i
−
1
c
2
j
+
1
[
2
(
i
−
j
)
+
1
]
!
{\displaystyle c_{2i+1}=-\sum _{j=0}^{i-1}{\frac {c_{2j+1}}{[2(i-j)+1]!}}}
(
E
q
.6
.6
.2
)
{\displaystyle (Eq.6.6.2)\ }
p
2
i
(
t
)
=
∑
j
=
0
i
c
2
j
+
1
t
2
(
i
−
j
)
[
2
(
i
−
j
)
]
!
{\displaystyle p_{2i}(t)=\sum _{j=0}^{i}c_{2j+1}{\frac {t^{2(i-j)}}{[2(i-j)]!}}}
(
E
q
.6
.6
.3
)
{\displaystyle (Eq.6.6.3)\ }
function [c,p]= traperror ( n)
The input of this function
n
{\displaystyle n}
is a integer that ranges from 1 to
∞
{\displaystyle \infty }
. The function returns two arrays
c
{\displaystyle c}
and
p
{\displaystyle p}
, corresponding to the constants
c
n
{\displaystyle c_{n}}
and the coefficients
p
2
n
{\displaystyle p_{2n}}
.
Set initial values for f, g, cn and cd, where the initial value of cn and cd are the numerator and denominator of
c
1
=
−
1
{\displaystyle c_{1}=-1}
, respectively.
Start iteration from
k
=
1
{\displaystyle k=1}
to
k
=
n
{\displaystyle k=n}
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
f
i
k
=
f
i
k
−
1
g
i
k
−
1
(
g
i
k
−
1
+
1
)
{\displaystyle f_{i}^{k}=f_{i}^{k-1}g_{i}^{k-1}(g_{i}^{k-1}+1)}
, where
f
i
k
{\displaystyle f_{i}^{k}}
is the denominator of the factors of
c
i
{\displaystyle c_{i}}
in (Eqn 6.6.2) as
i
=
k
{\displaystyle i=k}
[ 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
c
2
k
−
1
{\displaystyle c_{2k-1}}
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
g
i
{\displaystyle g_{i}}
and
f
i
{\displaystyle f_{i}}
are respectively the numerator and denominator of the factors of
c
2
i
+
1
{\displaystyle c_{2i+1}}
in (Eqn 6.6.2) as
i
=
k
{\displaystyle i=k}
[ 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
p
2
k
{\displaystyle p_{2k}}
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
c
2
k
+
1
{\displaystyle c_{2k+1}}
and
p
2
k
{\displaystyle p_{2k}}
to the output arrays c and p
function [nsum,dsum]= fracsum ( n,d)
This function is to calculate the summation of
n
i
/
d
i
{\displaystyle n_{i}/d_{i}}
from vectors
n
{\displaystyle n}
and
d
{\displaystyle d}
. The outputs nsum and dsum returns the numerator and denominator,respectively.
n
s
u
m
d
s
u
m
=
∑
i
n
i
d
i
{\displaystyle {\frac {nsum}{dsum}}=\sum _{i}{\frac {n_{i}}{d_{i}}}}
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
n
i
{\displaystyle n_{i}}
and
d
i
{\displaystyle d_{i}}
, while the ratio
n
i
/
d
i
{\displaystyle n_{i}/d_{i}}
is preserved.
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
div = gcd ( round ( nsum ), round ( dsum ));
Egm6341.s11.team3.ren 06:31, 30 March 2011 (UTC)
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg31.djvu&page=3
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg32.djvu&page=3
REMARK: We did the solution on our own
Given: Bifolium and Ellipse Expressions [ edit | edit source ]
The expression of a bifolium:
r
(
θ
)
=
2
s
i
n
θ
c
o
s
2
θ
{\displaystyle r(\theta )=2sin\theta cos^{2}\theta \quad }
with
θ
∈
[
0
,
π
]
{\displaystyle \theta \in [0,\pi ]\quad }
(
E
q
.6
.7
.1
)
{\displaystyle \displaystyle (Eq.6.7.1)}
And an ellipse with
a
=
10
,
b
=
3
{\displaystyle a=10,b=3\quad }
:
Compute the circumferences of the bifolium and the ellipse by integrations:
C
=
∫
d
l
=
∫
[
r
2
+
(
d
r
d
θ
)
2
]
1
2
d
θ
{\displaystyle C=\int dl=\int [r^{2}+({\frac {dr}{d\theta }})^{2}]^{\frac {1}{2}}d\theta \quad }
(
E
q
.6
.7
.2
)
{\displaystyle \displaystyle (Eq.6.7.2)}
with the ellipse represented by:
r
(
θ
)
=
a
(
1
−
e
2
)
1
−
e
c
o
s
θ
{\displaystyle r(\theta )={\frac {a(1-e^{2})}{1-ecos\theta }}\quad }
(
E
q
.6
.7
.3
)
{\displaystyle \displaystyle (Eq.6.7.3)}
Also compute the circumferences of the ellipse using:
C
=
4
a
E
(
e
)
=
∫
0
π
/
2
4
a
[
1
−
e
2
s
i
n
2
θ
]
1
2
d
θ
{\displaystyle C=4aE(e)=\int _{0}^{\pi /2}4a[1-e^{2}sin^{2}\theta ]^{\frac {1}{2}}d\theta \quad }
(
E
q
.6
.7
.4
)
{\displaystyle \displaystyle (Eq.6.7.4)}
Computer the integrations using the following methods:
a) Comp. Trap.
b) Romberg.
c) Clenshaw-Curtis.
d) Fast Clenshaw-Curtis.
(1) Circumference for Bifolium: [ edit | edit source ]
From (Eq. 6.7.1), (Eq. 6.7.2), we derive the Integration as:
C
=
I
=
∫
0
π
[
(
2
s
i
n
θ
c
o
s
2
θ
)
2
+
(
2
c
o
s
3
θ
−
4
s
i
n
2
θ
c
o
s
θ
)
2
]
1
2
d
θ
{\displaystyle C=I=\int _{0}^{\pi }[(2sin\theta cos^{2}\theta )^{2}+(2cos^{3}\theta -4sin^{2}\theta cos\theta )^{2}]^{\frac {1}{2}}d\theta \quad }
(
E
q
.6
.7
.5
)
{\displaystyle \displaystyle (Eq.6.7.5)}
Using WolframAlpha, the exact solution is calculated to be [4] :
C
e
x
a
c
t
=
I
e
x
a
c
t
=
3.57777
{\displaystyle C_{exact}=I_{exact}=3.57777\quad }
(
E
q
.6
.7
.6
)
{\displaystyle \displaystyle (Eq.6.7.6)}
a)
The Comp. Trap. rule [5] is written as:
I
n
=
h
[
1
2
f
(
x
0
)
+
f
(
x
1
)
+
.
.
.
+
f
(
x
n
−
1
)
+
1
2
f
(
x
n
)
]
{\displaystyle I_{n}=h[{\frac {1}{2}}f(x_{0})+f(x_{1})+...+f(x_{n-1})+{\frac {1}{2}}f(x_{n})]\quad }
(
E
q
.6
.7
.7
)
{\displaystyle \displaystyle (Eq.6.7.7)}
b)
The Romberg method is implemented following the Richardson extrapolation [6] :
I
n
=
T
k
(
n
)
=
2
2
k
T
k
−
1
(
2
n
)
−
T
k
−
1
(
n
)
2
2
k
−
1
{\displaystyle I_{n}=T_{k}(n)={\frac {2^{2k}T_{k-1}(2n)-T_{k-1}(n)}{2^{2k}-1}}\quad }
(
E
q
.6
.7
.8
)
{\displaystyle \displaystyle (Eq.6.7.8)}
c)
The Clenshaw-Curtis quadrature method could be realized with the Trefethen's clencurt code [7] . The original integration could be written as:
I
n
=
∫
0
π
F
(
θ
)
d
θ
=
∑
i
=
0
n
F
(
θ
i
)
w
i
{\displaystyle I_{n}=\int _{0}^{\pi }F(\theta )d\theta =\sum _{i=0}^{n}F(\theta _{i})w_{i}\quad }
(
E
q
.6
.7
.9
)
{\displaystyle \displaystyle (Eq.6.7.9)}
where
θ
i
=
π
2
(
x
i
+
1
)
{\displaystyle \theta _{i}={\frac {\pi }{2}}(x_{i}+1)\quad }
, and
F
(
θ
)
{\displaystyle F(\theta )\quad }
determined by,
F
(
θ
)
=
[
c
o
s
2
θ
+
1
2
(
(
s
i
n
2
θ
)
2
+
(
3
c
o
s
2
θ
−
1
)
2
)
]
1
2
{\displaystyle F(\theta )=[{\frac {cos2\theta +1}{2}}((sin2\theta )^{2}+(3cos2\theta -1)^{2})]^{\frac {1}{2}}\quad }
(
E
q
.6
.7
.10
)
{\displaystyle \displaystyle (Eq.6.7.10)}
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:
Matlab code (1)
%Authors: Xi Xia
%--------------------------------------------------------------------------
clc ; close all ; clear all ;
%--------------------------------------------------------------------------
%1(a) Bifolium Comp. Trap.
%--------------------------------------------------------------------------
tic ;
it = 1 ;
E ( 1 ) = 1 ;
Out = fopen ( '7(1)a.txt' , 'wt' );
while ( abs ( E ( it )) > 1e-5 )
n = 2 ^it ;
h = pi / n ;
theta = [ 0 : h : pi ];
f = 2 * abs ( cos ( theta )) .* ( cos ( theta ) .^ 4 - 3 * ( cos ( theta ) .* sin ( theta )) .^ 2 + 4 * sin ( theta ) .^ 4 ) .^ ( 0.5 );
I ( it ) = h * ( f ( 1 ) + f ( n + 1 )) / 2 ;
for i = 2 : n
I ( it ) = I ( it ) + f ( i ) * h ;
end
T ( 1 , it ) = I ( it );
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 3.57777 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.5f%9.5f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
clear I ;
clear E ;
clear it ;
clear f ;
toc ;
t1 = toc ;
%--------------------------------------------------------------------------
%1(b) Bifolium Romberg
%--------------------------------------------------------------------------
Out = fopen ( '7(1)b.txt' , 'wt' );
for k = 2 : 8
fprintf ( Out , 'k=%6.1f\n' , k - 1 );
for n = 1 : 9 - k
T ( k , n ) = ( 2 ^( 2 * ( k - 1 )) * T ( k - 1 , n + 1 ) - T ( k - 1 , n )) / ( 2 ^( 2 * ( k - 1 )) - 1 );
E = abs ( T ( k , n ) - 3.57777 );
fprintf ( Out , 'n=%6.1f' , 2 ^n );
fprintf ( Out , '%9.5f%9.5f\n' , T ( k , n ), E );
end
fprintf ( Out , '\n' );
end
clear E ;
clear T ;
toc ;
t2 = toc - t1 ;
%--------------------------------------------------------------------------
%1(c) Bifolium Clenshaw-Curtis
%--------------------------------------------------------------------------
Out = fopen ( '7(1)c.txt' , 'wt' );
it = 1 ;
E ( 1 ) = 1 ;
while ( abs ( E ( it )) > 1e-5 )
n = 2 ^it ;
[ x , w ] = clencurt ( n );
theta = ( 1 + x ) * pi / 2 ;
f = 2 * abs ( cos ( theta )) .* ( cos ( theta ) .^ 4 - 3 * ( cos ( theta ) .* sin ( theta )) .^ 2 + 4 * sin ( theta ) .^ 4 ) .^ ( 0.5 );
I ( it ) = pi / 2 * w * f ;
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 3.57777 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.5f%9.5f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
clear I ;
clear E ;
clear it ;
clear theta ;
clear x ;
clear w ;
clear f ;
toc ;
t3 = toc - t1 - t2 ;
%--------------------------------------------------------------------------
%1(d) Bifolium fast CC
%--------------------------------------------------------------------------
Out = fopen ( '7(1)d.txt' , 'wt' );
it = 1 ;
E ( 1 ) = 1 ;
while ( abs ( E ( it )) > 1e-5 )
n = 2 ^it ;
[ theta , w ] = fclencurt ( n , 0 , pi );
f = 2 * abs ( cos ( theta )) .* ( cos ( theta ) .^ 4 - 3 * ( cos ( theta ) .* sin ( theta )) .^ 2 + 4 * sin ( theta ) .^ 4 ) .^ ( 0.5 );
I ( it ) = f '* w ;
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 3.57777 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.5f%9.5f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
toc ;
t4 = toc - t1 - t2 - t3 ;
The integration results for a), c) and d) are shown and compared in the following table:
The Romberg table for b) is:
(2) Circumference for Ellipse: Method1 [ edit | edit source ]
From definition of eccentricity of ellipse [9] :
e
=
1
−
b
2
a
2
=
91
10
{\displaystyle e={\sqrt {1-{\frac {b^{2}}{a^{2}}}}}={\frac {\sqrt {91}}{10}}\quad }
(
E
q
.6
.7
.11
)
{\displaystyle \displaystyle (Eq.6.7.11)}
From (Eq. 6.7.2), (Eq. 6.7.3), we derive the Integration as:
C
=
I
=
∫
0
2
π
a
(
1
−
e
2
)
(
1
−
e
c
o
s
θ
)
−
2
(
1
+
e
2
−
2
e
c
o
s
θ
)
d
θ
{\displaystyle C=I=\int _{0}^{2\pi }a(1-e^{2})(1-ecos\theta )^{-2}{\sqrt {(1+e^{2}-2ecos\theta )}}d\theta \quad }
=
∫
0
2
π
9
10
(
1
−
91
10
c
o
s
θ
)
−
2
(
191
100
−
91
5
c
o
s
θ
)
1
2
d
θ
{\displaystyle =\int _{0}^{2\pi }{\frac {9}{10}}(1-{\frac {\sqrt {91}}{10}}cos\theta )^{-2}({\frac {191}{100}}-{\frac {\sqrt {91}}{5}}cos\theta )^{\frac {1}{2}}d\theta \quad }
(
E
q
.6
.7
.12
)
{\displaystyle \displaystyle (Eq.6.7.12)}
Using WolframAlpha, the exact solution is calculated to be [10] :
C
e
x
a
c
t
=
I
e
x
a
c
t
=
43.8591
{\displaystyle C_{exact}=I_{exact}=43.8591\quad }
(
E
q
.6
.7
.13
)
{\displaystyle \displaystyle (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:
Matlab code (2)
%Authors: Xi Xia
%--------------------------------------------------------------------------
clc ; close all ; clear all ;
%--------------------------------------------------------------------------
%2(a) Ellipse Comp. Trap.
%--------------------------------------------------------------------------
tic ;
it = 1 ;
E ( 1 ) = 1 ;
Out = fopen ( '7(2)a.txt' , 'wt' );
while ( abs ( E ( it )) > 1e-4 )
n = 2 ^it ;
h = 2 * pi / n ;
theta = [ 0 : h : 2 * pi ];
f = 0.9 * (( 1 - ( sqrt ( 91 ) / 10 ) * cos ( theta )) .^ ( - 2 )) .* ( 1.91 - ( sqrt ( 91 ) / 5 ) * cos ( theta )) .^ ( 0.5 );
I ( it ) = h * ( f ( 1 ) + f ( n + 1 )) / 2 ;
for i = 2 : n
I ( it ) = I ( it ) + f ( i ) * h ;
end
T ( 1 , it ) = I ( it );
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 43.8591 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.4f%9.4f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
clear I ;
clear E ;
clear it ;
clear f ;
toc ;
t1 = toc ;
%--------------------------------------------------------------------------
%2(b) Ellipse Romberg
%--------------------------------------------------------------------------
Out = fopen ( '7(2)b.txt' , 'wt' );
for k = 2 : 7
fprintf ( Out , 'k=%6.1f\n' , k - 1 );
for n = 2 : 9 - k
T ( k , n ) = ( 2 ^( 2 * ( k - 1 )) * T ( k - 1 , n + 1 ) - T ( k - 1 , n )) / ( 2 ^( 2 * ( k - 1 )) - 1 );
E = abs ( T ( k , n ) - 43.8591 );
fprintf ( Out , 'n=%6.1f' , 2 ^n );
fprintf ( Out , '%9.4f%9.4f\n' , T ( k , n ), E );
end
fprintf ( Out , '\n' );
end
clear E ;
clear T ;
toc ;
t2 = toc - t1 ;
%--------------------------------------------------------------------------
%2(c) Ellipse Clenshaw-Curtis
%--------------------------------------------------------------------------
Out = fopen ( '7(2)c.txt' , 'wt' );
it = 1 ;
E ( 1 ) = 1 ;
while ( abs ( E ( it )) > 1e-4 )
n = 2 ^it ;
[ x , w ] = clencurt ( n );
theta = ( 1 + x ) * pi ;
f = 0.9 * (( 1 - ( sqrt ( 91 ) / 10 ) * cos ( theta )) .^ ( - 2 )) .* ( 1.91 - ( sqrt ( 91 ) / 5 ) * cos ( theta )) .^ ( 0.5 );
I ( it ) = pi * w * f ;
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 43.8591 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.4f%9.4f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
clear I ;
clear E ;
clear it ;
clear theta ;
clear x ;
clear w ;
clear f ;
toc ;
t3 = toc - t1 - t2 ;
%--------------------------------------------------------------------------
%2(d) Ellipse fast CC
%--------------------------------------------------------------------------
Out = fopen ( '7(2)d.txt' , 'wt' );
it = 1 ;
E ( 1 ) = 1 ;
while ( abs ( E ( it )) > 1e-4 )
n = 2 ^it ;
[ theta , w ] = fclencurt ( n , 0 , 2 * pi );
f = 0.9 * (( 1 - ( sqrt ( 91 ) / 10 ) * cos ( theta )) .^ ( - 2 )) .* ( 1.91 - ( sqrt ( 91 ) / 5 ) * cos ( theta )) .^ ( 0.5 );
I ( it ) = f '* w ;
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 43.8591 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.4f%9.4f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
toc ;
t4 = toc - t1 - t2 - t3 ;
The integration results for a), c) and d) are shown and compared in the following table:
The Romberg table for b) is:
(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:
f
(
θ
)
=
4
a
[
1
−
e
2
s
i
n
2
θ
]
1
2
=
40
1
−
0.91
s
i
n
2
θ
{\displaystyle f(\theta )=4a[1-e^{2}sin^{2}\theta ]^{\frac {1}{2}}=40{\sqrt {1-0.91sin^{2}\theta }}\quad }
(
E
q
.6
.7
.14
)
{\displaystyle \displaystyle (Eq.6.7.14)}
Using WolframAlpha, the exact solution is calculated to be [11] :
C
e
x
a
c
t
=
I
e
x
a
c
t
=
43.8591
{\displaystyle C_{exact}=I_{exact}=43.8591\quad }
(
E
q
.6
.7
.13
)
{\displaystyle \displaystyle (Eq.6.7.13)}
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:
Matlab code (3)
%Authors: Xi Xia
%--------------------------------------------------------------------------
clc ; close all ; clear all ;
%--------------------------------------------------------------------------
%3(a) Ellipse Comp. Trap.
%--------------------------------------------------------------------------
tic ;
it = 1 ;
E ( 1 ) = 1 ;
Out = fopen ( '7(3)a.txt' , 'wt' );
while ( abs ( E ( it )) > 1e-4 )
n = 2 ^it ;
h = pi / 2 / n ;
theta = [ 0 : h : 2 * pi ];
f = 40 * ( 1 - 0.91 * ( sin ( theta )) .^ 2 ) .^ ( 0.5 );
I ( it ) = h * ( f ( 1 ) + f ( n + 1 )) / 2 ;
for i = 2 : n
I ( it ) = I ( it ) + f ( i ) * h ;
end
T ( 1 , it ) = I ( it );
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 43.8591 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.4f%9.4f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
clear I ;
clear E ;
clear it ;
clear f ;
toc ;
t1 = toc ;
%--------------------------------------------------------------------------
%3(b) Ellipse Romberg
%--------------------------------------------------------------------------
Out = fopen ( '7(3)b.txt' , 'wt' );
for k = 2 : 3
fprintf ( Out , 'k=%6.1f\n' , k - 1 );
for n = 1 : 4 - k
T ( k , n ) = ( 2 ^( 2 * ( k - 1 )) * T ( k - 1 , n + 1 ) - T ( k - 1 , n )) / ( 2 ^( 2 * ( k - 1 )) - 1 );
E = abs ( T ( k , n ) - 43.8591 );
fprintf ( Out , 'n=%6.1f' , 2 ^n );
fprintf ( Out , '%9.4f%9.4f\n' , T ( k , n ), E );
end
fprintf ( Out , '\n' );
end
clear E ;
clear T ;
toc ;
t2 = toc - t1 ;
%--------------------------------------------------------------------------
%3(c) Ellipse Clenshaw-Curtis
%--------------------------------------------------------------------------
Out = fopen ( '7(3)c.txt' , 'wt' );
it = 1 ;
E ( 1 ) = 1 ;
while ( abs ( E ( it )) > 1e-4 )
n = 2 ^it ;
[ x , w ] = clencurt ( n );
theta = ( 1 + x ) * pi / 4 ;
f = 40 * ( 1 - 0.91 * ( sin ( theta )) .^ 2 ) .^ ( 0.5 );
I ( it ) = pi / 4 * w * f ;
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 43.8591 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.4f%9.4f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
clear I ;
clear E ;
clear it ;
clear theta ;
clear x ;
clear w ;
clear f ;
toc ;
t3 = toc - t1 - t2 ;
%--------------------------------------------------------------------------
%3(d) Ellipse fast CC
%--------------------------------------------------------------------------
Out = fopen ( '7(3)d.txt' , 'wt' );
it = 1 ;
E ( 1 ) = 1 ;
while ( abs ( E ( it )) > 1e-4 )
n = 2 ^it ;
[ theta , w ] = fclencurt ( n , 0 , pi / 2 );
f = 40 * ( 1 - 0.91 * ( sin ( theta )) .^ 2 ) .^ ( 0.5 );
I ( it ) = f '* w ;
it = it + 1 ;
E ( it ) = abs ( I ( it - 1 ) - 43.8591 );
fprintf ( Out , 'n=%6.1f\n' , n );
fprintf ( Out , '%9.4f%9.4f\n\n' , I ( it - 1 ), E ( it ));
end
fclose ( Out );
toc ;
t4 = toc - t1 - t2 - t3 ;
The integration results for a), c) and d) are shown and compared in the following table:
The Romberg table for b) is:
Egm6341.s11.team3.Xia 18:00, 31 March 2011 (UTC)
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg33.djvu&page=1
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg33.djvu&page=2
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg33.djvu&page=3
↑ 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
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg7.djvu&page=4
↑ http://en.wikiversity.org/w/index.php?title=File:Nm1.s11.mtg29.djvu&page=5
↑ http://people.maths.ox.ac.uk/trefethen/clencurt.m
↑ http://www.mathworks.com/matlabcentral/fileexchange/6911-fast-clenshaw-curtis-quadrature
↑ http://mathworld.wolfram.com/Ellipse.html
↑ 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
↑ http://www.wolframalpha.com/input/?i=integrate+40sqrt%281-0.91sin%5E2x%29+for+x+from+0+to+pi%2F2
REMARK: We did the solution on our own
A
r
c
L
e
n
g
t
h
(
P
Q
)
=
∫
θ
P
θ
Q
[
r
2
+
(
d
r
d
θ
)
2
]
1
/
2
d
θ
{\displaystyle ArcLength(PQ)=\int _{\theta _{P}}^{\theta _{Q}}\left[r^{2}+\left({\frac {dr}{d\theta }}\right)^{2}\right]^{1/2}d\theta }
(
E
q
.6
.8
.1
)
{\displaystyle \displaystyle (Eq.6.8.1)}
Find: Derive the Given Expression [ edit | edit source ]
Use Triangle AOB and the Law of Cosines to derive Equation 6.8.1
The Law of Cosines [3] is as follows:
c
2
=
a
2
+
b
2
−
2
a
b
∗
c
o
s
(
γ
)
{\displaystyle c^{2}=a^{2}+b^{2}-2ab*cos(\gamma )\ }
(
E
q
.6
.8
.2
)
{\displaystyle \displaystyle (Eq.6.8.2)}
Where a, b and c are the sides of a triangle and
γ
{\displaystyle \gamma }
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):
(
A
B
¯
)
2
=
(
A
A
′
¯
)
2
+
(
A
′
B
¯
)
2
−
2
(
A
A
′
¯
)
(
A
′
B
¯
)
c
o
s
(
90
)
=
(
A
A
′
¯
)
2
+
(
A
′
B
¯
)
2
{\displaystyle ({\bar {AB}})^{2}=({\bar {AA'}})^{2}+({\bar {A'B}})^{2}-2({\bar {AA'}})({\bar {A'B}})cos(90)=({\bar {AA'}})^{2}+({\bar {A'B}})^{2}}
(
E
q
.6
.8
.3
)
{\displaystyle \displaystyle (Eq.6.8.3)}
In order to achieve the desired results, we must denote the triangular components in terms of elliptical components:
O
A
¯
=
r
(
θ
)
O
B
¯
=
r
(
θ
+
d
θ
)
A
A
′
¯
=
r
(
θ
)
d
θ
A
′
B
¯
=
d
r
(
θ
)
→
A
′
B
¯
=
d
r
(
θ
)
d
θ
d
θ
A
B
¯
=
d
l
{\displaystyle {\bar {OA}}=r(\theta )\quad {\bar {OB}}=r(\theta +d\theta )\quad {\bar {AA'}}=r(\theta )d\theta \quad {\bar {A'B}}=dr(\theta )\rightarrow {\bar {A'B}}={\frac {dr(\theta )}{d\theta }}d\theta \quad {\bar {AB}}=dl}
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:
(
d
l
)
2
=
(
r
(
θ
)
d
θ
)
2
+
(
d
r
(
θ
)
d
θ
d
θ
)
2
=
(
d
θ
)
2
[
r
2
(
θ
)
+
(
d
r
(
θ
)
d
θ
)
2
]
{\displaystyle (dl)^{2}=(r(\theta )d\theta )^{2}+\left({\frac {dr(\theta )}{d\theta }}d\theta \right)^{2}=(d\theta )^{2}\left[r^{2}(\theta )+\left({\frac {dr(\theta )}{d\theta }}\right)^{2}\right]}
(
E
q
.6
.8
.4
)
{\displaystyle \displaystyle (Eq.6.8.4)}
Solving for our desired parameter dl yields:
d
l
=
[
(
d
θ
)
2
[
(
r
(
θ
)
)
2
+
(
d
r
(
θ
)
d
θ
)
2
]
]
1
/
2
=
(
d
θ
)
[
(
r
(
θ
)
)
2
+
(
d
r
(
θ
)
d
θ
)
2
]
1
/
2
{\displaystyle dl=\left[(d\theta )^{2}\left[(r(\theta ))^{2}+\left({\frac {dr(\theta )}{d\theta }}\right)^{2}\right]\right]^{1/2}=(d\theta )\left[(r(\theta ))^{2}+\left({\frac {dr(\theta )}{d\theta }}\right)^{2}\right]^{1/2}}
(
E
q
.6
.8
.5
)
{\displaystyle \displaystyle (Eq.6.8.5)}
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
θ
Q
{\displaystyle \theta _{Q}}
and
θ
P
{\displaystyle \theta _{P}}
, respectively, the arc length may be found:
A
r
c
L
e
n
g
t
h
(
P
Q
)
=
∫
θ
P
θ
Q
d
l
=
∫
θ
P
θ
Q
(
(
d
θ
)
[
(
r
(
θ
)
)
2
+
(
d
r
(
θ
)
d
θ
)
2
]
1
/
2
)
{\displaystyle ArcLength(PQ)=\int _{\theta _{P}}^{\theta _{Q}}dl=\int _{\theta _{P}}^{\theta _{Q}}\left((d\theta )\left[(r(\theta ))^{2}+\left({\frac {dr(\theta )}{d\theta }}\right)^{2}\right]^{1/2}\right)}
(
E
q
.6
.8
.6
)
{\displaystyle \displaystyle (Eq.6.8.6)}
Which is exactly equivalent to the expression given in the problem statement by Equation 6.8.1.
∴
A
r
c
L
e
n
g
t
h
(
P
Q
)
=
∫
θ
P
θ
Q
[
r
2
+
(
d
r
d
θ
)
2
]
1
/
2
d
θ
{\displaystyle \therefore ArcLength(PQ)=\int _{\theta _{P}}^{\theta _{Q}}\left[r^{2}+\left({\frac {dr}{d\theta }}\right)^{2}\right]^{1/2}d\theta }
Egm6341.s11.team3.russo 22:38, 28 March 2011 (UTC)
↑ (Lecture 33 page 2)
↑ (Lecture 32 page 4)
↑ (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
N
i
(
s
)
{\displaystyle {N_{i}(s)}\,}
,
i
=
1
,
2
,
3
,
4
{\displaystyle i=1,2,3,4\,}
. Plot
N
i
¯
{\displaystyle {\bar {N_{i}}}\,}
's= [ edit | edit source ]
[1] ,[2]
6.9.a) The matrix is given to be:
[
1
0
0
0
0
1
0
0
1
1
1
1
0
1
2
3
]
.
[
C
0
C
1
C
2
C
3
]
[
z
i
z
i
′
z
i
+
1
z
i
+
1
′
]
{\displaystyle {{\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}}}
.
{\displaystyle \displaystyle }
d
3
¯
=
z
i
+
1
=
z
(
s
=
1
)
{\displaystyle \displaystyle {\bar {d_{3}}}=z_{i+1}=z(s=1)}
.
and
{\displaystyle \displaystyle }
d
4
¯
=
z
i
+
1
′
=
z
′
(
s
=
1
)
{\displaystyle \displaystyle {\bar {d_{4}}}=z'_{i+1}=z'(s=1)}
.
We know that
{\displaystyle \displaystyle }
z
(
s
)
=
∑
i
=
0
3
C
i
s
i
=
C
0
+
C
1
s
+
C
2
s
2
+
C
3
s
3
{\displaystyle \displaystyle z(s)=\sum _{i=0}^{3}C_{i}s^{i}=C_{0}+C_{1}s+C_{2}s^{2}+C_{3}s^{3}}
.
and
{\displaystyle \displaystyle }
z
′
(
s
)
=
∑
i
=
1
3
i
C
i
s
i
−
1
=
C
1
+
2
C
2
s
+
3
C
3
s
2
{\displaystyle \displaystyle z'(s)=\sum _{i=1}^{3}iC_{i}s^{i-1}=C_{1}+2C_{2}s+3C_{3}s^{2}}
.
thus using the above relations where
s
=
1
{\displaystyle s=1}
in both cases we then have
{\displaystyle \displaystyle }
z
i
+
1
=
C
0
+
C
1
+
C
2
+
C
3
{\displaystyle \displaystyle z_{i+1}=C_{0}+C_{1}+C_{2}+C_{3}}
.
and
{\displaystyle \displaystyle }
z
i
+
1
′
=
C
1
+
2
C
2
+
3
C
3
{\displaystyle \displaystyle z'_{i+1}=C_{1}+2C_{2}+3C_{3}}
.
Hence we can put forth that the 3rd and 4th rows are
[
1111
]
{\displaystyle [1111]}
and
[
0123
]
{\displaystyle [0123]}
respectively.
6.9.b)
From the matrix
A
=
[
1
0
0
0
0
1
0
0
1
1
1
1
0
1
2
3
]
{\displaystyle A={\begin{bmatrix}1&0&0&0\\0&1&0&0\\1&1&1&1\\0&1&2&3\end{bmatrix}}}
.
Verify the inverse of the above given matrix
A
−
1
=
[
1
0
0
0
0
1
0
0
−
3
−
2
3
−
1
2
1
−
2
1
]
{\displaystyle A^{-1}={\begin{bmatrix}1&0&0&0\\0&1&0&0\\-3&-2&3&-1\\2&1&-2&1\end{bmatrix}}}
.
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
N
i
(
s
)
{\displaystyle {N_{i}(s)}\,}
,
i
=
1
,
2
,
3
,
4
{\displaystyle i=1,2,3,4\,}
. Plot
N
i
¯
{\displaystyle {\bar {N_{i}}}\,}
's [ edit | edit source ]
∑
i
=
0
3
C
i
s
i
=
∑
i
=
1
4
N
i
(
s
)
d
i
{\displaystyle \sum _{i=0}^{3}C_{i}s^{i}=\sum _{i=1}^{4}N_{i}(s)d_{i}\,}
C
0
+
C
1
s
+
C
2
s
2
+
C
3
s
3
=
N
1
(
s
)
d
1
+
N
2
(
s
)
d
2
+
N
3
(
s
)
d
3
+
N
4
(
s
)
d
4
{\displaystyle C_{0}+C_{1}s+C_{2}s^{2}+C_{3}s^{3}=N_{1}(s)d_{1}+N_{2}(s)d_{2}+N_{3}(s)d_{3}+N_{4}(s)d_{4}\,}
C
0
=
z
i
{\displaystyle C_{0}=z_{i}\,}
C
1
=
z
i
′
{\displaystyle C_{1}=z_{i}'\,}
C
2
=
−
3
z
i
−
2
z
i
′
+
3
z
i
+
1
−
z
i
+
1
′
{\displaystyle C_{2}=-3z_{i}-2z_{i}'+3z_{i+1}-z_{i+1}'\,}
C
3
=
2
z
i
+
z
i
′
−
2
z
i
+
1
+
z
i
+
1
′
{\displaystyle C_{3}=2z_{i}+z_{i}'-2z_{i+1}+z_{i+1}'\,}
d
1
=
z
i
{\displaystyle d_{1}=z_{i}\,}
d
2
=
z
i
˙
=
z
i
′
1
h
{\displaystyle d_{2}={\dot {z_{i}}}=z_{i}'{\frac {1}{h}}\,}
d
3
=
z
i
+
1
{\displaystyle d_{3}=z_{i+1}\,}
d
4
=
z
i
+
1
˙
=
z
i
+
1
′
1
h
{\displaystyle d_{4}={\dot {z_{i+1}}}=z_{i+1}'{\frac {1}{h}}\,}
z
i
+
z
i
′
s
+
(
−
3
z
i
−
2
z
i
′
+
3
z
i
+
1
−
z
i
+
1
′
)
s
2
+
(
2
z
i
+
z
i
′
−
2
z
i
+
1
+
z
i
+
1
′
)
s
3
=
N
1
(
s
)
z
i
+
N
2
(
s
)
z
i
′
1
h
+
N
3
(
s
)
z
i
+
1
+
N
4
(
s
)
z
i
+
1
′
1
h
{\displaystyle z_{i}+z_{i}'s+(-3z_{i}-2z_{i}'+3z_{i+1}-z_{i+1}')s^{2}+(2z_{i}+z_{i}'-2z_{i+1}+z_{i+1}')s^{3}=N_{1}(s)z_{i}+N_{2}(s)z_{i}'{\frac {1}{h}}+N_{3}(s)z_{i+1}+N_{4}(s)z_{i+1}'{\frac {1}{h}}\,}
z
i
(
1
−
3
s
2
+
2
s
3
)
+
z
i
′
(
s
−
2
s
2
+
s
3
)
+
z
i
+
1
(
3
s
2
−
2
s
3
)
+
z
i
+
1
′
(
−
s
2
+
s
3
)
=
z
i
(
N
1
(
s
)
)
+
z
i
′
(
N
2
(
s
)
1
h
)
+
z
i
+
1
(
N
3
(
s
)
)
+
z
i
+
1
′
(
N
4
(
s
)
1
h
)
{\displaystyle z_{i}(1-3s^{2}+2s^{3})+z_{i}'(s-2s^{2}+s^{3})+z_{i+1}(3s^{2}-2s^{3})+z_{i+1}'(-s^{2}+s^{3})=z_{i}(N_{1}(s))+z_{i}'(N_{2}(s){\frac {1}{h}})+z_{i+1}(N_{3}(s))+z_{i+1}'(N_{4}(s){\frac {1}{h}})\,}
Note for Verifying that :
N
1
¯
(
0
)
=
1
,
N
1
¯
˙
(
0
)
=
0
,
N
1
¯
(
1
)
=
0
,
N
1
¯
˙
(
1
)
=
0
{\displaystyle {\bar {N_{1}}}(0)=1,{\dot {\bar {N_{1}}}}(0)=0,{\bar {N_{1}}}(1)=0,{\dot {\bar {N_{1}}}}(1)=0\,}
N
2
¯
(
0
)
=
0
,
N
2
¯
˙
(
0
)
=
1
,
N
2
¯
(
1
)
=
0
,
N
2
¯
˙
(
1
)
=
0
{\displaystyle {\bar {N_{2}}}(0)=0,{\dot {\bar {N_{2}}}}(0)=1,{\bar {N_{2}}}(1)=0,{\dot {\bar {N_{2}}}}(1)=0\,}
N
3
¯
(
0
)
=
0
,
N
3
¯
˙
(
0
)
=
0
,
N
3
¯
(
1
)
=
1
,
N
3
¯
˙
(
1
)
=
0
{\displaystyle {\bar {N_{3}}}(0)=0,{\dot {\bar {N_{3}}}}(0)=0,{\bar {N_{3}}}(1)=1,{\dot {\bar {N_{3}}}}(1)=0\,}
N
4
¯
(
0
)
=
0
,
N
4
¯
˙
(
0
)
=
0
,
N
4
¯
(
1
)
=
0
,
N
4
¯
˙
(
1
)
=
1
{\displaystyle {\bar {N_{4}}}(0)=0,{\dot {\bar {N_{4}}}}(0)=0,{\bar {N_{4}}}(1)=0,{\dot {\bar {N_{4}}}}(1)=1\,}
{\displaystyle \displaystyle }
Matlab code that was composed for numerical output values
clc ;
s = linspace ( 0 , 1 , 1000 );
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 )
plot ( s , N1 );
subplot ( 2 , 2 , 2 )
plot ( s , N2 );
subplot ( 2 , 2 , 3 )
plot ( s , N3 );
subplot ( 2 , 2 , 4 )
plot ( s , N4 );
Hv
Plot of N1, N2, N3, N4 as a function s
Egm6341.s11.team3.rakesh 19:21, 4 April 2011 (UTC)
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