Информатика: Движение тела в среде
Материал из Department of Theoretical and Applied Mechanics
Версия от 13:29, 18 декабря 2015; Станислав Лебедев (обсуждение | вклад)
Лебедев Станислав
Скачать можно тут.
Описание программы: программа записывает в четыре файла результаты вычисления:
- Координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- Координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости;
- Координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
- Координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
Визуализированный результат работы программы
- o1 - координаты, рассчитанные по формуле, при движении без сопротивления воздуха;
- o2 - координаты, полученные методом Верле при линейной зависимости силы сопротивлении воздуха от скорости;
- o3 - координаты, полученные из точного решения, при линейной зависимости силы сопротивлении воздуха от скорости;
- o4 - координаты, полученные методом Верле при квадратичной зависимости силы сопротивлении воздуха от скорости.
Для тела с массой 10,сопротивлением воздуха 1, угол бросания 30°, начальная скорость 30 м/с, ускорение свободного падения 9.8 м/c^2;
Примечание: графики o1 и o2 намеренно посчитаны с малой точностью, чтобы графики не сливались.
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 (VpV(VpV(r,VS(v,t)),VS(a,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 b2.WVTS(b3.positionReal(i));
184 b2.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 }
