RSS    

   Построение кубического сплайна функции - (реферат)

Построение кубического сплайна функции - (реферат)

Дата добавления: март 2006г.

    Построение
    кубического сплайна функции
    План:
    вывод расчётных формул;
    текст программы;
    тестирование.
    Текст программы.
    #include
    #include
    #include
    #include
    #include

#include "mat_vec. h" // классы для работы с матрицами и векторами #include "progonka. h" // решение системы ур-ний (для 3-х диагональных матриц) #include "funct. h" // второстепеннные функции программы (рисование и т. д. )

    // "корень" программы

void spline (float step, int dop, int n, double* &x, double* &y, double* &x1, double* &y1) {

    int k = 0;
    matrica Sp(n, n-1);
    for (int i = 1; i     Sp(i, n) = 3*(y[i-1] - 2*y[i] + y[i+1])/pow(step, 2);
    Sp(i, i) = 4;
    if (i < (n-1)) Sp(i, i+1) = 1;
    if (i > 1) Sp(i, i-1) = 1;
    }
    float *tmp;

progonka(Sp, tmp); // решение системы уравнений методом прогонки // (см. файл "progonka. h")

vector a(n), b(n+1), c(n), d(n); // вычисление коэф-тов многочленов b(1) = 0;

    b(n+1) = 0;
    for(int index = 0; index < n-1; index++)
    b(index+2) = tmp[index];
    delete [] tmp;
    for (i = 1; i     {
    d(i) = y[i-1];
    a(i) = (b(i+1)- b(i))/(-3*step);

c(i) = (y[i] - d(i) - pow(step, 2)*b(i) + pow(step, 3)*a(i) )/(-step); }

    i=0;

//построение графика сплайна при помощи полученный коэф-тов (см. выше) for (i=0; i < n; i++)

    for (int j=0; j < dop; j++)
    {
    x1[k] = x[i] + j*step / (dop);
    y1[k] = pow((x[i]-x1[k]), 3)*a(i+1)
    + pow((x[i]-x1[k]), 2)*b(i+1) + (x[i]-x1[k])*c(i+1)+d(i+1);
    k++;
    }
    x1[n*dop] = x[n];
    y1[n*dop] = y[n];
    }
    void main() {
    int n, dop; double step;
    cout > n;

cout > dop; cout > step;

    dop++;
    double *x, *y, *x1, *y1;
    initial(x, y, x1, y1, n, dop);

int i = 0; while (i < (n+1)) { // расчёт первоначальных значений функции x[i] = (i-n/2)*(step);

    y[i] = cos(x[i])*pow(x[i], 2);
    i++;
    }
    spline (step, dop, n, x, y, x1, y1);
    init(); interface(n, dop, x, y, x1, y1);
    delete x, y, x1, y1;
    closegraph();
    }
    #ifndef __FUNCT_H
    #define __FUNCT_H
    #include
    // инициализация графики
    void init() {
    int D, M; D = DETECT; M = 5;
    initgraph(&D, &M, "");
    }
    // рисование графика функции и сплайна

void paint(int Fx, int Fy, int key, int n, int dop, double* &x, double* &y, double* &x1, double* &y1) {

    int i = 0, a, b;
    a = getmaxx()/2; b = getmaxy()/2;
    setfillstyle(0, 0); bar(0, 0, a*2+1, b*2+1); setcolor(5);
    if ((key == 1) || (key == 3))
    while ( i < n ) {

line(x[i]*Fx + a, -y[i]*Fy + b, x[i+1]*Fx + a, -y[i+1]*Fy + b); i = i++;

    }
    if ((key == 2) || ( key == 3)) {
    i = 0;
    setcolor(3);
    while ( i < n*dop ) {

line(x1[i]*Fx + a, -y1[i]*Fy + b, x1[i+1]*Fx + a, -y1[i+1]*Fy + b); i = i++;

    }
    }
    setcolor(10); line(getmaxx()/2, 0, getmaxx()/2, getmaxy());
    line(0, getmaxy()/2, getmaxx(), getmaxy()/2);
    }

// функция для приближения (удаления) и масштабирования по осям графиков void interface(int n, int dop, double* &x, double* &y, double* &x1, double* &y1) {

    int c=16, z=16;
    char key='0';
    while (key ! = 27) {
    if (key == 75) c = c+4;
    if (key == 72) z = z+4;
    if (key == 77) c = c-4;
    if (key == 80) z = z-4;
    if (key == 45) { z = z-4; c = c-4; }
    if (key == 61) { z = z+4; c = c+4; }
    if (c < 0) c = 0;
    if (z < 0) z = 0;
    if (key == 's') paint(c, z, 2, n, dop, x, y, x1, y1);
    else if (key == 'f') paint(c, z, 1, n, dop, x, y, x1, y1);
    else paint(c, z, 3, n, dop, x, y, x1, y1);
    key = getch();
    }
    }
    // Инициализация динамических массивов

