Программная реализация модального управления для линейных стационарных систем
Курсовая работа:
«Программная реализация модального управления
для линейных стационарных систем»Постановка задачи:
1. Для объекта управления с математическим описанием
[pic], (1) [pic]- задано,
где [pic] — n-мерный вектор состояния, [pic],
[pic]- начальный вектор состояния,
[pic]- скалярное управление,
[pic]- матрица действительных коэффициентов,
[pic]- матрица действительных коэффициентов,
найти управление в функции переменных состояния объекта, т.е.
[pic], (2)
где[pic]- матрица обратной связи, такое, чтобы замкнутая система была
устойчивой.
2. Корни характеристического уравнения замкнутой системы
[pic] (3)
должны выбираться по усмотрению (произвольно) с условием устойчивости
системы (3).Задание:
1. Разработать алгоритм решения поставленной задачи.
2. Разработать программу решения поставленной задачи с интерактивным
экранным интерфейсом в системах Borland Pascal, Turbo Vision, Delphi — по
выбору.
3. Разработать программу решения систем дифференциальных уравнений (1) и
(3) с интерактивным экранным интерфейсом.
4. Разработать программу графического построения решений систем (1) и (3) с
интерактивным экранным интерфейсом.Введение
Наряду с общими методами синтеза оптимальных законов управления для
стационарных объектов всё большее применение находят методы, основанные на
решении задачи о размещении корней характеристического уравнения замкнутой
системы в желаемое положение. Этого можно добиться надлежащим выбором
матрицы обратной связи по состоянию. Решение указанной задачи является
предметом теории модального управления (термин связан с тем, что корням
характеристического уравнения соответствуют составляющие свободного
движения, называемые модами).Алгоритм модального управления.
Соглашения:
. Задаваемый объект управления математически описывается уравнением
[pic], (1)
где [pic] и [pic] — матрицы действительных коэффициентов,
[pic] — n-мерный вектор состояния
[pic]- скалярное управление,
[pic] — порядок системы (1).
. Обратная связь по состоянию имеет вид
[pic], (2)
где[pic]- матрица обратной связи.
. Система с введенной обратной связью описывается уравнением
[pic] (3)
. Характеристическое уравнение системы (1) имеет вид
[pic] (4)
. Характеристическое уравнение системы (3) с задаваемыми (желаемыми)
корнями [pic]имеет вид
[pic] (5)Алгоритм:
1. Для исходной системы (1) составляем матрицу управляемости
[pic]
2. Обращаем матрицу [pic], т.е. вычисляем [pic].
Если [pic] не существует (т.е. матрица [pic] — вырожденная), то
прекращаем вычисления: полное управление корнями характеристического
уравнения (5) не возможно.
3. Вычисляем матрицу [pic]
4. Составляем матрицу
[pic]
5. Вычисляем матрицу, обратную матрице [pic], т.е. [pic]
6. Вычисляем матрицу [pic] — матрицу [pic] в канонической форме фазовой
переменной:
[pic]
где [pic]- коэффициенты характеристического уравнения (4).
Матрица [pic] в канонической форме имеет вид
[pic]
7. Составляем вектор [pic] , элементам которого являются коэффициенты
характеристического уравнения (4), т.е. [pic], [pic],
где [pic] — элементы матрицы [pic].
8. Находим коэффициенты характеристического уравнения (5) (см. пояснения) и
составляем из них вектор [pic].
9. Вычисляем вектор [pic].
[pic] — искомая матрица обратной связи системы (3), но она вычислена для
системы, матрицы которой заданы в канонической форме фазовой переменной
([pic] и [pic]).
10. Для исходной системы (3) матрица обратной связи получается по формуле
[pic]
Матрица [pic] — искомая матрица обратной связи.Пояснения к алгоритму:
В данной работе рассматривается случай, когда управление единственно и
информация о переменных состояния полная. Задача модального управления
тогда наиболее просто решается, если уравнения объекта заданы в
канонической форме фазовой переменной.
Так как управление выбрано в виде линейной функции переменных состояния
[pic], где [pic] является матрицей строкой [pic]. В таком случае уравнение
замкнутой системы приобретает вид [pic]. Здесь
[pic]
[pic]
Характеристическое уравнение такой замкнутой системы будет следующим
[pic]
Поскольку каждый коэффициент матрицы обратной связи [pic] входит только в
один коэффициент характеристического уравнения, то очевидно, что выбором
коэффициентов [pic] можно получить любые коэффициенты характеристического
уравнения, а значит и любое расположение корней.
Если же желаемое характеристическое уравнение имеет вид
[pic],
то коэффициенты матрицы обратной связи вычисляются с помощью соотношений:
[pic]
Если при наличии одного управления нормальные уравнения объекта заданы не
в канонической форме (что наиболее вероятно), то, в соответствии с пунктами
№1-6 алгоритма, от исходной формы с помощью преобразования [pic] или [pic]
нужно перейти к уравнению [pic] в указанной канонической форме.
Управление возможно, если выполняется условие полной управляемости (ранг
матрицы управляемости M должен быть равен n). В алгоритме об управляемости
системы судится по существованию матрицы [pic]: если она существует, то
ранг матрицы равен ее порядку (n). Для объекта управления с единственным
управлением матрица [pic] оказывается также единственной.
Для нахождения коэффициентов [pic] характеристического уравнения (5), в
работе используется соотношения между корнями [pic] и коэффициентами
[pic] линейного алгебраического уравнения степени n:
[pic], (k = 1, 2, … , n)
где многочлены [pic]- элементарные симметрические функции, определяемые
следующим образом:
[pic],
[pic],
[pic],
…
[pic]
где Sk — сумма всех [pic] произведений, каждое из которых содержит 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.