Математика в исходниках (2 стр)

Автор: DROnik

Интегрирование дифуров

Решение дифференциальных уравнений любого порядка численными методами.

немного теории
Все численные методы могут решать только систему ОДУ (однородных дифференциальных уравнений) .
Так что, если у вас имеется дифур с производной выше первой, то его нужно привести к требуемому виду
путем понижения порядка.

Например, имеем диф уравнение:

X''' = sin(X') + tg(X) - X''

вводим обзначения

X   =  X
X'  =  A
X'' =  B

приводим уравнение к системе ОДУ

 X' = A
 A' = B
 B' = sin(A) + tg(X) - B

для чего
Например для физики.
Как известно,

m*X''(t)= F(t,X'',X',X);

Можно это уравнение решать как в школе для частных случаев, но в общем случае аналитически его решить нельзя.
Приходится прибегать к численным методам.

Исходник
Интегрирование методом Dormand-Prince 5(4)

Использование интегратора

program testIntegrator;

{$APPTYPE CONSOLE}
                                                   
uses
  SysUtils,
  Integrator in 'Integrator.pas';

//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>//
//     тестовая функция -   X = sin(t)           //
//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>//
//                                               //
//  дифур:  X'' = -X   // (sin(t))''= -sin(t)    //
//                                               //
//  приводим к Системе ОДУ   X' =  B             //
//                           B' = -X             //
//                                               //
//  вектор состояния  (0) (1)                    //
//                     X   B                     //
//                                               //
//                                               //
//  н.у.  X(0) = 0   //sin(0)                    //
//        B(0) = 1   //cos(0)                    //
//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>//
procedure Diffur(const Ti: IntegrFloat; const aY: array of IntegrFloat;
                 var dY: array of IntegrFloat);
begin
  //Собственно дифур
  //aY: (X , B )
  //dY: (X', B')
  dy[0]:=   aY[1]; {X' =  B}
  dy[1]:=  -aY[0]; {B' = -X}
end;    

var
 int: TRKDP54;
 aY: array [0..1] of single;
 Cur: array [0..1] of single;

const
 to_the = pi/7;

begin
  //задаем вектор начальных условий
  aY[0] := 0;  //sin(0)
  aY[1] := 1;  //cos(0)

  //создаем интегратор
  int:= TRKDP54.Create;

  //задаем начальные значения
  int.SetIniData(   2,     //aDim  - размерность
                    0,     //aT0   - начальное время
                    0.001, //ah    - шаг интегрирования
                    0.01,  //ahMax - максимальный шаг
                    1e-5,  //aTol  - оцениваемая точность
                    aY,    //aY    - вектор начальных условий
                    Diffur // Дифференциальное уравнение
                 );
  //ah, ahMax, aTol  - отвечают за точность интегрирования
  //К сожалению, в использованном методе интегрирования шаг
  //непостоянен, из-за этого дать однозначную рекомендацию
  //по соотношению производительность / качество я не могу.
  //Нужно в каждом случае подбирать эксперементально.

  //вычисляем вектор состояния
  // к моменту времени to_the+1
  int.RunTo(to_the+1,Cur);

  //следующее вычисление будет от предыдущей точки
  //причем интегрировать можно в любом направлении
  int.RunTo(to_the,Cur);

  //не забываем очистить память
  int.Free;

  //Выводим результат и ошибку....
  writeln( formatFloat('0.00000',Cur[0]) );
  writeln( formatFloat('0.00000',Cur[1]) );
  writeln;
  writeln;
  writeln( 'Error by X : ',  formatFloat('0.00000', abs(Cur[0]-sin(to_the)) ));
  writeln( 'Error by X'': ', formatFloat('0.00000', abs(Cur[1]-cos(to_the)) ));
       
  readln;
end.

Страницы: 1 2 3 4 5 6 Следующая »

19 декабря 2005 (Обновление: 26 дек 2005)

Комментарии [1]