HW SIX

1. 6.2 Exact value of the integral is given by
   $$I = \int\limits_1^1 e^{-2 x } d x = -\frac{e^{-2}-e^2}{2} = 3.626860407847019$$

   (a)

   Trapezoid

   $$T = \frac{1}{4}\left(\frac{e^{2}}{2} + e+ 1 + e^{-1}+\frac{e^{-2}}{2} \right) = 3.924178480357059,$$
   or, in Matlab,

   (1/2)*(exp(2)/2+exp(1)+ 1+exp(-1)+exp(-2)/2)
   

   The absolute relative error is $8.1976706869326$ percent

   (b) Simpson's method:

   $$S = \frac{1}{6}\left( e^{2}+ 4 e^{1}+ 2 + 4*e^{-1} +e^{-2}\right)= 3.633561469845151,$$
   or, in Matlab,

   (1/6)*(exp(2)+4*exp(1)+ 2+4*exp(-1)+exp(-2)/2)
   

   with the relative error of $7.40580511224$ percent.

   (c) Hermit's rule, or corrected trapezoid rule gives

   $$H = \frac{1}{4}\left(\frac{e^{2}}{2} + e+ 1 + e^{-1}+\frac{e^{-2}... ...\right) + \frac{1}{4\cdot 12}\left(-2e^{2}+2e^{-2}\right) = 3.621940113036475 ,$$
   or, in Matlab,

   (1/2)*(exp(2)/2+exp(1)+ 1+exp(-1)+exp(-2)/2)+1/(4*12)*(-2*exp(2)+2*exp(-2))
   

   with the absolute relative error of $0.1356626464$ percent.

   (d)

   Since $f'(x) = - 2 e^{-2 x}$, and $f''(x) = 4 e^{-2 x}$ on an interval from $-1$ to $1$ we have $\vert\vert f''(x)\vert\vert _\infty=4 e^2 = 29.55622439572260$, we obtain

   $$h\simeq\sqrt{\frac{6\cdot 10^{-6}}{\vert\vert f''(x)\vert\vert _\infty}}\simeq 4 4.505584\cdot 10^{-4}.$$

   (c) Since $f''''(x) = 16 e^{-2 x}$ we have, for this problem,

   $$\vert\vert f''(x)\vert\vert _\infty\simeq 118.2248,$$
   so that
   $$h\simeq \left(\frac{10^{-6}\cdot 180}{2\cdot 118.22}\right)^{\frac{1}{4}}\simeq 0.029538156.$$

2. 6.10

   (a)

   This problem was solved in class, see the lecture notes.

   There is a misprint in the book. The correct formulation of this problem is

   $$S = \frac{1}{3}M(n)+\frac{2}{3}T(n).$$

   Suppose $n=1$. Then the trapezoid method will give for $a$,$b$ and $2 m = a+ b$ with $h=(b-a)$:

   $$T = h (\frac{1}{2}f(a)+ \frac{1}{2}f(b)).$$
   Then the midpoint method for $[a,b]$ will give
   $$M = h f(m).$$
   Then
   $$\frac{2}{3}M+ \frac{1}{3}T = h (\frac{1}{6}f(a)+ \frac{2}{3} f(m)+ +\frac{1}{6}f(b),$$
   which is the Simpson's method.

   Results for arbitrary $n$ are obtained by trivial generalization.

   (b) There is a misprint in this problem as well, it should be $\frac{4}{3}T(n) - \frac{2}{3}T(n/2).$

   Now consider $n=2$ with the two intervals $[a,m]$ and $[m,b]$ with $2 m = a+b.$ We now have composite trapezoid method $T_2$ for two intervals to be given by

   $$T_2 = h (\frac{1}{2}f(a)+ f(m)+\frac{1}{2} f(b)),$$
   and the trapezoid method for the interval $[a,b]$ given by
   $$T_1 = h (\frac{1}{2}f(a)+ \frac{1}{2} f(b)).$$
   Then
   $$\frac{4}{3}T_2-\frac{1}{3}T_1 = h (\frac{1}{3}f(a)+ \frac{4}{3}f(m)+\frac{1}{3} f(b)),$$
   which is the Simpson's method.

   Results for arbitrary $n$ are obtained by trivial generalization.

3. 6.14

   (a)

   Assume $h = (b-a)/n$. Since

   $$\int\limits_a^b f(x) d x - I_s(n) =\alpha h^4 + \beta h^6,$$ (48)

   it follows that

   $$\int\limits_a^b f(x) d x - I_s(2 n) =\alpha \frac{h^4}{2^4} + \beta \frac{h^6}{2^6},$$ (49)

   because for $2n$ the value of $h$ is twice smaller.

   Then

   $$\int\limits_a^b f(x) d x - \frac{1}{15}(16 I_s(2 n) - I_s(n))= -\frac{\beta}{30} h^6.$$

   (b) Similarly, if

   $$\int\limits_a^b f(x) d x - I(n) =\alpha h^2 + \beta h^3+\gamma h^4,$$
   then
   $$\int\limits_a^b f(x) d x - I(2n) =\alpha \frac{h^2}{4} + \beta \frac{h^3}{8}+\gamma \frac{h^4}{16},$$
   and
   $$\int\limits_a^b f(x) d x - I(4n) =\alpha \frac{h^2}{16} + \beta \frac{h^3}{64}+\gamma \frac{h^4}{256}.$$
   Consequently,
   $$\int\limits_a^b f(x) d x -\frac{1}{21}(32 I(4n) - 12 I (2n)+ I (n)) = \gamma \frac{3}{8}h^4.$$

   (c)

   The above equations are now supplemented with

   $$\int\limits_a^b f(x) d x - I(3n) =\alpha \frac{h^2}{9} + \beta \frac{h^3}{7}+\gamma \frac{h^4}{81}.$$
   We now obtain
   $$\int\limits_a^b f(x) d x -\frac{1}{12}(27 I(3n) - 16 I (2n)+ I (n)) = \gamma \frac{1}{4}h^4.$$

4. 6.16 (a,b,e)

   (a) The value of $erf(2)$ obtained by Matlab is given by $0.995322265018952$. Interestingly, the “correct” result is
   

   $$0.99532226501895273416206925636725292861089179704006007673835232620043728071$$

   $$\hskip 2cm \dots\dots 99951773676290080196806805.$$

   
   

   To use trapezoid rule, we need to calculate $f''(x).$ We have

   $$f(x)=\frac{2}{\sqrt{\pi}}e^{-x^2},$$
   which gives
   $$f''(x) = \frac{8 x^2 - 4}{\sqrt{\pi}e^{x^2}}.$$
   The maximum value of $\vert f''(x)\vert$ is achieved at $x = 0$ and it is equal to $4/\sqrt{\pi}\simeq 2.26$. Therefore to achieve accuracy of six digits, i.e. $10^{-6}$ we have
   $$h = \sqrt{ \frac{12\cdot 10^{-6}}{2\cdot 2.26}} \simeq 0.001629376.$$

   (b)We need to know $f''''(x)$, which is equal to

   $$f''''(x) = \frac{(8*(3 + 4*x^2*(-3 + x^2)))}{e^{x^2}*\sqrt{\pi}}.$$
   Its maximum value is again achieved at $x = 0$, so that
   $$f''''(x=0)\simeq 13.54.$$
   We now can calculate $h$ for the composite Simpson's method to be equal to
   $$h=\left(\frac{180\cdot 10^{-6}}{2\cdot 13.54}\right)^{\frac{1}{4}}\simeq 0.0507757.$$

   (e) The following code in Mathematica plots the required curve:

   Plot[Erf[x], {x, 0, 10}, PlotRange -> All]
   

   To obtain this in Matlab, one can use composite Simson's method, as explained above, with $h\simeq 0.05.$

5. 6.21

   (a) Since the problem is symmetric with respect to $x\to -x$ and $y\to -y$, the curves are symmetric with respect to mirrowing over both $x$ and $y$ axis, so therefore it is sufficient to consider only positive $x$ and $y$ and multiply the result by $4$.

   (b)

   To evaluate the area under the curve

   $$x^4+2 y^4 =1,$$
   we need to evaluate numerically the integral
   $$I = \int\limits_{-1}^1 \Big(\frac{1-x^4}{2}\Big)^{\frac{1}{4}} d x.$$

   We are therefore integrating

   $$f(x) = \Big(\frac{1-x^4}{2}\Big)^{\frac{1}{4}} .$$

   We now obtain

   $$f''(x) = \frac{-3 x^2}{2^{\frac{1}{4}}(1-x^4)^{\frac{7}{4}}},$$
   and
   $$f''''(x) = \frac{-3 (2+45 x^4 +30 x^8} {2^{\frac{1}{4}}(1-x^4)^{\frac{15}{4}}}.$$
   Since we are integrating from $-1$ to $1$, both of these derivatives diverge at $x=1$. This is the reflection of the fact that the integrand is singular at $x=1$. Consequently, one can not obtain any useful error bound on these integrals.

   Solving for $x$ and integrating over $y$ will lead to the same problems. This can be observed by changing variables $x =2^{\frac{1}{4}} x'$ and $y=2^{-\frac{1}{4}} y'$.

   (c)

   After using adaptive quadrature we otain the result

   $$I\simeq 0.779542374877720.$$

   (d)

   We now calculate

   $$I_1 = \int\limits_0^{\frac{1}{2}} \Big(\frac{1-x^4}{2}\Big)^{\frac{1}{4}} d x \simeq 0.419116752331519,$$
   and
   $$I_2 = \int\limits_0^{(\frac{15}{32})^{\frac{1}{4}}} \Big((1-2 y^4)^{\frac{1}{4}} d y \simeq 0.774144487502060.$$

   To compare this result with the (c) we need to subtract $(15/32)^{(1/4)}/2$ from it

   $$\tilde{I} = 0.419116752331519 + 0.77414448750206 - (15/32)^{(1/4)}/2 = 0.779542374877720.$$
   Remarkably, this coincides with (c).

   (e) In the point (d) the issue of infinite second derivative was avoided by reformulating a problem in a way that no infinite second derivative was present. This is a result of cutting the domain in two, as described in the Holmes 2023.

6. 6.29

   (a)

   We studied in class how to construct a system of equations to obain the solution to this problem. We require that this quadrature rule integrates correctly $f(x)=1$, $f(x)=x$ and $f(x)=x^2$. The resulting equations are
   

   $$b - a = \omega_1+\omega_2,$$

   $$b^2/2 - a^2/2= \omega_1 a + \omega_2 z,$$

   $$b^3/3 - a^3 /3=\omega_1 a^2 +\omega_2 z^2.$$

   
   Its solutions is given by
   $$\omega_1 = \frac{b-a}{4},\omega_2 = \frac{3(b-a)}{4},z=\frac{a+ 2 b}{3}.$$

   (b)

   To find $K$ we need to integrate the $f(x)=x^3$ and calculate the error. We therefore calculate

   $$I = \int\limits_a^b x^3 = \frac{b^4}{4}-\frac{a^4}{4},$$
   and
   $$E = \left( \frac{b^4}{4}-\frac{a^4}{4}\right)- \left(\frac{b-a}{... ...+ \frac{3(b-a)}{4}\left(\frac{(a+2 b)}{3}\right)^3=\frac{(b-a)^4}{36} \right).$$

   Since $f'''(x)=6$ we obtain that the constant $K$ is equal to $K = ( 6 \cdot 36)^{-1}= (216)^{-1}.$

7. The error terms in teh Simpson's method were derived in class. To derive results for the Gauss method, please use the method which was outlined in the lectures (see lecture notes).