Исследование и решение задачи Коши для обыкновенных дифференциальных уравнений

double ratio_f (int i,double *r, double *f,double a,double h);

double p_dr (int i,double *r,double *f,double a);

double p_df (int i,double *r,double *f,double a);

double u_dr (int i,double *r,double *f,double a);

double u_df (int i,double *r,double *f,double a);

void graphik (int i,int key,double *r, double *f);

void file_write (int step,double *r,double *f);

 

 

//_________________________________________________//

//_______________Запись в файл____________________//

void file_write (int step,double *r,double *f) {

     ofstream file("Результат.txt");

      int k=0,i=0,b=0;

       while (i<step) {

            

                       k=0;

                        file<<"_______________________________";

                         file<<"_______________________________";

                          file<<"____________________________"<<endl;

                           b=0;

                          

                            while (b<5) {

                                         file<<"|"<<setw(8)<<"r"<<"|"<<setw(8)<<"f";

                                          b++;

                                           }

                                          

                              file<<endl<<"_______________________________";

                               file<<"_______________________________";

                                file<<"___________________________"<<endl;

                                   

                                 while (k<41) {

                                               file<<"|"<<setw(8)<<r[i]<<"|"<<setw(8)<<f[i];

                                                file<<"|"<<setw(8)<<r[i+40]<<"|"<<setw(8)<<f[i+40];

                                                 file<<"|"<<setw(8)<<r[i+80]<<"|"<<setw(8)<<f[i+80];

                                                  file<<"|"<<setw(8)<<r[i+120]<<"|"<<setw(8)<<f[i+120];

                                                   file<<"|"<<setw(8)<<r[i+160]<<"|"<<setw(8)<<f[i+160]<<endl;

                                                    k++;

                                                     i++;

                                                      }

                                   i+=160;           

                                  }

                                  file.close();

                   }

                  

//_________________________________________________//

//_______________График____________________________//

void graphik (int i,int key,double *r, double *f) {

    int m=30,j=0,n=28;          

    setcolor(15);

    setbkcolor(15);

    bar(0,0,1000,200);

    moveto(30,20);

    setcolor(0);

    lineto(30,150);

    lineto(830,150);

    setcolor(2);

    line(5,25,12,25);

    setcolor(12);

    line(5,35,12,35);

    setcolor(0);

    outtextxy(13,16,"r(t)");

    outtextxy(13,30,"f(t)");

    outtextxy(28,151,"t0");

    outtextxy(825,151,"t");

   // setbkcolor(15);

     while (j<key) {

           m+=2;

           n+=2;

            setcolor(2);

            // outtextxy(m,450-r[i],".");

              moveto(n,150-r[i]/10);

               lineto(m,150-r[i+1]/10);

              setcolor(12);

             //  outtextxy(m,450-f[i],".");

                moveto(n,150-f[i]/10);

                 lineto(m,150-f[i+1]/10);

             j++;

             i++;

           }

    

     }

 

//_________________________________________________//

//_______________Функция I____________________//

double func_dr (double r, double f,double a) {

 

       return (2 * r - a * r * f);

        }

 

//_________________________________________________//

//_______________Функция II____________________//       

double func_df (double r, double f,double a) {

     

       return (- f + a * r * f);

        }

       

//_________________________________________________//

//_______________Коэффициенты для I________________//

double ratio_r (int i,double *r, double *f,double a,double h) {

       double k1,k2,k3,k4;

       k1=func_dr (r[i],f[i],a);

        k2=func_dr (r[i]+(h/2),(f[i]+h/2*k1),a);

         k3=func_dr (r[i]+(h/2),(f[i]+h/2*k2),a);

          k4=func_dr ((r[i]+h),(f[i]+h*k3),a);

           return ((h/6)*(k1+(2*k2)+(2*k3)+k4));

       }

 

//_________________________________________________//

//_______________Коэффициенты для II_______________//       

