Home > Jejak langkah untuk bersekolah lagi, Python, Semester 2, Tutorial > Runge Kutta Implisit dengan python (3)

Runge Kutta Implisit dengan python (3)

Setelah melakukan iterasi dengan metode Newton, sekarang saatnya menggunakan metode Picard [1]. Metode ini jauh lebih sederhana dibanding metode Newton.

Seperti telah dijelaskan sebelumnya, persamaan yang mendasari perhitungan Runge Kutta ada dua, masing-masing adalah persamaan untuk menghitung nilai antara dan persamaan untuk menghitung nilai fungsi. Pada metode Picard, tidak diperlukan mencari matriks Jacobian dari persamaan nilai antara. Cukup menggunakan persamaan nilai antara saja.

Tahap pertama, nilai tebakan awal (yang juga merupakan initial value) dijadikan sebagai nilai antara. Selanjutnya, nilai tebakan awal itu digunakan untuk menghitung nilai antara yang diperbarui. Persamaan nilai antara adalah seperti persamaan (1)

y_{ik} = x_{k} + h \Sigma_{j=i}^{s} a_{ij}f(t_{k}+c_{j}h,y_{jk})      (1)

Dari persamaan (1) kita sekarang memiliki dua nilai antara. Yang pertama adalah nilai antara yang di-assigned dari nilai tebakan awal. Dan yang kedua adalah nilai antara dari hasil memasukkan nilai tebakan awal ke persamaan (1). Selisih nilai antara kedua dan pertama dibandingkan. Jika lebih kecil atau sama dengan nilai toleransi, maka perhitungan nilai antara untuk step pertama dapat dilanjutkan dengan perhitungan nilai fungsi. Jika masih lebih besar, maka nilai antara kedua menjadi nilai antara pertama. Sedangkan nilai antara kedua dihitung ulang dari nilai antara pertama yang diperbarui. Hal ini dilakukan berulang hingga seluruh rentang waktu diakomodasi. Berikut adalah program perhitungan Runge Kutta dengan metode iterasi Picard. Seperti pada iterasi Newton, program ini dibuat dengan asumsi jika iterasi pada suatu step mencapai 5000, sedangkan selisih nilai antara pertama dan kedua masih lebih besar dari nilai toleransi, maka kondisi tersebut mengindikasikan tidak tercapainya konvergensi, sehingga perhitungan harus dihentikan.


from numpy import *
from numpy import linalg as la
import os.path
from time import *
from Equations import *
from Butcher import *
from Symbols import *
class PicardRKSolver():
	def __init__(self,Eq,Symbols,Butcher,Initial,Stop,StepSize,TimeTarget):
		self.errorStatus=0
		
		"""lenEq: banyaknya persamaan"""
		self.lenEq=len(Eq)
		
		"""persamaan yang akan diselesaikan"""
		self.equations=Eq
		
		"""simbol/variabel yang digunakan"""
		self.symbols=Symbols
		
		"""orde: orde penyelesaian"""
		self.orde=len(Butcher[1])
		
		"""x: nilai target yang akan dihitung pada setiap step size"""
		self.x1=array(Initial)

		"""y: nilai antara dalam bentuk sub array untuk menghitung derivatif"""
		"""y1: nilai antara, yang tebakan awalnya disamakan dengan nilai awal"""
		tempY=[]
		for i in range(self.orde):
			tempY.append(Initial)
		self.y=array(tempY)
		self.y1=array(tempY)
		self.y1=reshape(self.y1,(self.orde*self.lenEq,1))
		
		"""y2: nilai antara yang diperbarui"""
		self.y2=arange(1)
		
		"""delta: perbaikan nilai antara"""
		self.delta=self.y1
		
		"""a,b,c: komponen array Butcher"""
		self.a=Butcher[0]
		self.b=Butcher[1]
		self.c=Butcher[2]
		
		"""
		self.stop: stopping criteria
		self.h=step size
		"""
		self.stop=Stop
		self.h=StepSize
		self.t=TimeTarget
		self.selisih=self.t/self.h-int(self.t/self.h)
		self.stop=Stop
		self.h=StepSize
		self.t=TimeTarget
		self.selisih=self.t/self.h-int(self.t/self.h)
		
		self.__setY()
		
		self.elapsed=0
		self.start=time()
		if self.selisih==float(0):
			i=0
			while(i<int(self.t/self.h) and self.errorStatus==0):
				self.__iterate()
				self.__setY()
				i+=1
			self.elapsed=time()-self.start
		else:
			i=0
			while(i<int(self.t/self.h) and self.errorStatus==0):
				self.__iterate()
				self.__setY()
				i+=1
			self.h=self.selisih
			i=int(self.t/self.h)+1
			self.__extraIterate()
			self.elapsed=time()-self.start

		self.namafile.close()
		if self.errorStatus==1:
			print 'Quit unseccessful in ' + str(self.elapsed) + 'seconds ... '
		else:
			print 'Computation done ...'	

	def __setY(self):
		l=0
		for i in range(self.orde):
			for j in range(self.lenEq):
				self.y[i][j]=self.y1[l]
				l+=1
		p=[]
		l=0
		for i in range(self.orde):
			for j in range(self.lenEq):
				q=0
				for k in range(self.orde):
					q=q+(self.a[i][k]*self.y[k][j])
				p.append(self.x1[j]+q*self.h)
				l+=1
		self.y2=array(p)
		self.y2=reshape(self.y2,(self.orde*self.lenEq,1))

		
	def __iterate(self):
		l=0
		self.delta=self.y2-self.y1
		self.y1=self.y2
		l+=1
		while(self.stop<linalg.norm(self.delta) and l=5000:
			self.errorStatus=1

	def __setX(self):
		self.y1=reshape(self.y1,(self.orde,self.lenEq))
		for i in range(self.lenEq):
			tempY=0
			for j in range(self.orde):
				tempY=tempY+(self.b[j]*self.y1[j][i])
			self.x1[i]=self.x1[i]+(self.h*tempY)
		self.y1=reshape(self.y1,(self.orde*self.lenEq,1))

Referensi
[1] R. Holsapple, R. Iyer, and D. Doman, Variable step-size selection methods for implicit integration schemes for odes, International Journal of Numerical Analysis and Modeling, vol. 4, no. 2, pp. 212–242, 2007.

About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: