Метод Эйлера и метод Рунге-Кутта: проверить код

@xSimplex · Регистрация: 24.07.2013

Студворк — интернет-сервис помощи студентам

Доброго времени суток. Хотел бы обратится к вам за помощью. Я написал программку, которая решает сит. диф. ур-й, двумя методами. Но метод Рунге-Кутты вроде не правильный, не могли бы вы проверить его, а так же математику, может быть где-то допустил ошибку. Так же до написание этих метод, я разобрал их на элементарных примерах и там вроде все работает, прикреплю.

Спасибо всем

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
unit uIntegrator;
 
interface
 
uses
  Math;
 
const
  MU = 3.986e5;
 
  type
    TData = array [1..4]of Extended;
    TOrbitMean = array[1..5] of Extended;
 
//var
//TStat_Vec : TData;
 
// Äåéñâèÿ ñ ìàòðèöàìè
function N_M(N: Extended; M: TData): TData;
function M_M(M1,M2: TData): TData;
 
// Ñèñòåìà äèô-ûõ óðàâíåíèé
function Eq(T: Extended; StatVec: TData): TData;
function Eq_1(T: Extended; StatVec: TData; Lam,F: Extended): TData;
 
// Ìåòîäû èíòåãðèðîâàíèÿ
function EL(T,Step: Extended; StatVec: TData): TData;
function EL_1(T,Step: Extended; StatVec: TData; Lam,F: Extended): TData;
 
function RK(T,Step: Extended; StatVec: TData): TData;
function RK_1(T,Step: Extended; StatVec: TData; Lam,F: Extended): TData;
 
// Íàõîæäåíèå ýëåìåíòîâ îðáèòû
function OrbitMean(T: Integer; StatVec: TData): TOrbitMean;
 
implementation
 
 
// Äåéñâèÿ ñ ìàòðèöàìè
 
// ×èñëî íà âåêòîð-ñòðîêó
function N_M(N: Extended; M: TData): TData;
var
  i: Integer;
begin
    for i:= 1 to 4 do
    begin
      Result[i]:= N*M[i];
    end;
end;
 
// Ìàññèâ ïëþñ Ìàññèâ
function M_M(M1, M2: TData): TData;
var
  i: Integer;
begin
    for i:= 1 to 4 do
    begin
      Result[i]:= M1[i]+M2[i];
    end;
end;
 
 
// Ñèñòåìà äèô-ûõ óðàâíåíèé
 
// Áåç òÿãè äâèãàòåëÿ
function Eq(T: Extended; StatVec: TData): TData;
begin
  Result[1]:= StatVec[3];
  Result[2]:= StatVec[4]/StatVec[1];
  Result[3]:= StatVec[4]*StatVec[4]/StatVec[1] - MU/(StatVec[1]*StatVec[1]);
  Result[4]:=  -StatVec[3]*StatVec[4]/StatVec[1];
end;
 
// Ñ òÿãîé äâèãàòåëÿ
function Eq_1(T: Extended; StatVec: TData; Lam,F: Extended): TData;
begin
  Result[1]:= StatVec[3];
  Result[2]:= StatVec[4]/StatVec[1];
  Result[3]:= StatVec[4]*StatVec[4]/StatVec[1] - MU/(StatVec[1]*StatVec[1]) + F*cos(Lam);
  Result[4]:= F*sin(Lam) - StatVec[3]*StatVec[4]/StatVec[1];
end;
 
// Ìåòîäû èíòåãðèðîâàíèÿ
 
// Ìåòîä Ýéëåðà
function EL(T,Step: Extended; StatVec: TData): TData;
begin
  Result:= M_M(StatVec,N_M(Step,Eq(T,StatVec)));
end;
 
function EL_1(T,Step: Extended; StatVec: TData; Lam,F: Extended): TData;
begin
  Result:= M_M(StatVec,N_M(Step,Eq_1(T,StatVec,Lam,F)));;
end;
 
 
// Ìåòîä Ðóíãå-Êóòòû
function RK(T,Step: Extended; StatVec: TData): TData;
var
  k1,k2,k3,k4: TData;
 
  eq1,eq2,eq3,eq4: TData;
