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

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

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

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

アルゴリズム

アルゴリズムとは、定められた目的を達成するために、プログラムの理論的な動作を定義するものです。

プログラミング言語

プログラミング言語はパソコン上で実行することができるソースコードを記述する為に扱う言語の総称です。

文字コード

文字コードとは、文字や記号をコンピュータ上で使用するために用いられるバイト表現を指します。

Python

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

解決済

オイラー法に関して、厳密解との誤差グラフで予想と違う挙動を示します

fin1920
fin1920

総合スコア0

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

アルゴリズム

アルゴリズムとは、定められた目的を達成するために、プログラムの理論的な動作を定義するものです。

プログラミング言語

プログラミング言語はパソコン上で実行することができるソースコードを記述する為に扱う言語の総称です。

文字コード

文字コードとは、文字や記号をコンピュータ上で使用するために用いられるバイト表現を指します。

Python

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

1回答

0評価

3クリップ

712閲覧

投稿2021/04/29 06:48

編集2021/05/01 13:20

###目的と予想していた実行結果
定数係数の二階微斉次・非斉次分方程式の数値解をオイラー法で求め、厳密解と比較を使用としています
比較方法として、数値解の差を絶対値を取り、プロットしました。
オイラー法は積み残し誤差が発生する手法なので、プロットされたグラフは徐々に発散することが予想されます。
今回は問題を簡単にするため、以下のような単振動を想定した定数に設定しています。
1d^2x/dt^2 =-1x
t0=0,x0=2,v0=0

###問題点
プロットされたグラフは徐々に発散するようなものを想定していたのですが、途中まで収束する方向に進み、ある地点から発散へと転換していきました。
実際のプロットされたグラフは下の方に添付しています。
また、質問のついでですが、プログラム初心者のたclassやファイル間でのやり取りを練習するためにこのような複雑なコードになってしまっています。こうすればもっと良いプログラムになるなどのアドバイスがあれば合わせてお願いいたします。

###プログラムの大まかな流れ
最終目的は先述した通りです。各行で行っている事や変数の意味はコメントアウトしているのでそちらをご覧ください。
「class numerical」について(二階微分方程式を数値的に解くことを目的としたクラスです)
・「def init」:扱う微分方程式の係数などの条件を設定しています。
・「def derivative」:ある時刻における二階の微分係数をで求めています。引数はx,z,tで、返り値は一つの定数です
・「def euler」:オイラー法で微分方程式の数値解を求めています。引数は無く、返り値は従属変数の数値解のみが配列で返っています。
・「def exact」:引数で与えられた関数(微分方程式の厳密解を想定)をオイラー法で解いた時と同じ条件(=独立変数のレンジと刻み幅)で解いています。引数は関数となっており、その関数の条件は一つの変数、返り値は定数としています。

「class exact」について(二階微分方程式の厳密解を初期条件に合わせた特解を求めることを目的としています、ただし、非斉次方程式の場合の計算についてはもう少し汎用性が出るように現在制作中です)
・「def init」:扱う微分方程式の係数等の条件を設定しています
・「def oed_exact」:初期条件に合わせた特解を求め、その関数自体を返り値に設定しています。また、特解の決定には、特性方程式が重解の時とそうでないときで場合分けをしており、各係数を行列を用いて計算しています。

「class visualize」について(与えられた二つの配列を視覚化することを目的としたクラスです)
・「def init」:扱う配列とデータの名前を取得しています。
・「def plot」:シンプルな二次元のグラフをプロットしています。
・「stakplot」:新たにもう一つ配列(y成分のデータ)を取得し、事前に取得したデータと合わせてプロットしています。
・「def complot」:新たにもう一つ配列(y成分のデータ)を取得し、事前に取得したデータとの差の値に絶対値を取ってプロットしています。
・「def animation」:事前に取得したy軸のデータを利用してアニメーション(詳細にはパラパラ漫画的なもの)を作ってい表示しています。

以上です。実際のコードは次のようになっています。

該当のソースコード(classのファイル)

python3

