Форум программистов, компьютерный форум, киберфорум
Delphi для начинающих
Войти
Регистрация
Восстановить пароль
 
Рейтинг 4.95/19: Рейтинг темы: голосов - 19, средняя оценка - 4.95
0 / 0 / 3
Регистрация: 24.07.2013
Сообщений: 51
1

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

18.08.2014, 15:10. Показов 3740. Ответов 11
Метки нет (Все метки)

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

Спасибо всем

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.
0
Вложения
Тип файла: rar Пример.rar (244.3 Кб, 75 просмотров)
Тип файла: docx Задание.docx (33.9 Кб, 60 просмотров)
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
18.08.2014, 15:10
Ответы с готовыми решениями:

Графики - Метод Рунге-Кутта и Метод Адамса
Пожалуйста помогите!!! очень срочно нужно сделать графики (два в одном): метод Рунге-Кутта и метод...

Графики - Метод Рунге-Кутта и Метод Адамса
Вот задание:

Метод Рунге-Кутта
Ребят, помогите разобраться почему полученный данные не выводятся в указанный файл(если что, то это...

Метод Рунге-Кутта
Здравствуйте, пожалуйста помогите в программе Delphi решить систему обыкновенных дифференциальных...

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

Решение

Ошибку нашел
0
174 / 161 / 71
Регистрация: 22.02.2013
Сообщений: 1,770
Записей в блоге: 2
24.08.2014, 13:28 4
Цитата Сообщение от xSimplex Посмотреть сообщение
Ошибку нашел
выкладывай рабочий код
0
0 / 0 / 3
Регистрация: 24.07.2013
Сообщений: 51
24.08.2014, 19:27  [ТС] 5
Вот тут все работает. Можете кстати проверить, вдруг что-то упустил, прошу проверить сами методы. Там мне еще надо провести пару манипуляций для сохранения значений в файл и т.п. , ну главное что скелет есть, а мясо нарастет

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, как бы вот так
0
0 / 0 / 1
Регистрация: 17.02.2014
Сообщений: 9
25.08.2014, 12:23 6
А сам проект с исхлодниками можете кинуть
0
0 / 0 / 3
Регистрация: 24.07.2013
Сообщений: 51
25.08.2014, 14:43  [ТС] 7
Kitsani, Я там doc файл прикрепил в начале поста
0
0 / 0 / 3
Регистрация: 24.07.2013
Сообщений: 51
26.08.2014, 21:01  [ТС] 8
Может, кто-то знает, что не так с графиком, почему он строиться так, когда должна быть окружность?
0
Миниатюры
Метод Эйлера и метод Рунге-Кутта: проверить код  
723 / 475 / 130
Регистрация: 24.12.2008
Сообщений: 3,924
26.08.2014, 21:08 9
TChart Сортирует координаты, по X. Рисуй на канве
0
Модератор
3475 / 2599 / 740
Регистрация: 19.09.2012
Сообщений: 7,966
26.08.2014, 21:15 10
Лучший ответ Сообщение было отмечено xSimplex как решение

Решение

Delphi
1
Series1.XValues.Order := loNone;
или в Инспекторе Объектов: Series1 -> XValues -> Order -> loNone
1
723 / 475 / 130
Регистрация: 24.12.2008
Сообщений: 3,924
26.08.2014, 21:16 11
Ух ты
0
0 / 0 / 3
Регистрация: 24.07.2013
Сообщений: 51
26.08.2014, 21:32  [ТС] 12
Всем спасибо
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
26.08.2014, 21:32

Заказываю контрольные, курсовые, дипломные работы и диссертации здесь.

Метод Рунге-Кутта
Здравствуйте, подскажите пожалуйста где я ошибся тут только для два значения икса и игрека...

Метод Рунге-Кутта 4-го порядка
Добрый день. Скачал программку, которая решает систему диф.уравнений методом Рунге-Кутты 4-го...

Метод Рунге-Кутта 4-го порядка
Доброго времени суток Господа. Прошу помочь с написанием кода для интегрирования сист. диф. ура-ий...

Метод Рунге-Кутта для 5 порядка
нужна ваша помощь! Нужно написать программу для Рунге-Кутта для дифференциального уравнения 5...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
12
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.