質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
87.20%
Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

解決済

python ルンゲクッタ法

Ackngawe.-
Ackngawe.-

総合スコア27

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

1回答

0評価

0クリップ

324閲覧

投稿2021/12/22 02:30

編集2021/12/22 03:22

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になってしまいます。どうすればよいでしょうか。![イメージ説明](3b0313782abbebc7a36cff90e99e922b.png)

良い質問の評価を上げる

以下のような質問は評価を上げましょう

  • 質問内容が明確
  • 自分も答えを知りたい
  • 質問者以外のユーザにも役立つ

評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

気になる質問をクリップする

クリップした質問は、後からいつでもマイページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

  • プログラミングに関係のない質問
  • やってほしいことだけを記載した丸投げの質問
  • 問題・課題が含まれていない質問
  • 意図的に内容が抹消された質問
  • 過去に投稿した質問と同じ内容の質問
  • 広告と受け取られるような投稿

評価を下げると、トップページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

jbpb0

2021/12/22 02:42 編集

plt.scatter (i, v, marker='.') のすぐ上に print(i) print(v) を(インデントを「plt.scatter」に合わせて)追加して実行してみて、値が正しく計算されてるのかを確認してみたらいかがでしょうか
Ackngawe.-

2021/12/22 02:51

確認してみたところ、iは正しくループできていましたが、vが発散してnanとなってしまっていました。 iをループせずに固定して、vの点をひとつプロットしたときはできたのですが、なぜiをループさせたらできないのでしょうか、、
jbpb0

2021/12/22 03:03 編集

> vが発散してnanとなってしまっていました。 それでしたら、 > なぜか点がグラフにプロットされません。 という質問ではなく、「v」がうまく計算できません、という質問ですね 上記分かったことを、質問を編集して追記することをお勧めします
Ackngawe.-

2021/12/22 03:19

わかりました。 ありがとうございます!
melian

2021/12/22 03:41

N = 100 としていますが、これですと np.average(fai2[6000:]) で範囲外アクセスになってしまいます。
jbpb0

2021/12/22 04:32

> np.average(fai2[6000:]) で範囲外アクセス は、 > iをループせずに固定して、その値ひとつをプロットしたときはできた の場合もダメなはずですよね
Ackngawe.-

2021/12/22 04:37

N =10000にしたらできました! ありがとうございます!

まだ回答がついていません

会員登録して回答してみよう

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
87.20%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問

同じタグがついた質問を見る

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。