#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.