Программная реализация модального управления для линейных стационарных систем

Курсовая работа
«Программная реализация модального управления для линейных стационарных систем»

Постановка задачи

1. Для объекта управления с математическим описанием
, (1) — задано,
где — n-мерный вектор состояния, ,
— начальный вектор состояния,
— скалярное управление,
— матрица действительных коэффициентов,
— матрица действительных коэффициентов,
найти управление в функции переменных состояния объекта, т.е.
, (2)
где- матрица обратной связи, такое, чтобы замкнутая система была устойчивой.
2. Корни характеристического уравнения замкнутой системы
(3)
должны выбираться по усмотрению (произвольно) с условием устойчивости системы (3).

Задание

1. Разработать алгоритм решения поставленной задачи.
2. Разработать программу решения поставленной задачи с интерактивным экранным интерфейсом в системах Borland Pascal, Turbo Vision, Delphi — по выбору.
3. Разработать программу решения систем дифференциальных уравнений (1) и (3) с интерактивным экранным интерфейсом.
4. Разработать программу графического построения решений систем (1) и (3) с интерактивным экранным интерфейсом.

Введение
Наряду с общими методами синтеза оптимальных законов управления для стационарных объектов всё большее примене­ние находят методы, основанные на решении задачи о размеще­нии корней характеристического уравнения замкнутой системы в желаемое положение. Этого можно добиться надлежащим выбором матрицы обратной связи по состоянию. Решение ука­занной задачи является предметом теории модального управ­ления (термин связан с тем, что корням характеристического уравнения соответствуют составляющие свободного движения, называемые модами).

Алгоритм модального управления.

Соглашения

Задаваемый объект управления математически описывается уравнением , (1) где и — матрицы действительных коэффициентов, — n-мерный вектор состояния — скалярное управление, — порядок системы (1).
Обратная связь по состоянию имеет вид , (2) где- матрица обратной связи.
Система с введенной обратной связью описывается уравнением (3)
Характеристическое уравнение системы (1) имеет вид (4)
Характеристическое уравнение системы (3) с задаваемыми (желаемыми) корнями имеет вид (5)Алгоритм 1. Для исходной системы (1) составляем матрицу управляемости 2. Обращаем матрицу , т.е. вычисляем . Если не существует (т.е. матрица — вырожденная), то прекращаем вычисления полное управление корнями характеристического уравнения (5) не возможно. 3. Вычисляем матрицу 4. Составляем матрицу 5. Вычисляем матрицу, обратную матрице , т.е. 6. Вычисляем матрицу — матрицу в канонической форме фазовой переменной где — коэффициенты характеристического уравнения (4). Матрица в канонической форме имеет вид 7. Составляем вектор , элементам которого являются коэффициенты характеристического уравнения (4), т.е. , , где — элементы матрицы . 8. Находим коэффициенты характеристического уравнения (5) (см. пояснения) и составляем из них вектор . 9. Вычисляем вектор . — искомая матрица обратной связи системы (3), но она вычислена для системы, матрицы которой заданы в канонической форме фазовой переменной ( и ). 10. Для исходной системы (3) матрица обратной связи получается по формуле Матрица — искомая матрица обратной связи.Пояснения к алгоритму В данной работе рассматривается случай, когда управление единственно и информация о переменных состояния полная. Задача модального управления тогда наиболее просто решается, если уравнения объекта заданы в канонической форме фазовой переменной. Так как управление выбрано в виде линейной функции переменных состояния , где является матрицей строкой . В таком случае уравнение замкнутой системы приобретает вид . Здесь Характеристическое уравнение такой замкнутой системы будет следующим Поскольку каждый коэффициент матрицы обратной связи входит только в один коэффициент характеристического уравнения, то очевидно, что выбором коэффициентов можно получить любые коэффициенты характеристического уравнения, а значит и любое расположение корней. Если же желаемое характеристическое уравнение имеет вид , то коэффициенты матрицы обратной связи вычисляются с помощью соотношений Если при наличии одного управления нормальные уравнения объекта заданы не в канонической форме (что наиболее вероятно), то, в соответствии с пунктами №1-6 алгоритма, от исходной формы с помощью преобразования или нужно перейти к уравнению в указанной канонической форме. Управление возможно, если выполняется условие полной управляемости (ранг матрицы управляемости M должен быть равен n). В алгоритме об управляемости системы судится по существованию матрицы если она существует, то ранг матрицы равен ее порядку (n). Для объекта управления с единственным управлением матрица оказывается также единственной. Для нахождения коэффициентов характеристического уравнения (5), в работе используется соотношения между корнями и коэффициентами линейного алгебраического уравнения степени n , (k = 1, 2, … , n) где многочлены — элементарные симметрические функции, определяемые следующим образом , , , … где Sk — сумма всех произведений, каждое из которых содержит k сомножителей xj с несовпадающими коэффициентами. Программная реализация алгоритма. Текст программной реализации приведен в ПРИЛОЖЕНИИ №1. Вот несколько кратких пояснений.
Программа написана на языке Object Pascal при помощи средств Delphi 2.0, и состоит из следующих основных файлов KursovayaWork.dpr MainUnit.pas SubUnit.pas Matrix.pas Operates.pas HelpUnit.pas OptsUnit.pas
KursovayaWork.dpr — файл проекта, содержащий ссылки на все формы проекта и инициализирующий приложение.
В модуле MainUnit.pas находится описание главной формы приложения, а также сконцентрированы процедуры и функции, поддерживаюшие нужный интерфейс программы.
Модули SubUnit.pas и Operates.pas содержат процедуры и функции, составляющие смысловую часть программной реализации алгоритма, т.е. процедуры решения задачи модально управления, процедуры решения систем дифференциальных уравнений, процедуры отображения графиков решений систем и т.д. Там также находятся процедуры отображения результатов расчетов на экран.
В модуле Matrix.pas расположено описание класса TMatrix — основа матричных данных в программе.
Модули HelpUnit.pas и OptsUnit.pas носят в программе вспомогательный характер.
Для решения систем дифференциальных уравнений использован метод Рунге-Кутта четвертого порядка точности с фиксированным шагом. Метод был позаимствован из пакета программ NumToolBox и адаптирован под новую модель матричных данных.
Обращение матриц производится методом исключения по главным диагональным элементам (метод Гаусса). Этот метод так же был позаимствован из NumToolBox и соответствующе адаптирован.

Пориложение.

program KursovayaWork;

uses
Forms,
MainUnit in ‘MainUnit.pas’ {Form_Main},
OptsUnit in ‘OptsUnit.pas’ {Form_Options},
SubUnit in ‘SubUnit.pas’,
Matrix in ‘Matrix.pas’,
Operates in ‘Operates.pas’,
HelpUnit in ‘HelpUnit.pas’ {Form_Help};

{$R *.RES}