double ratio_f (int i,double *r, double *f,double a,double h) {

       double k1,k2,k3,k4;

       k1=func_df (r[i],f[i],a);

        k2=func_df (r[i]+(h/2),(f[i]+h/2*k1),a);

         k3=func_df (r[i]+(h/2),(f[i]+h/2*k2),a);

          k4=func_df ((r[i]+h),(f[i]+h*k3),a);

           return ((h/6)*(k1+(2*k2)+(2*k3)+k4));

       } 

 

//_________________________________________________//

//_______________ Функция экстраполяционная I ________________//       

double p_dr (int i,double *r,double *f,double a) {

 

       return ((23*func_dr(r[i],f[i],a)-16*func_dr(r[i-1],f[i-1],a)+5*func_dr(r[i-2],f[i-2],a))/12);

      

       }

 

//_________________________________________________//

//_______________Функция экстраполяционная II________________//       

double p_df (int i,double *r,double *f,double a) {

      

       return ((23*func_df(r[i],f[i],a)-16*func_df(r[i-1],f[i-1],a)+5*func_df(r[i-2],f[i-2],a))/12);

      

       }

                    

//_______________ Функция интерполяционная I __________________________________//

//_______________________________//       

double u_dr (int i,double *r,double *f,double a) {

      

       return ((5*func_dr(r[i+1],f[i+1],a)+8*func_dr(r[i],f[i],a)-2*func_dr(r[i-1],f[i-1],a))/12);

      

       }

 

//_________________________________________________//

//_______________Функция интерполяционная II_______________//       

double u_df (int i,double *r,double *f,double a) {

      

       return ((5*func_df(r[i+1],f[i+1],a)+8*func_df(r[i],f[i],a)-2*func_df(r[i-1],f[i-1],a))/12);

      

       }      

                               

 

 

 

//_________________________________________________//

//_______________Основная функция__________________//

int main() {

   

    setlocale(LC_CTYPE, "Russian");

      

    double a=0.01,h=0.1,k1,k2,k3,k4,tmp1,tmp2;

    double r[10000], f[10000],_r,_f,t=0;

     cout<<"//_______________Модель хищник-жертва________________//\n";

      cout<<"//______Ввод начальных данных(решение задачи Коши)___//\n";

       cout<<"В данный момент времени:";

        cout<<"\nR(t0) = "; cin>>r[0];

         cout<<"\nF(t0) = "; cin>>f[0];

          cout<<"\nШаг H = "; cin>>h;

           cout<<"\nКоэффициент а= "; cin>>a;

            cout<<"\nРасчет будет происходить по нажатию Enter (Esc - выход)";

             cout<<"\nЧтобы продолжить нажмите любую клавишу";

              _getch();

             

   //________Первые 2 точки  методом Рунге-Кутта 4 порядка

           int i=0;

           while (i<2) {

                        _r=ratio_r (i,r,f,a,h);

                         _f=ratio_f (i,r,f,a,h); 

                          r[i+1]=r[i]+_r;;

                           f[i+1]=f[i]+_f;

                            i++;

                         }

 

    i=2; 

//_______Метод Адамса 

    int ch,step=0;

    initwindow(1000,200);

     while (ch=13) { 

                    step+=400;   

                     while (i<step) {

                                     r[i+1]=r[i]+h*p_dr(i,r,f,a);

                                      f[i+1]=f[i]+h*p_df(i,r,f,a);

                                       tmp1=r[i]+h*u_dr(i,r,f,a);

                                        tmp2=f[i]+h*u_df(i,r,f,a);

                                         r[i+1]=tmp1;  

                                          f[i+1]=tmp2;

                                           i++;    

                                            }   

     graphik (i-400,400,r,f);

      up: ch=_getch();

       if (ch!=13) {

                    if (ch=27) break;

                     goto up;

                      }

         }

     

  //______Запись               

    file_write (step,r,f);

    system("notepad.exe Результат.txt");

     

     

    cout<<endl;

    system("pause");

    return 0;

    }