Пример#16

Пример#16 Интегрирование системы дифференциальных уравнений

[integro.h]

#include<iostream>

#include<vector>

using namespace std;

typedef void (*Tracer)(double t, const vector<double> & x);

typedef vector<double> (*RightPart)(const vector<double> & x, double t);

typedef vector<double> (*StepMaker)(RightPart f, const vector<double> & x, double t, double h);

vector<double> F(const vector<double> & x, double t);

vector<double> Integration (RightPart f, const vector<double> & xo, double t, double tk, double h, StepMaker sm, Tracer P = 0);

vector<double> Eiler (RightPart f, const vector<double> & x, double t, double h);

vector<double> Eiler2 (RightPart f, const vector<double> & x, double t, double h);

vector<double> Hoin (RightPart f, const vector<double> & x, double t, double h);

vector<double> M3 (RightPart f, const vector<double> & x, double t, double h);

vector<double> M4 (RightPart f, const vector<double> & x, double t, double h);

vector<double> M5 (RightPart f, const vector<double> & x, double t, double h);

vector<double> RK4 (RightPart f, const vector<double> & x, double t, double h);

template < class T >

inline std::ostream& operator << (std::ostream& os, const std::vector<T>& v)

{

    for (std::vector<T>::const_iterator ii = v.begin(); ii != v.end(); ++ii)

    {

        os<<*ii ;

    }

    return os;

}

template < class T >

inline vector<T> operator + (const vector<T> &l, const vector<T> &r)

{

    vector<double> res = l;

    for (int i=0;i<res.size();i++)

        res[i]=l[i]+r[i];

    return res;

}

template < class T >

inline vector<T> operator - (const vector<T> &l, const vector<T> &r)

{

    vector<double> res = l;

    for (int i=0;i<res.size();i++)

        res[i]=l[i]-r[i];

    return res;

}

template < class T >

inline vector<T> operator * (double k, vector<T> &r)

{

    for (int i=0;i<r.size();i++)

        r[i]*=k;

    return r;

}

template < class T >

inline vector<T> operator / ( vector<T> &r, double k)

{

    for (int i=0;i<r.size();i++)

        r[i]/=k;

    return r;

}

[integro.cpp]

#include "integro.h"

#include<vector>

#include<iomanip>

vector<double> Integration (RightPart f, const vector<double> & xo, double t, double tk, double h, StepMaker sm, Tracer P){

    vector<double> x = xo;

    while (t <= tk){

        if (P!=0) {

            P(t,x);

            cout<<setiosflags(ios::left | ios::fixed)<<setw(16)<<t<<setw(16)<<setprecision(3)<<x<<endl;

        }

        x = sm(f, x, t, h);

        if (t==tk) break;

        if( t+h >= tk ) h = tk - t;

        t += h;

    }

    return x;

}

//грубый метод Эйлера

vector<double> Eiler (RightPart f, const vector<double> & x, double t, double h) {

    return x + h*f(x, t);

}

//метод Эйлера усовершенствованный

vector<double> Eiler2 (RightPart f, const vector<double> & x, double t, double h){

    vector<double> k1 = h/2.0 * f(x, t);

    return x + h*f(k1, t+h/2.0);

}

//метод второго порядка точости - коректирующий значение

vector<double> Hoin (RightPart f, const vector<double> & x, double t, double h){

    vector<double> k1 = h * f(x, t);

    vector<double> k2 = h * f(x+k1, t+h);

    return x + (k1 + k2)/2.0;

}

vector<double> M3 (RightPart f, const vector<double> & x, double t, double h){

    vector<double> k1 = h * f(x, t);

    vector<double> k2 = h * f(x+k1/2.0, t+h/2.0);

    return x + k2;

}

vector<double> M4 (RightPart f, const vector<double> & x, double t, double h){

    vector<double> k1 = h * f(x,t);

    vector<double> k2 = h * f(x+2*k1/3.0, t+2*h/3.0);

    return x + (k1+3*k2)/4.0;

}

vector<double> M5 (RightPart f, const vector<double> & x, double t, double h){

    vector<double> k1 = h * f(x, t);

    vector<double> k2 = h * f(x+k1/2.0, t+h/2.0);

    vector<double> k3 = h * f(x-k1+2*k2, t+h);

    return x + (k1+4*k2+k3)/6.0;

}

vector<double> RK4 (RightPart f, const vector<double> & x, double t, double h){

    vector<double> k1 = h * f(x, t);

    vector<double> k2 = h * f(x+k1/2.0, t+h/2.0);

    vector<double> k3 = h * f(x+k2/2.0, t+h/2.0);

    vector<double> k4 = h * f(x+k3, t+h);

    return x + (k1+2*k2+2*k3+k4)/6.0;

}