begin
  k1:= N_M(Step,Eq(T,StatVec));
  k2:= N_M(Step,Eq(T+Step/2,N_M(0.5,k1)));
  k3:= N_M(Step,Eq(T+Step/2,N_M(0.5,k2)));
  k4:= N_M(Step,Eq(T+Step,k3));
 
  //Result:= M_M(  StatVec, N_M(Step/6, M_M( M_M(k1, N_M(2,k2) ), M_M(N_M(2,k3),k4) ))  ) ;
 
  eq1:= M_M(k1,N_M(2,k2));
  eq2:= M_M(k4,N_M(2,k3));
  eq3:= M_M(eq1,eq2);
  eq4:= N_M(step/6,eq3);
 
  Result:= M_M(StatVec,eq4);
 
  
end;
 
function RK_1(T,Step: Extended; StatVec: TData; Lam,F: Extended): TData;
var
  k1,k2,k3,k4: TData;
begin
  k1:= N_M(Step,Eq_1(T,StatVec,Lam,F));
  k2:= N_M(Step,Eq_1(T+Step/2,N_M(0.5,k1),Lam,F));
  k3:= N_M(Step,Eq_1(T+Step/2,N_M(0.5,k2),Lam,F));
  k4:= N_M(Step,Eq_1(T+Step,k3,Lam,F));
 
  Result:= M_M(  StatVec, N_M(Step/6, M_M( M_M(k1, N_M(2,k2) ), M_M(N_M(2,k3),k4) ))  ) ;
end;
 
// Íàõîæäåíèå ýëåìåíòîâ îðáèòû
 
function OrbitMean(T: Integer; StatVec: TData): TOrbitMean;
var
  CosAnom: Extended;
begin
    Result[1]:= 1/MU * 
    sqrt(
    (StatVec[1] * StatVec[4] * StatVec[4] - MU) * (StatVec[1] * StatVec[4] * StatVec[4] - MU)
    + 
    (StatVec[1]*StatVec[3] * StatVec[4]) * (StatVec[1] * StatVec[3] * StatVec[4])
    );
    
    CosAnom:= (StatVec[1] * StatVec[4] * StatVec[4] - 1)
    /
    (
    (StatVec[1] * StatVec[4] * StatVec[4] - MU) * (StatVec[1] * StatVec[4] * StatVec[4] - MU)
    +
    (StatVec[1] * StatVec[3] * StatVec[4]) * (StatVec[1] * StatVec[3] * StatVec[4])
    );
    
    Result[2]:= StatVec[1] * (1 - Result[1] * CosAnom);
    
    Result[3]:= Result[2] / (1 - Result[1] * Result[1]);
    
    Result[4]:= Result[2] / (1 + Result[1]);
    
    Result[5]:= Result[2] / (1 - Result[1]);
end;
 
 
end.

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
unit uMain;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ExtCtrls, TeeProcs, TeEngine, Chart, ComCtrls, Series,Math,
   uIntegrator, Grids;
 
const
  g = 9.81;
 
type
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    GroupBox2: TGroupBox;
    PageControl1: TPageControl;
    TabSheet1: TTabSheet;
    TabSheet2: TTabSheet;
    TabSheet3: TTabSheet;
    TabSheet4: TTabSheet;
    TabSheet5: TTabSheet;
    TabSheet6: TTabSheet;
    Chart1: TChart;
    Chart2: TChart;
    Chart3: TChart;
    Chart4: TChart;
    Chart5: TChart;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Label5: TLabel;
    Label6: TLabel;
    Label7: TLabel;
    Label8: TLabel;
    Label9: TLabel;
    Label10: TLabel;
    ComboBox1: TComboBox;
    Button1: TButton;
    Button2: TButton;
    Edit1: TEdit;
    Edit2: TEdit;
    Edit3: TEdit;
    Edit4: TEdit;
    Label11: TLabel;
    Edit5: TEdit;
    Edit6: TEdit;
    Edit7: TEdit;
    Edit8: TEdit;
    Series1: TLineSeries;
    Series2: TLineSeries;
    Series3: TLineSeries;
    Series4: TLineSeries;
    Series5: TLineSeries;
    Memo1: TMemo;
    
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
 
  end;
 
var
  Form1: TForm1;
 
implementation
 
{$R *.dfm}
 
