2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


Посмотреть правила форума



Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Интерполяция сплайнами
Сообщение03.06.2016, 07:03 


28/02/15
52
Требуется написать программу для интерполяции точек сплайнами. Написал что-то такое:
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
//STL
#include <iostream>
#include <fstream>
#include <iomanip>
//C
#include <clocale>
#include <cmath>

using namespace std;

double *swype (double *a, double *b, double *c, double *d, int n){
    /**
    * Метод прогонки
    */

    double *x, *A, *B, e;
    x = new double [n];
    A = new double [n-1];
    B = new double [n-1];
     //прямая прогонка
    A[0]=-c[0]/b[0];
    B[0]=d[0]/b[0];
    for (int i=1; i<n-1; i++){
        e=a[i-1]*A[i-1]+b[i];
        B[i]=(d[i]-a[i-1]*B[i-1])/e;
        A[i]=-c[i]/e;
    }
    //обратная прогонка
    x[n-1]=(d[n-1]-a[n-2]*B[n-2])/(b[n-1]+a[n-2]*A[n-2]);
    for (int i=n-2; i>=0; i--) x[i]=A[i]*x[i+1]+B[i];
    delete [] A;
    delete [] B;
    return x;
}

double powr (double x, int n){
    /**
    * Возведение в степень
    */

    if (n==0) return 1;
    double y=1;
    for (int i=1; i<=n; i++) y*=x;
    return y;
}

int main()
{
    /**
    * Интерполяция сплайнами
    */

    setlocale(LC_ALL, "rus");
    double *a, *b, *c, *d, *aa, *bb, *cc, *dd, *x, *y, A=-0.5, B=0.5;
    int N;
    ofstream f1("graph.txt"), f2("plot.plt");
    cout<<"Введите N"<<endl;
    cin>>N;
    //инициализация
    a = new double [N];
    aa = new double [N-2];
    b = new double [N];
    bb = new double [N-1];
    c = new double [N];
    cc = new double [N-2];
    d = new double [N];
    dd = new double [N-1];
    x = new double [N+1];
    y = new double [N+1];
    //вычисления
    double h = (B-A)/(N);
    for (int i=0; i<=N; i++) {
        x[i]= A+i*h;
        if (fabs(x[i])<0.000001) x[i]=0;
        y[i]= powr(x[i],3)-3*x[i]+2+10*powr(x[i],10);
        if (i<N) a[i]=y[i];
        if (i<N-1) {
                bb[i]=4*h;
                dd[i]=3*(y[i+1]-2*y[i]+y[i-1])/h;
        }
        if (i<N-2){
            aa[i]=h;
            cc[i]=h;
        }
    }
    c[0]=0;
    c++;
    double * cd = swype(aa,bb,cc,dd,N-1);
    copy(cd,cd+N-1,c);
    c--;
    delete [] cd;
    d[N-1]=-c[N-1]/(3*h);
    b[N-1]=(y[N]-y[N-1])/h-2*c[N-1]/3*h;
    for (int i=0; i<N-1; i++){
        d[i]=(c[i+1]-c[i])/(3*h);
        b[i]=(y[i+1]-y[i])/h-(2*c[i]+c[i+1])/3*h;
    }
    //вывод на экран и запись в файл
    cout<<0<<" "<<0<<" "<<0<<" "<<0<<" "<<x[0]<<" "<<y[0]<<endl;
    f2<<showpos;
    f1<<x[0]<<" "<<y[0]<<endl;
    f2<<"set terminal pngcairo"<<endl<<"set output \"graph.png\""<<endl<<
    "set xrange ["<<x[0]<<":"<<x[N]<<"]"<<endl<<"f(x)=";
    for (int i=1; i<=N; i++){
        f1<<x[i]<<" "<<y[i]<<endl;
        cout<<a[i-1]<<" "<<b[i-1]<<" "<<c[i-1]<<" "<<d[i-1]<<" "<<x[i]<<" "<<y[i]<<endl;
        f2<<"x>="<<x[i-1]<<"&&x<"<<x[i]<<" ? "<<a[i-1]<<b[i-1]<<"*(x"<<-x[i-1]<<")"
        <<c[i-1]<<"*(x"<<-x[i-1]<<")**2"<<d[i-1]<<"*(x"<<-x[i-1]<<")**3 : ";
    }
    f2<<"1/0"<<endl<<"plot f(x), x**3-3*x+2+10*x**10 w l";
    delete [] a;
    delete [] aa;
    delete [] b;
    delete [] bb;
    delete [] c;
    delete [] cc;
    delete [] d;
    delete [] dd;
    delete [] x;
    delete [] y;
    return 0;
}

По методичке вроде бы всё правильно. Но как только я начну интерполировать набор точек
Используется синтаксис Text
-0.5 3.38477
-0.4 3.13705
-0.3 2.87306
-0.2 2.592
-0.1 2.299
0 2
0.1 1.701
0.2 1.408
0.3 1.12706
0.4 0.865049
0.5 0.634766
 

, сплайн получается каким-то некрасивым:
Изображение
Почему это так и что с этим делать?

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение03.06.2016, 08:59 
Заслуженный участник
Аватара пользователя


01/03/06
13626
Москва
Обычно интерполируют кусками кубических парабол, чтобы обеспечить непрерывность производной. Похоже, вы интерполировали кусками квадратичных функций, вот и вышло некрасиво.

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение03.06.2016, 13:36 


28/02/15
52
Brukvalub в сообщении #1128457 писал(а):
Обычно интерполируют кусками кубических парабол, чтобы обеспечить непрерывность производной. Похоже, вы интерполировали кусками квадратичных функций, вот и вышло некрасиво.