begin
Application.Initialize;
Application.Title = ‘Модальное управление’;
Application.CreateForm(TForm_Main, Form_Main);
Application.CreateForm(TForm_Options, Form_Options);
Application.CreateForm(TForm_Help, Form_Help);
Application.Run;
end.
unit MainUnit;

interface

uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
ComCtrls, Tabnotbk, Menus, StdCtrls, Spin, ExtCtrls, Buttons, Grids,
OleCtrls, VCFImprs, GraphSvr, ChartFX {, ChartFX3};

type
TForm_Main = class(TForm)
BevelMain TBevel;
TabbedNotebook_Main TTabbedNotebook;
SpinEdit_Dim TSpinEdit;
BitBtn_Close TBitBtn;
BitBtn_Compute TBitBtn;
StringGrid_Ap0 TStringGrid;
StringGrid_Anp0 TStringGrid;
StringGrid_Roots TStringGrid;
StringGrid_Kpp0 TStringGrid;
StringGrid_Bp0 TStringGrid;
RadioGroup_RootsType TRadioGroup;
Label_A1p0 TLabel;
Label_Ap0 TLabel;
Label_mBp0 TLabel;
Label_Roots TLabel;
Label_Kpp0 TLabel;
BevelLine TBevel;
Label_Dim TLabel;
StringGrid_Ap1 TStringGrid;
StringGrid_Bp1 TStringGrid;
Label_Ap1 TLabel;
Label_Bp1 TLabel;
StringGrid_Kpp1 TStringGrid;
Label_Kpp1 TLabel;
StringGrid_InCond TStringGrid;
Label_InCond TLabel;
Label_U TLabel;
Edit_U TEdit;
BitBtn_Options TBitBtn;
BitBtn_Help TBitBtn;
StringGrid_ABKpp1 TStringGrid;
Label_ABKpp1 TLabel;
Edit_W TEdit;
Label_w TLabel;
RadioGroupChart TRadioGroup;
ChartFX TChartFX;
LabelW1 TLabel;
StringGrid_Solve1 TStringGrid;
StringGrid_Solve2 TStringGrid;
Label1 TLabel;
Label2 TLabel;
Label3 TLabel;
procedure BitBtn_CloseClick(Sender TObject);
procedure BitBtn_OptionsClick(Sender TObject);
procedure BitBtn_ComputeClick(Sender TObject);
procedure FormCreate(Sender TObject);
procedure SpinEdit_DimChange(Sender TObject);
procedure StringGrid_RootsSetEditText(Sender TObject; ACol,
ARow Longint; const Value string);
procedure RadioGroup_RootsTypeClick(Sender TObject);
procedure TabbedNotebook_MainChange(Sender TObject; NewTab Integer;
var AllowChange Boolean);
procedure StringGrid_SetEditText(Sender TObject; ACol,
ARow Longint; const Value string);
procedure BitBtn_HelpClick(Sender TObject);
procedure RadioGroupChartClick(Sender TObject);
private
procedure FillFixedCellsInAllGrids;
procedure FillCellsInAllGrids;
public
procedure BindGrids;
procedure UnBindGrids;
end;

var
Form_Main TForm_Main;

implementation

uses Matrix, SubUnit, OptsUnit, Operates, CFXOCX2, HelpUnit;

const
DefOptions = [goFixedVertLine, goFixedHorzLine,
goVertLine, goHorzLine,
goColSizing, goEditing,
goAlwaysShowEditor, goThumbTracking];
{$R *.DFM}

procedure TForm_Main.FillFixedCellsInAllGrids;
var
Order TOrder;
i byte;
Str string;
begin
Order = SpinEdit_Dim.Value;
for i = 1 to Order do
begin
Str = IntToStr(i);
StringGrid_Ap0.Cells[0, i] = Str;
StringGrid_Ap0.Cells[i, 0] = Str;
StringGrid_Bp0.Cells[0, i] = Str;
StringGrid_ANp0.Cells[i, 0] = Str;
StringGrid_ANp0.Cells[0, i] = Str;
StringGrid_Roots.Cells[i, 0] = Str;
StringGrid_Kpp0.Cells[i, 0] = Str;
StringGrid_Ap1.Cells[0, i] = Str;
StringGrid_Ap1.Cells[i, 0] = Str;
StringGrid_Bp1.Cells[0, i] = Str;
StringGrid_ABKpp1.Cells[i, 0] = Str;
StringGrid_ABKpp1.Cells[0, i] = Str;
StringGrid_InCond.Cells[i, 0] = Str;
StringGrid_Kpp1.Cells[i, 0] = Str;
StringGrid_Solve1.Cells[i, 0] = ‘X’ + IntToStr(i);
StringGrid_Solve2.Cells[i, 0] = ‘X’ + IntToStr(i);
StringGrid_Solve1.Cells[0, 0] = ‘Время’;
StringGrid_Solve2.Cells[0, 0] = ‘Время’;
end;
end;

procedure TForm_Main.FillCellsInAllGrids;
var
Order TOrder;
i, j byte;
begin
Order = SpinEdit_Dim.Value;
for i = 1 to Order do
for j = 1 to Order do
begin
StringGrid_Ap0.Cells[j, i] = ‘0’;
StringGrid_Ap0.Cells[i, i] = ‘1’;
StringGrid_Bp0.Cells[1, i] = ‘0’;
StringGrid_Roots.Cells[i, 1] = ‘-1’;
StringGrid_Roots.Cells[i, 2] = ‘0’;
StringGrid_Kpp0.Cells[i, 1] = ‘0’;
StringGrid_Ap1.Cells[j, i] = ‘0’;
StringGrid_Ap1.Cells[i, i] = ‘1’;
StringGrid_Bp1.Cells[1, i] = ‘0’;
StringGrid_ABKpp1.Cells[j, i] = ‘0’;
StringGrid_ABKpp1.Cells[i, i] = ‘1’;
StringGrid_InCond.Cells[i, 1] = ‘0’;
StringGrid_Kpp1.Cells[i, 1] = ‘0’;
end;
FillFixedCellsInAllGrids;
StringGrid_Roots.Cells[0, 1] = ‘Re’;
StringGrid_Roots.Cells[0, 2] = ‘Im’;
StringGrid_Bp1.Cells[1, 0] = ‘1’;
StringGrid_Bp0.Cells[1, 0] = ‘1’;
end;

procedure TForm_Main.BindGrids;
begin
CopyGrid(StringGrid_Ap1, StringGrid_Ap0);
CopyGrid(StringGrid_Bp1, StringGrid_Bp0);
CopyGrid(StringGrid_Kpp1, StringGrid_Kpp0);
StringGrid_Ap1.Options = DefOptions — [goEditing];
StringGrid_Bp1.Options = DefOptions — [goEditing];
StringGrid_Kpp1.Options = DefOptions — [goEditing];
end;

