движение самолета в Delphi

[cMep4[21]]

Здравствуйте!
Мне нужно смоделировать движение сверхзвукового самолета. Я составил программу в Delphi. У меня все считает, но я не смог сделать, чтобы вычисления производились с заданной точностью (проблемы со счетчиком итераций). Не получилось вывести все решения в таблицу. (Может выводить в какое-нибудь другое место?).
Решения в Delphi мне надо проверить в MathCAD'е. Там я тоже че-то намудрил ). Пожалуйста помогите исправить, доделать программу.
Еще есть файл в паскале, но там ошибка "Program1.pas".
К теме прикладываю программу в делфи, файлы маткада, скан книжки -теории.
Спасибо.

[cMep4[21]]

var
  Form1: TForm1;
  V0,//начальная скорость
  V1,//скорость V1
  z0,//начальная высота
  zmax,//конечная высота
  theta0,//начальный угол
  theta1,// угол (1)
  W,//вес самолета. В проекте используется британская система мер
  S,//площадь несущих крыльев
  Cd,//коэффициент сопротивления воздуха
  Cl,//Коэффициент подъемной силы
  h,//шаг изменения значения высоты
  E,//точность
  grav,//ускорение свободного падения
  D,//Сила сопротивления воздуха
  L,//Подъемная сила
  T,//Тяга двигателей
  ro0,//нормальная плотность воздуха
  DiffV,//Ошибка ограничения DiffV
  DiffA,//Ошибка ограничения DiffA
  A0,//начальный угол в радианах
  F,// FF(V0,A0) - первое уравнение движения при значениях V0 и A0
  G,// GF(V0,A0) - второе уравнение движения при значениях V0 и A0
  VV,// скорость
  AA,// угол в радианах
  A1,//угол (1) в радианах
  A2,//угол (2) в радианах
  V20,//предсказанная скорость (2)
  A20,//предсказанный угол (2)
  itn,// счетчик количества итераций
  V2new,// скорректированная скорость (2)
  A2new,// скорректированный угол в радианах (2)
  theta,// угол
  rho,// плотность
  z1,// высота (1)
  z,//  высота
  z2,// высота (1)
  theta2,//угол (2)
  V2:real;//скорость (2)
  m,i:integer;// количество строк таблицы
implementation

{$R *.dfm}



  procedure TForm1.Button1Click(Sender: TObject);
function DenstyF(z:Double):Double;
begin
  DenstyF:=0.0023769*exp((-z)/23600);//плотность воздуха,зависит от высоты
  end;
  function ThrustF(V:Double;z:Double):Double;
begin
ThrustF:=(1+0.5*exp(-sqr(0.001*V-1.5)))*(10000-0.27*z);//тяга двигателей
end;
  function FF(V:Double;A:Double):Double;
  begin
  FF:=(grav*(T-(W*sin(A))-0.5*Cd*rho*S*sqr(V)))/(W*V*sin(A));// Первое уравнение движения самолета. rho-плотность
  end;
function GF(V:Double;A:Double):Double;
begin
GF:=(grav*(0.5*Cl*rho*S*sqr(V)-(W*cos(A))))/(W*sqr(V)*sin(A)); // Второе уравнение движения самолета.
end;
begin
grav:=32.174;// ускорение свободного падения. Используется британская система мер
    z0:=strtofloat(edit1.Text); //начальная высота
    zmax:=strtofloat(edit2.Text);// конечная высота
    theta0:=strtofloat(edit3.Text); //начальный угол
    h:=strtofloat(edit4.Text);  // шаг изменения высоты
    E:=strtofloat(edit5.Text);  // точность