procedure TForm1.Button1Click(Sender: TObject);
var
  Metod_Name: Integer;
  a,m,s,h,i: Integer;
  u,w,e,p,anom,Vr,Vt,r,F,lam: Extended;
 
  T,T_0: Integer;
  T1: Extended;
 
  Iter: Extended;
  Iter1: Integer;
 
  EqM: TData;
  OrbM: TOrbitMean;
begin
  T_0:= 0;
 
  a:= StrToInt(Edit1.Text);
  m:= StrToInt(Edit5.Text);
  s:= StrToInt(Edit6.Text);
  h:= StrToInt(Edit8.Text);
 
  u:= StrToFloat(Edit4.Text);
  w:= StrToFloat(Edit3.Text);
  lam:= StrToFloat(Edit7.Text);
  e:= StrToFloat(Edit2.Text);
 
  u:= DegToRad(u);
  w:= DegToRad(w);
  lam:= DegToRad(lam);
 
 
  T1:= 2*pi*a*sqrt(a/mu);
    p:= a*(1-sqr(e));
    anom:= u-w;
    r:= p/(1-e*cos(anom));
    Vr:= sqrt(mu/p)*sin(anom);
    Vt:= sqrt(mu/p)*(1-e*cos(anom));
    F:= (m*g-s)/m;
 
  EqM[1]:= r;
  EqM[2]:= u;
  EqM[3]:= Vr;
  EqM[4]:= Vt;
 
 
  // Óçíàåì êîë-âî èòåðàöèé
  Iter:= (Round(T1)-T_0)/h;
  Iter1:= Round(Iter);
 
  if ComboBox1.ItemIndex = 0 then
  begin
    Memo1.Clear;
    Series1.Clear;
    Series2.Clear;
    Series3.Clear;
    Series4.Clear;
    Series5.Clear;
 
    Memo1.Lines.Add('Âûáðàí ìåòîä Ýéëåðà');
 
    for i:= 0 to (Iter1-1) do
    begin
      T_0:= T_0+h;
 
     // EqM:= EL_1(T_0, h, EqM, lam, F);  // Ñ òÿãîé äâèãàòåëÿ
      EqM:= EL(T_0, h, EqM); // Áåç òÿãè äâèãàòåëÿ
 
      OrbM:= OrbitMean(T_0, EqM);
 
 
      Memo1.Lines.Add( 'e(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[1]) + ';  ' +
                       'p(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[2]) + ';  ' +
                       'a(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[3]) + ';  ' +
                       'Rï(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[4]) + ';  ' +
                       'Rà(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[5]) + ';  '
                      );
 
      Series1.AddXY(T_0,OrbM[1]);
      Series2.AddXY(T_0,OrbM[2]);
      Series3.AddXY(T_0,OrbM[3]);
      Series4.AddXY(T_0,OrbM[4]);
      Series5.AddXY(T_0,OrbM[5]);
    end;
  end;
 
  if ComboBox1.ItemIndex = 1 then
  begin
    Memo1.Clear;
    Series1.Clear;
    Series2.Clear;
    Series3.Clear;
    Series4.Clear;
    Series5.Clear;
 
    Memo1.Lines.Add('Âûáðàí ìåòîä Ðóíãå-Êóòòû');
 
    for i:= 0 to (Iter1-1) do
    begin
      T_0:= T_0+h;
 
     // EqM:= RK_1(T_0, h, EqM, lam, F);  // Ñ òÿãîé äâèãàòåëÿ
      EqM:= RK(T_0, h, EqM); // Áåç òÿãè äâèãàòåëÿ
 
      OrbM:= OrbitMean(T_0, EqM);
 
 
      Memo1.Lines.Add( 'e(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[1]) + ';  ' +
                       'p(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[2]) + ';  ' +
                       'a(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[3]) + ';  ' +
                       'Rï(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[4]) + ';  ' +
                       'Rà(' + IntToStr(T_0)+ ') = ' + FloatToStr(OrbM[5]) + ';  '
                      );
 
      Series1.AddXY(T_0,OrbM[1]);
      Series2.AddXY(T_0,OrbM[2]);
      Series3.AddXY(T_0,OrbM[3]);
      Series4.AddXY(T_0,OrbM[4]);
      Series5.AddXY(T_0,OrbM[5]);
    end;
  end;
 