procedure TForm_Main.UnBindGrids;
begin
StringGrid_Ap1.Options = DefOptions;
StringGrid_Bp1.Options = DefOptions;
StringGrid_Kpp1.Options = DefOptions;
end;

procedure TForm_Main.BitBtn_CloseClick(Sender TObject);
begin
Close;
end;

procedure TForm_Main.BitBtn_OptionsClick(Sender TObject);
var
V0, V1, V2, V3 LongInt;
LS TCheckBoxState;
begin
with Form_Options do
begin
V0 = SpinEdit0.Value;
V1 = SpinEdit1.Value;
V2 = SpinEdit2.Value;
V3 = SpinEdit3.Value;
LS = CheckBox_Link.State;
ShowModal;
if ModalResult = mrCancel then
begin
SpinEdit0.Value = V0;
SpinEdit1.Value = V1;
SpinEdit2.Value = V2;
SpinEdit3.Value = V3;
CheckBox_Link.State = LS;
end
else
if ((SpinEdit0.Value <> V0) or (SpinEdit1.Value <> V1)) or
((SpinEdit2.Value <> V2) or (SpinEdit3.Value <> V3)) then
begin
BitBtn_Compute.Enabled = True;
case BitBtn_Compute.Tag of
4, 5 BitBtn_Compute.Tag = BitBtn_Compute.Tag — 4;
6, 7 BitBtn_Compute.Tag = BitBtn_Compute.Tag — 4;
8, 9 BitBtn_Compute.Tag = BitBtn_Compute.Tag — 8;
10, 11 BitBtn_Compute.Tag = BitBtn_Compute.Tag — 8;
12, 13 BitBtn_Compute.Tag = BitBtn_Compute.Tag — 12;
14, 15 BitBtn_Compute.Tag = BitBtn_Compute.Tag — 12;
end;
end;
end;
end;

procedure TForm_Main.BitBtn_ComputeClick(Sender TObject);
begin
BitBtn_Compute.Enabled = False;
if Form_Options.CheckBox_Link.State = cbChecked then BindGrids;
case TabbedNotebook_Main.PageIndex of
0 begin
ComputeFromPage0;
BitBtn_Compute.Tag = BitBtn_Compute.Tag + 1;
end;
1 begin
ComputeFromPage1;
ShowChart(Succ(RadioGroupChart.ItemIndex));
BitBtn_Compute.Tag = BitBtn_Compute.Tag + 14;
end;
2 begin
ComputeFromPage2;
BitBtn_Compute.Tag = BitBtn_Compute.Tag + 4;
end;
3 begin
ComputeFromPage3;
BitBtn_Compute.Tag = BitBtn_Compute.Tag + 8;
end;
end;
end;

procedure TForm_Main.FormCreate(Sender TObject);
const
FirstColWidth = 20;
begin
StringGrid_Ap0.ColWidths [0] = FirstColWidth;
StringGrid_Anp0.ColWidths [0] = FirstColWidth;
StringGrid_Bp0.ColWidths [0] = FirstColWidth;
StringGrid_Roots.ColWidths [0] = FirstColWidth;
StringGrid_Ap1.ColWidths [0] = FirstColWidth;
StringGrid_ABKpp1.ColWidths [0] = FirstColWidth;
StringGrid_Bp1.ColWidths [0] = FirstColWidth;
StringGrid_Kpp0.ColWidths [0] = FirstColWidth;
StringGrid_Kpp1.ColWidths [0] = FirstColWidth;
StringGrid_InCond.ColWidths [0] = FirstColWidth;
FillCellsInAllGrids;
BindGrids;
end;

procedure TForm_Main.SpinEdit_DimChange(Sender TObject);
var
Order byte;
begin
Order = Succ(SpinEdit_Dim.Value);
StringGrid_Ap0.ColCount = Order;
StringGrid_Ap0.RowCount = Order;
StringGrid_Anp0.ColCount = Order;
StringGrid_Anp0.RowCount = Order;
StringGrid_Bp0.RowCount = Order;
StringGrid_Roots.ColCount = Order;
StringGrid_Kpp0.ColCount = Order;
StringGrid_Ap1.ColCount = Order;
StringGrid_Ap1.RowCount = Order;
StringGrid_Bp1.RowCount = Order;
StringGrid_ABKpp1.ColCount = Order;
StringGrid_ABKpp1.RowCount = Order;
StringGrid_InCond.ColCount = Order;
StringGrid_Kpp1.ColCount = Order;
FillFixedCellsInAllGrids;
BitBtn_Compute.Enabled = True;
end;

procedure TForm_Main.StringGrid_RootsSetEditText(Sender TObject; ACol,
ARow Longint; const Value string);
var
Val string;
begin
if (ARow = 2) and (Value <> ») then
begin
Val = StringGrid_Roots.Cells [ACol, ARow];
if StrToFloat (Value) <> 0 then
StringGrid_Roots.Cells[Succ(ACol),ARow] =FloatToStr(- StrToFloat(Value));
if StrToFloat (Value) = 0 then
StringGrid_Roots.Cells [Succ(ACol),ARow] = FloatToStr(0);
end;
end;

procedure TForm_Main.RadioGroup_RootsTypeClick(Sender TObject);
var
Order TOrder;
j byte;
NHalf byte;
StartAlfa, NAlfa, dAlfa Float;
W Float;
begin
Order = SpinEdit_Dim.Value;
W = StrToFloat (Edit_W.Text);
case RadioGroup_RootsType.ItemIndex of
0 StringGrid_Roots.Options = DefOptions;
1 begin
for j = 1 to Order do
begin
StringGrid_Roots.Cells [j, 1] = FloatToStr (-W);
StringGrid_Roots.Cells [j, 2] = ‘0’;
StringGrid_Roots.Options = DefOptions — [goEditing];
end
end;
2 begin
dAlfa = Pi / Order;
StartAlfa = Pi/2 — dAlfa/2;
NHalf = Order div 2;
for j = 1 to NHalf do
begin
NAlfa = StartAlfa + dAlfa * j;
StringGrid_Roots.Cells [j, 1] = FloatToStr (Cos (NAlfa) * W);
StringGrid_Roots.Cells [Order — Pred (j), 1] = FloatToStr (Cos (-NAlfa) * W);
StringGrid_Roots.Cells [j, 2] = FloatToStr (Sin (NAlfa) * W);
StringGrid_Roots.Cells [Order — Pred (j), 2] = FloatToStr (Sin (-NAlfa) * W);
end;
if Odd (Order) then
begin
StringGrid_Roots.Cells [NHalf +1, 1] = FloatToStr (-W);
StringGrid_Roots.Cells [NHalf +1, 2] = ‘0’;
end;
StringGrid_Roots.Options = DefOptions — [goEditing];
end;
end;

end;

