From 9164518de7acc0a4124ab78872fc376c1bb4e625 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andreas=20M=C3=BCller?= <andreas.mueller@ost.ch>
Date: Sat, 8 Jan 2022 17:42:29 +0100
Subject: add jacobi illustration

---
 .../chapters/070-orthogonalitaet/images/jacobi.cpp | 179 +++++++++++++++++++++
 1 file changed, 179 insertions(+)
 create mode 100644 buch/chapters/070-orthogonalitaet/images/jacobi.cpp

(limited to 'buch/chapters/070-orthogonalitaet/images/jacobi.cpp')

diff --git a/buch/chapters/070-orthogonalitaet/images/jacobi.cpp b/buch/chapters/070-orthogonalitaet/images/jacobi.cpp
new file mode 100644
index 0000000..e010fd6
--- /dev/null
+++ b/buch/chapters/070-orthogonalitaet/images/jacobi.cpp
@@ -0,0 +1,179 @@
+/*
+ * jacobi.cpp
+ *
+ * (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+ */
+#include <cstdlib>
+#include <vector>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <memory>
+#include <getopt.h>
+#include <cmath>
+
+/**
+ * \brief Pfad-Klasse
+ */
+class	pfad {
+	std::ostringstream	*_inhalt;
+	std::string	_name;
+	int	_punkte;
+public:
+	pfad(const pfad&) = delete;
+	pfad(const std::string name) : _name(name), _punkte(0) {
+		_inhalt = new std::ostringstream();
+	}
+	~pfad() {
+		if (_inhalt) {
+			delete _inhalt;
+		}
+	}
+	void	add(const std::pair<double,double>& XY) {
+		if (_punkte == 0) {
+			*_inhalt << "\\def\\" + _name + "{\n\t";
+		} else {
+			*_inhalt << "\n\t-- ";
+		}
+		*_inhalt << "({";
+		*_inhalt << XY.first;
+		*_inhalt << "*\\dx},{";
+		*_inhalt << XY.second;
+		*_inhalt << "*\\dy})";
+		_punkte++;
+	}
+	void	operator()(const std::pair<double,double>& XY) {
+		add(XY);
+	}
+	std::string	macro() const {
+		std::string	result = _inhalt->str();
+		result = result + ";\n}\n";
+		return result;
+	}
+};
+typedef	std::shared_ptr<pfad>	pfad_ptr;
+
+#define	c(k) ((k)+a+b)
+
+std::vector<double>	jacobiP(double a, double b, double x, int n) {
+	std::vector<double>	result;
+	double	p0 = 1.;
+	result.push_back(p0);
+	double	p1 = (a - b)/2. + (1. + (a + b)/2.) * x;
+	result.push_back(p1);
+	int	i = 1;
+	while (i++ < n) {
+		double	p;
+		p = c(2*i-1) * (c(2*i-2) * c(2*i) * x + a*a - b*b) * p1
+			- 2. * (i-1.+a) * (i-1.+b) * c(2*i) * p0;
+		p /= 2 * i * c(i) * c(2*i-2);
+		result.push_back(p);
+		p0 = p1;
+		p1 = p;
+	}
+	return result;
+}
+
+std::string	jacobiPlots(const std::string& prefix,
+			double a, double b, int n, int N) {
+	std::vector<pfad_ptr>	pfade;
+	for (int j = 0; j <= n; j++) {
+		char	c = 'a' + j;
+		std::string	name = prefix;
+		name.append(1, c);
+		pfad_ptr	p(new pfad(name));
+		pfade.push_back(p);
+	}
+	for (int i = -N; i <= N; i++) {
+		double	x = i / (double)N;
+		std::vector<double>	r = jacobiP(a, b, x, n);
+		for (int j = 0; j <= n; j++) {
+			std::pair<double,double>	XY
+				= std::make_pair(x, r[j]);
+			pfade[j]->add(XY);
+		}
+	}
+	std::ostringstream	result;
+	for (int j = 0; j <= n; j++) {
+		result << pfade[j]->macro();
+	}
+	return result.str();
+}
+
+std::string	jacobiWeight(const std::string& prefix,
+			double a, double b, int N) {
+	std::ostringstream	out;
+	out << "\\def\\" << prefix << "weight{" << std::endl;
+	int	_punkte = 0;
+	int	starti = -N;
+	int	endi = N;
+	if (b < 0) {
+		starti += 1;
+	}
+	if (a < 0) {
+		endi -= 1;
+	}
+	for (int i = starti; i <= endi; i++) {
+		double	x = i / (double)N;
+		double	y = pow(1-x, a) * pow(1+x, b);
+		if (_punkte > 0) {
+			out << std::endl << "\t-- ";
+		}
+		out << "({" << x << "*\\dx},{" << y << "*\\dwy})";
+		_punkte++;
+	}
+	out << std::endl << "}" << std::endl;
+	return out.str();
+}
+
+struct option 	options[] = {
+{ "a",		required_argument,	NULL,		'a' },
+{ "b",		required_argument,	NULL,		'b' },
+{ "degree",	required_argument,	NULL,		'n' },
+{ "points",	required_argument,	NULL,		'N' },
+{ "prefix",	required_argument,	NULL,		'p' },
+{ "outfile",	required_argument,	NULL,		'o' },
+{ NULL,		0,			NULL,		 0  }
+};
+
+int	main(int argc, char *argv[]) {
+	int	n = 10;
+	int	N = 100;
+	double	a = 1;
+	double	b = 1;
+	std::string	outfilename;
+	std::string	prefix("prefix");
+	int	c;
+	while (EOF != (c = getopt_long(argc, argv, "a:b:n:N:o:p:", options,
+		NULL)))
+		switch (c) {
+		case 'a':
+			a = std::stod(optarg);
+			break;
+		case 'b':
+			b = std::stod(optarg);
+			break;
+		case 'n':
+			n = std::stoi(optarg);
+			break;
+		case 'N':
+			N = std::stoi(optarg);
+			break;
+		case 'o':
+			outfilename = std::string(optarg);
+			break;
+		case 'p':
+			prefix = std::string(optarg);
+			break;
+		}
+	std::string	result = jacobiPlots(prefix, a, b, n, N);
+	result.append(jacobiWeight(prefix, a, b, N));
+	if (outfilename.size() == 0) {
+		std::cout << result;
+		return EXIT_SUCCESS;
+	}
+	std::ofstream	out(outfilename.c_str());
+	out << result;
+	out.close();
+	return EXIT_SUCCESS;
+}
-- 
cgit v1.2.1