Методы решения обыкновенных дифференциальных уравнений

Таблица 4. Правило Рунге для оценки погрешности.

,

.

За оценку погрешности решения, вычисленного с шагом , принимаем величину 0.113.

Метод прогонки

Дано:

,

,

,

,

,

.

Найти численное решение дифференциального уравнения.

Решение:

Необходимые формулы:

, ,

,

,

,

,

,

,

,

,

,

.

Находим :

,

,

,

.

Исходя из полученных данных, можно найти и все остальное.

Результаты соответствующих итераций приведены в таблице 5:

Таблица 5. Метод прогонки.

Реализация на языках программирования

Таблица идентификаторов

Реализация на С++

Метод Эйлера

Файл eiler.cpp

#include<iostream>

#include<cmath>

#include<iomanip>

using namespace std;

double function(double x, double y){

              return (1+y)/(1+pow(x,3)*y);

}

double x0, xN, h, y;

double dy = 0.0;

int main(){

  cout << "Method eilera" << endl;

  cout << "y'= exp(x*y)" << endl;

  cout << "Vvedite x0: " ;

  cin >> x0;

  cout << "Vvedite xN: " ;

  cin >> xN;

  cout << "Vvedite shag: " ;

  cin >> h;

  cout << "Vvedite y(0): " ;

  cin >> y;

  cout << "Result:" << endl;

  cout << setw(3) << "x" << setw(10) << "y"<< setw(10) << "dy" << endl;

  while(x0 < xN+h){

    dy=h*function(x0,y);

    cout << setw(3) <<  x0 << setw(10) << y << setw(10) << dy << endl; 

    y+=dy;

    x0+=h;

  }

  return 0;

}

Результат работы программы:

Метод Рунге-Кутты 4-го порядка

Файл RK4.cpp

#include<iostream>

#include<cmath>

using namespace std;

double function(double x, double y){

  return exp(x*y);

}

int main()

{

  double x0,xN,h,y,k1,k2,k3,k4;

  cout << "Vvedite a: ";

  cin >> x0;

  cout << "Vvedite b: ";

  cin >> xN;

  cout << "Vvedite shag: ";

  cin >> h;

  cout << "Vvedite y(0): ";

  cin >> y;

  double x=x0;

  cout << "x\ty" << endl;

  while(x < xN+h)

  {

    cout << x << "\t" << y << endl;

    k1 = h*function(x,y);

    k2 = h*function(x + h/2, y + k1/2);

    k3 = h*function(x + h/2, y + k2/2);

    k4 = h*function(x + h, y + k3);

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

    x+=h;

  }

  return 0;

}

Результат работы программы:

Метод Рунге-Кутты 5-го порядка

Файл RK5.cpp

#include<iostream>

#include<cmath>

using namespace std;

double function(double x,double y)

{return exp(x*y);}

int main()

{

  double x0,xN,h,y,k1,k2,k3,k4,k5;

  cout << "Vvedite x0: ";

  cin >> x0;

  cout << "Vvedite xN: ";

  cin >> xN;

  cout << "Vvedite shag: ";

  cin >> h;

  cout << "Vvedite y(0): ";

  cin >> y;

  double x=x0;

  cout << "x\ty"<< endl;

  while(x < xN+h)

  {

    cout << x << "\t" << y << endl;

    k1 = h*function(x,y)/3;

    k2 = h*function(x + h/3, y + k1)/3;

    k3 = h*function(x + h/3, y + k1/2 + k2/2)/3;

    k4 = h*function(x + h/2, y + 3*k1/8 + 9*k3/8)/3;

    k5 = h*function(x + h, y + 3*k1/2 - 9*k3/2 + 6*k4)/3;

    y+=(k1 + 4*k4 + k5)/2;

    x+=h;

  }

  return 0;

}

Результат работы программы:

Метод Адамса

Файл adams.cpp

#include<iostream>

#include<cmath>

using namespace std;

double function(double x,double y)

{return exp(x*y);}

int main()

