»» Scientific Computing

Scientific Computing

The following section presupposes Set Theory and Nonstandard Analysis.

Definition: Let $$f_n^*(z) = f(\eta_nz)$$ sisters of the TS $$f(z) \in \mathcal{O}(D)$$ centred on 0 on the domain $$D \subseteq {}^{\omega}\mathbb{C}$$ where $$m, n \in {}^{\omega}\mathbb{N}^{*}$$ and $$\eta_n^m := i^{2^{\lceil m/n \rceil}}$$. Then let $$\delta_n^*f = \tilde{2}(f – f_n^*)$$ the halved sister distances of $$f$$. Let $$\mu_n^m := m!n!/(m + n)!$$ and $$\widetilde{(-n)!} := 0$$. For $$j, k \in {}^{\omega}\mathbb{N}, f \in \mathcal{C}_{\pi}^{j+2}$$ and Fourier coefficients $$c_k := \tilde{\hat{\pi}}{\uparrow}_{-\pi}^{\pi}{f(t)\tilde{e}^{ikt}{\downarrow t}}$$1cf. Walter, Wolfgang: Analysis 2; 5., erw. Aufl.; 2002; Springer; Berlin, p. 358 ff., the series $${+}_{k=-\omega}^{\omega}{c_ke^{ikt}}$$ is called Fourier series.$$\triangle$$

Remark: A conversion to other period lengths than $$\hat{\pi }$$ is easily possible2see loc. cit., p. 364. On the level of TS3cf. Remmert, Reinhold: Funktionentheorie 1; 3., verb. Aufl.; 1992; Springer; Berlin, p. 165 f., $$\mu$$ and $$\eta$$ form a $$\mu$$-$$\eta$$-calculus, which allows an easy and finite closed representation of integrals and derivatives.

