Версия для печати


Компоненты HIntegrator, HDiffObject (решение систем обыкновенных дифференциальных уравнений)
http://www.delphikingdom.com/asp/viewitem.asp?catalogID=1010

Николай
дата публикации 24-05-2004 19:00

Компоненты HIntegrator, HDiffObject (решение систем обыкновенных дифференциальных уравнений)Два компонента для решения и визуализации решений систем обыкновенных дифференциальных уравнений и функционалов одной независимой переменной (очень краткое описание).

HDiffObject:

Это - формализованное описание системы обыкновенных дифференциальных уравнений или функционала одной независимой переменной. Если посмотреть с другой стороны, то это - формализованное описание любого физического объекта или явления функционалом или системой ДУ. Естественно, что такие объекты могут быть совершенно не похожие друг на друга. Их объединяет лишь то, что все существуют в одно и тоже время. Они могут создаваться и разрушаться, взаимодействовать между собой. Каждый из них может проявлять себя в различных формах(например, изображать себя на экране). Привлекает еще то, что большинство физических явлений описывать системой ДУ гораздо проще, чем функционалами.

Любой HDiffObject можно изобразить следующей схемой:
HDiffObject
Инициализация Вычисление Проверка контакта с другой системой Изображение себя

Каждый из блоков схемы описывается обработчиком события в понятиях компонента Delphi. В блоке "Инициализация" определяется начальное состояние объекта, например, его положение, размеры и т. п., то есть начальные значения переменных системы ДУ. В блоке "Вычисление" реализуется вычисление производных или значений переменных как функций времени. В следующем блоке можно обрабатывать и обсчитывать сразу два объекта(вероятность тройных соударений очень мала). И на конец объект может как-то проявить себя во внешнем мире - вынести на обозрение свое состояние, нарисоваться на экране в виде таблицы, графика, рисунка, или проявиться в виде звука и т.п.

Свойства HDiffObject
  • Alias — Пользовательское имя;
  • Divider — Делитель производных;
  • Integrator — Ссылка на HIntegrator;
  • Name — Имя компонента;
  • SystemSize — Размерность системы.
События HDiffObject
  • OnCalculate — Программа вычисления производных или переменных;
  • OnContactCheck — Проверка контакта с другим объектом;
  • OnCreate — Инициализация переменных;
  • OnDraw — Вывод.

Переменные системы ДУ хранятся в динамических массивах Y и dY - значения переменных и производных по времени, соответственно. Индексация переменных начинается с единицы. Элемент с индексом равным нулю - не используется.

HDiffObject, содержит всего один единственный метод - CreateNew. Это - функция, которая возвращает указатель на вновь созданную копию объекта. Только у копии не установлено свойство интегратор.

Множество объектов живут в строго синхронизированном времени, и это обеспечивает центральный компонент - HIntegrator.

HIntegrator:

Для того, что бы все дифференциальные объекты могли существовать и функционировать в реальном времени каждый из них должен иметь ссылку на интегратор. Интегратор, как правило, один на множество ДиффОбъектов. Его назначение - рассчитать следующее состояние всех объектов на заданный квант времени вперед. Он собирает все объекты, которые ссылаются на него, объединяет их в одну большую систему ДУ и продвигает ее по времени. Он же и вызывает все события, которые программируются в HDiffObject, кроме OnCreate. Соответственно состояние всех объектов одного интегратора изменяются в одном и том же временном измерении. Одним из основных свойств интегратора, как компонента, является - величина кванта времени. Она задается в миллисекундах. Величина кванта времени определяется требованиями визуализации. Чем больше квант времени, тем реже объекты будут давать о себе знать на экране и тем ниже будет нагрузка на сам интегратор. Например, если мы хотим видеть положение шарика на столе каждые 50 миллисекунд(20 раз в секунду), то интегратор с этой задачей справится довольно легко. Метод численного интегрирования выполнит шаг намного быстрее даже для системы большого размера. Но если задать квант времени 1 мс то процессор простаивать не будет совсем, а сам интегратор может не успеть обеспечить такую скорость.

В одном приложении можно использовать больше одного интегратора. Но тогда конечно, придется использовать многозадачность. Сам интегратор хорошо уживается с компонентом HProcess. И если он определяет, что работает не в VCL-потоке, то будет использовать программы синхронизации компонента HProcess. Между интеграторами нет никакой связи, хотя любой ДиффОбъект по своему усмотрению может динамически выбирать свой интегратор. Это дает большие возможности при моделировании очень сложных нелинейных систем.

Свойства HIntegrator

Простой пример( "шар на столе" ):
  1. Шар катается по биллиардному столу без трения, и с идеальным отскоком от стенки.
  2. Такую систему ДУ можно описать двумя уравнениями и структурой постоянных величин характеризующих сам объект :
Type
  Param = Record
    FW, FH : integer;
    VX, VY : real;
    iX, iY : real;
    SH : TShape;
  end;

Примечание:
В данном примере столом будет сама форма, а шаром Shape круглой формы. FW и FH - соответственно ширина и высота формы, VX и VY - скорость шара по осям координат, iX и iY - положение шара на столе.

Все постоянные параметры, кроме размеров формы получаются генератором случайных чисел. И самой длинной программой в этом примере является программа построения случайного шара на столе. Сама система уравнений не представляет ни какой сложности:
Здесь X и Y - координаты шара(положение), Vx и Vy - его скорость по осям координат.

Когда шар ударяется в борт стола, то величина его скорости по соответствующей координате меняет знак на противоположный. Вот так выглядит программы вычисления производных для объекта "ШАР":

