Форум программистов, компьютерный форум CyberForum.ru

Ошибка в алгоритме Рунге-Кутты - C++

Восстановить пароль Регистрация
Другие темы раздела
C++ Как в файле перейти на новую строчку? http://www.cyberforum.ru/cpp-beginners/thread1004953.html
Как в файле перейти на новую строчку? Пробовал так fstream sc("Save\\1.txt");//Открыл файл sc.getline(infscore,'endl'); //Считал инфу до конца строки ......... lenscore=strlen(infscore)+1; //Жалкая попытка перейти sc.seekg(lenscore); // на новую строку
C++ Вызов функции проверки Доброго времени суток! 2-ой день мучаюсь с задачей. Есть массив прямоугольников, вершины которых я считал с файла. Задача состоит в том, чтобы проверить пересекаются ли они. Сделал соответствующую функцию : bool Intersects(Rect Obj1, Rect Obj2) { int ax,ay,ax1,ay1,bx,by,bx1,by1; ax = Obj1.ItsUpperLeftGetX(); ay = Obj1.ItsUpperLeftGetY(); ax1 = Obj1.ItsLowerRightGetX(); http://www.cyberforum.ru/cpp-beginners/thread1004940.html
Поправка программы C++
Есть вот такой вот код #include <iostream> #include <conio.h> #include <stdio.h> #include <math.h> using namespace std; int main () { system("COLOR 0A"); double ob, x1 , x2 , y1 , y2, p1 , p2 , pi , skolko , status;
C++ Считывание и запись определителя матрицы
в файле записан массив , его надо считать и найти его определитель, и записать его в другой файл. программа не хочет считывать данные. #include <stdio.h> #include <stdlib.h> #include <math.h> #include <iostream> #include <fstream> double **M(int m, int im, int jm, double ** arr)
C++ Массивы http://www.cyberforum.ru/cpp-beginners/thread1004904.html
. Даны действительные числа х1, ..., х17. Найти сумму значений |xi-xj| (1<= i< j <=17) (1 меньше или равно i меньше j меньше или равно 17)
C++ Сравнить элементы char написан код нужно рядом диогональ отсортировать пузырьковым методом код программы есть вот только не сортирует помогите: Добавлено через 34 секунды #include <cstdlib> #include <iostream> #include <string.h> #include <cstring> using namespace std; int main(){ подробнее

Показать сообщение отдельно
SpideR13
 Аватар для SpideR13
0 / 0 / 0
Регистрация: 02.06.2010
Сообщений: 7
11.11.2013, 22:06     Ошибка в алгоритме Рунге-Кутты
Здравствуйте! При выполнении курсовой работы по вычислению координат положения спутника ГЛОНАСС столкнулся с ошибкой, что неправильно вычисляются координаты. Числа похожие, но не те. Есть рабочая программа написанная на языке Pascal. Никак не могу найти, чем мой код не соответствует коду, который рабочий. Просьба помочь найти ошибку. Код C++ полностью готов к компилированию, код Pascal не мой и был написан в Borland Delphi, привожу только необходимую часть кода.

Код C++
C++
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
#include <math.h>
#include <iostream>
#include <stdio.h>
using namespace std;
 
// Задаём исходные данные для расчёта
int Ns[4] = {36, 45, 47, 54}, // номера спутников
Litera[4] = {12, 6, 4, 10}, // литеры частота
Day[4] = {1031, 1031, 1031, 1031}; // номер суток внутри четырехлетнего интервала от високосного года
double Tk[4] = {56640, 56640, 56640, 56640}, // показания часов спутника в момент начала излучения кадра
Tb[4] = {56700, 56700, 56700, 56700}, // МДВ для заданных ниже координат X, Y, Z; составляющих вектора скорости Xt, Yt, Zt и составляющих ускорения Xtt, Ytt, Ztt
X[4] = {2.3958496093750000e+004, -3.9609389648437500e+003, 1.7244073242187500e+004, -1.1319433105468750e+004},
Y[4] = {-6.1954052734375000e+002, 1.8279998535156250e+004, -9.5787133789062500e+003, 1.5059692382812500e+003},
Z[4] = {8.7343925781250000e+003, 1.7395937500000000e+004, 1.6202982910156250e+004, 2.2819439453125000e+004},
Xt[4] = {1.2139244079589840e+000, -1.3721923828125000e+000, -1.4764909744262700e+000, -1.2335987091064450e+000},
Yt[4] = {-3.2507896423339840e-002, 1.9431800842285160e+000, 1.6627349853515620e+000, -2.8666810989379880e+000},
Zt[4] = {-3.3318939208984370e+000, -2.3511285781860350e+000, 2.5435619354248050e+000, -4.2124843597412110e-001},
Xtt[4] = {0.0000000000000000e+000, 0.0000000000000000e+000, 0.0000000000000000e+000, -9.3132257461547850e-010},
Ytt[4] = {0.0000000000000000e+000, 0.0000000000000000e+000, 0.0000000000000000e+000, -9.3132257461547850e-010},
Ztt[4] = {-1.8626451492309570e-009, -1.8626451492309570e-009, -9.3132257461547850e-010, -9.3132257461547850e-010},
Gn[4] = {7.2759576141834260e-012, 0.0000000000000000e+000, -4.5474735088646410e-012, 0.0000000000000000e+000}, // относительное отклонение прогнозируемого значения несущей
Tn[4] = {-3.4392718225717540e-004, -1.4806631952524190e-004, -7.0374459028244020e-005, 3.2387673854827880e-005}, // сдвиг шкалы времени спутника относительно часов системы ГЛОНАСС
Tc[4] = {-3.8743019104003910e-007, -3.8743019104003910e-007, -3.8743019104003910e-007, -3.8743019104003910e-007}, // поправка к шкале времени системы ГЛОНАСС к МДВ
PD[4] = {12.0723405158, 12.0715064653, 12.0716056073, 12.0760442303}, // псевдозадержки из MCA сообщений
Tp = 218960; // PBN
 
// Вводим необходимые для расчётов переменные
double Tj[4], Tbp[4], h, s[6][4], arg[6][4], dSj[6], k1[6], k2[6], k3[6], k4[6];
 
// Необходимые константы
const double mu = 398600.44, Re = 6378.136, C20 = -1082.63e-6, we = 0.7292115e-4, c = 299792458;
 
// Задаем функцию abs
double abs(double Chislo)
{
    int m;
    if (Chislo < 0) m = -1; 
    else m = 1;
    Chislo*= m;
    return Chislo;
}
 
// Процедура вычисления функции f[i], где i - номер используемой функции
double solvef(int p, int m)
{
    double r, A, f[6];
        f[0] = arg[3][m];
        f[1] = arg[4][m];
        f[2] = arg[5][m];
        r = sqrt ( arg[0][m]*arg[0][m] +  arg[1][m]*arg[1][m] + arg[2][m]*arg[2][m] );
        A = mu/(r*r*r);
        f[3] = (we*we - A)*arg[0][m] + 2*we*arg[4][m] + 1.5*C20 * (mu*Re*Re)/(r*r*r*r*r) * arg[0][m] * (1 - ( (5*arg[2][m]*arg[2][m])/(r*r) )) + Xtt[m];
        f[4] = (we*we - A)*arg[1][m] + 2*we*arg[3][m] + 1.5*C20 * (mu*Re*Re)/(r*r*r*r*r) * arg[1][m] * (1 - ( (5*arg[2][m]*arg[2][m])/(r*r) )) + Ytt[m];
        f[5] = (-1)*A*arg[2][m] + 1.5*C20 * (mu*Re*Re)/(r*r*r*r*r) * arg[2][m] * (3 - ( (5*arg[2][m]*arg[2][m])/(r*r) )) + Ztt[m];
        return f[p];
}
 
// Главная программа
int main()
{
    int i, j, n;
    
    // задаём вектор с наименованием выводимого числа
    string coord[6] = {"  X", "  Y", "  Z", " Xt", " Yt", " Zt"};
 
    // цикл для вычисления показаний спутниковых часов на момент предшествия
    for (j = 0; j < 4; j++) {
        Tj[j] = Tp + 10800 - 2*86400 - PD[j];
        
        cout<<"№ "<<Ns[j]<<endl;
        
        cout<<" Tj = ";
        if ( Tj[j] > 0 ) cout<<" ";
        printf("%.20f\n", Tj[j]);
 
        // вычисляем МДВ для спутника
        Tbp[j] = Tj[j] + Tc[j] + Tn[j] - Gn[j]*(Tj[j] - Tb[j]);
 
        cout<<"Tbp = ";
        if ( Tbp[j] > 0 ) cout<<" ";
        printf("%.20f\n", Tbp[j]);
 
        // задаем вектор вычисляемых значений координат спутника и составляющих Vx s[i][j], и arg[i][j]; где i - номер значения (e.g. 0 - X, 1 - Y и т.п.), j - номер спутника
        s[0][j]=X[j];
        s[1][j]=Y[j];
        s[2][j]=Z[j];
        s[3][j]=Xt[j];
        s[4][j]=Yt[j];
        s[5][j]=Zt[j];
        
        // цикл с while с вычислением конечного h
        double Tbc = Tb[j];
        if ( Tbc < Tbp[j] ) h = 10; else h = -10;
        n = 0;
        while ( n == 0 ) {
            if ( abs(Tbp[j] - Tbc) < abs(h) ) {
                n = 1;
                h = Tbp[j] - Tbc;
            }            
            Tbc+= h; // вычисление Tbc для цикла while //
        
            for ( i = 0; i < 6; i++ ) arg[i][j] = s[i][j];
            for ( i = 0; i < 6; i++ ) k1[i] = solvef(i, j);
            for ( i = 0; i < 6; i++ ) k1[i] = k1[i]*h;
            for ( i = 0; i < 6; i++ ) arg[i][j] = s[i][j] + k1[i]/2;
            for ( i = 0; i < 6; i++ ) k2[i] = solvef(i, j);
            for ( i = 0; i < 6; i++ ) k2[i] = k2[i]*h;
            for ( i = 0; i < 6; i++ ) arg[i][j] = s[i][j] + k2[i]/2;
            for ( i = 0; i < 6; i++ ) k3[i] = solvef(i, j);
            for ( i = 0; i < 6; i++ ) k3[i] = k3[i]*h;
            for ( i = 0; i < 6; i++ ) arg[i][j] = s[i][j] + k3[i];
            for ( i = 0; i < 6; i++ ) k4[i] = solvef(i, j);
            for ( i = 0; i < 6; i++ ) k4[i] = k4[i]*h;
            for ( i = 0; i < 6; i++ ) dSj[i] = (k1[i] + 2*k2[i] + 2*k3[i] + k4[i])/6;
            for ( i = 0; i < 6; i++ ) s[i][j]+= dSj[i];
            } // конец цикла while
        // вывод s[i][j]
        for ( i = 0; i < 6; i++ ) {
            cout<<coord[i]<<" = ";
            if ( s[i][j] > 0 ) cout<<" ";
            printf("%.20f\n", s[i][j]);
        }
        cout<<endl;
    } // конец цикла по спутникам
return 0;
}


Код Pascal
Pascal
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
//f6 - процедура-функция, вычисляющая правые части системы дифференциальных уравнений
    function f6 (ArrayofArg: array of extended): TMyArray;
      var
      r,A:extended;
 
      begin
         Result[1]:=arg[m,4];
         Result[2]:=arg[m,5];
         Result[3]:=arg[m,6];
            r:=sqrt(arg[m,1]*arg[m,1]+arg[m,2]*arg[m,2]+arg[m,3]*arg[m,3]);
            A:=(mu)/(r*r*r);
         Result[4]:=(wz*wz-A)*arg[m,1]+2*wz*arg[m,5]+
         +1.5*C20*((mu*Re*Re)/(IntPower(r,5)))*arg[m,1]*(1-((5*arg[m,3]*arg[m,3])/(r*r)))+xtt[m];
         Result[5]:=(wz*wz-A)*arg[m,2]-2*wz*arg[m,4]+
         +1.5*C20*((mu*Re*Re)/(IntPower(r,5)))*arg[m,2]*(1-((5*arg[m,3]*arg[m,3])/(r*r)))+ytt[m];
         Result[6]:=-A*arg[m,3]+
         +1.5*C20*((mu*Re*Re)/(IntPower(r,5)))*arg[m,3]*(3-((5*arg[m,3]*arg[m,3])/(r*r)))+ztt[m];
      end;
 
{ В Ы Ч И С Л Е Н И Я   П О   А Л Г О Р И Т М У
(вычисление коор-т и составляющих векторов скорости всех спутников) } 
 
//нумерация пунктов согласно выданному материалу практического задания
procedure TForm1.Button2Click(Sender: TObject);
//процедура обработки нажатия кнопки "Вычисление результатов"
  var
    Trtm2,h,hnew,ti:extended;
    Hs: array of array of extended;//градиентная матрица Hs (двумерный массив)
    Sr,Nev,Popr,tau: array [1..4] of extended;//вектора приближений, невязок, поправок и задержек
    PR:integer;
    Matrix, Matrix2 : MatrixPtr;//структура типа "матрица" - исп-ся для поиска обратной матрицы)
  begin
    h:=StrToFloat(Edit1.Text);
    if (h<=0) or (h>60) then
      begin
        Showmessage('Введен некорректный шаг h');
        Exit;
      end;
    //определение 6-мерного вектора S (двумерный массив):
    setlength(S,5);
    for m:=1 to 4 do
      begin
        setlength(S[m],7);
    //для удобства индексации нулевой элемент массива не используем, поэтому 7, а не 6 элементов
      end;
    {определение 6-мерного вектора arg, компоненты которого используются как аргументы для
    вычисления функций, стоящий в правых частях системы дифф-ых уравнений: arg[i,j],
    где i - порядковый номер спутника (1-4), j - номер аргумента (1-6):}
    setlength(arg,5);
    for m:=1 to 4 do
      begin
        setlength(arg[m],7);
      end;
    //определение 6-мерного вектора Ds:
    setlength(Ds,5);
    for m:=1 to 4 do
      begin
        setlength(Ds[m],7);
      end;
    //начало цикла по четырем спутникам (m=1..4):
    for m:=1 to 4 do
      begin
        //1.Вычисление показаний спутниковых часов Tjtjprec на момент предшествия tjprec:
        Trtm2:=Trtm+10800-2*86400;//+3 часа и вычли целое числое суток (2), чтоб остаток был меньше 86400с
        Tjtjprec[m]:=Trtm2-PD[m];//результат
        //2. Вычисление показаний часов МДВ TMDVtjprec на момент предшествия tjprec:
        TMDVtjprec[m]:=Tjtjprec[m]+Tc[m]+Tn[m]-Gn[m]*(Tjtjprec[m]-Tb[m]);
        //3. Вычисление координат x,y,z и составляющих вектора скорости vx,vy,vz m-го спутника с момента,
        //когда показания часов МДВ равны Tb, на момент предшествия tjprec, когда показания этих же часов
        //равны TMDVtjprec - численное интегрирование системы дифференциальных уравнений:
        //заполнение начальных значений компонент вектора S:
        S[m,1]:=x[m];
        S[m,2]:=y[m];
        S[m,3]:=z[m];
        S[m,4]:=xt[m];
        S[m,5]:=yt[m];
        S[m,6]:=zt[m];
        ti:=Tb[m];//текущий момент времени при интегрировании
        //непосредственно итерационный цикл интегрирования
        while (abs(TMDVtjprec[m]-ti)>h) do
          begin
            //1 строка (алгоритм содержит 10 строк - блоку присваивается номер строки алгоритма,
            //которой он соответствует)
            for k:=1 to 6 do
              begin
                arg[m,k]:=S[m,k];
              end;
            //2 строка
            K1:=f6(arg[m]);
            for k:=1 to 6 do
              begin
                K1[k]:=K1[k]*h;
              end;
            //3 строка
            for k:=1 to 6 do
              begin
                arg[m,k]:=S[m,k]+0.5*K1[k];
              end;
            //4 строка
            K2:=f6(arg[m]);
            for k:=1 to 6 do
              begin
                K2[k]:=K2[k]*h;
              end;
            //5 строка
            for k:=1 to 6 do
              begin
                arg[m,k]:=S[m,k]+0.5*K2[k];
              end;
            //6 строка
            K3:=f6(arg[m]);
            for k:=1 to 6 do
              begin
                K3[k]:=K3[k]*h;
              end;
            //7 строка
            for k:=1 to 6 do
              begin
                arg[m,k]:=S[m,k]+K3[k];
              end;
            //8 строка
            K4:=f6(arg[m]);
            for k:=1 to 6 do
              begin
                K4[k]:=K4[k]*h;
              end;
            //9 строка
            for k:=1 to 6 do
              begin
                Ds[m,k]:=(1/6)*(K1[k]+2*K2[k]+2*K3[k]+K4[k]);
              end;
            //10 строка - вычисление (i+1)-ых значений для компонент вектора s
            for k:=1 to 6 do
              begin
                S[m,k]:=S[m,k]+Ds[m,k];
              end;
            ti:=ti+h;
          end;//конец цикла интегрирования по итерациям. Осталась последняя точка (TMDVtjprec)
        //последняя итерация с новым шагом:
        if (abs(TMDVtjprec[m]-ti)<=h) then
          begin
            hnew:=abs(TMDVtjprec[m]-ti);
            //1 строка
            for k:=1 to 6 do
              begin
                arg[m,k]:=S[m,k];
              end;
            //2 строка
            K1:=f6(arg[m]);
            for k:=1 to 6 do
              begin
                K1[k]:=K1[k]*hnew;
              end;
            //3 строка
            for k:=1 to 6 do
              begin
                arg[m,k]:=S[m,k]+0.5*K1[k];
              end;
            //4 строка
            K2:=f6(arg[m]);
            for k:=1 to 6 do
              begin
                K2[k]:=K2[k]*hnew;
              end;
            //5 строка
            for k:=1 to 6 do
              begin
                arg[m,k]:=S[m,k]+0.5*K2[k];
              end;
            //6 строка
            K3:=f6(arg[m]);
            for k:=1 to 6 do
              begin
                K3[k]:=K3[k]*hnew;
              end;
            //7 строка
            for k:=1 to 6 do
              begin
                arg[m,k]:=S[m,k]+K3[k];
              end;
            //8 строка
            K4:=f6(arg[m]);
            for k:=1 to 6 do
              begin
                K4[k]:=K4[k]*hnew;
              end;
            //9 строка
            for k:=1 to 6 do
              begin
                Ds[m,k]:=(1/6)*(K1[k]+2*K2[k]+2*K3[k]+K4[k]);
              end;
            //10 строка - вычисление последнего значения вектора S на интервале интегрирования
            for k:=1 to 6 do
              begin
                S[m,k]:=S[m,k]+Ds[m,k];
              end;
          end;//закончили весь итерационный процесс
        //выводим вычисленные коор-ты и сост-ие вектора скорости спутников в таблицу результатов
        for k:=1 to 6 do
          begin
            StringGrid2.Cells[k,m]:=FloatToStr(S[m,k]);
          end;
      end;//конец цикла по m=1..4 спутникам
// далее идут другие вычисления, которые никак не влияют на результат выше


Вывод результата программы C++
Код
№ 36
 Tj =  56947.92765948420128552243
Tbp =  56947.92731516778439981863
  X =  24245.55882723983449977823
  Y = -621.96948941025243584590
  Z =  7902.08078379823564318940
 Xt =  1.10182357112657003384
 Yt =  0.01224006963147607728
 Zt = -3.38143983822364857517

№ 45
 Tj =  56947.92849353470228379592
Tbp =  56947.92834508095256751403
  X = -4290.30597361219133745180
  Y =  18745.20837423249759012833
  Z =  16800.38049823398614535108
 Xt = -1.28530772213540189775
 Yt =  1.80974429849611628818
 Zt = -2.45256645089889735445

№ 47
 Tj =  56947.92839439270028378814
Tbp =  56947.92832363193883793429
  X =  16875.63835070742425159551
  Y = -9167.71195700164935260545
  Z =  16821.51713725931404042058
 Xt = -1.49538996039706617580
 Yt =  1.65231698870472998308
 Zt =  2.44545297703395414501

№ 54
 Tj =  56947.92395576970011461526
Tbp =  56947.92398776994377840310
  X = -11631.64480262850884173531
  Y =  788.91681057900791529391
  Z =  22698.20635907495307037607
 Xt = -1.28507509437282840814
 Yt = -2.91752988306132898799
 Zt = -0.55661749205162081022


Вывод программы Pascal
Ошибка в алгоритме Рунге-Кутты


Результаты соответственно должны совпадать. На Delphi точно написано правильно, результаты сравниваются с теми, что у преподавателя. Входные данные для получения данных результатов ввожу такие же, проверял неоднократно.
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
 
Текущее время: 03:53. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2016, vBulletin Solutions, Inc.
Рейтинг@Mail.ru