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.