[main.cpp]

#define _USE_MATH_DEFINES

#include<iostream>

#include<conio.h>

#include<fstream>

#include<vector>

#include<cmath>

#include <sys/timeb.h>

#include "integro.h"

#include "GL/glut.h"

#define xmax  180

#define xmin -180

#define ymax  90

#define ymin -90

vector<double> x1;

vector<double> x2;

using namespace std;

void fout (double t, const vector<double> & x)

{

    static ofstream File_Out1 ("result.txt");

    File_Out1 << t<<'\t'<< x;

    x1.push_back(x[0]);

    x2.push_back(x[1]);

}

int GetMilliCount()

{

    timeb tb;

    ftime( &tb );

    int nCount = tb.millitm + (tb.time & 0xfffff) * 1000;

    return nCount;

}

int GetMilliSpan( int nTimeStart )

{

    int nSpan = GetMilliCount() - nTimeStart;

    if ( nSpan < 0 )

        nSpan += 0x100000 * 1000;

    return nSpan;

}

vector<double> F(const vector<double>  & x, double t)

{

    vector<double> F(x.size());

    F[0] = x[1]+t;

    F[1] = sin(x[1]+x[0])+t/2;

    return  F;

}

void draw()

{

    // очистка экрана цветом из буфера

    glClear(GL_COLOR_BUFFER_BIT);

    glLineWidth(1);

    //вертикальная сетка

    for (int y=xmin;y<=xmax;y+=5)

    { 

        glBegin(GL_LINES);

        glColor3f(0, 1, 0);

        glVertex2f(y, ymax);

        glVertex2f(y, ymin);

        glEnd();

    }

    glBegin(GL_LINES);

    glColor3f(0, 0, 0);

    glVertex2f(0, ymax);

    glVertex2f(0, ymin);

    glEnd();

    //горизонтальная сетка

    for (int x=ymin;x<=ymax;x+=5)

    { 

        glBegin(GL_LINES);

        glColor3f(0, 1, 0);

        glVertex2f(xmax, x);

        glVertex2f(xmin, x);

        glEnd();

    }

    glBegin(GL_LINES);

    glColor3f(0, 0, 0);

    glVertex2f(xmax, 0);

    glVertex2f(xmin, 0);

    glEnd();

    glColor3f(1, 0, 0);

    glPointSize(2);

    float tt=0;

    glBegin(GL_LINE_STRIP);

    for(int i=0;i<x1.size()-1;i+=2)

    {

        glVertex2f( tt,     x1[i]   );

        glVertex2f( tt+0.1, x1[i+1] );

        tt+=0.2;

    }

    glEnd();

    glColor3f(0, 0, 1);

    tt=0;

    glBegin(GL_LINE_STRIP);

    for(int i=0;i<x2.size()-1;i+=2)

    {

        glVertex2f( tt,     x2[i]   );

        glVertex2f( tt+0.1, x2[i+1] );

        tt+=0.2;

    }   

    glEnd();

    glColor3f(0.4, 0.1, 0.5);

    glLineWidth(2);

    glBegin(GL_LINE_STRIP);   

    for(int i=0;i<x2.size();i++)

    {

        glVertex2f( x2[i], x1[i] );

    }   

    glEnd();

    glutSwapBuffers();

}

int main(int argc, char **argv)

{

    setlocale(0,"");

    vector <double> x0(2);

    x0[0]=1;

    x0[1]=0;

    vector <double> yo(1);

    yo[0]=-1;

    cout<<"НУ: "<<x0;

    double tk=100;

    //cout<<"Введите tk = ";

    //cin>>tk;

    double h=0.05;

    //cout<<"Введите шаг интегрирования h = ";

    //cin>>h;

    //cout<<"Для запуска расчета нажмите любую клавишу";

    //getch();

    //system("cls");

    int nTimeStart = GetMilliCount();

    glutInit(&argc, argv);

    glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGBA |GLUT_DEPTH);//( GLUT_SINGLE | GLUT_RGB);

    glutInitWindowPosition(glutGet (GLUT_SCREEN_WIDTH)/4, glutGet (GLUT_SCREEN_HEIGHT)/4);

    glutInitWindowSize(800,600);

    glutCreateWindow("Решение ДУ");

    glClearColor(1.0, 1.0, 1.0, 0.0);

    //ортогональная проекция

    gluOrtho2D(-10, 20, -10, 20);

    vector<double> xk1 = Integration (F, x0, 0, tk, h, Eiler, fout);

    cout<<"\nВремя расчета = "<<GetMilliSpan( nTimeStart )/1000<<" (c)"<<endl;

    glutDisplayFunc(draw);

    glutMainLoop();

    return 0;

}