python　ルンゲクッタ法

python
1import numpy as np
2import matplotlib.pyplot as plt
3from numpy.lib.type_check import isrealobj
4
5def f1(φ2):
6    return φ2
7
8def f2(i,φ1,φ2,Bc):
9    return (i-np.sin(φ1)-φ2)/Bc
10
11a = 0.0
12b = 2.0
13N = 200
14h = (b-a)/N
15j = 0
16Bc = 4.0
17
18t = np.arange(a,b,h)
19φ1_0 = 0.0
20φ2_0 = 0.0
21i_0 = 0.5
22
23φ1 = np.empty(N)
24φ2 = np.empty(N)
25i = np.empty(N)
26φ1[0] = φ1_0
27φ2[0] = φ2_0
28i[0] = i_0
29
30
31while j < N-1:
32    
33    print(f'{j=}\t{φ2[j]=:.6f}\t{t[j]=:.3f}')
34
35    k_11 = h * f1(φ2[j])
36    k_12 = h * f2(φ1[j] , φ2[j] , i[j] , Bc)
37   
38    k_21 = h * f1(φ2[j] + (h/2) * k_12)
39    k_22 = h * f2(φ1[j] + (h/2) * k_11,  φ2[j] + (h/2) * k_12 , i[j] , Bc)
40   
41    k_31 = h * f1(φ2[j] + (h/2) * k_22)
42    k_32 = h * f2(φ1[j] + (h/2) * k_21 , φ2[j] + (h/2) * k_22 , i[j] , Bc)
43
44    k_41 = h * f1(φ2[j] + h * k_32)
45    k_42 = h * f2(φ1[j] + h * k_31 , φ2[j] + h * k_32 , i[j] , Bc)     
46    
47    φ1[j+1] = φ1[j] + 1/6 * (k_11 + 2 * k_21 + 2 * k_31 + k_41) 
48    φ2[j+1] = φ2[j] + 1/6 * (k_12 + 2 * k_22 + 2 * k_32 + k_42)   
49
50    if j == N-1:
51      break
52
53    j += 1
54
55plt.plot (t,φ2)
56plt.xlabel("t")
57plt.ylabel("φ2")
58plt.show()
59```ルンゲクッタ法を用いたコードです。
60Φ2とｔについてのグラフをプロットしたいのですが、結果は交流電流のようなグラフになるはずのところがぜんぜん異なる値とグラフが出てきてしまいます。どこを修正すればよいでしょうか。
61![イメージ説明](e0048c9bcbab3030a4a545fce94e12aa.png)

行動規範の内容に同意します

回答1件

ベストアンサー

i[j+1] が設定されていません。例えば以下の様に +1 すると、-x^2 状のカーブになります。

python
1    φ1[j+1] = φ1[j] + 1/6 * (k_11 + 2 * k_21 + 2 * k_31 + k_41) 
2    φ2[j+1] = φ2[j] + 1/6 * (k_12 + 2 * k_22 + 2 * k_32 + k_42)
3    i[j+1] = i[j] + 1  ## 追加