end;
 
procedure TForm1.Button2Click(Sender: TObject);
begin
 Form1.Close;
end;
 
 
end.

@xSimplex · 19.08.2014, 23:08 **[ТС]**

Если кто-то читал, то могу сказать, что я еще раз проверил и у меня получаются слишком большие значения при расчете коэф. в рунге-кутте. Из-за чего это может быть? Ошибка в исходных ур-иях или я что-то криво подал??

@xSimplex · 24.08.2014, 04:41 **[ТС]**

Ошибку нашел

@NotBeginner · 24.08.2014, 13:28

Сообщение от xSimplex

Ошибку нашел

выкладывай рабочий код

@xSimplex · 24.08.2014, 19:27 **[ТС]**

Вот тут все работает. Можете кстати проверить, вдруг что-то упустил, прошу проверить сами методы. Там мне еще надо провести пару манипуляций для сохранения значений в файл и т.п. , ну главное что скелет есть, а мясо нарастет

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
unit uIntegrator;
 
interface
 
uses
  Math;
 
const
  MU = 3.986e5;
  G = 9.81;
 
type
  TData = array[1..4] of Extended;
  TOrbitMean = array[1..5] of Extended;
 
// Дифура
function Eq(T: Extended; StatVec: TData):TData;  overload;
function Eq(T: Extended; StatVec: TData; Lam,F: Extended):TData;  overload;
// Методы интегрирования
function EL(T,Step: Extended; StatVec: TData): TData;  overload;
function EL(T,Step: Extended; StatVec: TData; Lam,F: Extended): TData;  overload;
 
function RK(T,Step: Extended; StatVec: TData): TData;  overload;
function RK(T,Step: Extended; StatVec: TData; Lam,F: Extended): TData;  overload;
// Действия с матрицами
function N_M(N: Extended; M: TData): TData;
function M_M(M1, M2: TData): TData;
// Нахождение элементов орбиты
function OrbitMean(T: Integer; StatVec: TData): TOrbitMean;
 
implementation
// Дейсвия с матрицами
 
// Число на вектор-строку
function N_M(N: Extended; M: TData): TData;
var
  i: Integer;
begin
    for i:= 1 to 4 do
    begin
      Result[i]:= N*M[i];
    end;
end;
 
// Массив плюс Массив
function M_M(M1, M2: TData): TData;
var
  i: Integer;
begin
    for i:= 1 to 4 do
    begin
      Result[i]:= M1[i]+M2[i];
    end;
end;
 
// Дифура
 
// Без тяги
function Eq(T: Extended; StatVec: TData): TData; overload;
begin
   Result[1]:= StatVec[3];
   Result[2]:= StatVec[4]/StatVec[1];
   Result[3]:= StatVec[4]*StatVec[4]/StatVec[1] - MU/(StatVec[1]*StatVec[1]);
   Result[4]:=  -StatVec[3]*StatVec[4]/StatVec[1];
end;
 
// С тягой
function Eq(T: Extended; StatVec: TData; Lam,F: Extended): TData; overload;
begin
  Result[1]:= StatVec[3];
  Result[2]:= StatVec[4]/StatVec[1];
  Result[3]:= StatVec[4]*StatVec[4]/StatVec[1] - MU/(StatVec[1]*StatVec[1]) + F*cos(Lam);
  Result[4]:= F*sin(Lam) - StatVec[3]*StatVec[4]/StatVec[1];
end;
 
// Методы итегрирования
 
// Метод Эйлера
function EL(T,Step: Extended; StatVec: TData): TData; overload;
begin
   Result:= M_M(StatVec,N_M(Step,Eq(T,StatVec)));
end;
 
function EL(T,Step: Extended; StatVec: TData; Lam,F: Extended): TData; overload;
begin
   Result:= M_M(StatVec,N_M(Step,Eq(T,StatVec,Lam,F)));
end;
 
 
// Метод Рунге-Кутты
function RK(T,Step: Extended; StatVec: TData): TData; overload;
var
  k1,k2,k3,k4,
  eq1,eq2,eq3,
  eq11,eq12,eq13,eq14: TData;