procedure TForm_Main.TabbedNotebook_MainChange(Sender TObject;
NewTab Integer; var AllowChange Boolean);
begin
with BitBtn_Compute do
case NewTab of
0 begin
SpinEdit_Dim.Enabled = True;
if Tag in [1, 3, 5, 7, 9, 11, 13, 15] then Enabled = False
else Enabled = True;
BitBtn_Compute.Caption = ‘Рассчитать модальное управление’;
end;
1 begin
SpinEdit_Dim.Enabled = True;
if Tag in [2, 3, 6, 7, 10, 11, 14, 15] then Enabled = False
else Enabled = True;
BitBtn_Compute.Caption = ‘Решить системы дифф. уравнений ‘;
if Form_Options.CheckBox_Link.State = cbChecked then BindGrids;
end;
2 begin
SpinEdit_Dim.Enabled = False;
if Tag in [4, 5, 6, 7, 12, 13, 14, 15] then Enabled = False
else Enabled = True;
BitBtn_Compute.Caption = ‘Обновить результаты решений ‘;
end;
3 begin
SpinEdit_Dim.Enabled = False;
if Tag in [8, 9, 10, 11, 12, 13, 14, 15] then Enabled = False
else Enabled = True;
BitBtn_Compute.Caption = ‘Обновить диаграмму решения ‘;
end;
end;
end;

procedure TForm_Main.StringGrid_SetEditText(Sender TObject; ACol,
ARow Longint; const Value string);
begin
if not BitBtn_Compute.Enabled then
case TabbedNotebook_Main.PageIndex of
0 if Form_Options.CheckBox_Link.State = cbChecked then
BitBtn_Compute.Tag = BitBtn_Compute.Tag — 3
else
BitBtn_Compute.Tag = BitBtn_Compute.Tag — 1;
1 BitBtn_Compute.Tag = BitBtn_Compute.Tag — 2;
end;
BitBtn_Compute.Enabled = True;
end;

procedure TForm_Main.BitBtn_HelpClick(Sender TObject);
begin
Form_Help.ShowModal;
end;

procedure TForm_Main.RadioGroupChartClick(Sender TObject);
begin
case RadioGroupChart.ItemIndex of
0 ShowChart(1);
1 ShowChart(2);
end;
end;

end.
unit SubUnit;

interface

uses
SysUtils, Matrix, Operates, Grids;

procedure CopyGrid(AGrid, BGrid TStringGrid);
procedure LoadMatrixSolveFromStrGrd (AMatrix TMatrix; AGrid TStringGrid);
procedure ComputeFromPage0;
procedure ComputeFromPage1;
procedure ComputeFromPage2;
procedure ComputeFromPage3;
procedure ShowChart(NumberOfChart Byte);

implementation

uses
MainUnit, OptsUnit, CFXOCX2;

procedure CopyGrid(AGrid, BGrid TStringGrid);
var
i, j LongInt;
begin
AGrid.ColCount = BGrid.ColCount;
AGrid.RowCount = BGrid.RowCount;
for j = 0 to AGrid.ColCount do
for i = 0 to AGrid.RowCount do
AGrid.Cells[j, i] = BGrid.Cells[j, i];
end;

function CropStr (Str String) String;
var
i Byte;
Str_1 String;
Begin
for i = Length(Str) downto 1 do
if Str [i] = ‘ ‘ then Str = Copy(Str, 1, i-1)
else Break;
Str_1 = Str;
for i = 1 to Length(Str) do
if Str[i] = ‘ ‘ then Str_1 = Copy(Str, i+1, Length(Str) — i)
else Break;
CropStr = Str_1;
End;

procedure LoadMatrixFromStrGrd (AMatrix TMatrix; AGrid TStringGrid);
var
i, j Word;
begin
AMatrix.Resize (Pred(AGrid.ColCount), Pred(AGrid.RowCount));
for i = 1 to AMatrix.RowCount do
for j = 1 to AMatrix.ColCount do
begin
if CropStr(AGrid.Cells[j, i]) = » then AGrid.Cells[j, i] = ‘0’;
AMatrix[j ,i] = StrToFloat(AGrid.Cells[j, i])
end
end;

procedure OutPutMatrixToStrGrd (AMatrix TMatrix; AGrid TStringGrid);
var
i, j Word;
begin
AGrid.ColCount = Succ(AMatrix.ColCount);
AGrid.RowCount = Succ(AMatrix.RowCount);
for i = 1 to AMatrix.RowCount do
for j = 1 to AMatrix.ColCount do
begin
AGrid.Cells[j, 0] = IntToStr (j);
AGrid.Cells[0, i] = IntToStr (i);
AGrid.Cells[j, i] = FloatToStrF(AMatrix[j ,i],ffGeneral,5,3);
end
end;

procedure OutPutMatrixSolveToStrGrd (AMatrix TMatrix; AGrid TStringGrid);
var
i, j, k Word;
begin
AGrid.ColCount = AMatrix.ColCount;
AGrid.RowCount = Succ(AMatrix.RowCount);
for i = 1 to AMatrix.RowCount do
for j = 1 to AMatrix.ColCount do
begin
if j = AMatrix.ColCount then k = 0 else k = j;
AGrid.Cells[j, 0] = ‘X’ + IntToStr (j);
AGrid.Cells[k, i] = FloatToStrF(AMatrix[j ,i],ffGeneral,5,3);
end;
AGrid.Cells[0, 0] = ‘Время’;
end;

procedure LoadMatrixSolveFromStrGrd (AMatrix TMatrix; AGrid TStringGrid);
var
i, j, k Word;
begin
AMatrix.Resize (AGrid.ColCount, Pred(AGrid.RowCount));
for i = 1 to AMatrix.RowCount do
for j = 0 to AMatrix.ColCount do
begin
if j = 0 then k = AMatrix.ColCount else k = j;
if CropStr(AGrid.Cells[j, i]) = » then AGrid.Cells[j, i] = ‘0’;
AMatrix[k ,i] = StrToFloat(AGrid.Cells[j, i])
end
end;

procedure ComputeFromPage0;
var
Order TOrder;
i, j byte;
K ShortInt;
mDummy1, mDummy2,
mA, mB, mKp,
mM, mN, mN1 TMatrix;
cvRoots TComplexVector;
begin
with Form_Main do
begin
Order = SpinEdit_Dim.Value;

mA = TMatrix.Create(Order, Order);
mB = TMatrix.Create(1, Order);
mM = TMatrix.Create(Order, Order);
mDummy1 = TMatrix.Create(Order, Order);
mN1 = TMatrix.Create(Order, 1);
mN = TMatrix.Create(Order, Order);
mDummy2 = TMatrix.Create(Order, Order);
mKp = TMatrix.Create(Order, 1);

LoadMatrixFromStrGrd (mA, StringGrid_Ap0);
LoadMatrixFromStrGrd (mB, StringGrid_Bp0);

