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