Representation theorem for integrals: The TS (see below) $$f(z) \in \mathcal{O}(D)$$ centred on 0 on $$D \subseteq {}^{\omega}\mathbb{C}$$ gives for $$\grave{m}, n \in {}^{\omega}\mathbb{N}^*$$${\uparrow}_0^z…{\uparrow}_0^{\zeta_2}{f(\zeta_1){\downarrow}\zeta_1\;…\;{\downarrow}\zeta_n} = \widetilde{n!} f(z\mu_n) z^n.\square$Example: For the TS $$f(x), g(x) \in {}^{\omega}\mathbb{R}$$, it holds that${\uparrow}_0^x{f(v){\downarrow}v}{\uparrow}_0^x{\uparrow}_0^{y}{g(v){\downarrow}v{\downarrow}y} = \tilde{2}f(x\mu_1)g(x\mu_2)x^3.$Representation theorem for derivatives: For $$\mathbb{B}_{\tilde{\nu}}(0) \subset D \subseteq {}^{\omega}\mathbb{C}, n$$-th unit roots, TS$f(z):=f(0) + {+}_{m=1}^{\omega }{\widetilde{m!}\,{{f}^{(m)}}(0){z^m}},$$$b_n := \tilde{\varepsilon}^{n}\,\acute{n}! = 2^j, j, n \in {}^{\omega}\mathbb{N}^{*}, \varepsilon \in ]0, r[, u := e^{\tilde{n}\hat{\pi}i}$$ and $$f$$’s radius of convergence $$r \in {}^{\nu}{\mathbb{R}}_{>0}$$ imply${{f}^{(n)}}(0)=b_n{+}_{k=1}^{n}{\delta_n^* f(\varepsilon u^k)}.$Fundamental theorem of algebra: Every non-constant polynomial $$p \in {}^{(\omega)}\mathbb{C}$$ has at least one complex root.

Indirect proof: By performing an affine substitution of variables, reduce to the case $$\widetilde{p(0)} \ne \mathcal{O}(\iota)$$. Suppose that $$p(z) \ne 0$$ for all $$z \in {}^{(\omega)}\mathbb{C}$$. Since $$f(z) := \widetilde{p(z)}$$ is holomorphic, it holds that $$f(\tilde{\iota}) = \mathcal{O}(\iota)$$. By the mean value inequality $$|f(0)| \le {|f|}_{\gamma}$$4see a.a.O., p. 160 for $$\gamma = \partial\mathbb{B}_{r}(0)$$ and arbitrary $$r \in {}^{(\omega)}\mathbb{R}_{>0}$$, and hence $$f(0) = \mathcal{O}(\iota)$$, which is a contradiction.$$\square$$

Universal multistep theorem: For $$n \in {}^{\nu}\mathbb{N}_{\le p}, k, m, p \in {}^{\nu}\mathbb{N}^{*}, d_{\curvearrowright B} x \in\, ]0, 1[, x \in [a, b] \subseteq {}^{\omega}\mathbb{R}, y : [a, b] \rightarrow {}^{\omega}\mathbb{R}^q, f : [a, b]\times{}^{\omega}\mathbb{R}^{q \times n} \rightarrow {}^{\omega}\mathbb{R}^q, g_k(\curvearrowright B x) := g_{\acute{k}}(x)$$, and $$g_0(a) = f((\curvearrowleft B)a, y_0, … , y_{\acute{n}})$$, the TS of the initial value problem $$y^\prime(x) = f(x, y((\curvearrowright B)^0 x), … , y((\curvearrowright B)^{\acute{n}} x))$$ of order $$n$$ implies$y(\curvearrowright B x) = y(x) + {\downarrow}_{\curvearrowright B}x{\pm}_{k=1}^{p}{\left (g_{p-k}((\curvearrowright B) x){+}_{m=k}^{p}{\widetilde{m!}\tbinom{\acute{m}}{\acute{k}}}\right )} + \mathcal{O}(({\downarrow}_{\curvearrowright B} x)^{\grave{p}}).\square$Theorem for (anti-) derivatives of TS: For $$j \in {}^{\omega}\mathbb{Z}$$, $$q = \tilde{\varepsilon}(z-a)$$, $$a \in D$$ and $$k, m \in \mathbb{N}_{<n}$$, modular arithmetic5cf. Knuth, Donald Ervin: The Art of Computer Programming Volume 2; 3rd Ed.; 1997; Addison Wesley; Reading, p. 302 – 311 and $$n$$-th unit roots result in the corresponding DFT form:${\updownarrow}^jf_n(z) := \tilde{n}(q^k)^T(\delta_{km}\widetilde{\varepsilon q}^j\widetilde{(k-j)!}k!)({\tilde{u}}^{km})(f({\varepsilon u}^m+a))+\mathcal{O}(\varepsilon^n).\square$Remark: The identity instead of $$\delta_n^*$$ provides arbitrarily precise approximations for the $$f^{(n)}$$. The last theorems are equally valid for multidimensional TS (with several sums) and Laurent series. The error is $$\mathcal{O}(\varepsilon^n)$$ instead of $$\mathcal{O}(\varepsilon)$$ for analogously defined $$m$$-dimensional DFT forms of the TS with $$\binom{m+n}{n}$$ derivatives where the effort is comparable. A good choice is $$\varepsilon = \tilde{2}$$ and $$n = 64$$ resp. $$n = 16$$ in the next conclusion.

Conclusion: DFT-zero methods iterate zeros $$a \in {}^{\omega}\mathbb{C}$$ for every function $$f(z) \in {}^{\omega}\mathbb{C}$$ that can be developed into a TS at defaults $$z_0 \in {}^{\omega}\mathbb{C}$$ with each time the same convergence as in methods similar to Simpson’s, if $${\updownarrow}^0f_n(z)$$ is differentiated (several times) for sufficiently precise $$|z_{\grave{m}} – z_m|$$ and $$m \in {}^{\omega}\mathbb{N}.\square$$

Theorem for derivatives of Fourier series: For $$f \in \mathcal{C}_{\pi}^{j+2}$$, $$j \in {}^{\omega}\mathbb{N}, k \in [-\acute{n}, n] \cap \mathbb{Z}, t \in [-\pi, \pi]$$ and $$m \in \mathbb{N}_{<\hat{n}}$$, $$\hat{n}$$-th unit roots result in the following DFT form:${\downarrow}^jf_n(t):=(u^{\tilde{\pi}knt})^T(\delta_{(k+\acute{n})m}(ik)^j)(\tilde{u}^{km})(f(\pi m/n – \pi))/{\hat{n}}+\mathcal{O}(\tilde{n}).\square$Conclusion: Supporting points $$mr := \pi m/n$$ of the smooth $$f(mr)$$ yield for $$k$$ like $$m$$ interpolating in $$\mathcal{O}(_en n)$$:${\downarrow}^jf_n(t):=(u^{\tilde{r}kt})^T(\delta_{km}(ik)^j)(\tilde{u}^{km})(f(mr))/{\hat{n}}.\square$Theorem (DFT fixed-point method): Every $$z \in {}^{\omega}\mathbb{C}$$ of an arbitrary $$m$$-polynomial $$p(z) = 0$$ with $$m \in [2, n] \cap \mathbb{N}$$ for $$n := 2^r, r \in {}^{\nu}\mathbb{N}^*$$ and coefficients from $${}^{\nu}\mathbb{C}$$ can be (initiating) determined also in $$\mathcal{O}(_en n)$$.

Proof and algorithm: Let $$U = (\tilde{u}^{jk})$$ for $$j, k \in \mathbb{N}_{<n}, u :=e^{\tilde{n} \hat{\pi} i}, q := 2z$$ and $$s_k := p(\tilde{2}u^k)$$. A simple transform achieves $$|q| < \tilde{2}$$ for all zeros $$\zeta$$ of $$p(z)$$ and $$p(0) = 1$$. It follows from $$p(z) = \tilde{n}(q^j)^TUs = \tilde{n}\mu^Ts = 0$$ the simplified iteration $$\mu^* = U_{1}^{-T}\mu U((\delta_{jk}\tilde{u}^j)U^{-1}\mu-(U_{\acute{n}}^{-T}\mu+\beta s^T\mu, \beta s^T\mu, …, \beta s^T\mu)^T)$$ with Kronecker delta $$\delta_{jk}$$, starting point $$q := \tilde{2}$$ and $$\beta \in {}^{\nu}\mathbb{C}^*$$ such that each time $$||\mu^*-\mu||$$ is roughly halved and $$\mu^Ts = 0$$ holds. Finish by polynomial division where $$m > 2.\square$$

Theorem (series product): For $${a}_{m}, {b}_{n} \in {}^{(\omega)}\mathbb{K}$$, the next equation replaces the Cauchy product6see Walter, Wolfgang: Analysis 1; 3., verb. Aufl.; 1992; Springer; Berlin, p. 103:${+}_{m=1}^{\omega }{{{a}_{m}}}{+}_{n=1}^{\omega }{{{b}_{n}}}={+}_{m=1}^{\omega }{\left( {+}_{n=1}^{m}{\left( {{a}_{n}}{{b}_{m-\acute{n}}}+{{a}_{\omega -\acute{n}}}{{b}_{\omega -m+n}} \right)}-{{a}_{m}}{{b}_{\omega -\acute{m}}} \right)}.\square$Example: The following series product has the finite value7cf. Gelbaum, Bernard R.; Olmsted, John M. H.: Counterexamples in Analysis; Republ., unabr., slightly corr.; 2003; Dover Publications; Mineola, New York, p. 61 f.:$\left({\pm}_{m=1}^{\mathrm{\omega}}{{\tilde{m}}^{\tilde{2}}}\right)^2={+}_{m=1}^{\mathrm{\omega}}{\left(\left(\frac{\tilde{m}}{\mathrm{\omega}-\acute{m}}\right)^{\tilde{2}}-{i^{\hat{m}}}{+}_{n=1}^{m}\left(\left(\frac{\tilde{n}}{m-\acute{n}}\right)^{\tilde{2}}+\left(\frac{\widetilde{\mathrm{\omega}-\acute{n}}}{\mathrm{\omega}-m\ \mathrm{+\ }n}\right)^{\tilde{2}}\right)\right)}=0,36590…\ \ \ \ll\frac{{\zeta\left(\tilde{2}\right)}^2}{3+2\;2^{\tilde{2}}}.$Example: The signum function sgn yields the following series product8cf. loc. cit., p. 62: ${+}_{m=0}^{\omega }{{2}^{{{m}^{\text{sgn}(m)}}}}{+}_{n=0}^{\omega}{\text{sgn}(n-\gamma)} = \acute{\omega}{2}^{\grave{\omega}}\gg -2.$Remarks: If the moduli of $$x \in \mathbb{C}$$, $${\downarrow}x$$ or $$\widetilde{{\downarrow}x}$$ have different orders of magnitude, the identity${{s}^{(0)}}(x):={\pm}_{m=0}^{n}{x^m}=(1-{{(-x)}^{\grave{n}}})/\grave{x}$yields by differentiating${{s}^{(1)}}(x)={\mp}_{m=1}^{n}{m{x^{\acute{m}}}}=(\grave{n}{{(-x)}^{n}}-n{{(-x)}^{\grave{n}}}-1)/{{{\grave{x}}^{2}}}.$The formulas above were sometimes miscalculated. For sufficiently small $$x$$, and sufficiently, but not excessively large $$n$$, the latter can be further simplified to $${-\grave{x}}^{-2}$$, and remains valid when $$x \ge 1$$ is not excessively large. By successively multiplying $${s}^{(j)}(x)$$ by $$x$$ for $$j \in {}^{\omega}\mathbb{N}^{*}$$ and subsequently differentiating, other formulas can be derived for $${s}^{(j+1)}(x)$$, providing an example of divergent series. However, if $${s}^{(0)}(-x)$$ is integrated from 0 to 1 and set $$n := \omega$$, an integral expression for $${_e}\omega + \gamma$$ is obtained for Euler’s constant $$\gamma$$.

L’Hôpital’s rule solves the case of $$x = -1$$. Substituting $$y := -\acute{x}$$, by the binomial series a series is obtained with infinite coefficients (if $${_e}\omega$$ is also expressed as a series, even an expression for $$\gamma$$ is obtained). If the numerator of $${s}^{(0)}(x)$$ is illegitimately simplified, finding incorrect results is risked, especially when $$|x| \ge 1$$. So $${s}^{(0)}(-{e}^{\pi i})$$ is e.g. 0 for odd $$n$$, and 1 for even $$n$$, but not $$\tilde{2}$$.

Series theorem for integrals: DFT forms of TS show for $$c, x \in [0, a], f \in {}^{\nu}\mathcal{C}^{n}$$ and $$n, \hat{p}, r := 2^p \in {}^{\nu}2\mathbb{N}^*$$${\uparrow}_0^1 g(y){\downarrow}y = 2{+}_{\check{q}=1}^{\check{r}}{{+}_{\check{m}=0}^{\check{n}-1}{\widetilde{\grave{m}!}g^{(m)}\left (\tilde{r}\acute{q} \right){\tilde{r}}^{\grave{m}}}}+\mathcal{O}\left(\widetilde{\grave{n}!}g^{(n)}(\tilde{a}c){\tilde{r}}^n\right).\square$Remarks: DFT equivalents may replace any $$g^{(m)}$$ for $$g(y) := af(ay)$$ and $$a > 0$$. Transforming $$z := e^{\pm x}$$ makes here infinite bounds of integration finite. A finite integral decreases the remainder’s modulus sufficiently. The midpoint rule is more advantageous than the trapezoidal rule.

Example: Complete elliptic integrals of the first and second kind are for $$g(\vartheta) := \check{\pi}(1 – {\varepsilon}^2 {\sin}^2(\check{\pi}\vartheta))^{\mp\tilde{2}}$$ and $$\varepsilon \in \; ]0, 1[$$${\uparrow}_0^{1} g(\vartheta) {\downarrow}\vartheta = g\left (\tilde{2} \right) + \widetilde{24} g^{(2)}\left (\tilde{2} \right) + \widetilde{1920} g^{(4)}\left (\tilde{2} \right) + \mathcal{O}\left(\widetilde{9!}g^{(6)}\left (\tilde{2} \right)\right).$Sum formula: The preceding theorem implies for $$f(\check{q}) := g(\tilde{r}\acute{q}), a := r$$ and $$H_{\acute{m}} := 0$$${+}_{\check{q}=1}^{\check{r}} f(\check{q}) = {\uparrow}_{\tilde{2}}^{\check{r}+\tilde{2}} f(x){\downarrow}x – {+}_{\check{m}=0}^{\check{n}-1} H_m \left (f^{(\acute{m})}(\check{r} + \tilde{2}) – f^{(\acute{m})}(\tilde{2}) \right ) + \mathcal{O}\left (H_n \left (f^{(\acute{n})}(\check{r} + \tilde{2}) – f^{(\acute{n})}(\tilde{2}) \right )\right ).$Proof: From${+}_{\check{q}=1}^{\check{r}} g(\tilde{r}\acute{q}) = \tilde{2}r{\uparrow}_0^1 g(y){\downarrow}y – {+}_{\check{q}=1}^{\check{r}}{+}_{\check{m}=0}^{\check{n}-1} {\widetilde{\grave{m}!}g^{(m)}\left (\tilde{r}\acute{q} \right)} + \mathcal{O}\left(\widetilde{\grave{n}!}g^{(n)}(\tilde{a}c){\tilde{r}}^n\right),$the alternating auxiliary numbers $$H_m$$ for the even derivatives $$g^{(m)}$$ are given by inserting and combining the factorials that are counted depending on partitions.$$\square$$

Remark: It is $$H_2 = -1/24, H_4 = 7/5760, H_6 = -31/2017280, H_8 = 127/154828800$$ and so on. Every numerator of uncancelled $$H_m$$ can shown to be divisible by $$2^{\acute{m}} – 1$$ and the limit of $$H_m$$ to be precisely $$2/\hat{\pi}^m$$. The sequence given by $$|H_m|$$ is strictly monotonically decreasing, contrary to Bernoulli numbers.

Corollary: For $$f(n) = {}_en$$, the approximation $$n! \sim \check{5}(\tilde{e}n,5)^{n,5}$$ is more precise than Stirling’s formula.

Proof: Putting $$\check{r} := n$$ and $$n,5 := n + \tilde{2}$$ shows where the limit of the fix part is roughly $$\check{5}$$ that:\begin{aligned}{}_en! &={\uparrow}_{\check{3}}^{n,5} f(x){\downarrow}x – {+}_{\check{m}=1}^n H_m \left (f^{(\acute{m})} (n,5) – f^{(\acute{m})} (\check{3}) \right ) + \mathcal{O}\left( H_{n+2} \left (f^{(\grave{n})} (n,5) – f^{(\grave{n})} (\check{3})\right)\right) \\ &= n,5(_e{(n,5)}- 1) – \check{3}\left (_e{\check{3}} – 1 \right ) – \left (\widetilde{n,5} – \check{3}^{-1} \right )/24 + 7\left (\widetilde{n,5}^3 – \check{3}^{-3} \right )/2880 + \ldots\square\end{aligned}

code of the FFT form