W:=20000;// вес самолета (все величины в английской системе)
S:=160; //  площадь крыльев
Cd:=0.02;// коэффициент сопротивления воздуха
Cl:=0.1; // коэффициент подъемной силы
V0:=1500; //начальная скорость
theta0:=15; //начальный угол
    z1:=z0+h;
    rho:=DenstyF(z0);
    T:=ThrustF(V0,z0);
    A0:=theta0/57.29578;//в радианы для внутренних вычислений
    F:=FF(V0,A0);
    G:=GF(V0,A0);
    VV:=V0+h*F;
    AA:=A0+H*G;
    rho:=DenstyF(z1);
    T:=ThrustF(VV,z1);
    V1:=V0+0.5*h*(F+FF(VV,AA));//Исправленный метод Эйлера
    A1:=A0+0.5*h*(G+GF(VV,AA))
    theta:=A1*57.29578; //перевод в градусы
    theta1:=A1*57.29578;
     edit6.Text:=floattostr(V1); 
     edit7.Text:=floattostr(z1);
     edit8.Text:=floattostr(theta);
   z2:=z1+h;
   rho:=DenstyF(z1);
   T:=ThrustF(V1,z1);
   F:=FF(V1,A1);
   G:=GF(V1,A1);
   V2:=V0+2*h*F;
   A2:=A0+2*h*G;
   V20:=V2;
   A20:=A2;
   rho:=DenstyF(z2);
  T:=ThrustF(V2,z2);
   V2new:=V1+0.5*h*(F+FF(V2,A2));
   A2new:=A1+0.5*h*(G+GF(V2,A2));
   theta:=57.29578*A2;
   while (abs((V2-V2new)/V2new)<e) and (abs((A2-A2new)/A2new)<e) do
begin
end;
if (itn<10) then
V2:=V2new;
A2:=A2new;
T:=ThrustF(V2,z2);
   V2new:=V1+0.5*h*(F+FF(V2,A2));//скорректированное значение 
   A2new:=A1+0.5*h*(G+GF(V2,A2));
   theta:=57.29578*A2;
begin
end;
V2:=V2new;
A2:=A2new;
DiffV:=0.2*(V20-V2new);
DiffA:=0.2*(A20-A2new);
V2:=V2new+DiffV;
A2:=A2new+DiffA;
theta:=57.29578*A2;
DiffA:=57.29578*DiffA;
     edit9.Text:=floattostr(V2);
     edit10.Text:=floattostr(z2);
     edit11.Text:=floattostr(theta);
     edit13.Text:=floattostr(DiffV);
     edit14.Text:=floattostr(DiffA);
     m:=round((zmax-z)/h)+2;
  StringGrid1.RowCount:=m; //количество строк
  StringGrid1.Cells[0,0]:='Высота';
  StringGrid1.Cells[1,0]:='Скорость';
  StringGrid1.Cells[2,0]:='Угол';
  StringGrid1.Cells[3,0]:='Ошибка ограничения DiffV';
  StringGrid1.Cells[4,0]:='Ошибка ограничения DiffA';
  StringGrid1.Cells[5,0]:='Счетчик итераций';
  for i:=1 to m do
   begin
   StringGrid1.Cells[0,1]:=FloatToStr(z0);
   StringGrid1.Cells[1,1]:=FloatToStr(V0);         {высота,скорость,угол,.......}
   StringGrid1.Cells[2,1]:=FloatToStr(theta0);
   StringGrid1.Cells[0,2]:=FloatToStr(z1);
   StringGrid1.Cells[0,3]:=FloatToStr(z2);
   StringGrid1.Cells[1,2]:=FloatToStr(V1);
   StringGrid1.Cells[1,3]:=FloatToStr(V2);
   StringGrid1.Cells[2,2]:=FloatToStr(theta1);
   StringGrid1.Cells[2,3]:=FloatToStr(theta);
   StringGrid1.Cells[3,3]:=FloatToStr(DiffV);
   StringGrid1.Cells[4,3]:=FloatToStr(DiffA);
end;
end;
end.

[Rin]

У меня все считает, но я не смог сделать, чтобы вычисления производились с заданной точностью (проблемы со счетчиком итераций).

Если правильно понимаю, то нужно повторять вычисление, пока не достигните определенного эпсилон, что будет выполняться условие |F1-F2|< epsilon. Если так, то просто введите epsilon="нужная вам точность". Нужна точность с 10 значащими цифрами после запятой? Epsilon=0.0000000001 . Проверяете на условие, которое написано выше и все.

Не получилось вывести все решения в таблицу. (Может выводить в какое-нибудь другое место?).

