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

Question

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

###クラスなどをまとめているファイルです
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のファイル)

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の組み方にこうしたほうがもっと見やすくなる、責任を分散のさせ方が悪い部分などのご指摘もあれば合わせてコメントいただけると助かります。よろしくお願いします。

Comments

[ppaul]

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

[hayataka2049]

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

[fin1920]

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

[ppaul]

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

[fin1920]

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

Answer 1

解決しました、オイラー法の部分が間違えていました。
修正箇所は以下の通りです

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-1],z,self.variable[i-1])*self.step#この行が修正部分
return x

コメント頂い方たありがとうございました。