###クラスなどをまとめているファイルです import matplotlib.pyplot as plt import numpy as np import math class numerical:#定数係数の2階微分方程式を数値的に解くクラス def __init__(self,constant,inhomo,initial,step,Max): #微分方程式の係数,非斉次項,初期条件(独立変数→従属変数の階数が高い順),刻み幅,独立変数の最大値 self.constant=constant self.initial=initial self.step=step self.variable=[round((initial[0]+i*step),5) for i in range(int((Max-initial[0])/step)+1)]#数値解と解析解で扱う微小時間を共通にしたいので最初のうちに計算しておく[内包表記] self.inhomo=inhomo#非斉次項の関数を代入 def derivative(self,x,z,t):#微分係数を微分方程式から求める return -(self.constant[1]*z+self.constant[2]*x-self.inhomo(t))/self.constant[0] def euler(self):#オイラー法 x=[] x.append(self.initial[1])#初期条件の代入 z=self.initial[2] for i in range(1,len(self.variable)): x.append(x[i-1]+z*self.step) z=z+self.derivative(x[i],z,self.variable[i])*self.step return x def exact(self,fanction):#入力された解析解を数値解と同じ条件(レンジや刻み幅)で解く exact=[fanction(i) for i in self.variable] return exact class exact:#定数係数二階微分方程式を初期条件に合わせた特解を求めるクラス、返り値は関数で設定(高階関数:クロージャ)、返す関数の引数と返り値はともに一つのみ def __init__(self,constant,initial,special,derivative):#微分方程式の係数、初期条件、非斉次方程式の特解、特解の微分係数 self.constant=constant self.initial=initial self.special=special self.derivative=derivative def ode_exact(self):#非斉次の二階微分方程式 m,c,k=self.constant[0],self.constant[1],self.constant[2] t0,x0,v0=self.initial[0],self.initial[1],self.initial[2] epsilon=c/pow(4*m*k,0.5) alpha=self.special(t0) beta=self.derivative(t0) if (epsilon!=1):#特性方程式が重解でないとき s1=-epsilon+pow(epsilon**2-1,0.5)#一般解の係数 s2=-epsilon-pow(epsilon**2-1,0.5) A=np.array([[x0-alpha],[v0-beta]]) H=np.array([[1,1],[s1,s2]]) Ad=np.dot(np.linalg.inv(H),A) a=Ad[0][0] b=Ad[1][0] def f(t): return (a*pow(math.e,s1*t)+b*pow(math.e,s2*t)+self.special(t)).real return f else: A=np.array([[x0-alpha],[v0-beta]]) H=np.array([[1,0],[-epsilon*k/m,1]]) Ad=np.dot(np.linalg.inv(H),A) a=Ad[0] b=Ad[1] def f(t): return ((a+b*t)*pow(math.e,-epsilon*k*t/m)+self.special(t)).real return f class visualize:#視覚化するクラス(アニメーションなどを追加予定) def __init__(self,x,y,label):#視覚化したいデータを取得 self.x=x self.y=y self.label=label def plot(self,Xlim,Ylim):#取得したデータをグラフ化 fig=plt.figure() ax=plt.axes(xlim=Xlim,ylim=Ylim) ax.plot(self.x,self.y,label=self.label) plt.legend() plt.show() def stakplot(self,Xlim,Ylim,x2,Label):#既存のデータと入力されたデータを重ねて比較 fig=plt.figure() ax=plt.axes(xlim=Xlim,ylim=Ylim) ax.plot(self.x,self.y,label=self.label) ax.plot(self.x,x2,label=Label) plt.legend() plt.show() def complot(self,Xlim,Ylim,x2,Label):#二つのデータを差を取って比較し、プロット if (len(self.y)!=len(x2)): print("比較する要素数が違います") return c=[] for i in range(len(self.x)): c.append(abs(x2[i]-self.y[i])) fig=plt.figure() ax=plt.axes(xlim=Xlim,ylim=Ylim) ax.plot(self.x,c,label=Label) plt.legend() plt.show() def animetion(self,Xlim,Ylim): fig=plt.figure() ax=plt.axes(xlim=Xlim,ylim=Ylim) ax.set_aspect('equal') circle=plt.Circle((self.y[0],0),0.1,fc="r") ax.add_patch(circle) for i in range(1,len(self.y)-1): circle.center=self.y[i],0.0 plt.pause(0.001) plt.show()

該当のソースコード(mainのファイル)

python3

import odeclass import matplotlib.pyplot as plt import numpy as np import math m,c,k=1,0,1 constant=[m,c,k]#微分方程式の係数:階数が高い変数にかかっている係数から並べる initial=[0,2.0,0]#初期条件:t0,x0,v0 def zero(t): return 0 f=odeclass.exact(constant,initial,zero,zero) fanc=f.ode_exact()#厳密解を初期条件に合わせた特解を決定 step=0.01 z=odeclass.numerical(constant,zero,initial,step,100) t=z.variable x1=z.euler() x2=z.exact(fanc) v=odeclass.visualize(t,x1,"euler") v.complot((0,100),(0,0.1),x2,"error")

###実行結果
オイラー法で解いた数値解と解析解との差の値に対して絶対値を取り、時間との関係を表示しています。

試したこと

実際に手で計算した結果と比較した時も同じ結果が出ました。ここにおける比較対象はオイラー法で解いた数値解以外に、プログラムで決定した特解の数値解両方確認しました。

補足情報

プログラム初心者であり、classの使い方とファイル間のやり取りを練習をしたかったので見にくいプログラムになってしまっています。classの組み方にこうしたほうがもっと見やすくなる、責任を分散のさせ方が悪い部分などのご指摘もあれば合わせてコメントいただけると助かります。よろしくお願いします。

良い質問の評価を上げる

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

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

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

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

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

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

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

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

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

ppaul

2021/04/29 07:05

出力結果の画像を貼り付けていただけませんか。
hayataka2049

2021/04/29 07:23

やっている計算の説明と、どういう結果が期待されているのかを書いて頂いた方が回答がつきやすくなるかな……と思います。
fin1920

2021/04/29 07:25

了解しました。補足します
ppaul

2021/04/29 07:47

うなりというのは、一定周期で周波数が変化しているように見えるという意味ですか?
fin1920

2021/04/29 07:50

厳密にいえばうなりでも何でもないと思うのですが、誤差がある時刻まで収束しそうになって、その後発散するという挙動を見せています。実行結果を添付しましたので確認してもらえるでしょうか

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

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

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

ただいまの回答率
87.20%

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

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

質問する

関連した質問

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

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

アルゴリズム

アルゴリズムとは、定められた目的を達成するために、プログラムの理論的な動作を定義するものです。

プログラミング言語

プログラミング言語はパソコン上で実行することができるソースコードを記述する為に扱う言語の総称です。

文字コード

文字コードとは、文字や記号をコンピュータ上で使用するために用いられるバイト表現を指します。

Python

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