{

  double x0, xN, h, k1, k2, k3, k4, k5, y, y0, x;

  cout << "Vvedite x0: "; cin >> x0;

  cout << "Vvedite xN: "; cin >> xN;

  cout << "Vvedite shag: "; cin >> h;

  cout << "Vvedite y(" << x0 << "): "; cin >> y0;

  int n = ceil(fabs(xN-x0)/h);

  double *yh = new double[n];

  double *yh2 = new double[2*n];

  y = y0; x = x0;

  for(int i = 0; x <= xN; i++){

    yh[i] = y;

    k1 = h*function(x,y);

    k2 = h*function(x + h/2, y + k1/2);

    k3 = h*function(x + h/2, y + k2/2);

    k4 = h*function(x + h, y + k3);

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

    x += h;

  }

  int i = 4;

  x = x0 + i*h;

  while(x <= xN){

yh[i] = yh[i-1] + h*(55*function(x,yh[i-1]) - 59*function(x-h,yh[i-2]) + 37*function(x-2*h,yh[i-3])-9*function(x-3*h,yh[i-4]))/24;

    x = x+h;

    i++;

  }

  cout << "Result: " << endl;

  cout << "y(" << xN << ") = " << yh[n] << endl;

  y = y0; x = x0; h = h/2;

  for(int i = 0; x <= xN; i++){

    yh2[i] = y;

    k1 = h*function(x,y);

    k2 = h*function(x + h/2, y + k1/2);

    k3 = h*function(x + h/2, y + k2/2);

    k4 = h*function(x + h, y + k3);

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

    x += h;

  }

  i = 8;

  x = x0 + i*h;

  while(x <= xN){

    yh2[i] = yh2[i-1] + h*(55*function(x,yh2[i-1])-59*function(x-h,yh2[i-2]) + 37*function(x-2*h,yh2[i-3])-9*function(x-3*h,yh2[i-4]))/24;

    x = x+h;

    i++;

  }

  cout << "eps = " << (yh[n] - yh2[2*n])/(pow(2, 4)-1) << endl;

  return 0;

}

Результат работы программы:

где eps – погрешность вычисленного решения.

Метод прогонки

Файл progon.cpp

#include <iostream>

#include <cmath>

using namespace std;

const double pi = 3.1416;

double funA(double x)

{

              return 0.4+2*x+x*x/2;

}

double funB(double x)

{

              return pow(x,0.7);

double funC(double x)

{

              return -0.95*x*x;

}

double funF(double x)

{

              return (pow(x,0.5) + cos(2*x/pi))/funA(x);

}

int main()

{

              double a, b, x, h;

              int N;

              cout << "Vvedite a = "; cin >> a;

              cout << "Vvedite b = "; cin >> b;

              if (a>=b){ cout << "bad interval" << endl; return -1; }

              cout << "Vvedite N = "; cin >> N;

              double y[N+1] ,c[N+1], d[N+1], m[N+1], n[N+1];

              cout << "Vvedite y(" << a << ") = "; cin >> y[0];

              cout << "Vvedite y(" << b << ") = "; cin >> y[N];

              h = (b-a)/N;

              m[0]=-2+funB(a)*h/funA(a);

              n[0]=1-funB(a)*h/funA(a)+funC(a)*h*h/funA(a);

              c[0]=1/m[0];

              d[0]=-n[0]*y[0]+funF(a)*h*h;

              for (int i=1;i<N+1;i++)

              {

                            x=a+i*h;

                            m[i]=-2+funB(x)*h/funA(x);

                            n[i]=1-funB(x)*h/funA(x)+funC(x)*h*h/funA(x);

                            c[i]=1/(m[i]-n[i]*c[i-1]);

                            d[i]=funF(x)*h*h-n[i]*c[i-1]*d[i-1];

              }

              for (int i=N;i>1; i--)

              {

                            y[i-1]=c[i-2]*(d[i-2]-y[i]);

              }

              cout << "x\ty" << endl;

              for (int i=0;i<N+1;i++)

              {

                cout << a+i*h << '\t' <<y[i] << endl;

              }

    return 0;

}

Результат работы программы:

Реализация на Fortran

Метод Эйлера

Файл eiler.f

Результат работы программы:

Метод Рунге-Кутты 4-го порядка

Файл RK4.f

Результат работы программы:

Метод Рунге-Кутты 5-го порядка

Файл RK5.f