Информатика: Движение тела в среде
Материал из Department of Theoretical and Applied Mechanics
Версия от 00:12, 17 января 2016; 178.162.33.3 (обсуждение)
Описание программы: программа записывает в четыре файла результаты вычисления:
- Координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- Координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости;
- Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
- Координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
Скачать можно тут.
Визуализированный результат работы программы
- o1 - координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- o2 - координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости;
- o3 - координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
- o4 - координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
Для тела с массой 10,сопротивлением воздуха 1, угол бросания 30°, начальная скорость 30 м/с, ускорение свободного падения 9.8 м/c^2;
Примечание: графики o1 и o2 намеренно посчитаны с малой точностью, чтобы графики не сливались.
Файл "main.cpp"
1 #include <iostream>
2 #include <math.h>
3 #include "Vector.h"
4 #include <cstring>
5 #include <cmath>
6 #include <malloc.h>
7 #include <fstream>
8
9 using namespace std;
10
11 int n = 100;
12 ofstream outfile;
13
14 class Ball //класс бросаемого тела
15 {
16 private:
17 double angle,m,k; //угол броска,масса,коэффицент сопротивления воздуха
18 Vector3D r,v,a; //радиус-вектор,вектор скорости,ускорения
19 public:
20
21 //задание начальных параметров через угол,начальное положение,скорость и ускорение,с которым движется тело. Без сопротивления воздуха
22 Ball(double _angle, Vector3D _r, Vector3D _v, Vector3D _a)
23 {
24 angle = _angle;
25 r = _r;
26 v = _v;
27 a = _a;
28 }
29
30 //задание начальных параметров через угол,начальное положение,скорость и ускорение,с которым движется тело. Без сопротивления воздуха
31 Ball(double _angle, double _m, double _k, Vector3D _r, Vector3D _v, Vector3D _a)
32 {
33 angle = _angle;
34 r = _r;
35 v = _v;
36 a = _a;
37 m = _m;
38 k = _k;
39 }
40
41 //точная формула зависимости координаты от времени
42 Vector3D positionReal(double t)
43 {
44 double c1 = m/k,c2 = fabs(a.y)*c1, c3 = exp(-t/c1), c4 = c2*t;
45 return MakeVector(v.x*c1*(1 - c3), c1*(v.y + c2)*(1 - c3) - c4 , 0 );
46 }
47
48 //вывод положения на экран
49 void writePosToScreen()
50 {
51 cout << r.x << " " << r.y << " " << r.z << endl;
52 }
53
54 //вывод положения в файл
55 void writePosToFile(char s[])
56 {
57 outfile.open(s,ios :: app);
58 outfile << r.x << " " << r.y << endl;
59 outfile.close();
60 }
61
62 //вывод произвольного вектора на экран
63 void WVTS(Vector3D v)
64 {
65 cout.width(15);
66 cout << v.x;
67 cout.width(15);
68 cout << v.y;
69 cout.width(15);
70 cout << v.z << endl;
71 }
72
73 //вывод произвольного вектора в файл
74 void WVTF(Vector3D v,char s[])
75 {
76 outfile.open(s,ios :: app);
77 outfile << v.x << " " << v.y << endl;
78 outfile.close();
79 }
80
81 //"пересчет" координаты по Верле(Линейная зависмость)
82 void changeR(Vector3D r1, double dt)
83 {
84 r = MakeVector(2 * r.x - r1.x - k/m*v.x*dt*dt,2*r.y - r1.y - (abs(a.y) + k/m*v.y)*dt*dt, 0 );
85 }
86
87 //"пересчет" координаты по Верле(Квадратичная зависимость)
88 void changeRSQ(Vector3D r1, double dt)
89 {
90 r = MakeVector(2 * r.x - r1.x - k/m*Length(v)*v.x*dt*dt,2*r.y - r1.y - (abs(a.y) + k/m*Length(v)*v.y)*dt*dt, 0 );
91 }
92 //пересчет скорости по Верле
93 void changeV(Vector3D r1,double dt)
94 {
95 v =VS((VmV(r,r1)),1/(2*dt));
96 }
97
98 //рассчет предыдущегт к 0ому элементу
99 Vector3D MR1(double dt)
100 {
101 return MakeVector(r.x - v.x * dt,r.y - v.y * dt,0);
102 }
103
104 //возращает координату тела
105 Vector3D getR()
106 {
107 return r;
108 }
109
110 //рассчет времени полета
111 double TimeOfFly()
112 {
113 return (2*Length(v)*sin(angle)/Length(a));
114 }
115
116 //рассчет координаты по точной формуле. без сопротивления воздуха.
117 Vector3D position(double t)
118 {
119 return MakeVector(r.x + v.x*t + a.x*t*t/2,r.y + v.y*t + a.y*t*t/2,r.z + v.z*t + a.z*t*t/2);
120 }
121
122 };
123
124 int main()
125 {
126 //задание начальных параметров
127 Vector3D g = {0,-9.8,0};
128 double a,dt = 0;
129 char s[20];
130
131 // cin >> dt;
132
133 dt = 0.1;
134 a = (M_PI * 30)/180;
135 Ball b1(a, MakeVector(0,0,0),MakeVector(30,a),g);
136
137 double tof = b1.TimeOfFly()+1; //единичка прибавлена,чтобы график красивым был
138
139 //Без сопротивления возлуха
140 strcpy(s,"");
141 strcat(s, "o1.txt");
142 outfile.open(s, ios :: trunc);
143 outfile.close();
144 for (double i = 0; i <= tof; i += dt)
145 {
146 b1.WVTS(b1.position(i));
147 b1.WVTF(b1.position(i), s);
148 }
149
150
151 //Верле(Линейная зависимость)
152 dt = 0.1;
153 a = (M_PI * 30)/180;
154 Ball b2(a,10 , 1, MakeVector(0,0,0),MakeVector(30,a),g);
155
156 strcpy(s,"");
157 strcat(s, "o2.txt");
158 outfile.open(s,ios :: trunc);
159 outfile.close();
160 Vector3D r1 = b2.MR1(dt),rp;
161 for (double i = 0; i <= 20; i += dt)
162 {
163 rp = b2.getR();
164 b2.writePosToFile(s);
165 b2.writePosToScreen();
166 b2.changeR(r1,dt);
167 b2.changeV(r1,dt);
168 r1.x = rp.x;
169 r1.y = rp.y;
170 }
171
172 //Точное решение (Линейная зависимость)
173 dt = 0.1;
174 a = (M_PI * 30)/180;
175 Ball b3(a,10 , 1, MakeVector(0,0,0),MakeVector(30,a),g);
176
177 strcpy(s,"");
178 strcat(s, "o3.txt");
179 outfile.open(s, ios :: trunc);
180 outfile.close();
181 for (double i = 0; i <= 20; i += dt)
182 {
183 b3.WVTS(b3.positionReal(i));
184 b3.WVTF(b3.positionReal(i), s);
185 }
186
187
188 //Верле (Квадратичная зависимость)
189 dt = 0.1;
190 a = (M_PI * 30)/180;
191 Ball b4(a,10 , 1, MakeVector(0,0,0),MakeVector(30,a),g);
192
193 strcpy(s,"");
194 strcat(s, "o4.txt");
195 outfile.open(s, ios :: trunc);
196 outfile.close();
197 r1 = b4.MR1(dt);
198 for (double i = 0; i <= 20; i += dt)
199 {
200 rp = b4.getR();
201 b4.writePosToFile(s);
202 b4.writePosToScreen();
203 b4.changeRSQ(r1,dt);
204 b4.changeV(r1,dt);
205 r1.x = rp.x;
206 r1.y = rp.y;
207 }
208
209 return 0;
210 }
Файл "Vector.h"
1 #ifndef VECTOR_H_INCLUDED
2 #define VECTOR_H_INCLUDED
3
4 struct Vector3D
5 {
6 double x,y,z;
7 };
8
9 Vector3D VmV(Vector3D v1,Vector3D v2) //векторное вычитание
10 {
11 Vector3D v = {v1.x - v2.x,v1.y - v2.y,v1.z - v2.z };
12 return v;
13 };
14 Vector3D VpV(Vector3D v1,Vector3D v2) //векторное сложение
15 {
16 Vector3D v = {v1.x + v2.x,v1.y + v2.y,v1.z + v2.z };
17 return v;
18 }
19
20 double VV(Vector3D v1,Vector3D v2) //скалярное умножение
21 {
22 return (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);
23 }
24
25 Vector3D VxV(Vector3D v1,Vector3D v2) //векторное умножение
26 {
27 Vector3D v = {v1.y*v2.z - v1.z*v2.y, v1.z*v2.x - v1.x*v2.z,v1.x*v2.y - v1.y*v2.x};
28 return v;
29 }
30
31 bool Kol(Vector3D v1,Vector3D v2)
32 {
33 return ((v1.x/v2.x == v1.y/v2.y)&&(v1.z/v2.z == v1.y/v2.y))? true:false;
34 }
35
36 Vector3D VS(Vector3D v1, double s)
37 {
38 Vector3D v = {v1.x*s, v1.y*s, v1.z*s};
39 return v;
40 }
41
42 double Length(Vector3D v1)
43 {
44 return sqrt(VV(v1,v1));
45 }
46
47 Vector3D MakeVector(double x,double y,double z)
48 {
49 Vector3D v = {x,y,z};
50 return v;
51 }
52
53 Vector3D MakeVector(double length,double angle)
54 {
55 Vector3D v = {length * cos(angle), length * sin(angle),0};
56 return v;
57 }
58
59 double Proection(Vector3D base, Vector3D dir)
60 {
61 return (VV(base,dir)/Length(base));
62 }
63 #endif // VECTOR_H_INCLUDED
Описание программы: пользователь вводит начальную скорость полета, угол бросания и шаг, с которым будут рассчитаны точки.
Программа записывает в один файл результаты вычисления:
- Координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
- Координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
- Координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости;
Скачать можно тут.
Визуализированный результат работы программы
Для тела с массой 1 кг,сопротивлением воздуха 0.05, угол бросания 30°, начальная скорость 30 м/с, ускорение свободного падения 9.8 м/c^2, шаг 0.01;
- "Zapis.txt" using 1 : 2 - координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- "Zapis.txt" using 3 : 4 - координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
- "Zapis.txt" using 5 : 6 - координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости;
- "Zapis.txt" using 7 : 8 - координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости.
1 #include <iostream>
2 #include <locale.h>
3 #include <math.h>
4 #include <fstream>
5 #include<iomanip>
6 #include <cmath>
7
8 using namespace std;
9
10 class fly ///создаем класс полета тела
11 {
12
13 private: ///объявляем тип переменных в привате
14 double Vo, Agrad, Brad, time, step, amountdouble; ///Vo-начальная скорость тела; Agrad-угол, под которым летит тело, в градусах;
15 ///Brad-угол, под которым летит тело, в радианах; time-время полета тела; step-шаг;
16 ///amountdouble-количество точек (типа double)
17
18 public: ///объявляем переменные в паблике
19 int amount; ///amoun-количество точек (типа int)
20
21 fly (double _step, double _Agrad, double _Vo):step(_step),Agrad(_Agrad),Vo(_Vo) ///создаем конструктор функции с объявлением переменных
22 {
23
24 double g=9.8; ///объявляем тип и значение переменной g(ускорение свободного падения)
25 Brad=3.14159*Agrad/180.0; ///переводим значение угла из градусов в радианы
26
27 time=2*Vo*sin(Brad)/g; ///рассчитываем время полета тела
28 amountdouble=(round(time/step)+1); ///подсчитываем количество точек с заданым шагом
29 amount=static_cast<int>(amountdouble); ///преобразуем количество из типа double к типу int
30
31 }
32 void zapis (char Zapis[]) ///создаем функцию записи
33 {
34 double g=9.8, m=1, n=0.05; ///объявляем тип и значения переменных g (ускорение свободного падения), m (масса тела), n(коэффициэнт сопротивления)
35 double Xb, Yb; ///объявляем тип переменных для полёта тела без сопротивления ветра Xb(координата тела по Х), Yb(координата тела по У)
36 double V=Vo, Vxv=Vo*cos(Brad), Vyv=Vo*sin(Brad), Xo=0, Yo=0, Xv=0, Yv=0, Y_1=Yo-Vo*sin(Brad)*step, X_1=Xo-Vo*cos(Brad)*step, Y=0, X=0;
37 ///объявляем тип переменных для метода ВерлеI V (скорость тела по модулю), Vxv (составляющая скорости по Х),
38 ///Vyv (составляющая скорости по У), Xo (начальное положение тела на оси Х), Yo (начальное положение тела на оси У),
39 ///Xv (координата тела на оси Х), Yv (координата тела на оси У), Y_1 (координата тела на (n-1)ом шаге на оси У),
40 ///X_1 (координата тела на (n-1)ом шаге на оси Х), Y (координата тела на n-ом шаге на оси У),
41 ///X (координата тела на n-ом шаге на оси Х);
42 double Vxv2=Vo*cos(Brad), Vyv2=Vo*sin(Brad), Xo2=0, Yo2=0, Xv2=0, Yv2=0, Y_12=Yo2-Vo*sin(Brad)*step, X_12=Xo-Vo*cos(Brad)*step, Y2=0, X2=0;
43 ///объявляем тип переменных для метода ВерлеII V (скорость тела по модулю), Vxv (составляющая скорости по Х),
44 ///Vyv (составляющая скорости по У), Xo (начальное положение тела на оси Х), Yo (начальное положение тела на оси У),
45 ///Xv (координата тела на оси Х), Yv (координата тела на оси У), Y_1 (координата тела на (n-1)ом шаге на оси У),
46 ///X_1 (координата тела на (n-1)ом шаге на оси Х), Y (координата тела на n-ом шаге на оси У),
47 ///X (координата тела на n-ом шаге на оси Х);
48 double Yt=0, Xt=0, Yot=0, Xot=0, Voxt=Vo*cos(Brad), Voyt=Vo*sin(Brad);
49
50 ofstream outfile("Zapis.txt"); ///запись элементов функции в фаил "Zapis.txt"
51 outfile<<setw(20)<<"Xb"<<setw(20)<<"Yb"<<setw(20)<<"Xt"<<setw(20)<<"Yt"<<setw(20)<<"Xv"<<setw(20)<<"Yv"<<setw(20)<<"Xv2"<<setw(20)<<"Yv2"<<" \n"; ///вывод на экран по столбцам
52 ///X (координата тела на оси Х без ветра),
53 ///Y (координата тела на оси У без ветра),
54 ///Xv (координата тела на оси Х с ветром для метода Верле),
55 ///Yv (координата тела на оси У с ветром для метода Верле)
56 ///setw() размер столбиков
57
58 for (int l=0; Yb>=0; ++l) ///создаем цикл от 0 до тех пор пока У больше нуля
59 {
60 outfile<<setw(20)<<Xb<<setw(20)<<Yb<<setw(20)<<Xt<<setw(20)<<Yt<<setw(20)<<Xv<<setw(20)<<Yv<<setw(20)<<Xv2<<setw(20)<<Yv2<<" \n";
61 ///вывод на экран по столбцам Xv, Yv;
62
63 ///полёт без ветра
64 Xb=Vo*cos(Brad)*l*step;
65 Yb=Vo*sin(Brad)*l*step-(9.8*l*step*l*step*0.5);
66
67 ///точный метод
68 Xt=Xot+(m/n)*Voxt*(1.0 - exp((-n*l*step)/m));
69 Yt=Yot+(m/n)*(Voyt+g*(m/n))*(1.0 - exp((-n*l*step)/m))-g*l*step*(m/n);
70
71 ///метод Верле I
72 Xv=2*X-X_1-(n/m)*V*Vxv*step*step; ///расчитываем координату Х в момент времени t для метода Верле
73 Yv=2*Y-Y_1-(g+(n/m)*V*Vyv)*step*step; ///расчитываем координату У в момент времени t для метода Верле
74 Vxv=(Xv-X_1)/(2.0*step); ///расчитываем скорость тела по оси Х в момент времени t для метода Верле
75 Vyv=(Yv-Y_1)/(2.0*step); ///расчитываем скорость тела по оси У в момент времени t для метода Верле
76 V=sqrt(Vxv*Vxv+Vyv*Vyv); ///рассчитываем скорость тела по модулю
77 X_1=X; ///присваиваем значению координаты Х на (n-1)ом шаге значение координаты Х на n-ом шаге
78 X=Xv; ///присваиваем значению координаты Х на n-ом шаге значение координаты Х
79 Y_1=Y; ///присваиваем значению координаты У на (n-1)ом шаге значение координаты У на n-ом шаге
80 Y=Yv; ///присваиваем значению координаты У на n-ом шаге значение координаты У
81
82 ///метод Верле II
83 Xv2=2*X2-X_12-(n/m)*Vxv2*step*step; ///расчитываем координату Х в момент времени t для метода Верле
84 Yv2=2*Y2-Y_12-(g+(n/m)*Vyv2)*step*step; ///расчитываем координату У в момент времени t для метода Верле
85 Vxv2=(Xv2-X_12)/(2.0*step); ///расчитываем скорость тела по оси Х в момент времени t для метода Верле
86 Vyv2=(Yv2-Y_12)/(2.0*step); ///расчитываем скорость тела по оси У в момент времени t для метода Верле
87 X_12=X2; ///присваиваем значению координаты Х на (n-1)ом шаге значение координаты Х на n-ом шаге
88 X2=Xv2; ///присваиваем значению координаты Х на n-ом шаге значение координаты Х
89 Y_12=Y2; ///присваиваем значению координаты У на (n-1)ом шаге значение координаты У на n-ом шаге
90 Y2=Yv2; ///присваиваем значению координаты У на n-ом шаге значение координаты У
91
92 }
93
94 outfile.close();
95
96 }
97
98 };
99
100 int main()
101 {
102
103 setlocale(LC_ALL,"RUS"); ///функция, позволяющая с++ распознавать русский язык
104
105 double Vo, Agrad, Brad, time, step; ///объявляем тип переменных Vo (начальная скорость тела), Agrad (угол, под которым летит тело, в градусах);
106 ///Brad (угол, под которым летит тело, в радианах); time (время полета тела); step (шаг)
107 cout<<"Задайте начальную скорость тела в м/с: Vo="; ///на экран выводится сообщение с просьюой задать начальную скорость тела
108 cin>>Vo; ///пользователь вводит начальную скорость тела
109 cout<<'\n'<<"Задайте в градусах угол, под которым брошено тело (угол должен принимать значения от 0 до 90): a=";
110 ///на экран выводится сообщение с просьбой задать угол, под которым летит тело, в градусах
111 cin>>Agrad; ///пользователь вводит угол, под которым летит тело
112 cout<<'\n'<<"Задайте шаг (шаг должен быть очень маленьким): шаг="; ///на экран выводится сообщение с просьбой ввести шаг
113 cin>>step; ///пользователь вводит шаг
114
115 fly X(step,Agrad,Vo); ///объявление коструктора, создание функции Х, зависящей от step,Agrad,Vo
116 X.zapis("координаты.txt"); ///запись элементов функции в файл
117
118 return 0; ///конец программы
119
120 }
Описание программы: пользователь вводит начальную скорость полета, угол бросания и шаг, с которым будут рассчитаны точки.
Программа записывает в один файл результаты вычисления:
- Координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- Координаты, полученные при численном интегрировании - метод Эйлера;
- Координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
- Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
Скачать можно тут.
Визуализированный результат работы программы
Для тела с массой 0.5 кг,сопротивлением воздуха 0.1, угол бросания 30°, начальная скорость 30 м/с, ускорение свободного падения 9.8 м/c^2, шаг 0.001;
- "output.txt" using 1 : 2 - координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- "output.txt" using 3 : 4 - координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости (численное интегрирование - метод Эйлера);
- "output.txt" using 5 : 6 - координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости;
- "output.txt" using 7 : 8 - координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости.
1 #include <iostream>
2 #include <fstream>
3 #include <iomanip>
4 #include <time.h>
5 #include <conio.h>
6 #include <stdlib.h>
7 #include <math.h>
8 #include <cstring>
9 using namespace std;
10
11 class pad ///создаем класс
12 {
13 private: ///в закрытом доступе
14 double *X, *Y, *E, *G, *Z, *S, *V, *U; ///координаты по х и у; *E, *G : численное интегрирование метод Эйлера, *Z, *S- метод верле, *V, *U- точный метод
15 double a, g , pi;///ускорение, коэфф свободного падения, значение числа пи
16 int Size; ///размер массива- сколько точек считать в первом способе
17 double v, shag, b, vx, vy;///скорость, шаг по времени, угол в градусах скорость по х, скорость по у
18
19 void SetX() ///создаем функцию для вычисления значенй по х для простого падения и по эйлеру
20 {
21 g = 9.8;///коэфф свободного падения
22 float t = 0; ///время
23 X = new double [Size];///создаем массив для координат по х для простого падения
24 E = new double [Size];///создаем массив для координат по х для интегрирования по эйлеру
25 X[0] = 0; ///задаем значение по х в нуле
26 E[0] = 0 ; ///задаем значение по х в нуле по эйлеру
27 for (int i = 1; i < Size; i++) ///задаем цикл от 1 (для нуля мы задали), чтобы считать координаты
28 {
29 t += shag; ///каждый раз прибавляем по времени шаг
30 X[i] = v * cos(a) * t; ///координаты по х для простого падения
31 E[i] = E[i-1] + v * cos(a) * shag; ///х из интегрирования по эйлеру
32 }
33 }
34
35 void SetY()///создаем функцию для вычисления значенй по у для простого падения и по эйлеру
36 {
37 g = 9.8; ///коэфф свободного падения
38 double Vy; /// переменная для значение скорости по у для метода эйлера
39 float t = 0; ///время
40 Y = new double [Size];///создаем массив для координат по у для простого падения
41 G = new double [Size];///создаем массив для координат по х для интегрирования по эйлеру
42 Vy = v * sin (a);/// значение скорости по у для метода эйлера
43 Y[0] = 0;///задаем значение по у в нуле
44 G[0] = 0;///задаем значение по у в нуле по эйлеру
45 for (int i = 1; i < Size; i++)///задаем цикл от 1 (для нуля мы задали), чтобы считать координаты
46 {
47 t += shag; ///каждый раз прибавляем по времени шаг
48 Y[i] = v * sin (a) *t - (g * t * t) * 0.5; ///координаты по у для простого падения
49 Vy -= g * shag; ///значение скорости по у для метода эйлера
50 G[i] = G[i-1] + Vy * shag;///у из интегрирования по эйлеру
51 }
52 }
53
54 void SetVerle() ///функция для метода верле
55 {
56 double k = 0.1, m = 0.5; ///коэфф сопротивления водуха, масса тела
57 g = 9.8; /// коэфф свободного падения
58 uint32_t Size1 = 1000000.0; ///размер массива
59 S = new double [Size1]; ///создаем массив для значений по х для метода верле
60 Z = new double [Size1]; ///создаем массив для значений по у для метода верле
61 vx = v * cos(a); ///формулы для вычисления скорости по оси х
62 vy = v * sin(a); ///формулы для вычисления скорости по оси у
63 S[1] = 0; ///значение х метод верле
64 S[0] = -vx * shag; ///значение в нуле
65 Z[1] = 0; ///значение у метод верле
66 Z[0] = -vy * shag; ///значение в нуле
67 for (int i = 0; i < Size1-2; i++) ///задаем цикл
68 {
69 S[i+2] = 2.0 * S[i+1] - S[i] - (k / m) * v * vx * shag * shag;///значения по х для верле
70 vx = 0.5 * ( 1.0 / shag )* ( S[i+2] - S[i]);///считаем значения скорости по оси х
71 Z[i+2] = 2.0 * Z[i+1] - Z[i] - ( g + (k / m) * v * vy ) * shag * shag;///значения по х для верле
72 vy = 0.5 * ( 1.0 / shag )* ( Z[i+2] - Z[i]);///считаем значения скорости по оси х
73 v = sqrt (vx * vx + vy * vy); ///модуль общей скорости
74 }
75 }
76 void SetVerleLast() ///функция для точного метода верле
77 {
78 double k = 0.1, m = 0.5;///коэфф сопротивления водуха, масса тела
79 g = 9.8; /// коэфф свободного падения
80 uint32_t Size2 = 1000000.0; ///размер массива
81 float t = 0; ///время
82 V = new double [Size2]; ///создаем массив для значений по х для точного метода верле
83 U = new double [Size2]; ///создаем массив для значений по у для точного метода верле
84 vx = v * cos(a); ///формулы для вычисления скорости по оси х
85 vy = v * sin(a); ///формулы для вычисления скорости по оси у
86 ///double e = 2.7 ;///значение экспоненты
87 V[0] = 0; ///значение х точный метод верле
88 U[0] = 0; ///значение у точный метод верле
89 for (int i = 1; i < Size2; i++)
90 {
91 t += shag; ///увеличиваем время на шаг
92 V[i] = vx * (m / k) * (1.0 - exp(((-k) / m) * t)); ///значения по х для точного верле
93 U[i] = (m / k) * (vy + g * (m / k)) * (1.0 - (exp(((-k) / m) * t))) - g * t * (m / k);///значения по х для точного верле
94 }
95 }
96
97 public: ///в открытом
98 pad()
99 {
100 X = 0; ///зануляем значения
101 Y = 0;
102 Size = 0;
103 v = 0;
104 shag = 0;
105 b = 0;
106 }
107 pad(double _v, double _shag, double _b) ///конструктор с параметрами
108 {
109 pi = M_PI; ///значение числа пи
110 g = 9.8; ///коэфф свободного падения
111 v = _v;/// присваиваем значения переменных значению параметров в конструкторе
112 shag = _shag;
113 b = _b;
114 a = (pi * b) / 180.0 ; ///вычисляем значение угла в радианах
115 double t = (2.0 * v * sin(a)) / g; /// считаем значение времени
116 Size = abs( t / shag )+1;///ищем значение размера массива
117 SetX(); ///вызываем функции зависящие от параметров конструктора
118 SetY();
119 SetVerle();
120 SetVerleLast();
121 }
122
123 void FilePrint() ///функция записи в файл
124 {
125 ofstream fout("output.txt"); ///открываем файл уже созданный в папке с программой
126 fout << "X: " << " Y: " << " E: " << " G: " << " S: " << " Z: "<< " V: "<< " U: "<<"\n" ; ///выводим стоку с разными названиями массивов, соотв. координатам по х и у различных методов
127 for (int i = 0; i < Size; i++) ///цикл
128 fout << X[i] << " " << Y[i] << " " << E[i] << " " << G[i] << " " << S[i] << " " << Z[i] << " " << V[i] <<" "<< U[i] <<"\n"; ///забивает сами значения массивов
129 fout.close();///закрываем файл
130 };
131 };
132
133 int main()/// основная функция
134 {
135 double shag, b, v; ///шаг, угол в градусах, скорость начальная
136 cout << "vvedite v "; ///просим пользователя ввести значение скорости начальной
137 cin >> v; ///считываем начальную скорость
138 cout << "vvedite ygol ";///просим пользователя ввести угол в градусах
139 cin >> b;/// считываем угол
140 cout << "vvedite shag ";///просим пользователя ввести шаг по времени
141 cin >> shag; ///считываем значение шага
142 pad F1(v, shag, b); ///объявление коструктора, создание функции F1 с переменными v, shag, b
143 F1.FilePrint(); ///вызываем функцию для записи файла
144 }
Краткое описание алгоритма: в классе находятся координаты по формулам и записываются в файл.
Инструкция : Пользователь должен ввести начальную скорость, угол и шаг, с которым будут рассчитываться координаты. В файл координаты записываются в таком порядке: 1, 2 столбики - Координаты, рассчитанные по формуле, при движении без сопротивления воздуха; 3, 4 - Координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости; 5,6 - Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости; 7,8 - Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости.
Скачать можно тут.
1
Иванова Яна
Описание программы: в программе выполняются четыре метода подсчета координат тела, брошенного под углом к горизонту. Координаты записываются в файл, строятся четыре графика, иллюстрирующие поведение тела при полете. Код написан для определенных начальных условий (для примера), если Вы хотите выполнить расчет для другой конфигурации, внесите изменения в начальные данные программы в самом коде.
Начальная скорость: 40 м/с, угол бросания: 45 градусов, коэффициент сопротивления воздуха: 0.023, шаг по времени : 0.1 секунды.
Скачать программу можно здесь
Визуализированный результат работы программы
File:graph.png
1 #include <iostream>
2 #include <math.h>
3 #include <iomanip>
4 #include <fstream>
5 #include <conio.h>
6
7
8 using namespace std;
9
10 ofstream outfile;
11
12 double perevod (double angle) //перевод из градусов в радианы
13 {
14 return (angle * M_PI / 180 );
15 }
16
17
18 int main()
19 {
20 //объявление переменных и задание начальных значений
21 double X, Xnext, Xprev, Y, Ynext, Yprev, Vx, Vy, V,
22 m = 1 , dt = 0.1 , g = 9.8,t = 0,
23 ugol = 45, alpha, R = 0.023, Xo = 0, Yo = 0, Vo = 40;
24
25 alpha = perevod (ugol);
26
27 //точное решение для случая движения без сопротивления воздуха
28 Y = Yo;
29 X = Xo;
30
31 outfile.open("1.txt");
32
33 while (Y >= Yo)
34 {
35 X = Xo + Vo * cos(alpha) * t;
36 Vx = Vo * cos(alpha);
37 Y = Yo + Vo * sin(alpha) * t - 0.5 * g * t * t;
38 Vy = Vo * sin(alpha) - g * t;
39 t += dt;
40
41 outfile << X << ' ' << Y << endl;
42 }
43 outfile.close();
44
45 //начальные условия для квадратичной зависимости (метод Верле)
46 Yprev = Yo - Vo*sin(alpha)*dt;
47 Xprev = Xo - Vo*cos(alpha)*dt;
48 X = Xo;
49 Y = Yo;
50 V = Vo;
51 Vx = Vo * cos(alpha);
52 Vy = Vo * sin(alpha);
53
54 outfile.open("2.txt");
55
56 while (Y >= Yo)
57 {
58 Xnext = 2.0 * X - Xprev - (R / m) * V * Vx * (dt * dt);
59 Vx = ( Xnext - Xprev )/ (2.0 * dt);
60 Ynext = 2.0 * Y - Yprev - (g + (R / m) * V * Vy) * (dt * dt);
61 Vy = (Ynext - Yprev)/ (2.0 * dt);
62 V = sqrt(Vy*Vy + Vx*Vx );
63 outfile << X << ' ' << Y << endl;
64
65 Xprev = X;
66 X = Xnext;
67 Yprev = Y;
68 Y = Ynext;
69 }
70 outfile.close();
71
72 //начальные условия для линейной зависимости (метод Верле)
73 Yprev = Yo - Vo*sin(alpha)*dt;
74 Xprev = Xo - Vo*cos(alpha)*dt;
75 X = Xo;
76 Y = Yo;
77 V = Vo;
78 Vx = Vo * cos(alpha);
79 Vy = Vo * sin(alpha);
80
81 outfile.open("3.txt");
82
83 while (Y >= Yo)
84 {
85 Xnext = 2.0 * X - Xprev - (R / m) * Vx * (dt * dt);
86 Vx = ( Xnext - Xprev )/ (2.0 * dt);
87 Ynext = 2.0 * Y - Yprev - (g + (R / m) * Vy) * (dt * dt);
88 Vy = (Ynext - Yprev)/ (2.0 * dt);
89 V = sqrt(Vy*Vy + Vx*Vx );
90 outfile << X << ' ' << Y << endl;
91
92 Xprev = X;
93 X = Xnext;
94 Yprev = Y;
95 Y = Ynext;
96 }
97 outfile.close();
98
99 //точное решения для линейной зависимости
100 Y = Yo;
101 X = Xo;
102 t = 0;
103
104 outfile.open("4.txt");
105
106 while (Y >= Yo)
107 {
108 Vx = Vo * cos(alpha);
109 Vy = Vo * sin(alpha);
110 X = (m * Vx / R)* (1 - exp(-1 * R * t / m));
111 Y = (m/R)*((Vy + g * m / R)*(1 - exp(-1 * R * t / m)) - g * t);
112 t += dt;
113 outfile << X << ' ' << Y << endl;
114 }
115 outfile.close();
116
117 return 0;
118
119 }
Описание программы: программа записывает в четыре файла результаты вычисления:
- Координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- Координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости;
- Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
- Координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
Описание: Пользователя попросят ввести начальную скорость,угол бросания,массу тела и коэф.сопротивления воздуха,тогда программа запишет в 4 разных файла результаты следующих вычислений:
- Координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- Координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости;
- Координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
- Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
Графики полученные при скорости =10 m/c;угле = 30 градусам;массе=10 кг;коэф.сопротивления=1;
.png)
1 #include<iostream>
2 using namespace std;
3 #define N 1
4 #define PI 3.14159265
5 #include <fstream>
6 #include<cmath>
7 double g=9.8;
8 double step=0.01;
9 #include<math.h>
10
11 void Func(double v,double r)//1.Координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
12
13 {
14 double x,y;
15
16 ofstream fout;//открытие файла
17 fout.open("C:\\Users\\Light\\Desktop\\1.txt");//указываем путь записи
18 cout<<"method1"<<endl;
19
20 for(double t=0.01;t<N;t=t+0.01)
21 {
22 y=v*t*sin(r*PI/ 180)-g*t*t/2;//координата y
23 x=v*t*cos(r*PI / 180);//координата х
24
25 fout<<x<<" ";//запись в файл х
26 cout<<"X="<<x<<endl;//вывод на экран
27 fout<<y<<" ";//запись в файл
28 cout<<"Y="<<y<<endl<<endl;//вывод на экран
29 fout<<endl;
30 }
31 }
32
33 void Verle1( double n,double m ,double v0,double r)//Координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости;
34
35 {
36 double x,y,x0=0,y0=0,xn_1,yn_1;//x0,y0=Xn,Yn начальные значения;xn_1=X(n-1);yn_1=Y(n-1);
37
38 double vx=v0*cos(r*PI / 180);//рассчитваем Vn для первого случая n=0 для x
39 double vy=v0*sin(r*PI/ 180)-g*step;//рассчитваем Vn для первого случая n=0 для y
40
41 xn_1=x0-vx*step;//X(n-1) для первого случая n=0
42 yn_1=y0-vy*step;//Y(n-1) для первого случая n=0
43 ofstream fout;//открытие файла
44 fout.open("C:\\Users\\Light\\Desktop\\2.txt");//путь записи в файл
45 cout<<"Verle1"<<endl<<endl;
46 for(double t=0.02;t<N;t=t+step)
47 {
48 x=2*x0-xn_1-(n*vx*step*step)/m;//считаем Хn+1
49 vx=(x-xn_1)/(2*step);
50 xn_1=x0;//для следущего шага Xn-1=Xn
51 x0=x;//для следущего шага Xn=Xn+1
52
53 y=2*y0-yn_1-(g+(n*vy)/m)*step*step;//Yn+1
54 vy=(y-yn_1)/(2*0.01);//скорость
55 yn_1=y0;//для следущего шага Yn-1=Yn
56 y0=y;//для следущего шага Yn=Yn+1
57 cout<<"X="<<x<<endl;
58 cout<<"Y="<<y<<endl<<endl;
59 fout<<x<<" ";
60 fout<<y<<" ";
61 fout<<endl;
62
63
64 }
65 }
66
67 void Verle2( double n,double m ,double v0,double r)//3.Координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
68
69 {
70 double x,y,x0=0,y0=0,xn_1,yn_1,v;//x0,y0=Xn,Yn начальные значения;xn_1=X(n-1);yn_1=Y(n-1);
71
72 double vx=v0*cos(r*PI / 180);//рассчитваем Vn для первого случая n=0 для x
73 double vy=v0*sin(r*PI/ 180)-g*step;//рассчитваем Vn для первого случая n=0 для y
74
75 xn_1=x0-vx*step;//X(n-1) для первого случая n=0
76 yn_1=y0-vy*step;//Y(n-1) для первого случая n=0
77 ofstream fout;//открытие файла
78 fout.open("C:\\Users\\Light\\Desktop\\3.txt");//путь записи
79 cout<<"Verle2"<<endl<<endl;
80 for(double t=0.02;t<N;t=t+step)
81 {
82
83 v=sqrt(vx*vx+vy*vy);//скорость V
84
85 x=2*x0-xn_1-(n*v*vx*step*step)/m;//Xn+1
86 vx=(x-xn_1)/(2*step);//скорость Vx
87 xn_1=x0;//для следущего шага Xn-1=Xn
88 x0=x;//для следущего шага Xn=Xn+1
89
90 y=2*y0-yn_1-(g+(n*vy*v)/m)*step*step;//Yn+1
91 vy=(y-yn_1)/(2*0.01);//скорость Vy
92 yn_1=y0;//для следущего шага Yn-1=Yn
93 y0=y;//для следущего шага Yn=Yn+1
94 cout<<"X="<<x<<endl;
95 cout<<"Y="<<y<<endl<<endl;
96 fout<<x<<" ";
97 fout<<y<<" ";
98 fout<<endl;
99 }
100
101 }
102
103 void method4(double n,double m ,double v0,double r)//Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
104 {
105 double x,y,x0=0,y0=0;//x0,y0 начальные значения;
106
107 double vx=v0*cos(r*PI / 180);//рассчитваем Vx
108 double vy=v0*sin(r*PI/ 180)-g*step;//рассчитваем Vy
109
110 ofstream fout;//открытие файла
111 fout.open("C:\\Users\\Light\\Desktop\\4.txt");
112 cout<<"4"<<endl<<endl;
113 for(double t=0.01;t<N;t=t+step)
114 {
115 x=x0+m*vx*(1-(exp((-1*n)*t/m)))/n;//координата х
116
117 y = y0+(m/n)*((vy + g * m / n)*(1 - exp(-1 * n * t / m))) - g * t*m/n;//координата у
118
119 //вывод в файл и на экран
120 cout<<"X="<<x<<endl;
121 cout<<"Y="<<y<<endl<<endl;
122 fout<<x<<" ";
123 fout<<y<<" ";
124 fout<<endl;
125
126 }
127 }
128
129 int main(void)
130 {
131
132 double v0,r,m,n;//v0-начальная скорость,r-угол в градусах,m-масса;n-коэф.сопротивления ветра
133
134 cout<<"Enter start speed:"<<endl;
135 cin>>v0;
136 cout<<"Enter angle less than 90 deg:"<<endl;
137 cin>>r;
138 cout<<"Enter mass:"<<endl;
139 cin>>m;
140 cout<<"Coefficient of resistance:"<<endl;
141 cin>>n;
142
143 Func(v0,r);
144 Verle1(n,m,v0,r);
145 Verle2(n,m,v0,r);
146 method4(n,m,v0,r);
147
148
149
150 int k;
151 cin>>k;
152 return 0;
153 }
