Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <cstdio>
- #include <cmath>
- #include <vector>
- #define a 0
- #define b 1.2
- double Function(double x) {
- return 10 * x * x * cosh(x) * sin(13 * x);
- }
- double Chebyshev(double x, int n) {
- if (n == 0) return 1;
- double tn = x;
- double tnm1 = 1;
- for (int i = 1; i < n; i++) {
- double tnp1 = 2 * x * tn - tnm1;
- tnm1 = tn;
- tn = tnp1;
- }
- return tn;
- }
- double F(double x, int k1, int k2, int N) {
- double t = (2 * x - (a + b)) / (b - a);
- if (k2 == N) {
- return Chebyshev(t, k1) * Function(x);
- }
- else {
- return Chebyshev(t, k1) * Chebyshev(t, k2);
- }
- }
- double Simpson(int k1, int k2, int n, int N) {
- double h = (b - a) / (double)n;
- double s1 = 0;
- double s2 = 0;
- for (int i = 1; i < n; i++) {
- double x = a + i * h;
- if (i % 2 == 0) {
- s2 += F(x, k1, k2, N);
- }
- else {
- s1 += F(x, k1, k2, N);
- }
- }
- return h / 3.0 * (2.0 * s2 + 4.0 * s1 + F(a, k1, k2, N) + F(b, k1, k2, N));
- }
- double Runge(int k1, int k2, int N) {
- int n = (b - a) * 100;
- double I2n = Simpson(k1, k2, n, N);
- for (;;) {
- double In = I2n;
- I2n = Simpson(k1, k2, n *= 2, N);
- if (fabs(In - I2n) < 15e-8) break;
- }
- return I2n;
- }
- void FillPolinom(std::vector<double> &p) {
- int N = p.size();
- double **M = new double*[N];
- for (int i = 0; i < N; i++) {
- M[i] = new double[N + 1]{ 0 };
- }
- for (int i = 0; i < N; i++) {
- for (int j = 0; j <= N; j++) {
- M[i][j] = Runge(i, j, N);
- }
- }
- for (int j = 0; j < N; j++) {
- double m = M[j][j];
- for (int i = j; i <= N; i++) {
- M[j][i] /= m;
- }
- for (int k = j + 1; k < N; k++) {
- double m = M[k][j];
- for (int l = j; l <= N; l++) {
- M[k][l] -= M[j][l] * m;
- }
- }
- }
- for (int i = N - 1; i >= 0; i--) {
- double s = 0;
- for (int j = i + 1; j < N; j++) {
- s += M[i][j] * p[j];
- }
- p[i] = M[i][N] - s;
- }
- for (int i = 0; i < N; i++) {
- delete[] M[i];
- }
- delete[] M;
- }
- double Polinom(double x, const std::vector<double> &p) {
- double t = (2 * x - b - a) / (b - a);
- double s = 0;
- for (int i = 0, N = p.size(); i < N; i++) {
- s += p[i] * Chebyshev(t, i);
- }
- return s;
- }
- std::vector<double> FindPolinom() {
- for (std::vector<double> p(4); ; p.resize(p.size() << 1)) {
- FillPolinom(p);
- double sqrError = 0;
- for (int i = 0; i <= 60; i++) {
- double x = a + i * 0.02;
- double d = Function(x) - Polinom(x, p);
- sqrError += d * d;
- }
- if (sqrError < 61e-6) return p;
- }
- }
- int main() {
- auto const &p = FindPolinom();
- printf("N = %i\n", p.size());
- for (int i = 0; i <= 60; i++) {
- double x = a + i * 0.02;
- printf("%-.2f\t%.8f\n", x, Polinom(x, p));
- }
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement