Приближение полиномами Чебышева

@Milkho · Регистрация: 24.07.2015

Студворк — интернет-сервис помощи студентам

Нужно програмно реализовать приближения (интерполяцию) функции полинонами Чебышева. Т.е. у нас есть n-ое количество пар x,y, и нужно по заданному x получить приблизительный y.

Приближающая функция ищется в виде суммы многочленов Чебышева, формулы во вложениях.

Пробовал реализовать, пока безуспешно:

C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
 public class Polynom
    {
        public double[] X { get; set; }
        public double[] Y { get; set; }
 
        public Polynom(double[] x, double[] y)
        {
            X = x;
            Y = y;
        }
 
        public Polynom(List<DoublePoint> list)
        {
            X = new double[list.Count];
            Y = new double[list.Count];
            int j = 0;
            foreach (var i in list)
            {
                X[j] = i.X;
                Y[j] = i.Y;
                j++;
            }
        }
 
        public double GetY(double x)
        {
            double res = GetRatio(0, X.Length) * (1 / Math.Pow(2, 0.5)); //отдельно считаем первое слогаемое (для Т с тильдой)
            for (int i = 1; i < X.Length - 1; i++)
            {
                res += (GetRatio(i, X.Length) * GetChebPolynom(x, i)); //сумма многочленов
            }
            return res; 
        }
 
        public double GetChebPolynom(double x, long n) //считаем многочлен
        {
            if (n == 1) return x;
            double result = Math.Cos(n * Math.Acos(x));
            return result;
        }
 
        public double GetRatio(double j, double n)
        {
            double res = Y[0] * 1; //отдельно считаем первое слогаемое (для Т без тильды)
        for (var i = 1; i < n - 1; i++)
            {
                res +=  Y[i] * GetChebPolynom(X[i], (long)j);
            }
            return res * 2 / n;
        }
 
    }

Класс DoublePoint:

C#
1
2
3
4
5
6
7
8
9
10
11
public class DoublePoint : IComparable
{
    public double X { get; set; }
    public double Y { get; set; }
 
    public DoublePoint(double x, double y)
    {
        X = x;
        Y = y;
    }
}

Пожалуйста, укажите мне на мои ошибки. Буду рад любой помощи.

Возможно, я не до конца понял этот метод. Ведь я так и не понял, для чего формула, которая вычисляет Xj. И не знаю, правильно ли я считаю полиномы.

Заранее спасибо.

@Milkho · 21.05.2016, 15:02 **[ТС]**

Работает, но неправильно считает. Знающие люди, подскажите, пожалуйста, где ошибка.

C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
  public class Polynom
    {
        public double[] X { get; set; }
        public double[] Y { get; set; }
 
        public Polynom(double[] x, double[] y)
        {
            X = x;
            Y = y;
        }
 
        public Polynom(List<DoublePoint> list)
        {
            X = new double[list.Count];
            Y = new double[list.Count];
            int j = 0;
            foreach (var i in list)
            {
                X[j] = i.X;
                Y[j] = i.Y;
                j++;
            }
        }
 
        public double GetY(double x)
        {
            double res = 0;
            for (int i = 0; i < X.Length - 1; i++)
            {
                res += (GetRatio(i, X.Length) * GetChebPolynom(x, i));
            }
            return res;
        }
 
        public double GetChebPolynom(double x, long n) 
        {
            if (n == 0) return 1;
            if (n == 1) return x;
            double result = 2 * x * x - 1;
            if (n == 2) return result;
            for (int j = 3; j <= n; j++)
                result = result * 2 * x;
            return result;
        }
        public double GetXCheb(int m)
        {
            double a = X.Min();
            double b = X.Max();
            int n = X.Length - 1;
            double res = ((a + b + (b - a) * Math.Cos(2 * Math.PI + m - 1) / 2 * n)) / 2;
            return res;
        }
        public double GetRatio(double j, double n)
        {
            double res = 0; 
            for (var k = 1; k < n - 1; k++)
            {
                //double xk = GetXCheb(k);
                res += Y[k] * GetChebPolynom(X[k], (long)j);
            }
            return res * 2 / n;
        }
 
    }

Формулы взяты остюда - http://e-science.ru/node/125491

@afront · 21.05.2016, 15:47

Пример на c#
https://code.msdn.microsoft.co... d-b10563c8

@Milkho · 21.05.2016, 16:50 **[ТС]**

afront, Cпасибо, но в том коде сложно разобраться.

Как я понял, там интерполируется заданная функция. Но что значит массив коэффициентов (ChebyshevCalculateCoefficients) ? Как их использовать для нахождения значения функции при заданном Х?

@afront · 21.05.2016, 17:28

мне кажется там находятся коэффициенты разложения для заданной функции
вот в этом куске, из баттонт1_клик

C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
            else if (method == 0)
            {
                double[] fy = new double[number + 1];
                double h = 3.2 / number;
                double x = -1.6;
 
                for (int i = 0; i <= number; i++)
                {
                    fy[i] = f(x);
                    x += h;
                }
 
                inter = new Interpolate();
                coeff = inter.ChebyshevCalculateCoefficients(number, fy);
            }

@Milkho · 21.05.2016, 17:37 **[ТС]**

afront, Ну вот, я о них и говорю. Это, значит, коэффициенты возле многочленов Чебышева в разложении?
Я подставил этот метод поиска коэффициентов в свою формулу. Считает по-другому, но все равно, неправильно.

C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
   public class Polynom
    {
        public double[] X { get; set; }
        public double[] Y { get; set; }
 
        public Polynom(double[] x, double[] y)
        {
            X = x;
            Y = y;
        }
 
        //public Polynom() { }
 
        public Polynom(List<DoublePoint> list)
        {
            X = new double[list.Count];
            Y = new double[list.Count];
            int j = 0;
            foreach (var i in list)
            {
                X[j] = i.X;
                Y[j] = i.Y;
                j++;
            }
        }
 
        //public long Fact(long x)
        //{
        //    return (x == 0) ? 1 : x * Fact(x - 1);
        //}
 
        public double GetY(double x)
        {
            // double res = GetRatio(0, X.Length) * (1 / Math.Pow(2, 0.5)); //отдельно считаем первое слогаемое (для Т с тильдой)
            double res = 0;
            for (int i = 0; i < X.Length - 1; i++)
            {
                res += (ChebyshevCalculateCoefficients(i, X.Length) * GetChebPolynom(x, i)) - ChebyshevCalculateCoefficients(0, X.Length) / 2; //сумма многочленов
            }
            return res;
        }
 
        public double GetChebPolynom(double x, long n) 
        {
            if (n == 0) return 1;//Math.Pow(2, 0.5);
            if (n == 1) return x;
            double result = 2 * x * x - 1;
            if (n == 2) return result;
            double temp = 1;
            double res = x;
            for (int j = 3; j < n - 1; j++)
            {
                result = res * 2 * x - temp;
                temp = res;
                res = result;
            }
            return result;
        }
 
        public double ChebyshevCalculateCoefficients(double k, double n)
        {
            double v =0;
 
                double sum = 0.0;
                double arg = k * Math.PI / n;
 
                for (int j = 0; j < n; j++)
                {
                    double dj = 1.0;
 
                    if (j % 2 == 0)
                        dj = 2.0;
 
                    sum += Math.Cos(j * arg) * X[j] / dj;
                }
 
                double dk = 1.0;
 
                if (k % 2 == 0)
                    dk = 2.0;
 
                v = 2.0 * sum / (n * dk);
 
            return v;
        }
 
    }

@afront · 21.05.2016, 17:42

а какую функцию вам нужно интерполировать вы не написали?

@Milkho · 21.05.2016, 17:50 **[ТС]**

afront, У меня нет функции. Есть просто дискретный набор пар точек: X и F(X). Их n штук.

@afront · 21.05.2016, 18:06

хотите на fortrane рабочий код?
там есть и С и С++
http://people.sc.fsu.edu/~jbur... yshev.html

Добавлено через 7 минут
единственно там функция задается аналитически а не таблично
может в этом дело, сделайте тест - задайте известную функцию аналитически и таблично и проверьте

@Milkho · 21.05.2016, 18:10 **[ТС]**

afront, спасибо еще раз.
Но и это не совсем то, что мне нужно. Интерполяция там происходит не по заданным точкам. Вместо этого, алгоритм использует "свои" узлы (нули полинома). К тому же мне придется очень долго изучать весь код этой библиотеки.

@afront · 21.05.2016, 18:12

Milkho, код там элементарный, хотите переведу на c#? с известной функцией?

@Milkho · 21.05.2016, 18:14 **[ТС]**

может в этом дело, сделайте тест - задайте известную функцию аналитически и таблично и проверьте

afront, Что проверить? На примере моей программы? Ну так я так и проверяю, на примере параболы.

Добавлено через 1 минуту

код там элементарный, хотите переведу на c#? с известной функцией?

afront, хочу. По видимому я увидел не тот код. Или не нашел нашу функцию.

@afront · 21.05.2016, 18:21

потребуется время))

@Milkho · 21.05.2016, 18:26 **[ТС]**

afront, я подожду, спасибо огромное вам за помощь.

@afront · 21.05.2016, 19:55

вот что то получилось

C#
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
 class Program
    {
        static double pi = 3.141592653589793;
        static void Main(string[] args)
        {
            int n = 5;
            double a = -1.0;
            double b = +1.0;
            double [] c = new double[n];
            double [] x = new double[n];
            double [] fc = new double[n];
 
            chebyshev_coefficients ( a, b, n, c );
            chebyshev_zeros ( n, x );
 
             for ( int i = 0; i < n; i++)
             {
               x[i] = 0.5 * ( a + b ) + x[i] * 0.5 * ( b - a );
             }
 
           int m = n;
           chebyshev_interpolant ( a, b, n, c, m, x, fc );
           for ( int i = 0; i < n; i++)
             {
                 Console.WriteLine("{0}  {1}  {2}  {3}  {4}", i, x[i] , c[i] , f(x[i]) , fc[i]);
             }
            Console.ReadKey();
        }
        public static double f(double x)
        {
            return (Math.Cos(2.0 * pi * x) * Math.Sin(3.0 * pi * x));
        }
 
        public static void chebyshev_coefficients(double a, double b, int n, double[] c)
        {
            double angle;
            double x;
            double[] fx = new double[n];
 
            for ( int i = 0; i < n; i++)
            {
              int j = i + 1;
              angle = (2.0 * j -1) * pi /  (2.0 * n );
              x = Math.Cos ( angle );
              x = 0.5 * ( a + b ) + x * 0.5 * ( b - a );
              fx[i] = f ( x );
            }
            for (int i = 0; i < n; i++)
            {
                c[i] = 0.0;
                for (int j = 0; j < n; j++)
                {
                 int k = j + 1;
                 angle =  i * (2.0 * k - 1.0) * pi / ( 2.0 * n );
                 c[i] = c[i] + fx[j] * Math.Cos(angle);
                }
            }
            for (int i = 0; i < n; i++)
            {
                c[i] = 2.0 * c[i] / n;
            }
        }
        public static void chebyshev_interpolant(double a, double b, int n, double[] c, int m, double[] x, double[] cf)
        {
            for (int j = 0; j < m; j++)
            { 
                double dip1 = 0.0;
                double di = 0.0;
                double y = (2.0 * x[j] - a - b) / (b - a);
                for (int i = n - 1; i > 0; i--)
                {
                    double dip2 = dip1;
                    dip1 = di;
                    di = 2.0 * y * dip1 - dip2 + c[i];
                }
                cf[j] = y * di - dip1 + 0.5 * c[0];
            }
        }
        public static void chebyshev_zeros(int n, double[] x)
        {
            for (int i = 0; i < n; i++)
            {
                int j = i + 1;
                double angle = (2.0 * (n - j) +1.0 )* pi / (2.0 * n);
                x[i] = Math.Cos(angle);
            }
        }
    }

Добавлено через 4 минуты
это для функции заданной аналитически - cos(x)

@Milkho · 21.05.2016, 23:35 **[ТС]**

Спасибо. Сейчас буду разбираться.

Добавлено через 3 часа 3 минуты
В общем, мне не удалось переделать код под интерполяцию по заданым точкам.

@afront · 22.05.2016, 08:40

Мне кажется, что интерполяция Чебышева может проводиться исключительно по точкам, являющихся корнями полинома
Чебышева, см. стр . 21 (внизу)
http://www.apmath.spbu.ru/ru/s... numet6.pdf

Добавлено через 12 минут
Для функции заданной таблично этот метод не подходит

@Milkho · 22.05.2016, 09:31 **[ТС]**

afront, да, интерполяционный полином должен быть разложен по n точкам - набору корней полинома n-го порядка.
Но где эти корни использовать и как по них строить интерполянт я не могу понять.

https://youtu.be/KpCJutTQ4WE?t=1h6m10s

@afront · 22.05.2016, 09:47

C#
1
chebyshev_zeros(int n, double[] x)

находит нули для заданного N

C#
1
2
3
4
          for ( int i = 0; i < n; i++)
             {
               x[i] = 0.5 * ( a + b ) + x[i] * 0.5 * ( b - a );
             }

определяет нули на интервале A,B

C#
1
chebyshev_coefficients(double a, double b, int n, double[] c)

находит коэффициенты для заданного интервала A,B

C#
1
chebyshev_coefficients(double a, double b, int n, double[] c)

находит значение функции

Добавлено через 4 минуты
вам нужно сразу было написать чему должно быть равно N, А, В и что у вас за фунция

@Milkho · 22.05.2016, 09:58 **[ТС]**

Так ведь функции нет. Значения функции вводятся с клавиатуры А = X.Min(), B = X.Max(), N = X.Length;

Добавлено через 6 минут
У меня есть, допустим, два массива - один со значениями x, другой со значениями y. Стоит задача: пользователь вводит с клавиатуры какое-то значение х, программа должна вывести приблизительный y.

Новые блоги и статьи Все статьи Все блоги /
делаю науч статью по влиянию грибов на сукцессию anaschu 13.03.2026 прикрепляю статью	SDL3 для Desktop (MinGW): Создаём пустое окно с нуля для 2D-графики на SDL3, Си и C++ 8Observer8 10.03.2026 Содержание блога Финальные проекты на Си и на C++: hello-sdl3-c. zip hello-sdl3-cpp. zip Результат:	Установка CMake и MinGW 13.1 для сборки С и C++ приложений из консоли и из Qt Creator в EXE 8Observer8 10.03.2026 Содержание блога MinGW - это коллекция инструментов для сборки приложений в EXE. CMake - это система сборки приложений. Здесь описаны базовые шаги для старта программирования с помощью CMake и. . .	Как дизайн сайта влияет на конверсию: 7 решений, которые реально повышают заявки Neotwalker 08.03.2026 Многие до сих пор воспринимают дизайн сайта как “красивую оболочку”. На практике всё иначе: дизайн напрямую влияет на то, оставит человек заявку или уйдёт через несколько секунд. Даже если у вас. . .
Модульная разработка через nuget packages DevAlt 07.03.2026 Сложившийся в . Net-среде способ разработки чаще всего предполагает монорепозиторий в котором находятся все исходники. При создании нового решения, мы просто добавляем нужные проекты и имеем. . .	Модульный подход на примере F# DevAlt 06.03.2026 В блоге дяди Боба наткнулся на такое определение: В этой книге («Подход, основанный на вариантах использования») Ивар утверждает, что архитектура программного обеспечения — это структуры,. . .	Управление камерой с помощью скрипта OrbitControls.js на Three.js: Вращение, зум и панорамирование 8Observer8 05.03.2026 Содержание блога Финальная демка в браузере работает на Desktop и мобильных браузерах. Итоговый код: orbit-controls-threejs-js. zip. Сканируйте QR-код на мобильном. Вращайте камеру одним пальцем,. . .	SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .

@afront 1498 / 1213 / 821 Регистрация: 29.02.2016 Сообщений: 3,631
	21.05.2016, 15:47
	Пример на c# https://code.msdn.microsoft.co... d-b10563c8 1

@Milkho 4 / 4 / 3 Регистрация: 24.07.2015 Сообщений: 75
	21.05.2016, 16:50 [ТС]
	afront, Cпасибо, но в том коде сложно разобраться. Как я понял, там интерполируется заданная функция. Но что значит массив коэффициентов (ChebyshevCalculateCoefficients) ? Как их использовать для нахождения значения функции при заданном Х? 0

@afront 1498 / 1213 / 821 Регистрация: 29.02.2016 Сообщений: 3,631
	21.05.2016, 17:42
	а какую функцию вам нужно интерполировать вы не написали? 0

@Milkho 4 / 4 / 3 Регистрация: 24.07.2015 Сообщений: 75
	21.05.2016, 17:50 [ТС]
	afront, У меня нет функции. Есть просто дискретный набор пар точек: X и F(X). Их n штук. 0

@afront 1498 / 1213 / 821 Регистрация: 29.02.2016 Сообщений: 3,631
	21.05.2016, 18:06
	хотите на fortrane рабочий код? там есть и С и С++ http://people.sc.fsu.edu/~jbur... yshev.html Добавлено через 7 минут единственно там функция задается аналитически а не таблично может в этом дело, сделайте тест - задайте известную функцию аналитически и таблично и проверьте 1

@Milkho 4 / 4 / 3 Регистрация: 24.07.2015 Сообщений: 75
	21.05.2016, 18:10 [ТС]
	afront, спасибо еще раз. Но и это не совсем то, что мне нужно. Интерполяция там происходит не по заданным точкам. Вместо этого, алгоритм использует "свои" узлы (нули полинома). К тому же мне придется очень долго изучать весь код этой библиотеки. 0

@afront 1498 / 1213 / 821 Регистрация: 29.02.2016 Сообщений: 3,631
	21.05.2016, 18:12
	Milkho, код там элементарный, хотите переведу на c#? с известной функцией? 0

@afront 1498 / 1213 / 821 Регистрация: 29.02.2016 Сообщений: 3,631
	21.05.2016, 18:21
	потребуется время)) 0

@Milkho 4 / 4 / 3 Регистрация: 24.07.2015 Сообщений: 75
	21.05.2016, 18:26 [ТС]
	afront, я подожду, спасибо огромное вам за помощь. 0

@Milkho 4 / 4 / 3 Регистрация: 24.07.2015 Сообщений: 75
	21.05.2016, 23:35 [ТС]
	Спасибо. Сейчас буду разбираться. Добавлено через 3 часа 3 минуты В общем, мне не удалось переделать код под интерполяцию по заданым точкам. 0

@afront 1498 / 1213 / 821 Регистрация: 29.02.2016 Сообщений: 3,631
	22.05.2016, 08:40
	Мне кажется, что интерполяция Чебышева может проводиться исключительно по точкам, являющихся корнями полинома Чебышева, см. стр . 21 (внизу) http://www.apmath.spbu.ru/ru/s... numet6.pdf Добавлено через 12 минут Для функции заданной таблично этот метод не подходит 0

@Milkho 4 / 4 / 3 Регистрация: 24.07.2015 Сообщений: 75
	22.05.2016, 09:31 [ТС]
	afront, да, интерполяционный полином должен быть разложен по n точкам - набору корней полинома n-го порядка. Но где эти корни использовать и как по них строить интерполянт я не могу понять. https://youtu.be/KpCJutTQ4WE?t=1h6m10s Изображения 0

@Milkho 4 / 4 / 3 Регистрация: 24.07.2015 Сообщений: 75
	22.05.2016, 09:58 [ТС]
	Так ведь функции нет. Значения функции вводятся с клавиатуры А = X.Min(), B = X.Max(), N = X.Length; Добавлено через 6 минут У меня есть, допустим, два массива - один со значениями x, другой со значениями y. Стоит задача: пользователь вводит с клавиатуры какое-то значение х, программа должна вывести приблизительный y. 0