質問編集履歴
3
プログラムの大まかな流れを追加しました
test
CHANGED
|
@@ -1 +1 @@
|
|
|
1
|
-
オイラー法に関して、厳密解との誤差グラフでう
|
|
1
|
+
オイラー法に関して、厳密解との誤差グラフで予想と違う挙動を示します
|
test
CHANGED
|
@@ -1,10 +1,10 @@
|
|
|
1
1
|
###目的と予想していた実行結果
|
|
2
2
|
|
|
3
|
-
定数係数の二階微斉次・非斉次分方程式の数値解をオイラー法で求め、厳密解と
|
|
3
|
+
定数係数の二階微斉次・非斉次分方程式の数値解をオイラー法で求め、厳密解と比較を使用としています
|
|
4
4
|
|
|
5
5
|
比較方法として、数値解の差を絶対値を取り、プロットしました。
|
|
6
6
|
|
|
7
|
-
オイラー法は積み残し誤差が発生する手法なので、プロットされたグラフは発散することが予想されます。
|
|
7
|
+
オイラー法は積み残し誤差が発生する手法なので、プロットされたグラフは徐々に発散することが予想されます。
|
|
8
8
|
|
|
9
9
|
今回は問題を簡単にするため、以下のような単振動を想定した定数に設定しています。
|
|
10
10
|
|
|
@@ -16,9 +16,53 @@
|
|
|
16
16
|
|
|
17
17
|
###問題点
|
|
18
18
|
|
|
19
|
-
プロットされたグラフ
|
|
19
|
+
プロットされたグラフは徐々に発散するようなものを想定していたのですが、途中まで収束する方向に進み、ある地点から発散へと転換していきました。
|
|
20
|
-
|
|
20
|
+
|
|
21
|
-
実
|
|
21
|
+
実際のプロットされたグラフは下の方に添付しています。
|
|
22
|
+
|
|
23
|
+
また、質問のついでですが、プログラム初心者のたclassやファイル間でのやり取りを練習するためにこのような複雑なコードになってしまっています。こうすればもっと良いプログラムになるなどのアドバイスがあれば合わせてお願いいたします。
|
|
24
|
+
|
|
25
|
+
|
|
26
|
+
|
|
27
|
+
###プログラムの大まかな流れ
|
|
28
|
+
|
|
29
|
+
最終目的は先述した通りです。各行で行っている事や変数の意味はコメントアウトしているのでそちらをご覧ください。
|
|
30
|
+
|
|
31
|
+
「class numerical」について(二階微分方程式を数値的に解くことを目的としたクラスです)
|
|
32
|
+
|
|
33
|
+
・「def __init__」:扱う微分方程式の係数などの条件を設定しています。
|
|
34
|
+
|
|
35
|
+
・「def derivative」:ある時刻における二階の微分係数をで求めています。引数はx,z,tで、返り値は一つの定数です
|
|
36
|
+
|
|
37
|
+
・「def euler」:オイラー法で微分方程式の数値解を求めています。引数は無く、返り値は従属変数の数値解のみが配列で返っています。
|
|
38
|
+
|
|
39
|
+
・「def exact」:引数で与えられた関数(微分方程式の厳密解を想定)をオイラー法で解いた時と同じ条件(=独立変数のレンジと刻み幅)で解いています。引数は関数となっており、その関数の条件は一つの変数、返り値は定数としています。
|
|
40
|
+
|
|
41
|
+
|
|
42
|
+
|
|
43
|
+
「class exact」について(二階微分方程式の厳密解を初期条件に合わせた特解を求めることを目的としています、ただし、非斉次方程式の場合の計算についてはもう少し汎用性が出るように現在制作中です)
|
|
44
|
+
|
|
45
|
+
・「def __init__」：扱う微分方程式の係数等の条件を設定しています
|
|
46
|
+
|
|
47
|
+
・「def oed_exact」:初期条件に合わせた特解を求め、その関数自体を返り値に設定しています。また、特解の決定には、特性方程式が重解の時とそうでないときで場合分けをしており、各係数を行列を用いて計算しています。
|
|
48
|
+
|
|
49
|
+
|
|
50
|
+
|
|
51
|
+
「class visualize」について(与えられた二つの配列を視覚化することを目的としたクラスです)
|
|
52
|
+
|
|
53
|
+
・「def __init__」:扱う配列とデータの名前を取得しています。
|
|
54
|
+
|
|
55
|
+
・「def plot」:シンプルな二次元のグラフをプロットしています。
|
|
56
|
+
|
|
57
|
+
・「stakplot」:新たにもう一つ配列(y成分のデータ)を取得し、事前に取得したデータと合わせてプロットしています。
|
|
58
|
+
|
|
59
|
+
・「def complot」:新たにもう一つ配列(y成分のデータ)を取得し、事前に取得したデータとの差の値に絶対値を取ってプロットしています。
|
|
60
|
+
|
|
61
|
+
・「def animation」:事前に取得したy軸のデータを利用してアニメーション(詳細にはパラパラ漫画的なもの)を作ってい表示しています。
|
|
62
|
+
|
|
63
|
+
|
|
64
|
+
|
|
65
|
+
以上です。実際のコードは次のようになっています。
|
|
22
66
|
|
|
23
67
|
|
|
24
68
|
|
|
@@ -100,7 +144,7 @@
|
|
|
100
144
|
|
|
101
145
|
self.derivative=derivative
|
|
102
146
|
|
|
103
|
-
def ode
|
|
147
|
+
def ode_exact(self):#非斉次の二階微分方程式
|
|
104
148
|
|
|
105
149
|
m,c,k=self.constant[0],self.constant[1],self.constant[2]
|
|
106
150
|
|
|
@@ -274,7 +318,7 @@
|
|
|
274
318
|
|
|
275
319
|
f=odeclass.exact(constant,initial,zero,zero)
|
|
276
320
|
|
|
277
|
-
fanc=f.ode
|
|
321
|
+
fanc=f.ode_exact()#厳密解を初期条件に合わせた特解を決定
|
|
278
322
|
|
|
279
323
|
step=0.01
|
|
280
324
|
|
|
@@ -298,7 +342,7 @@
|
|
|
298
342
|
|
|
299
343
|
###実行結果
|
|
300
344
|
|
|
301
|
-
オイラー法で解いた数値解と解析解との差
|
|
345
|
+
オイラー法で解いた数値解と解析解との差の値に対して絶対値を取り、時間との関係を表示しています。
|
|
302
346
|
|
|
303
347
|
|
|
304
348
|
|
2
説明を補足しました
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -1,14 +1,4 @@
|
|
|
1
|
-
###二階微分方程式をオイラー法で解き、厳密解との比較をしようとしています。
|
|
2
|
-
|
|
3
|
-
|
|
4
|
-
|
|
5
|
-
ここに質問の内容を詳しく書いてください。
|
|
6
|
-
|
|
7
|
-
（例）PHP(CakePHP)で●●なシステムを作っています。
|
|
8
|
-
|
|
9
|
-
■■な機能を実装中に以下のエラーメッセージが発生しました。
|
|
10
|
-
|
|
11
|
-
###
|
|
1
|
+
###目的と予想していた実行結果
|
|
12
2
|
|
|
13
3
|
定数係数の二階微斉次・非斉次分方程式の数値解をオイラー法で求め、厳密解との比較を使用としています
|
|
14
4
|
|
|
@@ -22,10 +12,14 @@
|
|
|
22
12
|
|
|
23
13
|
t0=0,x0=2,v0=0
|
|
24
14
|
|
|
15
|
+
|
|
16
|
+
|
|
25
17
|
###問題点
|
|
26
18
|
|
|
27
19
|
プロットされたグラフが発散せず、うなりのようなものが生じました。
|
|
28
20
|
|
|
21
|
+
実行結果は下の方に添付しています
|
|
22
|
+
|
|
29
23
|
|
|
30
24
|
|
|
31
25
|
|
1
実行結果を追加しました,実行目的の補足をしました
test
CHANGED
|
File without changes
|
test
CHANGED
|
@@ -8,9 +8,25 @@
|
|
|
8
8
|
|
|
9
9
|
■■な機能を実装中に以下のエラーメッセージが発生しました。
|
|
10
10
|
|
|
11
|
-
|
|
11
|
+
###最終目的と予想していた実行結果
|
|
12
|
+
|
|
12
|
-
|
|
13
|
+
定数係数の二階微斉次・非斉次分方程式の数値解をオイラー法で求め、厳密解との比較を使用としています
|
|
14
|
+
|
|
15
|
+
比較方法として、数値解の差を絶対値を取り、プロットしました。
|
|
16
|
+
|
|
17
|
+
オイラー法は積み残し誤差が発生する手法なので、プロットされたグラフは発散することが予想されます。
|
|
18
|
+
|
|
19
|
+
今回は問題を簡単にするため、以下のような単振動を想定した定数に設定しています。
|
|
20
|
+
|
|
21
|
+
1*d^2x/dt^2 =-1*x
|
|
22
|
+
|
|
23
|
+
t0=0,x0=2,v0=0
|
|
24
|
+
|
|
25
|
+
###問題点
|
|
26
|
+
|
|
13
|
-
|
|
27
|
+
プロットされたグラフが発散せず、うなりのようなものが生じました。
|
|
28
|
+
|
|
29
|
+
|
|
14
30
|
|
|
15
31
|
|
|
16
32
|
|
|
@@ -258,37 +274,17 @@
|
|
|
258
274
|
|
|
259
275
|
initial=[0,2.0,0]#初期条件:t0,x0,v0
|
|
260
276
|
|
|
261
|
-
f0=2.0#調和振動加振力の最大値
|
|
262
|
-
|
|
263
|
-
w=2#角振動数
|
|
264
|
-
|
|
265
|
-
def special(t):#非斉次方程式の特解
|
|
266
|
-
|
|
267
|
-
alpha=f0*(k-m*w*w+1j*c*w)/(pow(k-m*w*w,2)+pow(c*w,2))
|
|
268
|
-
|
|
269
|
-
return alpha*pow(math.e,-1j*w*t)
|
|
270
|
-
|
|
271
|
-
def derivative(t):#上記の特解の微分係数
|
|
272
|
-
|
|
273
|
-
alpha=f0*(k-m*w*w+1j*c*w)/(pow(k-m*w*w,2)+pow(c*w,2))
|
|
274
|
-
|
|
275
|
-
return -1j*w*alpha*pow(math.e,-1j*w*t)
|
|
276
|
-
|
|
277
277
|
def zero(t):
|
|
278
278
|
|
|
279
279
|
return 0
|
|
280
280
|
|
|
281
|
-
def inhom(t):#非斉次項
|
|
282
|
-
|
|
283
|
-
return f0*pow(math.e,-1j*w*t).real
|
|
284
|
-
|
|
285
281
|
f=odeclass.exact(constant,initial,zero,zero)
|
|
286
282
|
|
|
287
283
|
fanc=f.ode2_exact()#厳密解を初期条件に合わせた特解を決定
|
|
288
284
|
|
|
289
285
|
step=0.01
|
|
290
286
|
|
|
291
|
-
z=odeclass.numerical(constant,zero,initial,step,
|
|
287
|
+
z=odeclass.numerical(constant,zero,initial,step,100)
|
|
292
288
|
|
|
293
289
|
t=z.variable
|
|
294
290
|
|
|
@@ -298,9 +294,7 @@
|
|
|
298
294
|
|
|
299
295
|
v=odeclass.visualize(t,x1,"euler")
|
|
300
296
|
|
|
301
|
-
v.stakplot((0,80),(-2.5,2.5),x2,"exact2")
|
|
302
|
-
|
|
303
|
-
v.complot((0,
|
|
297
|
+
v.complot((0,100),(0,0.1),x2,"error")
|
|
304
298
|
|
|
305
299
|
|
|
306
300
|
|
|
@@ -308,6 +302,12 @@
|
|
|
308
302
|
|
|
309
303
|
|
|
310
304
|
|
|
305
|
+
###実行結果
|
|
306
|
+
|
|
307
|
+
オイラー法で解いた数値解と解析解との差を絶対値を取って表示しています。
|
|
308
|
+
|
|
309
|
+
|
|
310
|
+
|
|
311
311
|
### 試したこと
|
|
312
312
|
|
|
313
313
|
|