Метод численного решения систем ОДУ реализован с автоматическим выбором шага интегрирования и отлично справляется с ситуацией, когда производные резко меняют знак. Если квант времени довольно большой, то можно не увидеть столкновение шара с бортом стола, хотя его положение в следующий момент времени, будет вычислено, верно.

В демонстрационном примере, я постарался выбросить все лишнее, чтобы не усложнять понимание сути. Это, например, программу проверки соударений шаров. Но оставил компонент HProcess для демонстрации мультизадачности. К тому же если выполнять интегратор в VCL - потоке, то форму просто так закрыть не удастся. Еще HProcess выводит на форму время загрузки процессора.

Квант времени интегратора - 20 миллисекунд. Видно, что за 22-е секунды времени работы приложения, интегратор работал всего 160 миллисекунд, а процессор 9 секунд. В эти девять секунд входит время работы операционной системы на переключение подзадач, но основное время процессор затратил на рисование(VCL не отличается реентерабельностью, и шары перерисовываются последовательно).

Примечание:
Собственно, это то, для чего и предназначены компоненты - численным методом определить следующее состояние системы в реальном времени с минимальными вычислительными затратами.

Вообще, это вторая версия компонентов. В первой, была возможность выбора численного метода решения СОДУ. Таких алгоритмов было три. Классический RKF45 переведенный с Фортран-реализации программистов SANDIA LABORATORIES, DP54 в моей реализации и самый крутой из этой тройки - LSODE для "жестких" и "не жестких" систем (если найдутся желающие - могу прислать все). Во второй версии остался один DP54. В основном в этом виноват приведенный выше пример. Так алгоритм RKF45 из-за переоценки ошибок иногда просто замораживал пузырь около бортика, что совсем не вписывалось в реальное время. LSODE, часто выдавал плохие коды завершения, а разбираться не было ни желания, ни времени. В отношении же других задач они работают безупречно.

Еще один пример - "Шатунно-поршневая группа"

В этом примере показано как в программе события OnCalculate компонента HDiffObject вычисляются производные и переменные функционала. Вся программа вычисления занимает четыре строки исходного текста:

Procedure TForm1.HD1Calculate(T: Real; HD: THDiffObject);
begin
  HD.dY[1]:=-HD.Y[2];
  HD.dY[2]:=HD.Y[1];
  HD.Y[3]:=HD.Y[1] + Sqrt(Sqr(LS) - Sqr(HD.Y[2]));
  HD.Y[4]:=HD.Y[3] + PS;
end;
Здесь Y[1], Y[2] - координаты X и Y нижней головки шатуна, Y[3] - координата X верхней головки, Y[4] - координата X юбки поршня. LS - длина шатуна, PS - расстояние от юбки поршня до верхней головки шатуна.

Программа задания начальных значений переменных системы ДУ:
Procedure TForm1.HD1Create(T: Real; HD: THDiffObject);
begin
  Rm:=NG.Width / 2; // Радиус головки шатуна
  R:= KV.Width / 2;
  X:= KV.Left + R; // Нулевая точка системы
  Y:= KV.Top  + R; // координат
  Rt:= Sqrt(Sqr(X - NG.Left - Rm) + Sqr(Y - NG.Top - Rm));
  HD.Y[1]:=0;
  HD.Y[2]:=-Rt; // Радиус коленвала
  HD.Y[3]:=VG.Left + Rm;
  HD.Y[4]:=PR.Left;
  PS:=PR.Left - VG.Left;
  LS:=Sqrt(Sqr(HD.Y[1] - HD.Y[3]) + Sqr(HD.Y[2]));
  SpinEdit1.Value:=200;
end;

Время работы программы показано в миллисекундах это - время работы приложения / время работы процессора на все приложение / время работы интегратора. В дополнение еще можно показать программу визуализации, т. е. обработчик события OnDraw компонента HDiffObject:

Procedure TForm1.HD1Draw(HD: THDiffObject);
begin
  NG.Left:=Round((HD.Y[1] + X) - Rm);
  NG.Top :=Round((HD.Y[2] + Y) - Rm);
  VG.Left:=Round(HD.Y[3] - Rm);
  PR.Left:=Round(HD.Y[4] + Rm);
end;
В этой программе просто изменяются положения всех движущихся частей механизма(Shape).

Компоненты приводятся в исходных текстах. В дополнение, еще и компоненты, предназначенные для параллельного программирования, модули реализации матричных вычислений и решения систем ДУ. Компоненты проверены на D5, D6, D7, модули с незначительными изменениями пришли с 5-го Паскаля. Мне уже задавали вопросы по поводу некоторых методов в модуле матричных операций, которые взяты из классики. Скажу сразу: не в моих правилах выдавать за свое то, что написал не я. Все переписки с Фортрана я стараюсь, по мере возможности, сохранять в нетронутом виде. Программа перевода сделана так, что оставляет Фортран-стиль без изменения. К тому же многие алгоритмы более 80-ти лет не претерпели, практически, ни каких изменений. Две программы, по которым были ко мне вопросы, взяты мной из пакета 56-го года. То, что касается решения систем ДУ, то это алгоритм Дорманда-Принца 5-го 4-го порядка, и довольно свежий, реализован еще на 5-м Паскале десять лет назад. В программе вычисление производных не является необходимым атрибутом. Алгоритм построен так, что абсолютно нет разницы, вычисляете Вы производные переменных, либо сами переменные, как функции времени. В таком виде он используется в моей программе HCalc, предназначенной для инженерных вычислений (решения ДУ, краевых задач, моделирования, аппроксимации, оптимизации, расчета регуляторов и многого другого).

Скачать HIntegrator.zip (48K)

С уважением, HNV.