python
import numpy as np import matplotlib.pyplot as plt from numpy.lib.type_check import isrealobj #関数の定義 def f1(fai2): return fai2 def f2(fai1,fai2): return (i-np.sin(fai1)-fai2)/bc #初期値 a = 0.0 b = 2.0 N = 100 h = (b-a)/N j = 0 bc = 2.0 t = np.arange(a,b,h) fai1_0 = 0.0 fai2_0 = 0.0 #配列 fai1 = np.empty(N) fai2 = np.empty(N) fai1[0] = fai1_0 fai2[0] = fai2_0 #iを増加させたとき for i in np.arange(0,2.0,0.01): #ルンゲクッタ法 for j in range (N-2): k1 = h * f1(fai2[j]) d1 = h * f2(fai1[j] , fai2[j]) k2 = h * f1(fai2[j] + d1/2) d2 = h * f2(fai1[j] + k1/2, fai2[j] + d1/2) k3 = h * f1(fai2[j] + d2/2) d3 = h * f2(fai1[j] + k2/2 , fai2[j] +d2/2) k4 = h * f1(fai2[j] + d3) d4 = h * f2(fai1[j] + k3 , fai2[j] + d3) fai1[j+1] = fai1[j] + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) fai2[j+1] = fai2[j] + 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) if j == N-1: break j += 1 #足した値を平均化する v = np.average(fai2[6000:]) plt.scatter (i, v, marker='.') #iを減少させたとき for i in np.arange(2.0,0,-0.01): #ルンゲクッタ法 for j in range (N-2): k1 = h * f1(fai2[j]) d1 = h * f2(fai1[j] , fai2[j]) k2 = h * f1(fai2[j] + d1/2) d2 = h * f2(fai1[j] + k1/2, fai2[j] + d1/2) k3 = h * f1(fai2[j] + d2/2) d3 = h * f2(fai1[j] + k2/2 , fai2[j] +d2/2) k4 = h * f1(fai2[j] + d3) d4 = h * f2(fai1[j] + k3 , fai2[j] + d3) fai1[j+1] = fai1[j] + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) fai2[j+1] = fai2[j] + 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) if j == N-1: break j += 1 #足した値を平均化する v = np.average(fai2[6000:]) plt.scatter (i, v, marker='.') plt.xlabel("v") plt.ylabel("i") plt.show() ```ルンゲクッタ法を用いたコードなのですが、横軸をv、縦軸をiとして,ルンゲクッタ法を用いて計算した6000以降の値を平均した値を、iを0から2.0までの間で繰り返し、それぞれ得た値をプロットしていくプログラムなのですが、実行してもがvが発散してnanとなり、グラフにプロットされません。 iをループせずに固定して、その値ひとつをプロットしたときはできたのですが、i をループさせた途端にできなくなってしまいました。 また、iの範囲は0から2.0までと指定しているのに、結果として出てくるグラフは―0.04から0.04になってしまいます。どうすればよいでしょうか。
plt.scatter (i, v, marker='.')
のすぐ上に
print(i)
print(v)
を(インデントを「plt.scatter」に合わせて)追加して実行してみて、値が正しく計算されてるのかを確認してみたらいかがでしょうか
確認してみたところ、iは正しくループできていましたが、vが発散してnanとなってしまっていました。
iをループせずに固定して、vの点をひとつプロットしたときはできたのですが、なぜiをループさせたらできないのでしょうか、、
> vが発散してnanとなってしまっていました。
それでしたら、
> なぜか点がグラフにプロットされません。
という質問ではなく、「v」がうまく計算できません、という質問ですね
上記分かったことを、質問を編集して追記することをお勧めします
わかりました。
ありがとうございます!
N = 100 としていますが、これですと np.average(fai2[6000:]) で範囲外アクセスになってしまいます。
> np.average(fai2[6000:]) で範囲外アクセス
は、
> iをループせずに固定して、その値ひとつをプロットしたときはできた
の場合もダメなはずですよね
N =10000にしたらできました!
ありがとうございます!
まだ回答がついていません
会員登録して回答してみよう