Поиск численных решений ОДУ (методы Эйлера, Эйлера-Коши и Рунге – Кутты ) C++

Метод Эйлера

// Метод Эйлера.cpp
//

#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <conio.h>

typedef double typeOfFloat;

typeOfFloat f( typeOfFloat x, typeOfFloat y );

typeOfFloat integral( typeOfFloat a, typeOfFloat b, typeOfFloat y, int n );

typeOfFloat integralOnPiece( typeOfFloat a, typeOfFloat b, typeOfFloat y, typeOfFloat E );

int N = 0;

int main( void )
{

  typeOfFloat a,b;
  typeOfFloat h;
  typeOfFloat e;
  typeOfFloat y0;

  using namespace std;

  cout << "Input a: ";
  cin >> a;
   
  printf("Input b: ");
  cin >> b;

  printf("Input h: ");
  cin >> h;

  printf("Input e: ");
  cin >> e;

  printf("Input y0: ");
  cin >> y0;

  cout << endl;
  cout.precision(9);

  for( typeOfFloat x = a + h; x <= b; x = x + h)
  {

  y0 = y0 + integralOnPiece( x-h, x, y0, e);

  cout << "y = " << y0 << " N = " << N << endl;
   
  }

  getch();

  return 0;

}

// Нахождение интеграла на участке от a до b с погрешностью E
typeOfFloat integralOnPiece( typeOfFloat a, typeOfFloat b, typeOfFloat y, typeOfFloat E )
{

  typeOfFloat s;
  typeOfFloat prS = 0;
  int k;

  prS = integral(a,b,y,2);
  s = integral(a,b,y,3);

   
  for( k = 4; fabs( prS - s ) > E; k++ )
  {
  prS = s;

  s = integral(a,b,y,k);

  }

  return s;
}

// Найти интеграл функции f(x,y) на [a,b] при разбиении на n отрезков методом левых прямоугольников
typeOfFloat integral( typeOfFloat a, typeOfFloat b, typeOfFloat y, int n )
{

  typeOfFloat x, xb, dx;
  typeOfFloat s = 0;

  N=0;
  dx = (b - a)/n;

  xb = a;

  for ( int i = 0; i < n; i++ )
  {
  x = xb + i*dx;
  s = s + f(x,y)*dx;
  N++;
  }

  return s;

}

// Нахождение f(x,y)
typeOfFloat f( typeOfFloat x, typeOfFloat y )  
{  
  return ( x*x*x + y*y*y );  
}  

<b>Эйлера-Коши</b>

// Метод Эйлера-коши.cpp
//

#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <conio.h>

typedef double typeOfFloat;

typeOfFloat f( typeOfFloat x, typeOfFloat y );

typeOfFloat integral( typeOfFloat a, typeOfFloat b, typeOfFloat y, int n );

typeOfFloat integralOnPiece( typeOfFloat a, typeOfFloat b, typeOfFloat y, typeOfFloat E );

   

int N;

int main( void )
{

  typeOfFloat a,b,y3;
  typeOfFloat h;
  typeOfFloat e;
  typeOfFloat y0;
  typeOfFloat Y;

  using namespace std;

  cout << "Input a: ";
  cin >> a;
   
  printf("Input b: ");
  cin >> b;

  printf("Input h: ");
  cin >> h;

  printf("Input e: ");
  cin >> e;

  printf("Input y0: ");
  cin >> y0;

  cout << endl;
  cout.precision(8);

  for( typeOfFloat x = a + h; x <= b; x = x + h)
  {

  Y = y0 + h*integralOnPiece( x-h, x, y0, e);

  y3 = y0 + h * 0.5 * ( f(x-h,y0) + f(x,Y) );

  y0 = Y;

  cout << "Y = " << y3 << " n = " << N << endl;
   
  }

  getch();

  return 0;

}

// Нахождение интеграла на участке от a до b с погрешностью E
typeOfFloat integralOnPiece( typeOfFloat a, typeOfFloat b, typeOfFloat y, typeOfFloat E )
{

  typeOfFloat s;
  typeOfFloat prS = 0;
  int k;

  prS = integral(a,b,y,2);
  s = integral(a,b,y,3);

   
  for( k = 4; fabs( prS - s ) > E; k++ )
  {
  prS = s;

  s = 0;

  s = integral(a,b,y,k);

  }

  return s;
}

// Найти интеграл функции f(x,y) на [a,b] при разбиении на n отрезков методом трапеций
typeOfFloat integral( typeOfFloat a, typeOfFloat b, typeOfFloat y, int n )
{

  typeOfFloat x,dx;
  typeOfFloat s = 0;
  typeOfFloat a1,b1;

  N=0;

  dx = (b - a)/n;

  a1 = a;
  b1 = a + dx;

  for ( int i = 0; i < n; i++ )
  {
  N++;
  s = s + dx/2*(f(a1,y)+f(b1,y));
  a1+=dx;
  b1+=dx;

  }

  return s;

}

// Нахождение f(x,y)
typeOfFloat f( typeOfFloat x, typeOfFloat y )  
{  
  return ( x*x*x + y*y*y );  
}

Метод Рунге-Кутты
// Метод Рунге-Кутта.cpp
//

#include "stdafx.h"
#include <conio.h>
#include <math.h>
#include <iostream>

typedef double typeOfFloat;

typeOfFloat f( typeOfFloat x, typeOfFloat y );

int _tmain(int argc, _TCHAR* argv[])
{

  using namespace std;

  typeOfFloat a,b,y0,h,e,y,y1,y2,H,y3;
  typeOfFloat k1,k2,k3,k4,x ;
  typeOfFloat k11,k22,k33,k44;

  printf("Input a ");
  cin >> a;

  printf("Input b ");
  cin >> b;

  printf("Input h ");
  cin >> h;

  printf("Input e ");
  cin >> e;

  printf("Input y0 ");
  cin >> y0;
   
  y1 = y2 = y0;

  for( x = a + h ; x < b; x += h)
  {  
  int n = 1;

  k1 = h*( f(x-h,y0) );

  k2 = h*( f( x-h/2, y0+k1/2));

  k3 = h*( f( x-h/2, y0+k2/2));

  k4 = h*( f( x, y0 + k3) );

  y2 = y0 + (k1 + 2*k2 + 2*k3 + k4)/6;  

  while(abs(y1 - y2) > e)
  {
  y1 = y2;

  y2 = y0;

  n++;
   
  H = h/n;
  for( int i = 0; i < n; i++)
  {
  k11 = H*( f(x-h,y2) );
  k22 = H*( f(x-h/2, y2+k11/2 ));  
  k33 = H*( f(x-h/2, y2+k22/2 ));
  k44 = H*( f(x, y2+k33 ));  
  y2 = y2 + (k11 + 2*k22 + 2*k33 + k44)/6;  
  }

  }
  printf("y = %f, n = %d\n", y2,n);
  y0 = y2;
  }

  getch();
  return 0;
}

// Нахождение f(x,y)
typeOfFloat f( typeOfFloat x, typeOfFloat y )  
{  
  return ( x*x*x + y*y*y );  
}</div>