summaryrefslogtreecommitdiff
path: root/src/basis.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/basis.c')
-rw-r--r--src/basis.c236
1 files changed, 220 insertions, 16 deletions
diff --git a/src/basis.c b/src/basis.c
index 5ea043e..20e9f12 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -1,11 +1,34 @@
+/*
+ * Basis sets for pseudospectral methods
+ * Copyright (C) 2016 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
#include <math.h>
-#include "maximal_slicing_axi.h"
+#include "basis.h"
+#include "common.h"
+/*
+ * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9)
+ * SB(x, n) = sin((n + 1) arccot(|x| / L))
+ * They are symmetric wrt origin and decay as 1/x in infinity.
+ */
static double sb_even_eval(double coord, int idx)
{
- double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+ double val = atan2(SCALE_FACTOR, coord);
idx *= 2; // even only
@@ -14,20 +37,20 @@ static double sb_even_eval(double coord, int idx)
static double sb_even_eval_diff1(double coord, int idx)
{
- double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+ double val = atan2(SCALE_FACTOR, coord);
idx *= 2; // even only
- return - SCALE_FACTOR * (idx + 1) * SGN(coord) * cos((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
+ return - SCALE_FACTOR * (idx + 1) * cos((idx + 1) * val) / (SQR(SCALE_FACTOR) + SQR(coord));
}
static double sb_even_eval_diff2(double coord, int idx)
{
- double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
+ double val = atan2(SCALE_FACTOR, coord);
idx *= 2; // even only
- return SCALE_FACTOR * (idx + 1) * SGN(coord) * (2 * fabs(coord) * cos((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
+ return SCALE_FACTOR * (idx + 1) * (2 * coord * cos((idx + 1) * val) - SCALE_FACTOR * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(SCALE_FACTOR) + SQR(coord));
}
static double sb_even_colloc_point(int order, int idx)
@@ -38,11 +61,15 @@ static double sb_even_colloc_point(int order, int idx)
//order *= 2;
//t = (idx + 2) * M_PI / (order + 4);
+#if MS_POLAR
+ t = (idx + 2) * M_PI / (2 * order + 3);
+#else
t = (idx + 2) * M_PI / (2 * order + 2);
+#endif
return SCALE_FACTOR / tan(t);
}
-const BasisSet msa_sb_even_basis = {
+const BasisSet ms_sb_even_basis = {
.eval = sb_even_eval,
.eval_diff1 = sb_even_eval_diff1,
.eval_diff2 = sb_even_eval_diff2,
@@ -87,17 +114,77 @@ static double tb_even_colloc_point(int order, int idx)
//order *= 2;
//t = (idx + 2) * M_PI / (order + 4);
- t = (idx + 3) * M_PI / (2 * (order + 2));
+ t = (idx + 2) * M_PI / (2 * order + 4);
return SCALE_FACTOR / tan(t);
}
-const BasisSet msa_tb_even_basis = {
+const BasisSet ms_tb_even_basis = {
.eval = tb_even_eval,
.eval_diff1 = tb_even_eval_diff1,
.eval_diff2 = tb_even_eval_diff2,
.colloc_point = tb_even_colloc_point,
};
+static double tl_eval(double coord, int idx)
+{
+ double t = 2 * atan2(sqrt(SCALE_FACTOR), sqrt(fabs(coord)));
+
+ //idx++;
+
+ return cos(idx * t);// - 1.0;
+}
+
+static double tl_eval_diff1(double coord, int idx)
+{
+ double y = fabs(coord);
+ double y12 = sqrt(y);
+ double t = 2 * atan2(sqrt(SCALE_FACTOR), y12);
+
+ //idx++;
+
+ if (y < 1e-10)
+ return -2.0 * SQR(idx) * pow(-1.0, idx) / SCALE_FACTOR;
+
+ return idx * sqrt(SCALE_FACTOR) * sin(idx * t) / (y12 * (y + SCALE_FACTOR));
+}
+
+static double tl_eval_diff2(double coord, int idx)
+{
+ double y = fabs(coord);
+ double y12 = sqrt(y);
+ double t = 2 * atan2(sqrt(SCALE_FACTOR), y12);
+
+ //idx++;
+
+ if (y < 1e-10)
+ return 4.0 * SQR(idx) * (SQR(idx) + 2) * pow(-1.0, idx) / (3.0 * SQR(SCALE_FACTOR));
+
+ return -(SCALE_FACTOR * SQR(idx) * cos(idx * t) / (y * SQR(y + SCALE_FACTOR)) + sqrt(SCALE_FACTOR) * idx * sin(idx * t) / (y12 * SQR(y + SCALE_FACTOR)) + sqrt(SCALE_FACTOR) * idx * sin(idx * t) / (2 * y * y12 * (y + SCALE_FACTOR)));
+}
+
+static double tl_colloc_point(int order, int idx)
+{
+ double t;
+
+ if (!idx)
+ return 0.0;
+
+ idx = order - idx - 1;
+
+ order++;
+ idx++;
+
+ t = (2 * idx + 1) * M_PI / (2 * order);
+ return SCALE_FACTOR / SQR(tan(0.5 * t));
+}
+
+const BasisSet ms_tl_basis = {
+ .eval = tl_eval,
+ .eval_diff1 = tl_eval_diff1,
+ .eval_diff2 = tl_eval_diff2,
+ .colloc_point = tl_colloc_point,
+};
+
static double full_eval(double coord, int idx)
{
double val = (coord == 0.0) ? M_PI_2 : atan(SCALE_FACTOR / fabs(coord));
@@ -153,7 +240,7 @@ static double full_colloc_point(int order, int idx)
}
-const BasisSet msa_full_basis = {
+const BasisSet ms_full_basis = {
.eval = full_eval,
.eval_diff1 = full_eval_diff1,
.eval_diff2 = full_eval_diff2,
@@ -162,28 +249,145 @@ const BasisSet msa_full_basis = {
static double cheb_eval(double coord, int idx)
{
- return cos(2 * idx * acos(coord / SCALE_FACTOR));
+ double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));
+
+ //idx += 1;
+ return cos(idx * acos(x));
}
static double cheb_eval_diff1(double coord, int idx)
{
- return 2 * idx * sin(2 * idx * acos(coord / SCALE_FACTOR)) / sqrt(SQR(SCALE_FACTOR) - SQR(coord));
+ double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));
+ double t = acos(x);
+
+ //idx += 1;
+
+ if (fabs(x - 1.0) < 1e-8)
+ return (2.0 / SCALE_FACTOR) * SQR(idx);
+ if (fabs(x + 1.0) < 1e-8)
+ return -(2.0 / SCALE_FACTOR) * SQR(idx) * pow(-1.0, idx);
+ return (2.0 / SCALE_FACTOR) * idx * sin(idx * t) / sin(t);
}
static double cheb_eval_diff2(double coord, int idx)
{
- double t = acos(coord / SCALE_FACTOR);
- return 2 * idx * (cos(2 * idx * t) * 2 * idx / (SQR(SCALE_FACTOR) - SQR(coord)) + sin(2 * idx * t) * coord / pow(SQR(SCALE_FACTOR) - SQR(coord), 1.5));
+ double x = MIN(1.0, MAX(2.0 * coord / SCALE_FACTOR - 1.0, -1.0));
+ double t = acos(x);
+ double st = sin(t);
+
+ //idx += 1;
+
+ if (fabs(x - 1.0) < 1e-8)
+ return SQR(2.0 / SCALE_FACTOR) * SQR(idx) * (SQR(idx) - 1.0) / 3.;
+ if (fabs(x + 1.0) < 1e-8)
+ return SQR(2.0 / SCALE_FACTOR) * pow(-1.0, idx) * SQR(idx) * (SQR(idx) - 1.0) / 3.;
+
+ return SQR(2.0 / SCALE_FACTOR) * (-SQR(idx) * cos(idx * t) / SQR(st) + idx * cos(t) * sin(idx * t) / (st * SQR(st)));
}
static double cheb_colloc_point(int order, int idx)
{
- return SCALE_FACTOR * cos((idx + 0.01) * M_PI / (4 * order + 4));
+ double y;
+
+ idx = order - idx - 1;
+
+ //order += 1;
+ //idx += 1;
+ y = cos(idx * M_PI / (order - 1));
+ //y = cos((2 * idx + 1) * M_PI / (2 * order));
+
+ return 0.5 * SCALE_FACTOR * (y + 1.0);
}
-const BasisSet msa_cheb_basis = {
+const BasisSet ms_cheb_basis = {
.eval = cheb_eval,
.eval_diff1 = cheb_eval_diff1,
.eval_diff2 = cheb_eval_diff2,
.colloc_point = cheb_colloc_point,
};
+
+static double cheb_even_eval(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));
+
+ idx *= 2;
+ //idx += 1;
+ return cos(idx * acos(x));
+}
+
+static double cheb_even_eval_diff1(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));
+ double t = acos(x);
+
+ idx *= 2;
+ //idx += 1;
+
+ if (fabs(fabs(x) - 1.0) < 1e-8)
+ return (1.0 / SCALE_FACTOR) * SQR(idx);
+ return (1.0 / SCALE_FACTOR) * idx * sin(idx * t) / sin(t);
+}
+
+static double cheb_even_eval_diff2(double coord, int idx)
+{
+ double x = MIN(1.0, MAX(coord / SCALE_FACTOR, -1.0));
+ double t = acos(x);
+ double st = sin(t);
+
+ idx *= 2;
+ //idx += 1;
+
+ if (fabs(fabs(x) - 1) < 1e-8)
+ return SQR(1.0 / SCALE_FACTOR) * SQR(idx) * (SQR(idx) - 1.0) / 3.;
+
+ return SQR(1.0 / SCALE_FACTOR) * (-SQR(idx) * cos(idx * t) / SQR(st) + idx * cos(t) * sin(idx * t) / (st * SQR(st)));
+}
+
+static double cheb_even_colloc_point(int order, int idx)
+{
+ double y;
+
+ idx = order - idx - 1;
+
+ //order += 1;
+ //idx += 1;
+ order *= 2;
+ y = cos(idx * M_PI / (order - 1));
+ y = cos((2 * idx + 2) * M_PI / (2 * order));
+
+ return SCALE_FACTOR * y;
+}
+
+const BasisSet ms_cheb_even_basis = {
+ .eval = cheb_even_eval,
+ .eval_diff1 = cheb_even_eval_diff1,
+ .eval_diff2 = cheb_even_eval_diff2,
+ .colloc_point = cheb_even_colloc_point,
+};
+
+static double cos_even_eval(double coord, int idx)
+{
+ return cos(2 * idx * coord);
+}
+
+static double cos_even_eval_diff1(double coord, int idx)
+{
+ return -2 * idx * sin(2 * idx * coord);
+}
+
+static double cos_even_eval_diff2(double coord, int idx)
+{
+ return -4 * SQR(idx) * cos(2 * idx * coord);
+}
+
+static double cos_even_colloc_point(int order, int idx)
+{
+ return M_PI * idx / (2 * order - 0);
+}
+
+const BasisSet ms_cos_even_basis = {
+ .eval = cos_even_eval,
+ .eval_diff1 = cos_even_eval_diff1,
+ .eval_diff2 = cos_even_eval_diff2,
+ .colloc_point = cos_even_colloc_point,
+};
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 12:22:54 +0000