for j = 1 to Order do
begin
mDummy1.Assign(mA);
mDummy1.NthPower(j — 1);
mDummy1.MultFromRight(mB);
for i = 1 to Order do
mM[j, i] = mDummy1[1, i];
end;

if not mM.Inverse then
Raise ESingularMatrix.Create(‘Система неполностью управляема ‘ +
‘матрица M — вырожденная !!!’#10 +
‘Измените значения коэффициентов матриц А и B’);

mN1.SetNull;
mN1[Order, 1] = 1;
mN1.MultFromRight(mM);

for i = 1 to Order do
begin
mDummy2.Assign(mA);
mDummy2.NthPower(i-1);
mDummy1.Assign(mN1);
mDummy1.MultFromRight(mDummy2);
for j = 1 to Order do mN[j, i] = mDummy1[j, 1];
end;

mDummy1.Assign(mN);
if not mDummy1.Inverse then
Raise ESingularMatrix.Create(‘Не могу обратить матрицу N !!!’#10 +
‘(не разбрасывайтесь порядками коэффициентов матриц)’);
mA.MultFromLeft(mN);
mA.MultFromRight(mDummy1);
OutPutMatrixToStrGrd(mA, StringGrid_Anp0);

cvRoots.Dim = Order;
for j = 1 to Order do
begin
cvRoots.Data[j].Re = StrToFloat(StringGrid_Roots.Cells[j, 1]);
cvRoots.Data[j].Im = StrToFloat(StringGrid_Roots.Cells[j, 2]);
end;

for j = 1 to Order do
begin
if Odd (j) then K = -1 else K = +1;
mKp[Order-Pred(j), 1] = — mA[Order-Pred(j), Order] —
K * SymmetricalFunction(cvRoots, j);
end;
mKp.MultFromRight(mN);
OutPutMatrixToStrGrd (mKp, StringGrid_Kpp0);

mDummy1.Free;
mDummy2.Free;
mA.Free;
mB.Free;
mKp.Free;
mM.Free;
mN.Free;
mN1.Free;
end;
end;

procedure ComputeFromPage1;
var
Order TOrder;
mA, mB, mABKp, mInCond, mKp TMatrix;
mSolutionValues TMatrix;
LowerLimit, UpperLimit, NumReturn, NumIntervals Word;
begin
with Form_Main do
begin
Order = SpinEdit_Dim.Value;

mA = TMatrix.Create(Order, Order);
mB = TMatrix.Create(1, Order);
mKp = TMatrix.Create(Order, 1);
mInCond = TMatrix.Create(Order, 1);

LoadMatrixFromStrGrd(mA, StringGrid_Ap1);
LoadMatrixFromStrGrd(mB, StringGrid_Bp1);
LoadMatrixFromStrGrd(mKp, StringGrid_Kpp1);
LoadMatrixFromStrGrd(mInCond, StringGrid_InCond);

mABKp = TMatrix.Create(Order, Order);
mABKp.Assign(mB);
mABKp.MultFromRight(mKp);
mABKp.AddMatrix(mA);

OutPutMatrixToStrGrd(mABKp, StringGrid_ABKpp1);

mB.MultConst(StrToFloat(Edit_U.Text));

with Form_Options do
begin
LowerLimit = SpinEdit0.Value;
UpperLimit = SpinEdit1.Value;
NumReturn = SpinEdit2.Value;
NumIntervals = SpinEdit3.Value;
end;

mSolutionValues = TMatrix.Create(1, 1);

try
DiffSystemSolve (mA, mB,
LowerLimit, UpperLimit,
mInCond,
NumReturn, NumIntervals,
mSolutionValues);
OutPutMatrixSolveToStrGrd(mSolutionValues, StringGrid_Solve1);

mSolutionValues.ReSize(1, 1);
DiffSystemSolve (mABKp, mB,
LowerLimit, UpperLimit,
mInCond,
NumReturn, NumIntervals,
mSolutionValues);
OutPutMatrixSolveToStrGrd(mSolutionValues, StringGrid_Solve2);

except
on EO EOverflow do
begin
EO.Message = ‘Не буду считать !!!’#10 +
‘С уменьшите разброс коэффициентов в матрицах’#10 +
‘либо измените опции (уменьшите их pls.)’;
Raise;
end;
end;

mA.Free;
mB.Free;
mABKp.Free;
mInCond.Free;
mKp.Free;
mSolutionValues.Free;
end;
end;

procedure ShowChart(NumberOfChart Byte);
var
Order, Serie TOrder;
NumReturn, Point Word;
mSolutionValues TMatrix;

procedure SetAdm;
const
Divisor = 3.4E+38;
var
i, j LongInt;
Greatest, Least Float;
begin
Greatest = mSolutionValues[1, 1];
Least = Greatest;
for j = 1 to Order do
for i = 1 to NumReturn do
begin
if mSolutionValues[j, i] > Greatest then Greatest = mSolutionValues[j, i];
if mSolutionValues[j, i] < Least then Least = mSolutionValues[j, i];
end;
Form_Main.ChartFX.Adm[CSA_MAX] = Greatest;
Form_Main.ChartFX.Adm[CSA_MIN] = Least;
Form_Main.ChartFX.Title[CHART_TOPTIT] = ‘Y = Y » * ‘;
end;

begin
with Form_Main do
begin
Order = SpinEdit_Dim.Value;
NumReturn = Form_Options.SpinEdit2.Value;
mSolutionValues = TMatrix.Create(1, 1);
ComputeFromPage1;
case NumberOfChart of
1 begin
LoadMatrixSolveFromStrGrd(mSolutionValues, StringGrid_Solve1);
SetAdm;
ChartFX.OpenDataEx(Cod_Values, Order, Pred(NumReturn));
for Serie = 1 to Order do
begin
ChartFX.SerLeg[Pred(Serie)] = ‘X ‘ + IntToStr(Serie);
ChartFX.ThisSerie = Pred(Serie);
for Point = 0 to Pred(NumReturn) do
ChartFX.Value[Point] = mSolutionValues[Serie, Succ(Point)];
end;
ChartFX.CloseData(Cod_Values);
{
ChartFX.OpenDataEx(Cod_XValues, Order, Pred(NumReturn));
for Serie = 1 to Order do
begin
ChartFX.ThisSerie = Pred(Serie);
for Point = 0 to Pred(NumReturn) do
ChartFX.XValue[Point] = mSolutionValues[1, Succ(Point)];
end;
ChartFX.CloseData(Cod_XValues);
}
end;
2 begin
LoadMatrixSolveFromStrGrd(mSolutionValues, StringGrid_Solve2);
SetAdm;
ChartFX.OpenDataEx(Cod_Values, Order, Pred(NumReturn));
for Serie = 1 to Order do
begin
ChartFX.SerLeg[Pred(Serie)] = ‘X ‘ + IntToStr(Serie);
ChartFX.ThisSerie = Pred(Serie);
for Point = 0 to Pred(NumReturn) do
ChartFX.Value[Point] = mSolutionValues[Serie, Succ(Point)];
end;
ChartFX.CloseData(Cod_Values);
end;
end;
mSolutionValues.Free;
end;
end;