Ну вообще-то там вычисляются 4 коэффициента для, соответственно, кубической функции.

(Оффтоп)

(Я же весь код дал!)
Но выглядит это вот так.

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение03.06.2016, 13:48 
Заслуженный участник


06/07/11
5627
кран.набрать.грамота
byulent в сообщении #1128553 писал(а):
Ну вообще-то там вычисляются 4 коэффициента для, соответственно, кубической функции.
У вас же там икс в десятой степени! Или я что-то не понимаю?
byulent в сообщении #1128449 писал(а):
Используется синтаксис C++
 y[i]= powr(x[i],3)-3*x[i]+2+10*powr(x[i],10);


 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение03.06.2016, 13:54 
Заслуженный участник
Аватара пользователя


23/08/07
5501
Нов-ск
rockclimber в сообщении #1128560 писал(а):
У вас же там икс в десятой степени! Или я что-то не понимаю?
Это интерполируемая функция.

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение03.06.2016, 14:12 
Заслуженный участник
Аватара пользователя


01/03/06
13626
Москва
byulent в сообщении #1128553 писал(а):
Ну вообще-то там вычисляются 4 коэффициента для, соответственно, кубической функции.

Приношу свои извинения, я неправильно понял, где исходная функция, а где - сплайн.

(Оффтоп)

(Код я не разбирал, многабукаффниасилю).

Тогда можно предложить для улучшения картинки изменять шаг сплайна, других предложений у меня тет. :oops:

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение03.06.2016, 23:33 


28/02/15
52
Brukvalub в сообщении #1128570 писал(а):
Тогда можно предложить для улучшения картинки изменять шаг сплайна

Лучше не стало.
Изображение
У меня такое чувство, что каждый второй "сплайн" (или, наоборот, каждый первый) должен быть перевёрнут относительно интерполируемой функции, но как это исправить, я не знаю.

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение03.06.2016, 23:41 
Заслуженный участник
Аватара пользователя


01/03/06
13626
Москва
Все-таки, каким цветом на графике изображен сплайн? Если ему отвечает фиолетовый цвет, то косяк где-то в вычислении коэффициентов сплайна, поскольку его первая производная явно разрывна!

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение04.06.2016, 00:56 
Заслуженный участник


20/08/14
11926
Россия, Москва
Я бы предположил что коэффициенты сплайнов вычисляются независимо и сплайны не стыкуются друг с другом, потому все в одну сторону и вылезли. Для всех кроме первого надо брать первую производную от предыдущего куска, а не произвольно.

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение04.06.2016, 01:01 
Заслуженный участник
Аватара пользователя


23/07/08
10910
Crna Gora
Да, ошибка какая-то такая, что непрерывность сплайна обеспечивается, а его первой производной — нет.

Кстати, почему для интерполяции выбран самый неинтересный, почти линейный участок? Стоило ли ради такой унылой зависимости вводить полином десятой степени?

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение04.06.2016, 01:26 
Заслуженный участник


06/07/11
5627
кран.набрать.грамота
Рискну предположить, что коэффициент при третьей степени получается везде нулевым, как следствие - получается "сплайн" из кусков парабол...

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение04.06.2016, 02:16 
Заслуженный участник


20/08/14
11926
Россия, Москва
Ещё подозрение: а почему вообще сплайн так сильно отличается от исходной функции? Уж с такой гладкой функцией он должен вообще практически совпадать. На всех кусочках. Хотя бы первый (который крайний), но вообще говоря все.

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение04.06.2016, 03:15 
Заслуженный участник
Аватара пользователя


23/07/08
10910
Crna Gora
byulent
В функции main в самом первом цикле по i есть оператор
dd[i]=3*(y[i+1]-2*y[i]+y[i-1])/h;
Цикл начинается с i=0, при этом в y[i-1] происходит выход за границы массива: обращение к несуществующему элементу с индексом -1.
Когда Вы это исправите, напишите, как именно исправили, чтобы я синхронизировал свою программу с Вашей.

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение04.06.2016, 06:46 


28/02/15
52
svv в сообщении #1128752 писал(а):
byulent
В функции main в самом первом цикле по i есть оператор
dd[i]=3*(y[i+1]-2*y[i]+y[i-1])/h;
Цикл начинается с i=0, при этом в y[i-1] происходит выход за границы массива: обращение к несуществующему элементу с индексом -1.
Когда Вы это исправите, напишите, как именно исправили, чтобы я синхронизировал свою программу с Вашей.

Исправил:
Используется синтаксис C++
dd[i]=3*(y[i+2]-2*y[i+1]+y[i])/h;

Но это не помогло.
Изображение

 Профиль  
                  
 
 Re: Интерполяция сплайнами
Сообщение04.06.2016, 11:33 
Заслуженный участник
Аватара пользователя


23/07/08
10910
Crna Gora
Ничего, всё равно это надо было исправить.
Ещё одна ошибка в том же операторе (даже исправленном)
dd[i]=3*(y[i+2]-2*y[i+1]+y[i])/h;
Когда значение индекса равно i, для вычисления dd[i] используются значения y[i+1] и y[i+2], но они ещё не определены! Вы определяете их только на следующих итерациях цикла.

Выход: написать два последовательных цикла. В одном инициализировать массивы x и y. А во втором всё остальное, когда x и y уже сформированы.

Вероятно, в этих ошибках и дело: у меня после исправления кривая стала намного красивее.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 23 ]  На страницу 1, 2  След.

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: gris


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group