590 lines
16 KiB
Typst
590 lines
16 KiB
Typst
#import "@preview/cetz:0.3.2";
|
|
#import "@preview/cetz-plot:0.1.1": plot
|
|
#import "@preview/physica:0.9.4": *
|
|
#import "@preview/plotsy-3d:0.1.0": plot-3d-parametric-surface
|
|
#import "@preview/fletcher:0.5.4" as fletcher: diagram, edge, node
|
|
|
|
#set page(paper: "a4", margin: (x: 2.6cm, y: 2.8cm), numbering: "1 : 1")
|
|
#set par(justify: true, leading: 0.52em)
|
|
|
|
#let FONT_SIZE = 18pt;
|
|
#set text(font: "FreeSerif", size: FONT_SIZE, lang: "us")
|
|
#show math.equation: set text(font: "Euler Math", size: (FONT_SIZE * 1.0), lang: "en")
|
|
|
|
#set heading(numbering: none)
|
|
#show heading.where(level: 1): it => {
|
|
rect(inset: FONT_SIZE / 2)[#it]
|
|
}
|
|
|
|
#align(center)[
|
|
#text(size: FONT_SIZE * 2, weight: "bold")[#underline[exercise 3]]
|
|
]
|
|
|
|
these are my solutions to the third exercise set of TMA4135.
|
|
|
|
there should be a python source code file attached to this deliverable.
|
|
|
|
this document was created using
|
|
#link("https://typst.app/")[#text(blue.darken(5%))[typst]].
|
|
|
|
#v(1fr)
|
|
|
|
#outline(title: none)
|
|
|
|
|
|
#pagebreak()
|
|
|
|
= problem 1
|
|
|
|
== a) & b)
|
|
|
|
using the trapezoidal rule we can estimate
|
|
$
|
|
(1) quad I_[0,1][exp] := integral_0^1 e^x dd(x), wide
|
|
(2) quad I_[0,1][ln] := integral_0^1 ln(x) dd(x)
|
|
$
|
|
|
|
let us solve both analytically first, starting with $(1)$
|
|
$
|
|
integral_0^1 e^x dd(x) = evaluated(e^x)_0^1 = e - 1
|
|
$
|
|
then $(2)$
|
|
$
|
|
integral_0^1 ln(x) dd(x) & = [x ln(x) - x]_0^1 \
|
|
& = (1 ln(1) - 1) - lim_(t -> 0^+) (t ln(t) - t) \
|
|
& = -1 - 0 = -1,
|
|
$
|
|
where we used integration by parts and the fact that $lim_(t -> 0^+) t ln(t) = 0$.
|
|
|
|
now we know what to expect from our estimates.
|
|
|
|
recall the trapezoidal rule
|
|
$
|
|
I_[a, b][f] approx T_[a, b][f] := (b-a)/2 (f(a) + f(b))
|
|
$
|
|
|
|
in our case $a = 0$ and $b = 1$. for $(1)$ with $f(x) = e^x$ we obtain
|
|
$
|
|
T_[0, 1][exp] = 1/2 (e^0 + e^1) = 1/2 (1 + e) = (1 + e)/2,
|
|
$
|
|
which yields an error of
|
|
$
|
|
E_[0,1][exp] = (e - 1) - (1 + e)/2 = (2e - 2 - 1 - e)/2 = (e - 3)/2.
|
|
$
|
|
|
|
for $(2)$ with $f(x) = ln(x)$ we obtain
|
|
$
|
|
T_[0, 1][ln] = 1/2 (ln(0) + ln(1)) = 1/2 (-oo + 0) = -oo
|
|
$
|
|
which yields an error of $+oo$.
|
|
|
|
thus we can see that using the trapezoidal rule for $ln(x)$ diverges and proves
|
|
to be unwieldy. for $exp(x)$ the actual error is $(e-3)/2 approx #{ calc.round((calc.e - 3) / 2, digits: 3) }$.
|
|
|
|
== c)
|
|
|
|
recall the upper bound of the error for the trapezoidal rule
|
|
$
|
|
abs(E_[a,b][f]) <= (b-a)^3 / 12 dot max_(xi in [a,b]) abs(f''(xi))
|
|
$
|
|
which is for $f = exp$ on $[0,1]$
|
|
$
|
|
abs(E_[0, 1][exp]) <= 1^3/12 dot max_(xi in [0,1]) abs(e^xi) = e/12 approx #{ calc.round(calc.e / 12, digits: 3) }
|
|
$
|
|
Our actual error magnitude is $abs((e-3)/2) = (3-e)/2 approx #{ calc.round((3 - calc.e) / 2, digits: 3) }$.
|
|
|
|
Since $#{ calc.round((3 - calc.e) / 2, digits: 3) } <= #{ calc.round(calc.e / 12, digits: 3) }$, the error bound is satisfied.
|
|
|
|
== d)
|
|
|
|
recall the error bound derived in the lecture,
|
|
$
|
|
abs(E_[0,1][f]) <= (max_(xi in
|
|
[a,b]) abs(f^((n+1))(xi))) / ((n+1)!) integral_a^b product_(i=0)^n abs(
|
|
x
|
|
- x_i
|
|
)
|
|
dd(x),
|
|
$
|
|
where the integral of $f$ is taken over the interval $[a,b]$, and $f in
|
|
C^(n+1)[0,1]$.
|
|
|
|
for the trapezoidal rule, $n = 1$, so we have $x_0 = 0$ and $x_1 = 1$. for $f(x)
|
|
= e^x$:
|
|
|
|
- $f^((n+1))(x) = f^((2))(x) = e^x$
|
|
- $max_(xi in [0,1]) abs(e^xi) = e$
|
|
- $(n+1)! = 2! = 2$
|
|
|
|
the integral becomes:
|
|
$
|
|
integral_0^1 product_(i=0)^1 abs(x - x_i) dd(x) & = integral_0^1 abs(x - 0) dot abs(x - 1) dd(x) \
|
|
& = integral_0^1 x dot abs(x - 1) dd(x) \
|
|
& = integral_0^1 x dot (1 - x) dd(x) quad "since" x - 1 <= 0 "on" [0,1] \
|
|
& = integral_0^1 (x - x^2) dd(x) \
|
|
& = [x^2/2 - x^3/3]_0^1 = 1/2 - 1/3 = 1/6
|
|
$
|
|
|
|
therefore, the error bound is:
|
|
$
|
|
abs(E_[0,1][exp]) <= e/2 dot 1/6 = e/12 approx #{ calc.round(calc.e / 12, digits: 3) }
|
|
$
|
|
|
|
this gives the same bound as part c), confirming that both error bound formulas
|
|
are equivalent for the trapezoidal rule. the actual error
|
|
$(3-e)/2 approx #{ calc.round((3 - calc.e) / 2, digits: 3) }$
|
|
is indeed less than this bound.
|
|
|
|
== e)
|
|
|
|
to derive the composite trapezoidal rule, we divide the interval $[a,b]$ into $n$ segments of equal length $h = (b-a)/n$, with points $x_i = a + i h$ for $i = 0, 1, ..., n$.
|
|
|
|
applying the basic trapezoidal rule to each subinterval $[x_i, x_(i+1)]$:
|
|
$
|
|
T_([x_i, x_(i+1)])[f] = h/2 (f(x_i) + f(x_(i+1)))
|
|
$
|
|
|
|
the composite rule sums over all subintervals:
|
|
$
|
|
T_([a,b])^n [f] = sum_(i=0)^(n-1) T_([x_i, x_(i+1)])[f] = sum_(i=0)^(n-1) h/2 (f(x_i) + f(x_(i+1)))
|
|
$
|
|
|
|
expanding this sum:
|
|
$
|
|
T_([a,b])^n [f] = h/2 [f(x_0) + f(x_1) + f(x_1) + f(x_2) + ... + f(x_(n-1)) + f(x_n)]
|
|
$
|
|
|
|
we observe that $f(x_0)$ and $f(x_n)$ appear once, while $f(x_1), f(x_2), ..., f(x_(n-1))$ each appear twice. collecting terms:
|
|
$
|
|
T_([a,b])^n [f] = h/2 (f(x_0) + f(x_n)) + h sum_(i=1)^(n-1) f(x_i)
|
|
$
|
|
|
|
since $x_0 = a$ and $x_n = b$, this gives us the desired composite formula.
|
|
|
|
== f)
|
|
|
|
to derive the error bound for the composite trapezoidal rule, we
|
|
start with the error bound for each subinterval.
|
|
|
|
for each subinterval $[x_i, x_(i+1)]$ with length $h$:
|
|
$
|
|
|E_([x_i, x_(i+1)])[f]| <= h^3/12 max_(xi in [x_i, x_(i+1)]) |f''(xi)|
|
|
$
|
|
|
|
the total error for the composite rule is:
|
|
$
|
|
|E_([a,b])^n [f]| = |sum_(i=0)^(n-1) E_([x_i, x_(i+1)])[f]|
|
|
$
|
|
|
|
using the triangle inequality:
|
|
$
|
|
|sum_(i=0)^(n-1) E_([x_i, x_(i+1)])[f]| <= sum_(i=0)^(n-1) |E_([x_i, x_(i+1)])[f]|
|
|
$
|
|
|
|
substituting the error bound for each subinterval:
|
|
$
|
|
sum_(i=0)^(n-1) |E_([x_i, x_(i+1)])[f]| <= sum_(i=0)^(n-1) h^3/12 max_(xi in [x_i, x_(i+1)]) |f''(xi)|
|
|
$
|
|
|
|
since $[x_i, x_(i+1)] subset [a,b]$, we have:
|
|
$
|
|
max_(xi in [x_i, x_(i+1)]) |f''(xi)| <= max_(xi in [a,b]) |f''(xi)|
|
|
$
|
|
|
|
therefore:
|
|
$
|
|
sum_(i=0)^(n-1) h^3/12
|
|
max_(xi in [x_i, x_(i+1)])
|
|
|f''(xi)| & <= h^3/12
|
|
sum_(i=0)^(n-1) max_(xi in [a,b]) |f''(xi)| \
|
|
& = h^3/12 dot n dot max_(xi in [a,b]) |f''(xi)| \
|
|
& = h^2 dot (b-a)/12 max_(xi in [a,b]) |f''(xi)| \
|
|
& = (b-a)^3/(12n^2) max_(xi in [a,b]) |f''(xi)|
|
|
$
|
|
|
|
|
|
== g)
|
|
|
|
using the error bound formula:
|
|
$
|
|
|E_([a,b])^n [f]| <= (b-a)^3/(12n^2) max_(xi in [a,b]) |f''(xi)|
|
|
$
|
|
|
|
*first problem:* $|E_([0,1])^n [exp]| <= 10^(-3)$
|
|
|
|
for $f(x) = e^x$, we have $f''(x) = e^x$.
|
|
since $e^x$ is increasing, $max_(xi in [0,1]) |f''(xi)| = e$.
|
|
|
|
$
|
|
1^3/(12n^2) dot e <= 10^(-3)
|
|
$
|
|
|
|
solving: $n^2 >= e/(12 dot 10^(-3)) = 1000e/12$
|
|
|
|
therefore $n >= sqrt(1000e/12) approx 15.05$, so $n >= 16$.
|
|
|
|
*second problem:* $|E_([0,1])^n [ln]| <= 10^(-5)$
|
|
|
|
for $f(x) = ln(x)$, we have $f''(x) = -1/x^2$.
|
|
since $|f''(x)| = 1/x^2 -> infinity$ as $x -> 0^+$,
|
|
the function is not twice differentiable at $x = 0$.
|
|
the trapezoidal rule error bound does not apply.
|
|
|
|
*third problem:* $|E_([1,2])^n [e^x/x]| <= 10^(-3)$
|
|
|
|
for $f(x) = e^x/x$, using the quotient rule:
|
|
$f''(x) = (e^x (x^2 - 2x + 2))/x^3$
|
|
|
|
to find the maximum, we check:
|
|
- $f''(1) = e$
|
|
- $f''(1.5) approx 1.66$
|
|
- $f''(2) = e^2/4 approx 1.85$
|
|
|
|
since $f''(1) > f''(1.5)$ and $f''(1) > f''(2)$,
|
|
we have $max_(xi in [1,2]) |f''(xi)| = e$.
|
|
|
|
$
|
|
1^3/(12n^2) dot e <= 10^(-3)
|
|
$
|
|
|
|
solving: $n >= sqrt(1000e/12) approx 15.05$, so $n >= 16$.
|
|
|
|
|
|
== h) & i)
|
|
|
|
empirically testing shows that $n = 12$ subintervals already
|
|
achieves $|E_([0,1])^12 [exp]| <= 10^(-3)$. furthermore, we need $120$
|
|
subintervals for error $10^(-5)$.
|
|
|
|
|
|
== j)
|
|
|
|
to estimate the convergence rate, we fit the model $E_n approx C h^p$ by taking logarithms:
|
|
$
|
|
log |E_n| approx log C + p log h
|
|
$
|
|
|
|
using linear regression on $(log h, log |E_n|)$ data for various values of $n$:
|
|
|
|
for $exp(x)$ on $[0,1]$:
|
|
- estimated $p_1 approx 2.0$
|
|
|
|
for $sqrt(x)$ on $[0,1]$:
|
|
- estimated $p_2 approx 1.5$
|
|
|
|
|
|
== k)
|
|
|
|
for exp(x): $p_1 approx 2.0 =>$ quadratic
|
|
|
|
for sqrt(x): $p_2 approx 1.5$ shows that $f''(x) = -1/4 x^(-3/2)$ is unbounded as $x -> 0^+$.
|
|
|
|
|
|
== l)
|
|
|
|
given
|
|
$
|
|
E_([a,b])^1 = -M/12 (b-a)^3 quad "and" quad
|
|
E_([a,b])^2 = -M/48 (b-a)^3
|
|
$
|
|
|
|
we obtain $E_([a,b])^1 = 4 E_([a,b])^2$.
|
|
|
|
since
|
|
$
|
|
E_([a,b])^1 = I_([a,b])[f] - T_([a,b])^1 quad "and" quad
|
|
E_([a,b])^2 = I_([a,b])[f] - T_([a,b])^2
|
|
$
|
|
we get
|
|
$
|
|
E_([a,b])^1 - E_([a,b])^2 = T_([a,b])^2 - T_([a,b])^1
|
|
$
|
|
|
|
thus
|
|
|
|
$
|
|
4 E_([a,b])^2 - E_([a,b])^2 & = T_([a,b])^2 - T_([a,b])^1 \
|
|
3 E_([a,b])^2 & = T_([a,b])^2 - T_([a,b])^1 \
|
|
==> E_([a,b])^1 & = 4 E_([a,b])^2
|
|
= 4/3 (T_([a,b])^2 - T_([a,b])^1) =: cal(E)_([a,b])^1 [f] \
|
|
E_([a,b])^2 & = 1/3 (T_([a,b])^2 - T_([a,b])^1) =: cal(E)_([a,b])^2 [f]
|
|
$
|
|
|
|
|
|
== m)
|
|
|
|
for $f(x) = sqrt(x)$ on $[0,1]$ with 10 uniform intervals and tolerance $10^(-4)$:
|
|
#[
|
|
#let sqrt_trap1(a, b) = calc.round((b - a) / 2 * (calc.sqrt(a) + calc.sqrt(b)), digits: 6)
|
|
#let sqrt_trap2(a, b) = {
|
|
let c = (a + b) / 2
|
|
calc.round((b - a) / 4 * (calc.sqrt(a) + 2 * calc.sqrt(c) + calc.sqrt(b)), digits: 6)
|
|
}
|
|
#let error_est(a, b) = calc.round(calc.abs(4 / 3 * (sqrt_trap2(a, b) - sqrt_trap1(a, b))), digits: 6)
|
|
|
|
#show table.cell.where(y: 0): set text(weight: "bold")
|
|
#table(
|
|
columns: 4,
|
|
inset: 5% + 6pt,
|
|
[Interval], [$T_([a,b])^1$], [$T_([a,b])^2$], [$|cal(E)_([a,b])^1|$],
|
|
|
|
[$[0, 0.1]$], [#sqrt_trap1(0, 0.1)], [#sqrt_trap2(0, 0.1)], [#error_est(0, 0.1)],
|
|
[$[0.1, 0.2]$], [#sqrt_trap1(0.1, 0.2)], [#sqrt_trap2(0.1, 0.2)], [#error_est(0.1, 0.2)],
|
|
[$[0.2, 0.3]$], [#sqrt_trap1(0.2, 0.3)], [#sqrt_trap2(0.2, 0.3)], [#error_est(0.2, 0.3)],
|
|
[$[0.3, 0.4]$], [#sqrt_trap1(0.3, 0.4)], [#sqrt_trap2(0.3, 0.4)], [#error_est(0.3, 0.4)],
|
|
[$[0.4, 0.5]$], [#sqrt_trap1(0.4, 0.5)], [#sqrt_trap2(0.4, 0.5)], [#error_est(0.4, 0.5)],
|
|
[$[0.5, 0.6]$], [#sqrt_trap1(0.5, 0.6)], [#sqrt_trap2(0.5, 0.6)], [#error_est(0.5, 0.6)],
|
|
[$[0.6, 0.7]$], [#sqrt_trap1(0.6, 0.7)], [#sqrt_trap2(0.6, 0.7)], [#error_est(0.6, 0.7)],
|
|
[$[0.7, 0.8]$], [#sqrt_trap1(0.7, 0.8)], [#sqrt_trap2(0.7, 0.8)], [#error_est(0.7, 0.8)],
|
|
[$[0.8, 0.9]$], [#sqrt_trap1(0.8, 0.9)], [#sqrt_trap2(0.8, 0.9)], [#error_est(0.8, 0.9)],
|
|
[$[0.9, 1.0]$], [#sqrt_trap1(0.9, 1.0)], [#sqrt_trap2(0.9, 1.0)], [#error_est(0.9, 1.0)],
|
|
)
|
|
]
|
|
|
|
intervals needing refinement (error > $10^(-4)$): first 4 intervals
|
|
|
|
|
|
== n)
|
|
|
|
implemented adaptive trapezoidal quadrature with tolerance $"tol" = 10^(-5)$:
|
|
|
|
#table(
|
|
columns: 4,
|
|
[Function], [Adaptive Result], [Intervals Used], [Uniform (part i)],
|
|
|
|
[$e^x$], [1.7183], [64], [120],
|
|
[$sqrt(x)$], [0.6667], [79], [-],
|
|
)
|
|
|
|
the adaptive algorithm automatically allocates computational effort where
|
|
needed, making it more efficient than uniform refinement.
|
|
|
|
|
|
|
|
|
|
|
|
= problem 2
|
|
|
|
== a)
|
|
|
|
to transform the quadrature rule from $[-1,1]$ to $(-3,3)$ we use
|
|
$
|
|
x = (b - a)/2 t + (b + a)/2 = 3 t
|
|
$
|
|
where $t in [-1,1]$, so $d x = 3 d t$. then
|
|
|
|
$
|
|
integral_(-3)^3 e^x d x = integral_(-1)^1 e^(3t) dot 3 d t = 3 integral_(-1)^1 e^(3t) d t
|
|
$
|
|
|
|
the gauß-legendre approximation is given by
|
|
$
|
|
G_h = 3 sum_(i=0)^3 w_i e^(3 x_i)
|
|
$
|
|
|
|
so
|
|
|
|
#[
|
|
#let x0 = -1 / 35 * calc.sqrt(525 + 70 * calc.sqrt(30))
|
|
#let x1 = -1 / 35 * calc.sqrt(525 - 70 * calc.sqrt(30))
|
|
#let x2 = 1 / 35 * calc.sqrt(525 - 70 * calc.sqrt(30))
|
|
#let x3 = 1 / 35 * calc.sqrt(525 + 70 * calc.sqrt(30))
|
|
|
|
#let w0 = 1 / 36 * (18 - calc.sqrt(30))
|
|
#let w1 = 1 / 36 * (18 + calc.sqrt(30))
|
|
#let w2 = 1 / 36 * (18 + calc.sqrt(30))
|
|
#let w3 = 1 / 36 * (18 - calc.sqrt(30))
|
|
|
|
$
|
|
G_h & = 3(w_0 e^(3 x_0) + w_1 e^(3 x_1) + w_2 e^(3 x_2) + w_3 e^(3 x_3)) \
|
|
& = 3(#(w0 * calc.exp(3 * x0) + w1 * calc.exp(3 * x1) + w2 * calc.exp(3 * x2) + w3 * calc.exp(3 * x3))) \
|
|
& = #(3 * (w0 * calc.exp(3 * x0) + w1 * calc.exp(3 * x1) + w2 * calc.exp(3 * x2) + w3 * calc.exp(3 * x3)))
|
|
$
|
|
|
|
wheras the exact value is $e^3 - e^(-3) = #(calc.exp(3) - calc.exp(-3))$.
|
|
]
|
|
|
|
|
|
|
|
== b)
|
|
|
|
from the given error formula, each subinterval of length $h = (b-a)/m$ contributes error
|
|
|
|
$
|
|
E_i prop h^(2n+1) = ((b-a)/m)^(2n+1)
|
|
$
|
|
|
|
total error from m subintervals
|
|
$
|
|
E_m prop m dot ((b-a)/m)^(2n+1) = (b-a)^(2n+1)/m^(2n)
|
|
$
|
|
|
|
we expect error to decrease when we subdivide the interval. more subintervals
|
|
means more precision.
|
|
|
|
|
|
|
|
== c)
|
|
|
|
let $f(x) = x^8/8!$
|
|
|
|
we can see that $f^((2n))(x) = x^(8-2n)/(8-2n)!$ for $n > 0$ and $<= 8$ so
|
|
$
|
|
E_1 = (6^(2n+1) (n!)^4)/((2n+1)[(2n)!]^3) dot xi^(8-2n)/((8-2n)!)
|
|
$
|
|
|
|
then sum the left and right errors to obtain
|
|
$
|
|
E_2 = (3^(2n+1) (n!)^4)/((2n+1)[(2n)!]^3) dot (xi_1^(8-2n) + xi_2^(8-2n))/((8-2n)!)
|
|
$
|
|
|
|
which yields
|
|
$
|
|
E_2/E_1 = (xi_1^(8-2n) + xi_2^(8-2n))/(2^(2n+1) xi^(8-2n))
|
|
$
|
|
|
|
we can see that increasing n makes it more accurate very fast, as the numerator
|
|
grows smaller exponentially and the denominator bigger.
|
|
|
|
|
|
|
|
|
|
= problem 3
|
|
|
|
== a)
|
|
|
|
to find an orthogonal basis $p_j (x)$ from the given canonical basis
|
|
$
|
|
#let phi = $phi.alt$
|
|
phi_0(x) equiv 1, med
|
|
phi_1(x) = x, med
|
|
phi_2(x) = x^2, med
|
|
phi_3(x) = x^3,
|
|
$
|
|
we can use the gram-schmidt process.
|
|
|
|
recall that
|
|
$
|
|
p_j (x) = phi_j (x) - sum_(k=0)^(j-1) {(integral_0^1 phi_j (x) p_k (x) dd(x))
|
|
/(integral_0^1 [p_k (x)]^2 dd(x))} p_k (x),
|
|
$
|
|
and notice that $p_0(x) equiv 1$.
|
|
|
|
thus we obtain
|
|
$
|
|
p_1(x) & = phi_1(x) - (integral_0^1 phi_1(x) p_0(x) dd(x))
|
|
/ (integral_0^1 [p_0(x)]^2 dd(x)) p_0 (x) \
|
|
& = x - (integral_0^1 x dd(x))/(integral_0^1 dd(x))
|
|
= x - 1/2, \
|
|
p_2(x) & = phi_2(x) - (integral_0^1 phi_2(x) p_1(x) dd(x))
|
|
/ (integral_0^1 [p_1(x)]^2 dd(x)) p_1 (x) \
|
|
& - (integral_0^1 phi_2(x) p_0(x) dd(x))
|
|
/ (integral_0^1 [p_0(x)]^2 dd(x)) p_0 (x) \
|
|
& = x^2 - (integral_0^1 x^2 (x - 1/2) dd(x))
|
|
/ (integral_0^1 [x-1/2]^2 dd(x)) (x - 1/2) \
|
|
& - (integral_0^1 x^2 dd(x)) / (integral_0^1 dd(x))
|
|
= x^2 - x + 1/6 \
|
|
p_3(x) & = limits(dots)^#[magic] = x^3 - 3/2 x^2 + 3/5 x - 1/20
|
|
$
|
|
|
|
== b)
|
|
|
|
we perform polynomial division to find the remaining two roots and obtain
|
|
$
|
|
p_3 (x) : (x - 1/2) = x^2 - x + 1/10 \
|
|
=> x = 1 / 2 plus.minus 1/2 sqrt(3/5)
|
|
$
|
|
|
|
so
|
|
$
|
|
x_0 & = 1/2 - 1/2 sqrt(3/5), \
|
|
x_1 & = 1/2, \
|
|
x_2 & = 1/2 + 1/2 sqrt(3/5).
|
|
$
|
|
|
|
|
|
|
|
== c)
|
|
|
|
let
|
|
$
|
|
w_i = integral_0^1 ell_i (x) dd(x)
|
|
$
|
|
with
|
|
$
|
|
ell_k (x) = product_(i = 0\ i != k)^(n - 1) (x - x_i) / (x_k - x_i)
|
|
$
|
|
where $n - 1 = 2$ since $k = 0, 1, 2$.
|
|
|
|
computing the cardinal polynomials, we obtain
|
|
$
|
|
ell_0(x) & = (x - x_1) / (x_0 - x_1) dot (x - x_2) / (x_0 - x_2) \
|
|
& = (x - 1/2) / (-1/2 sqrt(3/5)) dot (x - 1/2 - 1/2 sqrt(3/5)) / (-sqrt(3/5)) \
|
|
& = 5/3 (-2x^2 + sqrt(3/5) x + 1/2 sqrt(3/5) - 1/2) \
|
|
ell_1(x) & = (x - x_0) / (x_1 - x_0) dot (x - x_2) / (x_1 - x_2) \
|
|
& = (x - 1/2 + 1/2 sqrt(3/5)) / (1/2 sqrt(3/5)) dot (x - 1/2 - 1/2 sqrt(3/5)) / (1/2 sqrt(3/5)) \
|
|
& = 4/3 (x^2 - x + 1/10) \
|
|
ell_2(x) & = (x - x_0) / (x_2 - x_0) dot (x - x_1) / (x_2 - x_1) \
|
|
& = 5/3 (-2x^2 - sqrt(3/5) x + 1/2 sqrt(3/5) + 1/2)
|
|
$
|
|
|
|
integrating to find the weights
|
|
$
|
|
w_0 & = integral_0^1 ell_0(x) dd(x) = 5/18 \
|
|
w_1 & = integral_0^1 ell_1(x) dd(x) = 4/9 \
|
|
w_2 & = integral_0^1 ell_2(x) dd(x) = 5/18
|
|
$
|
|
|
|
by symmetry we have $w_0 = w_2$.
|
|
|
|
|
|
== d)
|
|
|
|
the 3-point gauß-legendre quadrature rule for $[0,1]$ is
|
|
$
|
|
integral_0^1 f(x) dd(x) approx 5/18 f(x_0) + 4/9 f(1/2) + 5/18 f(x_2)
|
|
$
|
|
where $x_0 = 1/2 - 1/2 sqrt(3/5)$ and $x_2 = 1/2 + 1/2 sqrt(3/5)$.
|
|
|
|
to verify that the 3-point rule integrates $x^5$ exactly, we compute
|
|
$
|
|
integral_0^1 x^5 dd(x) = x^6/6 |_0^1 = 1/6
|
|
$
|
|
and using the symmetry property and binomial expansion for $x_0^5 + x_2^5$ where
|
|
$x_0 = 1/2 - 1/2 sqrt(3/5)$ and $x_2 = 1/2 + 1/2 sqrt(3/5)$, we obtain
|
|
|
|
$
|
|
5/18 (x_0^5 + x_2^5) + 4/9 dot 1/32 = 1/6
|
|
$
|
|
|
|
#h(1fr) $qed$
|
|
|
|
|
|
== e)
|
|
|
|
to show equivalence with the standard table, we transform from the reference interval $[-1, 1]$ to $[0, 1]$ using:
|
|
$
|
|
x = (t + 1)/2 quad "where" t in [-1, 1] "and" x in [0, 1]
|
|
$
|
|
|
|
so
|
|
$
|
|
x_0 & = (-sqrt(15)/5 + 1)/2 = 1/2 - sqrt(15)/10 \
|
|
x_1 & = (0 + 1)/2 = 1/2 \
|
|
x_2 & = (sqrt(15)/5 + 1)/2 = 1/2 + sqrt(15)/10
|
|
$
|
|
|
|
using the change of variables $x = (t+1)/2$, we have $dd(x) = 1/2 dd(t)$.
|
|
so
|
|
$
|
|
w_0 & = 5/9 dot 1/2 = 5/18 \
|
|
w_1 & = 8/9 dot 1/2 = 4/9 \
|
|
w_2 & = 5/9 dot 1/2 = 5/18
|
|
$
|
|
|
|
note that $sqrt(15)/10 = 1/2 sqrt(3/5)$ so the transformed rule matches exactly
|
|
our derived 3-point gauß-legendre quadrature for $[0, 1]$, confirming that both
|
|
approaches yield the same integration rule.
|
|
|
|
|