procedure ComputeFromPage2;
begin
ComputeFromPage1;
end;

procedure ComputeFromPage3;
begin
case Form_Main.RadioGroupChart.ItemIndex of
0 ShowChart(1);
1 ShowChart(2);
end;
end;

end.
unit Matrix;

interface

uses SysUtils;

type
Float = Extended;
EMatrixOperatingError = class (Exception);

const
NearlyZero = 1E-15;

type
TMatrix = class (TObject)
private
DataPtr Pointer;
FCols, FRows Word;
function GetCell (ACol, ARow Word) Float;
procedure SetCell (ACol, ARow Word; AValue Float);
function GetItem (NumItem LongInt) Float;
procedure SetItem (NumItem LongInt; AValue Float);
procedure SwitchRows (FirstRow, SecondRow Word);
public
constructor Create (NCols, NRows Word);
destructor Destroy; override;
procedure Assign (AMatrix TMatrix);
procedure ReSize (NewCols, NewRows Word);
procedure SetNull;
procedure SetSingle;
procedure SetNegative;
procedure AddConst (AConst Float);
procedure AddMatrix (AMatrix TMatrix);
procedure MultConst (MConst Float);
procedure MultFromRight (MMatrix TMatrix);
procedure MultFromLeft (MMatrix TMatrix);
procedure NthPower (Power Word);
procedure Transpose;
function Inverse Boolean;
function Determinant Float;
function Rang Float;
property ColCount Word read FCols;
property RowCount Word read FRows;
property Cells [ACol, ARow Word] Float read GetCell write SetCell; default;
property Items [NumItem LongInt] Float read GetItem write SetItem;
end;

implementation

uses Windows;

function IncPtr (p Pointer; i LongInt) Pointer;
asm
push EBX
mov EBX,EAX
add EBX,EDX
mov EAX,EBX
pop EBX
end;

function TMatrix.GetCell (ACol, ARow Word) Float;
var
CellPtr ^Float;
begin
CellPtr = IncPtr(DataPtr, (FRows * Pred(ACol) + Pred(ARow)) * SizeOf(Float));
Result = CellPtr^;
end;

procedure TMatrix.SetCell (ACol, ARow Word; AValue Float);
var
CellPtr ^Float;
begin
CellPtr = IncPtr(DataPtr, (FRows * Pred(ACol) + Pred(ARow)) * SizeOf(Float));
CellPtr^ = AValue;
end;

function TMatrix.GetItem (NumItem LongInt) Float;
var
CellPtr ^Float;
begin
CellPtr = IncPtr(DataPtr, Pred(NumItem) * SizeOf(Float));
Result = CellPtr^;
end;

procedure TMatrix.SetItem (NumItem LongInt; AValue Float);
var
CellPtr ^Float;
begin
CellPtr = IncPtr(DataPtr, Pred(NumItem) * SizeOf(Float));
CellPtr^ = AValue;
end;

procedure TMatrix.SwitchRows (FirstRow, SecondRow Word);
var
i Word;
Buffer Float;
begin
for i = 1 to FCols do
begin
Buffer = GetCell(i, FirstRow);
SetCell(i, FirstRow, GetCell(i, SecondRow));
SetCell(i, SecondRow, Buffer);
end;
end;

constructor TMatrix.Create (NCols, NRows Word);
begin
inherited Create;
FCols = NCols;
FRows = NRows;
DataPtr = AllocMem(FCols * FRows * SizeOf(Float));
end;

destructor TMatrix.Destroy;
begin
FreeMem(DataPtr);
inherited Destroy;
end;

procedure TMatrix.Assign (AMatrix TMatrix);
var
NewMatrixSize LongInt;
begin
NewMatrixSize = AMatrix.ColCount * AMatrix.RowCount * SizeOf(Float);
ReAllocMem(DataPtr, NewMatrixSize);
CopyMemory(DataPtr, AMatrix.DataPtr, NewMatrixSize);
FCols = AMatrix.ColCount;
FRows = AMatrix.RowCount
end;

procedure TMatrix.ReSize (NewCols, NewRows Word);
var
NewMatrixSize LongInt;
begin
NewMatrixSize = NewCols * NewRows * SizeOf(Float);
ReAllocMem(DataPtr, NewMatrixSize);
FCols = NewCols;
FRows = NewRows;
end;

procedure TMatrix.SetNull;
begin
ZeroMemory (DataPtr, FCols * FRows * SizeOf(Float));
end;

procedure TMatrix.SetSingle;
var
i Word;
begin
if FCols <> FRows then
Raise EMatrixOperatingError.Create (‘Единичная матрица должна быть ‘+
‘квадратной’)
else
begin
SetNull;
for i = 1 to FCols do SetCell (i, i, 1);
end;
end;

procedure TMatrix.SetNegative;
var
i LongInt;
begin
for i = 1 to FCols * FRows do SetItem(i, — GetItem(i));
end;

procedure TMatrix.AddConst (AConst Float);
var
i LongInt;
begin
for i = 1 to FCols * FRows do SetItem (i, GetItem(i) + AConst);
end;

procedure TMatrix.AddMatrix (AMatrix TMatrix);
var
i LongInt;
begin
for i = 1 to FCols * FRows do SetItem (i, GetItem(i) + AMatrix.Items [i]);
end;

procedure TMatrix.MultConst (MConst Float);
var
i LongInt;
begin
for i = 1 to FCols * FRows do SetItem (i, GetItem(i) * MConst);
end;

procedure TMatrix.MultFromRight (MMatrix TMatrix);
var
j, i, k Word;
DummyRes Float;
DummyMatrix TMatrix;
begin
DummyMatrix = TMatrix.Create (MMatrix.ColCount, FRows);
if FCols <> MMatrix.RowCount then
Raise EMatrixOperatingError.Create (‘Перемножаемые матрицы должны быть ‘+
‘соответствующей размерности’)
else
for i = 1 to FRows do
for j = 1 to MMatrix.ColCount do
begin
DummyRes = 0;
for k = 1 to FCols do
DummyRes = DummyRes + Cells[k, i] * MMatrix[j, k];
DummyMatrix[j, i] = DummyRes;
end;
Assign(DummyMatrix);
DummyMatrix.Free;
end;

