|
63 | 63 | "\n",
|
64 | 64 | "However, for our purposes here it's somewhat more convenient to use the *Newton Polynomial* Basis\n",
|
65 | 65 | "\n",
|
66 |
| - "$$n_j(x) = \\prod^{j-1}_{i=0} (x - x_i)$$" |
| 66 | + "\\begin{align}\n", |
| 67 | + " n_0 &= 1\\\\\n", |
| 68 | + " n_j(x) &= \\prod^{j-1}_{i=0} (x - x_i)\\quad\\mathrm{for}\\quad j>0 \\\\\n", |
| 69 | + "\\end{align}\n" |
67 | 70 | ]
|
68 | 71 | },
|
69 | 72 | {
|
|
96 | 99 | "source": [
|
97 | 100 | "Note: The Newton polynomials have some features common to both the monomials and the Lagrange polynomials\n",
|
98 | 101 | "\n",
|
99 |
| - "* Like the monomials they are of monic polynomials of increasing degree in $x$, i.e. $n_2(x)$ is quadratic\n", |
| 102 | + "* Like the monomials they are monic polynomials of increasing degree in $x$, i.e. $n_2(x)$ is quadratic\n", |
100 | 103 | "* Like the Lagrange polynomials, some Newton polynomials vanish identically at the interpolation points $x_i$. In particular\n",
|
101 | 104 | "\n",
|
102 | 105 | "$$\n",
|
|
118 | 121 | "\n",
|
119 | 122 | "$$P_N(x) = \\sum^N_{j=0} a_j n_j(x)$$\n",
|
120 | 123 | "\n",
|
121 |
| - "as the Linear system\n", |
| 124 | + "as the $(N+1)\\times (N+1)$ linear system\n", |
122 | 125 | "\n",
|
123 | 126 | "$$\n",
|
124 |
| - " P_N(x_i) = y_i\n", |
| 127 | + " P_N(x_i) = y_i\\quad\\mathrm{for}\\quad i=0,\\ldots,N\n", |
125 | 128 | "$$"
|
126 | 129 | ]
|
127 | 130 | },
|
|
279 | 282 | },
|
280 | 283 | "outputs": [],
|
281 | 284 | "source": [
|
282 |
| - "\n", |
283 | 285 | "# Calculate a polynomial in Newton Form\n",
|
284 | 286 | "data = numpy.array([[-2.0, 1.0], [-1.5, -1.0], [-0.5, -3.0], [0.0, -2.0], [1.0, 3.0], [2.0, 1.0]])\n",
|
285 | 287 | "x = numpy.linspace(-2.0, 2.0, 100)\n",
|
286 | 288 | "basis = newton_basis(x, data[:,0])\n",
|
287 |
| - "P = P_newton(x, data)" |
| 289 | + "P = P_newton(x, data)\n", |
| 290 | + "N = data.shape[0] - 1" |
288 | 291 | ]
|
289 | 292 | },
|
290 | 293 | {
|
|
358 | 361 | "cell_type": "markdown",
|
359 | 362 | "metadata": {
|
360 | 363 | "slideshow": {
|
361 |
| - "slide_type": "subslide" |
| 364 | + "slide_type": "fragment" |
362 | 365 | }
|
363 | 366 | },
|
364 | 367 | "source": [
|
365 |
| - "We know from Lagrange's Theorem that the remainder term looks like\n", |
366 |
| - "\n", |
367 |
| - "$$R_N(x) = (x - x_0)(x - x_1)\\cdots (x - x_{N})(x - x_{N+1}) \\frac{f^{(N+1)}(c)}{(N+1)!}$$\n", |
| 368 | + "And we know from Lagrange's Theorem that the remainder term looks like\n", |
368 | 369 | "\n",
|
| 370 | + "$$R_N(x) = (x - x_0)(x - x_1)\\cdots (x - x_{N})(x - x_{N+1}) \\frac{f^{(N+1)}(c)}{(N+1)!}$$" |
| 371 | + ] |
| 372 | + }, |
| 373 | + { |
| 374 | + "cell_type": "markdown", |
| 375 | + "metadata": { |
| 376 | + "slideshow": { |
| 377 | + "slide_type": "subslide" |
| 378 | + } |
| 379 | + }, |
| 380 | + "source": [ |
369 | 381 | "noting that we need to require that $f(x) \\in C^{N+1}$ on the interval of interest. Taking the derivative of the interpolant $P_N(x)$ then leads to \n",
|
370 | 382 | "\n",
|
371 | 383 | "$$P_N'(x) = [y_0, y_1] + ((x - x_1) + (x - x_0)) [y_0, y_1, y_2] + \\cdots + \\left(\\sum^{N-1}_{i=0}\\left( \\prod^{N-1}_{j=0,~j\\neq i} (x - x_j) \\right )\\right ) [y_0, y_1, \\ldots, y_N]$$"
|
|
375 | 387 | "cell_type": "markdown",
|
376 | 388 | "metadata": {
|
377 | 389 | "slideshow": {
|
378 |
| - "slide_type": "subslide" |
| 390 | + "slide_type": "fragment" |
379 | 391 | }
|
380 | 392 | },
|
381 | 393 | "source": [
|
|
401 | 413 | "cell_type": "markdown",
|
402 | 414 | "metadata": {
|
403 | 415 | "slideshow": {
|
404 |
| - "slide_type": "subslide" |
| 416 | + "slide_type": "fragment" |
405 | 417 | }
|
406 | 418 | },
|
407 | 419 | "source": [
|
408 | 420 | "If we let $\\Delta x = \\max_i |x_k - x_i|$ we then know that the remainder term will be $\\mathcal{O}(\\Delta x^N)$ as $\\Delta x \\rightarrow 0$ thus showing that this approach converges and we can find arbitrarily high order approximations (ignoring floating point error)."
|
409 | 421 | ]
|
410 | 422 | },
|
| 423 | + { |
| 424 | + "cell_type": "markdown", |
| 425 | + "metadata": { |
| 426 | + "slideshow": { |
| 427 | + "slide_type": "fragment" |
| 428 | + } |
| 429 | + }, |
| 430 | + "source": [ |
| 431 | + "### <font color='red'>Caution</font>\n", |
| 432 | + "\n", |
| 433 | + "High order does not necessarily imply high-accuracy!" |
| 434 | + ] |
| 435 | + }, |
411 | 436 | {
|
412 | 437 | "cell_type": "code",
|
413 | 438 | "execution_count": null,
|
|
486 | 511 | {
|
487 | 512 | "cell_type": "markdown",
|
488 | 513 | "metadata": {
|
| 514 | + "hide_input": true, |
489 | 515 | "slideshow": {
|
490 | 516 | "slide_type": "slide"
|
491 | 517 | }
|
|
720 | 746 | "source": [
|
721 | 747 | "$$\\begin{aligned}\n",
|
722 | 748 | " P_2'(x) &= [f(x_n), f(x_{n+1})] + ((x - x_n) + (x - x_{n+1})) [f(x_n), f(x_{n+1}), f(x_{n-1})] \\\\\n",
|
723 |
| - " &= \\frac{f(x_{n+1}) - f(x_n)}{x_{n+1} - x_n} + ((x - x_n) + (x - x_{n+1})) \\left ( \\frac{f(x_{n-1}) - f(x_{n+1})}{(x_{n-1} - x_{n+1})(x_{n-1} - x_n)} - \\frac{f(x_{n+1}) - f(x_n)}{(x_{n+1} - x_n)(x_{n-1} - x_n)} \\right )\n", |
| 749 | + " &= \\frac{f(x_{n+1}) - f(x_n)}{x_{n+1} - x_n} + ((x - x_n) + (x - x_{n+1}))\\\\ &* \\left ( \\frac{f(x_{n-1}) - f(x_{n+1})}{(x_{n-1} - x_{n+1})(x_{n-1} - x_n)} - \\frac{f(x_{n+1}) - f(x_n)}{(x_{n+1} - x_n)(x_{n-1} - x_n)} \\right )\n", |
724 | 750 | "\\end{aligned}$$"
|
725 | 751 | ]
|
726 | 752 | },
|
|
777 | 803 | "cell_type": "code",
|
778 | 804 | "execution_count": null,
|
779 | 805 | "metadata": {
|
| 806 | + "hide_input": true, |
780 | 807 | "slideshow": {
|
781 |
| - "slide_type": "skip" |
| 808 | + "slide_type": "subslide" |
782 | 809 | }
|
783 | 810 | },
|
784 | 811 | "outputs": [],
|
|
794 | 821 | "\n",
|
795 | 822 | "# Compute derivative\n",
|
796 | 823 | "f_prime_hat = numpy.empty(x_hat.shape)\n",
|
797 |
| - "f_prime_hat[1:-1] = (f(x_hat[2:]) - f(x_hat[:-2])) / (2 * delta_x)\n", |
| 824 | + "f_prime_hat[1:-1] = (f(x_hat[2:]) - f(x_hat[:-2])) / (2. * delta_x)\n", |
798 | 825 | "\n",
|
799 | 826 | "# Use first-order differences for points at edge of domain\n",
|
800 | 827 | "f_prime_hat[0] = (f(x_hat[1]) - f(x_hat[0])) / delta_x # Forward Difference at x_0\n",
|
|
807 | 834 | "cell_type": "code",
|
808 | 835 | "execution_count": null,
|
809 | 836 | "metadata": {
|
810 |
| - "hide_input": false, |
| 837 | + "hide_input": true, |
811 | 838 | "slideshow": {
|
812 | 839 | "slide_type": "fragment"
|
813 | 840 | }
|
|
818 | 845 | "plt.rcdefaults()\n",
|
819 | 846 | "axes = fig.add_subplot(1, 1, 1)\n",
|
820 | 847 | "\n",
|
821 |
| - "axes.plot(x, f_prime(x), 'k')\n", |
| 848 | + "axes.plot(x, f_prime(x), 'k',label=\"$f'(x)$\")\n", |
822 | 849 | "axes.plot(x_hat, f_prime_hat, 'ro')\n",
|
| 850 | + "axes.plot(x, f(x),'g--',label=\"$f(x)$\")\n", |
| 851 | + "axes.plot(x_hat, f(x_hat), 'go')\n", |
823 | 852 | "axes.set_xlim((x[0], x[-1]))\n",
|
824 | 853 | "# axes.set_ylim((-1.1, 1.1))\n",
|
825 | 854 | "axes.grid()\n",
|
| 855 | + "axes.legend()\n", |
826 | 856 | "\n",
|
827 | 857 | "plt.show()"
|
828 | 858 | ]
|
|
884 | 914 | {
|
885 | 915 | "cell_type": "markdown",
|
886 | 916 | "metadata": {
|
887 |
| - "hide_input": false, |
| 917 | + "hide_input": true, |
888 | 918 | "slideshow": {
|
889 | 919 | "slide_type": "subslide"
|
890 | 920 | }
|
|
907 | 937 | "source": [
|
908 | 938 | "Say we want to derive the second order accurate, first derivative approximation that we just did, this requires the values $(x_{n+1}, f(x_{n+1})$ and $(x_{n-1}, f(x_{n-1})$. We can express these values via our Taylor series approximation above as\n",
|
909 | 939 | "\n",
|
910 |
| - "$$\\begin{aligned}\n", |
| 940 | + "\\begin{aligned}\n", |
911 | 941 | " f(x_{n+1}) &= f(x_n) + (x_{n+1} - x_n) f'(x_n) + \\frac{(x_{n+1} - x_n)^2}{2!} f''(x_n) + \\frac{(x_{n+1} - x_n)^3}{3!} f'''(x_n) + \\mathcal{O}((x_{n+1} - x_n)^4) \\\\\n",
|
912 |
| - " &= f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\n", |
913 |
| - "\\end{aligned}$$\n", |
914 |
| - "\n", |
| 942 | + "\\end{aligned}" |
| 943 | + ] |
| 944 | + }, |
| 945 | + { |
| 946 | + "cell_type": "markdown", |
| 947 | + "metadata": { |
| 948 | + "slideshow": { |
| 949 | + "slide_type": "fragment" |
| 950 | + } |
| 951 | + }, |
| 952 | + "source": [ |
| 953 | + "or\n", |
| 954 | + "\\begin{aligned}\n", |
| 955 | + "&= f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\n", |
| 956 | + "\\end{aligned}" |
| 957 | + ] |
| 958 | + }, |
| 959 | + { |
| 960 | + "cell_type": "markdown", |
| 961 | + "metadata": { |
| 962 | + "slideshow": { |
| 963 | + "slide_type": "subslide" |
| 964 | + } |
| 965 | + }, |
| 966 | + "source": [ |
915 | 967 | "and\n",
|
916 | 968 | "\n",
|
917 |
| - "$$f(x_{n-1}) = f(x_n) + (x_{n-1} - x_n) f'(x_n) + \\frac{(x_{n-1} - x_n)^2}{2!} f''(x_n) + \\frac{(x_{n-1} - x_n)^3}{3!} f'''(x_n) + \\mathcal{O}((x_{n-1} - x_n)^4) $$\n", |
918 |
| - "\n", |
919 |
| - "$$ = f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4) $$" |
| 969 | + "\\begin{align}\n", |
| 970 | + "f(x_{n-1}) &= f(x_n) + (x_{n-1} - x_n) f'(x_n) + \\frac{(x_{n-1} - x_n)^2}{2!} f''(x_n) + \\frac{(x_{n-1} - x_n)^3}{3!} f'''(x_n) + \\mathcal{O}((x_{n-1} - x_n)^4) \n", |
| 971 | + "\\end{align}" |
| 972 | + ] |
| 973 | + }, |
| 974 | + { |
| 975 | + "cell_type": "markdown", |
| 976 | + "metadata": { |
| 977 | + "slideshow": { |
| 978 | + "slide_type": "fragment" |
| 979 | + } |
| 980 | + }, |
| 981 | + "source": [ |
| 982 | + "\\begin{align} \n", |
| 983 | + "&= f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\n", |
| 984 | + "\\end{align}" |
| 985 | + ] |
| 986 | + }, |
| 987 | + { |
| 988 | + "cell_type": "markdown", |
| 989 | + "metadata": { |
| 990 | + "slideshow": { |
| 991 | + "slide_type": "subslide" |
| 992 | + } |
| 993 | + }, |
| 994 | + "source": [ |
| 995 | + "Or all together (for regularly spaced points),\n", |
| 996 | + "\\begin{align} \n", |
| 997 | + "f(x_{n+1}) &= f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\\\\\n", |
| 998 | + "f(x_{n-1})&= f(x_n) - \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) - \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\n", |
| 999 | + "\\end{align}" |
920 | 1000 | ]
|
921 | 1001 | },
|
922 | 1002 | {
|
|
933 | 1013 | " f'(x_n) + R(x_n) = A f(x_{n+1}) + B f(x_n) + C f(x_{n-1})\n",
|
934 | 1014 | "$$\n",
|
935 | 1015 | "\n",
|
936 |
| - "where $R(x_n)$ is our error. Plugging in the Taylor series approximations we find\n", |
| 1016 | + "where $R(x_n)$ is our error. " |
| 1017 | + ] |
| 1018 | + }, |
| 1019 | + { |
| 1020 | + "cell_type": "markdown", |
| 1021 | + "metadata": { |
| 1022 | + "slideshow": { |
| 1023 | + "slide_type": "fragment" |
| 1024 | + } |
| 1025 | + }, |
| 1026 | + "source": [ |
| 1027 | + "Plugging in the Taylor series approximations we find\n", |
937 | 1028 | "\n",
|
938 | 1029 | "$$\\begin{aligned}\n",
|
939 | 1030 | " f'(x_n) + R(x_n) &= A \\left ( f(x_n) + \\Delta x f'(x_n) + \\frac{\\Delta x^2}{2!} f''(x_n) + \\frac{\\Delta x^3}{3!} f'''(x_n) + \\mathcal{O}(\\Delta x^4)\\right ) \\\\\n",
|
|
985 | 1076 | "cell_type": "markdown",
|
986 | 1077 | "metadata": {
|
987 | 1078 | "slideshow": {
|
988 |
| - "slide_type": "subslide" |
| 1079 | + "slide_type": "skip" |
989 | 1080 | }
|
990 | 1081 | },
|
991 | 1082 | "source": [
|
|
1072 | 1163 | "cell_type": "code",
|
1073 | 1164 | "execution_count": null,
|
1074 | 1165 | "metadata": {
|
| 1166 | + "hide_input": true, |
1075 | 1167 | "slideshow": {
|
1076 |
| - "slide_type": "skip" |
| 1168 | + "slide_type": "fragment" |
1077 | 1169 | }
|
1078 | 1170 | },
|
1079 | 1171 | "outputs": [],
|
|
0 commit comments