begin
  k1:= Eq(T,StatVec);
 
  eq1:= M_M(StatVec,N_M(0.5,k1));
  k2:= Eq(T+Step/2,eq1);
 
  eq2:=  M_M(StatVec,N_M(0.5,k2));
  k3:= Eq(T+Step/2,eq2);
 
  eq3:=  M_M(StatVec,k3);
  k4:= Eq(T+Step,eq3);
 
  eq11:= M_M(k1,N_M(2,k2));
  eq12:= M_M(k4,N_M(2,k3));
  eq13:= M_M(eq11,eq12);
  eq14:= N_M(Step/6,eq13);
 
  Result:= M_M(StatVec,eq14);
end;
 
function RK(T,Step: Extended; StatVec: TData; Lam,F: Extended): TData; overload;
var
  k1,k2,k3,k4,
  eq1,eq2,eq3,
  eq11,eq12,eq13,eq14: TData;
begin
  k1:= Eq(T,StatVec,Lam,F);
 
  eq1:= M_M(StatVec,N_M(0.5,k1));
  k2:= Eq(T+Step/2,eq1,Lam,F);
 
  eq2:=  M_M(StatVec,N_M(0.5,k2));
  k3:= Eq(T+Step/2,eq2,Lam,F);
 
  eq3:=  M_M(StatVec,k3);
  k4:= Eq(T+Step,eq3,Lam,F);
 
  eq11:= M_M(k1,N_M(2,k2));
  eq12:= M_M(k4,N_M(2,k3));
  eq13:= M_M(eq11,eq12);
  eq14:= N_M(Step/6,eq13);
 
  Result:= M_M(StatVec,eq14);
end;
 
 
// Нахождение элементов орбиты
 
function OrbitMean(T: Integer; StatVec: TData): TOrbitMean;
var
  CosAnom: Extended;
begin
    Result[1]:= 1/MU * 
    sqrt(
    (StatVec[1] * StatVec[4] * StatVec[4] - MU) * (StatVec[1] * StatVec[4] * StatVec[4] - MU)
    + 
    (StatVec[1]*StatVec[3] * StatVec[4]) * (StatVec[1] * StatVec[3] * StatVec[4])
    );
    
    CosAnom:= (StatVec[1] * StatVec[4] * StatVec[4] - 1)
    /
    (
    (StatVec[1] * StatVec[4] * StatVec[4] - MU) * (StatVec[1] * StatVec[4] * StatVec[4] - MU)
    +
    (StatVec[1] * StatVec[3] * StatVec[4]) * (StatVec[1] * StatVec[3] * StatVec[4])
    );
    
    Result[2]:= StatVec[1] * (1 - Result[1] * CosAnom);
    
    Result[3]:= Result[2] / (1 - Result[1] * Result[1]);
    
    Result[4]:= Result[2] / (1 + Result[1]);
    
    Result[5]:= Result[2] / (1 - Result[1]);
end;
end.

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
unit uMain;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, TeEngine, Series, ExtCtrls, TeeProcs, Chart, StdCtrls, ComCtrls,
  Math, uIntegrator;
 
type
 
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    GroupBox2: TGroupBox;
    PageControl1: TPageControl;
    TabSheet1: TTabSheet;
    TabSheet2: TTabSheet;
    PageControl2: TPageControl;
    PageControl3: TPageControl;
    TabSheet3: TTabSheet;
    TabSheet4: TTabSheet;
    TabSheet5: TTabSheet;
    TabSheet6: TTabSheet;
    TabSheet7: TTabSheet;
    TabSheet8: TTabSheet;
    TabSheet9: TTabSheet;
    TabSheet10: TTabSheet;
    TabSheet11: TTabSheet;
    TabSheet12: TTabSheet;
    TabSheet13: TTabSheet;
    Memo1: TMemo;
    Memo2: TMemo;
    Chart1: TChart;
    Chart2: TChart;
    Chart3: TChart;
    Chart4: TChart;
    Chart5: TChart;
    Chart6: TChart;
    Chart7: TChart;
    Chart8: TChart;
    Chart9: TChart;
    Series1: TLineSeries;
    Series2: TLineSeries;
    Series3: TLineSeries;
    Series4: TLineSeries;
    Series5: TLineSeries;
    Series6: TLineSeries;
    Series7: TLineSeries;
    Series8: TLineSeries;
    Series9: TLineSeries;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Label5: TLabel;
    Label6: TLabel;
    Label7: TLabel;
    Label8: TLabel;
    Label10: TLabel;
    Label11: TLabel;
    ComboBox1: TComboBox;
    Label12: TLabel;
    Button1: TButton;
    Button2: TButton;
    Edit1: TEdit;
    Edit2: TEdit;
    Edit3: TEdit;
    Edit4: TEdit;
    Edit5: TEdit;
    Edit6: TEdit;
    Edit8: TEdit;
    Edit9: TEdit;
    Label9: TLabel;
    Label13: TLabel;
    Edit7: TEdit;
    Edit10: TEdit;
    Label14: TLabel;
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
 
  end;
  
