Полезное для программистов:

Фриланс
Новости
Статьи
   
Рубрики:


Решение дифференциальных уравнений

Поиск:
Метод Эйлера

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Eiler method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y=h*f(x,y(i))
//
//    In this example y'=cos(y)+3x, [a;b]=[0;1], x0=0, y0=1.3
//
//////////////////////////////////////////////////////////////////////////////

#include <math.h>
#include <stdio.h>

const float a=0,b=1; // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=1.3; // initial conditions
int M=1;

float f(float x, float y)
{
  return (cos(y)+3*x); // y'=cos(y)+3*x
}

void calculate(int m,float *y)
{
  float x,yi,h;

  h=(b-a)/(num_points*m);
  yi=y0; x=x0;
  for (int i=0;i<num_points;i++)
  {
    for (int k=0;k<m;k++)
    {
      yi+=h*f(x,yi);
      x+=h;
    }
    *(y+i)=yi;
  }
}

void main(void)
{

  float yh[num_points],yh2[num_points];

  calculate(M,yh);
  calculate(2*M,yh2);  // doubled step for better accuracy
  // epsilon is difference between solutions with single
  // and double steps
  printf("XttYHttYH2ttEPSILONn");
  for (int i=0;i<num_points;i++)
    printf("%ft%ft%ft%fn",(x0+((i+1)*(b-a))/num_points),
                                      yh[i],yh2[i],yh2[i]-yh[i]);
}



Метод Адамса

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Adams method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    We use Runge-Kutta method to process first 'num_starting_points'
//    points and then use Adams extrapolation to process other points.
//    In this example y'=x+y, [a;b]=[0;2], x0=0, y0=1
//
//////////////////////////////////////////////////////////////////////////////

#include <math.h>
#include <stdio.h>

const float a=0,b=2;             // bounds of the interval
const int num_points=10,         // number of points to solve
          num_starting_points=4; // number of points to solve with Runge-Kutta method
float x0=0,y0=1;                 // starting conditions

float f(float x, float y)
{
  return x+y;  // y'=x+y
}

// this function realises Runge-Kutta method for n starting points
void calculate(float *y)
{
  float k1,k2,k3,k4,x,yi,h;

  h=(b-a)/num_points;  // step
  yi=y0; x=x0;
  for (int i=0;i<num_starting_points;i++)
  {
    k1=h*f(x,yi);
    k2=h*f(x+h/2,yi+k1/2);
    k3=h*f(x+h/2,yi+k2/2);
    k4=h*f(x+h,yi+k3);
    yi+=(k1+2*k2+2*k3+k4)/6;
    x+=h;
    *(y+i+1)=yi;
  }
}

void main(void)
{

  float y[num_points+1],h;

  y[0]=y0;
  // apply Runge-Kutta method
  calculate(y);
  h=(b-a)/num_points;
  // extrapolating
  for (int i=num_starting_points;i<num_points;i++)
    y[i] = y[i-1]+h/24*(55*f(x0+(i-1)*h,y[i-1])-
             59*f(x0+(i-2)*h,y[i-2])+
             37*f(x0+(i-3)*h,y[i-3])-
             9*f(x0+(i-4)*h,y[i-4]));

  printf("XttYttExact solutionn");
  for (i=0;i<num_points;i++)
    printf("%ft%ft%fn",(x0+i*h),y[i],(2*exp(x0+i*h)-(x0+i*h)-1));
}



