Календарь
Ноябрь
Пн   5 12 19 26
Вт   6 13 20 27
Ср   7 14 21 28
Чт 1 8 15 22 29
Пт 2 9 16 23 30
Сб 3 10 17 24  
Вс 4 11 18 25  

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



Скачать: Построение кубического сплайна функции

Текст программы.

#include <iostream.h>

#include <fstream.h>

#include <conio.h>

#include <math.h>

#include <dos.h>

#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 <= (n-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 <= n; 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 << "Введите количество интервалов: "; cin >> n;

 cout << "Введите количество доп. т. на интервале: "; cin >> dop;

 cout << "Введите шаг интервала: "; cin >> 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 <graphics.h>

// инициализация графики

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 <stdlib.h>

#include <iostream.h>

// класс матриц

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<<(ostream& out, const matrica& matr);

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

};

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

matrica :: matrica(int column, int string) : Column(column), String(string) {

   WW = new float*[string];

   if(!WW) {

     cout << "\n !!! Не хватает памяти конструктору matrica\n";

     exit(EXIT_FAILURE);

     }

  for(int i = 0; i < string; i++) {

     WW[i] = new float[column];

     if(!WW[i]) {

                   cout << "\n !!! Не хватает памяти конструктору matrica\n";

                   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 <= String) && (j > 0) && (j <= Column))

           return WW[i - 1][j - 1];

    else {

    cout << "\n Ошибка доступа к элементу (" << i << ", " << j << ") ! \n";

    exit(EXIT_FAILURE);

    }

}

// вывод матрицы в поток

ostream& operator<<(ostream& out,  matrica& WW) {

          for(int i = 1; i <= WW.String; i++) {

             for(int j = 1; j <= WW.Column; j++)

                     out << WW(i, j) << " ";

             out << endl;

             }

          return out << "";

}

// ввод матрицы из потока

istream& operator>>(istream& in, matrica& WW) {

          for(int i = 1; i <= WW.String; i++)

             for(int j = 1; j <= WW.Column; 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<<(ostream& out, const vector& vec);

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

};

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

vector :: vector(int column) : Column(column) {

  vect = new float[column];

   if(!vect) {

   cout << endl << "\n !!!Не хватает памяти конструктору vector! \n";

           exit(EXIT_FAILURE);

   }

   for(int i = 0; i < Column; i++)

           vect[i] = 0;

}

// деструктор

vector :: ~vector() {

   delete [] vect;

}

// операция доступа к эелементу

float& vector :: operator()(int i) {

   if((i > 0) && (i <= Column))

                return vect[i - 1];

          else {

            cout << "\n !!!Ошибка доступа к элементу вектора - " << i;

            exit(EXIT_FAILURE);

            }

}

// вывод вектора в поток

ostream& operator << (ostream& out, vector& vec) {

   for(int i = 1; i <= vec.Column; i++)

            out << vec(i) << ' ';

   return out << endl;

}

// ввод вектора из потока

istream& operator>>(istream& in, vector& vec) {

   for(int i = 1; i <= vec.Column; 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 <= n - 1; 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 точки на каждом интервале.



  © Реферат плюс


Поиск

  © REFERATPLUS.RU  

Яндекс.Метрика