質問編集履歴
4
修正
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -298,7 +298,7 @@
|
|
|
298
298
|
|
|
299
299
|
```
|
|
300
300
|
|
|
301
|
-
リストxiの中には要素が１０個あるのでjを0から10まで増やしながら対応する_v[-1]とsum(_v)を返してもらうコードを書きました。
|
|
301
|
+
追記(訂正） リストxiの中には要素が１０個あるのでjを0から10まで増やしながら対応する_v[-1]とsum(_v)を返してもらうコードを書きました。これをどうにかしたいと思っています。
|
|
302
302
|
|
|
303
303
|
```
|
|
304
304
|
|
|
@@ -322,7 +322,7 @@
|
|
|
322
322
|
|
|
323
323
|
|
|
324
324
|
|
|
325
|
-
while j <= 9 :
|
|
325
|
+
while j <= 9 :
|
|
326
326
|
|
|
327
327
|
t = 0
|
|
328
328
|
|
|
@@ -336,7 +336,7 @@
|
|
|
336
336
|
|
|
337
337
|
|
|
338
338
|
|
|
339
|
-
while x <= 0.60: #境界面がL以下の位置にあるとき
|
|
339
|
+
while x <= 0.60: #境界面がL以下の位置にあるとき これを超えると水が完全に切れたことに
|
|
340
340
|
|
|
341
341
|
|
|
342
342
|
|
|
@@ -412,7 +412,7 @@
|
|
|
412
412
|
|
|
413
413
|
|
|
414
414
|
|
|
415
|
-
Hi = sum(_vi)*dt + (_vi[-1])**2/2/g
|
|
415
|
+
Hi = sum(_vi)*dt + (_vi[-1])**2/2/g #最高到達点を求める式
|
|
416
416
|
|
|
417
417
|
|
|
418
418
|
|
|
@@ -486,7 +486,17 @@
|
|
|
486
486
|
|
|
487
487
|
問題はペットボトルロケットの水の量に対する最高到達点を求めようというものです。
|
|
488
488
|
|
|
489
|
+
ペットボトルロケットを真上に発射すると水が入っている状態ではルンゲクッタで速度を求めることが出来ますが、水がなくなるとその時の速度を初速度とする投げ上げ運動に変わります。
|
|
490
|
+
|
|
491
|
+
最高到達点は、水が切れた時の速度と到達点が分かれば求めることが出来ます。
|
|
492
|
+
|
|
493
|
+
|
|
494
|
+
|
|
495
|
+
つまり、最初の座標ｘを変えながら速度と到達点を返してくれるコードを書くことを目的としています。
|
|
496
|
+
|
|
497
|
+
|
|
498
|
+
|
|
489
|
-
この問題に至るまでに時間と速度の関係等を求める問題があったのですが（水の量は固定されていた）、この問題は時間は関係なく今度は初期値を変える必要があり頭が混乱しています。
|
|
499
|
+
この問題に至るまでに時間と速度の関係等を求める問題があったのですが（水の量は固定されていた）、この問題は時間は関係なく今度は初期値（最初にどれくらい水が入っているのか）を変える必要があり頭が混乱しています。
|
|
490
500
|
|
|
491
501
|
力技で自分で水の量を毎回変えながら最高到達点を求めてグラフを描くことには成功したのですが...
|
|
492
502
|
|
3
コード中に出てくる文字の意味を追記しました。
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -230,13 +230,77 @@
|
|
|
230
230
|
|
|
231
231
|
|
|
232
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
|
+
import math
|
|
260
|
+
|
|
261
|
+
import numpy as np
|
|
262
|
+
|
|
263
|
+
import matplotlib.pyplot as plt
|
|
264
|
+
|
|
265
|
+
```
|
|
266
|
+
|
|
267
|
+
文字の意味です。
|
|
268
|
+
|
|
269
|
+
```
|
|
270
|
+
|
|
271
|
+
pi = math.pi
|
|
272
|
+
|
|
273
|
+
S1 = pi*(0.08/2)**2 #ペットボトルの断面積 m^2
|
|
274
|
+
|
|
275
|
+
S2 = pi*(0.02/2)**2 #噴出口の断面積 m^2
|
|
276
|
+
|
|
277
|
+
SS = (0.02/ 0.08)**2 # 断面積の比 S1/S2
|
|
278
|
+
|
|
233
|
-
|
|
279
|
+
L = 0.60 # ペットボトルの長さ m
|
|
280
|
+
|
|
234
|
-
|
|
281
|
+
rho = 1000.0 # 水の密度 kg/m^3
|
|
282
|
+
|
|
235
|
-
|
|
283
|
+
atm = 101325 # atm → N/m^3
|
|
284
|
+
|
|
236
|
-
|
|
285
|
+
P0 = 4.0*atm # ペットボトル内部の圧力 N/m^2
|
|
286
|
+
|
|
287
|
+
PA = 1.0*atm # 大気圧 N/m^2
|
|
288
|
+
|
|
289
|
+
X0 = 0.48 # 水と空気の境界面 水の量 m
|
|
290
|
+
|
|
291
|
+
gm = 7/5 # Cp/Cv 気体の比熱比
|
|
292
|
+
|
|
293
|
+
dt = 0.000001 #時間幅 sec
|
|
294
|
+
|
|
295
|
+
g = 9.8 #重力加速度 m/sec^2
|
|
296
|
+
|
|
297
|
+
M0 = 0.12 # ロケットの質量（水抜き） kg
|
|
298
|
+
|
|
299
|
+
```
|
|
300
|
+
|
|
237
|
-
|
|
301
|
+
リストxiの中には要素が１０個あるのでjを0から10まで増やしながら対応する_v[-1]とsum(_v)を返してもらうコードを書きました。
|
|
302
|
+
|
|
238
|
-
|
|
303
|
+
```
|
|
239
|
-
|
|
240
304
|
|
|
241
305
|
j = 0
|
|
242
306
|
|
2
書式の改善
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -52,9 +52,7 @@
|
|
|
52
52
|
|
|
53
53
|
while i <= 0.55 :
|
|
54
54
|
|
|
55
|
-
|
|
56
|
-
|
|
57
|
-
|
|
55
|
+
water = S1 * rho * (0.60 - i )
|
|
58
56
|
|
|
59
57
|
|
|
60
58
|
|
|
@@ -124,7 +122,7 @@
|
|
|
124
122
|
|
|
125
123
|
```
|
|
126
124
|
|
|
127
|
-
|
|
125
|
+
|
|
128
126
|
|
|
129
127
|
_water = []
|
|
130
128
|
|
|
@@ -230,7 +228,7 @@
|
|
|
230
228
|
|
|
231
229
|
i += di;
|
|
232
230
|
|
|
233
|
-
|
|
231
|
+
|
|
234
232
|
|
|
235
233
|
while以降右にずれていないのですが実際上はwhileの下はずれています。
|
|
236
234
|
|
|
@@ -238,127 +236,151 @@
|
|
|
238
236
|
|
|
239
237
|
追記 リストxiには１０個の要素があるのでその要素を1つずつ取り出すべくjを導入して各jに対して_v[-1]とmax(_v)を配列に追加するコードを作成してみましたがずっと実行したまま進みません...
|
|
240
238
|
|
|
239
|
+
|
|
240
|
+
|
|
241
|
+
j = 0
|
|
242
|
+
|
|
243
|
+
dj = 1
|
|
244
|
+
|
|
245
|
+
|
|
246
|
+
|
|
247
|
+
_t=[]
|
|
248
|
+
|
|
249
|
+
_x=[]
|
|
250
|
+
|
|
251
|
+
_vf=[]
|
|
252
|
+
|
|
253
|
+
_vsum = []
|
|
254
|
+
|
|
255
|
+
_u=[]
|
|
256
|
+
|
|
257
|
+
|
|
258
|
+
|
|
259
|
+
|
|
260
|
+
|
|
261
|
+
while j <= 9 : #jを定義したい
|
|
262
|
+
|
|
263
|
+
t = 0
|
|
264
|
+
|
|
265
|
+
v = 0
|
|
266
|
+
|
|
267
|
+
x = _xi[j]#xの初期値 これを変えていきたい。
|
|
268
|
+
|
|
269
|
+
dt = 0.000001
|
|
270
|
+
|
|
271
|
+
|
|
272
|
+
|
|
273
|
+
|
|
274
|
+
|
|
275
|
+
while x <= 0.60: #境界面がL以下の位置にあるとき
|
|
276
|
+
|
|
277
|
+
|
|
278
|
+
|
|
279
|
+
#ルンゲクッタで時間tでの水の座標を求める。
|
|
280
|
+
|
|
281
|
+
k1 = dxdt(x)*dt;
|
|
282
|
+
|
|
283
|
+
k2 = dxdt((x+k1/2.0))*dt;
|
|
284
|
+
|
|
285
|
+
k3 = dxdt((x+k2/2.0))*dt;
|
|
286
|
+
|
|
287
|
+
k4 = dxdt((x+k3))*dt;
|
|
288
|
+
|
|
289
|
+
k = (k1 + 2.0*(k2 + k3 )+ k4)/6.0;
|
|
290
|
+
|
|
291
|
+
|
|
292
|
+
|
|
293
|
+
u = math.sqrt((2/rho)*(P0*((i/x)**gm) - PA)); # m/s
|
|
294
|
+
|
|
295
|
+
|
|
296
|
+
|
|
297
|
+
#ルンゲクッタで速度を求める。
|
|
298
|
+
|
|
299
|
+
k1v = dvdt(x)*dt;
|
|
300
|
+
|
|
301
|
+
k2v = dvdt((x+k1v/2.0))*dt;
|
|
302
|
+
|
|
303
|
+
k3v = dvdt((x +k2v/2.0))*dt;
|
|
304
|
+
|
|
305
|
+
k4v = dvdt((x +k3v))*dt;
|
|
306
|
+
|
|
307
|
+
kv = (k1v + 2.0*(k2v + k3v) *k4v)/6.0;
|
|
308
|
+
|
|
309
|
+
|
|
310
|
+
|
|
311
|
+
vf = _v[-1] #最終速度
|
|
312
|
+
|
|
313
|
+
vsum = sum(_v) #vの和
|
|
314
|
+
|
|
315
|
+
|
|
316
|
+
|
|
317
|
+
|
|
318
|
+
|
|
319
|
+
#配列へ
|
|
320
|
+
|
|
321
|
+
_t.append(t)
|
|
322
|
+
|
|
323
|
+
_x.append(x)
|
|
324
|
+
|
|
325
|
+
_u.append(u)
|
|
326
|
+
|
|
327
|
+
_v.append(v)
|
|
328
|
+
|
|
329
|
+
_vf.append(vf)
|
|
330
|
+
|
|
331
|
+
_vsum.append(vsum)
|
|
332
|
+
|
|
333
|
+
|
|
334
|
+
|
|
335
|
+
|
|
336
|
+
|
|
337
|
+
|
|
338
|
+
|
|
339
|
+
t += dt;
|
|
340
|
+
|
|
341
|
+
x += k;
|
|
342
|
+
|
|
343
|
+
v += kv;
|
|
344
|
+
|
|
345
|
+
|
|
346
|
+
|
|
347
|
+
|
|
348
|
+
|
|
349
|
+
|
|
350
|
+
|
|
351
|
+
Hi = sum(_vi)*dt + (_vi[-1])**2/2/gm #最高到達点を求める式
|
|
352
|
+
|
|
353
|
+
|
|
354
|
+
|
|
355
|
+
_Hi.append(Hi)
|
|
356
|
+
|
|
357
|
+
|
|
358
|
+
|
|
359
|
+
j += dj;
|
|
360
|
+
|
|
361
|
+
```
|
|
362
|
+
|
|
363
|
+
こういうエラーが出てきてしまいます。
|
|
364
|
+
|
|
241
365
|
```ここに言語を入力
|
|
242
366
|
|
|
243
|
-
j = 0
|
|
244
|
-
|
|
245
|
-
dj = 1
|
|
246
|
-
|
|
247
|
-
|
|
248
|
-
|
|
249
|
-
_t=[]
|
|
250
|
-
|
|
251
|
-
_x=[]
|
|
252
|
-
|
|
253
|
-
_vf=[]
|
|
254
|
-
|
|
255
|
-
_vsum = []
|
|
256
|
-
|
|
257
|
-
_u=[]
|
|
258
|
-
|
|
259
|
-
|
|
260
|
-
|
|
261
|
-
|
|
262
|
-
|
|
263
|
-
while j <= 9 : #jを定義したい
|
|
264
|
-
|
|
265
|
-
t = 0
|
|
266
|
-
|
|
267
|
-
v = 0
|
|
268
|
-
|
|
269
|
-
x = _xi[j]#xの初期値 これを変えていきたい。
|
|
270
|
-
|
|
271
|
-
dt = 0.000001
|
|
272
|
-
|
|
273
|
-
|
|
274
|
-
|
|
275
|
-
|
|
276
|
-
|
|
277
|
-
while x <= 0.60: #境界面がL以下の位置にあるとき
|
|
278
|
-
|
|
279
|
-
|
|
280
|
-
|
|
281
|
-
#ルンゲクッタで時間tでの水の座標を求める。
|
|
282
|
-
|
|
283
|
-
k1 = dxdt(x)*dt;
|
|
284
|
-
|
|
285
|
-
k2 = dxdt((x+k1/2.0))*dt;
|
|
286
|
-
|
|
287
|
-
k3 = dxdt((x+k2/2.0))*dt;
|
|
288
|
-
|
|
289
|
-
k4 = dxdt((x+k3))*dt;
|
|
290
|
-
|
|
291
|
-
|
|
367
|
+
IndexErrorTraceback (most recent call last)
|
|
292
|
-
|
|
293
|
-
|
|
294
|
-
|
|
368
|
+
|
|
295
|
-
|
|
369
|
+
<ipython-input-17-14b49cbf9f76> in <module>()
|
|
296
|
-
|
|
297
|
-
|
|
298
|
-
|
|
299
|
-
|
|
370
|
+
|
|
300
|
-
|
|
301
|
-
k1v = dvdti(x)*dt;
|
|
302
|
-
|
|
303
|
-
k2v = dvdti((x+k1v/2.0))*dt;
|
|
304
|
-
|
|
305
|
-
k3v = dvdti((x +k2v/2.0))*dt;
|
|
306
|
-
|
|
307
|
-
k4v = dvdti((x +k3v))*dt;
|
|
308
|
-
|
|
309
|
-
kv = (k1v + 2.0*(k2v + k3v) *k4v)/6.0;
|
|
371
|
+
35 kv = (k1v + 2.0*(k2v + k3v) *k4v)/6.0;
|
|
372
|
+
|
|
310
|
-
|
|
373
|
+
36
|
|
311
|
-
|
|
312
|
-
|
|
374
|
+
|
|
313
|
-
vf = _v[-1]
|
|
375
|
+
---> 37 vf = _v[-1] #最終速度
|
|
314
|
-
|
|
376
|
+
|
|
315
|
-
vsum = sum(_v)
|
|
377
|
+
38 vsum = sum(_v) #vの和
|
|
316
|
-
|
|
317
|
-
|
|
318
|
-
|
|
319
|
-
|
|
320
|
-
|
|
378
|
+
|
|
321
|
-
|
|
379
|
+
39
|
|
322
|
-
|
|
323
|
-
|
|
380
|
+
|
|
324
|
-
|
|
325
|
-
|
|
381
|
+
|
|
326
|
-
|
|
327
|
-
|
|
382
|
+
|
|
328
|
-
|
|
329
|
-
_v.append(v)
|
|
330
|
-
|
|
331
|
-
_vf.append(vf)
|
|
332
|
-
|
|
333
|
-
_vsum.append(vsum)
|
|
334
|
-
|
|
335
|
-
|
|
336
|
-
|
|
337
|
-
|
|
338
|
-
|
|
339
|
-
|
|
340
|
-
|
|
341
|
-
t += dt;
|
|
342
|
-
|
|
343
|
-
x += k;
|
|
344
|
-
|
|
345
|
-
v += kv;
|
|
346
|
-
|
|
347
|
-
|
|
348
|
-
|
|
349
|
-
|
|
350
|
-
|
|
351
|
-
|
|
352
|
-
|
|
353
|
-
|
|
383
|
+
IndexError: list index out of range
|
|
354
|
-
|
|
355
|
-
|
|
356
|
-
|
|
357
|
-
_Hi.append(Hi)
|
|
358
|
-
|
|
359
|
-
|
|
360
|
-
|
|
361
|
-
j += dj;
|
|
362
384
|
|
|
363
385
|
```
|
|
364
386
|
|
1
viはvに変更しましたがコードをコピーしてここに張り付けるとwhile文があってもずれてくれないのでご了承ください。新しいコードを作成しました
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -110,7 +110,7 @@
|
|
|
110
110
|
|
|
111
111
|
|
|
112
112
|
|
|
113
|
-
def dvdt
|
|
113
|
+
def dvdt(x):
|
|
114
114
|
|
|
115
115
|
u_square = (2.0/rho)*(P0*((X0/x)**gm) - PA)
|
|
116
116
|
|
|
@@ -128,17 +128,15 @@
|
|
|
128
128
|
|
|
129
129
|
_water = []
|
|
130
130
|
|
|
131
|
-
#_Hi = []
|
|
132
|
-
|
|
133
131
|
|
|
134
132
|
|
|
135
133
|
while i <= 0.55:
|
|
136
134
|
|
|
135
|
+
|
|
136
|
+
|
|
137
137
|
t = 0
|
|
138
138
|
|
|
139
|
-
|
|
139
|
+
v = 0
|
|
140
|
-
|
|
141
|
-
vi = 0
|
|
142
140
|
|
|
143
141
|
x = i #xの初期値
|
|
144
142
|
|
|
@@ -158,7 +156,7 @@
|
|
|
158
156
|
|
|
159
157
|
_x=[]
|
|
160
158
|
|
|
161
|
-
_v
|
|
159
|
+
_v=[]
|
|
162
160
|
|
|
163
161
|
_u=[]
|
|
164
162
|
|
|
@@ -212,7 +210,7 @@
|
|
|
212
210
|
|
|
213
211
|
_u.append(u)
|
|
214
212
|
|
|
215
|
-
_v
|
|
213
|
+
_v.append(v)
|
|
216
214
|
|
|
217
215
|
|
|
218
216
|
|
|
@@ -224,15 +222,9 @@
|
|
|
224
222
|
|
|
225
223
|
x += k;
|
|
226
224
|
|
|
227
|
-
v
|
|
225
|
+
v += kv;
|
|
228
|
-
|
|
229
|
-
|
|
230
|
-
|
|
231
|
-
|
|
226
|
+
|
|
232
|
-
|
|
233
|
-
|
|
234
|
-
|
|
235
|
-
|
|
227
|
+
|
|
236
228
|
|
|
237
229
|
|
|
238
230
|
|
|
@@ -244,7 +236,137 @@
|
|
|
244
236
|
|
|
245
237
|
|
|
246
238
|
|
|
247
|
-
|
|
239
|
+
追記 リストxiには１０個の要素があるのでその要素を1つずつ取り出すべくjを導入して各jに対して_v[-1]とmax(_v)を配列に追加するコードを作成してみましたがずっと実行したまま進みません...
|
|
240
|
+
|
|
241
|
+
```ここに言語を入力
|
|
242
|
+
|
|
243
|
+
j = 0
|
|
244
|
+
|
|
245
|
+
dj = 1
|
|
246
|
+
|
|
247
|
+
|
|
248
|
+
|
|
249
|
+
_t=[]
|
|
250
|
+
|
|
251
|
+
_x=[]
|
|
252
|
+
|
|
253
|
+
_vf=[]
|
|
254
|
+
|
|
255
|
+
_vsum = []
|
|
256
|
+
|
|
257
|
+
_u=[]
|
|
258
|
+
|
|
259
|
+
|
|
260
|
+
|
|
261
|
+
|
|
262
|
+
|
|
263
|
+
while j <= 9 : #jを定義したい
|
|
264
|
+
|
|
265
|
+
t = 0
|
|
266
|
+
|
|
267
|
+
v = 0
|
|
268
|
+
|
|
269
|
+
x = _xi[j]#xの初期値 これを変えていきたい。
|
|
270
|
+
|
|
271
|
+
dt = 0.000001
|
|
272
|
+
|
|
273
|
+
|
|
274
|
+
|
|
275
|
+
|
|
276
|
+
|
|
277
|
+
while x <= 0.60: #境界面がL以下の位置にあるとき
|
|
278
|
+
|
|
279
|
+
|
|
280
|
+
|
|
281
|
+
#ルンゲクッタで時間tでの水の座標を求める。
|
|
282
|
+
|
|
283
|
+
k1 = dxdt(x)*dt;
|
|
284
|
+
|
|
285
|
+
k2 = dxdt((x+k1/2.0))*dt;
|
|
286
|
+
|
|
287
|
+
k3 = dxdt((x+k2/2.0))*dt;
|
|
288
|
+
|
|
289
|
+
k4 = dxdt((x+k3))*dt;
|
|
290
|
+
|
|
291
|
+
k = (k1 + 2.0*(k2 + k3 )+ k4)/6.0;
|
|
292
|
+
|
|
293
|
+
|
|
294
|
+
|
|
295
|
+
u = math.sqrt((2/rho)*(P0*((i/x)**gm) - PA)); # m/s
|
|
296
|
+
|
|
297
|
+
|
|
298
|
+
|
|
299
|
+
#ルンゲクッタで速度を求める。
|
|
300
|
+
|
|
301
|
+
k1v = dvdti(x)*dt;
|
|
302
|
+
|
|
303
|
+
k2v = dvdti((x+k1v/2.0))*dt;
|
|
304
|
+
|
|
305
|
+
k3v = dvdti((x +k2v/2.0))*dt;
|
|
306
|
+
|
|
307
|
+
k4v = dvdti((x +k3v))*dt;
|
|
308
|
+
|
|
309
|
+
kv = (k1v + 2.0*(k2v + k3v) *k4v)/6.0;
|
|
310
|
+
|
|
311
|
+
|
|
312
|
+
|
|
313
|
+
vf = _v[-1] #最終速度
|
|
314
|
+
|
|
315
|
+
vsum = sum(_v) #vの和
|
|
316
|
+
|
|
317
|
+
|
|
318
|
+
|
|
319
|
+
|
|
320
|
+
|
|
321
|
+
#配列へ
|
|
322
|
+
|
|
323
|
+
_t.append(t)
|
|
324
|
+
|
|
325
|
+
_x.append(x)
|
|
326
|
+
|
|
327
|
+
_u.append(u)
|
|
328
|
+
|
|
329
|
+
_v.append(v)
|
|
330
|
+
|
|
331
|
+
_vf.append(vf)
|
|
332
|
+
|
|
333
|
+
_vsum.append(vsum)
|
|
334
|
+
|
|
335
|
+
|
|
336
|
+
|
|
337
|
+
|
|
338
|
+
|
|
339
|
+
|
|
340
|
+
|
|
341
|
+
t += dt;
|
|
342
|
+
|
|
343
|
+
x += k;
|
|
344
|
+
|
|
345
|
+
v += kv;
|
|
346
|
+
|
|
347
|
+
|
|
348
|
+
|
|
349
|
+
|
|
350
|
+
|
|
351
|
+
|
|
352
|
+
|
|
353
|
+
Hi = sum(_vi)*dt + (_vi[-1])**2/2/gm #最高到達点を求める式
|
|
354
|
+
|
|
355
|
+
|
|
356
|
+
|
|
357
|
+
_Hi.append(Hi)
|
|
358
|
+
|
|
359
|
+
|
|
360
|
+
|
|
361
|
+
j += dj;
|
|
362
|
+
|
|
363
|
+
```
|
|
364
|
+
|
|
365
|
+
|
|
366
|
+
|
|
367
|
+
|
|
368
|
+
|
|
369
|
+
|
|
248
370
|
|
|
249
371
|
|
|
250
372
|
|