Метод Эйлера-Коши

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Eiler-Coshi method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y=h*(f(x,y(i))+f(x,z))/2, z=y(i)+h*(f(x(i),y(i))
//
//    In this example y'=cos(y)+3x, [a;b]=[0;1], x0=0, y0=1.3
//
//////////////////////////////////////////////////////////////////////////////

#include <math.h>
#include <stdio.h>

const float a=0,b=1;     // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=1.3;       // initial conditions
int M=1;

float f(float x, float y)
{
  return (cos(y)+3*x); // y'=cos(y)+3*x
}

void calculate(int m,float *y)
{
  float x,yi,h,z;

  h=(b-a)/(num_points*m);
  yi=y0; x=x0;
  for (int i=0;i<num_points;i++)
  {
    for (int k=0;k<m;k++)
    {
      z=yi+h*f(x,yi);
      yi+=h*(f(x,yi)+f(x,z))/2;
      x+=h;
    }
    *(y+i)=yi;
  }
}

void main(void)
{

  float yh[num_points],yh2[num_points];

  calculate(M,yh);
  calculate(2*M,yh2);  // doubled step for better accuracy
  // epsilon is difference between solutions with single
  // and double steps
  printf("XttYHttYH2ttEPSILONn");
  for (int i=0;i<num_points;i++)
    printf("%ft%ft%ft%fn",(x0+((i+1)*(b-a))/num_points),
                                       yh[i],yh2[i],yh2[i]-yh[i]);
}



Метод Рунге-Кутта

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving differential equations (Runge-Kutta method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given differential equation y'=f(x,y) and starting conditions
//    x=x0, y(x0)=y0. We should find solution of this equation
//    at [a;b] interval.
//    y(i+1) calculated as y(i) + delta y(i)
//    delta y =1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4
//    k1=hf(x,y)
//    k2=hf(x+h/2,y+k1/2)
//    k3=hf(x+h/2,y+k2/2)
//    k4=hf(x+h,y+k3)
//
//    In this example y'=cos(x+y)+x-y, [a;b]=[0;1], x0=0, y0=0
//
//////////////////////////////////////////////////////////////////////////////

#include <math.h>
#include <stdio.h>

const float a=0,b=1;     // bound of the interval
const int num_points=10; // number of point to solve
float x0=0,y0=0;         // initial conditions
int    M=1;

float f(float x, float y)
{
  return (cos(x+y)+x-y); // y'=cos(x+y)+x-y
}

void calculate(int m,float *y)
{
  float k1,k2,k3,k4,x,yi,h;

  h=(b-a)/(num_points*m);
  yi=y0; x=x0;
  for (int i=0;i<num_points;i++)
  {
    for (int k=0;k<m;k++)
    {
      // calculatin coefficients k
      k1=h*f(x,yi);
      k2=h*f(x+h/2,yi+k1/2);
      k3=h*f(x+h/2,yi+k2/2);
      k4=h*f(x+h,yi+k3);
      yi+=(k1+2*k2+2*k3+k4)/6;
      x+=h;
   }
   *(y+i)=yi;
  }
}

void main(void)
{
  float yh[num_points],yh2[num_points];

  calculate(M,yh);
  calculate(2*M,yh2);  // doubled step for better accuracy
  // epsilon is difference between solutions with single
  // and double steps
  printf("XttYHttYH2ttEPSILONn");
  for (int i=0;i<num_points;i++)
    printf("%ft%ft%ft%fn",(x0+((i+1)*(b-a))/num_points),
                                    yh[i],yh2[i],(yh2[i]-yh[i]));
}
Сайт: forum.vingrad.ru






Просмотров: 3659

 

 

Новые статьи:


Популярные:
  1. Как сделать цикличным проигрывание MIDI-файла?
  2. Создание AVI файла из рисунков
  3. Как устройство "отключить в данной конфигурации"?
  4. Kто в данный момент присоединен через Сеть?
  5. Как узнать количество доступной памяти?
  6. Как реализовать в RichEdit разноцветный текст?
  7. Как скрыть свое приложение от ProcessViewer
  8. Как программно нажать/скрыть/показ кнопку "Start"?
  9. Модуль работы с ресурсами в PE файлах
10. Функции вызова диалоговых окон выбора
11. Проверка граматики средствами Word'а из Delphi.
12. Модуль для упрощенного вызова сообщений
13. Функции для записи и чтение своих данных в, ЕХЕ- файле
14. Рекурсивный просмотр директорий
15. Network Traffic Monitor
16. Разные модули
17. Универсальная функция для обращения к любым экспортируем функциям DLL
18. Библиотека от VladS
19. Протектор для UPX'а
20. Еще об ICQ, сообщения по контакт листу?
21. Использование открытых интерфейсов
22. Теория и практика использования RTTI
23. Работа с TApplication
24. Примеры использования Drag and Drop для различных визуальных компонентов
25. Что такое порт? Правила для работы с портами
26. Симфония на клавиатуре
27. Загрузка DLL
28. Исправление автоинкремента
29. Взаимодействие с чужими окнами
30. Проверить дубляжи в столбце


 

 

 
 
На главную