var
  Form1: TForm1;
 
implementation
 
{$R *.dfm}
 
 
procedure TForm1.Button1Click(Sender: TObject);
var
  a,m,s,h: Integer;
  e,w,u,lam: Extended;
 
  T1,p,anom,r,Vr,Vt,F: Extended;
  T0,T,i: Integer;
 
  EqM: TData;
  OrM: TOrbitMean;
 
  Iter: Extended;
  Iter1: Integer;
 
  PowerOn, PowerOff,
  TOn, TOff: Integer;
begin
  a:= StrToInt(Edit1.Text);
  e:= StrToFloat(Edit2.Text);
 
  w:= StrToFloat(Edit3.Text);
  w:= DegToRad(w);
 
  u:= StrToFloat(Edit4.Text);
  u:= DegToRad(u);
 
  m:= StrToInt(Edit5.Text);
  s:= StrToInt(Edit6.Text);
 
  lam:= StrToFloat(Edit8.Text);
  lam:= DegToRad(lam);
 
  h:= StrToInt(Edit9.Text);
 
  T0:= 0;
  T1:= 2*pi*a*sqrt(a/mu);
  T:= Round(T1);
 
    p:= a*(1-sqr(e));
    anom:= u-w;
    r:= p/(1-e*cos(anom));
    Vr:= sqrt(mu/p)*sin(anom);
    Vt:= sqrt(mu/p)*(1-e*cos(anom));
    F:= s/m;
 
  EqM[1]:= r;
  EqM[2]:= u;
  EqM[3]:= Vr;
  EqM[4]:= Vt;
 
  Iter:= (T-T0)/h;
  Iter1:= Round(Iter);
 
  if ComboBox1.ItemIndex = 0 then
  begin
    Memo1.Clear;
    Memo2.Clear;
    Series1.Clear;
    Series2.Clear;
    Series3.Clear;
    Series4.Clear;
    Series5.Clear;
    Series6.Clear;
    Series7.Clear;
    Series8.Clear;
    Series9.Clear;
 
    Memo1.Lines.Add('Выбран метод Эйлера');
    Memo2.Lines.Add('Выбран метод Эйлера');
      for i:= 0 to Iter1 do
      begin
         T0:= T0+h;
 
         if i < PowerOff then
           begin
             if i >= PowerOn then
                EqM:= EL(T0, h, EqM, lam, F)
             else
                EqM:= EL(T0, h, EqM);
           end
         else
            EqM:= EL(T0, h, EqM);
            
            
         Memo1.Lines.Add('r('+ IntToStr(T0) +')= '+ FloatToStr(EqM[1])+';  '+
                         'U('+ IntToStr(T0) +')= '+ FloatToStr(EqM[2])+';  '+
                         'Vr('+ IntToStr(T0) +')= '+ FloatToStr(EqM[3])+';  '+
                         'Vt('+ IntToStr(T0) +')= '+ FloatToStr(EqM[4])+';'
                        );
 
         Series1.AddXY(T0,EqM[1]);
         Series2.AddXY(T0,EqM[2]);
         Series3.AddXY(T0,EqM[3]);
         Series4.AddXY(T0,EqM[4]);
 
         OrM:= OrbitMean(T0, EqM);
 
         Memo2.Lines.Add('e(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[1]) + ';  ' +
                         'p(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[2]) + ';  ' +
                         'a(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[3]) + ';  ' +
                         'Rп(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[4]) + ';  ' +
                         'Rа(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[5]) + ';  '
                        );
 
         Series5.AddXY(T0,OrM[1]);
         Series6.AddXY(T0,OrM[2]);
         Series7.AddXY(T0,OrM[3]);
         Series8.AddXY(T0,OrM[4]);
         Series9.AddXY(T0,OrM[5]);
      end;
  end;
 
    if ComboBox1.ItemIndex = 1 then
  begin
    Memo1.Clear;
    Memo2.Clear;
    Series1.Clear;
    Series2.Clear;
    Series3.Clear;
    Series4.Clear;
    Series5.Clear;
    Series6.Clear;
    Series7.Clear;
    Series8.Clear;
    Series9.Clear;
 
    Memo1.Lines.Add('Выбран метод Рунге-Кутты');
    Memo2.Lines.Add('Выбран метод Рунге-Кутты');
      for i:= 0 to Iter1 do
      begin
         T0:= T0+h;
 
         if i < PowerOff then
           begin
             if i >= PowerOn then
                EqM:= RK(T0, h, EqM, lam, F)
             else
                EqM:= RK(T0, h, EqM);
           end
         else
            EqM:= RK(T0, h, EqM);
 
         Memo1.Lines.Add('r('+ IntToStr(T0) +')= '+ FloatToStr(EqM[1])+';  '+
                         'U('+ IntToStr(T0) +')= '+ FloatToStr(EqM[2])+';  '+
                         'Vr('+ IntToStr(T0) +')= '+ FloatToStr(EqM[3])+';  '+
                         'Vt('+ IntToStr(T0) +')= '+ FloatToStr(EqM[4])+';'
                        );
 
         Series1.AddXY(T0,EqM[1]);
         Series2.AddXY(T0,EqM[2]);
         Series3.AddXY(T0,EqM[3]);
         Series4.AddXY(T0,EqM[4]);
 
         OrM:= OrbitMean(T0, EqM);
 
         Memo2.Lines.Add('e(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[1]) + ';  ' +
                         'p(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[2]) + ';  ' +
                         'a(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[3]) + ';  ' +
                         'Rп(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[4]) + ';  ' +
                         'Rа(' + IntToStr(T0) + ') = ' + FloatToStr(OrM[5]) + ';  '
                        );
 
         Series5.AddXY(T0,OrM[1]);
         Series6.AddXY(T0,OrM[2]);
         Series7.AddXY(T0,OrM[3]);
         Series8.AddXY(T0,OrM[4]);
         Series9.AddXY(T0,OrM[5]);
      end;
  end;
 
end;
 
procedure TForm1.Button2Click(Sender: TObject);
begin
   Form1.Close;
end;
 
 
 
 
end.

Добавлено через 3 минуты
NotBeginner, как бы вот так

@Kitsani · 25.08.2014, 12:23

А сам проект с исхлодниками можете кинуть

@xSimplex · 25.08.2014, 14:43 **[ТС]**

Kitsani, Я там doc файл прикрепил в начале поста

@xSimplex · 26.08.2014, 21:01 **[ТС]**

Может, кто-то знает, что не так с графиком, почему он строиться так, когда должна быть окружность?

@Игорь[Igor] · 26.08.2014, 21:08

TChart Сортирует координаты, по X. Рисуй на канве

@FIL · 26.08.2014, 21:15

Delphi
1
Series1.XValues.Order := loNone;

или в Инспекторе Объектов: Series1 -> XValues -> Order -> loNone

@Игорь[Igor] · 26.08.2014, 21:16

Ух ты

@xSimplex · 26.08.2014, 21:32 **[ТС]**

Всем спасибо

Новые блоги и статьи Все статьи Все блоги /
Как писать чистый, тестируемый и качественный код на Python py-thonny 12.07.2025 Помню свой первый проект на Python. Работал тогда быстро, грязно, лишь бы работало. Код был похож на запутанный клубок - переменные по одной букве, функции на 200 строк, комментарии отсутствовали как. . .	Blazor и контроллер сервопривода IoT Meadow Maple Wired 11.07.2025 Я решил разобраться, как можно соединить современные веб-технологии с миром "железа". Интересная комбинация получилась из Blazor в качестве веб-интерфейса и микроконтроллера Meadow с его веб-сервером. . .	Генерация OpenQASM из кода Q# EggHead 10.07.2025 Летом 2024-го я начал эксперименты с библиотекой Q# Bridge, и знаете что? Она оказалась просто находкой для тех, кто работает на стыке разных квантовых экосистем. Основная фишка этой библиотеки -. . .	Изучаем новый шаблон ИИ-чата .NET AI Chat Web App stackOverflow 10.07.2025 В . NET появилось интересное обновление - новый шаблон ИИ-чата под названием . NET AI Chat Web App. Когда я впервые наткнулся на анонс этого шаблона, то сразу понял, что Microsoft наконец-то. . .	Результаты исследования от команды ARP (июль 2025 г.) Programma_Boinc 10.07.2025 Результаты исследования от команды ARP (июль 2025 г. ) Африканский проект по дождям (ARP) World Community Grid снова запущен! Мы рады поделиться обновленной информацией о нашем прогрессе с осени. . .
Angular vs Svelte - что лучше? Reangularity 09.07.2025 Сегодня рынок разделился на несколько четких категорий: тяжеловесы корпоративного уровня (Angular), гибкие универсалы (React), прогрессивные решения (Vue) и новая волна компилируемых фреймворков. . .	Code First и Database First в Entity Framework UnmanagedCoder 09.07.2025 Entity Framework дает нам свободу выбора, предлагая как Code First, так и Database First подходы. Но эта свобода порождает вечный вопрос — какой подход выбрать? Entity Framework — это. . .	Как использовать Bluetooth-модуль HC-05 с Arduino Wired 08.07.2025 Bluetooth - это технология, созданная чтобы заменить кабельные соединения. Обычно ее используют для связи небольших устройств: мобильных телефонов, ноутбуков, наушников и т. д. Работает она на частоте. . .	Руководство по структурам данных Python AI_Generated 08.07.2025 Я отчетливо помню свои первые серьезные проекты на Python - я писал код, он работал, заказчики были относительно довольны. Но однажды мой наставник, взглянув на мою реализацию поиска по огромному. . .	Тестирование энергоэффективности и скорости вычислений видеокарт в BOINC проектах Programma_Boinc 08.07.2025 Тестирование энергоэффективности и скорости вычислений видеокарт в BOINC проектах Опубликовано: 07. 07. 2025 Рубрика: Uncategorized Автор: AlexA Статья размещается на сайте с разрешения. . .