void initial (double* &x, double* &y, double* &x1, double* &y1, int n, int dop) { x = new double [n+1];

    y = new double[n+1];
    for (int i = 0 ; i < (n+1); i++) {
    y[i] = 0;
    x[i] = 0; }
    x1 = new double[n*dop+1];
    y1 = new double[n*dop+1];
    for ( i = 0 ; i < (n*dop+1); i++) {
    x1[i] = 0;
    y1[i] = 0; }
    }
    #endif
    #ifndef __MAT_VEC_H
    #define __MAT_VEC_H
    #include
    #include
    // класс матриц
    class matrica {
    public:

const int Column, String; //кол-во столбцов и строк матрицы matrica(int column, int string);

    ~matrica();
    private:
    float **WW;
    matrica(const matrica& rhs);
    matrica& operator=(const matrica& rhs);
    public:
    float& operator()(int i, int j);

friend ostream& operator>(istream& in, const matrica& matr); };

    // конструктор

matrica : : matrica(int column, int string) : Column(column), String(string) { WW = new float*[string];

    if(! WW) {
    cout     exit(EXIT_FAILURE);
    }
    for(int i = 0; i < string; i++) {
    WW[i] = new float[column];
    if(! WW[i]) {
    cout     exit(EXIT_FAILURE);
    }
    for(int j = 0; j < column; j++)
    WW[i][j] = 0;
    }
    }
    // деструктор
    matrica : : ~matrica() {
    for(int i = 0; i < String; i++)
    delete [] WW[i];
    delete [] WW;
    }
    // операция доступа к элементу
    float& matrica : : operator()(int i, int j) {
    if((i > 0) && (i 0) && (j     return WW[i - 1][j - 1];
    else {

cout
    }
    }
    // вывод матрицы в поток
    ostream& operator    for(int i = 1; i     for(int j = 1; j     out     out     }
    return out     }
    // ввод матрицы из потока
    istream& operator>>(istream& in, matrica& WW) {
    for(int i = 1; i     for(int j = 1; j     in >> WW(i, j);
    return in;
    }
    // класс векторов
    class vector {
    public:
    vector(int column);
    ~vector();
    const int Column; // кол-во элементов вектора
    private:
    float *vect;
    vector(const vector& rhs);
    vector& operator=(const vector& rhs);
    public:
    float& operator()(int i);

friend ostream& operator>(istream& in, const vector& vec); };

    // кнструктор vector
    vector : : vector(int column) : Column(column) {
    vect = new float[column];
    if(! vect) {

cout
    }
    for(int i = 0; i < Column; i++)
    vect[i] = 0;
    }
    // деструктор
    vector : : ~vector() {
    delete [] vect;
    }
    // операция доступа к эелементу
    float& vector : : operator()(int i) {
    if((i > 0) && (i     return vect[i - 1];
    else {

cout
    exit(EXIT_FAILURE);
    }
    }
    // вывод вектора в поток
    ostream& operator     for(int i = 1; i     out     return out     }
    // ввод вектора из потока
    istream& operator>>(istream& in, vector& vec) {
    for(int i = 1; i     in >> vec(i);
    return in;
    }
    #endif
    #ifndef __PROGONKA_H
    #define __PROGONKA_H
    #include "mat_vec. h"
    int progonka(matrica &mat, float* &x) {
    x = new float[mat. String];
    if(! x)
    return 0;

int i, y = mat. Column, n = mat. String; vector h(n), d(n); d(1) = - mat(1, 2) / mat(1, 1);

    h(1) = mat(1, y) / mat(1, 1);
    for(i = 2; i     d(i) = mat(i, i+1) / (mat(i, i-1) * d(i-1) - mat(i, i));

h(i) =(mat(i, y)-mat(i, i-1) * h(i-1))/(mat(i, i-1) * d(i-1) + mat(i, i)); }

h(n) =(mat(n, y)-mat(n, n-1) * h(n-1))/(mat(n, n-1) * d(n-1) + mat(n, n)); x[n-1] = h(n); for ( i=n - 1; i >= 1; i--)

    x[i - 1] = d(i) * x[i] + h(i);
    return 1;
    }
    #endif
    Тестирование:

Зеленым цветом – график функции построенный в пределе от –5 до 5, с шагом = 1. Красным цветом –график сплайна, полученный при интерполировании исходного графика, причём дополнительно построено всего 3 точки на каждом интервале.


Новости


Быстрый поиск

Группа вКонтакте: новости

Пока нет

Новости в Twitter и Facebook

                   

Новости

© 2010.