procedure TMatrix.MultFromLeft (MMatrix TMatrix);
var
j, i, k Word;
DummyRes Float;
DummyMatrix TMatrix;
begin
DummyMatrix = TMatrix.Create (FCols, MMatrix.RowCount);
if MMatrix.ColCount <> FRows then
Raise EMatrixOperatingError.Create (‘Перемножаемые матрицы должны быть ‘+
‘соответствующей размерности’)
else
for i = 1 to MMatrix.ColCount do
for j = 1 to FCols do
begin
DummyRes = 0;
for k = 1 to MMatrix.ColCount do
DummyRes = DummyRes + MMatrix[k, i] * Cells[j, k];
DummyMatrix[j, i] = DummyRes;
end;
Assign(DummyMatrix);
DummyMatrix.Free;
end;

procedure TMatrix.NthPower (Power Word);
var
i Word;
DummyMatrix TMatrix;
begin
DummyMatrix = TMatrix.Create (FCols, FRows);
DummyMatrix.Assign (Self);
if FCols <> FRows then
Raise EMatrixOperatingError.Create (‘Возводимая в степень матрица должна ‘+
‘быть квадратной’)
else
case Power of
0 SetSingle;
1 begin end;
else
for i = 2 to Power do MultFromRight (DummyMatrix);
end;
DummyMatrix.Free;
end;

procedure TMatrix.Transpose;
var
i, j Word;
Dummy Float;
begin
if FCols <> FRows then
Raise EMatrixOperatingError.Create (‘Транспонируемая матрица должна быть ‘+
‘квадратной’)
else
for i = 1 to FCols do
for j = 1 to FRows do
if j > i then
begin
Dummy = GetCell(j, i);
SetCell(j, i, GetCell(i, j));
SetCell(i, j, Dummy);
end
end;

function TMatrix.Inverse Boolean;
var
DummyMatrix TMatrix;
Divisor, Multiplier Float;
Row, RefRow, NewRow, Term Word;
Singular Boolean;
begin
Singular = False;
DummyMatrix = TMatrix.Create (FCols, FRows);
if (FCols <> FRows) or (FCols = 0) then
Raise EMatrixOperatingError.Create (‘Инвертируемая матрица должна быть ‘+
‘квадратной и ненулевого размера’);
if FCols = 1 then
if ABS(GetItem(1)) < NearlyZero then Singular = True
else DummyMatrix.Items[1] = 1 / GetItem(1);
if FCols > 1 then
begin
DummyMatrix.SetSingle;
RefRow = 0;
repeat
Inc(RefRow);
if ABS(Cells[RefRow, RefRow]) < NearlyZero then
begin
Singular = TRUE;
NewRow = RefRow;
repeat
Inc(NewRow);
if ABS(Cells[RefRow, NewRow]) > NearlyZero then
begin
SwitchRows(NewRow, RefRow);
DummyMatrix.SwitchRows(NewRow, RefRow);
Singular = False;
end;
until (not Singular) or (NewRow >= FCols);
end;
if not Singular then
begin
Divisor = Cells[RefRow, RefRow];
for Term = 1 to FCols do
begin
SetCell(Term, RefRow, GetCell(Term, RefRow)/Divisor);
DummyMatrix[Term, RefRow] = DummyMatrix[Term, RefRow]/Divisor;
end;
for Row = 1 to FCols do
if (Row <> RefRow) and (ABS(Cells[RefRow, Row]) > NearlyZero) then
begin
Multiplier = — Cells[RefRow, Row] / Cells[RefRow, RefRow];
for Term = 1 to FCols do
begin
SetCell(Term, Row, GetCell(Term, Row) +
Multiplier * GetCell(Term, RefRow));
DummyMatrix[Term, Row] = DummyMatrix[Term, Row] +
Multiplier * DummyMatrix[Term, RefRow];
end
end;
end;
until Singular or (RefRow >= FCols);
end;
Assign(DummyMatrix);
DummyMatrix.Free;
if not Singular then Result = True
else Result = False;
end;

function TMatrix.Determinant Float;
begin
Result = 0;
end;

function TMatrix.Rang Float;
begin
Result = 0;
end;

end.
unit Operates;

interface

uses Matrix, Grids, SysUtils;

const
MaxArraySize = 30;

type
Float = Extended;
TOrder = 1..MaxArraySize;
ESingularMatrix = class (Exception);

type
TComplex = record
Re, Im Float;
end;

TComplexVector = record
Data array [1..MaxArraySize] of TComplex;
Dim TOrder;
end;

function SymmetricalFunction (Roots TComplexVector; K byte) Float;
procedure DiffSystemSolve (matrixA,
matrixB TMatrix;
LowerLimit,
UpperLimit Float;
InitialValues TMatrix;
NumReturn,
NumIntervals Word;
SolutionValues TMatrix);

implementation

function SymmetricalFunction (Roots TComplexVector; K byte) Float;
var
Z TComplex;

function SummComplex (FirstNC, SecondNC TComplex) TComplex;
begin
Result.Re = FirstNC.Re + SecondNC.Re;
Result.Im = FirstNC.Im + SecondNC.Im;
end;

function MultComplex (FirstNC, SecondNC TComplex) TComplex;
begin
Result.Re = FirstNC.Re * SecondNC.Re — FirstNC.Im * SecondNC.Im;
Result.Im = FirstNC.Re * SecondNC.Im + FirstNC.Im * SecondNC.Re;
end;

function DivComplex (FirstNC, SecondNC TComplex) TComplex;
var
Z Float;
begin
Z = Sqr(SecondNC.Re) + Sqr(SecondNC.Im);
Result.Re = (FirstNC.Re * SecondNC.Re + FirstNC.Im * SecondNC.Im) / Z;
Result.Im = (FirstNC.Im * SecondNC.Re — FirstNC.Re * SecondNC.Im) / Z;
end;

function CombinationSumm (LowLimit, HighLimit, K byte) TComplex;
var i byte;
begin
Result.Re = 0;
Result.Im = 0;
if LowLimit = HighLimit then Result = Roots.Data[LowLimit]
else
for i = LowLimit to HighLimit — K + 1 do
if K = 1 then Result = SummComplex(Result, Roots.Data [i])
else Result = SummComplex(Result,
MultComplex(Roots.Data [i],
CombinationSumm(i + 1,
HighLimit, K-1)));

end;
begin
Z = CombinationSumm(1, Roots.Dim, K);
Result = Z.Re;
end;

procedure DiffSystemSolve (matrixA, matrixB TMatrix;
LowerLimit, UpperLimit Float;
InitialValues TMatrix;
NumReturn, NumIntervals Word;
SolutionValues TMatrix);
type
Ptr = ^Data;
Data = record
Values TMatrix;
Next Ptr;
end;
var
ValuesStack Ptr;
Spacing, HalfSpacing Float;
Index, Term Word;
F1, F2, F3, F4,
CurrentValues,
TempValues TMatrix;
NumEquations, NumTimeCol Word;