Метод Эйлера и метод Рунге-Кутта: проверить код

Решение

Решение

@xSimplex 0 / 0 / 3 Регистрация: 24.07.2013 Сообщений: 51
	19.08.2014, 23:08 [ТС]
	Если кто-то читал, то могу сказать, что я еще раз проверил и у меня получаются слишком большие значения при расчете коэф. в рунге-кутте. Из-за чего это может быть? Ошибка в исходных ур-иях или я что-то криво подал?? 0

@xSimplex 0 / 0 / 3 Регистрация: 24.07.2013 Сообщений: 51
	24.08.2014, 04:41 [ТС]
	Сообщение было отмечено droider как решение Решение Ошибку нашел 0

@Kitsani 0 / 0 / 1 Регистрация: 17.02.2014 Сообщений: 12
	25.08.2014, 12:23
	А сам проект с исхлодниками можете кинуть 0

@xSimplex 0 / 0 / 3 Регистрация: 24.07.2013 Сообщений: 51
	25.08.2014, 14:43 [ТС]
	Kitsani, Я там doc файл прикрепил в начале поста 0

@xSimplex 0 / 0 / 3 Регистрация: 24.07.2013 Сообщений: 51
	26.08.2014, 21:01 [ТС]
	Может, кто-то знает, что не так с графиком, почему он строиться так, когда должна быть окружность? Миниатюры 0

@Игорь[Igor] 726 / 478 / 130 Регистрация: 24.12.2008 Сообщений: 3,924
	26.08.2014, 21:08
	TChart Сортирует координаты, по X. Рисуй на канве 0

@Игорь[Igor] 726 / 478 / 130 Регистрация: 24.12.2008 Сообщений: 3,924
	26.08.2014, 21:16
	Ух ты 0

@xSimplex 0 / 0 / 3 Регистрация: 24.07.2013 Сообщений: 51
	26.08.2014, 21:32 [ТС]
	Всем спасибо 0

Опции темы