From: Martin Kronbichler Date: Fri, 12 Jul 2019 07:56:37 +0000 (+0200) Subject: Switch Hermite-like interpolation to Jacobi polynomials P^{4,4} X-Git-Tag: v9.2.0-rc1~1378^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6b877bf84eb2f94a500aa2ab50071b74a52f8e83;p=dealii.git Switch Hermite-like interpolation to Jacobi polynomials P^{4,4} --- diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index 4ba120336b..3d99c7010b 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -617,11 +617,12 @@ namespace Polynomials * respectively. Furthermore, the extension of these polynomials to higher * degrees $n>3$ is constructed by adding additional nodes inside the unit * interval, again ensuring better conditioning. The nodes are computed as - * the roots of the Jacobi polynomials for $\alpha=\beta=2$, which are - * orthogonal against the generating function $x^2(1-x)^2$ with the Hermite + * the roots of the Jacobi polynomials for $\alpha=\beta=4$, which are + * orthogonal against the square of the generating function $x^2(1-x)^2$ + * with the Hermite * property. Then, these polynomials are constructed in the usual way as - * Lagrange polynomials with double roots at $x=0$ and $x=1$. For example at - * $n=4$, all of $p_0, p_1, p_3, p_4$ get an additional root at $x=0.5$ + * Lagrange polynomials with double roots at $x=0$ and $x=1$. For example + * with $n=4$, all of $p_0, p_1, p_3, p_4$ get an additional root at $x=0.5$ * through the factor $(x-0.5)$. In summary, this basis is dominated by * nodal contributions, but it is not a nodal one because the second and * second to last polynomials that are non-nodal, and due to the presence of @@ -629,12 +630,11 @@ namespace Polynomials * set such that the sum of all polynomials with unit weight represents the * constant function 1, similarly to Lagrange polynomials. * - * The basis only contains Hermite information at degree>=3, + * The basis only contains Hermite information for degree>=3, * but it is also implemented for degrees between 0 and two. For the linear * case, the usual hat functions are implemented, whereas the polynomials * for degree=2 are $p_0(x)=(1-x)^2$, $p_1(x)=2x(x-1)$, and - * $p_2(x)=x^2$, in accordance with the construction principle for degree 3 - * that allows a non-zero of $p_0$ and $p_2$. + * $p_2(x)=x^2$, in accordance with the construction principle for degree 3. * * These two relaxations improve the condition number of the mass matrix * (i.e., interpolation) significantly, as can be seen from the following @@ -663,32 +663,32 @@ namespace Polynomials * * n=5 * 1.875e+04 - * 19.37 + * 15.99 * * * n=6 * 6.033e+04 - * 18.99 + * 16.34 * * * n=10 * 9.756e+05 - * 25.65 + * 20.70 * * * n=15 * 9.431e+06 - * 36.47 + * 27.91 * * * n=25 * 2.220e+08 - * 62.28 + * 43.54 * * * n=35 * 2.109e+09 - * 91.50 + * 59.51 * * * diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index 0ec454b3b2..3eacc200a2 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -1389,22 +1389,24 @@ namespace Polynomials // // | x 0 x x x x x | // | 0 x x x . . . x x x | - // | x x x x x x x | - // | x x x x x x x | + // | x x x 0 0 x x | + // | x x 0 x 0 x x | // | . . . | // M = | . . . | // | . . . | - // | x x x x x x x | + // | x x 0 0 x x x | // | x x x x . . . x x 0 | // | x x x x x 0 x | // // We find the inner points as the zeros of the Jacobi polynomials - // with alpha = beta = 2 which is the polynomial with the kernel - // (1-x)^2 (1+x)^2, the two polynomials achieving zero value and zero - // derivative at the boundary. + // with alpha = beta = 4 which is the polynomial with the kernel + // (1-x)^4 (1+x)^4. Since polynomials (1-x)^2 (1+x)^2 are contained + // in every interior polynomial (bubble function), their product + // leads us to the orthogonality condition of the Jacobi(4,4) + // polynomials. std::vector jacobi_roots = - jacobi_polynomial_roots(degree - 3, 2, 2); + jacobi_polynomial_roots(degree - 3, 4, 4); AssertDimension(jacobi_roots.size(), degree - 3); // iteration from variable support point N with secant method diff --git a/tests/base/polynomial_hermite_like.output b/tests/base/polynomial_hermite_like.output index 6e48f5b9a5..423b2d3799 100644 --- a/tests/base/polynomial_hermite_like.output +++ b/tests/base/polynomial_hermite_like.output @@ -62,9 +62,9 @@ DEAL::0.875000 0.0498047 -0.102539 0.191406 0.717773 0.143555 DEAL::1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 DEAL:: DEAL::degree 6 -DEAL::0 0: 1.00000 -21.1429 +DEAL::0 0: 1.00000 -21.7000 DEAL::0 1: 0.00000 0.00000 -DEAL::1 0: 0.00000 21.1429 +DEAL::1 0: 0.00000 21.7000 DEAL::1 1: 0.00000 0.00000 DEAL::2 0: 0.00000 0.00000 DEAL::2 1: 0.00000 0.00000 @@ -73,23 +73,23 @@ DEAL::3 1: 0.00000 0.00000 DEAL::4 0: 0.00000 0.00000 DEAL::4 1: 0.00000 0.00000 DEAL::5 0: 0.00000 0.00000 -DEAL::5 1: 0.00000 -21.1429 +DEAL::5 1: 0.00000 -21.7000 DEAL::6 0: 0.00000 0.00000 -DEAL::6 1: 1.00000 21.1429 +DEAL::6 1: 1.00000 21.7000 DEAL::0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -DEAL::0.125000 -0.0775452 0.521667 0.643097 -0.131592 0.0836483 -0.0745239 0.0352478 -DEAL::0.250000 0.0627790 -0.185826 1.02264 0.140625 -0.0734223 0.0619420 -0.0287388 -DEAL::0.375000 0.126103 -0.314549 0.613544 0.714111 -0.242755 0.188729 -0.0851833 +DEAL::0.125000 -0.139293 0.671701 0.646751 -0.275146 0.141723 -0.0959573 0.0502213 +DEAL::0.250000 -0.0149414 0.0381445 1.00964 -0.0468750 0.0202009 -0.0127148 0.00654297 +DEAL::0.375000 0.0787735 -0.173836 0.587642 0.640869 -0.185359 0.104301 -0.0523911 DEAL::0.500000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 -DEAL::0.625000 -0.0851833 0.188729 -0.242755 0.714111 0.613544 -0.314549 0.126103 -DEAL::0.750000 -0.0287388 0.0619420 -0.0734223 0.140625 1.02264 -0.185826 0.0627790 -DEAL::0.875000 0.0352478 -0.0745239 0.0836483 -0.131592 0.643097 0.521667 -0.0775452 +DEAL::0.625000 -0.0523911 0.104301 -0.185359 0.640869 0.587642 -0.173836 0.0787735 +DEAL::0.750000 0.00654297 -0.0127148 0.0202009 -0.0468750 1.00964 0.0381445 -0.0149414 +DEAL::0.875000 0.0502213 -0.0959573 0.141723 -0.275146 0.646751 0.671701 -0.139293 DEAL::1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 DEAL:: DEAL::degree 9 -DEAL::0 0: 1.00000 -43.0000 +DEAL::0 0: 1.00000 -46.0000 DEAL::0 1: 0.00000 0.00000 -DEAL::1 0: 0.00000 43.0000 +DEAL::1 0: 0.00000 46.0000 DEAL::1 1: 0.00000 0.00000 DEAL::2 0: 0.00000 0.00000 DEAL::2 1: 0.00000 0.00000 @@ -104,15 +104,15 @@ DEAL::6 1: 0.00000 0.00000 DEAL::7 0: 0.00000 0.00000 DEAL::7 1: 0.00000 0.00000 DEAL::8 0: 0.00000 0.00000 -DEAL::8 1: 0.00000 -43.0000 +DEAL::8 1: 0.00000 -46.0000 DEAL::9 0: 0.00000 0.00000 -DEAL::9 1: 1.00000 43.0000 +DEAL::9 1: 1.00000 46.0000 DEAL::0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -DEAL::0.125000 0.0734501 -0.287123 1.00390 0.270272 -0.0928510 0.0555242 -0.0438505 0.0446656 -0.0410176 0.0170339 -DEAL::0.250000 -0.0360489 0.103340 -0.157805 1.00328 0.120509 -0.0544771 0.0391483 -0.0382418 0.0344467 -0.0141525 -DEAL::0.375000 -0.0402864 0.106060 -0.136357 0.215809 0.941854 -0.131662 0.0793256 -0.0726297 0.0636361 -0.0257504 -DEAL::0.500000 0.0664062 -0.167969 0.200133 -0.244967 0.646397 0.646397 -0.244967 0.200133 -0.167969 0.0664062 -DEAL::0.625000 -0.0257504 0.0636361 -0.0726297 0.0793256 -0.131662 0.941854 0.215809 -0.136357 0.106060 -0.0402864 -DEAL::0.750000 -0.0141525 0.0344467 -0.0382418 0.0391483 -0.0544771 0.120509 1.00328 -0.157805 0.103340 -0.0360489 -DEAL::0.875000 0.0170339 -0.0410176 0.0446656 -0.0438505 0.0555242 -0.0928510 0.270272 1.00390 -0.287123 0.0734501 +DEAL::0.125000 -0.0171645 0.0438649 1.01661 -0.0656195 0.0365983 -0.0236670 0.0151330 -0.00863548 0.00626642 -0.00338620 +DEAL::0.250000 0.0176958 -0.0370004 0.0912862 1.00058 -0.107265 0.0550291 -0.0321169 0.0175116 -0.0123335 0.00661360 +DEAL::0.375000 -0.0273213 0.0538620 -0.0970813 0.310409 0.891621 -0.193262 0.0945173 -0.0479195 0.0323172 -0.0171422 +DEAL::0.500000 0.0267857 -0.0513393 0.0815481 -0.184165 0.627170 0.627170 -0.184165 0.0815481 -0.0513393 0.0267857 +DEAL::0.625000 -0.0171422 0.0323172 -0.0479195 0.0945173 -0.193262 0.891621 0.310409 -0.0970813 0.0538620 -0.0273213 +DEAL::0.750000 0.00661360 -0.0123335 0.0175116 -0.0321169 0.0550291 -0.107265 1.00058 0.0912862 -0.0370004 0.0176958 +DEAL::0.875000 -0.00338620 0.00626642 -0.00863548 0.0151330 -0.0236670 0.0365983 -0.0656195 1.01661 0.0438649 -0.0171645 DEAL::1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000