function TargetALL (matrixA, mayrixB TMatrix; Values TMatrix; KRow Word) Float;
var
j Word;
begin
try
Result = matrixB.Items[KRow];
for j = 1 to NumEquations do
Result = Result + matrixA[j, KRow] * Values.Items[j];
except
on EO EOverflow do EO.Message = ‘Не буду считать !!!’#10 +
‘С уменьшите разброс коэффициентов в матрице А’#10 +
‘либо измените опции (уменьшите их pls.)’;
end;
end;

procedure Push (var ValuesStack Ptr;
CurrentValues TMatrix);
var
NewNode Ptr;
begin
New(NewNode);
NewNode^.Values = TMatrix.Create(NumTimeCol, 1);
NewNode^.Values.Assign(CurrentValues);
NewNode^.Next = ValuesStack;
ValuesStack = NewNode;
end; { procedure Push }

procedure Pop (var ValuesStack Ptr;
CurrentValues TMatrix);
var
OldNode Ptr;
begin
OldNode = ValuesStack;
ValuesStack = OldNode^.Next;
CurrentValues.Assign(OldNode^.Values);
OldNode^.Values.Free;
Dispose(OldNode);
end; { procedure Pop }

procedure GetValues(NumReturn, NumIntervals Word; var ValuesStack Ptr;
SolutionValues TMatrix);
var
Index, Term Integer;
j Word;
CurrValues TMatrix;
begin
SolutionValues.ReSize(NumTimeCol, Succ(NumReturn));
CurrValues = TMatrix.Create(NumTimeCol, 1);
Term = NumIntervals;
for Index = NumReturn downto 0 do
begin
Pop(ValuesStack, CurrValues);
Dec(Term);
while (Term / NumIntervals >= Index / NumReturn) and (Term >= 0) do
begin
Pop(ValuesStack, CurrValues);
Dec(Term);
end;
for j = 1 to NumTimeCol do
SolutionValues[j, Succ(Index)] = CurrValues.Items[j];
end;
CurrValues.Free;
end; { procedure GetValues }

procedure Step(Spacing Float; CurrentValues TMatrix; F TMatrix);
var
i byte;
begin
for i = 1 to NumEquations do
F.Items[i] = Spacing * TargetALL (matrixA, matrixB, CurrentValues, i);
end; { procedure Step }

begin
NumEquations = matrixA.RowCount;
NumTimeCol = Succ(NumEquations);
ValuesStack = nil;
Spacing = (UpperLimit — LowerLimit) / NumIntervals;
CurrentValues = TMatrix.Create(1, 1);
CurrentValues.Assign(InitialValues);
CurrentValues.ReSize(NumTimeCol, 1);
CurrentValues.Items[NumTimeCol] = LowerLimit;
TempValues = TMatrix.Create(NumTimeCol, 1);
F1 = TMatrix.Create(NumTimeCol, 1);
F2 = TMatrix.Create(NumTimeCol, 1);
F3 = TMatrix.Create(NumTimeCol, 1);
F4 = TMatrix.Create(NumTimeCol, 1);
Push(ValuesStack, CurrentValues);
HalfSpacing = Spacing / 2;
for Index = 1 to NumIntervals do
begin
{ First step — calculate F1 }
Step(Spacing, CurrentValues, F1);
TempValues.Items[NumTimeCol] = CurrentValues.Items[NumTimeCol] + HalfSpacing;
for Term = 1 to NumEquations do
TempValues.Items[Term] = CurrentValues.Items[Term] + 0.5 * F1.Items[Term];

{ 2nd step — calculate F2 }
Step(Spacing, TempValues, F2);
for Term = 1 to NumEquations do
TempValues.Items[Term] = CurrentValues.Items[Term] + 0.5 * F2.Items[Term];

{ Third step — calculate F3 }
Step(Spacing, TempValues, F3);
TempValues.Items[NumTimeCol] = CurrentValues.Items[NumTimeCol] + Spacing;
for Term = 1 to NumEquations do
TempValues.Items[Term] = CurrentValues.Items[Term] + F3.Items[Term];

{ Fourth step — calculate F4[1]; first equation }
Step(Spacing, TempValues, F4);

{ Combine F1, F2, F3, and F4 to get }
{ the solution at this mesh point }
CurrentValues.Items[NumTimeCol] = CurrentValues.Items[NumTimeCol] + Spacing;
for Term = 1 to NumEquations do
CurrentValues.Items[Term] = CurrentValues.Items[Term] +
(F1.Items[Term] + 2 * F2.Items[Term] +
2 * F3.Items[Term] + F4.Items[Term]) /6;
Push(ValuesStack, CurrentValues);
end;
GetValues(NumReturn, NumIntervals, ValuesStack, SolutionValues);
F1.Free;
F2.Free;
F3.Free;
F4.Free;
CurrentValues.Free;
TempValues.Free;
end;

end.
unit HelpUnit;

interface

uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls, ExtCtrls, Buttons;

type
TForm_Help = class(TForm)
BitBtn1 TBitBtn;
Bevel1 TBevel;
Label1 TLabel;
Label2 TLabel;
Label3 TLabel;
Label4 TLabel;
Label5 TLabel;
Label6 TLabel;
Label7 TLabel;
private
{ Private declarations }
public
{ Public declarations }
end;

var
Form_Help TForm_Help;

implementation

{$R *.DFM}

end.
unit OptsUnit;

interface

uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls, Spin, Buttons, ExtCtrls;

type
TForm_Options = class(TForm)
CheckBox_Link TCheckBox;
SpinEdit1 TSpinEdit;
SpinEdit2 TSpinEdit;
SpinEdit3 TSpinEdit;
Label_UpLimit TLabel;
Label_PointsNumber TLabel;
Label_Intervals TLabel;
Label_1 TLabel;
BitBtn_Ok TBitBtn;
BitBtn_Cancel TBitBtn;
SpinEdit0 TSpinEdit;
Label1 TLabel;
Bevel1 TBevel;
procedure SpinEdit0Change(Sender TObject);
procedure SpinEdit2Change(Sender TObject);
procedure CheckBox_LinkClick(Sender TObject);
private
{ Private declarations }
public
{ Public declarations }
end;

var
Form_Options TForm_Options;

implementation

uses MainUnit, SubUnit;

{$R *.DFM}

procedure TForm_Options.SpinEdit0Change(Sender TObject);
begin
SpinEdit1.MinValue = Succ(SpinEdit0.Value);
if SpinEdit1.Value < SpinEdit1.MinValue then SpinEdit1.Value = SpinEdit1.MinValue;
end;

procedure TForm_Options.SpinEdit2Change(Sender TObject);
begin
SpinEdit3.MinValue = SpinEdit2.Value;
if SpinEdit3.Value < SpinEdit3.MinValue then SpinEdit3.Value = SpinEdit3.MinValue;
end;

procedure TForm_Options.CheckBox_LinkClick(Sender TObject);
begin
if CheckBox_Link.State = cbChecked then Form_Main.BindGrids
else Form_Main.UnBindGrids
end;

end.