Проблема, кажется, в этом:

  for i:=1 to m do
   begin
   StringGrid1.Cells[0,1]:=FloatToStr(z0);
   StringGrid1.Cells[1,1]:=FloatToStr(V0);         {высота,скорость,угол,.......}
   StringGrid1.Cells[2,1]:=FloatToStr(theta0);
   StringGrid1.Cells[0,2]:=FloatToStr(z1);
   StringGrid1.Cells[0,3]:=FloatToStr(z2);
   StringGrid1.Cells[1,2]:=FloatToStr(V1);
   StringGrid1.Cells[1,3]:=FloatToStr(V2);
   StringGrid1.Cells[2,2]:=FloatToStr(theta1);
   StringGrid1.Cells[2,3]:=FloatToStr(theta);
   StringGrid1.Cells[3,3]:=FloatToStr(DiffV);
   StringGrid1.Cells[4,3]:=FloatToStr(DiffA);
end;

[cMep4[21]]

Спасибо за ответ.
По поводу эпсилон: по условию задачи e=0,0001. Я почему-то сделал, что его вводит пользователь. Мне нужно, чтобы вычисления были с точностью e=0,0001 и при этом программа считала количество итераций и выводила.
Это практическое задание из книги "Мак-Кракен Д. , Дорн У. Численные методы и программирование на ФОРТРАНе", 1977 год. Там текст программы приведен на ФОРТРАНе. У меня не полностью получилось "переделать" этот код под Delphi.
Мне нужно примерно как на скриншоте вывести результаты.
По поводу вывода результатов: это я сам вручную написал первые 3 результата из 40. Как мне сделать так, чтобы программа сама заполняла таблицы? (Может есть способ полегче, чем выводить в таблицу StrihgGrid?)

[Rin]

Чтобы программа сама заполняла таблицу Тебе нужно вместо цифры писать переменную "m" :

for m:= 1 to 'сколько-то там' do
begin
// функция вычисления полета

  StringGrid1.Cells[0,m]:='переменная под высоту';
  StringGrid1.Cells[1,m]:='переменная под скорость';
  StringGrid1.Cells[2,m]:='переменная под угол';
  StringGrid1.Cells[3,m]:='переменная под ошибку ограничения DiffV';
  StringGrid1.Cells[4,m]:='переменная под ошибку ограничения DiffA';
  StringGrid1.Cells[5,m]:='переменная под счетчик итераций';

[cMep4[21]]

Да, действительно) Таблица заполняется). Спасибо.
А как мне выводить не только первые значения, а все от первого до i-того?
i получается при таком выражении:

i:=(zmax-z0)/h

[Rin]

Оформите вычисление i-тых значений полета как процедуру. И запускайте её в цикле. Таким образом у вас будет вычисляться I-ые параметры полета и записываться рез-ты в stringgrid. И так по кругу столько раз, сколько вам нужно.

for m:= 1 to 'сколько-то там' do
begin
// процедура вычисления полета

  StringGrid1.Cells[0,m]:='переменная под высоту';
  StringGrid1.Cells[1,m]:='переменная под скорость';
  StringGrid1.Cells[2,m]:='переменная под угол';
  StringGrid1.Cells[3,m]:='переменная под ошибку ограничения DiffV';
  StringGrid1.Cells[4,m]:='переменная под ошибку ограничения DiffA';
  StringGrid1.Cells[5,m]:='переменная под счетчик итераций';
end;

[cMep4[21]]

Если честно, я не понимаю как написать процедуру

[s-andriano]

Но Вы понимаете, что она должна делать?

[cMep4[21]]

Сначала надо найти значения V1 и theta1(угол) исправленным методом Эйлера. Потом методом прогноза и коррекции найти V2 и theta2. Это я сделал.
Вычисления следующих основано на известных предыдущих значениях. Например: находим V1 исправленным методом Эйлера .Чтобы найти V2 используем V1. Чтобы найти V3 используем V2, и т.д.
Нужно как-то сделать, чтобы программа вычисляла эти значения и выводила.
Я пока не знаю куда выводить. Надо сделать, чтобы вывелось первое значение, сохранилось, туда-же вывелось второе.
А вычисление следующих можно заменять:
z0:=z1;
z1:=z2;
v0:=v1;
v1:=v2;
theta1:=theta2;
Чтобы так заменять, надо сначала куда-то выводить надо.
Как и куда выводить я пока не знаю.